Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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
share|improve this question
    
I forgot to say that I checked the results against the Ordnance Survey Ireland on-line converter at osi.ie/calculators/converter_index.asp?alias=/services/… –  aacinelli Aug 29 '12 at 11:37
    
How does your output compare to the worked example on p. 36? –  whuber Aug 29 '12 at 12:09
3  
Not 100% on topic, but have you tried using the pyproj library to compare your results also? –  om_henners Aug 29 '12 at 12:35
2  
@aacinelli When porting fiddly numeric code you have to reproduce any worked examples you can find: it's not enough just to get the right answers in a few test cases. A single typographical error can cause the results to be very wrong for just a small subset of possible inputs and only a tiny bit wrong for all other inputs, making it hard to detect, test, and debug: that's why the example is there showing lots of intermediate results. –  whuber Aug 29 '12 at 13:06
1  
I found the bug in the script and it has been amended. When calculating variables M1, M2, M3 and M4 the script now uses floating numbers divisions and results obtained match the ones from Ordnance Survey Ireland on-line converter and the pyproj Python lib. Thanks everyone for your comments and help. Andrea –  aacinelli Sep 1 '12 at 15:09
show 4 more comments

1 Answer 1

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)
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.