8

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

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

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

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

问题是结果的准确性有点偏离（大约 1%）。有什么可以改善的吗？

谢谢。

4

7 回答 7

8

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

现在我已经编写了代码片段来计算在地球表面绘制的多边形的面积。多边形可以有 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;
}
于 2015-10-02T17:33:30.913 回答
3

我正在修改谷歌地图，以便用户可以通过单击顶点来计算多边形的面积。在我确定 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 以度为单位。现在就像一个魅力。

于 2011-04-28T14:42:43.660 回答
2

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

如果您想要一种更准确的方法来执行此操作，那么这个问题有一个很好的答案。如果您想要更快的方法，您可以使用 WGS84 大地水准面而不是参考球体来转换为笛卡尔坐标 (ECEF)。这是该转换的wiki 链接。

于 2010-05-18T21:59:49.807 回答
0

基于 Risky Pathak 的解决方案，这里是 SQL (Redshift) 计算GeoJSON多面体面积的解决方案（假设线串 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
;
于 2019-02-27T11:28:33.863 回答
0

将 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;
}
于 2019-03-21T19:27:46.720 回答
0

谢谢风险帕塔克！

本着分享的精神，这是我在 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;
于 2019-06-03T15:20:35.677 回答
0

将 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
于 2020-03-20T10:28:19.157 回答