Shapely计算面积不对?投影需要转吗?

“我算的面积怎么差了十倍？”——Shapely 的投影陷阱

你是不是也遇到过这种情况：用 Shapely 算出来的多边形面积，跟 ArcGIS 或 QGIS 里显示的结果完全对不上？甚至差出几个数量级？别慌，这不是代码写错了，而是你忘了最关键一步——投影转换。

我在参与某省耕地保护项目时，就曾因为没转投影，导致计算的地块面积比实际小了整整87%。领导拿着报告问我：“Dr. Gis，咱们省的耕地缩水了吗？”——那一刻，我真想找个地缝钻进去。

为什么经纬度不能直接算面积？

想象一下，你手里拿着一个橘子（地球），上面画了个圈。如果你直接在橘子皮上量这个圈的“长×宽”，得到的是弧面距离。但当你把橘子皮剥下来、压平到桌面上（投影平面），这个圈的形状和尺寸都变了——这才是我们能在地图上测量的真实面积。

Shapely 默认在笛卡尔平面坐标系下工作。如果你传入的是 WGS84 经纬度（EPSG:4326），它会天真地把经度当“X”、纬度当“Y”，然后套用平面几何公式计算面积——结果当然是错的，而且错得离谱。

实战：三步搞定正确面积计算

解决方法其实很简单：先转投影，再算面积。以下是标准操作流程：

1. 判断原始数据坐标系：确认你的 GeoJSON 或 Shapefile 是不是 WGS84（大多数在线数据默认都是）。
2. 选择合适的投影坐标系：根据区域位置选 UTM 分区，或 Albers 等面积投影。中国全域推荐 EPSG:4547（CGCS2000 / Gauss-Kruger CM 105E）或 EPSG:9807（Albers China）。
3. 用 PyProj + Shapely 联合处理：先投影变换，再调用 .area 方法。

from shapely.geometry import Polygon
from pyproj import Transformer

# 假设你有一个WGS84下的多边形
gps_poly = Polygon([(116.3, 39.9), (116.4, 39.9), (116.4, 40.0), (116.3, 40.0)])

# 创建投影转换器：从WGS84转到UTM Zone 50N（北京地区）
transformer = Transformer.from_crs("EPSG:4326", "EPSG:32650", always_xy=True)

# 对多边形每个点进行投影转换
projected_poly = transform(transformer.transform, gps_poly)

# 现在计算的面积单位是平方米！
area_sq_m = projected_poly.area
print(f"正确面积：{area_sq_m:.2f} 平方米")

常见误区与避坑指南

总结：投影不是可选项，是必选项

记住这个口诀：“经纬度不算面积，投影之后再算”。Shapely 本身没有错，错的是我们忘了给数据“穿对鞋”。地理计算中，坐标系就是数据的“语境”，脱离语境谈数值，必然南辕北辙。

你在项目中踩过类似的坑吗？或者有更好的投影选择经验？欢迎在评论区分享你的“血泪史”——让我们一起少走弯路，多产准数！