那么，我该如何在Python中实现这个功能呢？如果您想要更具“GIS”风格的方法，则需要选择用于测量面积的单位，并找到一个保留面积（并非所有投影都能）的适当投影。由于您要计算任意多边形，因此建议使用诸如Lambert等面积方位投影。将投影的中心/原点设置为多边形的中心，将多边形投影到新坐标系统中，然后使用标准平面技术计算其面积。
假设你拥有一个以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
这是与测量区域非常接近的近似值。 对于更复杂的功能，您需要沿边缘，顶点之间进行采样，以获得精确的值。 所有关于日期线等的警告都适用。 如果您只对面积感兴趣，可以在投影之前将要素移开日期线。
[-102.05, 41.0]。规范有时候对这些事情描述得不太明确。 - Brad Kocharea=np.abs(0.5*np.sum(y[:-1]*np.diff(x) - x[:-1]*np.diff(y))) 就可以让你不需要使用 shapely 模块。 - Jasondef 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
shapely包来计算面积，方法是shapely.geometry.Polygon(list(zip(x, y))).area，以得到平方米为单位的面积。 - Ali或者简单使用一个库：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
...以平方米为单位返回面积。
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
139013699.103) 与我使用“网格”方法得到的值 (809339.212) 相差很大。该方法是从一个矩形墨卡托网格中获取多边形内所有网格单元，并将网格单元面积相加。每个单元格的面积是其纵向和经向间隔的乘积（从纬度、经度转换为公里）。 - Jasoncontour函数在一个轮廓上应用了它，事实上我有几十个这样的轮廓，每个都使用您的方法报告负面积。 - Jasonx=x[::5]; y=y[::5]），我可以得到一个接近“真实”值的正值。 - Jasonbasemap而不是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平方千米。
更新：Basemap已被弃用，因此您可能需要首先考虑其他替代方案。
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平方米。我没有找到美国人口普查局实际方法的来源，但我认为它基于在地面上测量或使用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投影进行比较是没有意义的。 直观地说，需要进行更少的计算。 另一方面，三角函数可能会占用大量计算资源。
结论：使用直接积分更加精确。
我知道，回答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.
根据Yellows的说法，直接积分更精确。
但是Yellows使用的地球半径为6378 137米，这是WGS-84椭球体的半长轴，而Sulkeh使用的是6371 000米。
在Sulkeh的方法中使用半径为6378 137米，则得到269533625893平方米。
假设科罗拉多州面积的真实值（来自美国人口普查局）为269601367661平方米，则Sulkeh方法与实际值的偏差为：-0.025％，比线积分法的-0.07％更好。
因此，Sulkeh的提议似乎是迄今为止更精确的。
为了能够对解决方案进行数字比较，并假设地球是一个球体，所有计算必须使用相同的地球半径。
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平方米。