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 = 213371DECLARE @a real = 6377563.396,@b real = 6356256.910,@F0 real = 0.9996012717,@lat0 real = 49*PI()/180,@lon0 real = -2*PI()/180DECLARE @N0 real = -100000,@E0 real = 400000,@e2 real = 1 - (@b*@b)/(@a*@a),@n real = (@[email protected])/(@[email protected])DECLARE @n2 real = @n*@n,@n3 real = @n*@n*@nDECLARE @lat real = @lat0,@M real = 0WHILE (@[email protected]@M >= 0.00001)BEGINSET @lat = ((@[email protected]@M)/(@a*@F0)) + @latDECLARE @Ma real = (1 + @n + (5/4)*@n2 + (5/4)*@n3) * (@[email protected]),@Mb real = (3*@n + 3*@n*@n + (21/8)*@n3) * SIN(@[email protected]) * COS(@[email protected]),@Mc real = ((15/8)*@n2 + (15/8)*@n3) * SIN(2*(@[email protected])) * COS(2*(@[email protected])),@Md real = (35/24)*@n3 * SIN(3*(@[email protected])) * COS(3*(@[email protected]))SET @M = @b * @F0 * (@Ma - @Mb + @Mc - @Md)ENDDECLARE @cosLat real = COS(@lat),@sinLat real = SIN(@lat)DECLARE @nu real = @a*@F0/sqrt([email protected]*@sinLat*@sinLat)DECLARE @rho real = @a*@F0*([email protected])/POWER([email protected]*@sinLat*@sinLat, 1.5)DECLARE @eta2 real = @nu/@rho-1DECLARE @tanLat real = tan(@lat)DECLARE @tan2lat real = @tanLat*@tanLatDECLARE @tan4lat real = @tan2lat*@tan2latDECLARE @tan6lat real = @tan4lat*@tan2latDECLARE @secLat real = 1/@cosLatDECLARE @nu3 real = @nu*@nu*@nuDECLARE @nu5 real = @nu3*@nu*@nuDECLARE @nu7 real = @nu5*@nu*@nuDECLARE @VII real = @tanLat/(2*@rho*@nu)DECLARE @VIII real = @tanLat/(24*@rho*@nu3)*(5+3*@[email protected]*@tan2lat*@eta2)DECLARE @IX real = @tanLat/(720*@rho*@nu5)*(61+90*@tan2lat+45*@tan4lat)DECLARE @X real = @secLat/@nuDECLARE @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 = (@[email protected])DECLARE @dE2 real = @dE*@dEDECLARE @dE3 real = @dE2*@dEDECLARE @dE4 real = @dE2*@dE2,@dE5 real = @dE3*@dE2DECLARE @dE6 real = @dE4*@dE2,@dE7 real = @dE5*@dE2SET @lat = @lat - @VII*@dE2 + @VIII*@dE4 - @IX*@dE6DECLARE @lon real = @lon0 + @X*@dE - @XI*@dE3 + @XII*@dE5 - @XIIA*@dE7SELECT @lon, @lat``

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.

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

--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 @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
``````

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()

};

}

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

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

);

}
``````

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);
``````

Top