0
我試圖自動計算OpenStreetMaps上顯示的建築物的面積。
爲了做到這一點，我通過立交橋API獲取了建築物多邊形的座標。與PostGis ST_Area相比，爲什麼我會得到一個稍微不同的區域？
[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方面myselfs，我把它們放到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平方米，橢圓形爲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會 - 得到更高的結果。
我真的好奇我爲什麼會得到這樣不同的結果。
使用谷歌，我發現這裏 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()
（但只有當polygon用左手法則定義[這就是爲什麼陣列在這裏顛倒]，否則我得到System.ArgumentException: 24200: The specified input does not represent a valid geography instance.）
參考：
ST_Area
Source Code
