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

Question

标题基本上已经说明了一切。我需要用Python计算地球表面一个多边形内部的面积。计算地球表面任意多边形包围的面积提到了一些相关内容，但在技术细节上说得比较模糊：

如果你想用更“地理信息系统（GIS）”的方式来做这个，那么你需要选择一个面积的单位，并找到一个合适的投影方式来保持面积（并不是所有的投影都能做到这一点）。因为你要计算的是一个任意的多边形，我建议使用类似于“兰伯特方位等面积投影”的方法。将投影的原点/中心设置为你的多边形的中心，然后把这个多边形投影到新的坐标系统中，最后用标准的平面技术来计算面积。

那么，我该如何在Python中实现这个呢？

Answer 1

或者直接使用一个库： https://github.com/scisco/area

from area import area
>>> obj = {'type':'Polygon','coordinates':[[[-180,-90],[-180,90],[180,90],[180,-90],[-180,-90]]]}
>>> area(obj)
511207893395811.06

...返回的结果是以平方米为单位的面积。

Answer 2

我觉得最简单的方法就是把东西投影到一个非常简单的等面积投影上，然后用一些常见的平面技术来计算面积。

首先，我假设你问这个问题是因为你觉得地球是个球体，这样就足够了。如果不是，那你需要用合适的椭球体重新投影你的数据，这种情况下你可能需要用到实际的投影库（现在很多东西都是在后台使用proj4的），比如Python的GDAL/OGR或者更友好的pyproj。

不过，如果你可以接受地球是个球体，那其实不需要任何专业的库就能很简单地做到这一点。

最简单的等面积投影是正弦投影。基本上，你只需要把纬度乘以一个纬度的长度，把经度乘以一个纬度的长度再乘以纬度的余弦值。

def reproject(latitude, longitude):
    """Returns the x & y coordinates in meters using a sinusoidal projection"""
    from math import pi, cos, radians
    earth_radius = 6371009 # in meters
    lat_dist = pi * earth_radius / 180.0

    y = [lat * lat_dist for lat in latitude]
    x = [long * lat_dist * cos(radians(lat)) 
                for lat, long in zip(latitude, longitude)]
    return x, y

好了……现在我们只需要计算一个任意多边形在平面上的面积。

有很多方法可以做到这一点。我将使用可能是最常见的一种。

def area_of_polygon(x, y):
    """Calculates the area of an arbitrary polygon given its verticies"""
    area = 0.0
    for i in range(-1, len(x)-1):
        area += x[i] * (y[i+1] - y[i-1])
    return abs(area) / 2.0

希望这能给你指明方向……

Answer 3

假设你有一个关于科罗拉多州的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

这个计算结果与实际测量的面积非常接近。对于更复杂的形状，你需要在边缘和顶点之间取样，以获得更准确的值。之前提到的关于国际日期线等的注意事项仍然适用。如果你只关心面积，可以在投影之前把你的特征移离国际日期线。