My ultimate goal:
Write a TSQL function to convert NAD83 State Plane (specifically ESRI Projection 102738 http://spatialreference.org/ref/esri/102738/) to WGS84 Lat/Lng so I can, in near real-time, plot the original coordinates on a Google Map.
I do not want to use an external library such as Proj4 or GDAL because I want to avoid external, non-sql dependencies that I am afraid won't get reinstalled by future admins when our database servers are upgraded.
I found javascript code (http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html) to convert NAD83 UTM to WGS84 Lat/Long (see below). I have ported/converted this to SQL successfully.
I don't know how to find the math to write a function that can convert NAD83 State Plane coordinates to NAD 83 UTM (so I can then use my existing function to convert UTM to Lat/Lng)
TSQL I've ported from JS to convert NAD83 UTM to WGS84:
/*Adapted from: http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html*/
-----------------------------------------
--USER INPUTS
DECLARE @Operation VARCHAR(50) = 'FromUTMToGeo'
DECLARE @Input_X_or_Lat FLOAT = 686847.2227879118
DECLARE @Input_Y_or_Lng FLOAT = 3653063.504053675
DECLARE @Input_Geog_Zone INT = 14
DECLARE @Input_Geog_Hemi VARCHAR(1) = 'N'
------------------------------------------
--CONSTANTS
DECLARE @gbl_pi FLOAT = 3.14159265358979
DECLARE @gbl_sm_a FLOAT = 6378137.0;
DECLARE @gbl_sm_b FLOAT = 6356752.314;
DECLARE @gbl_sm_EccSquared FLOAT = 6.69437999013e-03;
DECLARE @gbl_UTMScaleFactor FLOAT = 0.9996;
--RadToDeg(rad): (rad / pi * 180.0)
--DegToRad(deg): (deg / 180.0 * pi)
------------------------------------------
IF @Operation = 'FromUTMToGeo'
BEGIN
----------------------------------
DECLARE @x float = 0
DECLARE @y float = 0
DECLARE @southhemi bit = 0
IF @Input_Geog_Hemi = 'S' SET @southhemi = 1
----------------------------------
--X and Y
SET @x = (@Input_X_or_Lat - 500000.0)
SET @x = @x / @gbl_UTMScaleFactor
IF @southhemi = 1 BEGIN SET @y = (@Input_Y_or_Lng - 10000000.0) END
SET @y = @Input_Y_or_Lng / @gbl_UTMScaleFactor
----------------------------------
--cmeridian
-- DegToRad (-183.0 + (zone * 6.0));
DECLARE @cmeridian FLOAT = (-183.0 + (@Input_Geog_Zone * 6.0)) / 180 * @gbl_pi
DECLARE @lambda0 FLOAT = @cmeridian
----------------------------------
--Calculate Latitude and Longtitude
DECLARE @phif FLOAT
DECLARE @Nf FLOAT
DECLARE @Nfpow FLOAT
DECLARE @nuf2 FLOAT
DECLARE @ep2 FLOAT
DECLARE @tf FLOAT
DECLARE @tf2 FLOAT
DECLARE @tf4 FLOAT
DECLARE @cf FLOAT
DECLARE @x1frac FLOAT
DECLARE @x2frac FLOAT
DECLARE @x3frac FLOAT
DECLARE @x4frac FLOAT
DECLARE @x5frac FLOAT
DECLARE @x6frac FLOAT
DECLARE @x7frac FLOAT
DECLARE @x8frac FLOAT
DECLARE @x2poly FLOAT
DECLARE @x3poly FLOAT
DECLARE @x4poly FLOAT
DECLARE @x5poly FLOAT
DECLARE @x6poly FLOAT
DECLARE @x7poly FLOAT
DECLARE @x8poly FLOAT
DECLARE @Lat_AsRad FLOAT
DECLARE @Lng_AsRad FLOAT
DECLARE @Latitude FLOAT
DECLARE @Longitude FLOAT
-------------------------------------------------------------
/* Get the value of phif, the footpoint latitude. */
DECLARE @y_ FLOAT
DECLARE @alpha_ FLOAT
DECLARE @beta_ FLOAT
DECLARE @gamma_ FLOAT
DECLARE @delta_ FLOAT
DECLARE @epsilon_ FLOAT
DECLARE @n FLOAT
DECLARE @FootpointLatitude FLOAT
/* Precalculate n (Eq. 10.18) */
SET @n = (@gbl_sm_a - @gbl_sm_b) / (@gbl_sm_a + @gbl_sm_b);
/* Precalculate alpha_ (Eq. 10.22) */
/* (Same as alpha in Eq. 10.17) */
SET @alpha_ = ((@gbl_sm_a + @gbl_sm_b) / 2.0) * (1 + (POWER(@n, 2.0) / 4) + (POWER(@n, 4.0) / 64));
/* Precalculate y_ (Eq. 10.23) */
SET @y_ = @y / @alpha_;
/* Precalculate beta_ (Eq. 10.22) */
SET @beta_ = (3.0 * @n / 2.0) + (-27.0 * POWER(@n, 3.0) / 32.0) + (269.0 * POWER(@n, 5.0) / 512.0);
/* Precalculate gamma_ (Eq. 10.22) */
SET @gamma_ = (21.0 * POWER(@n, 2.0) / 16.0) + (-55.0 * POWER(@n, 4.0) / 32.0);
/* Precalculate delta_ (Eq. 10.22) */
SET @delta_ = (151.0 * POWER(@n, 3.0) / 96.0) + (-417.0 * POWER(@n, 5.0) / 128.0);
/* Precalculate epsilon_ (Eq. 10.22) */
SET @epsilon_ = (1097.0 * POWER(@n, 4.0) / 512.0);
/* Now calculate the sum of the series (Eq. 10.21) */
SET @FootpointLatitude = @y_ + (@beta_ * SIN(2.0 * @y_))
+ (@gamma_ * SIN(4.0 * @y_))
+ (@delta_ * SIN(6.0 * @y_))
+ (@epsilon_ * SIN(8.0 * @y_));
SET @phif = @FootpointLatitude
-------------------------------------------------------------
/* Precalculate ep2 */
SET @ep2 = (POWER(@gbl_sm_a, 2.0) -POWER(@gbl_sm_b, 2.0)) / POWER(@gbl_sm_b, 2.0);
/* Precalculate cos (phif) */
SET @cf = COS(@phif);
/* Precalculate nuf2 */
SET @nuf2 = @ep2 * POWER(@cf, 2.0);
/* Precalculate Nf and initialize Nfpow */
SET @Nf = POWER(@gbl_sm_a, 2.0) / (@gbl_sm_b * SQRT(1 + @nuf2));
SET @Nfpow = @Nf;
/* Precalculate tf */
SET @tf = TAN(@phif);
SET @tf2 = @tf * @tf;
SET @tf4 = @tf2 * @tf2;
/* Precalculate fractional coefficients for x**n in the equations
below to simplify the expressions for latitude and longitude. */
SET @x1frac = 1.0 / (@Nfpow * @cf);
SET @Nfpow *= @Nf; /* now equals Nf**2) */
SET @x2frac = @tf / (2.0 * @Nfpow);
SET @Nfpow *= @Nf; /* now equals Nf**3) */
SET @x3frac = 1.0 / (6.0 * @Nfpow * @cf);
SET @Nfpow *= @Nf; /* now equals Nf**4) */
SET @x4frac = @tf / (24.0 * @Nfpow);
SET @Nfpow *= @Nf; /* now equals Nf**5) */
SET @x5frac = 1.0 / (120.0 * @Nfpow * @cf);
SET @Nfpow *= @Nf; /* now equals Nf**6) */
SET @x6frac = @tf / (720.0 * @Nfpow);
SET @Nfpow *= @Nf; /* now equals Nf**7) */
SET @x7frac = 1.0 / (5040.0 * @Nfpow * @cf);
SET @Nfpow *= @Nf; /* now equals Nf**8) */
SET @x8frac = @tf / (40320.0 * @Nfpow);
/* Precalculate polynomial coefficients for x**n.
-- x**1 does not have a polynomial coefficient. */
SET @x2poly = -1.0 - @nuf2;
SET @x3poly = -1.0 - 2 * @tf2 - @nuf2;
SET @x4poly = 5.0 + 3.0 * @tf2 + 6.0 * @nuf2 - 6.0 * @tf2 * @nuf2 - 3.0 * (@nuf2 *@nuf2) - 9.0 * @tf2 * (@nuf2 * @nuf2);
SET @x5poly = 5.0 + 28.0 * @tf2 + 24.0 * @tf4 + 6.0 * @nuf2 + 8.0 * @tf2 * @nuf2;
SET @x6poly = -61.0 - 90.0 * @tf2 - 45.0 * @tf4 - 107.0 * @nuf2 + 162.0 * @tf2 * @nuf2;
SET @x7poly = -61.0 - 662.0 * @tf2 - 1320.0 * @tf4 - 720.0 * (@tf4 * @tf2);
SET @x8poly = 1385.0 + 3633.0 * @tf2 + 4095.0 * @tf4 + 1575 * (@tf4 * @tf2);
/* Calculate latitude */
SET @Lat_AsRad = @phif + @x2frac * @x2poly * (@x * @x)
+ @x4frac * @x4poly * POWER(@x, 4.0)
+ @x6frac * @x6poly * POWER(@x, 6.0)
+ @x8frac * @x8poly * POWER(@x, 8.0);
--Radians to Degrees
SET @Latitude = (@Lat_AsRad / @gbl_pi * 180.0)
/* Calculate longitude */
SET @Lng_AsRad = @lambda0 + @x1frac * @x
+ @x3frac * @x3poly * POWER(@x, 3.0)
+ @x5frac * @x5poly * POWER(@x, 5.0)
+ @x7frac * @x7poly * POWER(@x, 7.0);
--Radians to Degrees
SET @Longitude = (@Lng_AsRad / @gbl_pi * 180.0)
------------------------------------------------------------------
SELECT @Latitude As Lat, @Longitude As Lng
END
Thanks for any assistance!
