用于将英国操作系统坐标从东向/北向转换为经度和纬度的 SQL 函数

Question

请有人发布一个 SQL 函数来将东距/北距转换为经度/纬度。我知道它非常复杂，但我还没有找到任何人用 T-SQL 记录它。

此 javascript 代码 可以工作，但我无法将其转换为SQL。

我有 16,000 个坐标，需要将它们全部转换为纬度/经度。

这是我到目前为止所拥有的，但它没有超过 while 循环。

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

Answer 1

我已经在这个问题上挣扎了一段时间了。
我在 OSGB36 中有很多北向/东向点，必须定期进行动态转换。
请注意，下面的 UDF 将 OSGB36（地形测量）投影中的北距/东距转换为 WGS84 投影中的纬度/经度，以便它们可以在 Google 地图中使用。

/****** Object:  UserDefinedFunction [dbo].[NEtoLL]    Script Date: 09/06/2012 17:06:39 ******/
SET ANSI_NULLS ON
GO

SET QUOTED_IDENTIFIER ON
GO

CREATE FUNCTION [dbo].[NEtoLL] (@East INT, @North INT, @LatOrLng VARCHAR(3)) RETURNS FLOAT AS
BEGIN

--Author: Sandy Motteram
--Date:   06 September 2012

--UDF adapted from javascript at http://www.bdcc.co.uk/LatLngToOSGB.js
--found on page http://mapki.com/wiki/Tools:Snippets

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

--Sample values below
--DECLARE @East INT, @North INT, @LatOrLng VARCHAR(3)
--SELECT @East = 529000, @North = 183650 --that combo should be the corner of Camden High St and Delancey St


    DECLARE @Pi              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  @Pi = 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
    /*    , @e2
          , @ex
          , @n1
          , @n2
          , @n3
          , @OriginNorthings*/

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

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

    DECLARE @ReturnMe FLOAT
    SET @ReturnMe = 0

    IF @LatOrLng = 'Lat'
        SET @ReturnMe = @newLat
    IF @LatOrLng = 'Lng'
        SET @ReturnMe = @newLon

    RETURN @ReturnMe
END
GO

Answer 2

我最终使用以下 javascript 函数来转换值。我知道这不是一个 SQL 解决方案，但它为我完成了这项工作。

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

Answer 3

我需要相同的功能，而 JavaScript 使得与数据库交互变得困难。我已将您的 JS 转换为 PHP，这在更新数据库时可能更有用 - 即：查询表、循环结果集、调用函数、更新表。

function OSGridToLatLong($E, $N) {

  $a = 6377563.396; 
  $b = 6356256.910;              // Airy 1830 major & minor semi-axes
  $F0 = 0.9996012717;                             // NatGrid scale factor on central meridian
  $lat0 = 49*M_PI/180; 
  $lon0 = -2*M_PI/180;  // NatGrid true origin
  $N0 = -100000;
  $E0 = 400000;                     // northing & easting of true origin, metres
  $e2 = 1 - ($b*$b)/($a*$a);                          // eccentricity squared
  $n = ($a-$b)/($a+$b);
  $n2 = $n*$n;
  $n3 = $n*$n*$n;

  $lat=$lat0;
  $M=0;
  do {
    $lat = ($N-$N0-$M)/($a*$F0) + $lat;

    $Ma = (1 + $n + (5/4)*$n2 + (5/4)*$n3) * ($lat-$lat0);
    $Mb = (3*$n + 3*$n*$n + (21/8)*$n3) * sin($lat-$lat0) * cos($lat+$lat0);
    $Mc = ((15/8)*$n2 + (15/8)*$n3) * sin(2*($lat-$lat0)) * cos(2*($lat+$lat0));
    $Md = (35/24)*$n3 * sin(3*($lat-$lat0)) * cos(3*($lat+$lat0));
    $M = $b * $F0 * ($Ma - $Mb + $Mc - $Md);                // meridional arc

  } while ($N-$N0-$M >= 0.00001);  // ie until < 0.01mm

  $cosLat = cos($lat);
  $sinLat = sin($lat);
  $nu = $a*$F0/sqrt(1-$e2*$sinLat*$sinLat);              // transverse radius of curvature
  $rho = $a*$F0*(1-$e2)/pow(1-$e2*$sinLat*$sinLat, 1.5);  // meridional radius of curvature
  $eta2 = $nu/$rho-1;

  $tanLat = tan($lat);
  $tan2lat = $tanLat*$tanLat;
  $tan4lat = $tan2lat*$tan2lat;
  $tan6lat = $tan4lat*$tan2lat;
  $secLat = 1/$cosLat;
  $nu3 = $nu*$nu*$nu;
  $nu5 = $nu3*$nu*$nu;
  $nu7 = $nu5*$nu*$nu;
  $VII = $tanLat/(2*$rho*$nu);
  $VIII = $tanLat/(24*$rho*$nu3)*(5+3*$tan2lat+$eta2-9*$tan2lat*$eta2);
  $IX = $tanLat/(720*$rho*$nu5)*(61+90*$tan2lat+45*$tan4lat);
  $X = $secLat/$nu;
  $XI = $secLat/(6*$nu3)*($nu/$rho+2*$tan2lat);
  $XII = $secLat/(120*$nu5)*(5+28*$tan2lat+24*$tan4lat);
  $XIIA = $secLat/(5040*$nu7)*(61+662*$tan2lat+1320*$tan4lat+720*$tan6lat);

  $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;
  $lon = $lon0 + $X*$dE - $XI*$dE3 + $XII*$dE5 - $XIIA*$dE7;


  return array(

    'longitude' => $lon * 180 / M_PI,
    'latitude'  => $lat * 180 / M_PI

  );

}

Answer 4

如果有人对非 SQL 解决方案感兴趣，我强烈建议使用这个 http://www.howtocreate。 co.uk/php/gridref.php PHP/JavaScript 类。

这里要提到的一件重要的事情是该库支持 Helmert 转换。

PHP

$grutoolbox = Grid_Ref_Utils::toolbox();
$source_coords = Array(54.607720,-6.411990);
//get the ellipsoids that will be used
$Airy_1830 = $grutoolbox->get_ellipsoid('Airy_1830');
$WGS84 = $grutoolbox->get_ellipsoid('WGS84');
$Airy_1830_mod = $grutoolbox->get_ellipsoid('Airy_1830_mod');
//get the transform parameters that will be used
$UK_to_GPS = $grutoolbox->get_transformation('OSGB36_to_WGS84');
$GPS_to_Ireland = $grutoolbox->get_transformation('WGS84_to_Ireland65');
//convert to GPS coordinates
$gps_coords = $grutoolbox->Helmert_transform($source_coords,$Airy_1830,$UK_to_GPS,$WGS84);
//convert to destination coordinates
print $grutoolbox->Helmert_transform($source_coords,$WGS84,$GPS_to_Ireland,$Airy_1830_mod,$grutoolbox->HTML);

JavaScript

var grutoolbox = gridRefUtilsToolbox();
var sourceCoords = [54.607720,-6.411990];
//get the ellipsoids that will be used
var Airy1830 = grutoolbox.getEllipsoid('Airy_1830');
var WGS84 = grutoolbox.getEllipsoid('WGS84');
var Airy1830Mod = grutoolbox.getEllipsoid('Airy_1830_mod');
//get the transform parameters that will be used
var UKToGPS = grutoolbox.getTransformation('OSGB36_to_WGS84');
var GPSToIreland = grutoolbox.getTransformation('WGS84_to_Ireland65');
//convert to GPS coordinates
var gpsCoords = grutoolbox.HelmertTransform(sourceCoords,Airy1830,UKToGPS,WGS84);
//convert to destination coordinates
element.innerHTML = grutoolbox.HelmertTransform(sourceCoords,WGS84,GPSToIreland,Airy1830Mod,grutoolbox.HTML);

Answer 5

我在 .NET 中开发了一个库，可以从 transact sql 调用 将 WGS84/UTM 坐标转换为纬度和经度

它的作用恰恰相反，但由于它使用 CooperativeSharp，您可以下载代码并轻松更改它以从纬度/经度转换为改为 wgs84。

您可以从github下载它：

https://github.com/jv-garcia/UTM2LATITUDE

  usage:

    SELECT dbo.UTM2LATITUDE(723399.51,4373328.5,'S',30) AS Latitude, dbo.UTM2LONGITUDE(723399.51,4373328.5,'S',30) AS Longitude

    result: 

        39,4805657453054    -0,402592727245112

        <param name="XUTM">pos UTM X</param>
        <param name="YUTM">pos UTM Y</param>
        <param name="LatBand">Latitude band grid zone designation letter (see http://www.dmap.co.uk/utmworld.htm) </param>
        <param name="LongBand">Longitude band grid zone designation number (see http://www.dmap.co.uk/utmworld.htm) </param>

Answer 6

基于 @Niall 的答案，此 PHP 函数还包括从 OSGB36 到 WSG84 的最终转换。

当前的前 3 个答案（SQL、JavaScript 和 PHP）都使用相同的算法（我认为）来生成 OSGB36 纬度和经度。这是 1836 年的标准，基于当时对地球曲率的理解。 WSG84是1984年的标准，其中坐标精确定位的位置最多可能相差约150m。

资料来源：http://www.movable-type.co。英国/scripts/latlong-os-gridref.html

function OSGridToLatLong($E, $N)
{
  $a = 6377563.396;
  $b = 6356256.910;     // Airy 1830 major & minor semi-axes
  $F0 = 0.9996012717;   // NatGrid scale factor on central meridian
  $lat0 = 49*M_PI/180;
  $lon0 = -2*M_PI/180;  // NatGrid true origin
  $N0 = -100000;
  $E0 = 400000;         // northing & easting of true origin, metres
  $e2 = 1 - ($b*$b)/($a*$a);  // eccentricity squared
  $n = ($a-$b)/($a+$b);
  $n2 = $n*$n;
  $n3 = $n*$n*$n;

  $lat=$lat0;
  $M=0;
  do {
    $lat = ($N-$N0-$M)/($a*$F0) + $lat;

    $Ma = (1 + $n + (5/4)*$n2 + (5/4)*$n3) * ($lat-$lat0);
    $Mb = (3*$n + 3*$n*$n + (21/8)*$n3) * sin($lat-$lat0) * cos($lat+$lat0);
    $Mc = ((15/8)*$n2 + (15/8)*$n3) * sin(2*($lat-$lat0)) * cos(2*($lat+$lat0));
    $Md = (35/24)*$n3 * sin(3*($lat-$lat0)) * cos(3*($lat+$lat0));
    $M = $b * $F0 * ($Ma - $Mb + $Mc - $Md);  // meridional arc

  } while ($N-$N0-$M >= 0.00001); // ie until < 0.01mm

  $cosLat = cos($lat);
  $sinLat = sin($lat);
  $nu = $a*$F0/sqrt(1-$e2*$sinLat*$sinLat);               // transverse radius of curvature
  $rho = $a*$F0*(1-$e2)/pow(1-$e2*$sinLat*$sinLat, 1.5);  // meridional radius of curvature
  $eta2 = $nu/$rho-1;

  $tanLat = tan($lat);
  $tan2lat = $tanLat*$tanLat;
  $tan4lat = $tan2lat*$tan2lat;
  $tan6lat = $tan4lat*$tan2lat;
  $secLat = 1/$cosLat;
  $nu3 = $nu*$nu*$nu;
  $nu5 = $nu3*$nu*$nu;
  $nu7 = $nu5*$nu*$nu;
  $VII = $tanLat/(2*$rho*$nu);
  $VIII = $tanLat/(24*$rho*$nu3)*(5+3*$tan2lat+$eta2-9*$tan2lat*$eta2);
  $IX = $tanLat/(720*$rho*$nu5)*(61+90*$tan2lat+45*$tan4lat);
  $X = $secLat/$nu;
  $XI = $secLat/(6*$nu3)*($nu/$rho+2*$tan2lat);
  $XII = $secLat/(120*$nu5)*(5+28*$tan2lat+24*$tan4lat);
  $XIIA = $secLat/(5040*$nu7)*(61+662*$tan2lat+1320*$tan4lat+720*$tan6lat);

  $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;
  $lon = $lon0 + $X*$dE - $XI*$dE3 + $XII*$dE5 - $XIIA*$dE7;

  $f = 1/299.3249646; // Airy 1830
  // Convert to cartesian
  $h = 0;
  $sinLat = \sin($lat);
  $cosLat = \cos($lat);
  $sinLon = \sin($lon);
  $cosLon = \cos($lon);
  $eSq = 2 * $f - $f * $f; // 1st eccentricity squared ≡ (a²-b²)/a²
  $v = $a / \sqrt(1 - $eSq * $sinLat * $sinLat); // radius of curvature in prime vertical
  $x = ($v + $h) * $cosLat * $cosLon;
  $y = ($v + $h) * $cosLat * $sinLon;
  $z = ($v * (1 - $eSq) + $h) * $sinLat;

  // Transform from OSGB36 to WSG84
  $t = [446.448, -125.157, 542.060, -20.4894, 0.1502, 0.2470, 0.8421];
  $tx = $t[0];                       // x-shift in metres
  $ty = $t[1];                       // y-shift in metres
  $tz = $t[2];                       // z-shift in metres
  $s  = $t[3] / 1e6 + 1;             // scale: normalise parts-per-million to (s+1)
  $rx = ($t[4] / 3600) * M_PI / 180; // x-rotation: normalise arcseconds to radians
  $ry = ($t[5] / 3600) * M_PI / 180; // y-rotation: normalise arcseconds to radians
  $rz = ($t[6] / 3600) * M_PI / 180; // z-rotation: normalise arcseconds to radians
  $x2 = $tx + $x * $s - $y * $rz + $z * $ry;
  $y2 = $ty + $x * $rz + $y * $s - $z * $rx;
  $z2 = $tz - $x * $ry + $y * $rx + $z * $s;

  // Convert back to WSG84 latitude and longitude
  $a = 6378137;
  $b = 6356752.314245;
  $f = 1/298.257223563;

  $e2 = 2 * $f - $f * $f;      // 1st eccentricity squared ≡ (a²−b²)/a²
  $epsilon2 = $e2 / (1 - $e2); // 2nd eccentricity squared ≡ (a²−b²)/b²
  $p = \sqrt($x2 * $x2 + $y2 * $y2); // distance from minor axis
  $R = \sqrt($p * $p + $z2 * $z2); // polar radius
  // parametric latitude (Bowring eqn.17, replacing tanβ = z·a / p·b)
  $tanBeta = ($b * $z2) / ($a * $p) * (1 + $epsilon2 * $b / $R);
  $sinBeta = $tanBeta / \sqrt(1 + $tanBeta * $tanBeta);
  $cosBeta = $sinBeta / $tanBeta;
  // geodetic latitude (Bowring eqn.18: tanφ = z+epsilon²⋅b⋅sin³β / p−e²⋅cos³β)
  $lat = \is_nan($cosBeta) ? 0 :
    \atan2(
      $z2 + $epsilon2 * $b * $sinBeta * $sinBeta * $sinBeta,
      $p - $e2 * $a * $cosBeta * $cosBeta * $cosBeta
    );
  // longitude
  $lon = \atan2($y2, $x2);
  // height above ellipsoid (Bowring eqn.7)
  $sinLat = \sin($lat);
  $cosLat = \cos($lat);
  $v = $a / \sqrt(1 - $e2 * $sinLat * $sinLat); // length of the normal terminated by the minor axis
  $h = $p * $cosLat + $z2 * $sinLat - ($a * $a / $v);

  return array(
    'longitude' => $lon * 180 / M_PI,
    'latitude'  => $lat * 180 / M_PI,
    'height' => $h,
  );
}

Answer 7

基于之前的答案（谢谢！）我必须在 Google Bigquery 中执行此操作。

我使用以下命令在 BQ 中运行它，并在临时函数末尾添加一些 SQL：

CREATE TEMP FUNCTION OSGridToLatLong(E FLOAT64, N FLOAT64)
RETURNS STRUCT<latitude FLOAT64, longitude FLOAT64, height FLOAT64>
LANGUAGE js AS """
  var a = 6377563.396;
  var b = 6356256.910;     
  var F0 = 0.9996012717;   
  var lat0 = 49 * Math.PI / 180;
  var lon0 = -2 * Math.PI / 180;  
  var N0 = -100000;
  var E0 = 400000;         
  var e2 = 1 - (b*b)/(a*a);  
  var n = (a-b)/(a+b);
  var n2 = n*n;
  var n3 = n*n*n;

  var lat = lat0;
  var 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);  

  } while (N - N0 - M >= 0.00001);

  var cosLat = Math.cos(lat);
  var sinLat = Math.sin(lat);
  var nu = a * F0 / Math.sqrt(1 - e2 * sinLat * sinLat);               
  var rho = a * F0 * (1 - e2) / Math.pow(1 - e2 * sinLat * sinLat, 1.5);  
  var eta2 = nu / rho - 1;

  var tanLat = Math.tan(lat);
  var tan2lat = tanLat * tanLat;
  var tan4lat = tan2lat * tan2lat;
  var tan6lat = tan4lat * tan2lat;
  var secLat = 1 / cosLat;
  var nu3 = nu * nu * nu;
  var nu5 = nu3 * nu * nu;
  var 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);
  var dE2 = dE * dE;
  var dE3 = dE2 * dE;
  var dE4 = dE2 * dE2;
  var dE5 = dE3 * dE2;
  var dE6 = dE4 * dE2;
  var dE7 = dE5 * dE2;
  lat = lat - VII * dE2 + VIII * dE4 - IX * dE6;
  var lon = lon0 + X * dE - XI * dE3 + XII * dE5 - XIIA * dE7;

  var f = 1 / 299.3249646; 
  var h = 0;
  sinLat = Math.sin(lat);
  cosLat = Math.cos(lat);
  var sinLon = Math.sin(lon);
  var cosLon = Math.cos(lon);
  var eSq = 2 * f - f * f;
  var v = a / Math.sqrt(1 - eSq * sinLat * sinLat);
  var x = (v + h) * cosLat * cosLon;
  var y = (v + h) * cosLat * sinLon;
  var z = (v * (1 - eSq) + h) * sinLat;

  var t = [446.448, -125.157, 542.060, -20.4894, 0.1502, 0.2470, 0.8421];
  var tx = t[0];                      
  var ty = t[1];                       
  var tz = t[2];                       
  var s  = t[3] / 1e6 + 1;             
  var rx = (t[4] / 3600) * Math.PI / 180; 
  var ry = (t[5] / 3600) * Math.PI / 180; 
  var rz = (t[6] / 3600) * Math.PI / 180; 
  var x2 = tx + x * s - y * rz + z * ry;
  var y2 = ty + x * rz + y * s - z * rx;
  var z2 = tz - x * ry + y * rx + z * s;

  a = 6378137;
  b = 6356752.314245;
  f = 1 / 298.257223563;

  e2 = 2 * f - f * f;
  var epsilon2 = e2 / (1 - e2);
  var p = Math.sqrt(x2 * x2 + y2 * y2);
  var R = Math.sqrt(p * p + z2 * z2);
  var tanBeta = (b * z2) / (a * p) * (1 + epsilon2 * b / R);
  var sinBeta = tanBeta / Math.sqrt(1 + tanBeta * tanBeta);
  var cosBeta = sinBeta / tanBeta;
  lat = isNaN(cosBeta) ? 0 :
    Math.atan2(
      z2 + epsilon2 * b * sinBeta * sinBeta * sinBeta,
      p - e2 * a * cosBeta * cosBeta * cosBeta
    );
  lon = Math.atan2(y2, x2);
  sinLat = Math.sin(lat);
  cosLat = Math.cos(lat);
  v = a / Math.sqrt(1 - e2 * sinLat * sinLat);
  h = p * cosLat + z2 * sinLat - (a * a / v);

  return {latitude: lat * 180 / Math.PI, longitude: lon * 180 / Math.PI, height: h};
""";

WITH input_data AS (
  SELECT 
    574807.122806563 AS E, --Add you data here, can also be a whole table to convert
    221464.559972554 AS N
)

SELECT
  OSGridToLatLong(E, N).latitude AS latitude,
  OSGridToLatLong(E, N).longitude AS longitude

FROM
  input_data