【Python】RTKLIB POS文件转换至ENU误差并绘图

#!/usr/bin/python
# coding=utf-8
'''
RTKLIB POS文件转换至ENU误差并绘图
Version：1.0
Author:LZ_CUMT
脚本功能：
1：批量读取rtklib解算结果文件（.pos)，并将XYZ或LLH坐标转换为ENU形式
2：转换时选取的坐标基准通过snx文件进行读取，若snx文件中无当前测站或未找到snx文件，以定位结果的均值作为替代
3：转换完成后将ENU输出至csv格式文件保存，通过matplotlib库（需额外安装）绘制ENU定位误差图
'''import os
import re
import csv
from math import sin,cos,atan,atan2,pi,sqrt
import matplotlib.pyplot as plt
# from tkinter import Tk
# from tkinter import filedialog# 找目录下含该关键词的文件名及其路径
def findfile(filepath,key_word):filename = []for files in os.listdir(filepath):  # 遍历目录文件名if (key_word == 'pos' ):if re.match(r'[A-Za-z0-9\.]+\.pos',files):filename.append(files)  # 文件名及其路径添加到数组elif (key_word == 'snx'):if re.match(r'igs[0-9]+\.snx',files):filename.append(files)  # 文件名及其路径添加到数组return filename  # 返回数组# 经纬度转化为xyz(参考RTKLIB)
def llh2xyz(lat,lon,h):RE_WGS84 = 6378137.0FE_WGS84 = 1.0/298.257223563lat = lat * pi / 180lon = lon * pi / 180sinp = sin(lat)cosp = cos(lat)sinl = sin(lon)cosl = cos(lon)e2 = FE_WGS84*(2.0-FE_WGS84)v = RE_WGS84 / sqrt(1.0 - e2 * sinp * sinp)x = (v+h)*cosp*cosly = (v+h)*cosp*sinlz = (v*(1.0-e2)+h)*sinpreturn [x,y,z]# 从pos文件中读取坐标结果，并以xyz坐标形式输出
def readpos(posfilepath) :f = open(posfilepath,encoding='gb18030', errors='ignore')ln = f.readline()while ln:if '%  GPST' in ln:if 'x-ecef(m)' in ln:posmode = 'xyz'elif 'latitude(deg)' in ln:posmode = 'llh'breakln = f.readline()xyz = []while ln:ln = f.readline()if not ln:break#content = ln.split(' ')if posmode == 'xyz':x = float(ln[24:38])y = float(ln[39:53])z = float(ln[54:69])xyz.append([x,y,z])elif posmode == 'llh':lat = float(ln[24:38])lon = float(ln[38:53])h = float(ln[53:64])xyz.append(llh2xyz(lat,lon,h))f.close()return xyz# 根据测站id在snx文件中查找测站精确坐标
def getcrd(siteid,snxfilepath,xyz):snxcrd=[]if snxfilepath == '':print('[WARNING] Not find the snxfile, use the average pos as substitute')else:f = open(snxfilepath,encoding='gb18030', errors='ignore')ln = f.readline()while ln:ln = f.readline()if not ln:print('[WARNING] Not find the siteid',siteid,', use the average pos as substitute')breakif ln[14:18]==siteid:snxcrd.append(float(ln[47:68]))ln = f.readline()snxcrd.append(float(ln[47:68]))ln = f.readline()snxcrd.append(float(ln[47:68]))breakf.close()if snxcrd == []:x = 0y = 0z = 0for i in range(len(xyz)):x = x + xyz[i][0]y = y + xyz[i][1]z = z + xyz[i][2]snxcrd = [x/len(xyz),y/len(xyz),z/len(xyz)]else:print('[INFO] Find the sitecrd in the snx file')print('[INFO] The sitecrd is', snxcrd)return snxcrd# xyz转换为llh（经纬度）
def xyz2llh(ecef):aell = 6378137 # WGS - 84 not ITRFfell = 1 / 298.257223563deg = pi / 180u = ecef[0]v = ecef[1]w = ecef[2]esq = 2*fell-fell*fellif w == 0:lat = 0else :lat0 = atan(w/(1-esq)*sqrt(u*u+v*v))j = 0delta = 10^6limit = 0.000001/3600*degwhile delta > limit:N = aell / sqrt(1 - esq * sin(lat0)*sin(lat0))lat = atan((w / sqrt(u*u + v*v)) * (1 + (esq * N * sin(lat0) / w)))delta = abs(lat0 - lat)lat0 = latj = j + 1if j > 10:breaklong = atan2(v,u)h=(sqrt(u*u+v*v)/cos(lat))-Nllh=[lat,long,h]return llh# xyz转换至以测站精确坐标为基准的enu坐标
def xyz2enu(xyz,snxcrd):enu=[]llhcrd = xyz2llh(snxcrd)phi = llhcrd[0] * pi / 180lam = llhcrd[1] * pi / 180sinphi = sin(phi)cosphi = cos(phi)sinlam = sin(lam)coslam = cos(lam)for i in range(len(xyz)):difxyz=[xyz[i][0]-snxcrd[0],xyz[i][1]-snxcrd[1],xyz[i][2]-snxcrd[2]]e = -sinlam*difxyz[0]+coslam*difxyz[1]n = -sinphi*coslam*difxyz[0]-sinphi*sinlam*difxyz[1]+cosphi*difxyz[2]u =  cosphi*coslam*difxyz[0]+cosphi*sinlam*difxyz[1]+sinphi*difxyz[2]enu.append([e,n,u])return enu# 保存enu误差结果为CSV文件
def saveenu(enu,posfilepath):f = open(posfilepath[:-4] + '.csv', 'a', newline="")csv_write = csv.writer(f)for i in range(len(enu)):csv_write.writerow([enu[i][0],enu[i][1],enu[i][2]])f.close()return 0# 绘制enu误差时间序列图
def plotenu(enu,posfilepath):e = []n = []u = []for i in range(len(enu)):e.append(enu[i][0])n.append(enu[i][1])u.append(enu[i][2])fig = plt.figure()ax = fig.add_subplot(1, 1, 1)plt.axis([0, len(enu)+1, -1, 1])ax.plot(range(1,len(enu)+1), e)ax.plot(range(1,len(enu)+1), n)ax.plot(range(1,len(enu)+1), u)ax.grid(True, linestyle='--')plt.xlabel('[Epoch]')plt.ylabel('[Error]')plt.legend(['E','N','U'])#plt.show()plt.savefig(posfilepath[:-4] + '.png', dpi=400)plt.close()return 0if __name__ == '__main__':# window = Tk()# window.withdraw()# filepath = filedialog.askdirectory(title='选择数据所在文件夹')# 此处修改文件夹路径，将要处理的pos文件和下载的snx坐标文件统一放于此文件夹下filepath = 'D:\\data\\all'posfilelist = findfile(filepath,'pos')snxfilelist = findfile(filepath,'snx')if snxfilelist == []:snxfilepath = ''else:snxfilepath = filepath + '\\' + snxfilelist[0]print('[INFO] Find the snxfile :',snxfilelist[0])for posfile in posfilelist:fileindex = posfilelist.index(posfile) + 1posfilepath = filepath + '\\' + posfileprint('[INFO] Process the '+str(fileindex)+'th posfile :',posfile)siteid = posfile[0:4].upper()print('[INFO] The siteid is :',siteid)xyz = readpos(posfilepath)snxcrd = getcrd(siteid,snxfilepath,xyz)enu = xyz2enu(xyz,snxcrd)saveenu(enu,posfilepath)plotenu(enu,posfilepath)print('[INFO] Finish the '+str(fileindex)+'th posfile processing successfully')

【Python】RTKLIB POS文件转换至ENU误差并绘图相关推荐

python将doc文件转换docx
一:python将doc文件批量转换docx: # -*- coding:utf-8 -*- # @Author : yyzhang import os import time from win32c ...
基于Python的DICOM文件转换教程，使用pydicom将图片存为DICOM文件。
基于Python的DICOM文件转换教程,使用pydicom将图片存为DICOM文件. DICOM是医学图像和信息的数字化标准,可用于将医学影像数据.诊断报告等信息在医疗领域进行传输.分享和分析.而常 ...
Python中将字节流文件转换成图片文件
Python中将字节流文件转换成图片文件 import urllib3 import os #PIL图像处理标准库 from PIL import Image from io import Bytes ...
用Python将word文件转换成html（转）
用Python将word文件转换成html 序最近公司一个客户大大购买了一堆医疗健康方面的科普文章,希望能放到我们正在开发的健康档案管理软件上.客户大大说,要智能推送!要掌握节奏!要深度学习!要让用 ...
python word处理_妙用Python将word文件转换成html 方法超简单
什么方法可以将word文件转换成html,找了一圈,没有发现合适的应用可以把word或indd转化成干净的html.机缘巧合,无意间听说python很擅长文本处理,用Python将word文件转换成h ...
Python将PDF文件转换成PNG的方案
2019独角兽企业重金招聘Python工程师标准>>> 目前最靠谱的是基于 mupdf 的 Python 绑定: https://github.com/rk700/PyMuPDF ...
python把csv文件转换txt_Python实现txt文件转csv格式
码农公社 210.net.cn 210= 1024 10月24日一个重要的节日--码农(程序员)节把txt文件转成成csv文件格式,通过手动打开excel文件,然后导入txt来生产csv文件. ...
Python将.dat文件转换成.csv文件
在找数据的时候有时候会找到.dat文件,我发现了两种方式. 第一种很简单,是利用可以打开dat的软件进行转换,但是这种方式在只有一个数据的时候可以使用,数据集多时很不方便. 另一种就是使用python ...
用Python将word文件转换成html
序最近公司一个客户大大购买了一堆医疗健康方面的科普文章,希望能放到我们正在开发的健康档案管理软件上.客户大大说,要智能推送!要掌握节奏!要深度学习!要让用户留恋网站无法自拔! 话说符合以上特点的我也 ...

【Python】RTKLIB POS文件转换至ENU误差并绘图

【Python】RTKLIB POS文件转换至ENU误差并绘图相关推荐

最新文章

热门文章