与 PostGis ST_Area 相比，为什么我得到的区域略有不同？
Why exactly do I get a slightly different area compared to PostGis ST_Area?
我正在尝试自动计算 OpenStreetMaps 上显示的建筑物的面积。
为此，我通过立交桥 API 获取建筑物多边形的坐标。
[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。
阿尔伯斯投影：
圆柱等积投影：
我正在尝试自动计算 OpenStreetMaps 上显示的建筑物的面积。
为此，我通过立交桥 API 获取建筑物多边形的坐标。
[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。
阿尔伯斯投影：
圆柱等积投影：