User:The Anome/lat-long.py

Lat/long functions in Python
Converts OSGB and Irish grid references to WGS84 latitude/longitude
and formats the resulting lat/long data as a Wikipedia {{coord}} tag

Mostly auto-translated from a original PHP source code by
(various) original authors as listed below, with some manual cleanup.

Python conversion by The Anome, November 2007

You may use and redistribute this code under the terms of the GPL
see http://www.gnu.org/copyleft/gpl.html


ORIGINAL SOURCE COMMENTS BELOW
Lat Long functions in PHP
Note: from http://www.megalithia.com/search/llfuncshighlight.php

With thanks to Andy, G4JNT for inspiration in GEOG, and to the OSGB
for their white paper on coordinate transformation describing the
iterative method used thanks to the Ordnance survey of Ireland for
details of the true and false origins of the Irish grid

You may use and redistribute this code under the terms of the GPL
see http://www.gnu.org/copyleft/gpl.html


Written by Richard
www.megalithia.com
v0.something 27/2/2000
v 1.01 28 June 2004
v 1.02 6 Aug 2004 line 89 add "0" to chars in $ngr = stripcharsnotinbag Thx Andy
v 1.03 9 Mar 2005 Jan (Klingon) added conversion to WGS84 map datum and removed limitation of digits of the grid ref
v 1.04 10 Aug 2005 Richard correct error trapping (only manifest on malformed ngrs

This code is predicated on the assumption that your are ONLY feeding
it Irish or UK Grid references.  It uses the single char prefix of
Irish grid refs to tell the difference, UK grid refs have a two
letter prefix.  We would like an even number of digits for the rest
of the grid ref.  Anything in the NGR other than 0-9, A-Z, a-z is
eliminated.  WARNING this assumes there are no decimal points in
your NGR components.  (before v 1.03) you must have at least 4
numerical digits eg TM1234 and no more than eight eg TM12345678
otherwise you and your data are viewed with suspicion.  NEW (Jan):
at least the letter prefix of the Grid ref, no upper limit

The transformation from OSGB36/Ireland 1965 to WGS84 is more precise than 5 m.

The function NGR2LL(ngr) is case insensitive: its return value is
(country, lat, long) or an exception is raised for any errors


import string, math, re, sys

Get the following maths fns

sin = math.sin
cos = math.cos
tan = math.tan
atan = math.atan
atan2 = math.atan2
floor = math.floor
sqrt = math.sqrt
PI = math.pi
def modified_alpha_offset(ch):
n = ord(ch) - ord("A")
if n > 8: # Letter 'I' is not used
return n-1
else:
return n
def NGR2LL(text):
# Must have at least one letter and one digit, with optional spaces
matches = re.findall(r"^([A-Z]+) *([0-9]+) *([0-9]*)$", string.upper(text))
if not matches:
raise "Malformed NGR"
letters, left, right = matches[0]
# If digits are split,
# check for balanced digit groups
if len(right) != 0 and len(left) != len(right):
raise "Unbalanced digit groups in NGR"
digits = left + right
if len(digits) % 2:
raise "Odd number of digits in NGR"
halflen = len(digits)/2
multiplier = pow(10.0, (5.0-halflen))
n = int(digits[halflen:])*multiplier
e = int(digits[:halflen])*multiplier
if len(letters) == 1:
return IE_NGR2LL(letters, e, n)
if len(letters) == 2:
return GB_NGR2LL(letters, e, n)
raise "Wrong number of letters in NGR"
def IE_NGR2LL(lett, e, n):
# USE IRISH values
# now for the heavy stuff
T1 = modified_alpha_offset(lett[0])
e = 100000.0*(T1%5)+e
n = n+100000.0*(4.0-floor(T1/5.0))
# deg2rad
dr = 180.0/PI
# True Origin is 8 deg W
phi0ir = -8.0
# True Origin is 53.5 deg N
lambda0ir = 53.5
# scale factor @ central meridian
F0ir = 1.000035
# True origin in 200 km E of false origin
E0ir = 200000.0
#True origin is 250km N of false origin
N0ir = 250000.0
# semi-major axis (in line to equator) 1.000035 is yer scale @ central meridian
air = 6377340.189*F0ir
#semi-minor axis (in line to poles)
bir = 6356034.447*F0ir
# flatness = a1-b1/(a1 + b1)
n1ir = 0.001673220384152058651484728058385228837777
# first eccentricity squared = 2*f-f^2 where f = (a1-b1)/a1
e2ir = 0.006670540293336110419293763349975612794125
# radius of earth
#re = 6371.29
k = (n-N0ir)/air+lambda0ir/dr
nextcounter = 0
while 1:
nextcounter = nextcounter+1
k3 = k-lambda0ir/dr
k4 = k+lambda0ir/dr
j3 = (1.0+n1ir+1.25*pow(n1ir, 2.0)+1.25*pow(n1ir, 3.0))*k3
j4 = (3.0*n1ir+3.0*pow(n1ir, 2.0)+2.625*pow(n1ir, 3.0))*sin(k3)*cos(k4)
j5 = (1.875*pow(n1ir, 2.0)+1.875*pow(n1ir, 3.0))*sin(2.0*k3)*cos(2.0*k4)
j6 = 35.0/24.0*pow(n1ir, 3.0)*sin(3.0*k3)*cos(3.0*k4)
m = bir*(j3-j4+j5-j6)
k = k+(n-N0ir-m)/air
if not ((abs(n-N0ir-m)>0.000000000001) and (nextcounter<10000)):
break
v = air/sqrt(1.0-e2ir*pow(sin(k), 2.0))
r = v*(1.0-e2ir)/(1.0-e2ir*pow(sin(k), 2.0))
h2 = v/r-1.0
y1 = e-E0ir
j3 = tan(k)/(2.0*r*v)
j4 = tan(k)/(24.0*r*pow(v, 3.0))*(5.0+3.0*pow(tan(k), 2.0)+h2-9.0*pow(tan(k), 2.0)*h2)
j5 = tan(k)/(720.0*r*pow(v, 5.0))*(61.0+90.0*pow(tan(k), 2.0)+45.0*pow(tan(k), 4.0))
k9 = k-y1*y1*j3+pow(y1, 4.0)*j4-pow(y1, 6.0)*j5
j6 = 1.0/(cos(k)*v)
j7 = 1.0/(cos(k)*6.0*pow(v, 3.0))*(v/r+2.0*pow(tan(k), 2.0))
j8 = 1.0/(cos(k)*120.0*pow(v, 5.0))*(5.0+28.0*pow(tan(k), 2.0)+24.0*pow(tan(k), 4.0))
j9 = 1.0/(cos(k)*5040.0*pow(v, 7.0))
j9 = j9*(61.0+662.0*pow(tan(k), 2.0)+1320.0*pow(tan(k), 4.0)+720.0*pow(tan(k), 6.0))
long = phi0ir+dr*(y1*j6-y1*y1*y1*j7+pow(y1, 5.0)*j8-pow(y1, 7.0)*j9)
lat = k9*dr
# v1.04 this bracket moved to just before elsif # }
# convert long/lat to Cartesian coordinates
v = 6377340.189/sqrt(1.0-e2ir*pow(sin(k), 2.0))
cartxa = v*cos(k9)*cos(long/dr)
cartya = v*cos(k9)*sin(long/dr)
cartza = (1.0-e2ir)*v*sin(k9)
# Helmert-Transformation from Ireland 1965 to WGS84 map date
rotx = 1.042/3600.0*PI/180.0
roty = 0.214/3600.0*PI/180.0
rotz = 0.631/3600.0*PI/180.0
scale = 8.15/1000000.0
cartxb = 482.53+(1.0+scale)*cartxa+rotz*cartya-roty*cartza
cartyb = -130.596-rotz*cartxa+(1.0+scale)*cartya+rotx*cartza
cartzb = 564.557+roty*cartxa-rotx*cartya+(1.0+scale)*cartza
# convert Cartesian to long/lat
awgs84 = 6378137.0
bwgs84 = 6356752.3141
e2wgs84 = 0.00669438003551279089034150031998869922791
lambdaradwgs84 = atan(cartyb/cartxb)
long = lambdaradwgs84*180.0/PI
pxy = sqrt(pow(cartxb, 2.0)+pow(cartyb, 2.0))
phiradwgs84 = atan(cartzb/pxy/(1.0-e2wgs84))
nextcounter = 0
while 1:
nextcounter = nextcounter+1
v = awgs84/sqrt(1.0-e2wgs84*pow(sin(phiradwgs84), 2.0))
phinewwgs84 = atan((cartzb+e2wgs84*v*sin(phiradwgs84))/pxy)
# Now testing _before_ the assignment below...
if abs(phinewwgs84-phiradwgs84)<0.000000000001:
break
if nextcounter > 20:
raise "loop not converging"
phiradwgs84 = phinewwgs84
lat = phiradwgs84*180.0/PI
return ("IE", lat, long)
def GB_NGR2LL(lett, e, n):
# Use British values
# now for the heavy stuff
T1 = modified_alpha_offset(lett[0])
T2 = modified_alpha_offset(lett[1])
e = 500000.0*(T1%5)+100000.0*(T2%5)-1000000.0+e
n = 1900000.0-500000.0*floor(T1/5.0)-100000.0*floor(T2/5.0)+n
# deg2rad
dr = 180.0/PI
# True Origin is 2 deg W
phi0uk = -2.0
# True Origin is 49 deg N
lambda0uk = 49.0
# scale factor @ central meridian
F0uk = 0.9996012717
# True origin in 400 km E of false origin
E0uk = 400000.0
#True origin is 100 km S of false origin
N0uk = -100000.0
# semi-major axis (in line to equator) 0.996012717 is yer scale @ central meridian
auk = 6377563.396*F0uk
#semi-minor axis (in line to poles)
buk = 6356256.91*F0uk
# flatness = a1-b1/(a1+b1)
n1uk = 0.00167322025032508731869331280635710896296
# first eccentricity squared = 2*f-f^2where f = (a1-b1)/a1
e2uk = 0.006670539761597529073698869358812557054558
# radius of earth
#re = 6371.29
k = (n-N0uk)/auk+lambda0uk/dr
nextcounter = 0
while 1:
nextcounter = nextcounter+1
k3 = k-lambda0uk/dr
k4 = k+lambda0uk/dr
j3 = (1.0+n1uk+1.25*pow(n1uk, 2.0)+1.25*pow(n1uk, 3.0))*k3
j4 = (3.0*n1uk+3.0*pow(n1uk, 2.0)+2.625*pow(n1uk, 3.0))*sin(k3)*cos(k4)
j5 = (1.875*pow(n1uk, 2.0)+1.875*pow(n1uk, 3.0))*sin(2.0*k3)*cos(2.0*k4)
j6 = 35.0/24.0*pow(n1uk, 3.0)*sin(3.0*k3)*cos(3.0*k4)
m = buk*(j3-j4+j5-j6)
k = k+(n-N0uk-m)/auk
if not ((abs(n-N0uk-m)>0.000000000001) and (nextcounter<10000)):
break
v = auk/sqrt(1.0-e2uk*pow(sin(k), 2.0))
r = v*(1.0-e2uk)/(1.0-e2uk*pow(sin(k), 2.0))
h2 = v/r-1.0
y1 = e-E0uk
j3 = tan(k)/(2.0*r*v)
j4 = tan(k)/(24.0*r*pow(v, 3.0))*(5.0+3.0*pow(tan(k), 2.0)+h2-9.0*pow(tan(k), 2.0)*h2)
j5 = tan(k)/(720.0*r*pow(v, 5.0))*(61.0+90.0*pow(tan(k), 2.0)+45.0*pow(tan(k), 4.0))
k9 = k-y1*y1*j3+pow(y1, 4.0)*j4-pow(y1, 6.0)*j5
j6 = 1.0/(cos(k)*v)
j7 = 1.0/(cos(k)*6.0*pow(v, 3.0))*(v/r+2.0*pow(tan(k), 2.0))
j8 = 1.0/(cos(k)*120.0*pow(v, 5.0))*(5.0+28.0*pow(tan(k), 2.0)+24.0*pow(tan(k), 4.0))
j9 = 1.0/(cos(k)*5040.0*pow(v, 7.0))
j9 = j9*(61.0+662.0*pow(tan(k), 2.0)+1320.0*pow(tan(k), 4.0)+720.0*pow(tan(k), 6.0))
long = phi0uk+dr*(y1*j6-y1*y1*y1*j7+pow(y1, 5.0)*j8-pow(y1, 7.0)*j9)
lat = k9*dr
# v1.04 this bracket moved to just before elsif # }
# convert long/lat to Cartesian coordinates
v = 6377563.396/sqrt(1.0-e2uk*pow(sin(k), 2.0))
cartxa = v*cos(k9)*cos(long/dr)
cartya = v*cos(k9)*sin(long/dr)
cartza = (1.0-e2uk)*v*sin(k9)
# Helmert-Transformation from OSGB36 to WGS84 map date
rotx = -0.1502/3600.0*PI/180.0
roty = -0.2470/3600.0*PI/180.0
rotz = -0.8421/3600.0*PI/180.0
scale = -20.4894/1000000.0
cartxb = 446.448+(1.0+scale)*cartxa+rotz*cartya-roty*cartza
cartyb = -125.157-rotz*cartxa+(1.0+scale)*cartya+rotx*cartza
cartzb = 542.06+roty*cartxa-rotx*cartya+(1.0+scale)*cartza
# convert Cartesian to long/lat
awgs84 = 6378137.0
bwgs84 = 6356752.3141
e2wgs84 = 0.00669438003551279089034150031998869922791
lambdaradwgs84 = atan(cartyb/cartxb)
long = lambdaradwgs84*180.0/PI
pxy = sqrt(pow(cartxb, 2.0)+pow(cartyb, 2.0))
phiradwgs84 = atan(cartzb/pxy/(1.0-e2wgs84))
nextcounter = 0
while 1:
nextcounter = nextcounter+1
v = awgs84/sqrt(1.0-e2wgs84*pow(sin(phiradwgs84), 2.0))
phinewwgs84 = atan((cartzb+e2wgs84*v*sin(phiradwgs84))/pxy)
# Now testing _before_ the assignment below...
if abs(phinewwgs84-phiradwgs84)<0.000000000001:
break
if nextcounter > 20:
raise "loop not converging"
phiradwgs84 = phinewwgs84
lat = phiradwgs84*180.0/PI
# v 1.04 mod
return ("GB", lat, long)
def dirselect(x, options):
if x < 0:
return options[1]
else:
return options[0]
osref = string.join(sys.argv[1:])
country, lat, long = NGR2LL(osref)
print (
"{{coor title d|%1.5lf|%s|%1.5lf|%s|region:%s_source:enwiki-osgb36(%s)}}" %
(abs(lat), dirselect(lat, "NS"),
abs(long), dirselect(long, "EW"),
country, osref))