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