google 地图 Python3 的 ITM(爱尔兰横向坐标)转换为 GPS
ITM (Irish Transverse Coordinate) conversion to GPS for google maps Python3
我对坐标一无所知。我的问题是我有一个包含 ITM 格式坐标的数据集(Irish_X 和 Irish_Y)。 我想将 ITM 坐标转换为 Google 地图可读坐标。
在网上我找到了一个可能有用但我不知道如何使用的库,并且文档使用了我不习惯的非常具体的行话:
https://proj.org
我也在同一个库中评论 gitHub 回购寻找答案:
https://github.com/OSGeo/PROJ/issues/1687
非常感谢您的帮助!
要从 ITM 转换为 WGS84,只需调用 def itm2geo(x,y):
函数,如 Teste Values
部分代码末尾所示。
前两个函数(arcmer
和xy2geo
)是辅助函数,不需要显式调用(它们被itm2geo(x,y)
调用
from math import *
############################################################################
# Meridian Arc
############################################################################
def arcmer(a,equad,lat1,lat2):
b=a*sqrt(1-equad)
n=(a-b)/(a+b)
a0=1.+((n**2)/4.)+((n**4)/64.)
a2=(3./2.)*(n-((n**3)/8.))
a4=(15./16.)*((n**2)-((n**4)/4.))
a6=(35./48.)*(n**3)
s1=a/(1+n)*(a0*lat1-a2*sin(2.*lat1)+a4*sin(4.*lat1)-a6*sin(6.*lat1))
s2=a/(1+n)*(a0*lat2-a2*sin(2.*lat2)+a4*sin(4.*lat2)-a6*sin(6.*lat2))
return s2-s1
###############################################################################
#
# Transverse Mercator Inverse Projection
#
###############################################################################
def xy2geo(m,p,a,equad,lat0,lon0):
lat0=radians(lat0)
lon0=radians(lon0)
sigma1=p
fil=lat0+sigma1/(a*(1-equad))
deltafi=1
while deltafi > 0.0000000001:
sigma2=arcmer(a,equad,lat0,fil)
RO=a*(1-equad)/((1-equad*(sin(fil)**2))**(3./2.))
deltafi=(sigma1-sigma2)/RO
fil=fil+deltafi
N=a/sqrt(1-equad*(sin(fil))**2)
RO=a*(1-equad)/((1-equad*(sin(fil)**2))**(3./2.))
t=tan(fil)
psi=N/RO
lat=fil-(t/RO)*((m**2)/(2.*N))+(t/RO)*((m**4)/(24.*(N**3)))*(-4.*(psi**2)-9.*psi*(1.-t**2)+12.*(t**2))-(t/RO)*(m**6/(720.*(N**5)))*(8.*(psi**4)*(11.-24.*(t**2))-12.*(psi**3)*(21.-71.*(t**2))+15.*(psi**2)*(15.-98.*(t**2)+15.*(t**4))+180.*psi*(5.*(t**2)-3.*(t**4))-360.*(t**4))+(t/RO)*((m**8)/(40320.*(N**7)))*(1385.+3633.*(t**2)+4095.*(t**4)+1575.*(t**6))
lon=(m/(N))-((m**3)/(6.*(N**3)))*(psi+2.*(t**2))+((m**5)/(120.*(N**5)))*(-4.*(psi**3)*(1.-6.*(t**2))+(psi**2)*(9.-68.*(t**2))+72.*psi*(t**2)+24.*(t**4))-((m**7)/(5040.*(N**7)))*(61.+662.*(t**2)+1320.*(t**4)+720.*(t**6))
lon=lon0+lon/cos(fil)
lat=degrees(lat)
lon=degrees(lon)
return lat,lon
#############################################################################
# Irish Transverse Mercator - Inverse
#############################################################################
def itm2geo(x,y):
# GRS-80
a = 6378137.
equad =0.00669437999
# Natural Origin
lat0=53.5
lon0=-8.
k0=0.999820
p = (y - 750000.) /k0
m = (x - 600000.) /k0
lat,lon = xy2geo(m,p,a,equad,lat0,lon0)
return lat,lon
#############################################################################
# Test values
#############################################################################
#lat=53.5
#lon=-8.
test = itm2geo(600000.,750000.)
print ("latitude= %.16f" %test[0])
print ("longitude= %.16f" %test[1])
我对坐标一无所知。我的问题是我有一个包含 ITM 格式坐标的数据集(Irish_X 和 Irish_Y)。 我想将 ITM 坐标转换为 Google 地图可读坐标。
在网上我找到了一个可能有用但我不知道如何使用的库,并且文档使用了我不习惯的非常具体的行话: https://proj.org
我也在同一个库中评论 gitHub 回购寻找答案: https://github.com/OSGeo/PROJ/issues/1687
非常感谢您的帮助!
要从 ITM 转换为 WGS84,只需调用 def itm2geo(x,y):
函数,如 Teste Values
部分代码末尾所示。
前两个函数(arcmer
和xy2geo
)是辅助函数,不需要显式调用(它们被itm2geo(x,y)
from math import *
############################################################################
# Meridian Arc
############################################################################
def arcmer(a,equad,lat1,lat2):
b=a*sqrt(1-equad)
n=(a-b)/(a+b)
a0=1.+((n**2)/4.)+((n**4)/64.)
a2=(3./2.)*(n-((n**3)/8.))
a4=(15./16.)*((n**2)-((n**4)/4.))
a6=(35./48.)*(n**3)
s1=a/(1+n)*(a0*lat1-a2*sin(2.*lat1)+a4*sin(4.*lat1)-a6*sin(6.*lat1))
s2=a/(1+n)*(a0*lat2-a2*sin(2.*lat2)+a4*sin(4.*lat2)-a6*sin(6.*lat2))
return s2-s1
###############################################################################
#
# Transverse Mercator Inverse Projection
#
###############################################################################
def xy2geo(m,p,a,equad,lat0,lon0):
lat0=radians(lat0)
lon0=radians(lon0)
sigma1=p
fil=lat0+sigma1/(a*(1-equad))
deltafi=1
while deltafi > 0.0000000001:
sigma2=arcmer(a,equad,lat0,fil)
RO=a*(1-equad)/((1-equad*(sin(fil)**2))**(3./2.))
deltafi=(sigma1-sigma2)/RO
fil=fil+deltafi
N=a/sqrt(1-equad*(sin(fil))**2)
RO=a*(1-equad)/((1-equad*(sin(fil)**2))**(3./2.))
t=tan(fil)
psi=N/RO
lat=fil-(t/RO)*((m**2)/(2.*N))+(t/RO)*((m**4)/(24.*(N**3)))*(-4.*(psi**2)-9.*psi*(1.-t**2)+12.*(t**2))-(t/RO)*(m**6/(720.*(N**5)))*(8.*(psi**4)*(11.-24.*(t**2))-12.*(psi**3)*(21.-71.*(t**2))+15.*(psi**2)*(15.-98.*(t**2)+15.*(t**4))+180.*psi*(5.*(t**2)-3.*(t**4))-360.*(t**4))+(t/RO)*((m**8)/(40320.*(N**7)))*(1385.+3633.*(t**2)+4095.*(t**4)+1575.*(t**6))
lon=(m/(N))-((m**3)/(6.*(N**3)))*(psi+2.*(t**2))+((m**5)/(120.*(N**5)))*(-4.*(psi**3)*(1.-6.*(t**2))+(psi**2)*(9.-68.*(t**2))+72.*psi*(t**2)+24.*(t**4))-((m**7)/(5040.*(N**7)))*(61.+662.*(t**2)+1320.*(t**4)+720.*(t**6))
lon=lon0+lon/cos(fil)
lat=degrees(lat)
lon=degrees(lon)
return lat,lon
#############################################################################
# Irish Transverse Mercator - Inverse
#############################################################################
def itm2geo(x,y):
# GRS-80
a = 6378137.
equad =0.00669437999
# Natural Origin
lat0=53.5
lon0=-8.
k0=0.999820
p = (y - 750000.) /k0
m = (x - 600000.) /k0
lat,lon = xy2geo(m,p,a,equad,lat0,lon0)
return lat,lon
#############################################################################
# Test values
#############################################################################
#lat=53.5
#lon=-8.
test = itm2geo(600000.,750000.)
print ("latitude= %.16f" %test[0])
print ("longitude= %.16f" %test[1])