我正在使用 Python 3.4 和 shapely 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 不处理任何投影并实现笛卡尔空间中的所有几何对象，因此在该多边形上调用区域方法，如下所示：
poly.area
给出该多边形的面积（以平方度为单位）。为了获得像平方米这样的平面单位的面积，我想我必须使用不同的投影（哪一个？）来转换多边形的坐标。
我多次读到 pyproj 库应该提供执行此操作的方法。使用pyproj，有没有办法将整个形状良好的Polygon对象转换为另一个投影，然后计算面积？
我用我的多边形做了一些其他的事情（不是你现在想的那样），只有在某些情况下，我需要计算面积。
到目前为止，我只找到了这个例子： http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/
这意味着将每个 Polygon 对象分割为其外环和内环（如果存在），抓取坐标，将每对坐标转换为另一个投影并重建 Polygon 对象，然后计算其面积（那么它是什么单位？） 。这看起来是一个解决方案，但不太实用。
还有更好的想法吗？
计算一个测地面积，非常准确，只需要一个椭球体（而不是投影）。这可以使用 pyproj 2.3.0 或更高版本来完成。
from pyproj import Geod
from shapely import wkt
poly = wkt.loads("""\
POLYGON ((-116.904 43.371, -116.823 43.389, -116.895 43.407,
-116.908 43.375, -116.904 43.371))
""")
print(f"Cartesian area: {poly.area:.6f}")
# Cartesian area: 0.001467
# Specify a named ellipsoid
geod = Geod(ellps="WGS84")
geod_area = abs(geod.geometry_area_perimeter(poly)[0])
print(f"Geodesic area: {geod_area:.3f} m^2")
# Geodesic area: 13205034.647 m^2
abs()好吧，我终于用matplotlib库的Basemap工具包做到了。我会解释它是如何工作的，也许这有时会对某人有所帮助。
1. 在您的系统上下载并安装 matplotlib 库。 http://matplotlib.org/downloads.html
对于 Windows 二进制文件，我建议使用此页面： http://www.lfd.uci.edu/~gohlke/pythonlibs/#matplotlib 请注意以下提示：
需要 numpy、dateutil、pytz、pyparsing、六个，以及可选的枕头， pycairo、龙卷风、wxpython、pyside、pyqt、ghostscript、miktex、ffmpeg、 mencoder、avconv 或 imagemagick。
因此，如果您的系统上尚未安装，则必须下载并安装 numpy、dateutil、pytz、pyparsing 和 6 个软件才能使 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)
转换为弧度并假设地球是半径为 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