udf_OrdnanceSurvey_WGS84.sql

create function dbo.udf_OrdnanceSurvey_WGS84 (@east int, @north int)
returns geography
as
begin
    if @east is null
        or @north is null
		or (@east = 0 and @north = 0) -- 0,0 is a legitmate co-ordinate, but in our DB it's used instead of missing locations
        return null

	-- Adapted from https://stackoverflow.com/a/12317217/905
	-- JS reference: http://www.movable-type.co.uk/scripts/latlong-os-gridref.html
	-- See also: https://osedok.wordpress.com/2012/01/17/conversion-of-british-national-grid-wkid27700-to-wgs84wkid4326-and-then-to-webmercator-wkid102100/

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

    declare @π 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  @π = 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

		-- 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 * @π / 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 / @π, @lon = @lon * 180.0 / @π + @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 = @π / 180, @rad2deg = 180 / @π

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

    return geography::Point(@newLat, @newLon, 4326)
end