空間直角坐標(biāo)轉(zhuǎn)高斯克呂格平面坐標(biāo)

• 原理示意
• 源碼

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

原理示意

源碼

源碼中還有其他坐標(biāo)轉(zhuǎn)換的小函數(shù)，也可參考

# -*- 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'''以下為三角函數(shù)調(diào)用''' sqrt = mt.sqrt sin = mt.sin cos = mt.cos atan = mt.atan tan = mt.tan'''以下為WGS84橢球參數(shù)''' a = 6378137.0 # 參考橢球的長(zhǎng)半軸, 單位 m b = 6356752.31414 # 參考橢球的短半軸, 單位 m e = sqrt((a**2-b**2)/(a**2)) #橢球第一偏心率 e1 = sqrt((a**2-b**2)/(b**2)) #橢球第二偏心率'''以下為高斯投影中的參數(shù)''' 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''''以下為弧度轉(zhuǎn)換角度函數(shù)''' def rad2angle(r):"""該函數(shù)可以實(shí)現(xiàn)弧度到角度的轉(zhuǎn)換.:param r: 弧度:return: a, 對(duì)應(yīng)的角度"""a = r*180.0/mt.pireturn a'''以下為角度轉(zhuǎn)換弧度函數(shù)''' def angle2rad(a):"""該函數(shù)可以實(shí)現(xiàn)角度到弧度的轉(zhuǎn)換.:param a: 角度:return: r, 對(duì)應(yīng)的弧度"""r = a*mt.pi/180.0return r'''以下為XYZ轉(zhuǎn)BLH函數(shù)''' def XYZ2BLH(X, Y, Z, a, b):"""該函數(shù)實(shí)現(xiàn)把某點(diǎn)在參心空間直角坐標(biāo)系下的坐標(biāo)（X, Y, Z)轉(zhuǎn)為大地坐標(biāo)（B, L, H).:param X: X方向坐標(biāo)，單位 m:param Y: Y方向坐標(biāo), 單位 m:param Z: Z方向坐標(biāo), 單位 m:param a: 地球長(zhǎng)半軸，即赤道半徑，單位 m:param b: 地球短半軸，即大地坐標(biāo)系原點(diǎn)到兩級(jí)的距離, 單位 m:return: B, L, H, 大地緯度、經(jīng)度、海拔高度 (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、經(jīng)度L、海拔高度H (m)'''以下為高斯克呂格投影6度分帶法，經(jīng)度Lon與投影帶號(hào)Zone以及帶號(hào)Zone與投影帶中央子午線經(jīng)度L0的計(jì)算''' def GK(L):'''高斯克呂格投影6度分帶法，經(jīng)度Lon與投影帶號(hào)Zone以及帶號(hào)Zone與投影帶中央子午線經(jīng)度L0的計(jì)算:param L:弧度制:return L0:'''if L>=0:'''東半球時(shí)'''z = int(L/6) + 1L0 = (6 * (z) - 3) * mt.pielse:'''西半球時(shí)'''z = int(L/6) + 60L0 = ((6 * z - 3) - 360)* mt.pireturn L0'''以下為將大地坐標(biāo)轉(zhuǎn)換為高斯平面坐標(biāo)函數(shù)''' def GKxy(B,L,a,b,A1,B1,C1,D1,E1):'''將大地坐標(biāo)轉(zhuǎn)換為高斯平面坐標(biāo):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'''以下定義空間直角坐標(biāo)的存儲(chǔ)列表''' X = [] Y = [] Z = [] '''以下定義大地坐標(biāo)存儲(chǔ)列表''' B = [] L = [] H = [] '''以下定義高斯平面坐標(biāo)存儲(chǔ)列表''' x = [] y = []if __name__ == '__main__':path = input('請(qǐng)輸入空間直角坐標(biāo)源文件：') #C://Users//23646//Downloads//數(shù)據(jù)//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的弧度數(shù)寫(xiě)入到csv文件中'''savepath = 'C://Users//23646//Downloads//數(shù)據(jù)//IGSBLH-rad.csv'data = pd.DataFrame({'B':B,'L':L,'H':H})data.to_csv(savepath,index=False)'''以下為執(zhí)行大地坐標(biāo)轉(zhuǎn)高斯平面坐標(biāo)代碼塊'''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的高斯平面坐標(biāo)寫(xiě)入到csv文件中'''savepath1 = 'C://Users//23646//Downloads//數(shù)據(jù)//IGSxy.csv'data = pd.DataFrame({'x': x, 'y': y})data.to_csv(savepath1, index=False)

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。