事情是这样的,最近博主在做一个测量接收机RTK定位精度的实验,在野外采集到的数据是NMEA的GPGGA格式的定位信息,需要将其转换为处理软件需要用到的格式,因此编写了一个小的转换程序,下面详细介绍:
在野外采集到的数据如下所示:
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
$GPGGA,031851.00,1.235,N,1.235,E,4,38,0.4,1.235,M,11.785,M,1.0,1029*4C
需要将其转换为如下格式:
-2633658.0904 3758412.8295 4261486.8771
-2633658.0904 3758412.8295 4261486.8771
-2633658.0904 3758412.8295 4261486.8771
-2633658.0904 3758412.8295 4261486.8771
-2633658.0904 3758412.8295 4261486.8771
-2633658.0904 3758412.8295 4261486.8771
-2633658.0904 3758412.8295 4261486.8771
-2633658.0904 3758412.8295 4261486.8771
-2633658.0904 3758412.8295 4261486.8771
-2633658.0904 3758412.8295 4261486.8771
-2633658.0904 3758412.8295 4261486.8771
-2633658.0904 3758412.8295 4261486.8771
上面是测量点位的XYZ信息,需要将原始数据中的BLH数据信息提取出来,再利用BLH->XYZ的转换公式,将其变为目标格式。话不多说,上代码:
# -*- coding:utf-8 -*-
"""
created on Sat Jun 5 10:36:25 2021
@author: xymeng
"""
import math as mt
'''以下为WGS84椭球参数'''
a = 6378137.0 # 参考椭球的长半轴, 单位 m
b = 6356752.31414 # 参考椭球的短半轴, 单位 m
'''以下为三角函数调用'''
sqrt = mt.sqrt
sin = mt.sin
cos = mt.cos
atan = mt.atan
'''以下为弧度转换角度函数'''
def rad2angle(r):
"""
该函数可以实现弧度到角度的转换.
:param r: 弧度
:return: a, 对应的角度
"""
a = r*180.0/mt.pi
return a
'''以下为角度转换弧度函数'''
def angle2rad(a):
"""
该函数可以实现角度到弧度的转换.
:param a: 角度
:return: r, 对应的弧度
"""
r = a*mt.pi/180.0
return r
'''以下为BLH转XYZ函数'''
def BLH2XYZ(B, L, H, a, b):
"""
该函数实现把某点在参心空间直角坐标系下的坐标(X, Y, Z)转为大地坐标(B, L, H).
:param X: X方向坐标,单位 m
:param Y: Y方向坐标, 单位 m
:param Z: Z方向坐标, 单位 m
:param a: 地球长半轴,即赤道半径,单位 m
:param b: 地球短半轴,即大地坐标系原点到两级的距离, 单位 m
:return: B, L, H, 大地纬度、经度、海拔高度 (m)
"""
B = angle2rad(B)
L = angle2rad(L)
e = sqrt((a**2-b**2)/(a**2)) #椭球第一偏心率
N = a / sqrt((1 - e * e * sin(B) * sin(B)))
X = (N+H)*cos(B)*cos(L)
Y = (N+H)*cos(B)*sin(L)
Z = (N*(1-e *e)+H)*sin(B)
return X, Y, Z # 返回大地纬度B、经度L、海拔高度H (m)
if __name__ == '__main__':
path = input('请输入要转换的文本:')
with open(path,'r') as f:
data = f.readlines()
i = 0
while i<len(data):
string = data[i]
Bdeg = string[17:19]
Bdeg = int(Bdeg)
Bmin = string[19:29]
Bmin = float(Bmin)
Ldeg = string[32:35]
Ldeg = int(Ldeg)
Lmin = string[35:45]
Lmin = float(Lmin)
H = string[57:64]
H = float(H)
B = Bdeg+Bmin/60
L = Ldeg+Lmin/60
X,Y,Z = BLH2XYZ(B,L,H,a,b)
X = str(X)
Y = str(Y)
Z = str(Z)
with open(path+'geokit.txt','a+') as f1:
f1.write(X+' '+Y+' '+Z+'\n')
i = i+1
print('转换成功!')
代码中主要功能的实现便是如下几行代码:
path = input('请输入要转换的文本:')
with open(path,'r') as f:
data = f.readlines()
with open(path+'geokit.txt','a+') as f1:
f1.write(X+' '+Y+' '+Z+'\n')
文章到此结束啦,博主最近任务比较重,暂时没有很充足的时间更一些内容,待这一阵事情处理完,定会为朋友们带来更多文章!感谢阅读,感谢交流!