SQL function to convert UK OS coordinates from easting/northing to longitude and latitude

Question

Please can someone post a SQL function to convert easting/northing to longitude/latitude. I know it's incredibly complicated but I haven't found anyone who has documented it in T-SQL.

This javascript code works but I'm having trouble converting it to SQL.

I have 16,000 coordinates and need them all converted to lat/long.

This is what I have so far but it's not getting past the while loop.

DECLARE @east real = 482353, 
        @north real = 213371

DECLARE @a real = 6377563.396, 
        @b real = 6356256.910,
        @F0 real = 0.9996012717,

        @lat0 real = 49*PI()/180, 
        @lon0 real = -2*PI()/180

DECLARE @N0 real = -100000, 
        @E0 real = 400000,
        @e2 real = 1 - (@b*@b)/(@a*@a),
        @n real = (@a-@b)/(@a+@b)

DECLARE @n2 real = @n*@n, 
        @n3 real = @n*@n*@n

DECLARE @lat real = @lat0, 
        @M real = 0

WHILE (@north-@N0-@M >= 0.00001)
BEGIN

    SET @lat = ((@north-@N0-@M)/(@a*@F0)) + @lat

    DECLARE @Ma real = (1 + @n + (5/4)*@n2 + (5/4)*@n3) * (@lat-@lat0),
            @Mb real = (3*@n + 3*@n*@n + (21/8)*@n3) * SIN(@lat-@lat0) * COS(@lat+@lat0),
            @Mc real = ((15/8)*@n2 + (15/8)*@n3) * SIN(2*(@lat-@lat0)) * COS(2*(@lat+@lat0)),
            @Md real = (35/24)*@n3 * SIN(3*(@lat-@lat0)) * COS(3*(@lat+@lat0))

    SET @M = @b * @F0 * (@Ma - @Mb + @Mc - @Md)

END

DECLARE @cosLat real = COS(@lat), 
        @sinLat real = SIN(@lat)

DECLARE @nu real = @a*@F0/sqrt(1-@e2*@sinLat*@sinLat)
DECLARE @rho real = @a*@F0*(1-@e2)/POWER(1-@e2*@sinLat*@sinLat, 1.5)
DECLARE @eta2 real = @nu/@rho-1

DECLARE @tanLat real = tan(@lat)
DECLARE @tan2lat real = @tanLat*@tanLat
DECLARE @tan4lat real = @tan2lat*@tan2lat
DECLARE @tan6lat real = @tan4lat*@tan2lat
DECLARE @secLat real = 1/@cosLat
DECLARE @nu3 real = @nu*@nu*@nu
DECLARE @nu5 real = @nu3*@nu*@nu
DECLARE @nu7 real = @nu5*@nu*@nu
DECLARE @VII real = @tanLat/(2*@rho*@nu)
DECLARE @VIII real = @tanLat/(24*@rho*@nu3)*(5+3*@tan2lat+@eta2-9*@tan2lat*@eta2)
DECLARE @IX real = @tanLat/(720*@rho*@nu5)*(61+90*@tan2lat+45*@tan4lat)
DECLARE @X real = @secLat/@nu
DECLARE @XI real = @secLat/(6*@nu3)*(@nu/@rho+2*@tan2lat)
DECLARE @XII real = @secLat/(120*@nu5)*(5+28*@tan2lat+24*@tan4lat)
DECLARE @XIIA real = @secLat/(5040*@nu7)*(61+662*@tan2lat+1320*@tan4lat+720*@tan6lat)

DECLARE @dE real = (@east-@E0)
DECLARE @dE2 real = @dE*@dE
DECLARE @dE3 real = @dE2*@dE
DECLARE @dE4 real = @dE2*@dE2, 
        @dE5 real = @dE3*@dE2
DECLARE @dE6 real = @dE4*@dE2, 
        @dE7 real = @dE5*@dE2

SET @lat = @lat - @VII*@dE2 + @VIII*@dE4 - @IX*@dE6

DECLARE @lon real = @lon0 + @X*@dE - @XI*@dE3 + @XII*@dE5 - @XIIA*@dE7

SELECT @lon, @lat

Marinamarinade

What have you tried and where are you having trouble>? Your question is a bit "Give me the code" at present – Isobaric

@CL4NCY - Why on earth would you want to do such a complicated calculation in T-SQL? It seems like it would make a lot more sense to do it inside a script or another application! – Endotoxin

See above my SQL so far. Well I have a table full of coordinates which I need to convert to lat/long only once. – Marinamarinade

Hmm. Last time I did this, I ended up throwing the table onto disk, converting them using a batch converter, and throwing them back into the database. Though I now can't find the batch converter I used. Looks like there might be an OS-approved converter available here but you need to register and I don't know whether it's batch. EDIT: Aha! Found the one I used. It was a bit fiddly from what I I remember, but once I got it working, it produced good results. – Venezuela

Thanks Matt. I've seen a few software solutions but I thought someone might know an easier code way of doing this. – Marinamarinade

@CL4NCY It's hard to do in SQL. I seem to remember having to implement a couple of the odder trig functions on SQL Server because they weren't built in when I tried to do something like this years ago. I think for a one-off conversion, you'll be better off using a tool outside the DB. (Also, see my earlier comment, now updated with a link to the batch converter I definitely used last time.) – Venezuela

Thanks, seeing as it's a one-off conversion I've used the javascript code via ajax. Thanks for your help. – Marinamarinade

Answer 1

I ended up using the following javascript functions to convert the values. I know it's not a SQL solution but it did the job for me.

function OSGridToLatLong(E, N) {

  var a = 6377563.396, b = 6356256.910;              // Airy 1830 major & minor semi-axes
  var F0 = 0.9996012717;                             // NatGrid scale factor on central meridian
  var lat0 = 49*Math.PI/180, lon0 = -2*Math.PI/180;  // NatGrid true origin
  var N0 = -100000, E0 = 400000;                     // northing & easting of true origin, metres
  var e2 = 1 - (b*b)/(a*a);                          // eccentricity squared
  var n = (a-b)/(a+b), n2 = n*n, n3 = n*n*n;

  var lat=lat0, M=0;
  do {
    lat = (N-N0-M)/(a*F0) + lat;

    var Ma = (1 + n + (5/4)*n2 + (5/4)*n3) * (lat-lat0);
    var Mb = (3*n + 3*n*n + (21/8)*n3) * Math.sin(lat-lat0) * Math.cos(lat+lat0);
    var Mc = ((15/8)*n2 + (15/8)*n3) * Math.sin(2*(lat-lat0)) * Math.cos(2*(lat+lat0));
    var Md = (35/24)*n3 * Math.sin(3*(lat-lat0)) * Math.cos(3*(lat+lat0));
    M = b * F0 * (Ma - Mb + Mc - Md);                // meridional arc

  } while (N-N0-M >= 0.00001);  // ie until < 0.01mm

  var cosLat = Math.cos(lat), sinLat = Math.sin(lat);
  var nu = a*F0/Math.sqrt(1-e2*sinLat*sinLat);              // transverse radius of curvature
  var rho = a*F0*(1-e2)/Math.pow(1-e2*sinLat*sinLat, 1.5);  // meridional radius of curvature
  var eta2 = nu/rho-1;

  var tanLat = Math.tan(lat);
  var tan2lat = tanLat*tanLat, tan4lat = tan2lat*tan2lat, tan6lat = tan4lat*tan2lat;
  var secLat = 1/cosLat;
  var nu3 = nu*nu*nu, nu5 = nu3*nu*nu, nu7 = nu5*nu*nu;
  var VII = tanLat/(2*rho*nu);
  var VIII = tanLat/(24*rho*nu3)*(5+3*tan2lat+eta2-9*tan2lat*eta2);
  var IX = tanLat/(720*rho*nu5)*(61+90*tan2lat+45*tan4lat);
  var X = secLat/nu;
  var XI = secLat/(6*nu3)*(nu/rho+2*tan2lat);
  var XII = secLat/(120*nu5)*(5+28*tan2lat+24*tan4lat);
  var XIIA = secLat/(5040*nu7)*(61+662*tan2lat+1320*tan4lat+720*tan6lat);

  var dE = (E-E0), dE2 = dE*dE, dE3 = dE2*dE, dE4 = dE2*dE2, dE5 = dE3*dE2, dE6 = dE4*dE2, dE7 = dE5*dE2;
  lat = lat - VII*dE2 + VIII*dE4 - IX*dE6;
  var lon = lon0 + X*dE - XI*dE3 + XII*dE5 - XIIA*dE7;


  return {

    longitude: lon.toDeg(),
    latitude: lat.toDeg()

  };

}

Number.prototype.toRad = function() {  // convert degrees to radians
  return this * Math.PI / 180;
}
Number.prototype.toDeg = function() {  // convert radians to degrees (signed)
  return this * 180 / Math.PI;
}
Number.prototype.padLZ = function(w) {
  var n = this.toString();
  for (var i=0; i<w-n.length; i++) n = '0' + n;
  return n;
}

Marinamarinade

Answer 2

I've been struggling with this one for a while. I had a lot of northing/easting points in OSGB36 that have to be converted on the fly on a regular basis. Please note that the UDF below converts northings/eastings in OSGB36 (Ordnance Survey) projection to latitude/longitude in WGS84 projection so they can be used in Google Maps.

/****** Object:  UserDefinedFunction [dbo].[NEtoLL]    Script Date: 09/06/2012 17:06:39 ******/
SET ANSI_NULLS ON
GO

SET QUOTED_IDENTIFIER ON
GO

CREATE FUNCTION [dbo].[NEtoLL] (@East INT, @North INT, @LatOrLng VARCHAR(3)) RETURNS FLOAT AS
BEGIN

--Author: Sandy Motteram
--Date:   06 September 2012

--UDF adapted from javascript at http://www.bdcc.co.uk/LatLngToOSGB.js
--found on page http://mapki.com/wiki/Tools:Snippets

--Instructions:
--Latitude and Longitude are calculated based on BOTH the easting and northing values from the OSGB36
--This UDF takes both easting and northing values in OSGB36 projection and you must specify if a latitude or longitude co-ordinate should be returned.
--IT first converts E/N values to lat and long in OSGB36 projection, then converts those values to lat/lng in WGS84 projection

--Sample values below
--DECLARE @East INT, @North INT, @LatOrLng VARCHAR(3)
--SELECT @East = 529000, @North = 183650 --that combo should be the corner of Camden High St and Delancey St


    DECLARE @Pi              FLOAT
          , @K0              FLOAT
          , @OriginLat       FLOAT
          , @OriginLong      FLOAT
          , @OriginX         FLOAT
          , @OriginY         FLOAT
          , @a               FLOAT
          , @b               FLOAT
          , @e2              FLOAT
          , @ex              FLOAT
          , @n1              FLOAT
          , @n2              FLOAT
          , @n3              FLOAT
          , @OriginNorthings FLOAT
          , @lat             FLOAT
          , @lon             FLOAT
          , @Northing        FLOAT
          , @Easting         FLOAT

    SELECT  @Pi = 3.14159265358979323846
          , @K0 = 0.9996012717 -- grid scale factor on central meridean
          , @OriginLat  = 49.0
          , @OriginLong = -2.0
          , @OriginX =  400000 -- 400 kM
          , @OriginY = -100000 -- 100 kM
          , @a = 6377563.396   -- Airy Spheroid
          , @b = 6356256.910
    /*    , @e2
          , @ex
          , @n1
          , @n2
          , @n3
          , @OriginNorthings*/

    -- compute interim values
    SELECT  @a = @a * @K0
          , @b = @b * @K0

    SET     @n1 = (@a - @b) / (@a + @b)
    SET     @n2 = @n1 * @n1
    SET     @n3 = @n2 * @n1

    SET     @lat = @OriginLat * @Pi / 180.0 -- to radians

    SELECT  @e2 = (@a * @a - @b * @b) / (@a * @a) -- first eccentricity
          , @ex = (@a * @a - @b * @b) / (@b * @b) -- second eccentricity

    SET     @OriginNorthings = @b * @lat + @b * (@n1 * (1.0 + 5.0 * @n1 * (1.0 + @n1) / 4.0) * @lat
          - 3.0 * @n1 * (1.0 + @n1 * (1.0 + 7.0 * @n1 / 8.0)) * SIN(@lat) * COS(@lat)
          + (15.0 * @n1 * (@n1 + @n2) / 8.0) * SIN(2.0 * @lat) * COS(2.0 * @lat)
          - (35.0 * @n3 / 24.0) * SIN(3.0 * @lat) * COS(3.0 * @lat))

    SELECT  @northing = @north - @OriginY
         ,  @easting  = @east  - @OriginX

    DECLARE @nu       FLOAT
          , @phid     FLOAT
          , @phid2    FLOAT
          , @t2       FLOAT
          , @t        FLOAT
          , @q2       FLOAT
          , @c        FLOAT
          , @s        FLOAT
          , @nphid    FLOAT
          , @dnphid   FLOAT
          , @nu2      FLOAT
          , @nudivrho FLOAT
          , @invnurho FLOAT
          , @rho      FLOAT
          , @eta2     FLOAT

    /* Evaluate M term: latitude of the northing on the centre meridian */

    SET     @northing = @northing + @OriginNorthings

    SET     @phid  = @northing / (@b*(1.0 + @n1 + 5.0 * (@n2 + @n3) / 4.0)) - 1.0
    SET     @phid2 = @phid + 1.0

    WHILE (ABS(@phid2 - @phid) > 0.000001)
    BEGIN
        SET @phid = @phid2;
        SET @nphid = @b * @phid + @b * (@n1 * (1.0 + 5.0 * @n1 * (1.0 + @n1) / 4.0) * @phid
                   - 3.0 * @n1 * (1.0 + @n1 * (1.0 + 7.0 * @n1 / 8.0)) * SIN(@phid) * COS(@phid)
                   + (15.0 * @n1 * (@n1 + @n2) / 8.0) * SIN(2.0 * @phid) * COS(2.0 * @phid)
                   - (35.0 * @n3 / 24.0) * SIN(3.0 * @phid) * COS(3.0 * @phid))

        SET @dnphid = @b * ((1.0 + @n1 + 5.0 * (@n2 + @n3) / 4.0) - 3.0 * (@n1 + @n2 + 7.0 * @n3 / 8.0) * COS(2.0 * @phid)
                    + (15.0 * (@n2 + @n3) / 4.0) * COS(4 * @phid) - (35.0 * @n3 / 8.0) * COS(6.0 * @phid))

        SET @phid2 = @phid - (@nphid - @northing) / @dnphid
    END

    SELECT @c = COS(@phid)
         , @s = SIN(@phid)
         , @t = TAN(@phid)
    SELECT @t2 = @t * @t
         , @q2 = @easting * @easting

    SET    @nu2 = (@a * @a) / (1.0 - @e2 * @s * @s)
    SET    @nu = SQRT(@nu2)

    SET    @nudivrho = @a * @a * @c * @c / (@b * @b) - @c * @c + 1.0
    SET    @eta2 = @nudivrho - 1
    SET    @rho = @nu / @nudivrho;

    SET    @invnurho = ((1.0 - @e2 * @s * @s) * (1.0 - @e2 * @s * @s)) / (@a * @a * (1.0 - @e2))

    SET    @lat = @phid - @t * @q2 * @invnurho / 2.0 + (@q2 * @q2 * (@t / (24 * @rho * @nu2 * @nu) * (5 + (3 * @t2) + @eta2 - (9 * @t2 * @eta2))))
    SET    @lon = (@easting / (@c * @nu))
                - (@easting * @q2 * ((@nudivrho + 2.0 * @t2) / (6.0 * @nu2)) / (@c * @nu))
                + (@q2 * @q2 * @easting * (5 + (28 * @t2) + (24 * @t2 * @t2)) / (120 * @nu2 * @nu2 * @nu * @c))


    SELECT @lat = @lat * 180.0 / @Pi
         , @lon = @lon * 180.0 / @Pi + @OriginLong


--Now convert the lat and long from OSGB36 to WGS84

    DECLARE @OGlat  FLOAT
          , @OGlon  FLOAT
          , @height FLOAT

    SELECT  @OGlat  = @lat
          , @OGlon  = @lon
          , @height = 24 --London's mean height above sea level is 24 metres. Adjust for other locations.

    DECLARE @deg2rad  FLOAT
          , @rad2deg  FLOAT
          , @radOGlat FLOAT
          , @radOGlon FLOAT

    SELECT  @deg2rad = @Pi / 180
          , @rad2deg = 180 / @Pi

    --first off convert to radians
    SELECT  @radOGlat = @OGlat * @deg2rad
          , @radOGlon = @OGlon * @deg2rad
    --these are the values for WGS84(GRS80) to OSGB36(Airy) 

    DECLARE @a2       FLOAT
          , @h        FLOAT
          , @xp       FLOAT
          , @yp       FLOAT
          , @zp       FLOAT
          , @xr       FLOAT
          , @yr       FLOAT
          , @zr       FLOAT
          , @sf       FLOAT
          , @e        FLOAT
          , @v        FLOAT
          , @x        FLOAT
          , @y        FLOAT
          , @z        FLOAT
          , @xrot     FLOAT
          , @yrot     FLOAT
          , @zrot     FLOAT
          , @hx       FLOAT
          , @hy       FLOAT
          , @hz       FLOAT
          , @newLon   FLOAT
          , @newLat   FLOAT
          , @p        FLOAT
          , @errvalue FLOAT
          , @lat0     FLOAT

    SELECT  @a2 = 6378137             -- WGS84_AXIS
          , @e2 = 0.00669438037928458 -- WGS84_ECCENTRIC
          , @h  = @height             -- height above datum (from $GPGGA sentence)
          , @a  = 6377563.396         -- OSGB_AXIS
          , @e  = 0.0066705397616     -- OSGB_ECCENTRIC
          , @xp = 446.448
          , @yp = -125.157
          , @zp = 542.06
          , @xr = 0.1502
          , @yr = 0.247
          , @zr = 0.8421
          , @s  = -20.4894

    -- convert to cartesian; lat, lon are in radians
    SET @sf = @s * 0.000001
    SET @v = @a / (sqrt(1 - (@e * (SIN(@radOGlat) * SIN(@radOGlat)))))
    SET @x = (@v + @h) * COS(@radOGlat) * COS(@radOGlon)
    SET @y = (@v + @h) * COS(@radOGlat) * SIN(@radOGlon)
    SET @z = ((1 - @e) * @v + @h) * SIN(@radOGlat)

    -- transform cartesian
    SET @xrot = (@xr / 3600) * @deg2rad
    SET @yrot = (@yr / 3600) * @deg2rad
    SET @zrot = (@zr / 3600) * @deg2rad
    SET @hx = @x + (@x * @sf) - (@y * @zrot) + (@z * @yrot) + @xp
    SET @hy = (@x * @zrot) + @y + (@y * @sf) - (@z * @xrot) + @yp
    SET @hz = (-1 * @x * @yrot) + (@y * @xrot) + @z + (@z * @sf) + @zp

    -- Convert back to lat, lon
    SET @newLon = ATAN(@hy / @hx)
    SET @p = SQRT((@hx * @hx) + (@hy * @hy))
    SET @newLat = ATAN(@hz / (@p * (1 - @e2)))
    SET @v = @a2 / (SQRT(1 - @e2 * (SIN(@newLat) * SIN(@newLat))))
    SET @errvalue = 1.0;
    SET @lat0 = 0
    WHILE (@errvalue > 0.001)
    BEGIN
        SET @lat0 = ATAN((@hz + @e2 * @v * SIN(@newLat)) / @p)
        SET @errvalue = ABS(@lat0 - @newLat)
        SET @newLat = @lat0
    END

    --convert back to degrees
    SET @newLat = @newLat * @rad2deg
    SET @newLon = @newLon * @rad2deg

    DECLARE @ReturnMe FLOAT
    SET @ReturnMe = 0

    IF @LatOrLng = 'Lat'
        SET @ReturnMe = @newLat
    IF @LatOrLng = 'Lng'
        SET @ReturnMe = @newLon

    RETURN @ReturnMe
END
GO

Umbelliferous

Thank you! Like CL4NCY I have been looking for this function. Does the job perfectly. – Vegetable

@PeterMcMinn Hi Peter / SiO2y, Have you managed to test this SQL version in a production enviroment, is it accurate? – Ma

Bearing in mind the data I am working from is not that accurate to begin with, it seems to get within 200-300 meters of the actual point. This is good enough for what we need. – Vegetable

For any future visitors to this question, I've done some empirical tests and compared Niall's PHP version, and the OS OpenSpace Library - and there seems to be ~100 metre difference in the resulting location (but YMMV, depending on the accuracy of your original data, and what you want to do with it...) – Leialeibman

Doesn't appear to be accurate. According to – Accusatory

Sorry all. I never received any emails to say people had commented on this post. I found the solution to be extremely accurate when I used it in our production environment back in 2012. I wonder if any inaccuracies you've encountered are due to different projections being used? – Umbelliferous

You can avoid running this twice by returning SQL's built in type: return geography::Point(@newLat, @newLon, 4326) – Contrasty

The inaccuracy might be because it's providing OSGB36 latitude and longitude, rather than WSG84. There is a descrepancy of about 100-150m between them. – Village

Answer 3

I needed the same function, and javascript made it difficult to interact with the DB. I have converted your JS to PHP and this could be more useful when updating your database - ie: query table, loop through result set, call function, update table.

function OSGridToLatLong($E, $N) {

  $a = 6377563.396; 
  $b = 6356256.910;              // Airy 1830 major & minor semi-axes
  $F0 = 0.9996012717;                             // NatGrid scale factor on central meridian
  $lat0 = 49*M_PI/180; 
  $lon0 = -2*M_PI/180;  // NatGrid true origin
  $N0 = -100000;
  $E0 = 400000;                     // northing & easting of true origin, metres
  $e2 = 1 - ($b*$b)/($a*$a);                          // eccentricity squared
  $n = ($a-$b)/($a+$b);
  $n2 = $n*$n;
  $n3 = $n*$n*$n;

  $lat=$lat0;
  $M=0;
  do {
    $lat = ($N-$N0-$M)/($a*$F0) + $lat;

    $Ma = (1 + $n + (5/4)*$n2 + (5/4)*$n3) * ($lat-$lat0);
    $Mb = (3*$n + 3*$n*$n + (21/8)*$n3) * sin($lat-$lat0) * cos($lat+$lat0);
    $Mc = ((15/8)*$n2 + (15/8)*$n3) * sin(2*($lat-$lat0)) * cos(2*($lat+$lat0));
    $Md = (35/24)*$n3 * sin(3*($lat-$lat0)) * cos(3*($lat+$lat0));
    $M = $b * $F0 * ($Ma - $Mb + $Mc - $Md);                // meridional arc

  } while ($N-$N0-$M >= 0.00001);  // ie until < 0.01mm

  $cosLat = cos($lat);
  $sinLat = sin($lat);
  $nu = $a*$F0/sqrt(1-$e2*$sinLat*$sinLat);              // transverse radius of curvature
  $rho = $a*$F0*(1-$e2)/pow(1-$e2*$sinLat*$sinLat, 1.5);  // meridional radius of curvature
  $eta2 = $nu/$rho-1;

  $tanLat = tan($lat);
  $tan2lat = $tanLat*$tanLat;
  $tan4lat = $tan2lat*$tan2lat;
  $tan6lat = $tan4lat*$tan2lat;
  $secLat = 1/$cosLat;
  $nu3 = $nu*$nu*$nu;
  $nu5 = $nu3*$nu*$nu;
  $nu7 = $nu5*$nu*$nu;
  $VII = $tanLat/(2*$rho*$nu);
  $VIII = $tanLat/(24*$rho*$nu3)*(5+3*$tan2lat+$eta2-9*$tan2lat*$eta2);
  $IX = $tanLat/(720*$rho*$nu5)*(61+90*$tan2lat+45*$tan4lat);
  $X = $secLat/$nu;
  $XI = $secLat/(6*$nu3)*($nu/$rho+2*$tan2lat);
  $XII = $secLat/(120*$nu5)*(5+28*$tan2lat+24*$tan4lat);
  $XIIA = $secLat/(5040*$nu7)*(61+662*$tan2lat+1320*$tan4lat+720*$tan6lat);

  $dE = ($E-$E0);
  $dE2 = $dE*$dE;
  $dE3 = $dE2*$dE;
  $dE4 = $dE2*$dE2;
  $dE5 = $dE3*$dE2;
  $dE6 = $dE4*$dE2;
  $dE7 = $dE5*$dE2;
  $lat = $lat - $VII*$dE2 + $VIII*$dE4 - $IX*$dE6;
  $lon = $lon0 + $X*$dE - $XI*$dE3 + $XII*$dE5 - $XIIA*$dE7;


  return array(

    'longitude' => $lon * 180 / M_PI,
    'latitude'  => $lat * 180 / M_PI

  );

}

Teece

I had issues with the SQL UDF above (100% user error, not knowing how to add/use it). This PHP function worked perfectly, coords seem accurate when checking thru GMaps. Thanks! – Foreworn

Hoping this works, will try it tomorrow. Had extreme difficulty finding details of the transformation online, mostly due to useless search engines only finding online tools and websites which waffled on without providing any useful information. Even OS's own website only provides an online tool, but no details of the algorithm. This is not proprietary information. It should be much easier to find. – Village

Couldn't wait. Tested. Works like a dream. If I could give you a thousand upvotes, I would. Could possibly cache results, as I see there is a narrowing-down loop, but that would be up to the implementor. Thank you so much. – Village

Actually not quite like a dream. It provides OSGB36 latitude and longitude, which differ from WSG84 by around 100-150m. A final transform is needed to get WSG84 co-ordinates, which I'll have to post as a separate answer. – Village

Answer 4

If anyone's interested in non-SQL solution I strongly recommend using this http://www.howtocreate.co.uk/php/gridref.php PHP/JavaScript class.

One important thing to mention here is the library supports Helmert transformation.

PHP

$grutoolbox = Grid_Ref_Utils::toolbox();
$source_coords = Array(54.607720,-6.411990);
//get the ellipsoids that will be used
$Airy_1830 = $grutoolbox->get_ellipsoid('Airy_1830');
$WGS84 = $grutoolbox->get_ellipsoid('WGS84');
$Airy_1830_mod = $grutoolbox->get_ellipsoid('Airy_1830_mod');
//get the transform parameters that will be used
$UK_to_GPS = $grutoolbox->get_transformation('OSGB36_to_WGS84');
$GPS_to_Ireland = $grutoolbox->get_transformation('WGS84_to_Ireland65');
//convert to GPS coordinates
$gps_coords = $grutoolbox->Helmert_transform($source_coords,$Airy_1830,$UK_to_GPS,$WGS84);
//convert to destination coordinates
print $grutoolbox->Helmert_transform($source_coords,$WGS84,$GPS_to_Ireland,$Airy_1830_mod,$grutoolbox->HTML);

JavaScript

var grutoolbox = gridRefUtilsToolbox();
var sourceCoords = [54.607720,-6.411990];
//get the ellipsoids that will be used
var Airy1830 = grutoolbox.getEllipsoid('Airy_1830');
var WGS84 = grutoolbox.getEllipsoid('WGS84');
var Airy1830Mod = grutoolbox.getEllipsoid('Airy_1830_mod');
//get the transform parameters that will be used
var UKToGPS = grutoolbox.getTransformation('OSGB36_to_WGS84');
var GPSToIreland = grutoolbox.getTransformation('WGS84_to_Ireland65');
//convert to GPS coordinates
var gpsCoords = grutoolbox.HelmertTransform(sourceCoords,Airy1830,UKToGPS,WGS84);
//convert to destination coordinates
element.innerHTML = grutoolbox.HelmertTransform(sourceCoords,WGS84,GPSToIreland,Airy1830Mod,grutoolbox.HTML);

Tenuis

Answer 5

I have developed a library in .NET to be called from transact sql Converts WGS84/UTM coordinates to Latitude and Longitude

It does just the opposite, but as it uses CoordinateSharp you can download the code and change it easily to convert from lat/long to wgs84 instead.

You can download it from github:

https://github.com/j-v-garcia/UTM2LATITUDE

  usage:

    SELECT dbo.UTM2LATITUDE(723399.51,4373328.5,'S',30) AS Latitude, dbo.UTM2LONGITUDE(723399.51,4373328.5,'S',30) AS Longitude

    result: 

        39,4805657453054    -0,402592727245112

        <param name="XUTM">pos UTM X</param>
        <param name="YUTM">pos UTM Y</param>
        <param name="LatBand">Latitude band grid zone designation letter (see http://www.dmap.co.uk/utmworld.htm) </param>
        <param name="LongBand">Longitude band grid zone designation number (see http://www.dmap.co.uk/utmworld.htm) </param>

Keirakeiser

Please don't add the same answer to multiple questions. Answer the best one and flag the rest as duplicates. See Is it acceptable to add a duplicate answer to several questions? – Groningen

Answer 6

Building upon @Niall's answer, this PHP function also includes the final transform from OSGB36 to WSG84.

The current top 3 answers (SQL, JavaScript and PHP) all use the same algorithm (I think) which produces OSGB36 latitude and longitude. That is an 1836 standard, based on the understanding of the curvature of the Earth at the time. WSG84 is a 1984 standard, in which the location pinpointed by the co-ordinates may differ by up to about 150m.

Source: http://www.movable-type.co.uk/scripts/latlong-os-gridref.html

function OSGridToLatLong($E, $N)
{
  $a = 6377563.396;
  $b = 6356256.910;     // Airy 1830 major & minor semi-axes
  $F0 = 0.9996012717;   // NatGrid scale factor on central meridian
  $lat0 = 49*M_PI/180;
  $lon0 = -2*M_PI/180;  // NatGrid true origin
  $N0 = -100000;
  $E0 = 400000;         // northing & easting of true origin, metres
  $e2 = 1 - ($b*$b)/($a*$a);  // eccentricity squared
  $n = ($a-$b)/($a+$b);
  $n2 = $n*$n;
  $n3 = $n*$n*$n;

  $lat=$lat0;
  $M=0;
  do {
    $lat = ($N-$N0-$M)/($a*$F0) + $lat;

    $Ma = (1 + $n + (5/4)*$n2 + (5/4)*$n3) * ($lat-$lat0);
    $Mb = (3*$n + 3*$n*$n + (21/8)*$n3) * sin($lat-$lat0) * cos($lat+$lat0);
    $Mc = ((15/8)*$n2 + (15/8)*$n3) * sin(2*($lat-$lat0)) * cos(2*($lat+$lat0));
    $Md = (35/24)*$n3 * sin(3*($lat-$lat0)) * cos(3*($lat+$lat0));
    $M = $b * $F0 * ($Ma - $Mb + $Mc - $Md);  // meridional arc

  } while ($N-$N0-$M >= 0.00001); // ie until < 0.01mm

  $cosLat = cos($lat);
  $sinLat = sin($lat);
  $nu = $a*$F0/sqrt(1-$e2*$sinLat*$sinLat);               // transverse radius of curvature
  $rho = $a*$F0*(1-$e2)/pow(1-$e2*$sinLat*$sinLat, 1.5);  // meridional radius of curvature
  $eta2 = $nu/$rho-1;

  $tanLat = tan($lat);
  $tan2lat = $tanLat*$tanLat;
  $tan4lat = $tan2lat*$tan2lat;
  $tan6lat = $tan4lat*$tan2lat;
  $secLat = 1/$cosLat;
  $nu3 = $nu*$nu*$nu;
  $nu5 = $nu3*$nu*$nu;
  $nu7 = $nu5*$nu*$nu;
  $VII = $tanLat/(2*$rho*$nu);
  $VIII = $tanLat/(24*$rho*$nu3)*(5+3*$tan2lat+$eta2-9*$tan2lat*$eta2);
  $IX = $tanLat/(720*$rho*$nu5)*(61+90*$tan2lat+45*$tan4lat);
  $X = $secLat/$nu;
  $XI = $secLat/(6*$nu3)*($nu/$rho+2*$tan2lat);
  $XII = $secLat/(120*$nu5)*(5+28*$tan2lat+24*$tan4lat);
  $XIIA = $secLat/(5040*$nu7)*(61+662*$tan2lat+1320*$tan4lat+720*$tan6lat);

  $dE = ($E-$E0);
  $dE2 = $dE*$dE;
  $dE3 = $dE2*$dE;
  $dE4 = $dE2*$dE2;
  $dE5 = $dE3*$dE2;
  $dE6 = $dE4*$dE2;
  $dE7 = $dE5*$dE2;
  $lat = $lat - $VII*$dE2 + $VIII*$dE4 - $IX*$dE6;
  $lon = $lon0 + $X*$dE - $XI*$dE3 + $XII*$dE5 - $XIIA*$dE7;

  $f = 1/299.3249646; // Airy 1830
  // Convert to cartesian
  $h = 0;
  $sinLat = \sin($lat);
  $cosLat = \cos($lat);
  $sinLon = \sin($lon);
  $cosLon = \cos($lon);
  $eSq = 2 * $f - $f * $f; // 1st eccentricity squared ≡ (a²-b²)/a²
  $v = $a / \sqrt(1 - $eSq * $sinLat * $sinLat); // radius of curvature in prime vertical
  $x = ($v + $h) * $cosLat * $cosLon;
  $y = ($v + $h) * $cosLat * $sinLon;
  $z = ($v * (1 - $eSq) + $h) * $sinLat;

  // Transform from OSGB36 to WSG84
  $t = [446.448, -125.157, 542.060, -20.4894, 0.1502, 0.2470, 0.8421];
  $tx = $t[0];                       // x-shift in metres
  $ty = $t[1];                       // y-shift in metres
  $tz = $t[2];                       // z-shift in metres
  $s  = $t[3] / 1e6 + 1;             // scale: normalise parts-per-million to (s+1)
  $rx = ($t[4] / 3600) * M_PI / 180; // x-rotation: normalise arcseconds to radians
  $ry = ($t[5] / 3600) * M_PI / 180; // y-rotation: normalise arcseconds to radians
  $rz = ($t[6] / 3600) * M_PI / 180; // z-rotation: normalise arcseconds to radians
  $x2 = $tx + $x * $s - $y * $rz + $z * $ry;
  $y2 = $ty + $x * $rz + $y * $s - $z * $rx;
  $z2 = $tz - $x * $ry + $y * $rx + $z * $s;

  // Convert back to WSG84 latitude and longitude
  $a = 6378137;
  $b = 6356752.314245;
  $f = 1/298.257223563;

  $e2 = 2 * $f - $f * $f;      // 1st eccentricity squared ≡ (a²−b²)/a²
  $epsilon2 = $e2 / (1 - $e2); // 2nd eccentricity squared ≡ (a²−b²)/b²
  $p = \sqrt($x2 * $x2 + $y2 * $y2); // distance from minor axis
  $R = \sqrt($p * $p + $z2 * $z2); // polar radius
  // parametric latitude (Bowring eqn.17, replacing tanβ = z·a / p·b)
  $tanBeta = ($b * $z2) / ($a * $p) * (1 + $epsilon2 * $b / $R);
  $sinBeta = $tanBeta / \sqrt(1 + $tanBeta * $tanBeta);
  $cosBeta = $sinBeta / $tanBeta;
  // geodetic latitude (Bowring eqn.18: tanφ = z+epsilon²⋅b⋅sin³β / p−e²⋅cos³β)
  $lat = \is_nan($cosBeta) ? 0 :
    \atan2(
      $z2 + $epsilon2 * $b * $sinBeta * $sinBeta * $sinBeta,
      $p - $e2 * $a * $cosBeta * $cosBeta * $cosBeta
    );
  // longitude
  $lon = \atan2($y2, $x2);
  // height above ellipsoid (Bowring eqn.7)
  $sinLat = \sin($lat);
  $cosLat = \cos($lat);
  $v = $a / \sqrt(1 - $e2 * $sinLat * $sinLat); // length of the normal terminated by the minor axis
  $h = $p * $cosLat + $z2 * $sinLat - ($a * $a / $v);

  return array(
    'longitude' => $lon * 180 / M_PI,
    'latitude'  => $lat * 180 / M_PI,
    'height' => $h,
  );
}

Village

Answer 7

Building on the previous answer (thank you!) I had to do this in Google Bigquery.

I used the following to run it in BQ with some SQL on the end of a temp function:

CREATE TEMP FUNCTION OSGridToLatLong(E FLOAT64, N FLOAT64)
RETURNS STRUCT<latitude FLOAT64, longitude FLOAT64, height FLOAT64>
LANGUAGE js AS """
  var a = 6377563.396;
  var b = 6356256.910;     
  var F0 = 0.9996012717;   
  var lat0 = 49 * Math.PI / 180;
  var lon0 = -2 * Math.PI / 180;  
  var N0 = -100000;
  var E0 = 400000;         
  var e2 = 1 - (b*b)/(a*a);  
  var n = (a-b)/(a+b);
  var n2 = n*n;
  var n3 = n*n*n;

  var lat = lat0;
  var M = 0;
  do {
    lat = (N - N0 - M) / (a * F0) + lat;

    var Ma = (1 + n + (5/4)*n2 + (5/4)*n3) * (lat - lat0);
    var Mb = (3*n + 3*n*n + (21/8)*n3) * Math.sin(lat - lat0) * Math.cos(lat + lat0);
    var Mc = ((15/8)*n2 + (15/8)*n3) * Math.sin(2*(lat - lat0)) * Math.cos(2*(lat + lat0));
    var Md = (35/24)*n3 * Math.sin(3*(lat - lat0)) * Math.cos(3*(lat + lat0));
    M = b * F0 * (Ma - Mb + Mc - Md);  

  } while (N - N0 - M >= 0.00001);

  var cosLat = Math.cos(lat);
  var sinLat = Math.sin(lat);
  var nu = a * F0 / Math.sqrt(1 - e2 * sinLat * sinLat);               
  var rho = a * F0 * (1 - e2) / Math.pow(1 - e2 * sinLat * sinLat, 1.5);  
  var eta2 = nu / rho - 1;

  var tanLat = Math.tan(lat);
  var tan2lat = tanLat * tanLat;
  var tan4lat = tan2lat * tan2lat;
  var tan6lat = tan4lat * tan2lat;
  var secLat = 1 / cosLat;
  var nu3 = nu * nu * nu;
  var nu5 = nu3 * nu * nu;
  var nu7 = nu5 * nu * nu;
  var VII = tanLat / (2 * rho * nu);
  var VIII = tanLat / (24 * rho * nu3) * (5 + 3 * tan2lat + eta2 - 9 * tan2lat * eta2);
  var IX = tanLat / (720 * rho * nu5) * (61 + 90 * tan2lat + 45 * tan4lat);
  var X = secLat / nu;
  var XI = secLat / (6 * nu3) * (nu / rho + 2 * tan2lat);
  var XII = secLat / (120 * nu5) * (5 + 28 * tan2lat + 24 * tan4lat);
  var XIIA = secLat / (5040 * nu7) * (61 + 662 * tan2lat + 1320 * tan4lat + 720 * tan6lat);

  var dE = (E - E0);
  var dE2 = dE * dE;
  var dE3 = dE2 * dE;
  var dE4 = dE2 * dE2;
  var dE5 = dE3 * dE2;
  var dE6 = dE4 * dE2;
  var dE7 = dE5 * dE2;
  lat = lat - VII * dE2 + VIII * dE4 - IX * dE6;
  var lon = lon0 + X * dE - XI * dE3 + XII * dE5 - XIIA * dE7;

  var f = 1 / 299.3249646; 
  var h = 0;
  sinLat = Math.sin(lat);
  cosLat = Math.cos(lat);
  var sinLon = Math.sin(lon);
  var cosLon = Math.cos(lon);
  var eSq = 2 * f - f * f;
  var v = a / Math.sqrt(1 - eSq * sinLat * sinLat);
  var x = (v + h) * cosLat * cosLon;
  var y = (v + h) * cosLat * sinLon;
  var z = (v * (1 - eSq) + h) * sinLat;

  var t = [446.448, -125.157, 542.060, -20.4894, 0.1502, 0.2470, 0.8421];
  var tx = t[0];                      
  var ty = t[1];                       
  var tz = t[2];                       
  var s  = t[3] / 1e6 + 1;             
  var rx = (t[4] / 3600) * Math.PI / 180; 
  var ry = (t[5] / 3600) * Math.PI / 180; 
  var rz = (t[6] / 3600) * Math.PI / 180; 
  var x2 = tx + x * s - y * rz + z * ry;
  var y2 = ty + x * rz + y * s - z * rx;
  var z2 = tz - x * ry + y * rx + z * s;

  a = 6378137;
  b = 6356752.314245;
  f = 1 / 298.257223563;

  e2 = 2 * f - f * f;
  var epsilon2 = e2 / (1 - e2);
  var p = Math.sqrt(x2 * x2 + y2 * y2);
  var R = Math.sqrt(p * p + z2 * z2);
  var tanBeta = (b * z2) / (a * p) * (1 + epsilon2 * b / R);
  var sinBeta = tanBeta / Math.sqrt(1 + tanBeta * tanBeta);
  var cosBeta = sinBeta / tanBeta;
  lat = isNaN(cosBeta) ? 0 :
    Math.atan2(
      z2 + epsilon2 * b * sinBeta * sinBeta * sinBeta,
      p - e2 * a * cosBeta * cosBeta * cosBeta
    );
  lon = Math.atan2(y2, x2);
  sinLat = Math.sin(lat);
  cosLat = Math.cos(lat);
  v = a / Math.sqrt(1 - e2 * sinLat * sinLat);
  h = p * cosLat + z2 * sinLat - (a * a / v);

  return {latitude: lat * 180 / Math.PI, longitude: lon * 180 / Math.PI, height: h};
""";

WITH input_data AS (
  SELECT 
    574807.122806563 AS E, --Add you data here, can also be a whole table to convert
    221464.559972554 AS N
)

SELECT
  OSGridToLatLong(E, N).latitude AS latitude,
  OSGridToLatLong(E, N).longitude AS longitude

FROM
  input_data

Naphthalene