Groovy Script demonstrating how to convert from OSGB36 Eastings & Northings to Latitude/Longitude
I recently commented in my Contacts Domain Model post that Ordnance survey has released a free database of UK postal codes with precision down to the individual postcode, which is better than the other free alternative you’ve probably heard about from Nearby.org, which doesn’t add the level of precision down to the last two characters of a UK post code. This is the current database employed in the Grails Postcode plug-in.
Unfortunately the database is provided with data in the format of Eastings/Northings, rather than latitude & longitude.
After getting past the jargon and terminology and maths overload, I was able to conjure up a Groovy script to do the conversion, which I’ll later apply to the whole database, which consists of one CSV file per postcode area. There are 120 such CSV files in total in the download. As I mentioned earlier the post codes are in column A and the Eastings & Northings are in columns K & L respectively.
I’ve taken the Annexe section from A guide to coordinate systems in Great Britain and restructured it into a table that helped me summarize and get my head around the jargon like OSGB36 etc and shown the Ellipsoid & UK Transverse Mercator projections. These are defined as maps named airy1830, airy1830mod, i1924_ed50_utm30_31 & utm29 in my Groovy script below. If you look at the bottom of my script you will see the convert method gets invoked with the airy1830 version defined at the top of the script.
In essence you use the airy1830modified for Ireland and airy1830 for the UK. Other countries use different constants to represent the 3D shape of the earth.
Ellipsoid Airy 1830 Airy 1830 modified International 1924 aka Hayford 1909 GRS80 aka WGS84 ellipsoid Associated datums and projections OSGB36, National Grid Ireland 65, Irish National Grid ED50, UTM WGS84, ITRS, ETRS89 semi-major axis, a 6 377 563.396 6 377 340.189 6 378 388.000 6 378 137.000 semi-minor axis, b 6 356 256.910 6 356 034.447 6 356 911.946 6 356 752.3141 UK Transverse Mercator projections: Central Meridian Scale F0 0.9996012717 1.000035 0.9996 ? True origin φ0 lat 49° N lat 53°30′ N lat 0° True origin λ0 long 2° W long 8° W long 3° W (UTM Zone 30/31) long 9° W (UTM zone 29) Map coordinates of true origin (metres), E0 and E 400 000 E 200 000 E 500 000 N0 N -100 000 N 250 000 N0
Here is the script to do the conversion (en_to_ll.groovy):
package jgf import java.text.DecimalFormat // Ellipsoid/projection constants table A plus Transverse Mercator projections used in UK // Taken from P40 of Ordnance Survey : A guide to coordinate systems in Great Britain // E is positive, W negative. Since UK in Northern hemisphere lat always plus for UK Post codes..... // Airy 1830 for UK Mainland. Airy1830mod for Ireland... def airy1830 = [a: 6377563.396, b:6356256.910, f0:0.9996012717, lat0deg:49, long0deg:-2, e0:400000, n0:-100000] def airy1830mod = [a: 6377340.189, b:6356034.447, f0:1.000035, lat0deg:53.5, long0deg:-8, e0:200000, n0:250000] def i1924_ed50_utm30_31 = [a: 6378388.000, b:6356911.946, f0:0.9996, lat0deg:0, long0deg:-3, e0:500000, n0:0] def utm29 = [a: 6378388.000, b:6356911.946, f0:0.9996, lat0deg:0, long0deg:-9, e0:500000, n0:0] class east_north_to_lat_long { // Co-ordinates used in Code-Point CSV database. Eastings & Northings to be converted to latitude & longitude.... double easting double northing def convert(east, north, proj) { // Number formats for various intermediary results def sc10dp = new DecimalFormat("0.#########E0") def df7_8 = new DecimalFormat("#,###,###.########") def df6_3 = new DecimalFormat("###,###.###") def df1_9 = new DecimalFormat("0.#########") def df1_10 = new DecimalFormat("0.##########") def pr = proj // Always working in radians .... pr.lat0 = Math.toRadians(pr.lat0deg) pr.long0 = Math.toRadians(pr.long0deg) double a_p2 = Math.pow(pr.a, 2) double e_sq = (a_p2 - Math.pow(pr.b, 2)) / a_p2 double one_m_e_sq = 1.0 - e_sq double tol = 0.01 // 0.00001 double diff = 0.02 double n = (pr.a -pr.b)/(pr.a +pr.b) double n_p2 = Math.pow(n,2) double n_p3 = Math.pow(n,3) double af0 = (pr.a * pr.f0) double bf0 = (pr.b * pr.f0) def e_e0s = [] double e_e0 = east - pr.e0 (1..7).each{e_e0s << Math.pow(e_e0, it)} double n_n0 = north - pr.n0 double m_old = 0.0 double lat def latseed = pr.lat0 def latresults = [] while (diff > tol) { double n_n0_m_old = n_n0 - m_old double n_n0_m_old_o_af0 = n_n0_m_old / af0 lat = n_n0_m_old_o_af0 + latseed // In radians double latplat0 = lat + pr.lat0 double latmlat0 = lat - pr.lat0 double m_1 = (1 + n + ((5/4) * n_p2) + ((5/4) * n_p3)) * latmlat0 double m_2 = ((3 * n) +(3 * n_p2) +((21/8) * n_p3)) * Math.sin(latmlat0) * Math.cos(latplat0) double m_3 = ( ((15/8) * n_p2) + ((15/8) * n_p3) ) * Math.sin(2 *latmlat0) * Math.cos(2 * latplat0) double m_4 = ( (35/24) * n_p3) * Math.sin(3 *latmlat0) * Math.cos(3 * latplat0) double m_new = bf0 * (m_1 - m_2 + m_3 - m_4) double n_n0_m_new = n_n0 - m_new double abs_n_n0_m_new = Math.abs(n_n0_m_new) diff = abs_n_n0_m_new def latres = [:] latres.m_old = m_old latres.latseed = latseed latres.lat = lat latres.latplat0 = latplat0 latres.latmlat0 = latmlat0 latres.m_1 = m_1 latres.m_2 = m_2 latres.m_3 = m_3 latres.m_4 = m_4 latres.m_new = m_new latres.n_n0_m_new = n_n0_m_new latres.abs_n_n0_m_new = abs_n_n0_m_new latresults << latres latseed = lat m_old = m_new } lat = latresults[latresults.size() -1].lat double sin_sq_lat = Math.pow(Math.sin(lat), 2) double one_m_e_sq_sin_sq_lat = 1 - (e_sq * Math.pow(Math.sin(lat), 2)) double v = af0 * Math.pow(one_m_e_sq_sin_sq_lat,-0.5) double p = af0 * one_m_e_sq * Math.pow(one_m_e_sq_sin_sq_lat,-1.5) double n_sq = (v/p) -1 double tan_lat = Math.tan(lat) double tan_lat_p2 = Math.pow(tan_lat, 2) double tan_lat_p4 = Math.pow(tan_lat, 4) double tan_lat_p6 = Math.pow(tan_lat, 6) double sec_lat = 1 / (Math.cos(lat)) double vii = tan_lat / (2 * p * v) double viii = (tan_lat / (24 * p * Math.pow(v,3)) ) * (5 + (3 * tan_lat_p2) + n_sq - (9 * tan_lat_p2 * n_sq)) double ix = (tan_lat/(720 * p * Math.pow(v,5))) * (61 + (90 * tan_lat_p2) + 45 * tan_lat_p4) double x = sec_lat / v double xi = (sec_lat / (6 * Math.pow(v,3))) * ((v/p) + 2 * tan_lat_p2) double xii = (sec_lat/(120 * Math.pow(v,5))) * (5 + (28 * tan_lat_p2) + (24 * tan_lat_p4)) double xiia = (sec_lat/(5040 * Math.pow(v,7))) * (61 + (662 * tan_lat_p2) + (1320 * tan_lat_p4) + (720 * tan_lat_p6)) double latres = lat - (vii * e_e0s[1]) + (viii * e_e0s[3]) - (ix * e_e0s[5]) // =_lat1.2-(_vii*_E_E0_p2)+(_viii*_E_E0_p4)-(_ix*_E_E0_p6) double longres = pr.long0 + (x * e_e0s[0]) - (xi * e_e0s[2]) + (xii * e_e0s[4]) - (xiia * e_e0s[6]) // =_long0+(_x*_E_E0)-(_xi*_E_E0_p3)+(_xii*_E_E0_p5)-(_x11a*_E_E0_p7) double latresdegdec = Math.toDegrees(latres) double longresdegdec = Math.toDegrees(longres) def latresdms = converDegreesDecToDMS(latresdegdec) def longresdms = converDegreesDecToDMS(longresdegdec) // Dump out intermediary calculations... println "tol : ${df1_9.format(tol)}" println "pr.lat0 : ${df1_9.format(pr.lat0)}" println "pr.long0 : ${df1_9.format(pr.long0)}" println "one_m_e_sq : ${df1_10.format(one_m_e_sq)}" println "e_sq : ${sc10dp.format(e_sq)}" println "n : ${sc10dp.format(n)}" println "n_p2 : ${sc10dp.format(n_p2)}" println "n_p3 : ${sc10dp.format(n_p3)}" println "af0 : ${df7_8.format(af0)}" println "bf0 : ${df7_8.format(bf0)}" e_e0s.eachWithIndex{e0, ei -> print "e_e0^${ei+1} : ${(ei) ? sc10dp.format(e0) : df6_3.format(e0)} " } println " " println "n_n0 : ${df6_3.format(n_n0)}" println "---" latresults.eachWithIndex{it, i -> print "${i+1}: m_old: ${sc10dp.format(it.m_old)} " print "latseed: ${sc10dp.format(it.latseed)} " print "lat: ${sc10dp.format(it.lat)} " print "latplat0: ${df1_9.format(it.latplat0)} " print "latmlat0: ${df1_9.format(it.latmlat0)} " print "m_1: ${sc10dp.format(it.m_1)} m_2: ${sc10dp.format(it.m_2)} m_3: ${sc10dp.format(it.m_3)} m_4: ${sc10dp.format(it.m_4)} " print "m_new: ${sc10dp.format(it.m_new)} " print "n_n0_m_new: ${df7_8.format(it.n_n0_m_new)}" println "abs_n_n0_m_new: ${df7_8.format(it.abs_n_n0_m_new)}" } println '---' println "lat : ${sc10dp.format(lat)}" println "sin_sq_lat : ${sc10dp.format(sin_sq_lat)}" println "one_m_e_sq_sin_sq_lat: ${df1_10.format(one_m_e_sq_sin_sq_lat)}" println "v : ${sc10dp.format(v)}" println "p : ${sc10dp.format(p)}" println "n_sq : ${sc10dp.format(n_sq)}" print "tan_lat : ${sc10dp.format(tan_lat)} " print "tan_lat_p2: ${sc10dp.format(tan_lat_p2)} " print "tan_lat_p4: ${sc10dp.format(tan_lat_p4)} " print "tan_lat_p6: ${sc10dp.format(tan_lat_p6)} " println "sec_lat: ${sc10dp.format(sec_lat)}" print "vii : ${sc10dp.format(vii)} " print "viii : ${sc10dp.format(viii)} " print "ix : ${sc10dp.format(ix)} " print "x : ${sc10dp.format(x)} " print "xi : ${sc10dp.format(xi)} " print "xii : ${sc10dp.format(xii)} " println "xiia : ${sc10dp.format(xiia)} " print "latres: ${sc10dp.format(latres)} longres: ${sc10dp.format(longres)} " print "latresdegdec: ${df1_9.format(latresdegdec)} longresdegdec: ${df1_9.format(longresdegdec)} " println "latresdms: $latresdms longresdms: $longresdms" return null } def converDegreesDecToDMS(degreesIn) { def dms = [:] def degrees = degreesIn.toInteger() def fract = degreesIn - degrees def minutes = (fract * 60).toInteger() def seconds = fract *3600 - (minutes * 60) dms.degrees = degrees dms.minutes = minutes dms.seconds = seconds return dms } } def c = new east_north_to_lat_long() c.convert(651409.903, 313177.270, airy1830)
Here is the result of running the canned example as described in section C of conversion guide and my script.
(see second of ‘useful links’ below).
Basically the script hones in on a more and more precise value for the latitude in an iterative process, by applying delta adjustments.
You have to be careful though, because you could end up in an infinite loop if you try to get too accurate. For example I was able to tweak the tolerance field ‘tol‘ from 0.01. I tried 0.00001 and got another iteration through.
Your mileage may vary depending on the example you choose. You could try 485345 for Easting and 234472 for Northing and you’d know whereabouts I live. :-)
Version 2 (en_to_ll2.groovy):
I have since gone on to enhance this script to convert latitude & longitudes to the WGS84 datum format that most mapping systems use. This involved converting to Cartisian Co-ordinates, applying a Helmert Transformation and converting back to Latitiude & Longitude again. This facilitates a transformation from OSGB36 latitude & longitude to WGS84. Chris Veness’s blog post provided some of the missing pieces that didn’t jump out at me from the Ordnance survey guide. Namely how to approximate the level of precision.
- I’ve also made tweak to the converDegreesDecToDMS routine, renaming it convertRadtoDegrees and accomodated negative values with Compass points in a more robust fashion. See results below.
- I also factored out the Ellipsoid and TransverseMercatorProjection classes to make the script more O-O. I also hit upon some teething problems with incorrect assumptions about Groovy constructors and class initialization with Maps. See my other post for the ‘gotcha’.
package jgf import java.text.DecimalFormat // Transverse Mercator projection Constants used in UK // Taken from P40 of Ordnance Survey : A guide to coordinate systems in Great Britain // E is positive, W negative. Since UK in Northern hemisphere lat always plus for UK Post codes..... // Airy 1830 for UK Mainland. Airy1830mod for Ireland... def airy1830 = [a: 6377563.396, b: 6356256.910, f0:0.9996012717, lat0deg:49, long0deg:-2, e0:400000, n0:-100000] def airy1830mod = [a: 6377340.189, b: 6356034.447, f0:1.000035, lat0deg:53.5, long0deg:-8, e0:200000, n0:250000] def i1924_ed50_utm30_31 = [a: 6378388.000, b: 6356911.946, f0:0.9996, lat0deg:0, long0deg:-3, e0:500000, n0:0] def utm29 = [a: 6378388.000, b: 6356911.946, f0:0.9996, lat0deg:0, long0deg:-9, e0:500000, n0:0] // Ellipsoid constants def WGS84 = [a: 6378137 , b: 6356752.3142] // Helmert Transformation def OSGB36toWGS84 = [tx: 446.448, ty: -125.157, tz: 542.060, rx: 0.1502, ry: 0.2470, rz: 0.8421, s: -20.4894] enum LatLong { LATITUDE, LONGITUDE } class Ellipsoid { double a double b Ellipsoid(props) { this.a = props.a this.b = props.b } double a(pow) { return Math.pow(a, pow) } double b(pow) { return Math.pow(b, pow) } double n() { return n(1) } double n(pow) { return Math.pow(((a - b) / (a + b)), pow) } double esq() { def ap2 = a(2) return (ap2 - b(2))/ap2 } double one_m_esq() { return (1.0 - esq()) } } class TransverseMercatorProjection { Ellipsoid e double f0 // Scale factor on central meridian F0 double lat0deg double long0deg //double lat0rad //double long0rad int e0 int n0 TransverseMercatorProjection(props) { this.e = new Ellipsoid(props) this.f0 = props.f0 this.lat0deg = props.lat0deg this.long0deg = props.long0deg this.e0 = props.e0 this.n0 = props.n0 } double lat0rad() { return Math.toRadians(lat0deg) } double long0rad() { return Math.toRadians(long0deg) } double a() { return e.a } double b() { return e.b } double a(pow) { return e.a(pow) } double b(pow) { return e.b(pow) } double n() { return e.n() } double n(pow) { return e.n(pow) } double esq() { return e.esq() } double one_m_esq() { return e.one_m_esq() } double af0() { return e.a * f0 } double bf0() { return e.b * f0 } } class EastNorthToLatLong { def convertEN(east, north, tmpPar, trace) { // Number formats for various intermediary results def sc10dp = new DecimalFormat("0.#########E0") def df7_8 = new DecimalFormat("#,###,###.########") def df6_3 = new DecimalFormat("###,###.###") def df1_9 = new DecimalFormat("0.#########") def df1_10 = new DecimalFormat("0.##########") def df1_11 = new DecimalFormat("0.###########") def tmp = new TransverseMercatorProjection(tmpPar) def e_e0s = [] double e_e0 = east - tmp.e0 (0..7).each{e_e0s << Math.pow(e_e0, it)} // Don't bother with e_e0s[0].. will contain 1 (n^0 = 1) double n_n0 = north - tmp.n0 double m_old = 0.0 def latseed = tmp.lat0rad() def latresults = [] double lat double tol = 0.01 double diff = 0.02 while (diff > tol) { double n_n0_m_old = n_n0 - m_old double n_n0_m_old_o_af0 = n_n0_m_old / tmp.af0() lat = n_n0_m_old_o_af0 + latseed // In radians double latplat0 = lat + tmp.lat0rad() double latmlat0 = lat - tmp.lat0rad() double m_1 = (1 + tmp.n() + ((5/4) * tmp.n(2)) + ((5/4) * tmp.n(3))) * latmlat0 double m_2 = ((3 * tmp.n()) +(3 * tmp.n(2)) +((21/8) * tmp.n(3))) * Math.sin(latmlat0) * Math.cos(latplat0) double m_3 = ( ((15/8) * tmp.n(2)) + ((15/8) * tmp.n(3)) ) * Math.sin(2 *latmlat0) * Math.cos(2 * latplat0) double m_4 = ( (35/24) * tmp.n(3)) * Math.sin(3 *latmlat0) * Math.cos(3 * latplat0) double m_new = tmp.bf0() * (m_1 - m_2 + m_3 - m_4) double n_n0_m_new = n_n0 - m_new double abs_n_n0_m_new = Math.abs(n_n0_m_new) def latres = [:] latres.m_old = m_old latres.latseed = latseed latres.lat = lat latres.latplat0 = latplat0 latres.latmlat0 = latmlat0 latres.m_1 = m_1 latres.m_2 = m_2 latres.m_3 = m_3 latres.m_4 = m_4 latres.m_new = m_new latres.n_n0_m_new = n_n0_m_new latres.abs_n_n0_m_new = abs_n_n0_m_new latresults << latres diff = abs_n_n0_m_new m_old = m_new latseed = lat } lat = latresults[latresults.size() -1].lat double sin_sq_lat = Math.pow(Math.sin(lat), 2) double one_m_e_sq_sin_sq_lat = 1 - (tmp.esq() * Math.pow(Math.sin(lat), 2)) double v = tmp.af0() * Math.pow(one_m_e_sq_sin_sq_lat,-0.5) double p = tmp.af0() * tmp.one_m_esq() * Math.pow(one_m_e_sq_sin_sq_lat,-1.5) double n_sq = (v/p) -1 double tan_lat = Math.tan(lat) double tan_lat_p2 = Math.pow(tan_lat, 2) double tan_lat_p4 = Math.pow(tan_lat, 4) double tan_lat_p6 = Math.pow(tan_lat, 6) double sec_lat = 1 / (Math.cos(lat)) double vii = tan_lat / (2 * p * v) double viii = (tan_lat / (24 * p * Math.pow(v,3)) ) * (5 + (3 * tan_lat_p2) + n_sq - (9 * tan_lat_p2 * n_sq)) double ix = (tan_lat/(720 * p * Math.pow(v,5))) * (61 + (90 * tan_lat_p2) + 45 * tan_lat_p4) double x = sec_lat / v double xi = (sec_lat / (6 * Math.pow(v,3))) * ((v/p) + 2 * tan_lat_p2) double xii = (sec_lat/(120 * Math.pow(v,5))) * (5 + (28 * tan_lat_p2) + (24 * tan_lat_p4)) double xiia = (sec_lat/(5040 * Math.pow(v,7))) * (61 + (662 * tan_lat_p2) + (1320 * tan_lat_p4) + (720 * tan_lat_p6)) double latrad = lat - (vii * e_e0s[2]) + (viii * e_e0s[4]) - (ix * e_e0s[6]) // =_lat1.2-(_vii*_E_E0_p2)+(_viii*_E_E0_p4)-(_ix*_E_E0_p6) double longrad = tmp.long0rad() + (x * e_e0s[1]) - (xi * e_e0s[3]) + (xii * e_e0s[5]) - (xiia * e_e0s[7]) // =_long0+(_x*_E_E0)-(_xi*_E_E0_p3)+(_xii*_E_E0_p5)-(_x11a*_E_E0_p7) def res = [:] res.latitude = convertRadtoDegrees(latrad, LatLong.LATITUDE) res.longitude = convertRadtoDegrees(longrad, LatLong.LONGITUDE) if (trace) { // Dump out intermediary calculations... println "tol : ${df1_9.format(tol)}" println "pr.lat0 : ${df1_9.format(tmp.lat0rad())}" println "pr.long0 : ${df1_9.format(tmp.long0rad())}" println "one_m_e_sq : ${df1_10.format(tmp.one_m_esq())}" println "e_sq : ${sc10dp.format(tmp.esq())}" println "n : ${sc10dp.format(tmp.n())}" println "n_p2 : ${sc10dp.format(tmp.n(2))}" println "n_p3 : ${sc10dp.format(tmp.n(3))}" println "af0 : ${df7_8.format(tmp.af0())}" println "bf0 : ${df7_8.format(tmp.bf0())}" e_e0s.eachWithIndex{e0, ei -> print "e_e0^${ei+1} : ${(ei) ? sc10dp.format(e0) : df6_3.format(e0)} " } println " " println "n_n0 : ${df6_3.format(n_n0)}" println "---" latresults.eachWithIndex{it, i -> print "${i+1}: m_old: ${sc10dp.format(it.m_old)} " print "latseed: ${sc10dp.format(it.latseed)} " print "lat: ${sc10dp.format(it.lat)} " print "latplat0: ${df1_9.format(it.latplat0)} " print "latmlat0: ${df1_9.format(it.latmlat0)} " print "m_1: ${sc10dp.format(it.m_1)} m_2: ${sc10dp.format(it.m_2)} m_3: ${sc10dp.format(it.m_3)} m_4: ${sc10dp.format(it.m_4)} " print "m_new: ${sc10dp.format(it.m_new)} " print "n_n0_m_new: ${df7_8.format(it.n_n0_m_new)}" println "abs_n_n0_m_new: ${df7_8.format(it.abs_n_n0_m_new)}" } println '---' println "lat : ${sc10dp.format(lat)}" println "sin_sq_lat : ${sc10dp.format(sin_sq_lat)}" println "one_m_e_sq_sin_sq_lat: ${df1_10.format(one_m_e_sq_sin_sq_lat)}" println "v : ${sc10dp.format(v)}" println "p : ${sc10dp.format(p)}" println "n_sq : ${sc10dp.format(n_sq)}" print "tan_lat : ${sc10dp.format(tan_lat)} " print "tan_lat_p2: ${sc10dp.format(tan_lat_p2)} " print "tan_lat_p4: ${sc10dp.format(tan_lat_p4)} " print "tan_lat_p6: ${sc10dp.format(tan_lat_p6)} " println "sec_lat: ${sc10dp.format(sec_lat)}" print "vii : ${sc10dp.format(vii)} " print "viii : ${sc10dp.format(viii)} " print "ix : ${sc10dp.format(ix)} " print "x : ${sc10dp.format(x)} " print "xi : ${sc10dp.format(xi)} " print "xii : ${sc10dp.format(xii)} " println "xiia : ${sc10dp.format(xiia)} " println "latrad: ${sc10dp.format(latrad)} longres: ${sc10dp.format(longrad)} " println "res: $res" } return res } def convertENtoWGS84(enh, tmpPar, ePar, ht, trace) { // tmp Airy1830, e = WGS84, ht = Helmert Transformation if (trace) println tmpPar def tmp = new TransverseMercatorProjection(tmpPar) def e = new Ellipsoid(ePar) def osgb36 = convertEN(enh.e, enh.n, tmpPar, trace) def wgs84 = convertOSGB36toWGS84(osgb36, tmp, e, ht, enh.h, trace) def res = [:] res.e = enh.e res.n = enh.n res.hOSGB = enh.h res.hWGS = wgs84.wgs84ll.h res.latOSGBdms = osgb36.latitude.dms res.latOSGBdd = osgb36.latitude.degreesReal res.longOSGBdms = osgb36.longitude.dms res.longOSGBdd = osgb36.longitude.degreesReal res.latWGSdms = wgs84.wgs84ll.latitude.dms res.latWGSdd = wgs84.wgs84ll.latitude.degreesReal res.longWGSdms = wgs84.wgs84ll.longitude.dms res.longWGSdd = wgs84.wgs84ll.longitude.degreesReal prettyPrint(res) return res } def convertOSGB36toWGS84(osgb36,tmp,e,ht, h, trace) { def osgb36Cart = convertOSGB36ToCartesian(osgb36, tmp, h) def wgs84Cart = applyHelmertTransformation(osgb36Cart, ht) def wgs84ll = convertWGS84CartToLatLong(e, wgs84Cart, trace) def res = [:] res.osgb36Cart = osgb36Cart res.wgs84Cart = wgs84Cart res.wgs84ll = wgs84ll if (trace) { println "res.osgb36Cart: ${res.osgb36Cart}" println "res.wgs84Cart : ${res.wgs84Cart}" println "res.wgs84ll : ${res.wgs84ll}" } return res } def convertOSGB36ToCartesian(osgb36, tmp, h) { double sinlat = Math.sin(osgb36.latitude.rad) double coslat = Math.cos(osgb36.latitude.rad) double sinlong = Math.sin(osgb36.longitude.rad) double coslong = Math.cos(osgb36.longitude.rad) double osgb36v = tmp.a()*Math.pow((1 - (tmp.esq()*Math.pow(sinlat,2))),-0.5) double v_p_h = osgb36v + h double x = v_p_h * coslat * coslong double y = v_p_h * coslat * sinlong double z = ((tmp.one_m_esq() * osgb36v) + h) * sinlat def res = [:] res.sinlat = sinlat res.coslat = coslat res.sinlong = sinlong res.coslong = coslong res.v = osgb36v res.v_p_h = v_p_h res.x = x res.y = y res.z = z return res } def applyHelmertTransformation(osgb36Cart, ht) { double rx = convertHTCoord(ht.rx) double ry = convertHTCoord(ht.ry) double rz = convertHTCoord(ht.rz) double s = 1.0 +(ht.s / Math.pow(10,6)) double x = ht.tx + osgb36Cart.x * s - osgb36Cart.y * rz + osgb36Cart.z * ry // _tx+_x1*_s-_y1*_rz+_z1*_ry double y = ht.ty + osgb36Cart.x * rz + osgb36Cart.y * s - osgb36Cart.z * rx // _ty+_x1*_rz+_y1*_s-_z1*_rx double z = ht.tz - osgb36Cart.x * ry + osgb36Cart.y * rx + osgb36Cart.z * s // _tz-_x1*_ry+_y1*_rx+_z1*_s def res = [:] res.x = x res.y = y res.z = z return res } def convertHTCoord(coord) { double res = ((coord / 3600) * (Math.PI / 180)) return res } def convertWGS84CartToLatLong(e, wgs84Cart, trace) { double tol = 4 / e.a double p = Math.pow( (Math.pow(wgs84Cart.x, 2) + Math.pow(wgs84Cart.y, 2)), 0.5) double lat = Math.atan(wgs84Cart.z/p*e.one_m_esq()) double latPrior = 2 * Math.PI def res = [:] res.p = p res.latInitial = lat res.latPrior = latPrior res.tol = tol def latresults = [], latresult = [:], v latresult.lat = lat latresult.latPrior = latPrior updDiffCont(latresult, tol) while (latresult.cont) { v = e.a * Math.pow( (1 - e.esq()* Math.pow(Math.sin(lat), 2)) , -0.5) if (trace) println "v: $v lat: $lat" latresult.v = v latresults << latresult latresult = [:] latPrior = lat lat = Math.atan((wgs84Cart.z + e.esq() * v * Math.sin(latPrior))/p) latresult.lat = lat latresult.latPrior = latPrior updDiffCont(latresult, tol) } v = e.a * Math.pow( (1 - e.esq()* Math.pow(Math.sin(lat), 2)) , -0.5) if (trace) println "v: $v lat: $lat" latresult.v = v latresults << latresult res.latresults = latresults def lastIndex = latresults.size() -1 def latitude = latresults[lastIndex].lat def longitude = Math.atan(wgs84Cart.y/wgs84Cart.x) res.h = (p/Math.cos(latitude)) - latresults[lastIndex].v res.latitude = convertRadtoDegrees(latitude, LatLong.LATITUDE) res.longitude = convertRadtoDegrees(longitude, LatLong.LONGITUDE) return res } def updDiffCont(props, tol) { double diff = Math.abs(props.lat - props.latPrior) props.diff = diff props.cont = diff > tol return null } def convertRadtoDegrees(angleInRad, LatLong ll) { def df2dp = new DecimalFormat("0.00") def res = [:] res.rad = angleInRad res.degreesReal = Math.toDegrees(res.rad) def degreesRealAbs = Math.abs(res.degreesReal) res.degrees = degreesRealAbs.toInteger() def fract = degreesRealAbs - res.degrees res.minutes = (fract * 60).toInteger() res.seconds = fract *3600 - (res.minutes * 60) if (ll == LatLong.LATITUDE) { if (res.degreesReal < 0.0) { res.card = 'S' } else { res.card = 'N' } } else if (ll == LatLong.LONGITUDE) { if (res.degreesReal < 0.0) { res.card = 'W' } else { res.card = 'E' } } res.dms = """${res.degrees}°${res.minutes}"${df2dp.format(res.seconds)}'${res.card}""" return res } def prettyPrint(props) { println "E/N: [${props.e}, ${props.n}] H: (OSGB36) $props.hOSGB (WGS84) $props.hWGS" println "OSGB36 Lat/Long Dec : [${props.latOSGBdd}, ${props.longOSGBdd}]" println "OSGB36 Lat/Long DMS : [${props.latOSGBdms}, ${props.longOSGBdms}]" println "WGS84 Lat/Long Dec : [${props.latWGSdd}, ${props.longWGSdd}]" println "WGS84 Lat/Long DMS : [${props.latWGSdms}, ${props.longWGSdms}]" } } def c1 = new EastNorthToLatLong() def enh1 = [e: 651409.903, n:313177.270, h:24.00] def p1 = c1.convertENtoWGS84(enh1, airy1830, WGS84, OSGB36toWGS84, false) println "===" def c2 = new EastNorthToLatLong() def enh2 = [e: 485345, n:234472, h:24.700] def p2 = c2.convertENtoWGS84(enh2, airy1830, WGS84, OSGB36toWGS84, false) return null
Here is the result of running the script:
E/N: [651409.903, 313177.270] H: (OSGB36) 24.00 (WGS84) 68.70079326722771 OSGB36 Lat/Long Dec : [52.657570301518525, 1.7179215809789155] OSGB36 Lat/Long DMS : [52°39"27.25'N, 1°43"4.52'E] WGS84 Lat/Long Dec : [52.65797859738828, 1.716051953792222] WGS84 Lat/Long DMS : [52°39"28.72'N, 1°42"57.79'E] === E/N: [485345, 234472] H: (OSGB36) 24.700 (WGS84) 72.17237068805844 OSGB36 Lat/Long Dec : [52.001694746632644, -0.756629819124975] OSGB36 Lat/Long DMS : [52°0"6.10'N, 0°45"23.87'W] WGS84 Lat/Long Dec : [52.00213785138463, -0.7581842475369207] WGS84 Lat/Long DMS : [52°0"7.70'N, 0°45"29.46'W]
Version 3 (UKPostCodeAggregatorScript.groovy):
- I read through all the CSV files and appended the Latitude & Longitude to a new version of each file and wrote it out to an new TXT directory using ConfigSlurper to pick up the source/destination directories.
crawl.properties values: ukPostCodesCSVDir='/Users/JGF/Downloads/Code-Point Open/data/CSV' ukPostCodesTXTDir='/Users/JGF/Downloads/Code-Point Open/data/TXT'
- I used the SmartCsvParser class that I first described in my earlier post Groovy Script to parse Yahoo contacts in CSV format and export to XML based loosely on an enhanced version of one Scott Davis Groovy Recipes.
- There seemed to be an issue reading the first record of each file that Code-Point provided and it splitting the line up into an array of columnar data, so you’ll see a first time flag to stop an extra column appearing in first row of a CSV file.
- There is a Positional quality indicator that Code-Point uses in the CSV file, which when it has a value of ’90’ it won’t have an Easting and Northing to convert. So I also don’t do the complex conversion for such records either. (See Code-Point Technical Info/Spec links below and comments in code).
- Caveat:
- Opening in Excel doesn’t work well, because Microsoft didn’t think to give the option of opening a UTF-8 encoded file, so the degree symbol \u00B0 doesn’t look good. :-(
- But if you open it up with a more intelligent piece of software like TextMate it looks good. :-)
package jgf import java.text.DecimalFormat class UKPostCodeAggregator { def init() { def res = [:] res.with { h = System.getenv('HOME') // OS Shell var fs = System.getProperty('file.separator') // Java Sys Property nl = System.getProperty("line.separator") // Newline character d = "${h}${fs}Desktop" gsd = "${d}${fs}Groovy Scripts" def props = new ConfigSlurper().parse(new File("${gsd}${fs}crawl.properties").toURL()) ukPostCodesCSVDir = props.ukPostCodesCSVDir ukPostCodesTXTDir = props.ukPostCodesTXTDir encoding = 'UTF-8' numcols = 19 /* A-S 0-18 in array... 0 = pc, 10 = east, 11 = north, * 19 lat dec OSGB36 * 20 long dec OSGB36 * 21 lat dms OSGB36 * 22 long dms OSGB36 * 23 lat dec WGS84 * 24 long dec WGS84 * 25 lat dms WGS84 * 26 long dms WGS84 */ e = 11 n = 12 /* Transverse Mercator projection Constants used in UK * Taken from P40 of Ordnance Survey : A guide to coordinate systems in Great Britain * E is positive, W negative. Since UK in Northern hemisphere lat always plus for UK Post codes..... * Airy 1830 for UK Mainland. Airy1830mod for Ireland... */ airy1830 = [a: 6377563.396, b: 6356256.910, f0:0.9996012717, lat0deg:49, long0deg:-2, e0:400000, n0:-100000] airy1830mod = [a: 6377340.189, b: 6356034.447, f0:1.000035, lat0deg:53.5, long0deg:-8, e0:200000, n0:250000] i1924_ed50_utm30_31 = [a: 6378388.000, b: 6356911.946, f0:0.9996, lat0deg:0, long0deg:-3, e0:500000, n0:0] utm29 = [a: 6378388.000, b: 6356911.946, f0:0.9996, lat0deg:0, long0deg:-9, e0:500000, n0:0] // Ellipsoid constants WGS84 = [a: 6378137 , b: 6356752.3142] // Helmert Transformation OSGB36toWGS84 = [tx: 446.448, ty: -125.157, tz: 542.060, rx: 0.1502, ry: 0.2470, rz: 0.8421, s: -20.4894] latdegdecfmt = new DecimalFormat(" 00.0000000000;-00.0000000000") longdegdecfmt = new DecimalFormat(" 000.0000000000;-000.0000000000") linefmt = new DecimalFormat("00000") latdmsblank = ' '*13 longdmsblank = ' '*14 latdecblank = ' '*14 longdecblank = ' '*15 textQualifier = '' delimiter = ',' } return res } def convertFilesInUKPCDir() { def props = init() new File(props.ukPostCodesCSVDir).eachFile{file -> def fn = file.getName() if (file.isFile() && fn.endsWith('.csv')) { convertFile(file, props) } } } def convertFile(file, props) { def oldname = file.getName() def token = ' ' def newname = "${oldname[0..-5]}.txt" def newfull = "${props.ukPostCodesTXTDir}${props.fs}${newname}" println newfull def newFile = new File(newfull) newFile.delete() println newname def first = true def lineNo = 1 file.eachLine{line -> def fd, fd2 def fields = [] use(SmartCsvParser) { fd = line.smartSplit(token, props.nl, first) } first = false token = fd.token fd2 = fd.list fd2.each{fields << it} if (fields.size() == props.numcols) { processLine(fields, props, newFile, lineNo) } lineNo +=1 } } def processLine(fields, props, newFile, lineNo) { def res = new StringBuilder() def tqDel = "${props.textQualifier}${props.delimiter}" fields.each{field -> res += props.textQualifier res += field res += tqDel } def enll = new EastNorthToLatLong() def enh = [e: Double.parseDouble(fields[props.e]),n: Double.parseDouble(fields[props.n]), h:24.00] def p1 = enll.convertENtoWGS84(enh, props.airy1830, props.WGS84, props.OSGB36toWGS84) /* See 1: http://www.ordnancesurvey.co.uk/oswebsite/products/codepoint/techinfo.html * 2: http://www.ordnancesurvey.co.uk/oswebsite/products/codepoint/pdf/Code-Point_userguidev2-7.pdf P10 */ def pqi = fields[1] def noCoordinates = (pqi == '90') res += props.textQualifier res += (noCoordinates) ? props.latdecblank : props.latdegdecfmt.format(p1.latOSGBdd) res += tqDel res += props.textQualifier res += (noCoordinates) ? props.longdecblank : props.longdegdecfmt.format(p1.longOSGBdd) res += tqDel res += props.textQualifier res += (noCoordinates) ? props.latdmsblank : p1.latOSGBdms res += tqDel res += props.textQualifier res += (noCoordinates) ? props.longdmsblank : p1.longOSGBdms res += tqDel res += props.textQualifier res += (noCoordinates) ? props.latdecblank : props.latdegdecfmt.format(p1.latWGSdd) res += tqDel res += props.textQualifier res += (noCoordinates) ? props.longdecblank : props.longdegdecfmt.format(p1.longWGSdd) res += tqDel res += props.textQualifier res += (noCoordinates) ? props.latdmsblank : p1.latWGSdms res += tqDel res += props.textQualifier res += (noCoordinates) ? props.longdmsblank : p1.longWGSdms res += "${props.textQualifier}${props.nl}" newFile.append(res, props.encoding) println "${props.linefmt.format(lineNo)} ${fields[0]} $p1" } } class SmartCsvParser{ static def smartSplit(String line, String midToken, String nl, Boolean first){ def list = [] def thisToken def st = new StringTokenizer(line, ",") def ft = first while (st.hasMoreTokens()) { if (midToken) { thisToken = midToken midToken = '' } else { thisToken = st.nextToken() } while ( thisToken.startsWith("\"") && !thisToken.endsWith("\"") && st.hasMoreTokens() ) { thisToken += "," + st.nextToken() } if (thisToken.startsWith("\"") && !thisToken.endsWith("\"")) midToken = thisToken.replaceAll(',', nl) else if (!ft) list << thisToken.noQuote() ft = false } def res = [:] res.list = list res.token = midToken return res } static String noQuote(String token) { if(token.startsWith("\"")) { if (token.size() == 2) return '' else return token[1..-2] } else return token } } enum LatLong { LATITUDE, LONGITUDE } class Ellipsoid { double a double b Ellipsoid(props) { this.a = props.a this.b = props.b } double a(pow) { return Math.pow(a, pow) } double b(pow) { return Math.pow(b, pow) } double n() { return n(1) } double n(pow) { return Math.pow(((a - b) / (a + b)), pow) } double esq() { def ap2 = a(2) return (ap2 - b(2))/ap2 } double one_m_esq() { return (1.0 - esq()) } } class TransverseMercatorProjection { Ellipsoid e double f0 // Scale factor on central meridian F0 double lat0deg double long0deg //double lat0rad //double long0rad int e0 int n0 TransverseMercatorProjection(props) { this.e = new Ellipsoid(props) this.f0 = props.f0 this.lat0deg = props.lat0deg this.long0deg = props.long0deg this.e0 = props.e0 this.n0 = props.n0 } double lat0rad() { return Math.toRadians(lat0deg) } double long0rad() { return Math.toRadians(long0deg) } double a() { return e.a } double b() { return e.b } double a(pow) { return e.a(pow) } double b(pow) { return e.b(pow) } double n() { return e.n() } double n(pow) { return e.n(pow) } double esq() { return e.esq() } double one_m_esq() { return e.one_m_esq() } double af0() { return e.a * f0 } double bf0() { return e.b * f0 } } class EastNorthToLatLong { def convertEN(east, north, tmpPar) { def tmp = new TransverseMercatorProjection(tmpPar) def e_e0s = [] double e_e0 = east - tmp.e0 (0..7).each{e_e0s << Math.pow(e_e0, it)} // Don't bother with e_e0s[0].. will contain 1 (n^0 = 1) double n_n0 = north - tmp.n0 double m_old = 0.0 def latseed = tmp.lat0rad() def latresults = [] double lat double tol = 0.01 double diff = 0.02 while (diff > tol) { double n_n0_m_old = n_n0 - m_old double n_n0_m_old_o_af0 = n_n0_m_old / tmp.af0() lat = n_n0_m_old_o_af0 + latseed // In radians double latplat0 = lat + tmp.lat0rad() double latmlat0 = lat - tmp.lat0rad() double m_1 = (1 + tmp.n() + ((5/4) * tmp.n(2)) + ((5/4) * tmp.n(3))) * latmlat0 double m_2 = ((3 * tmp.n()) +(3 * tmp.n(2)) +((21/8) * tmp.n(3))) * Math.sin(latmlat0) * Math.cos(latplat0) double m_3 = ( ((15/8) * tmp.n(2)) + ((15/8) * tmp.n(3)) ) * Math.sin(2 *latmlat0) * Math.cos(2 * latplat0) double m_4 = ( (35/24) * tmp.n(3)) * Math.sin(3 *latmlat0) * Math.cos(3 * latplat0) double m_new = tmp.bf0() * (m_1 - m_2 + m_3 - m_4) double n_n0_m_new = n_n0 - m_new double abs_n_n0_m_new = Math.abs(n_n0_m_new) def latres = [:] latres.m_old = m_old latres.latseed = latseed latres.lat = lat latres.latplat0 = latplat0 latres.latmlat0 = latmlat0 latres.m_1 = m_1 latres.m_2 = m_2 latres.m_3 = m_3 latres.m_4 = m_4 latres.m_new = m_new latres.n_n0_m_new = n_n0_m_new latres.abs_n_n0_m_new = abs_n_n0_m_new latresults << latres diff = abs_n_n0_m_new m_old = m_new latseed = lat } lat = latresults[latresults.size() -1].lat double sin_sq_lat = Math.pow(Math.sin(lat), 2) double one_m_e_sq_sin_sq_lat = 1 - (tmp.esq() * Math.pow(Math.sin(lat), 2)) double v = tmp.af0() * Math.pow(one_m_e_sq_sin_sq_lat,-0.5) double p = tmp.af0() * tmp.one_m_esq() * Math.pow(one_m_e_sq_sin_sq_lat,-1.5) double n_sq = (v/p) -1 double tan_lat = Math.tan(lat) double tan_lat_p2 = Math.pow(tan_lat, 2) double tan_lat_p4 = Math.pow(tan_lat, 4) double tan_lat_p6 = Math.pow(tan_lat, 6) double sec_lat = 1 / (Math.cos(lat)) double vii = tan_lat / (2 * p * v) double viii = (tan_lat / (24 * p * Math.pow(v,3)) ) * (5 + (3 * tan_lat_p2) + n_sq - (9 * tan_lat_p2 * n_sq)) double ix = (tan_lat/(720 * p * Math.pow(v,5))) * (61 + (90 * tan_lat_p2) + 45 * tan_lat_p4) double x = sec_lat / v double xi = (sec_lat / (6 * Math.pow(v,3))) * ((v/p) + 2 * tan_lat_p2) double xii = (sec_lat/(120 * Math.pow(v,5))) * (5 + (28 * tan_lat_p2) + (24 * tan_lat_p4)) double xiia = (sec_lat/(5040 * Math.pow(v,7))) * (61 + (662 * tan_lat_p2) + (1320 * tan_lat_p4) + (720 * tan_lat_p6)) double latrad = lat - (vii * e_e0s[2]) + (viii * e_e0s[4]) - (ix * e_e0s[6]) // =_lat1.2-(_vii*_E_E0_p2)+(_viii*_E_E0_p4)-(_ix*_E_E0_p6) double longrad = tmp.long0rad() + (x * e_e0s[1]) - (xi * e_e0s[3]) + (xii * e_e0s[5]) - (xiia * e_e0s[7]) // =_long0+(_x*_E_E0)-(_xi*_E_E0_p3)+(_xii*_E_E0_p5)-(_x11a*_E_E0_p7) def res = [:] res.latitude = convertRadtoDegrees(latrad, LatLong.LATITUDE) res.longitude = convertRadtoDegrees(longrad, LatLong.LONGITUDE) return res } def convertENtoWGS84(enh, tmpPar, ePar, ht) { // tmp Airy1830, e = WGS84, ht = Helmert Transformation def tmp = new TransverseMercatorProjection(tmpPar) def e = new Ellipsoid(ePar) def osgb36 = convertEN(enh.e, enh.n, tmpPar) def wgs84 = convertOSGB36toWGS84(osgb36, tmp, e, ht, enh.h) def res = [:] res.latOSGBdms = osgb36.latitude.dms res.latOSGBdd = osgb36.latitude.degreesReal res.longOSGBdms = osgb36.longitude.dms res.longOSGBdd = osgb36.longitude.degreesReal res.latWGSdms = wgs84.wgs84ll.latitude.dms res.latWGSdd = wgs84.wgs84ll.latitude.degreesReal res.longWGSdms = wgs84.wgs84ll.longitude.dms res.longWGSdd = wgs84.wgs84ll.longitude.degreesReal return res } def convertOSGB36toWGS84(osgb36,tmp,e,ht, h) { def osgb36Cart = convertOSGB36ToCartesian(osgb36, tmp, h) def wgs84Cart = applyHelmertTransformation(osgb36Cart, ht) def wgs84ll = convertWGS84CartToLatLong(e, wgs84Cart) def res = [:] res.osgb36Cart = osgb36Cart res.wgs84Cart = wgs84Cart res.wgs84ll = wgs84ll return res } def convertOSGB36ToCartesian(osgb36, tmp, h) { double sinlat = Math.sin(osgb36.latitude.rad) double coslat = Math.cos(osgb36.latitude.rad) double sinlong = Math.sin(osgb36.longitude.rad) double coslong = Math.cos(osgb36.longitude.rad) double osgb36v = tmp.a()*Math.pow((1 - (tmp.esq()*Math.pow(sinlat,2))),-0.5) double v_p_h = osgb36v + h double x = v_p_h * coslat * coslong double y = v_p_h * coslat * sinlong double z = ((tmp.one_m_esq() * osgb36v) + h) * sinlat def res = [:] res.sinlat = sinlat res.coslat = coslat res.sinlong = sinlong res.coslong = coslong res.v = osgb36v res.v_p_h = v_p_h res.x = x res.y = y res.z = z return res } def applyHelmertTransformation(osgb36Cart, ht) { double rx = convertHTCoord(ht.rx) double ry = convertHTCoord(ht.ry) double rz = convertHTCoord(ht.rz) double s = 1.0 +(ht.s / Math.pow(10,6)) double x = ht.tx + osgb36Cart.x * s - osgb36Cart.y * rz + osgb36Cart.z * ry // _tx+_x1*_s-_y1*_rz+_z1*_ry double y = ht.ty + osgb36Cart.x * rz + osgb36Cart.y * s - osgb36Cart.z * rx // _ty+_x1*_rz+_y1*_s-_z1*_rx double z = ht.tz - osgb36Cart.x * ry + osgb36Cart.y * rx + osgb36Cart.z * s // _tz-_x1*_ry+_y1*_rx+_z1*_s def res = [:] res.x = x res.y = y res.z = z return res } def convertHTCoord(coord) { double res = ((coord / 3600) * (Math.PI / 180)) return res } def convertWGS84CartToLatLong(e, wgs84Cart) { double tol = 4 / e.a double p = Math.pow( (Math.pow(wgs84Cart.x, 2) + Math.pow(wgs84Cart.y, 2)), 0.5) double lat = Math.atan(wgs84Cart.z/p*e.one_m_esq()) double latPrior = 2 * Math.PI def res = [:] res.p = p res.latInitial = lat res.latPrior = latPrior res.tol = tol def latresults = [], latresult = [:], v latresult.lat = lat latresult.latPrior = latPrior updDiffCont(latresult, tol) while (latresult.cont) { v = e.a * Math.pow( (1 - e.esq()* Math.pow(Math.sin(lat), 2)) , -0.5) latresult.v = v latresults << latresult latresult = [:] latPrior = lat lat = Math.atan((wgs84Cart.z + e.esq() * v * Math.sin(latPrior))/p) latresult.lat = lat latresult.latPrior = latPrior updDiffCont(latresult, tol) } v = e.a * Math.pow( (1 - e.esq()* Math.pow(Math.sin(lat), 2)) , -0.5) latresult.v = v latresults << latresult res.latresults = latresults def lastIndex = latresults.size() -1 def latitude = latresults[lastIndex].lat def longitude = Math.atan(wgs84Cart.y/wgs84Cart.x) res.h = (p/Math.cos(latitude)) - latresults[lastIndex].v res.latitude = convertRadtoDegrees(latitude, LatLong.LATITUDE) res.longitude = convertRadtoDegrees(longitude, LatLong.LONGITUDE) return res } def updDiffCont(props, tol) { double diff = Math.abs(props.lat - props.latPrior) props.diff = diff props.cont = diff > tol return null } def convertRadtoDegrees(angleInRad, LatLong ll) { def longdegfmt = new DecimalFormat("000") def minfmt = new DecimalFormat("00") def secfmt = new DecimalFormat("00.00") def degfmt = (ll == LatLong.LATITUDE) ? minfmt : longdegfmt def res = [:] res.rad = angleInRad res.degreesReal = Math.toDegrees(res.rad) def degreesRealAbs = Math.abs(res.degreesReal) res.degrees = degreesRealAbs.toInteger() def fract = degreesRealAbs - res.degrees res.minutes = (fract * 60).toInteger() res.seconds = fract *3600 - (res.minutes * 60) if (ll == LatLong.LATITUDE) { if (res.degreesReal < 0.0) { res.card = 'S' } else { res.card = 'N' } } else if (ll == LatLong.LONGITUDE) { if (res.degreesReal < 0.0) { res.card = 'W' } else { res.card = 'E' } } res.dms = """${degfmt.format(res.degrees)}\u00B0${minfmt.format(res.minutes)}"${secfmt.format(res.seconds)}'${res.card}""" return res } } def ukpcc = new UKPostCodeAggregator() ukpcc.convertFilesInUKPCDir()
I found the following links useful when writing these scripts:
- Ordnance Survey Related:
- A guide to coordinate systems in Great Britain PDF
- Subsection on converting between co-ordinate systems
- Coordinate transformer page to try out conversions
- Excel 2002 Spreadsheet download page with VBScript that doesn’t work on a Mac with Office 2008!
- A guide to coordinate systems in Great Britain Web Page
- Ellipsoid and projection constants
- Historical Earth ellipsoids on Wiki (Cheers Chris Veness for this one)
- Quest Developer page
- Code-Point: The database with the 120 CSV’s, one per post code area, that contain Eastings/Northings per post code in it’s record layout. I got 120, not 121 CSV files in my download! They’ve probably includedD for Dublin in that number!!
- Miscellaneous:
- An excellent blog by Chris Veness on this and a lot more GIS related articles. Good place to start to get you head around this stuff
- A web page to test converted latitude & longitude looks right. (Remember E is positive, W is negative. We’re in Northern hemisphere here in UK, so all latitudes should be positive. I made this mistake initially!)
- Trigonometry on Wiki/Reciprocal functions
- COMPUTER GRAPHICS Systems & Concepts (ISBN: 0 201 14656 8 Publisher: Addison-Wesley (1987) Authors: Rod Salmon and Mel Slater). Section 13.4: 3D Transformations P390-397
- Java:
- Formatting:
- Maths:
About this entry
You’re currently reading “Groovy Script demonstrating how to convert from OSGB36 Eastings & Northings to Latitude/Longitude,” an entry on All things Grails and RIA
- Published:
- Saturday, April 24, 2010 / 4:53 pm
- Category:
- Groovy
3 Comments
Jump to comment form | comment rss [?] | trackback uri [?]