Convert feet to decimal degrees Florida East SRID=2236

252
3
Jump to solution
04-12-2021 09:18 AM
ScottMinter
New Contributor III

I have a shape column in a point Feature Class with projection srid = 2236, Florida East. I need the lat/long of the shape in decimal degrees in a 3rd party application. Using sql server I can get shape.STX and shape.STY but they are represented in feet.

Do you know how to take the STX in feet and convert it to decimal degrees? I found some vbscript and converted it to a sql function but I'm off by 10^3 on many points compared to the attribute table where I created a column to store x,y and ran calculate geometry to populate.

The programming must start and end in T-Sql, but I can call an outside program if necessary. I wouldn't mind calling an arcgis rest utility if I knew of one that would do the conversion for me.

One other alternative is to be able to automate the "calculate geometry" we perform manually to update the x,y of the shape. If I could run this nightly I could use the attribute table instead of trying to perform the math in t-sql.

0 Kudos
1 Solution

Accepted Solutions
ScottMinter
New Contributor III

This is the version that works. I had some stale "Calculate Geometry" values in my geodatabase. After updating these this solution works to get me to centimeter accuracy vs Esri's calculation for about 26,000 points I tested.

 

USE [ArcSDE]

GO

/****** Object:  UserDefinedFunction [dbo].[Langd_Convert_Feet_to_DecimalDegrees_Latitude]    Script Date: 4/13/2021 8:38:54 AM ******/

SET ANSI_NULLS ON

GO

SET QUOTED_IDENTIFIER ON

GO

-- =============================================

-- Author: Scott Minter

-- Create date: 9Apr21

-- Description: The purpose of this function is to convert STSRid 2236 lat and long from feet to decimal degrees.

--              This code was found in vb on Esri Community Website and I converted it to tsql

-- =============================================

ALTER FUNCTION [dbo].[Langd_Convert_Feet_to_DecimalDegrees_Latitude] 

(

-- Add the parameters for the function here

@y Numeric(38,8), -- latitude

@x Numeric(38,8)  -- longitude

)

RETURNS Numeric(38,8)

AS

BEGIN

 

 

 -- Constants below from NAD 1983 State Plane Florida East FIPS 0901(Feet).

 -- e_sqr is for NAD83 or GRS80 (http:en.wikipedia.org/wiki/NAD83#North_American_Datum_of_1983)

 /* From Esri WKT https://epsg.io/2236

 PROJCS["NAD83_Florida_East_ftUS",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID

 ["GRS_1980",6378137,298.257222101]],

 PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],

 PARAMETER["latitude_of_origin",24.33333333333333],PARAMETER["central_meridian",-81],

 PARAMETER["scale_factor",0.999941177],PARAMETER["false_easting",656166.667],PARAMETER["false_northing",0],

 UNIT["Foot_US",0.30480060960121924]]

 */

 Declare @false_easting Numeric(38,8) = 656166.667 -- original formula was (656166.6666666665 + (2/3)); -- false easting

 Declare @false_northing Numeric(38,8) = 0 -- original formula was (0 + (1/3))  -- Should be 0. this is a 'fudge' factor.

 Declare @a Numeric(38,8) =  6378137.0  --Equatorial Radius (meters) For GRS80

 Declare @e_sqr Numeric(38,8) =  0.006694300229 --Eccentricity of ellipsoid (GRS80)

 Declare @lat_o Numeric(38,8) = 24.33333333333333  --Origin Latitude FL state plane

 Declare @long_o Numeric(38,8) = 81.0 -- original formula was (81.0 + ((1/3) * 0.01)) -- Central Meridan FL state plane, no negative sign for this formula

 Declare @k_o Numeric(38,8) = 0.999941177 --0.9999411764705882  -- Scale factor

 

 Declare @pointx_ft Numeric(38,8) = @x    -- number   ft

 Set @pointx_ft = @x - @false_easting -- false easting accounted for.

 Declare @pointy_ft Numeric(38,8) = @y  --number      ft

 Set @pointy_ft = @y - @false_northing -- false northing accounted for.

 Declare @e_sqr_2 Numeric(38,8) = power(@e_sqr, 2)    -- need this to calculate M

 Declare @e_sqr_3 Numeric(38,8) = power(@e_sqr, 3)    -- need this to calculate M

 Declare @radianlat_o Numeric(38,8) = @lat_o/(180/PI())      --degrees converted to radians     57.2957

 Declare @radianLong_o Numeric(38,8) = @long_o/(180/PI())    --degrees converted to radians

 

 Declare @pointx_m Numeric(38,8) = @pointx_ft/3.2808 --@pointx_ft/3.2808399

 Declare @pointy_m Numeric(38,8) = @pointy_ft/3.2808 --@pointy_ft/3.2808399

 

 Declare @Mpart1_o Numeric(38,8) = ((1-(@e_sqr/4)-((3*@e_sqr_2)/64)-((5*@e_sqr_3)/256))*@radianlat_o)

 Declare @Mpart2_o Numeric(38,8) = ((((3*@e_sqr)/8)+((3*@e_sqr_2)/32)+((45*@e_sqr_3)/1024))*sin(2*@radianlat_o))

 Declare @Mpart3_o Numeric(38,8) = ((((15*@e_sqr_2)/256) + ((45*@e_sqr_3)/1024))*4*sin(4*@radianlat_o) )

 Declare @Mpart4_o Numeric(38,8) =  ((35*@e_sqr_3)/3072)*4*sin(6*@radianlat_o)

 Declare @M_o Numeric(38,8) = @a*(@Mpart1_o - @Mpart2_o + @Mpart3_o - @Mpart4_o)

 

  Declare @e_prim_sqr Numeric(38,8) = @e_sqr / (1 - @e_sqr)

  Declare @M Numeric(38,8) = @M_o + (@pointy_m / @k_o)

  Declare @e1 Numeric(38,8) = (1-sqrt((1 - @e_sqr))) / (1+sqrt((1 - @e_sqr)))             --note: constant only for @e_sqr of 0.00676866

  Declare @u Numeric(38,8) = @M / (   @a*( 1-(@e_sqr/4)- (3*(power(@e_sqr,2)/64) -   (5*(power(@e_sqr,3)/256) )))  )  --radians

 

 Declare @lat1_part1 Numeric(38,8) =  ((3 * (@e1/2)) - (27 * (power(@e1,3)/32) )) * sin(2*@u)

 Declare @lat1_part2 Numeric(38,8) =  ((21 * (power(@e1,2)/16)) - (55 * (power(@e1,4)/32) )) * sin(4*@u)

 Declare @lat1_part3 Numeric(38,8) =  ((151 * (power(@e1,3)/96))  * sin(6*@u)  )

 Declare @lat1_part4 Numeric(38,8) =  ((1097 * (power(@e1,4)/512))  * sin(8*@u)  )

 Declare @lat1 Numeric(38,8) = @u + @lat1_part1 + @lat1_part2 + @lat1_part3   + @lat1_part4

 Set @lat1 = @lat1 * (180/PI())

 Declare @radian_lat1 Numeric(38,8) = @lat1/(180/PI())    --degrees converted to radians

 

 Declare @C1 Numeric(38,8) = @e_prim_sqr * (0.5*(1 + cos(2*@radian_lat1)))

 Declare @T1 Numeric(38,8) = (1-cos(2*@radian_lat1))/(1+cos(2*@radian_lat1))

 Declare @N1 Numeric(38,8) = @a/(sqrt(1-(@e_sqr*0.5*(1-cos(2*@radian_lat1)))))

 Declare @R1 Numeric(38,8) = (@a*(1-@e_sqr))/ sqrt(power((1-@e_sqr * (0.5*(1 - cos(2*@radian_lat1))) ),3) )

 Declare @D Numeric(38,8) =  @pointx_m/(@N1*@k_o)

 

 Declare @Lat_final_part1 Numeric(38,8) = @lat1

 Declare @Lat_final_part2 Numeric(38,8) = (@N1 * (tan(@radian_lat1) / @R1))

 Declare @Lat_final_part3 Numeric(38,8) = power(@D,2)/2 - ((( 5 + (3*@T1) + (10*@C1) - (4* power(@C1,2)) - (9 * @e_prim_sqr)) * power(@D,4))/24)

 Declare @Lat_final_part4 Numeric(38,8) =  ((61 + (90 * @T1) + (298 * @C1) + (45 * power(@T1,2)) - (252 * @e_prim_sqr) - (3 * power(@C1,2)))* power(@D,6))/720

 Declare @Lat_final Numeric(38,8) = @Lat_final_part1 - ((@Lat_final_part2 *(@Lat_final_part3 + @Lat_final_part4))     * (180/PI()))

 

 Declare @Long_final_part1 Numeric(38,8) = @long_o

 Declare @Long_final_part2 Numeric(38,8) = @D   --D

 Declare @Long_final_part3 Numeric(38,8) = (((1 + (2*@T1) + @C1))* (power(@D,3))) /6

 Declare @Long_final_part4 Numeric(38,8) = 5 - (2 * @C1) + (28*@T1) - (3 * power(@C1,2)) + (8 * @e_prim_sqr) + (24 * power(@T1,2))

 Declare @Long_final_part4a Numeric(38,8) = @Long_final_part4 * (power(@D,5) )

 Declare @Long_final_part4b Numeric(38,8) = @Long_final_part4a /120

 Declare @Long_final Numeric(38,8) = (-@Long_final_part1 + (((@Long_final_part2 - (@Long_final_part3 + @Long_final_part4b))) /cos(@radian_lat1)) * 180/PI())

 

 RETURN @Lat_final 

 

END

 

 

/* just copy this and make a function for longitude like I did and then call it when needed. I use this to integrate with an import into Itron FCS software so we can read gas meters. Hope this helps someone out there. */

View solution in original post

0 Kudos
3 Replies
ScottMinter
New Contributor III

I found some code and converted it to T-sql. It works, but is not giving me the same result as Esri's calculate Geometry.  Here is the T-sql function (Note I did not use the degree/minute/second part but left it commented out in case someone else would want to use it):

 

FUNCTION [dbo].[Langd_Convert_Feet_to_DecimalDegrees_Latitude] 

(

-- Add the parameters for the function here

@y Numeric(24,12), -- latitude

@x Numeric(24,12)  -- longitude

)

RETURNS Numeric(24,12)

AS

BEGIN

-- Declare the return variable here

DECLARE @Result Numeric(24,12)

 

 -- Constants below from NAD 1983 State Plane Florida East FIPS 0901(Feet).

 -- e_sqr is for NAD83 or GRS80 (http:en.wikipedia.org/wiki/NAD83#North_American_Datum_of_1983)

 Declare @false_easting numeric(24,12) = (656166.6666666665 + (2/3)); -- false easting

 Declare @false_northing numeric(24,12) = (0 + (1/3))  -- Should be 0. this is a 'fudge' factor.

 Declare @a numeric(24,12) = 6378137.0  --Equatorial Radius (meters) For GRS80

 Declare @e_sqr numeric(24,12) =  0.006694300229 --Eccentricity of ellipsoid (GRS80)

 Declare @lat_o numeric(24,12) = 24.33333333333333   --Origin Latitude FL state plane

 Declare @long_o numeric(24,12) = (81.0 + ((1/3) * 0.01))    -- Central Meridan FL state plane, no negative sign for this formula

 Declare @k_o numeric(24,12) = 0.9999411764705882  -- Scale factor

 

 Declare @pointx_ft numeric(24,12) = @x    -- number   ft

 Set @pointx_ft = @x - @false_easting -- false easting accounted for.

 Declare @pointy_ft numeric(24,12) = @y  --number      ft

 Set @pointy_ft = @y - @false_northing -- false northing accounted for.

 Declare @e_sqr_2 numeric(24,12) = power(@e_sqr, 2)    -- need this to calculate M

 Declare @e_sqr_3 numeric(24,12) = power(@e_sqr, 3)    -- need this to calculate M

 Declare @radianlat_o numeric(24,12) = @lat_o/(180/PI())      --degrees converted to radians     57.2957

 Declare @radianLong_o numeric(24,12) = @long_o/(180/PI())    --degrees converted to radians

 

 Declare @pointx_m numeric(24,12) = @pointx_ft/3.2808399

 Declare @pointy_m numeric(24,12) = @pointy_ft/3.2808399

 

 Declare @Mpart1_o numeric(24,12) = ((1-(@e_sqr/4)-((3*@e_sqr_2)/64)-((5*@e_sqr_3)/256))*@radianlat_o)

 Declare @Mpart2_o numeric(24,12) = ((((3*@e_sqr)/8)+((3*@e_sqr_2)/32)+((45*@e_sqr_3)/1024))*sin(2*@radianlat_o))

 Declare @Mpart3_o numeric(24,12) = ((((15*@e_sqr_2)/256) + ((45*@e_sqr_3)/1024))*4*sin(4*@radianlat_o) )

 Declare @Mpart4_o numeric(24,12) =  ((35*@e_sqr_3)/3072)*4*sin(6*@radianlat_o)

 Declare @M_o numeric(24,12) = @a*(@Mpart1_o - @Mpart2_o + @Mpart3_o - @Mpart4_o)

 

  Declare @e_prim_sqr numeric(24,12) = @e_sqr / (1 - @e_sqr)

  Declare @M numeric(24,12) = @M_o + (@pointy_m / @k_o)

  Declare @e1 numeric(24,12) = (1-sqrt((1 - @e_sqr))) / (1+sqrt((1 - @e_sqr)))             --note: constant only for @e_sqr of 0.00676866

  Declare @u numeric(24,12) = @M / (   @a*( 1-(@e_sqr/4)- (3*(power(@e_sqr,2)/64) -   (5*(power(@e_sqr,3)/256) )))  )  --radians

 

 Declare @lat1_part1 numeric(24,12) =  ((3 * (@e1/2)) - (27 * (power(@e1,3)/32) )) * sin(2*@u)

 Declare @lat1_part2 numeric(24,12) =  ((21 * (power(@e1,2)/16)) - (55 * (power(@e1,4)/32) )) * sin(4*@u)

 Declare @lat1_part3 numeric(24,12) =  ((151 * (power(@e1,3)/96))  * sin(6*@u)  )

 Declare @lat1_part4 numeric(24,12) =  ((1097 * (power(@e1,4)/512))  * sin(8*@u)  )

 Declare @lat1 numeric(24,12) = @u + @lat1_part1 + @lat1_part2 + @lat1_part3   + @lat1_part4

 Set @lat1 = @lat1 * (180/PI())

 Declare @radian_lat1 numeric(24,12) = @lat1/(180/PI())    --degrees converted to radians

 

 Declare @C1 numeric(24,12) = @e_prim_sqr * (0.5*(1 + cos(2*@radian_lat1)))

 Declare @T1 numeric(24,12) = (1-cos(2*@radian_lat1))/(1+cos(2*@radian_lat1))

 Declare @N1 numeric(24,12) = @a/(sqrt(1-(@e_sqr*0.5*(1-cos(2*@radian_lat1)))))

 Declare @R1 numeric(24,12) = (@a*(1-@e_sqr))/ sqrt(power((1-@e_sqr * (0.5*(1 - cos(2*@radian_lat1))) ),3) )

 Declare @D numeric(24,12) =  @pointx_m/(@N1*@k_o)

 

 Declare @Lat_final_part1 numeric(24,12) = @lat1

 Declare @Lat_final_part2 numeric(24,12) = (@N1 * (tan(@radian_lat1) / @R1))

 Declare @Lat_final_part3 numeric(24,12) = power(@D,2)/2 - ((( 5 + (3*@T1) + (10*@C1) - (4* power(@C1,2)) - (9 * @e_prim_sqr)) * power(@D,4))/24)

 Declare @Lat_final_part4 numeric(24,12) =  ((61 + (90 * @T1) + (298 * @C1) + (45 * power(@T1,2)) - (252 * @e_prim_sqr) - (3 * power(@C1,2)))* power(@D,6))/720

 Declare @Lat_final numeric(24,12) = @Lat_final_part1 - ((@Lat_final_part2 *(@Lat_final_part3 + @Lat_final_part4))     * (180/PI()))

 

 Declare @Long_final_part1 numeric(24,12) = @long_o

 Declare @Long_final_part2 numeric(24,12) = @D   --D

 Declare @Long_final_part3 numeric(24,12) = (((1 + (2*@T1) + @C1))* (power(@D,3))) /6

 Declare @Long_final_part4 numeric(24,12) = 5 - (2 * @C1) + (28*@T1) - (3 * power(@C1,2)) + (8 * @e_prim_sqr) + (24 * power(@T1,2))

 Declare @Long_final_part4a numeric(24,12) = @Long_final_part4 * (power(@D,5) )

 Declare @Long_final_part4b numeric(24,12) = @Long_final_part4a /120

 Declare @Long_final numeric(24,12) = (-@Long_final_part1 + (((@Long_final_part2 - (@Long_final_part3 + @Long_final_part4b))) /cos(@radian_lat1)) * 180/PI())

 

 /*

 -- conversion to ddmmss

 Degreelat = Math.floor(Lat_final)

 Degreelong = (Math.floor(Long_final)  + 1) --adjustment

 minlat = Math.floor(((Lat_final - Degreelat)*60))

 minlong = Math.floor((((Long_final - Degreelong)*60))+ 1)

 seclat = Math.abs(((Lat_final - Degreelat - (minlat/60)) * 3600)) -- - 1)

 seclong = Math.abs(((Long_final - Degreelong - (minlong/60)) * 3600 ))

 minlong = Math.abs(minlong);

 

 

 if (@seclat >= 1) 

  Set @seclat = (@seclat -1);

 } else {

  seclat = (seclat + 59);

  minlat = (minlat -1)

 }

 

 seclat = seclat.toFixed(2);

 seclong = seclong.toFixed(2);

 

 

 --Declare latitude = Lat_final.toString();

 --Declare longitude = Long_final.toString();

 

 latitude = latitude.substr(0,10);

 longitude = longitude.substr(0,10);

 

 if (minlat<10) {

  minlat = "0" + minlat;

 }

 if (seclat<10) {

  seclat = "0" + seclat;

 }

 if (minlong<10) {

  minlong = "0" + minlong;

 }

 if (seclong<10) {

  seclong = "0" + seclong;

 }

 

 if (seclat == "60.00") {

  seclat = "00.00";

  minlat = minlat +1;

 }

 if (seclong == "60.00") {

  seclong = "00.00";

  minlong = minlong +1;

 }

 

 northing = uY.toFixed(0);

 easting = uX.toFixed(0);

 

 dojo.byId("coordDMS").innerHTML = " " + Degreelat + "&#176;" + minlat + "\'" + seclat + "\"N " + (-Degreelong) + "&#176;" + minlong + "\'" + seclong + "\"W ";

 dojo.byId("coordStPl").innerHTML = " Northing: " + northing + " Easting: " + easting + " ";

}

--END CUSTOM FUNCTION FOR LAT-LONG DISPLAY

*/

RETURN @Lat_final 

 

END

 

 

A comparison of Esri's calculate Geometry and the result from this function (and a duplicate function for longitude):

Esri Calculate Geometry Result for Shape: 28.55022776,          -81.70544149
This function results:                                       28.550685362983, -81.705464998695

0 Kudos
ScottMinter
New Contributor III

The conversion is not accurate enough for my purposes (yet). We are in the natural gas industry and the conversion difference between esri and the calculated position is off by a few houses (a hundred feet or more).

Point 2 (checked using google maps)
esri 28.56191545,-81.75219301
38,6 28.560284, -81.752160

 

So I'm close using tsql Shape.STX and Shape.STY and calling the custom function. But not matching up with Esri's calculation.

0 Kudos
ScottMinter
New Contributor III

This is the version that works. I had some stale "Calculate Geometry" values in my geodatabase. After updating these this solution works to get me to centimeter accuracy vs Esri's calculation for about 26,000 points I tested.

 

USE [ArcSDE]

GO

/****** Object:  UserDefinedFunction [dbo].[Langd_Convert_Feet_to_DecimalDegrees_Latitude]    Script Date: 4/13/2021 8:38:54 AM ******/

SET ANSI_NULLS ON

GO

SET QUOTED_IDENTIFIER ON

GO

-- =============================================

-- Author: Scott Minter

-- Create date: 9Apr21

-- Description: The purpose of this function is to convert STSRid 2236 lat and long from feet to decimal degrees.

--              This code was found in vb on Esri Community Website and I converted it to tsql

-- =============================================

ALTER FUNCTION [dbo].[Langd_Convert_Feet_to_DecimalDegrees_Latitude] 

(

-- Add the parameters for the function here

@y Numeric(38,8), -- latitude

@x Numeric(38,8)  -- longitude

)

RETURNS Numeric(38,8)

AS

BEGIN

 

 

 -- Constants below from NAD 1983 State Plane Florida East FIPS 0901(Feet).

 -- e_sqr is for NAD83 or GRS80 (http:en.wikipedia.org/wiki/NAD83#North_American_Datum_of_1983)

 /* From Esri WKT https://epsg.io/2236

 PROJCS["NAD83_Florida_East_ftUS",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID

 ["GRS_1980",6378137,298.257222101]],

 PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],

 PARAMETER["latitude_of_origin",24.33333333333333],PARAMETER["central_meridian",-81],

 PARAMETER["scale_factor",0.999941177],PARAMETER["false_easting",656166.667],PARAMETER["false_northing",0],

 UNIT["Foot_US",0.30480060960121924]]

 */

 Declare @false_easting Numeric(38,8) = 656166.667 -- original formula was (656166.6666666665 + (2/3)); -- false easting

 Declare @false_northing Numeric(38,8) = 0 -- original formula was (0 + (1/3))  -- Should be 0. this is a 'fudge' factor.

 Declare @a Numeric(38,8) =  6378137.0  --Equatorial Radius (meters) For GRS80

 Declare @e_sqr Numeric(38,8) =  0.006694300229 --Eccentricity of ellipsoid (GRS80)

 Declare @lat_o Numeric(38,8) = 24.33333333333333  --Origin Latitude FL state plane

 Declare @long_o Numeric(38,8) = 81.0 -- original formula was (81.0 + ((1/3) * 0.01)) -- Central Meridan FL state plane, no negative sign for this formula

 Declare @k_o Numeric(38,8) = 0.999941177 --0.9999411764705882  -- Scale factor

 

 Declare @pointx_ft Numeric(38,8) = @x    -- number   ft

 Set @pointx_ft = @x - @false_easting -- false easting accounted for.

 Declare @pointy_ft Numeric(38,8) = @y  --number      ft

 Set @pointy_ft = @y - @false_northing -- false northing accounted for.

 Declare @e_sqr_2 Numeric(38,8) = power(@e_sqr, 2)    -- need this to calculate M

 Declare @e_sqr_3 Numeric(38,8) = power(@e_sqr, 3)    -- need this to calculate M

 Declare @radianlat_o Numeric(38,8) = @lat_o/(180/PI())      --degrees converted to radians     57.2957

 Declare @radianLong_o Numeric(38,8) = @long_o/(180/PI())    --degrees converted to radians

 

 Declare @pointx_m Numeric(38,8) = @pointx_ft/3.2808 --@pointx_ft/3.2808399

 Declare @pointy_m Numeric(38,8) = @pointy_ft/3.2808 --@pointy_ft/3.2808399

 

 Declare @Mpart1_o Numeric(38,8) = ((1-(@e_sqr/4)-((3*@e_sqr_2)/64)-((5*@e_sqr_3)/256))*@radianlat_o)

 Declare @Mpart2_o Numeric(38,8) = ((((3*@e_sqr)/8)+((3*@e_sqr_2)/32)+((45*@e_sqr_3)/1024))*sin(2*@radianlat_o))

 Declare @Mpart3_o Numeric(38,8) = ((((15*@e_sqr_2)/256) + ((45*@e_sqr_3)/1024))*4*sin(4*@radianlat_o) )

 Declare @Mpart4_o Numeric(38,8) =  ((35*@e_sqr_3)/3072)*4*sin(6*@radianlat_o)

 Declare @M_o Numeric(38,8) = @a*(@Mpart1_o - @Mpart2_o + @Mpart3_o - @Mpart4_o)

 

  Declare @e_prim_sqr Numeric(38,8) = @e_sqr / (1 - @e_sqr)

  Declare @M Numeric(38,8) = @M_o + (@pointy_m / @k_o)

  Declare @e1 Numeric(38,8) = (1-sqrt((1 - @e_sqr))) / (1+sqrt((1 - @e_sqr)))             --note: constant only for @e_sqr of 0.00676866

  Declare @u Numeric(38,8) = @M / (   @a*( 1-(@e_sqr/4)- (3*(power(@e_sqr,2)/64) -   (5*(power(@e_sqr,3)/256) )))  )  --radians

 

 Declare @lat1_part1 Numeric(38,8) =  ((3 * (@e1/2)) - (27 * (power(@e1,3)/32) )) * sin(2*@u)

 Declare @lat1_part2 Numeric(38,8) =  ((21 * (power(@e1,2)/16)) - (55 * (power(@e1,4)/32) )) * sin(4*@u)

 Declare @lat1_part3 Numeric(38,8) =  ((151 * (power(@e1,3)/96))  * sin(6*@u)  )

 Declare @lat1_part4 Numeric(38,8) =  ((1097 * (power(@e1,4)/512))  * sin(8*@u)  )

 Declare @lat1 Numeric(38,8) = @u + @lat1_part1 + @lat1_part2 + @lat1_part3   + @lat1_part4

 Set @lat1 = @lat1 * (180/PI())

 Declare @radian_lat1 Numeric(38,8) = @lat1/(180/PI())    --degrees converted to radians

 

 Declare @C1 Numeric(38,8) = @e_prim_sqr * (0.5*(1 + cos(2*@radian_lat1)))

 Declare @T1 Numeric(38,8) = (1-cos(2*@radian_lat1))/(1+cos(2*@radian_lat1))

 Declare @N1 Numeric(38,8) = @a/(sqrt(1-(@e_sqr*0.5*(1-cos(2*@radian_lat1)))))

 Declare @R1 Numeric(38,8) = (@a*(1-@e_sqr))/ sqrt(power((1-@e_sqr * (0.5*(1 - cos(2*@radian_lat1))) ),3) )

 Declare @D Numeric(38,8) =  @pointx_m/(@N1*@k_o)

 

 Declare @Lat_final_part1 Numeric(38,8) = @lat1

 Declare @Lat_final_part2 Numeric(38,8) = (@N1 * (tan(@radian_lat1) / @R1))

 Declare @Lat_final_part3 Numeric(38,8) = power(@D,2)/2 - ((( 5 + (3*@T1) + (10*@C1) - (4* power(@C1,2)) - (9 * @e_prim_sqr)) * power(@D,4))/24)

 Declare @Lat_final_part4 Numeric(38,8) =  ((61 + (90 * @T1) + (298 * @C1) + (45 * power(@T1,2)) - (252 * @e_prim_sqr) - (3 * power(@C1,2)))* power(@D,6))/720

 Declare @Lat_final Numeric(38,8) = @Lat_final_part1 - ((@Lat_final_part2 *(@Lat_final_part3 + @Lat_final_part4))     * (180/PI()))

 

 Declare @Long_final_part1 Numeric(38,8) = @long_o

 Declare @Long_final_part2 Numeric(38,8) = @D   --D

 Declare @Long_final_part3 Numeric(38,8) = (((1 + (2*@T1) + @C1))* (power(@D,3))) /6

 Declare @Long_final_part4 Numeric(38,8) = 5 - (2 * @C1) + (28*@T1) - (3 * power(@C1,2)) + (8 * @e_prim_sqr) + (24 * power(@T1,2))

 Declare @Long_final_part4a Numeric(38,8) = @Long_final_part4 * (power(@D,5) )

 Declare @Long_final_part4b Numeric(38,8) = @Long_final_part4a /120

 Declare @Long_final Numeric(38,8) = (-@Long_final_part1 + (((@Long_final_part2 - (@Long_final_part3 + @Long_final_part4b))) /cos(@radian_lat1)) * 180/PI())

 

 RETURN @Lat_final 

 

END

 

 

/* just copy this and make a function for longitude like I did and then call it when needed. I use this to integrate with an import into Itron FCS software so we can read gas meters. Hope this helps someone out there. */

View solution in original post

0 Kudos