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 &amp; Northings to be converted to latitude &amp; 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).

Results from canned example : Converting OSGB36 Eastings/Northings to Latitude/Longitude

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'
```
UK Postcode Directory Structure
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
  Computer Graphics Systems & Concepts
Java:
- Formatting:
- Maths:
  - API
  - Tutorial

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

Tags:: Cartesian Co-Ordinate Conversion, Code-Point, ConfigSlurper, CSV, geo-coding UK Post Codes, Groovy, Helmert Transformation

gr8fanboy 4.24.10 / 5pm

I suspect the formula’s given in the PDF have a more precise variation. The results are within a minute or two walk of my house using the script and the Northing/Easting for my post code provided in the Code-Point database. I’ve noticed the OS Transformer page returns a more accurate result for latitude and longitude. I’ll scour the Excel VB Script some more to see if the results can be made more accurate. There’s some interesting reading material about how annual the UK moves 2.5cm on average what with plate tectonics, so map data from 1920’s say will actually have shifted over time…. It turns out the VB Script returns precisely the same results as my Groovy script. I’d seen that the VB Scripts calculated e squared differently, but it turns out the results was identical. (f0 could have been omitted completed in VBScript for e squared by applying a reduction). The only other tweak was the number of decimal places latitude and longitude was presented with in the Excel spreadsheet which I has to load up on my old Dell laptop. Will search elsewhere to see if there is a more precise formula. I found the answer back on Chris Veness’s blog. Have to transform latitude/longitude from OSGB36 to WGS-84 too. This is now superceded with Revision 2. Will now go on to write consolidation script for CSV’s and apply conversions. Will post this as revision 3.

Mark Harby 9.27.10 / 4pm

Hi Jeremy
Great blog but still battling with converting O/S N&E to Long & Lat.
Do you have a copy of the converted text/csv files after running them through your parser ? Would really dig me out of a hole.

Mark

gr8fanboy 9.27.10 / 5pm

Hi Mark
I just had a look at the size of the converted data.
It’s 342.6MB
Whilst the original is 161.2MB
So sending it to you would take twice as long.
You really are better off running the script and doing the conversion. I don’t want to get into any sort of license infringement by sending you the data either.
There is another reason it’s a good idea to get the OS data too. It’s continuously updated, so if you convert again now, you’ll have the most recently updated information from the new post codes that have come into being since I did the original conversion.