与 PostGis ST_Area 相比,为什么我得到的区域略有不同?

Why exactly do I get a slightly different area compared to PostGis ST_Area?

我正在尝试自动计算 OpenStreetMaps 上显示的建筑物的面积。
为此,我通过立交桥 API 获取建筑物多边形的坐标。

例如 https://overpass-turbo.eu/

[out:json];
way(29858799);
out ids geom;

哪个输出我

{
  "version": 0.6,
  "generator": "Overpass API",
  "osm3s": {
    "timestamp_osm_base": "2017-10-09T09:54:02Z",
    "copyright": "The data included in this document is from www.openstreetmap.org. The data is made available under ODbL."
  },
  "elements": [

{
  "type": "way",
  "id": 29858799,
  "bounds": {
    "minlat": 47.3604067,
    "minlon": 8.5342631,
    "maxlat": 47.3612503,
    "maxlon": 8.5352457
  },
  "geometry": [
    { "lat": 47.3612503, "lon": 8.5351944 },
    { "lat": 47.3612252, "lon": 8.5342631 },
    { "lat": 47.3610145, "lon": 8.5342755 },
    { "lat": 47.3610212, "lon": 8.5345227 },
    { "lat": 47.3606405, "lon": 8.5345451 },
    { "lat": 47.3606350, "lon": 8.5343411 },
    { "lat": 47.3604067, "lon": 8.5343545 },
    { "lat": 47.3604120, "lon": 8.5345623 },
    { "lat": 47.3604308, "lon": 8.5352457 },
    { "lat": 47.3606508, "lon": 8.5352328 },
    { "lat": 47.3606413, "lon": 8.5348784 },
    { "lat": 47.3610383, "lon": 8.5348551 },
    { "lat": 47.3610477, "lon": 8.5352063 },
    { "lat": 47.3612503, "lon": 8.5351944 }
  ]
}

  ]
}

现在,在我自己计算 JavaScript 中的面积之前,我将它们放入 PostGIS 中,然后查看 PostGIS 给我的面积:

SELECT     
    ST_Area
       (
       ST_GeomFromText('
POLYGON
(
       (
             47.3612503 8.5351944
             ,47.3612252 8.5342631,47.3610145 8.5342755,47.3610212 8.5345227,47.3606405 8.5345451
             ,47.3606350 8.5343411,47.3604067 8.5343545,47.3604120 8.5345623,47.3604308 8.5352457
             ,47.3606508 8.5352328,47.3606413 8.5348784,47.3610383 8.5348551,47.3610477 8.5352063
             ,47.3612503 8.5351944
       )
)'
             ,4326 -- WGS84
             )
             ,false -- 
    )
    -- false: 6379.25032051953
    -- true:  6350.65051177517

其中椭圆形为 6379.25032051953 m2,椭球体为 6350.65051177517。

现在我尝试计算 JavaScript
中的面积 所以我把这些坐标放到一个 JS 数组中:

var poly = [
        [47.3612503, 8.5351944],
        [47.3612252, 8.5342631],
        [47.3610145, 8.5342755],
        [47.3610212, 8.5345227],
        [47.3606405, 8.5345451],
        [47.3606350, 8.5343411],
        [47.3604067, 8.5343545],
        [47.3604120, 8.5345623],
        [47.3604308, 8.5352457],
        [47.3606508, 8.5352328],
        [47.3606413, 8.5348784],
        [47.3610383, 8.5348551],
        [47.3610477, 8.5352063],
        [47.3612503, 8.5351944]
    ];

并尝试计算 JavaScript 中的面积,使用 this gis post 作为参考。

Math.radians = function(degrees)
{
    return degrees * Math.PI / 180.0;
};


// https://gis.stackexchange.com/a/816/3997
function polygonArea()
{
    var poly = [
        [47.3612503, 8.5351944],
        [47.3612252, 8.5342631],
        [47.3610145, 8.5342755],
        [47.3610212, 8.5345227],
        [47.3606405, 8.5345451],
        [47.3606350, 8.5343411],
        [47.3604067, 8.5343545],
        [47.3604120, 8.5345623],
        [47.3604308, 8.5352457],
        [47.3606508, 8.5352328],
        [47.3606413, 8.5348784],
        [47.3610383, 8.5348551],
        [47.3610477, 8.5352063],
        [47.3612503, 8.5351944]
    ];


    var area = 0.0;
    var len = poly.length;

    if (len > 2)
    {

        var p1, p2;

        for (var i = 0; i < len - 1; i++)
        {

            p1 = poly[i];
            p2 = poly[i + 1];

            area += Math.radians(p2[0] - p1[0]) *
                (
                    2
                    + Math.sin(Math.radians(p1[1]))
                    + Math.sin(Math.radians(p2[1]))
                );
        }

        // https://en.wikipedia.org/wiki/Earth_radius#Equatorial_radius
        // https://en.wikipedia.org/wiki/Earth_ellipsoid
        // The radius you are using, 6378137.0 m corresponds to the equatorial radius of the Earth.
        var equatorial_radius = 6378137; // m
        var polar_radius = 6356752.3142; // m
        var mean_radius = 6371008.8; // m
        var authalic_radius = 6371007.2; // m (radius of perfect sphere with same surface as reference ellipsoid)
        var volumetric_radius = 6371000.8 // m (radius of a sphere of volume equal to the ellipsoid)
        // geodetic latitude φ 
        var siteLatitude = Math.radians(poly[0][0]);


        // https://en.wikipedia.org/wiki/Semi-major_and_semi-minor_axes
        // https://en.wikipedia.org/wiki/World_Geodetic_System
        var a = 6378137; // m 
        var b = 6356752.3142; // m 
        // where a and b are, respectively, the equatorial radius and the polar radius.

        var R1 = Math.pow(a * a * Math.cos(siteLatitude), 2) + Math.pow(b * b * Math.sin(siteLatitude), 2)
        var R2 = Math.pow(a * Math.cos(siteLatitude), 2) + Math.pow(b * Math.sin(siteLatitude), 2);

        // https://en.wikipedia.org/wiki/Earth_radius#Radius_at_a_given_geodetic_latitude
        // Geocentric radius
        var R = Math.sqrt(R1 / R2);
        // var merid_radius = ((a * a) * (b * b)) / Math.pow(Math.pow(a * Math.cos(siteLatitude), 2) + Math.pow(b * Math.sin(siteLatitude), 2), 3/2)



        // console.log(R);
        // var hrad = polar_radius + (90 - Math.abs(siteLatitude)) / 90 * (equatorial_radius - polar_radius);
        var radius = mean_radius;

        area = area * radius * radius / 2.0;
    } // End if len > 0

    // equatorial_radius: 6391.565558418869 m2
    // mean_radius:       6377.287126172337m2
    // authalic_radius:   6377.283923019292 m2
    // volumetric_radius: 6377.271110415153 m2
    // merid_radius:      6375.314923754325 m2
    // polar_radius:      6348.777989748668 m2
    // R:                 6368.48180842528 m2
    // hrad:              6391.171919886588 m2

    // http://postgis.net/docs/doxygen/2.2/dc/d52/geography__measurement_8c_a1a7c48d59bcf4ed56522ab26c142f61d.html
    // ST_Area(false)     6379.25032051953
    // ST_Area(true)      6350.65051177517

    // return area;
    return area.toFixed(2);
}

但无论我选择哪个半径,我都至少与 PostGis 输出相距 2 平方米。

更重要的是,当我使用纬度 X 处的精确半径计算面积时,我更接近 PostGIS 球形结果,而当我选择椭球体半径时,我没有像 PostGis 那样得到较低的结果 -我得到了更高的结果。

我真的很想知道为什么我会得到如此不同的结果。

使用 google 我在这里找到 http://postgis.net/docs/doxygen/2.2/dc/d52/geography__measurement_8c_a1a7c48d59bcf4ed56522ab26c142f61d.html PostGIS ST_Area 调用 geography_area,现在我想知道这两个结果集中哪个更错误...

这个计算有什么问题吗? 还是归咎于 PostGIS?

具有讽刺意味的是,当我在 SQL-Server 中计算面积时,我得到了 PostGIS 球形面积(实际上是 6350.65051472187)...

DECLARE @v_polygon_string varchar(1000);
DECLARE @g Geography;
SET @v_polygon_string = 'POLYGON((
47.3612503 8.5351944, 
47.3610477 8.5352063, 
47.3610383 8.5348551, 
47.3606413 8.5348784, 
47.3606508 8.5352328, 
47.3604308 8.5352457, 
47.3604120 8.5345623, 
47.3604067 8.5343545, 
47.3606350 8.5343411, 
47.3606405 8.5345451, 
47.3610212 8.5345227, 
47.3610145 8.5342755, 
47.3612252 8.5342631, 
47.3612503 8.5351944
)) ';
SET @g = Geography::STGeomFromText(@v_polygon_string,4326);
SELECT @g.STArea()

(但前提是多边形是使用左手法则定义的[这就是数组在这里反转的原因],否则我得到 System.ArgumentException: 24200: The specified input does not represent a valid geography instance.

参考:
ST_Area
Source Code

啊,自己找到答案了。
这是因为 Google-Maps 使用 Web-Mercator 投影。

坐标是WGS84,为了从中计算面积,需要将坐标转换成保留区域的地图投影。

否则,您将得到一个实际面积未必保留的区域(取决于地球上的哪个位置),并且计算出的面积将与实际面积不同。

真正的根本问题是:What does ST_Area actually do
我通过查看地理空间图书馆找到了这个问题的答案。

实际上 6377.28 的结果以 0.02 m2 的精度与圆柱等积球体 (6377.2695087222382) 或 EckertIV-Sphere (6377.26664171461) 拟合。

虽然 6350 的结果符合(非球形)圆柱等积世界或阿尔伯斯投影。

Details on the stackexchange-GIS-site

阿尔伯斯投影:

圆柱等积投影: