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

Question

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

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

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

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

有更好的主意吗？

Answer 1

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

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

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

Requires numpy, dateutil, pytz, pyparsing, six, and optionally pillow, pycairo, tornado, wxpython, pyside, pyqt, ghostscript, miktex, ffmpeg, mencoder, avconv, or imagemagick.

因此，如果您的系统上还没有安装，那么您必须下载并安装numpy、dateutil、pytz、pyparsing和6，matplotlib才能正常工作（对于Windows用户：所有这些都可以在页面上找到，对于Linux用户，python包管理器“pip”应该能做到这一点）。

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

三。 在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)

输出：

projected_coordinates (10587117.191355567, 14567974.064658936)

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

编辑：我努力不只是转换单个坐标，而是转换像多边形这样的整个几何体对象，因为这些对象已经被指定为形状对象，而不是通过它们的原始坐标。这意味着要写很多代码

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

然后我偶然发现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.)

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

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

Answer 2

转换成弧度，假设地球是半径为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 = Polygon(np.radians(p))

poly.area =4.468737548271707e-07

poly.area*6370**2 =18.132751662246623