Python绘制RTKLIB的POS文件中的XYZ误差曲线及计算RMSE

前面博客中有提到过，文中给出了部分代码，在这篇博客里面，给出剩下的代码

# -*- coding:utf-8 -*-
"""

created on Fri Apr 29 08:26:10 2022
@author: xymeng

"""
import numpy as np
import math as mt
import matplotlib.pyplot as plt
import sklearn.metrics as ms

sqrt = mt.sqrt
'''以下为列表定义'''
X = []
Y = []
Z = []#定位结果
X0 = []
Y0 = []
Z0 = []#参考真值
delX = []
delY = []
delZ = []#两者差值
XT=[]
YT=[]
ZT=[]#将真值放到列表中

if __name__ == '__main__':

    path = input('请输入pos定位结果文件：') #E://rtklib//rtklib-test1//不同方案对比数据//Huber权函数抗差最小二乘//MA36-x3-C-huber1.pos
    path0 = input('请输入pos参考真值文件：')#E://rtklib//rtklib-test1//data//MA36x3-ppp.pos

    '''以下代码是读取定位结果所有行的XYZ'''
    with open(path,'r') as f:
        lines = f.readlines()
        for line in lines[0:len(lines)]:
            if (line.strip('\n').split(' ')[0])!='%' :
                X.append(float(line.split('  ')[1]))
                Y.append(float(line.split('  ')[2]))
                Z.append(float(line.split('  ')[3]))

    '''以下代码是读取参考真值所有行的XYZ'''
    with open(path0, 'r') as f:
        lines = f.readlines()
        for line in lines[0:len(lines)]:
            if (line.strip('\n').split(' ')[0]) != '%':
                X0.append(float(line.split('  ')[1]))
                Y0.append(float(line.split('  ')[2]))
                Z0.append(float(line.split('  ')[3]))

    '''以下行是将定位均值设为XYZ真值'''
    Xt = np.mean(X0)
    Yt = np.mean(Y0)
    Zt = np.mean(Z0)

    '''以下行是将XYZ真值存入列表中，用于符合rmse计算结构'''
    for k in range(len(X)):
        XT.append(Xt)
        YT.append(Yt)
        ZT.append(Zt)

    '''以下是计算定位值与真值差值，作为绘制的纵坐标'''
    for i in range(len(X)):
        delX.append(X[i]-Xt)
        delY.append(Y[i]-Yt)
        delZ.append(Z[i]-Zt)

    '''下面为计算del的均值、标准差和均方根误差'''
    Xmean = round(np.mean(delX), 4)
    Ymean = round(np.mean(delY), 4)
    Zmean = round(np.mean(delZ), 4)

    Xstd = round(np.std(delX, ddof=1), 4)
    Ystd = round(np.std(delY, ddof=1), 4)
    Zstd = round(np.std(delZ, ddof=1), 4)

    Xrmse = round(sqrt(ms.mean_squared_error(XT, X)), 6)
    Yrmse = round(sqrt(ms.mean_squared_error(YT, Y)), 6)
    Zrmse = round(sqrt(ms.mean_squared_error(ZT, Z)), 6)

    '''以下为定义收敛图的x坐标范围'''
    xx = [i for i in range(1,len(X)+1)]

    '''以下为绘制误差曲线图'''
    #plt.style.use('seaborn-whitegrid')
    plt.rcParams['font.sans-serif'] = ['KaiTi']
    plt.rcParams['axes.unicode_minus'] = False
    plt.figure(dpi=500)#, figsize=(12, 6)
    plt.subplot(311)
    plt.plot(xx,delX,label='AVE={}m,STD={}m,RMSE={}m'.format(Xmean,Xstd,Xrmse),c='k')
    plt.xlabel('历元/s')
    plt.ylabel('X/m')
    y_ticks = np.arange(-10, 11 ,5)
    plt.yticks(y_ticks)
    #plt.legend()

    plt.subplot(312)
    plt.plot(xx, delY, label='AVE={}m,STD={}m,RMSE={}m'.format(Ymean, Ystd, Yrmse),c='k')
    plt.xlabel('历元/s')
    plt.ylabel('Y/m')
    y_ticks = np.arange(-10, 11, 5)
    plt.yticks(y_ticks)
    #plt.legend()

    plt.subplot(313)
    plt.plot(xx, delZ, label='AVE={}m,STD={}m,RMSE={}m'.format(Zmean, Zstd, Zrmse),c='k')
    plt.xlabel('历元/s')
    plt.ylabel('Z/m')
    #plt.legend()
    y_ticks = np.arange(-10, 11, 5)
    plt.yticks(y_ticks)
    plt.show()