空間直角坐標轉高斯克呂格平面坐標

• 原理示意
• 源碼

有時需要將點展開到地圖上，需要的是平面坐標，而手中只有三維的空間直角坐標XYZ或者BLH時該怎么辦呢？這時候就需要將坐標進行投影轉換，本文給出的程序代碼是用高斯克呂格投影進行坐標轉換

原理示意

源碼

源碼中還有其他坐標轉換的小函數，也可參考

# -*- coding:utf-8 -*- """created on Wed Apr 29 20:30:25 2021 @author: xymeng""" import numpy as np import math as mt import matplotlib.pyplot as plt import sklearn.metrics as ms import csv import pandas as pd'''以下為三角函數調用''' sqrt = mt.sqrt sin = mt.sin cos = mt.cos atan = mt.atan tan = mt.tan'''以下為WGS84橢球參數''' a = 6378137.0 # 參考橢球的長半軸, 單位 m b = 6356752.31414 # 參考橢球的短半軸, 單位 m e = sqrt((a**2-b**2)/(a**2)) #橢球第一偏心率 e1 = sqrt((a**2-b**2)/(b**2)) #橢球第二偏心率'''以下為高斯投影中的參數''' A1 = 1 + 3/4 * e**2 + 45/64 * e**4 +175/256 * e**6 +11025/16384 * e**8 +43659/65536 * e**10 B1 = 3/4 * e**2 + 15/16 * e**4 +525/512 * e**6 + 2205/2048 * e**8 + 72765/65536 * e**10 C1 = 15/64 * e**4 +105/256 * e**6 + 2205/4096 * e**8 + 10395/16384 * e**10 D1 = 35/512 * e**6 + 315/2048 * e**8 + 31185/131072 * e**10 E1 = 315/16384 * e**8 + 3645/65536 * e**10''''以下為弧度轉換角度函數''' def rad2angle(r):"""該函數可以實現弧度到角度的轉換.:param r: 弧度:return: a, 對應的角度"""a = r*180.0/mt.pireturn a'''以下為角度轉換弧度函數''' def angle2rad(a):"""該函數可以實現角度到弧度的轉換.:param a: 角度:return: r, 對應的弧度"""r = a*mt.pi/180.0return r'''以下為XYZ轉BLH函數''' def XYZ2BLH(X, Y, Z, 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)"""e = sqrt((a**2-b**2)/(a**2)) #橢球第一偏心率if X == 0 and Y > 0:L = 90elif X == 0 and Y < 0:L = -90elif X < 0 and Y >= 0:L = atan(Y/X)L = rad2angle(L)L = L+180elif X < 0 and Y <= 0:L = atan(Y/X)L = rad2angle(L)L = L-180else:L = atan(Y / X)L = rad2angle(L)b0 = atan(Z/sqrt(X**2+Y**2))N_temp = a/sqrt((1-e*e*sin(b0)*sin(b0)))b1 = atan((Z+N_temp*e*e*sin(b0))/sqrt(X**2+Y**2))while abs(b0-b1) > 1e-7:b0 = b1N_temp = a / sqrt((1 - e * e * sin(b0) * sin(b0)))b1 = atan((Z + N_temp * e * e * sin(b0)) / sqrt(X ** 2 + Y ** 2))B = b1N = a/sqrt((1-e*e*sin(B)*sin(B)))H = sqrt(X**2+Y**2)/cos(B)-NB = rad2angle(B)return B, L, H # 返回大地緯度B、經度L、海拔高度H (m)'''以下為高斯克呂格投影6度分帶法，經度Lon與投影帶號Zone以及帶號Zone與投影帶中央子午線經度L0的計算''' def GK(L):'''高斯克呂格投影6度分帶法，經度Lon與投影帶號Zone以及帶號Zone與投影帶中央子午線經度L0的計算:param L:弧度制:return L0:'''if L>=0:'''東半球時'''z = int(L/6) + 1L0 = (6 * (z) - 3) * mt.pielse:'''西半球時'''z = int(L/6) + 60L0 = ((6 * z - 3) - 360)* mt.pireturn L0'''以下為將大地坐標轉換為高斯平面坐標函數''' def GKxy(B,L,a,b,A1,B1,C1,D1,E1):'''將大地坐標轉換為高斯平面坐標:param B::param L::param a::param b::param A1::param B1::param C1::param D1::param E1::return x,y:'''X0 = a * (1 - e**2)*(A1 * B - (B1/2)*sin(2*B) + (C1/4)*sin(4*B) - (D1/6)*sin(6*B) + (E1/8)*sin(8*B))N = a / sqrt(1 - e*e *sin(B)*sin(B))t = tan(B)L0 = GK(L)l = L - L0y2 = e1*e1*cos(B)*cos(B)x = X0 + (N/2)*sin(B)*cos(B)*l*l +(N/24)*sin(B)*(cos(B)**3)*(5-t**2 +9*y2 + 4*y2*y2)*(l**4) + (N/720)*sin(B)*(cos(B)**5)*(61 - 58*t*t + (t)**4)*(l**6)y = N*cos(B)*l + (N/6)*(cos(B)**3)*(1-t**2 + y2)*(l**3) + (N/120)*(cos(B)**5)*(5-18*(t**2)+(t**4)+14*y2-58*y2*(t**2))*(l**5)return x,y'''以下定義空間直角坐標的存儲列表''' X = [] Y = [] Z = [] '''以下定義大地坐標存儲列表''' B = [] L = [] H = [] '''以下定義高斯平面坐標存儲列表''' x = [] y = []if __name__ == '__main__':path = input('請輸入空間直角坐標源文件：') #C://Users//23646//Downloads//數據//IGSNetwork.csvwith open(path,'r') as f:lines = csv.reader(f)lines = list(lines)for line in lines[1:]:X.append(float(line[1]))Y.append(float(line[2]))Z.append(float(line[3]))for i in range(len(X)):B0,L0,H0 = XYZ2BLH(X[i],Y[i],Z[i],a,b)B.append(B0 * mt.pi/180)L.append(L0 * mt.pi/180)H.append(H0 * mt.pi/180)'''以下是將BLH的弧度數寫入到csv文件中'''savepath = 'C://Users//23646//Downloads//數據//IGSBLH-rad.csv'data = pd.DataFrame({'B':B,'L':L,'H':H})data.to_csv(savepath,index=False)'''以下為執行大地坐標轉高斯平面坐標代碼塊'''for k in range(len(B)):x0 ,y0 = GKxy(B[k],L[k],a,b,A1,B1,C1,D1,E1)x.append(x0)y.append(y0)'''以下是將xy的高斯平面坐標寫入到csv文件中'''savepath1 = 'C://Users//23646//Downloads//數據//IGSxy.csv'data = pd.DataFrame({'x': x, 'y': y})data.to_csv(savepath1, index=False)

總結

以上是生活随笔為你收集整理的Python实现空间直角坐标转高斯克吕格平面坐标的全部內容，希望文章能夠幫你解決所遇到的問題。

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