代码中有详细的注释
根据《大地测量学（第三版）》武汉大学出版社的教材进行编写
import math as mimport sys as simport osclass Earth():# 建立椭球类，采用CGCS2000中国大地坐标系数据# 《大地测量学基础》P.106的表4-1def __init__(self): # 给定椭球数据# 长半轴self.a = 6378137.0# 短半轴self.b = 6356752.31414# 扁率self.f = 1 / 298.257222101# 第一偏心率的平方self.e1s = m.pow(0.0818191910428, 2)# 第二偏心率的平方self.e2s = self.e1s / (1 - self.e1s)class Transform():# 建立转换类def __init__(self): # 导入椭球数据self.earth = Earth()def BLH_XYZ(self, B, L, H): # 大地坐标转换为空间直角坐标# 将角度转换为弧度self.b = m.radians(B)self.l = m.radians(L)# 这里用个python的包就很简便的实现了角度到弧度间的转换# 计算辅助函数self.W = m.sqrt(1 - self.earth.e1s * m.pow(m.sin(self.b), 2))#《大地测量学基础》P.106的（4-5）# 转换为空间直角坐标# 《大地测量学基础》P.111的（4-30）self.N = self.earth.a / self.Wself.X = (self.N + H) * m.cos(self.b) * m.cos(self.l)self.Y = (self.N + H) * m.cos(self.b) * m.sin(self.l)self.Z = (self.N * (1 - self.earth.e1s) + H) * m.sin(self.b)def XYZ_BLH(self, X, Y, Z): # 空间直角坐标转换为大地坐标# 求出大地经度self.l = m.atan(Y / X) #《大地测量学基础》P.111的（4-31）# 求出大地纬度# 《大地测量学基础》P.112self.r = m.sqrt(X * X + Y * Y)self.tb1 = Z / self.r #《大地测量学基础》P.112的（4-32）中的sinB的B取0得到 self.e1s是第一偏心率的平方while True:'''self.tb2 = 1 / self.r * (Z + self.earth.a * self.earth.e1s * self.tb1 /m.sqrt(1 + self.tb1 * self.tb1 * (1 - self.earth.e1s))) #《大地测量学基础》P.112的（4-38）的变式'''# 这里的俩个方法的tb2的结果是一样，本质是公式的转变，下面这个是课本的self.tb2 = 1 / self.r * (Z + self.earth.a * self.earth.a * self.earth.e1s * self.tb1 /(self.earth.b * m.sqrt(1 + self.earth.e2s + self.tb1 * self.tb1)))# 《大地测量学基础》P.112的（4-38）if abs(self.tb2 - self.tb1) <= 5e-10: # 对于这个范围特地查找了下资料breakself.tb1 = self.tb2self.b = m.atan(self.tb2) # 利用迭代的方法求出B的值，直到符合要求时返回# 求出大地高self.W = m.sqrt(1 - self.earth.e1s * m.pow(m.sin(self.b), 2))self.N = self.earth.a / self.Wself.H = self.r / m.cos(self.b) - self.N #《大地测量学基础》P.112的（4-35）# 将弧度转换为角度self.B = m.degrees(self.b)self.L = m.degrees(self.l)# 这里用个python的包就很简便的实现了弧度到角度间的转换class Point():# 建立点类def BLH(self, B, L, H): # 建立大地坐标点类self.B = Bself.L = Lself.H = Hself.P = Transform()self.P.BLH_XYZ(B, L, H)self.X = self.P.Xself.Y = self.P.Yself.Z = self.P.Zdef XYZ(self, X, Y, Z): # 建立直角坐标点类self.X = Xself.Y = Yself.Z = Zself.P = Transform()self.P.XYZ_BLH(X, Y, Z)self.B = self.P.Bself.L = self.P.Lself.H = self.P.HSolveType = input('请输入解算类型：（键入A为大地坐标->直角坐标，B为直角坐标->大地坐标）\n')while True:if SolveType == 'A':# 当选择'大地坐标->直角坐标'模式时# 输入大地坐标I = input('请输入大地坐标B、L、H，中间用空格隔开：\n')while True:if I:passelse:s.exit(0)try:[B, L, H] = [float(n) for n in I.split()]breakexcept:I = input('输入错误，请重新输入大地坐标B、L、H，中间用空格隔开：\n')else:pass# 开始转换P = Point() # 创建一个实例对象P.BLH(B, L, H)B = str(B)L = str(L)H = str(H)X = str(P.X)Y = str(P.Y)Z = str(P.Z)massege = '\n您输入的大地坐标为：（' + B + ', ' + L + ', ' + H + '）\n转换得到直角坐标为：（' + X + ', ' + Y + ', ' + Z + '）'print(massege)breakelif SolveType == 'B':# 当选择'直角坐标->大地坐标'模式时# 输入直角坐标I = input('请输入直角坐标X、Y、Z，中间用空格隔开：\n')while True:if I:passelse:s.exit(0)try:[X, Y, Z] = [float(n) for n in I.split()]1/Xbreakexcept:I = input('输入错误，请重新输入直角坐标X、Y、Z，中间用空格隔开：\n')else:pass# 开始转换P = Point() # 创建一个实例对象P.XYZ(X, Y, Z)X = str(X)Y = str(Y)Z = str(Z)B = str(P.B)L = str(P.L)H = str(P.H)massege = '\n您输入的直角坐标为：（' + X + ', ' + Y + ', ' + Z + '）\n转换得到大地坐标为：（' + B + ', ' + L + ', ' + H + '）'print(massege)breakelse:SolveType = input('输入错误，请重新输入解算类型：（键入C为正算，D为反算）\n')os.system('pause') # 暂停命令
运行结果示例：
