代碼中有詳細的注釋

根據《大地測量學（第三版）》武漢大學出版社的教材進行編寫

import math as m import sys as s import 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') # 暫停命令

運行結果示例：

?

總結

以上是生活随笔為你收集整理的坐标转换（空间直角坐标系与大地坐标系）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。