ITM（爱尔兰横坐标）转换为谷歌地图 Python3 的 GPS

Question

我对坐标一无所知。我的问题是我有一个包含 ITM 格式（Irish_X 和 Irish_Y）坐标的数据集。 我想将 ITM 坐标转换为 Google 地图可读的坐标。

在线我发现了一个可能有用但我不知道如何使用的库，并且文档使用了我不习惯的非常具体的行话： https://proj.org

我还在同一个库 gitHub 存储库中评论了寻找答案： https://github.com/OSGeo/PROJ/issues/1687

非常感谢您的帮助！

Answer 1

要从 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])