ITM (Irish Transverse Coordinate) conversion to GPS for google maps Python3

Question

I don't know anything about coordinates. My problem is that I have a dataset that contains coordinates in the ITM format (Irish_X and Irish_Y). I want to convert ITM coordinates into Google Maps readable ones.

Online I found a library that may be useful but I don't know how to use and the documentation use a very specific jargon I'm not used to: https://proj.org

I also commented in the same library gitHub repo looking for answer: https://github.com/OSGeo/PROJ/issues/1687

Thank you very much for your help!

Answer 1

To convert from ITM to WGS84 just call the def itm2geo(x,y): function as exemplified in the end of the code in the Teste Values section.

The first two functions (arcmer and xy2geo) are auxiliary functions and don't need to be explicitly called (they are called by 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])