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

Question

我使用Python 3.4和shapely 1.3.2创建一个多边形对象的长/lat坐标对列表，我转换成一个众所周知的文本字符串，以便解析它们。这样的多边形可能看起来像：

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，有没有一种方法可以将整个形状优美的多边形对象转换为另一个投影，然后计算面积？

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

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

这意味着将每个多边形对象拆分为其外圈和内圈（如果存在的话），获取坐标，将每对坐标转换为另一个投影并重建多边形对象，然后计算其面积（它是什么单位？）。这看起来是一个解决方案，但不太实际。

还有更好的主意吗？

Answer 1

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

p=np。数组（[-116.904,43.371]、-116.823,43.389]、-116.895,43.407]、-116.908,43.375]、-116.904,43.371]）

Poly=多边形（np.radians（p））

聚。面积=4.468737548271707e-07

聚。面积*6370**2=18.132751662246623

Answer 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（）仅用于返回正区域。根据多边形的缠绕方向，可能会返回负区域。

Answer 3

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

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

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

需要numpy、dateutil、pytz、pyparsing、six，还可以选择枕头、pycairo、tornado、wxpython、pyside、pyqt、ghostscript、miktex、ffmpeg、mencoder、avconv或imagemagick。

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

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

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.19135556714567974.064658936）

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

编辑：我努力不只是变换单个坐标，而是变换整个几何体对象，比如多边形，而这些对象已经被赋予了形状优美的对象，而不是通过它们的原始坐标。这意味着要编写大量的代码

• 得到多边形外环的坐标
• 确定多边形是否有孔，如果有，分别处理每个孔
• 变换外环和任意孔的每对坐标
• 把整个东西放回一起，并创建一个多边形对象与投影坐标
• 这只适用于多边形...为多边形和几何集合添加更多的循环

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

当设置投影贴图时，例如如上所述：

m=Basemap（宽度=1，高度=1，分辨率=l'，投影=laea'，纬度=50，纬度0=50，经度0=-107。）

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

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