python - 在 Shapely 中以平面单位(例如平方米)计算多边形面积

Question

我正在使用 Python 3.4 和 shapely 1.3.2 从长/纬度坐标对列表中创建一个多边形对象，我将其转换为众所周知的文本字符串以便解析它们。这样的多边形可能看起来像:

POLYGON ((-116.904 43.371, -116.823 43.389, -116.895 43.407, -116.908 43.375, -116.904 43.371))

由于 shapely 不处理任何投影并在笛卡尔空间中实现所有几何对象，因此在该多边形上调用 area 方法，如下所示:

poly.area

以平方度为单位给出该多边形的面积。要获得像平方米这样的平面单位的面积，我想我必须使用不同的投影(哪个投影？)来转换多边形的坐标。

我多次读到 pyproj 库应该提供执行此操作的方法。使用pyproj，有没有办法将一个完整的形状多边形对象转换为另一个投影，然后计算面积？

我用我的多边形做了一些其他事情(不是你现在想的那样)，只有在某些情况下，我才需要计算面积。

到目前为止，我只找到了这个例子: http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/

这意味着将每个 Polygon 对象拆分为其外环和内环(如果存在)，获取坐标，将每对坐标转换为另一个投影并重建 Polygon 对象，然后计算其面积(无论如何单位是什么？)。这看起来像一个解决方案，但不是很实用。

有什么更好的主意吗？

Answer 1

计算测地面积，非常准确，只需要一个椭球体(不是投影)。这可以使用 pyproj 2.3.0 或更高版本来完成。

from pyproj import Geod
from shapely import wkt

# specify a named ellipsoid
geod = Geod(ellps="WGS84")

poly = wkt.loads('''\
POLYGON ((-116.904 43.371, -116.823 43.389, -116.895 43.407,
-116.908 43.375, -116.904 43.371))''')

area = abs(geod.geometry_area_perimeter(poly)[0])

print('# Geodesic area: {:.3f} m^2'.format(area))

# # Geodesic area: 13205034.647 m^2

abs() 用于仅返回正区域。根据多边形的缠绕方向，可能会返回负面积。