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

1%的误差似乎有点高,因为只有你的近似值。你是在比较实际的测量结果还是一些理想的计算结果?请记住,GPS中也存在可能导致故障的错误。

如果你想要一个更精确的方法来做这件事,有一个很好的答案 this 问题。如果你想要更快的方法,你可以使用wgs84大地水准面而不是你的参考球体来转换为笛卡尔坐标(ecef)。这里是 wiki link 为了那个转变。

我正在修改谷歌地图,以便用户可以计算面积 通过单击顶点来创建多边形。它没有给出正确的 直到我确定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、lonchor和latanchor以度为单位。现在工作得很有魅力。

我在网上查了各种多边形面积公式(或代码),但没有发现任何一个好的或容易实现的。

现在我编写了一个代码片段来计算在地球表面绘制的多边形的面积。多边形可以有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;
   }
   

基于风险Pathak的解决方案,这里是SQL(Redshift)计算区域的解决方案 GeoJSON multipolygons (假设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
;