SQL function to convert UK OS coordinates from east to north and longitude and latitude

Question

Please someone can send an SQL function to convert east / north to longitude / latitude. I know this is incredibly complicated, but I did not find anyone who registered it in T-SQL.

This javascript code works, but I am having problems converting it to SQL.

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

This is what I have so far, but it does not go 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 

Answer 1

As a result, I used the following javascript functions to convert the values. I know this is not a SQL solution, but he 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 = (ab)/(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<wn.length; i++) n = '0' + n; return n; } 

Answer 2

I struggled with this for a while. In OSGB36, I had a lot of north / east points that needed to be converted on the fly on a regular basis. Please note that the UDF below converts north / east directions in OSGB36 (Ordinance Poll) to latitude / longitude in the WGS84 projection so that 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 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 

Answer 3

I needed the same function, and javascript made it difficult to interact with the database. I converted your JS to PHP, and this can be more useful when updating the database, i.e.: Query table, loop through the 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 ); } 

Answer 4

If anyone is interested in a solution other than SQL, I highly recommend using this http://www.howtocreate.co.uk/php/gridref.php PHP / JavaScript class.

It is very important to note here that the library supports Helmert transform .

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