如何使用python计算地球表面上多边形的面积？

Question

标题基本上说明了一切。我需要使用Python计算地球表面上的多边形内部的面积。计算地球表面上任意多边形所围成的面积虽然可以说明这一点，但在技术细节上仍然含糊不清：

如果要使用更“
GIS”的样式来执行此操作，则需要为您的区域选择一个度量单位，并找到保留该区域的适当投影（并非全部如此）。由于您正在谈论计算任意多边形，因此我将使用类似Lambert
Azimuthal等面积投影的方法。将投影的原点/中心设置为多边形的中心，将多边形投影到新的坐标系，然后使用标准平面技术计算面积。

那么，如何在Python中做到这一点？

Answer 1

假设您以GeoJSON格式表示了科罗拉多州

{"type": "Polygon", 
 "coordinates": [[
   [-102.05, 41.0], 
   [-102.05, 37.0], 
   [-109.05, 37.0], 
   [-109.05, 41.0]
 ]]}

所有坐标均为经度，纬度。您可以使用pyproj投影坐标，并使用Shapely查找任何投影多边形的面积：

co = {"type": "Polygon", "coordinates": [
    [(-102.05, 41.0),
     (-102.05, 37.0),
     (-109.05, 37.0),
     (-109.05, 41.0)]]}
lon, lat = zip(*co['coordinates'][0])
from pyproj import Proj
pa = Proj("+proj=aea +lat_1=37.0 +lat_2=41.0 +lat_0=39.0 +lon_0=-106.55")

这是一个以关注区域为中心并包围其关注区域的等面积投影。现在制作新的投影GeoJSON表示形式，将其转换为Shapely几何对象，并采用以下区域：

x, y = pa(lon, lat)
cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
from shapely.geometry import shape
shape(cop).area  # 268952044107.43506

这是与被调查区域非常接近的近似值。对于更复杂的特征，您需要沿顶点之间的边缘进行采样，以获取准确的值。上面有关日期变更线等的所有警告均适用。如果您仅对区域感兴趣，则可以在投影之前将特征从日期线移开。