【问题标题】：How to calculate the area of a polygon on the earth's surface using python?如何使用python计算地球表面多边形的面积？
【发布时间】：2011-06-08 13:54:27
【问题描述】：

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

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

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

【问题讨论】：

如果选择区域保留投影，多边形的边将不再是直的

标签： python geometry geolocation geospatial

【解决方案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

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

【讨论】：

Strictly speaking，GeoJSON 应该有第五个结束点，[-102.05, 41.0]。有时，规范对这些事情有点含糊。
该区域结果的单位是什么？是平方米还是平方英寸？
如果其他人想知道，它以平方米为单位。你可以用谷歌搜索科罗拉多州的面积，得到 104,185 平方英里。将其转换为平方米得出大致相同的结果。
这对于可能不是简单四边形的多边形是否同样有效？
所以，基本上我可以使用任何等面积投影来转换坐标并使用格林定理来计算面积，对吗？如果是这样，那么可以替代地使用basemap 进行投影（例如'aea'、'cea'）。根据格林定理，单线 area=np.abs(0.5*np.sum(y[:-1]*np.diff(x) - x[:-1]*np.diff(y))) 让您无需使用 shapely 模块

【解决方案2】：

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

首先，如果您要问这个问题，我将假设球形地球足够接近您的目的。如果不是，那么您需要使用适当的椭球重新投影您的数据，在这种情况下，您将需要使用实际的投影库（现在一切都在幕后使用 proj4），例如 python 绑定到GDAL/OGR或（更友好的）pyproj。

但是，如果您对球形地球没问题，那么无需任何专门的库就可以很简单地做到这一点。

要计算的最简单的等面积投影是sinusoidal projection。基本上，你只需将纬度乘以一个纬度的长度，然后将经度乘以一个纬度的长度和纬度的余弦。

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

好的...现在我们要做的就是计算平面中任意多边形的面积。

有很多方法可以做到这一点。我将在这里使用probably the most common one。

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

无论如何，希望这将为您指明正确的方向......

【讨论】：

请注意，线条将被正弦投影扭曲。沿着点之间的大圆进行插值可能是明智的，这样您的投影多边形就不会失真。
@spacedman - 当然！我只是想展示一个简单的近似值。它并不打算完全准确。不过，只有当他试图计算一个边长非常且顶点很少的多边形时，多边形边的变形才有意义。
非常感谢您的重投影功能，您在迷路一天后救了我！旁注：重投影后，您还可以使用 Python 的 shapely 包计算面积，shapely.geometry.Polygon(list(zip(x, y))).area 以平方米为单位。

【解决方案3】：

也许有点晚了，但这里有一个不同的方法，使用 Girard 定理。它指出大圆多边形的面积是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)

现在我可以用它来找到角度，然后是面积（在下面，当然应该指定 lons 和 lats，并且它们应该按照正确的顺序。另外，应该指定球体的半径.)

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

【讨论】：

我试过你的方法，得到了一个负值。此外，它的绝对值（139013699.103）与我使用“网格”方法得到的值（809339.212）完全不同，该方法从矩形墨卡托网格中获取多边形内的所有网格单元，并总结网格单元区域。每个像元的面积是其纬向和经向间隔的乘积（从纬度、经度转换为公里）。
你能澄清你的评论吗？你提到的数字来自哪里？我再次尝试了我发布的代码，它给出了报告的结果......
我将它应用在使用matplotlib的contour函数找到的轮廓上，实际上我有几十个，每个人都用你的方法报告了一个负面区域。
好的。我仍然不明白你提到的数字来自哪里。如果您能更具体一点，我很乐意纠正任何错误。你能发布一个给出负面积的 lons/lats 的例子吗？
我在这里导出了一个轮廓坐标pastebin.com/2fMkFgvu也许你可以试一试。我注意到的一件事是，如果我通过每 5 个点 (x=x[::5]; y=y[::5]) 对 135 点轮廓进行二次采样，我可以获得一个接近“真实”值的正值。

【解决方案4】：

或者干脆使用库：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

...以平方米为单位返回面积。

【讨论】：

我喜欢这种方法让我可以直接使用纬度/经度点，并考虑地球的曲率，而无需将点投影到笛卡尔网格上。我将这种方法和标准格林定理的结果与谷歌上不同州的区域进行了比较，发现这个库给出了最准确的结果。谢谢！
AFAI 明白，这个库正在计算测地线面积，假设地球是一个完美的球体，半径为 defined 为 WGS84_RADIUS = 6378137 请记住这一点。

【解决方案5】：

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

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。

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

【讨论】：

为了澄清，用于计算面积的方程的来源是什么？
@blaylockbk 这是Green's theorem。

【解决方案6】：

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

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

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

将 lat-long 视为简单的 (x,y) 坐标系，甚至忽略任何坐标投影都会出现扭曲和中断的事实，可能会让您在球体上大失所望。

【讨论】：

非常真实！我只是在我的回答中给出了一个最基本情况的简单例子......它甚至没有远程尝试处理穿越本初子午线（0-360）或反子午线（-180到180）。
是的。当地球是平的时候，生活就轻松多了。

【解决方案7】：

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

此外，根据this discussion，似乎 Girard 定理（苏尔克的答案）在某些情况下并不能给出准确的结果，例如“由从极到极的 30º 月牙包围并以本初子午线和30ºE"（见here）。

更精确的解决方案是直接在球体上进行线积分。下面的比较表明这种方法更精确。

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

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 包there 中编写了一个更明确的版本（以及更多的参考和待办事项...）。

数值比较

科罗拉多州将作为参考，因为之前的所有答案都在其区域进行了评估。其精确的总面积为 104,093.67 平方英里（来自US Census Bureau，第 89 页，另见here），或 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 投影进行比较没有意义。直观上需要更少的计算。另一方面，三角函数可能是计算密集型的。

【讨论】：

我现在更喜欢这个解决方案而不是我自己的解决方案，因为它不依赖任何额外的包。底图已被弃用。
科罗拉多州的边界不是一个完美的正方形，在边界建立时，由于测量不完善，西侧有一个小“凸起”。人口普查区可能说明了这一点，这也是为什么这里所有地区都略微低估了它的原因。

【解决方案8】：

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

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%。

【讨论】：

【解决方案9】：

根据 Yellows 的说法，直接积分更精确。

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

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

假设科罗拉多地区的真实值（来自美国人口普查局）是 269601367661 平方米，那么 Sulkeh 方法与地面实况的变化为：-0,025%，优于线积分法的 -0.07。

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

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

【讨论】：

【解决方案10】：

这是一个 Python 3 实现，其中函数将获取 lats 和 long 的元组对列表，并返回投影多边形中包含的区域。它使用 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 Sq。山。

【讨论】：