[GIS] Python script to convert Lat Long to ITM (Irish Transverse Mercator)

coordinate systempythonwgs84

I am currently working on a python script that converts Lat Long (WGS84) to ITM. The formulas used are based on a paper published by OSGB (Appendix B):

http://www.leica-geosystems.co.uk/downloads123/gb/gps/gps1200/other/OS%20Transformations%20and%20OSGM02%20user%20guide_en.pdf.

I get good results in the Eastings but not in the Northings (from few mm to approx. 1.5m) and I can't understand why. Does anyone know the reasons?
Thank for your help.

Andrea

# this script convert Latitude and Longitude to ITM

# import modules
import math

# GRS80 Ellipsoid constants
a = 6378137.0
b = 6356752.3141
e2 = (a**2-b**2)/a**2

# ITM projection constants
F = 0.99982
Lat0 = 0.933751149817
Long0 = -0.13962634016
E0 = 600000.0
N0 = 750000.0

# input Lat degrees, minutes and seconds for latitude
dLat = float(raw_input('Enter latitude degrees: '))
mLat = float(raw_input('Enter latitude minutes: '))
sLat = float(raw_input('Enter latitude seconds: '))

# convert Lat to decimal degrees
ddLat = dLat + (mLat/60) + (sLat/3600)
print 'Latitude dd:', ddLat

# convert Lat to radians
rLat = (ddLat*math.pi)/180.0
Lat = rLat
print 'Latitude r:', Lat
print

# input Long degrees, minutes and seconds for longitude
dLong = float(raw_input('Enter longitude degrees: '))  #The degrees have to be entered with a posite sign
mLong = float(raw_input('Enter longitude minutes: '))
sLong = float(raw_input('Enter longitude seconds: '))

# convert Long to decimal degrees
ddLong = dLong + (mLong/60) + (sLong/3600)
print 'Longitude dd:', ddLong

# convert Long to radians
rLong = -(ddLong*math.pi)/180.0 # the negative sign account for west direction
Long = rLong
print 'Longitude r:', Long
print

# calculate constants for transformation
n = (a - b)/(a + b)
v = a*F*(1-e2*(math.sin(Lat)**2))**-0.5
p = a*F*(1-e2)*(1-e2*(math.sin(Lat)**2))**-1.5
n2 = v/p-1
print 'n:', n
print 'v:', v
print 'p:', p
print 'n2:', n2
print

M1 = (1+n+5.0/4.0*n**2+5.0/4.0*n**3)*(Lat-Lat0)
M2 = (3*n+3*n**2+21.0/8.0*n**3)*math.sin(Lat-Lat0)*math.cos(Lat+Lat0)
M3 = (15.0/8.0*n**2+15.0/8.0*n**3)*math.sin(2*(Lat-Lat0))*math.cos(2*(Lat+Lat0))
M4 = 35.0/24.0*n**3*math.sin(3*(Lat-Lat0))*math.cos(3*(Lat+Lat0))
M = b*F*(M1-M2+M3-M4)
print 'M1:', M1
print 'M2:', M2
print 'M3:', M3
print 'M4:', M4
print 'M:', M
print

I = M+N0
print 'I:', I
II = v/2*math.sin(Lat)*math.cos(Lat)
print 'II:', II
III = v/24*math.sin(Lat)*(math.cos(Lat))**3*(5-(math.tan(Lat)**2)+9*n2)
print 'III:', III
IIIA = v/720*math.sin(Lat)*(math.cos(Lat)**5)*(61-58*(math.tan(Lat)**2)+(math.tan(Lat)**4))
print 'IIIA:', IIIA
IV = v*math.cos(Lat)
print 'IV:', IV
V = v/6*(math.cos(Lat)**3)*(v/p-(math.tan(Lat)**2))
print 'V:', V
VI = v/120*(math.cos(Lat)**5)*(5-18*(math.tan(Lat)**2)+(math.tan(Lat)**4)+14*n2-58*(math.tan(Lat)**2)*n2)
print 'VI:', VI
print

# calculate Eastings and Northings
N = I+II*(Long-Long0)**2+III*(Long-Long0)**4+IIIA*(Long-Long0)**6
E = E0+IV*(Long-Long0)+V*(Long-Long0)**3+VI*(Long-Long0)**5
print 'Easting: ', E
print 'Northing:', N

Best Answer

As om_henners advised it is better to use available library for this purpose as it is already implemented and tested by many people...

So, take a look at pyproj Python lib.

Here is a sample code for reprojecting WGS-84 long/lat to ITM (EPSG:2157) x,y:

from pyproj import Proj, transform


def reproject_wgs_to_itm(longitude, latitude):
    prj_wgs = Proj(init='epsg:4326')
    prj_itm = Proj(init='epsg:2157')
    x, y = transform(prj_wgs, prj_itm, longitude, latitude)
    return x, y


print reproject_wgs_to_itm(-7.748108, 53.431628)
Related Question