10

我正在使用 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，有没有办法将整个形状优美的多边形对象转换为另一个投影，然后计算面积？

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

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

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

有更好的想法吗？

4

3 回答 3

11

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

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

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

需要 numpy、dateutil、pytz、pyparsing、6 和可选的pillow、pycairo、tornado、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)
于 2014-05-19T09:12:49.830 回答
9

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

于 2020-10-02T01:15:44.040 回答
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]])

多边形=多边形（np.弧度（p））

poly.area = 4.468737548271707e-07

poly.area*6370**2 =18.132751662246623

于 2014-06-12T17:42:47.473 回答