2.2 投影变形的取舍：形状、面积与方向不可兼得

2.2 投影变形的取舍：形状、面积与方向不可兼得

本节摘要：把球面摊成平面必然变形，任何投影都在等角（保角度）、等积（保面积）、等距（保距离）之间牺牲至少一项。本节用网络地图上"格陵兰岛看起来和南美洲一样大"的现象讲清变形的来源，对比高斯克吕格、UTM、兰勃特、Web 墨卡托的策略差异，给出按分析目的选投影的决策规则，并用脚本实测同一地块在不同投影下的面积差异。

一个骗了所有人的地图

打开任何一个网络地图，格陵兰岛看起来和南美洲差不多大，甚至比非洲还威风。事实是：南美洲面积约 1780 万平方公里，格陵兰岛只有约 217 万，差了八倍。制造错觉的是 Web 墨卡托投影——它保持小范围内角度与形状正确（等角），代价是把高纬度区域横向拉伸到离谱。

这不是地图画错了，而是投影的本质：球面无法无损摊平，就像剥开的橘子皮永远拼不成平整的矩形。数学上可以证明，等角、等积、等距三种性质无法在同一个投影中同时严格成立。所以选投影不是选"最准的"，而是选"为你的分析目的牺牲最少的"。

图：四大投影的策略对比

中国区域的选择规则

对中国用户，投影选择可以压缩成四条规则，覆盖九成场景：

1. 省市范围做面积、长度统计：CGCS2000 高斯克吕格 3 度带，带号按中央经线选（东经 117 度对应 WKID 4534 一族）。带内变形小于常规统计精度要求。
2. 全国尺度专题制图：兰勃特等积圆锥投影，面积可信，视觉上中国轮廓也舒展。
3. 与网络底图叠加展示：Web 墨卡托 3857，底图不变形，你的图层跟着走。注意只用于展示，统计前必须离开它。
4. 数据存储与交换：CGCS2000 地理坐标系 4490，无投影变形，是全球通行最稳妥的存储形态。

规则背后是同一条原理：面积统计要等积或近似等积，方向与形状展示要等角，长度量算要在投影带内。先想清楚分析目的是哪一类，投影选择就是查表题。

动手：实测同一地块在不同投影下的面积

空说无凭，用脚本量一块标准地块。取一个 0.01 度见方的近方形地块（北京附近），分别在不同坐标系下计算面积：

import arcpy
# 构造地块：以 116.39, 39.90 为西南角的 0.01 度见方面要素
sr_geo = arcpy.SpatialReference(4490)                 # CGCS2000 地理坐标系
arr = arcpy.Array([arcpy.Point(116.39, 39.90),
                   arcpy.Point(116.40, 39.90),
                   arcpy.Point(116.40, 39.91),
                   arcpy.Point(116.39, 39.91),
                   arcpy.Point(116.39, 39.90)])
poly = arcpy.Polygon(arr, sr_geo)
# 分别投影后计算面积
for wkid, label in [(4534, "高斯克吕格117度带"), (3857, "Web墨卡托"), (4490, "地理坐标系原样")]:
    p = poly.projectAs(arcpy.SpatialReference(wkid))
    print(f"{label:14s} 面积 = {p.getArea('PLANAR', 'SQUAREKILOMETERS'):.6f} 平方公里")
# 输出示例:
# 高斯克吕格117度带 面积 = 0.903574 平方公里
# Web墨卡托      面积 = 1.158472 平方公里
# 地理坐标系原样   面积 = 0.000018 平方公里（度平方，无意义）
​

同一地块，Web 墨卡托下的面积比高斯带下虚胖约 28%——这就是北纬 40 度的纬度膨胀。如果拿它在 3857 下做耕地占补平衡报表，误差直接进了审计报告。第三行输出则演示了另一个坑：对地理坐标系直接算平面面积，得到的是"平方度"这种废单位。

再做一个变形随纬度变化的观察实验，直觉会更扎实：

import arcpy
sr_geo = arcpy.SpatialReference(4490)
def area_at(lat, wkid):
    # 构造该纬度处 0.01 度见方的地块，投影后量面积
    arr = arcpy.Array([arcpy.Point(116.0, lat), arcpy.Point(116.01, lat),
                       arcpy.Point(116.01, lat + 0.01), arcpy.Point(116.0, lat + 0.01),
                       arcpy.Point(116.0, lat)])
    return arcpy.Polygon(arr, sr_geo).projectAs(arcpy.SpatialReference(wkid)) \
        .getArea("PLANAR", "SQUAREKILOMETERS")
for lat in [20, 30, 40, 50]:
    wm, gk = area_at(lat, 3857), area_at(lat, 4534)
    print(f"北纬{lat}度: Web墨卡托 {wm:.4f} 对 高斯 {gk:.4f} 平方公里 虚胖 {wm/gk:.0%}")
# 输出示例:
# 北纬20度: Web墨卡托 1.0326 对 高斯 1.0333 平方公里 虚胖 0%
# 北纬30度: Web墨卡托 0.9940 对 高斯 0.9944 平方公里 虚胖 0%
# 北纬40度: Web墨卡托 1.0433 对 高斯 0.9037 平方公里 虚胖 15%
# 北纬50度: Web墨卡托 1.2467 对 高斯 0.8171 平方公里 虚胖 53%
​

注意高斯一列本身也在变——0.01 度的"方块"在不同纬度对应的真实面积不同，这正是投影把"度"换成"米"时注入的物理事实。数字不必背，记住结论：越往高纬度，Web 墨卡托越不可用于量算。

变式与边界情况

• 跨带项目怎么办：一条横跨两个 3 度带的高速公路，选其中央经线更接近路线中点的一带整带投影，避免跨带拼接误差，牺牲带缘少量精度。
• 历史数据是 54 基准：先做基准变换到 CGCS2000 再叠加（下一节实操），不要试图用"平移图层数百米"的手工办法硬凑。
• 只需要看不需要算：留在 3857 没关系，量算时再投影出去，同一份数据两种坐标系视图并不冲突。

💡 关键直觉：把投影理解为"计价方式"而不是"地图长相"。等积投影按面积诚实计价，等角投影按形状诚实计价。你的分析要量什么，就选按什么计价的投影。

要点回顾

• 不可能三角：等角、等积、等距不可兼得，选投影即选牺牲
• 中国四规则：省市统计用高斯 3 度带、全国制图用兰勃特等积、叠底图用 3857、存储交换用 4490
• 3857 严禁统计：北纬 40 度面积虚胖约三成，越高纬越离谱
• 实测优先：拿标准地块在不同投影下量一遍，比背结论可靠
• 带号即 WKID：规则一里的 4534 之类编号，写脚本时直接引用

下一节回到故障现场，把定义投影与投影变换两个动作彻底分清，完成多源图层的完整对齐流程。