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

问题
标题基本上说明了一切。我需要在 Python 中计算地球表面多边形内的面积。计算由任意多边形包围的地球表面的面积说明了这一点，但在技术细节上仍然含糊不清：

那么，我怎样才能在 Python 中实现这一点呢？

回答
假设您有一个以 GeoJSON 格式表示的科罗拉多州

1. {"type": "Polygon",
   
2. "coordinates": [[
   
3.    [-102.05, 41.0],
   
4.    [-102.05, 37.0],
   
5.    [-109.05, 37.0],
   
6.    [-109.05, 41.0]
   
7. ]]}

复制代码

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

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

复制代码

这是一个以感兴趣区域为中心并用括号括起来的等面积投影。现在创建一个新的投影 GeoJSON 表示，将其转换为形状几何对象，并获取面积：

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

复制代码

这是一个非常接近的测量区域。对于更复杂的特征，您需要在顶点之间沿边进行采样，以获得准确的值。以上关于日期线等的所有警告都适用。如果您只对该区域感兴趣，则可以在投影之前将要素转换为日期线之外。