User:The Anome/lat-long.py


 * 1) Lat/long functions in Python
 * 2) Converts OSGB and Irish grid references to WGS84 latitude/longitude
 * 3) and formats the resulting lat/long data as a Wikipedia °N, °W tag
 * 4) Mostly auto-translated from a original PHP source code by
 * 5) (various) original authors as listed below, with some manual cleanup.
 * 6) Python conversion by The Anome, November 2007
 * 7) You may use and redistribute this code under the terms of the GPL
 * 8) see http://www.gnu.org/copyleft/gpl.html
 * 1) You may use and redistribute this code under the terms of the GPL
 * 2) see http://www.gnu.org/copyleft/gpl.html
 * 1) see http://www.gnu.org/copyleft/gpl.html


 * 1) ORIGINAL SOURCE COMMENTS BELOW
 * 2) Lat Long functions in PHP
 * 3) Note: from http://www.megalithia.com/search/llfuncshighlight.php
 * 4) With thanks to Andy, G4JNT for inspiration in GEOG, and to the OSGB
 * 5) for their white paper on coordinate transformation describing the
 * 6) iterative method used thanks to the Ordnance survey of Ireland for
 * 7) details of the true and false origins of the Irish grid
 * 8) You may use and redistribute this code under the terms of the GPL
 * 9) see http://www.gnu.org/copyleft/gpl.html
 * 10) Written by Richard
 * 11) www.megalithia.com
 * 12) v0.something 27/2/2000
 * 13) v 1.01 28 June 2004
 * 14) v 1.02 6 Aug 2004 line 89 add "0" to chars in $ngr = stripcharsnotinbag Thx Andy
 * 15) v 1.03 9 Mar 2005 Jan (Klingon) added conversion to WGS84 map datum and removed limitation of digits of the grid ref
 * 16) v 1.04 10 Aug 2005 Richard correct error trapping (only manifest on malformed ngrs
 * 17) This code is predicated on the assumption that your are ONLY feeding
 * 18) it Irish or UK Grid references.  It uses the single char prefix of
 * 19) Irish grid refs to tell the difference, UK grid refs have a two
 * 20) letter prefix.  We would like an even number of digits for the rest
 * 21) of the grid ref.  Anything in the NGR other than 0-9, A-Z, a-z is
 * 22) eliminated.  WARNING this assumes there are no decimal points in
 * 23) your NGR components.  (before v 1.03) you must have at least 4
 * 24) numerical digits eg TM1234 and no more than eight eg TM12345678
 * 25) otherwise you and your data are viewed with suspicion.  NEW (Jan):
 * 26) at least the letter prefix of the Grid ref, no upper limit
 * 27) The transformation from OSGB36/Ireland 1965 to WGS84 is more precise than 5 m.
 * 28) The function NGR2LL(ngr) is case insensitive: its return value is
 * 29) (country, lat, long) or an exception is raised for any errors
 * 1) otherwise you and your data are viewed with suspicion.  NEW (Jan):
 * 2) at least the letter prefix of the Grid ref, no upper limit
 * 3) The transformation from OSGB36/Ireland 1965 to WGS84 is more precise than 5 m.
 * 4) The function NGR2LL(ngr) is case insensitive: its return value is
 * 5) (country, lat, long) or an exception is raised for any errors
 * 1) The function NGR2LL(ngr) is case insensitive: its return value is
 * 2) (country, lat, long) or an exception is raised for any errors

import string, math, re, sys

sin = math.sin cos = math.cos tan = math.tan atan = math.atan atan2 = math.atan2 floor = math.floor sqrt = math.sqrt PI = math.pi
 * 1) Get the following maths fns

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 ( "" %         (abs(lat), dirselect(lat, "NS"), abs(long), dirselect(long, "EW"), country, osref))