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

Question

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

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

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

然后我使用 叉积< /a> 计算以确定面积。

问题是结果的准确性有点偏差（大约 1%）。我可以研究什么来改善这一点吗？

谢谢。

Answer 1

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

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

几个要点

1. 此函数的数组输入将包含“n + 1”个元素。最后一个元素的值与第一个元素的值相同。
2. 我编写了非常基本的 C# 代码，以便大家也可以将其改编为其他语言。
3. 6378137 是以米为单位的地球半径值。

4. 输出面积单位为平方米

   private static double CalculatePolygonArea(IList 坐标)
   {
       双倍面积=0；
   
       if (坐标.Count > 2)
       {
           for (var i = 0; i < 坐标.Count - 1; i++)
           {
               地图点 p1 = 坐标[i];
               地图点 p2 = 坐标[i + 1];
               面积 += ConvertToRadian(p2.Longitude - p1.Longitude) * (2 + Math.Sin(ConvertToRadian(p1.Latitude)) + Math.Sin(ConvertToRadian(p2.Latitude)));
           }
   
           面积=面积*6378137*6378137/2；
       }
   
       返回 Math.Abs​​(面积);
   }
   
   私有静态双ConvertToRadian（双输入）
   {
       返回输入 * Math.PI / 180;
   }
   

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

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

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

Answer 4

造成这种“1%”差异的原因是地球是非常轻微的椭球体，因此通过使用球形模型进行计算，给出或获取位置时，误差通常高达 0.3%。

Answer 5

基于 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
;

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

Answer 7

将 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

Answer 8

尝试在快速游乐场中执行此操作，但结果却相差甚远
示例坐标：(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)
    }
}

Answer 9

我不知道为什么，但我使用上面的公式得到的结果与谷歌地图测量同一区域时返回的结果相差甚远。

因此，我想出了这个 javascript 方法来计算由 GPS 坐标定义的多边形区域：

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

其中点存储如下值：{lng: -73.462556156587, lat: 45.48566183708046}

Answer 10

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

// mssing 关闭 )

Answer 11

Google Maps Utils 库提供了一种计算面积的方法。

您必须添加以下依赖项

implementation("com.google.maps.android:android-maps-utils:3.8.2")

，然后调用下面的方法并为其提供您的纬度和经度列表。它将返回以平方米为单位的面积。

fun surfaceArea(list: List<LatLng>): Double {
    if (list.size < 3) {
        return 0.0
    }
    return SphericalUtil.computeArea(list)
}