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

Question

标题基本上说明了一切。我需要使用 Python 计算地球表面多边形内的面积。 计算地球表面上任意多边形包围的面积说有关它的一些信息，但在技术细节上仍然含糊不清：

如果你想用更多的方式来做到这一点 “GIS”风味，那么你需要选择 您所在地区的计量单位以及 找到一个合适的投影 保留区域（并非全部保留）。自从你 正在谈论计算 任意多边形，我会使用 像兰伯特方位角的东西 等面积投影。设置 投影的原点/中心 多边形的中心，项目 多边形到新坐标 系统，然后使用计算面积 标准平面技术。

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

它非常接近测量区域。对于更复杂的特征，您需要沿着顶点之间的边缘进行采样，以获得准确的值。上述有关日期线等的所有注意事项均适用。如果您只对面积感兴趣，则可以在投影之前将要素从日期变更线平移。

Answer 2

最简单的方法（在我看来）是将事物投影到（一个非常简单的）等面积投影中，并使用一种常用的平面技术来计算面积。

首先，如果您问这个问题，我将假设球形地球足够接近您的目的。如果没有，那么您需要使用适当的椭球体重新投影数据，在这种情况下，您将需要使用实际的投影库（现在一切都在幕后使用 proj4），例如与 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

或者简单地使用一个库： 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 4

也许有点晚了，但这里有一种不同的方法，使用吉拉德定理。它指出大圆多边形的面积是 R**2 乘以多边形之间的角度之和减去 (N-2)*pi，其中 N 是角的数量。

我认为这值得发布，因为它不依赖于 numpy 之外的任何其他库，并且它是一种与其他方法完全不同的方法。当然，这只适用于球体，所以应用到地球上时会存在一些不准确的情况。

首先，我定义一个函数来计算从点 1 沿大圆到点 2 的方位角：

import numpy as np
from numpy import cos, sin, arctan2

d2r = np.pi/180

def greatCircleBearing(lon1, lat1, lon2, lat2):
    dLong = lon1 - lon2

    s = cos(d2r*lat2)*sin(d2r*dLong)
    c = cos(d2r*lat1)*sin(d2r*lat2) - sin(lat1*d2r)*cos(d2r*lat2)*cos(d2r*dLong)

    return np.arctan2(s, c)

现在我可以使用它来求角度，然后求面积（在下面，当然应该指定经度和纬度，并且它们应该以正确的顺序排列。另外，应该指定球体的半径。）

N = len(lons)

angles = np.empty(N)
for i in range(N):

    phiB1, phiA, phiB2 = np.roll(lats, i)[:3]
    LB1, LA, LB2 = np.roll(lons, i)[:3]

    # calculate angle with north (eastward)
    beta1 = greatCircleBearing(LA, phiA, LB1, phiB1)
    beta2 = greatCircleBearing(LA, phiA, LB2, phiB2)

    # calculate angle between the polygons and add to angle array
    angles[i] = np.arccos(cos(-beta1)*cos(-beta2) + sin(-beta1)*sin(-beta2))

area = (sum(angles) - (N-2)*np.pi)*R**2

根据另一个回复中给出的科罗拉多坐标，以及地球半径 6371 公里，我得到的面积是 268930758560.74808

Answer 5

我知道十年后回答有一些优势，但对于今天看这个问题的人来说，提供更新的答案似乎是公平的。

pyproj 直接计算面积，无需调用 shapely：

# Modules:
from pyproj import Geod
import numpy as np

# Define WGS84 as CRS:
geod = Geod('+a=6378137 +f=0.0033528106647475126')

# Data for Colorado (no need to close the polygon):
coordinates = np.array([
[-102.05, 41.0], 
[-102.05, 37.0], 
[-109.05, 37.0], 
[-109.05, 41.0]])
lats = coordinates[:,1]
lons = coordinates[:,0]

# Compute:
area, perim = geod.polygon_area_perimeter(lons, lats)

print(abs(area))  # Positive is counterclockwise, the data is clockwise.

结果为：269154.54988400977 km2，或报告正确值（269601.367661 km2）的-0.17%。

Answer 6

您可以直接在球体上计算面积，而不是使用等积投影。

此外，根据 这个讨论，似乎吉拉德定理（sulkeh 的答案）在某些情况下并没有给出准确的结果，例如“30°月牙所包围的区域从极点到极点，以本初子午线和 30°E 为界”（请参阅​​此处）。

更精确的解决方案是直接在球体上执行线积分。从下面的比较可以看出，这种方法更加精确。

与所有其他答案一样，我应该提到我们假设地球是球形的警告，但我认为对于非关键目的这已经足够了。

Python 实现

这是一个使用线积分和格林定理的 Python 3 实现：

def polygon_area(lats, lons, radius = 6378137):
    """
    Computes area of spherical polygon, assuming spherical Earth. 
    Returns result in ratio of the sphere's area if the radius is specified.
    Otherwise, in the units of provided radius.
    lats and lons are in degrees.
    """
    from numpy import arctan2, cos, sin, sqrt, pi, power, append, diff, deg2rad
    lats = np.deg2rad(lats)
    lons = np.deg2rad(lons)

    # Line integral based on Green's Theorem, assumes spherical Earth

    #close polygon
    if lats[0]!=lats[-1]:
        lats = append(lats, lats[0])
        lons = append(lons, lons[0])

    #colatitudes relative to (0,0)
    a = sin(lats/2)**2 + cos(lats)* sin(lons/2)**2
    colat = 2*arctan2( sqrt(a), sqrt(1-a) )

    #azimuths relative to (0,0)
    az = arctan2(cos(lats) * sin(lons), sin(lats)) % (2*pi)

    # Calculate diffs
    # daz = diff(az) % (2*pi)
    daz = diff(az)
    daz = (daz + pi) % (2 * pi) - pi

    deltas=diff(colat)/2
    colat=colat[0:-1]+deltas

    # Perform integral
    integrands = (1-cos(colat)) * daz

    # Integrate 
    area = abs(sum(integrands))/(4*pi)

    area = min(area,1-area)
    if radius is not None: #return in units of radius
        return area * 4*pi*radius**2
    else: #return in ratio of sphere total area
        return area

我在 sphericalgeometry 包 那里。

数值比较

科罗拉多州将作为参考，因为之前的所有答案都是在其面积上进行评估的。其精确总面积为 104,093.67 平方英里（来自美国人口普查局，第 89 页，另请参阅此处 ），或269601367661平方米。我没有找到 USCB 实际方法的来源，但我认为它是基于对地面实际测量值进行求和，或使用 WGS84/EGM2008 进行精确计算。

Method                 | Author     | Result       | Variation from ground truth
--------------------------------------------------------------------------------
Albers Equal Area      | sgillies   | 268952044107 | -0.24%
Sinusoidal             | J. Kington | 268885360163 | -0.26%
Girard's theorem       | sulkeh     | 268930758560 | -0.25%
Equal Area Cylindrical | Jason      | 268993609651 | -0.22%
Line integral          | Yellows    | 269397764066 | **-0.07%**

结论：使用直接积分更精确。

性能

我没有对不同的方法进行基准测试，将纯 Python 代码与编译的 PROJ 投影进行比较是没有意义的。直观上需要更少的计算。另一方面，三角函数可能需要大量计算。

Answer 7

这是一个使用 basemap 而不是 pyproj 和 shapely 进行坐标转换的解决方案。这个想法与@sgillies 所建议的相同。请注意，我添加了第五点，以便路径成为闭环。

import numpy
from mpl_toolkits.basemap import Basemap

coordinates=numpy.array([
[-102.05, 41.0], 
[-102.05, 37.0], 
[-109.05, 37.0], 
[-109.05, 41.0],
[-102.05, 41.0]])

lats=coordinates[:,1]
lons=coordinates[:,0]

lat1=numpy.min(lats)
lat2=numpy.max(lats)
lon1=numpy.min(lons)
lon2=numpy.max(lons)

bmap=Basemap(projection='cea',llcrnrlat=lat1,llcrnrlon=lon1,urcrnrlat=lat2,urcrnrlon=lon2)
xs,ys=bmap(lons,lats)

area=numpy.abs(0.5*numpy.sum(ys[:-1]*numpy.diff(xs)-xs[:-1]*numpy.diff(ys)))
area=area/1e6

print area

结果为 268993.609651，单位为 km^2。

更新：底图已被弃用，因此您可能需要首先考虑替代解决方案。

Answer 8

因为地球是一个封闭的表面，所以在其表面绘制的封闭多边形会创建两个多边形区域。您还需要定义哪一个在内部，哪一个在外部！

大多数时候人们会处理小多边形，所以这是“显而易见的”，但是一旦你有海洋或大陆大小的东西，你最好确保你以正确的方式得到它。

另外，请记住，线可以通过两种不同的方式从 (-179,0) 到 (+179,0)。一个比另一个长得多。同样，大多数情况下，您会假设这是一条从 (-179,0) 到 (-180,0) 的线，即 (+180,0)，然后到 (+179,0)，但是天……不会的。

将经纬度视为简单的 (x,y) 坐标系，或者甚至忽略任何坐标投影都会发生扭曲和断裂的事实，可能会让您在球体上失败。

Answer 9

根据黄氏的说法，直接积分更为精确。

但是 Yellows 使用地球半径 = 6378 137m，即 WGS-84 椭球半长轴，而 Sulkeh 使用 6371 000 m。

在 Sulkeh 方法中使用半径 = 6378 137 m，得出 269533625893 平方米。

假设科罗拉多州面积的真实值（来自美国人口普查局）为 269601367661 平方米，则 Sulkeh 方法与真实值的偏差为：-0,025%，优于线积分法的 -0.07。

所以到目前为止，Sulkeh 的建议似乎更为精确。

为了能够对解进行数值比较，假设地球是球形的，所有计算都必须使用相同的地球半径。

Answer 10

这是一个 Python 3 实现，其中该函数将获取纬度和经度的元组对列表，并返回投影多边形中包含的面积。它使用 pyproj 来投影坐标，然后使用 Shapely 来查找任何投影多边形的面积

def calc_area(lis_lats_lons):

import numpy as np
from pyproj import Proj
from shapely.geometry import shape


lons, lats = zip(*lis_lats_lons)
ll = list(set(lats))[::-1]
var = []
for i in range(len(ll)):
    var.append('lat_' + str(i+1))
st = ""
for v, l in zip(var,ll):
    st = st + str(v) + "=" + str(l) +" "+ "+"
st = st +"lat_0="+ str(np.mean(ll)) + " "+ "+" + "lon_0" +"=" + str(np.mean(lons))
tx = "+proj=aea +" + st
pa = Proj(tx)

x, y = pa(lons, lats)
cop = {"type": "Polygon", "coordinates": [zip(x, y)]}

return shape(cop).area 

对于一组纬度/经度样本，它给出的面积值接近于测量的近似值

calc_area(lis_lats_lons = [(-102.05, 41.0),
 (-102.05, 37.0),
 (-109.05, 37.0),
 (-109.05, 41.0)])

，输出面积为 268952044107.4342 平方。山。