将NAD83州平面坐标系转换为NAD83 UTM坐标系的数学公式(最终目标是转换为WGS84经纬度)

3

我的终极目标:

编写一个TSQL函数,将NAD83州平面坐标(特别是ESRI投影102738 http://spatialreference.org/ref/esri/102738/转换为WGS84经纬度,以便我可以在几乎实时的情况下,在Google地图上绘制原始坐标。

我不想使用外部库,如Proj4或GDAL,因为我想避免外部、非SQL依赖项,我担心未来管理员升级数据库服务器时不会重新安装这些依赖项。

  1. 我找到了JavaScript代码(http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html),用于将NAD83 UTM转换为WGS84经纬度(见下文)。我已成功地将其移植/转换为SQL。

  2. 我不知道如何找到数学公式,编写一个函数,可以将NAD83州平面坐标转换为NAD 83 UTM(然后我可以使用现有函数将UTM转换为经纬度)

我从JS移植的TSQL,用于将NAD83 UTM转换为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

感谢您的帮助!

你在这个存储过程上有进展了吗?我正在处理一个类似的存储过程,需要一些输入。 - crosan
不幸的是,为了使我们的数据更好地与Web地图配合使用,我们将SQL Server中的空间存储(无论它被称为什么)更改为SQL Server地理格式,并且现在使用WGS84坐标系。起初这并不受欢迎,但现在我们已经习惯了,现在将我们的数据放在Web地图上非常容易。 - user651745
你知道我可以用哪些库来实现这个吗? - Rotail
1个回答

2
我做了类似的事情,但不得不对公式进行修正,以创建最接近我的区域的逼近值。我知道有更好的处理修正的方法,但这对于我的预期目的已经足够接近了。
这个函数基于NOAA程序SPCS83(版本2.1) http://www.ngs.noaa.gov/PC_PROD/SPCS83/ 参考文献:http://www.softwright.com/faq/engineering/Coordinate%20Systems.html 以下是纬度函数:
DECLARE @X numeric(15,8)
DECLARE @Y numeric(15,8)
DECLARE @t numeric(15,8)
DECLARE @lon numeric(15,8)
DECLARE @x1 numeric(15,8)
DECLARE @lat0 numeric(15,8)
DECLARE @part1 numeric(15,8)
DECLARE @lat1 numeric(15,8)
DECLARE @lat numeric(15,8)

set @X =((@xCoord * .3048) - 500000)
set @Y =(@yCoord * .3048)
set @t = power(((Sqrt(power(@x,2) + Power((6184598.958 - @y),2))) / 11760085.63),1.376682)
set @lon = ((ATAN(@X / (6184598.958 - @Y))) / 0.726384231) + -2.103121749
set @x1 = @X + 500000
set @lat0 = 1.570796326 - (2 * Atan(power( ((Sqrt(power(@x,2) + power((6184598.958 - @y),2))) / 11760085.63) ,1.376682)))
set @part1 = (1 - (0.08227185422 * Sin(@lat0))) / (1 + (0.08227185422 * Sin(@lat0)))
set @lat1 = .570796326 - (2 * Atan(power(((Sqrt(Power(@x,2) + power((6184598.958 - @y),2))) / 11760085.63),1.376682) * power(@part1, 0.041136)))

While Abs(@lat1 - @lat0) > 0.000000003
    Begin
         Set @lat0 = @lat1
         Set @part1 = (1 - (0.08227185422 * Sin(@lat0))) / (1 + (0.08227185422 * Sin(@lat0)))
         Set @lat1 = 1.570796326 - (2 * Atan( @t * ( power(@part1 ,(0.08227185422 / 2) ) ) ))   
    END

set @lat = @lat1 / 0.01745329252
set @lon = @lon / 0.01745329252

set @lat=(((@lat * 100000) / 100000) + .00331)
--set @lon =  @lon * 100000 
-- exec StatePlane2Geodetic 1084101.380000, 114080.430000

--print  cast(@lat as varchar(20)) + ', '    + cast(@lon as varchar(20))
return @lat

这是 Lon 函数:
BEGIN
DECLARE @X numeric(15,8)
DECLARE @Y numeric(15,8)
DECLARE @t numeric(15,8)
DECLARE @lon numeric(15,8)
DECLARE @x1 numeric(15,8)
DECLARE @lat0 numeric(15,8)
DECLARE @part1 numeric(15,8)
DECLARE @lat1 numeric(15,8)
DECLARE @lat numeric(15,8)

set @X =((@xCoord * .3048) - 500000)
set @Y =(@yCoord * .3048)
set @t = power(((Sqrt(power(@x,2) + Power((6184598.958 - @y),2))) / 11760085.63),1.376682)
set @lon = ((ATAN(@X / (6184598.958 - @Y))) / 0.726384231) + -2.103121749
set @x1 = @X + 500000
set @lat0 = 1.570796326 - (2 * Atan(power( ((Sqrt(power(@x,2) + power((6184598.958 - @y),2))) / 11760085.63) ,1.376682)))
set @part1 = (1 - (0.08227185422 * Sin(@lat0))) / (1 + (0.08227185422 * Sin(@lat0)))
set @lat1 = .570796326 - (2 * Atan(power(((Sqrt(Power(@x,2) + power((6184598.958 - @y),2))) / 11760085.63),1.376682) * power(@part1, 0.041136)))

While Abs(@lat1 - @lat0) > 0.000000003
    Begin
         Set @lat0 = @lat1
         Set @part1 = (1 - (0.08227185422 * Sin(@lat0))) / (1 + (0.08227185422 * Sin(@lat0)))
         Set @lat1 = 1.570796326 - (2 * Atan( @t * ( power(@part1 ,(0.08227185422 / 2) ) ) ))   
    END

set @lat = @lat1 / 0.01745329252
set @lon = @lon / 0.01745329252

set @lat=(((@lat * 100000) / 100000) + .00331)
--set @lon =  @lon * 100000 


--print  cast(@lat as varchar(20)) + ', '    + cast(@lon as varchar(20))
return @lon

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接