标题基本上都说明了一切。我需要使用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为单位。