使用从笛卡尔空间和世界文件生成的纬度和经度计算多边形区域

Question

给定一系列GPS坐标对，我需要计算多边形的面积（n-gon）。这相对较小（不大于50,000平方英尺）。通过对来自世界文件的数据应用仿射变换来创建地理编码。

我试图通过将地理编码转换为笛卡尔坐标来使用两步法：

double xPos = (lon-lonAnchor)*( Math.toRadians( 6378137 ) )*Math.cos( latAnchor );
double yPos = (lat-latAnchor)*( Math.toRadians( 6378137 ) );

然后我使用cross product计算来确定面积。

问题是结果的准确性有点偏差（约1％）。我有什么可以改进的吗？

谢谢。

Answer 1

由于您的近似值，1％的误差似乎有点高。您是在与实际测量值或理想计算值进行比较吗？请记住，GPS中也存在可能有所贡献的错误。

如果你想要一个更准确的方法来做这个，那么在this问题上有一个很好的答案。如果您想要更快的方式，可以使用WGS84大地水准面而不是参考球来转换为笛卡尔坐标（ECEF）。这是转换的wiki link。

Answer 2

我正在修改Google地图，以便用户可以通过单击顶点来计算多边形的面积。在确定Math.cos（latAnchor）首先是弧度之前，它没有提供正确的区域

所以：

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的度数。现在就像一个魅力。

Answer 3

我在互联网上检查了各种多边形区域公式（或代码），但没有找到任何好的或易于实现的。

现在我编写了代码片段来计算地球表面上绘制的多边形的面积。多边形可以有n个顶点，每个顶点都有自己的纬度经度。

几点重要

1. 输入此函数的数组将具有“n + 1”个元素。最后一个元素的值与第一个元素的值相同。
2. 我编写了非常基本的C＃代码，这样人们也可以用其他语言进行调整。
3. 6378137是地球半径的值，以米为单位。
4. 输出区域的单位为平方米

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

Answer 4

基于Risky Pathak的解决方案，这里是SQL（Redshift）计算GeoJSON multipolygons区域的解决方案（假设线串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
;

Answer 5

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

Answer 6

谢谢Risky Pathak！

本着共享的精神，这是我在德尔福的改编：

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;