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

Posted 2023-02-14

【中文标题】使用从笛卡尔空间和世界文件生成的纬度和经度计算多边形面积

【英文标题】：Polygon area calculation using Latitude and Longitude generated from Cartesian space and a world file

【发布时间】：2011-02-21 02:54:55

【问题描述】：

给定一系列 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%）。有什么可以改进的吗？

谢谢。

【参考方案1】：

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

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

几个要点

此函数的数组输入将有“n + 1”个元素。最后一个元素的值与第一个元素的值相同。 我已经编写了非常基本的 C# 代码，因此人们也可以将其改编成其他语言。 6378137 是以米为单位的地球半径值。

输出面积单位为平方米

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;

【讨论】：

我试过你的代码，但有问题。有任何想法吗？见：code

你已经输入了“面积 = 面积 * R * R / 2;”在 for 循环内，而它应该在循环外。

我认为您也应该将p1.Longitude 和p2.Longitude 转换为弧度。进行此修改后，我得到了与 google.maps.geometry.spherical.computeArea 函数非常相似的结果。

修正后这似乎很好。并且与开放层中的getGeodesicArea 非常相似（减去投影部分）。见：github.com/openlayers/openlayers/blob/v2.13.1/lib/OpenLayers/…

【参考方案2】：

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

【参考方案3】：

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

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

【讨论】：

我正在与已知区域的实际测量值进行比较。一个有趣的旁注是，如果我通过 Haversine 方法运行 GPS 坐标，我会得到非常准确的距离计算，从而产生准确的周长值。

抱歉回复晚了，最终使用了带有 proj4 java 库的 WGs84 大地水准面。效果很好，感谢您的帮助。

【参考方案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
;

【参考方案5】：

将 RiskyPathak 的 sn-p 改编为 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;

【参考方案6】：

谢谢你有风险的帕塔克！

本着分享的精神，这是我在 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;

【参考方案7】：

将 RiskyPathak 的 sn-p 改编为 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