使用Shapely计算平面单位（例如平方米）中的多边形面积

Question

我正在使用Python 3.4和shapely 1.3.2将一系列经纬度坐标对转换为一个well-known-text字符串，以便进行解析并创建一个多边形对象。这样的多边形可能如下所示：

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，是否有一种方法可以将整个Shapely Polygon对象转换为另一个投影，并计算其面积？

我对我的多边形进行一些其他操作（不是你现在想的那个），只有在某些情况下，我才需要计算面积。

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

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

有更好的想法吗？

- Dirk

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() 用于返回仅为正面积。根据多边形的绕向方向，可能会返回负数面积。

- Mike T

---

不错。除了pyproj文档之外：对于多边形，孔应该使用与外部相反的遍历顺序（如果外部是逆时针，则孔/内部应该是顺时针）。 - Arjan

非常有帮助。作为一个新手，我一直在想为什么.area返回0。 - Kevin-Prichard

Answer 2

好的，我最终使用了matplotlib库的Basemap工具包。我将解释它的工作原理，也许对某人会有所帮助。

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

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

需要numpy、dateutil、pytz、pyparsing、six以及可选的pillow、pycairo、tornado、wxpython、pyside、pyqt、ghostscript、miktex、ffmpeg、mencoder、avconv或imagemagick。

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

2.从 http://matplotlib.org/basemap/users/installing.html 或者对于Windows二进制文件，也可以从这里下载并安装“basemap”工具包： 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的面积方法计算面积，您将得到平方米为单位的面积（根据所用的投影）。要获得平方千米，请除以10^6。

编辑： 我很难只转换单个坐标，而是像多边形那样处理整个几何对象，特别是当它们已经作为Shapely对象而不是通过其原始坐标给出时。这意味着要编写大量代码来

• 获取多边形外环的坐标
• 确定多边形是否具有孔，如果是，单独处理每个孔
• 转换外环和任何孔的每对坐标
• 重新组合所有内容并使用投影后的坐标创建多边形对象
• 对于多面体和几何集，需要添加更多的循环...

然后我偶然发现了Shapely文档中的这一部分，它使事情变得容易得多： 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.)

然后，可以使用此投影通过以下方式转换任何Shapely几何对象：

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

- Dirk

---

请参考此处的投影变换示例：https://dev59.com/S3vaa4cB1Zd3GeqPGbGw#21420950。请注意，区域计算的准确性取决于投影方式。例如，如果数据都很好地适合一个UTM区域，请使用它。 - Mike T

Answer 3

将角度转换为弧度，并假设地球是半径为6370公里的完美球体：

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

poly = Polygon(np.radians(p))

poly.area
# 4.468737548271707e-07

poly.area*6370**2
# 18.132751662246623

- DST

---

这个有很大的错误，它假定多边形非常小且位于赤道上。否则，数学计算将比这复杂得多！ (@DST, 我建议你删除或改进这篇文章，因为它极具误导性。) - zabop