在 Shapely 中以平面单位（例如平方米）计算多边形面积

Posted 2023-02-14

【中文标题】在 Shapely 中以平面单位（例如平方米）计算多边形面积

【英文标题】：Calculate Polygon area in planar units (e.g. square-meters) in Shapely

【发布时间】：2014-07-05 00:43:35

【问题描述】：

我正在使用 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 不处理任何投影并在 carthesian 空间中实现所有几何对象，因此在该多边形上调用 area 方法，如下所示：

poly.area

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

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

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

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

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

有更好的想法吗？

【参考方案1】：

好的，我终于用 matplotlib 库的 Basemap 工具包做到了。我会解释它是如何工作的，也许这对某人有帮助。

1。 在您的系统上下载并安装 matplotlib 库。 http://matplotlib.org/downloads.html

对于 Windows 二进制文件，我建议使用此页面： http://www.lfd.uci.edu/~gohlke/pythonlibs/#matplotlib 请注意以下提示：

需要 numpy、dateutil、pytz、pyparsing、6 和可选的枕头， pycairo，龙卷风，wxpython，pyside，pyqt，ghostscript，miktex，ffmpeg， mencoder、avconv 或 imagemagick。

因此，如果您的系统上尚未安装，您必须下载并安装 numpy、dateutil、pytz、pyparsing 和六个，以便 matplotlib 正常工作（对于 Windows 用户：所有这些都可以在页面上找到，对于 Linux 用户，python 包管理器“pip”应该可以解决问题）。

2。 下载并安装 matplotlib 的“底图”工具包。无论是从 http://matplotlib.org/basemap/users/installing.html 或者对于 Windows 二进制文件也来自这里： http://www.lfd.uci.edu/~gohlke/pythonlibs/#basemap

3。 在 Python 代码中进行投影：

#Import necessary libraries
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

#Coordinates that are to be transformed
long = -112.367
lat = 41.013

#Create a basemap for your projection. Which one you use is up to you.
#Some examples can be found at http://matplotlib.org/basemap/users/mapsetup.html
m = Basemap(projection='sinu',lon_0=0,resolution='c')

#Project the coordinates:
projected_coordinates = m(long,lat)

输出：

投影坐标 (10587117.191355567, 14567974.064658936)

就这么简单。现在，当使用投影坐标构建一个带有 shapely 的多边形，然后通过 shapely 的面积方法计算面积时，您将获得以平方米为单位的面积（根据您使用的投影）。要获得平方公里，除以 10^6。

编辑：我努力不只是变换单个坐标，而是像多边形这样的整个几何对象，因为这些对象已经作为匀称的对象而不是通过它们的原始坐标给出。这意味着要编写大量代码来

获取多边形外环坐标 确定多边形是否有孔，如果有，分别处理每个孔 变换外环和任意孔的每一对坐标 将整个东西重新组合在一起，并使用投影坐标创建一个多边形对象 这仅适用于多边形...为多多边形和几何集合添加更多循环

然后我偶然发现了这部分匀称的文档，它使事情变得容易多了： http://toblerity.org/shapely/manual.html#shapely.ops.transform

设置投影图时，比如上面做的：

m = Basemap(width=1,height=1, resolution='l',projection='laea',lat_ts=50,lat_0=50,lon_0=-107.)

然后，可以通过以下方式使用此投影变换任何形状匀称的几何对象：

from shapely.ops import transform
projected_geometry = transform(m,geometry_object)

【讨论】：

另见投影变换示例here。请注意，面积计算的准确性取决于投影。例如，如果数据都很好地放在一个中，则使用 UTM 区域。

【参考方案2】：

计算测地线面积，它非常准确，只需要一个椭球体（不是投影）。这可以使用 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() 用于仅返回正区域。根据多边形的缠绕方向，可能会返回负面积。

【讨论】：

不错。另外，来自the pyproj documentation：对于多边形，孔应该使用与外部相反的遍历（如果外部是逆时针，孔/内部应该是顺时针）。

【参考方案3】：

转换为弧度并假设地球是半径为 6370Km 的完美球体：

p = np.array([[-116.904,43.371],[-116.823, 43.389],[-116.895,43.407],[-116.908,43.375],[-116.904,43.371]])

poly = 多边形(np.radians(p))

poly.area =4.468737548271707e-07

poly.area*6370**2 =18.132751662246623

【讨论】：

这个非常错误，它假设多边形非常小&在赤道上。否则，数学比这复杂得多！ （@DST，我建议你删除或改进这篇文章，因为它具有很大的误导性。）