Python/shapely:以平面单位计算多边形面积(例如平方米)

Question

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

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对象,然后计算其面积(无论如何它是什么单位？).这看起来像一个解决方案,但不太实用.

有更好的想法吗？

Answer 1

好吧,我终于使用了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以及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

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)

输出:

projected_coordinates(10587117.191355567,14567974.064658936)

就那么简单.现在,当使用投影坐标来构建具有匀称的多边形,然后通过匀称的面积方法计算面积时,您将获得以平方米为单位的面积(根据您使用的投影).要达到平方公里,除以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)

Answer 2

Calculate a geodesic area, which is very accurate and only requires an ellipsoid (not a projection). This can be done with pyproj 2.3.0 or later.

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() is used to return only positive areas. A negative area may be returned depending on the winding direction of the polygon.