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

How to calculate the area of a polygon on the earth's surface using python?

提问by andreas-h

The title basically says it all. I need to calculate the area inside a polygon on the Earth's surface using Python. Calculating area enclosed by arbitrary polygon on Earth's surfacesays something about it, but remains vague on the technical details:

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

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.

如果你想用更“GIS”的风格来做到这一点，那么你需要为你的区域选择一个度量单位，并找到一个合适的投影来保留区域（不是全部都这样做）。由于您在谈论计算任意多边形，因此我会使用诸如兰伯特方位角等面积投影之类的东西。将投影的原点/中心设置为多边形的中心，将多边形投影到新坐标系，然后使用标准平面技术计算面积。

So, how do I do this in Python?

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

采纳答案by sgillies

Let's say you have a representation of the state of Colorado in GeoJSON format

假设您以 GeoJSON 格式表示科罗拉多州

{"type": "Polygon", 
 "coordinates": [[
   [-102.05, 41.0], 
   [-102.05, 37.0], 
   [-109.05, 37.0], 
   [-109.05, 41.0]
 ]]}

All coordinates are longitude, latitude. You can use pyprojto project the coordinates and Shapelyto find the area of any projected polygon:

所有坐标都是经度，纬度。您可以使用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")

That's an equal area projection centered on and bracketing the area of interest. Now make new projected GeoJSON representation, turn into a Shapely geometric object, and take the area:

这是一个以感兴趣区域为中心并包围感兴趣区域的等面积投影。现在制作新的投影 GeoJSON 表示，变成一个 Shapely 几何对象，并取区域：

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

It's a very close approximation to the surveyed area. For more complex features, you'll need to sample along the edges, between the vertices, to get accurate values. All caveats above about datelines, etc, apply. If you're only interested in area, you can translate your feature away from the dateline before projecting.

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

回答by Joe Kington

The easiest way to do this (in my opinion), is to project things into (a very simple) equal-area projection and use one of the usual planar techniques for calculating area.

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

First off, I'm going to assume that a spherical earth is close enough for your purposes, if you're asking this question. If not, then you need to reproject your data using an appropriate ellipsoid, in which case you're going to want to use an actual projection library (everything uses proj4 behind the scenes, these days) such as the python bindings to GDAL/OGRor (the much more friendly) pyproj.

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

However, if you're okay with a spherical earth, it quite simple to do this without any specialized libraries.

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

The simplest equal-area projection to calculate is a sinusoidal projection. Basically, you just multiply the latitude by the length of one degree of latitude, and the longitude by the length of a degree of latitude and the cosine of the latitude.

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

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

Okay... Now all we have to do is to calculate the area of an arbitrary polygon in a plane.

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

There are a number of ways to do this. I'm going to use what is probably the most common onehere.

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

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

Hopefully that will point you in the right direction, anyway...

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

回答by Spacedman

Because the earth is a closed surface a closed polygon drawn on its surface creates TWOpolygonal areas. You also need to define which one is inside and which is outside!

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

Most times people will be dealing with small polygons, and so it's 'obvious' but once you have things the size of oceans or continents, you better make sure you get this the right way round.

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

Also, remember that lines can go from (-179,0) to (+179,0) in two different ways. One is very much longer than the other. Again, mostly you'll make the assumption that this is a line that goes from (-179,0) to (-180,0) which is (+180,0) and then to (+179,0), but one day... it won't.

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

Treating lat-long like a simple (x,y) coordinate system, or even neglecting the fact that any coordinate projection is going to have distortions and breaks, can make you fail big-time on spheres.

将 lat-long 视为一个简单的 (x,y) 坐标系，甚至忽略任何坐标投影都会产生扭曲和断裂的事实，都会使您在球体上失败。

回答by sulkeh

A bit late perhaps, but here is a different method, using Girard's theorem. It states that the area of a polygon of great circles is R**2 times the sum of the angles between the polygons minus (N-2)*pi where N is number of corners.

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

I thought this would be worth posting, since it doesn't rely on any other libraries than numpy, and it is a quite different method than the others. Of course, this only works on a sphere, so there will be some inaccuracy when applying it to the Earth.

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

First, I define a function to compute the bearing angle from point 1 along a great circle to point 2:

首先，我定义了一个函数来计算从点 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)

Now I can use this to find the angles, and then the area (In the following, lons and lats should of course be specified, and they should be in the right order. Also, the radius of the sphere should be specified.)

现在我可以用它来找到角度，然后是面积（在下面，当然应该指定 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

With the Colorado coordinates given in another reply, and with Earth radius 6371 km, I get that the area is 268930758560.74808

使用另一个回复中给出的科罗拉多坐标，地球半径为 6371 公里，我得到该区域为 268930758560.74808

回答by Ikar Pohorsky

Or simply use a library: https://github.com/scisco/area

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

...returns the area in square meters.

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

回答by Jason

Here is a solution that uses basemap, instead of pyprojand shapely, for the coordinate conversion. The idea is the same as suggested by @sgillies though. NOTE that I've added the 5th point so that the path is a closed loop.

这是一个使用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

The result is 268993.609651 in km^2.

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

回答by Yellows

You can compute the area directly on the sphere, instead of using an equal-area projection.

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

Moreover, according to this discussion, it seems that Girard's theorem (sulkeh's answer) does not give accurate results in certain cases, for example "the area enclosed by a 30o lune from pole to pole and bounded by the prime meridian and 30oE" (see here).

此外，根据此讨论，似乎 Girard 定理（sulkeh 的答案）在某些情况下并没有给出准确的结果，例如“由一个 30° 的月线从极点到另一个极点包围，并以本初子午线和 30oE 为界的区域”（参见在这里）。

A more precise solution would be to perform line integral directly on the sphere. The comparison below shows this method is more precise.

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

Like all other answers, I should mention the caveat that we assume a spherical earth, but I assume that for non-critical purposes this is enough.

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

Python implementation

Python 实现

Here is a Python 3 implementation which uses line integral and Green's theorem:

这是一个使用线积分和格林定理的 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

I wrote a somewhat more explicit version (and with many more references and TODOs...) in the sphericalgeometrypackage there.

我在那里的球面几何包中写了一个更明确的版本（还有更多的参考资料和待办事项......）。

Numerical Comparison

数值比较

Colorado will be the reference, since all previous answers were evaluated on its area. Its precise total area is 104,093.67 square miles (from the US Census Bureau, p. 89, see also here), or 269601367661 square meters. I found no source for the actual methodology of the USCB, but I assume it is based on summing actual measurements on ground, or precise computations using WGS84/EGM2008.

科罗拉多州将作为参考，因为所有以前的答案都在其区域内进行了评估。它的精确总面积为 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%**

Conclusion: using direct integral is more precise.

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

Performance

表现

I have not benchmarked the different methods, and comparing pure Python code with compiled PROJ projections would not be meaningful. Intuitively less computations are needed. On the other hand, trigonometric functions may be computationally intensive.

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

回答by Styx

According to Yellows' assertion, direct integral is more precise.

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

But Yellows use an earth radius = 6378 137m, which is the WGS-84 ellipsoid, semi-major axis, while Sulkeh use 6371 000 m.

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

Using a radius = 6378 137 m in the Sulkeh' method, gives 269533625893 square meters.

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

Assuming that the true value of Colorado area (from the US Census Bureau) is 269601367661 square meters then the variation from the ground truth of Sulkeh' method is : -0,025%, better than -0.07 with the Line integral method.

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

So Sulkeh' proposal seems to be the more precise so far.

因此，到目前为止，Sulkeh 的提议似乎更为精确。

In order to be able to make a numerical comparison of the solutions, with the assumption of a spherical Earth, all calculations must use the same terrestrial radius.

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

回答by Kunal Verma

Here is a Python 3 implementation where the function would take a list of tuple-pairs of lats and longs and would return the area enclosed in the projected polygon.It uses pyproj to project the coordinates and then Shapely to find the area of any projected polygon

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

For a sample set of lats/longs, it gives an area value close to the surveyed approximation value

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

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

Which outputs an area of 268952044107.4342 Sq. Mts.

其输出面积为 268952044107.4342 Sq。山