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

wowwwai GIS研习社 · 工具流程与项目排障

引言：很多同学在 Python 里做矢量数据处理时，会遇到一个典型问题：Shapely计算面积不对?投影需要转吗? 简短回答是：如果你的几何坐标还是经纬度，也就是常见的 EPSG:4326，那么直接用 Shapely 的 .area 计算面积通常不符合真实平方米面积；在多数 GIS 面积统计场景下，需要先转到合适的投影坐标系，再计算面积。

背景：为什么 Shapely计算面积不对

Shapely 是 Python GIS 中常用的几何计算库，可以做缓冲区、相交、裁剪、面积、长度等操作。但 Shapely 本身并不知道你的坐标代表什么单位，它只会按照几何坐标的数值进行平面几何计算。

这就会导致一个常见误解：数据是 WGS84 经纬度，坐标看起来像 116.3, 39.9，然后直接写：

area = polygon.area

此时得到的不是平方米，也不是平方千米，而是“经纬度坐标数值构成的平面面积”，可以粗略理解为平方度。这个值通常不能直接用于土地面积、行政区面积、地块面积、缓冲区面积等正式统计。

如果你发现 Shapely 面积结果特别小、和 QGIS 或 ArcGIS Pro 里看到的面积差很多，优先检查坐标系单位，而不是先怀疑 Shapely 算错了。

原理：Shapely area 计算的到底是什么

Shapely 的 .area 属性计算的是平面笛卡尔坐标系中的几何面积。它不主动读取 CRS，也就是坐标参考系。CRS 是 GIS 里描述坐标含义的规则，例如 EPSG:4326 表示 WGS84 经纬度坐标系，单位是度；EPSG:3857 表示 Web Mercator 投影坐标系，单位通常是米。

关键点有三个：

• Shapely 只看坐标数值：它不知道 116 是经度、米还是其他单位。
• 面积单位来自坐标单位的平方：坐标单位是米，面积结果就是平方米；坐标单位是度，面积结果就是平方度。
• 经纬度不是等距平面坐标：同样 1 度经度在赤道和高纬地区对应的实际距离不同，不能直接当作米来算面积。

所以，Shapely计算面积不对的根本原因，通常不是算法错误，而是输入几何没有被转换到适合面积计算的投影坐标系。

步骤：用 GeoPandas 和 Shapely 正确计算面积

步骤 1：先确认数据坐标系

如果你的数据来自 Shapefile、GeoPackage、GeoJSON 或 PostGIS，建议先用 GeoPandas 读取，并检查 CRS：

import geopandas as gpd

gdf = gpd.read_file("parcels.shp")
print(gdf.crs)

如果输出类似下面这样，说明数据是经纬度坐标：

EPSG:4326

如果输出为空，说明数据没有声明坐标系。此时不要急着投影转换，应该先确认原始数据真实坐标系，再用 set_crs 设置。

gdf = gdf.set_crs("EPSG:4326")

注意：set_crs 只是声明坐标系，不会改变坐标数值；to_crs 才是真正进行投影转换。

步骤 2：选择合适的投影坐标系

要解决 Shapely计算面积不对，核心是把经纬度几何转换到合适的投影坐标系。常见选择如下：

• 小范围地块或城市级数据：优先使用当地 UTM 分带、高斯克吕格分带或地方坐标系。
• 中国区域分析：可根据项目规范使用 CGCS2000 高斯克吕格投影分带。
• 全球或跨区域面积统计：考虑等面积投影，例如 Albers Equal Area 或 Mollweide。
• 仅用于网页展示：EPSG:3857 可以显示方便，但不推荐作为严肃面积统计的默认选择。

例如，数据位于北京附近，可以根据实际项目选择合适的米制投影。这里仅用一个示例 EPSG 说明流程，实际项目请替换为你的目标投影：

gdf_proj = gdf.to_crs("EPSG:32650")

EPSG:32650 是 WGS84 UTM 50N，单位为米，适合部分位于该分带范围内的北半球区域。若数据不在该分带内，不应照抄。

步骤 3：计算平方米面积和平方千米面积

投影完成后，再调用 Shapely 的 .area，此时结果单位就是投影坐标单位的平方。若投影单位是米，面积就是平方米。

gdf_proj["area_m2"] = gdf_proj.geometry.area
gdf_proj["area_km2"] = gdf_proj["area_m2"] / 1_000_000

print(gdf_proj[["area_m2", "area_km2"]].head())

这一步得到的结果才适合用于常见的面积统计、报表汇总和专题制图。

步骤 4：如果只有 Shapely 几何对象，使用 pyproj 转换

有时你手里不是 GeoDataFrame，而是一个单独的 Shapely Polygon。这时可以结合 pyproj 和 shapely.ops.transform 做投影转换。

from shapely.geometry import Polygon
from shapely.ops import transform
from pyproj import Transformer

poly = Polygon([
    (116.30, 39.90),
    (116.31, 39.90),
    (116.31, 39.91),
    (116.30, 39.91),
    (116.30, 39.90)
])

transformer = Transformer.from_crs("EPSG:4326", "EPSG:32650", always_xy=True)
poly_proj = transform(transformer.transform, poly)

area_m2 = poly_proj.area
print(area_m2)

这里的 always_xy=True 很重要，它可以确保坐标顺序按 x, y 处理，也就是经度、纬度，避免轴顺序导致的投影错误。

步骤 5：把结果写回文件

如果你需要把面积字段保存下来，可以写出为 GeoPackage 或 Shapefile。更推荐 GeoPackage，因为字段名、编码和几何类型支持更稳定。

gdf_proj["area_m2"] = gdf_proj.geometry.area
gdf_proj.to_file("parcels_with_area.gpkg", layer="parcels", driver="GPKG")

如果必须保存为 Shapefile，字段名尽量控制在 10 个字符以内，避免被截断。

常见坑：Shapely投影转换后面积仍然不对

坑 1：把 set_crs 当成 to_crs

set_crs 只是在数据上贴一个 CRS 标签，不会改变坐标。很多面积错误来自把经纬度数据错误声明成米制投影。

# 错误示例：这不会把经纬度转换成米
gdf_wrong = gdf.set_crs("EPSG:32650", allow_override=True)

正确做法是：原始 CRS 正确声明后，再用 to_crs 转换。

gdf = gdf.set_crs("EPSG:4326")
gdf_proj = gdf.to_crs("EPSG:32650")

坑 2：直接用 EPSG:3857 计算正式面积

EPSG:3857 是 WebGIS 常用投影，适合地图瓦片显示，但面积变形明显。对于城市级粗略估算可能还能接受，但如果是地籍、生态红线、耕地保护、行政统计等正式场景，不建议直接用 EPSG:3857 面积作为最终结果。

坑 3：跨多个投影分带的数据用一个 UTM 分带

UTM 或高斯克吕格分带适合特定经度范围。如果数据横跨多个分带，用单一分带可能带来较大变形。跨省、全国或全球数据更适合选择等面积投影，或者分区投影后再汇总。

坑 4：几何无效导致面积异常

自相交、多部件异常、空几何等问题也会影响面积结果。计算前建议检查几何有效性：

invalid = gdf[~gdf.geometry.is_valid]
print(len(invalid))

如果存在无效几何，可根据数据情况修复。常见方式包括 Shapely 2.x 的 make_valid 或缓冲修复法：

from shapely.validation import make_valid

gdf["geometry"] = gdf.geometry.apply(make_valid)

坑 5：面积结果看起来差 10000 倍或 1000000 倍

这种情况通常是单位换算错误。平方米转公顷应除以 10000，平方米转平方千米应除以 1000000。

gdf_proj["area_ha"] = gdf_proj.geometry.area / 10000
gdf_proj["area_km2"] = gdf_proj.geometry.area / 1000000

方法比较：Shapely、GeoPandas、QGIS 怎么选

如果你是在 Python 中批量处理矢量文件，推荐使用 GeoPandas 管理 CRS，用 Shapely 执行几何计算。如果只是排查 Shapely计算面积不对，可以先在 QGIS 中加载同一份数据，对比图层 CRS、投影后面积字段和 Python 结果是否一致。

检查清单：计算前先确认这 8 件事

• 数据是否有 CRS？如果没有，是否知道真实坐标系？
• 当前 CRS 是否是 EPSG:4326 这类经纬度坐标系？
• 是否使用了 to_crs 进行真实投影转换，而不是只用 set_crs？
• 目标投影是否适合数据所在区域？
• 目标投影单位是否为米？
• 是否需要平方米、公顷、亩或平方千米之间的单位换算？
• 几何是否存在自相交、空几何或无效多边形？
• 结果是否与 QGIS、ArcGIS Pro 或已知样本面积做过抽查对比？

实务建议：不要只看代码能不能运行，要看坐标系、单位和业务尺度是否匹配。面积统计的错误，很多时候不是发生在 .area 这一行，而是发生在投影选择之前。

FAQ：Shapely计算面积不对的常见问题

Q1：Shapely 计算面积必须投影吗？

如果数据是经纬度坐标，通常需要先投影再计算面积。如果数据已经是合适的米制投影坐标系，可以直接用 Shapely 的 .area。判断依据不是文件格式，而是坐标单位和投影是否适合面积计算。

Q2：EPSG:4326 的 Shapely area 单位是什么？

EPSG:4326 的坐标单位是度，因此 Shapely 直接计算出来的面积可以理解为平方度，不是平方米。平方度不能简单乘一个固定系数转换成平方米，因为经纬度对应的实际距离会随纬度变化。

Q3：Shapely 能不能自动识别坐标系？

不能。Shapely 的几何对象只保存坐标和拓扑结构，不保存完整 CRS 信息。需要 GeoPandas、pyproj、Fiona、Rasterio 或 PostGIS 这类工具来管理和转换坐标系。

Q4：为什么 GeoPandas 的 area 和 Shapely 的 area 结果一样？

GeoPandas 的 gdf.geometry.area 底层也是对几何对象做平面面积计算。区别在于 GeoPandas 可以保存 CRS，并提供 to_crs 方法。只有先转换到合适投影后，面积结果才会变成有意义的平方米或其他投影单位平方。

Q5：用什么投影算面积最准确？

没有一个投影适合所有区域。小范围项目可选当地官方投影、UTM 分带或高斯克吕格分带；大范围统计优先考虑等面积投影。正式项目应遵循数据生产规范或行业标准，而不是随意选择 EPSG:3857。

Q6：Shapely 计算出来的面积能和 QGIS 完全一致吗？

在相同几何、相同 CRS、相同投影、相同单位换算条件下，结果通常应非常接近。如果不一致，优先检查 QGIS 图层 CRS、项目 CRS、字段计算表达式、椭球面积设置，以及 Python 中是否真的执行了 to_crs。

结论：先投影，再用 Shapely 算面积

遇到 Shapely计算面积不对时，最重要的排查思路是：先看坐标系，再看单位，最后再看几何本身。Shapely 的 .area 没有错，它只是按照当前坐标数值做平面面积计算。

如果你的数据是 EPSG:4326 经纬度，直接计算面积通常不适合业务使用。正确流程是：确认原始 CRS，选择适合研究区域的投影坐标系，用 GeoPandas 或 pyproj 完成投影转换，再用 Shapely 或 GeoPandas 计算面积，并完成平方米、公顷或平方千米换算。

一句话总结：Shapely 负责几何计算，投影和单位要由你负责。 只要把 CRS、投影和单位这三件事处理清楚，Python 中的面积统计就会稳定很多。