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

Posted 2021-04-18

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

如果你想用更“GIS”的味道来做这件事，那么你需要为你的区域选择一个度量单位，找到一个保留区域的适当投影（不是所有的）。既然你在谈论计算任意多边形，我会使用像Lambert Azimuthal Equal Area投影这样的东西。将投影的原点/中心设置为多边形的中心，将多边形投影到新坐标系，然后使用标准平面技术计算面积。

那么，我如何在Python中执行此操作？

答案

假设您以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

它与调查区域非常接近。对于更复杂的特征，您需要沿顶点之间的边缘采样，以获得准确的值。以上关于日期线等的所有警告均适用。如果您只对区域感兴趣，可以在投影前将您的功能从日期线转换出来。

另一答案

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

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

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

另一答案

也许有点晚了，但这是一种不同的方法，使用吉拉德定理。它指出，大圆形多边形的面积是多边形之间的角度之和减去（N-2）* pi的R ** 2倍，其中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 km，我得到的区域是268930758560.74808

另一答案

由于地球是封闭的表面，因此在其表面上绘制的闭合多边形会产生两个多边形区域。你还需要定义哪一个在里面，哪个在外面！

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

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

像简单的（x，y）坐标系一样处理lat-long，甚至忽略任何坐标投影都会产生扭曲和断裂的事实，可能会让你在球体上失败。

另一答案

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

...返回平方米的面积。

另一答案

这是一个使用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为单位。