给定一系列GPS坐标对，我需要计算多边形（n-gon）的面积。这个面积相对较小（不大于50,000平方英尺）。地理编码是通过使用来自世界文件的数据进行仿射变换创建的。
我已经尝试使用两步方法，将地理编码转换为笛卡尔坐标：
double xPos = (lon-lonAnchor)*( Math.toRadians( 6378137 ) )*Math.cos( latAnchor );
double yPos = (lat-latAnchor)*( Math.toRadians( 6378137 ) );
然后我使用叉积计算来确定面积。
问题是结果的准确度有些偏差（大约1%）。有什么方法可以改进吗？
谢谢。
我在互联网上查找了各种多边形面积公式(或代码)，但没有找到任何一个好的或容易实现的。
现在我已经编写了代码片段，可以计算绘制在地球表面上的多边形的面积。该多边形可以有n个顶点，每个顶点都有自己的纬度和经度。
几个重要的要点
输出面积的单位为平方米
private static double CalculatePolygonArea(IList<MapPoint> coordinates)
{
double area = 0;
if (coordinates.Count > 2)
{
for (var i = 0; i < coordinates.Count - 1; i++)
{
MapPoint p1 = coordinates[i];
MapPoint p2 = coordinates[i + 1];
area += ConvertToRadian(p2.Longitude - p1.Longitude) * (2 + Math.Sin(ConvertToRadian(p1.Latitude)) + Math.Sin(ConvertToRadian(p2.Latitude)));
}
area = area * 6378137 * 6378137 / 2;
}
return Math.Abs(area);
}
private static double ConvertToRadian(double input)
{
return input * Math.PI / 180;
}
double xPos = (lon-lonAnchor)*( Math.toRadians( 6378137 ) )*Math.cos( latAnchor );
成为：
double xPos = (lon-lonAnchor)*( 6378137*PI/180 ) )*Math.cos( latAnchor*PI/180 );
其中lon、lonAnchor和latAnchor的单位为度。现在完美运作。
我不确定为什么，但是使用上述公式得到的结果与Google地图测量同一区域时返回的结果相去甚远。
因此，我进行了搜索，并提出了这种由GPS坐标定义的多边形面积计算的JavaScript方法：
const PI = Math.PI;
const EARTH_RADIUS_IN_METERS = 6378137;
const EARTH_CIRCUMFERENCE_IN_METERS = 2*EARTH_RADIUS_IN_METERS*PI;
function areaClaculator(points) {
let area = null;
if (!isValueEmpty(points) && points.length > 2) {
let p0 = points[0]
let newPoints = [];
for (let i=1; i<points.length; i++) {
let p = points[i];
let y = (p.lat - p0.lat) / 360 * EARTH_CIRCUMFERENCE_IN_METERS;
let x = (p.lng - p0.lng) / 360 * EARTH_CIRCUMFERENCE_IN_METERS * Math.cos(rad(p.lat));
let entry = {};
entry.x = x;
entry.y = y;
newPoints.push(entry);
}
if (!isValueEmpty(newPoints) && newPoints.length > 1) {
area = 0;
for (let i=0;i< newPoints.length - 1; i++) {
let p1 = newPoints[i];
let p2 = newPoints[i+1];
area += ((p1.y * p2.x) - (p1.x*p2.y))/2;
}
area = Math.abs(area);
}
}
return area;
}
function rad(degrees) {
return degrees * PI / 180;
}
points 存储的值如下：{lng: -73.462556156587, lat: 45.48566183708046}
尝试在Swift Playground中执行此操作，但结果偏差很大。我正在将示例坐标（39.58571008386715，-104.94522892318253）输入函数中。
func deg2rad(_ number: Double) -> Double {
return number * .pi / 180
}
func areaCalc(lat: [Double]?, lon: [Double]?){
guard let lat = lat,
let lon = lon
else { return }
var area: Double = 0.0
if(lat.count > 2){
for i in stride(from: 0, to: lat.count - 1, by: 1) {
let p1lon = lon[i]
let p1lat = lat[i]
let p2lon = lon[i+1]
let p2lat = lat[i+1]
area = area + (deg2rad(p2lon - p1lon)) * (2 + sin(deg2rad(p1lat))) + (sin(deg2rad(p2lat)))
}
area = area * 6378137.0 * 6378137.0
area = abs(area / 2)
}
}
将RiskyPathak的片段改编为PHP
function CalculatePolygonArea($coordinates) {
$area = 0;
$coordinatesCount = sizeof($coordinates);
if ($coordinatesCount > 2) {
for ($i = 0; $i < $coordinatesCount - 1; $i++) {
$p1 = $coordinates[$i];
$p2 = $coordinates[$i + 1];
$p1Longitude = $p1[0];
$p2Longitude = $p2[0];
$p1Latitude = $p1[1];
$p2Latitude = $p2[1];
$area += ConvertToRadian($p2Longitude - $p1Longitude) * (2 + sin(ConvertToRadian($p1Latitude)) + sin(ConvertToRadian($p2Latitude)));
}
$area = $area * 6378137 * 6378137 / 2;
}
return abs(round(($area));
}
function ConvertToRadian($input) {
$output = $input * pi() / 180;
return $output;
}
基于Risky Pathak的解决方案，这里提供了SQL（Redshift）计算GeoJSON多边形面积的解决方案（假设linestring 0是最外层多边形）。
create or replace view geo_area_area as
with points as (
select ga.id as key_geo_area
, ga.name, gag.linestring
, gag.position
, radians(gag.longitude) as x
, radians(gag.latitude) as y
from geo_area ga
join geo_area_geometry gag on (gag.key_geo_area = ga.id)
)
, polygons as (
select key_geo_area, name, linestring, position
, x
, lag(x) over (partition by key_geo_area, linestring order by position) as prev_x
, y
, lag(y) over (partition by key_geo_area, linestring order by position) as prev_y
from points
)
, area_linestrings as (
select key_geo_area, name, linestring
, abs( sum( (x - prev_x) * (2 + sin(y) + sin(prev_y)) ) ) * 6378137 * 6378137 / 2 / 10^6 as area_km_squared
from polygons
where position != 0
group by 1, 2, 3
)
select key_geo_area, name
, sum(case when linestring = 0 then area_km_squared else -area_km_squared end) as area_km_squared
from area_linestrings
group by 1, 2
;
谢谢Risky Pathak！
为了分享精神，这是我在Delphi中的改编：
interface
uses
System.Math;
TMapGeoPoint = record
Latitude: Double;
Longitude: Double;
end;
function AreaInAcres(AGeoPoints: TList<TMapGeoPoint>): Double;
implementation
function AreaInAcres(AGeoPoints: TList<TMapGeoPoint>): Double;
var
Area: Double;
i: Integer;
P1, P2: TMapGeoPoint;
begin
Area := 0;
// We need at least 2 points
if (AGeoPoints.Count > 2) then
begin
for I := 0 to AGeoPoints.Count - 1 do
begin
P1 := AGeoPoints[i];
if i < AGeoPoints.Count - 1 then
P2 := AGeoPoints[i + 1]
else
P2 := AGeoPoints[0];
Area := Area + DegToRad(P2.Longitude - P1.Longitude) * (2 +
Sin(DegToRad(P1.Latitude)) + Sin(DegToRad(P2.Latitude)));
end;
Area := Area * 6378137 * 6378137 / 2;
end;
Area := Abs(Area); //Area (in sq meters)
// 1 Square Meter = 0.000247105 Acres
result := Area * 0.000247105;
end;
将 RiskyPathak 的片段改编为 Ruby
def deg2rad(input)
input * Math::PI / 180.0
end
def polygone_area(coordinates)
return 0.0 unless coordinates.size > 2
area = 0.0
coor_p = coordinates.first
coordinates[1..-1].each{ |coor|
area += deg2rad(coor[1] - coor_p[1]) * (2 + Math.sin(deg2rad(coor_p[0])) + Math.sin(deg2rad(coor[0])))
coor_p = coor
}
(area * 6378137 * 6378137 / 2.0).abs # 6378137 Earth's radius in meters
end
p1.Longitude和p2.Longitude也转换为弧度。在进行这个修改后，我得到了与google.maps.geometry.spherical.computeArea函数相似的结果。 - MC XgetGeodesicArea非常相似（减去投影部分）。请参见：https://github.com/openlayers/openlayers/blob/v2.13.1/lib/OpenLayers/Geometry/LinearRing.js#L251 - Nux