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

Question

标题基本上说明了一切。我需要使用 Python 计算地球表面多边形内的面积。 Calculating area enclosed by arbitrary polygon on Earth's surface对此有所说明，但在技术细节上仍含糊不清:

If you want to do this with a more "GIS" flavor, then you need to select an unit-of-measure for your area and find an appropriate projection that preserves area (not all do). Since you are talking about calculating an arbitrary polygon, I would use something like a Lambert Azimuthal Equal Area projection. Set the origin/center of the projection to be the center of your polygon, project the polygon to the new coordinate system, then calculate the area using standard planar techniques.

那么，我如何在 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

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