坐标与投影 | Knowledge

前情回顾：01 第一帧讲清了渲染循环的通用 5-6 步骨架与典型实现 8 步细化。其中”绘制”环节暗藏一个根本问题——引擎拿到的相机/对象位置到底是哪种坐标？为什么不能直接用经纬度喂给 GPU？本篇就打开这个黑盒。

直觉问题

打开任何一个 Web 地图，左下角都显示着类似 (116.40°E, 39.90°N) 的两个数字。这两个数字背后藏着两个让人困惑的事实：

同一对经纬度，在 3D 地球和 2D 平面地图上都指向”同一个地方”——但 3D 是球面，2D 是平面，球面和平面怎么可能用同一对坐标？
格陵兰岛在 Web 地图上看起来和非洲差不多大——但非洲实际面积是格陵兰的 14 倍。地图在骗你，骗你的方式还很有规律。

这两个问题的答案合起来就是 GIS 入门最大的认知坎：球面上的位置怎么变成屏幕上的像素。本篇要做四件事：

讲清”经纬度到底是什么”（大地测量学的 4 层抽象）
给出 WGS84 椭球的精确数学定义（不是完美球）
推导经纬度 → ECEF 笛卡尔坐标 → 投影米 → 屏幕像素这条变换链
解释 Web 墨卡托投影为什么在纬度 ±85.0511° 截断

读完这一篇，你能回答：“为什么 GPS 给的高度不能当海拔用？”、“为什么经纬度不能直接勾股定理算距离？”、“为什么 EPSG:4326 和 EPSG:3857 数值差 5 个数量级？“

核心概念白话讲

用”橘子皮”理解：为什么球面坐标必须用经纬度

想象把地球当一个橘子。剥橘子皮铺平在桌上——你发现做不到，橘子皮总会撕烂或皱起来。球面和平面在拓扑上不等价，这是高斯”绝妙定理”（Theorema Egregium）的推论。

所以描述球面上的位置，必须用”球面坐标系”：

经度 $\lambda$ ：东西方向位置（-180° 到 180°），本初子午线为 0°
纬度 $\varphi$ ：南北方向位置（-90° 到 90°），赤道为 0°
高度 $h$ ：离开椭球面的距离（米）

这三元组 $(\lambda, \varphi, h)$ 就是你看到的 GPS 给出的”经纬度”。

NOTE

关键认知：经纬度是球面坐标，不能直接做平面距离/面积计算。两个经纬度点的”距离”必须用球面公式（如 haversine），不是勾股定理。这是入门者最容易踩的第一个坑。

用”等高线”理解：地球形状的 4 层抽象

地球的真实形状是一个凹凸不平的土豆（山脉、海沟、平原）。直接拿土豆做数学计算太复杂，大地测量学家做了 4 层抽象：

层	名称	形状描述	用途
①	地形面 (Topography)	真实地表，含山脉海沟	高程数据（DEM）参考
②	大地水准面 (Geoid)	全球海平面延伸到大陆下方的等位面	海拔高的参考面（凹凸差 ±100m）
③	参考椭球面 (Ellipsoid)	数学旋转椭球（WGS84）	三维渲染与经纬度定义
④	球面近似 (Sphere)	完美球（半径 = 赤道半径）	Web 墨卡托等简化投影

Web GIS 主要用 ③ 参考椭球面（精确）和 ④ 球面近似（简化）。② 大地水准面只在涉及海拔数据时出现，本篇不展开（留 05 地形与 Worker 详讲）。

用”剥橘子”理解：投影必然失真

把球面变成平面叫投影。无论怎么投，都会失真——这是数学定理（没有完美投影）。三种主要失真方向：

角度失真：原来 90° 的角，投影后变了
面积失真：原来等大的两块区域，投影后不等了
距离失真：原来等距的两段，投影后不等了

这三者不可兼得。任何投影最多保一个，必然牺牲其他两个：

保角度 → Mercator 投影（面积失真，格陵兰看着和非洲一样大）
保面积 → Mollweide 投影（角度失真，形状扭曲）
保距离 → Equirectangular 投影（皆非完美，距离近似）

Web 地图普遍选 Mercator（保角）—— 因为航海/导航看方向，方向不能错；面积错了人眼不敏感。

为什么需要 EPSG 编码

世界上有上千种坐标系（不同椭球、不同投影、不同分带）。如果两份地图数据用不同坐标系，叠加就会偏移几公里。EPSG 编码（International Association of Oil & Gas Producers 维护）给每个坐标系分配唯一编号：

EPSG:4326 — WGS84 经纬度（单位度）
EPSG:3857 — Web Mercator（单位米）
EPSG:3395 — 标准墨卡托（椭球版）
EPSG:32650 — UTM Zone 50N（北京附近分带投影）

所有 GIS 数据必须附带 EPSG 编号，否则不可信。

原理与数学/机制

1. WGS84 椭球参数

WGS84（World Geodetic System 1984）是美国国防部定义、GPS 卫星使用的全球参考椭球。它的核心参数 5 个：

参数	符号	数值	含义
长半轴	$a$	$6\,378\,137.0$ m	赤道半径
短半轴	$b$	$6\,356\,752.314\,245$ m	极半径
扁率	$f$	$1/298.257\,223\,563$	$f = (a-b)/a$
第一偏心率平方	$e^2$	$6.694\,379\,990\,14 \times 10^{-3}$	$e^2 = 2f - f^2$
第二偏心率平方	$e'^2$	$6.739\,496\,742\,28 \times 10^{-3}$	$e'^2 = e^2/(1-e^2)$

关键观察：扁率 $f \approx 1/298.257$ ，意味着地球极半径比赤道半径短约 21 km。地球不是完美球，是”略扁的椭球”。这个 21 km 在三维渲染里就是相机的最大坐标精度瓶颈（推动 06 相机与拾取要用 RTE 高精度算法）。

2. 经纬度 → ECEF 变换公式

ECEF（Earth-Centered, Earth-Fixed）是以地球质心为原点的三维笛卡尔坐标系：

$x$ 轴：指向赤道与本初子午线交点
$z$ 轴：指向北极（地球自转轴）
$y$ 轴：右手系决定（指向赤道与东经 90° 交点）

为什么 GIS 引擎内部用 ECEF 而不是经纬度？因为线性代数（旋转、平移、投影矩阵）只在笛卡尔坐标系里成立。GPU 算的是矩阵乘法，必须先把经纬度变成 ECEF。

变换公式（从 $(\lambda, \varphi, h)$ 到 $(X, Y, Z)$ ）：

\begin{cases} N = \dfrac{a}{\sqrt{1 - e^2 \sin^2\varphi}} \\ X = (N + h)\cos\varphi\cos\lambda \\ Y = (N + h)\cos\varphi\sin\lambda \\ Z = \bigl(N(1 - e^2) + h\bigr)\sin\varphi \end{cases}

符号说明：

$\lambda$ — 经度（弧度，东经为正）
$\varphi$ — 纬度（弧度，北纬为正）
$h$ — 椭球面以上高度（米），不是海拔
$N$ — 卯酉圈曲率半径（Radius of curvature in the prime vertical），几何意义是”沿椭球面外法线方向到 z 轴的距离”

$N$ 的直观理解：在赤道， $N \approx a$ （椭球最鼓处，法线短）；在极点， $N \approx c = a/\sqrt{1-e^2} \approx 6\,399\,593$ m（椭球最扁处，法线长）。

3. Web 墨卡托投影（EPSG:3857）

Web 墨卡托是 Mercator 投影的球面近似版本——假设地球是完美球体（用赤道半径 $R = a$ ），简化计算。Google Maps / OpenStreetMap / Bing 都用它。

正变换（lonLat → 米）

\begin{cases} x = R\lambda \\ y = R\ln\tan\!\left(\dfrac{\pi}{4} + \dfrac{\varphi}{2}\right) \end{cases}

反变换（米 → lonLat）

\begin{cases} \lambda = \dfrac{x}{R} \\ \varphi = 2\arctan e^{y/R} - \dfrac{\pi}{2} \end{cases}

其中 $R = a = 6\,378\,137.0$ m。

±85.0511° 截断的数学根源

观察正变换的 $y$ 公式：当 $\varphi \to \pm 90°$ （极点）时， $\tan(\pi/4 + \varphi/2) \to \pm\infty$ ，所以 $y \to \pm\infty$ 。

极点在 Web 墨卡托上不存在——它会跑到无限远。实际工程上约定截断到让 $y_{max} = R\pi$ （与 $x_{max}$ 相等，保证正方形瓦片），反推截断纬度：

R\ln\tan\!\left(\dfrac{\pi}{4} + \dfrac{\varphi_{max}}{2}\right) = R\pi

\tan\!\left(\dfrac{\pi}{4} + \dfrac{\varphi_{max}}{2}\right) = e^\pi

\varphi_{max} = 2\arctan(e^\pi) - \dfrac{\pi}{2} \approx 85.051\,128\,78°

这就是为什么 Web 地图都限制在纬度 ±85.0511°。极地的部分用别的投影（如 UPS Polar Stereographic）补充。

4. haversine 公式：球面两点距离

经纬度不能直接勾股定理。两点 $(\lambda_1, \varphi_1)$ 与 $(\lambda_2, \varphi_2)$ 的球面距离用 haversine 公式：

d = 2R\arcsin\!\left(\sqrt{\sin^2\!\left(\dfrac{\Delta\varphi}{2}\right) + \cos\varphi_1\cos\varphi_2\sin^2\!\left(\dfrac{\Delta\lambda}{2}\right)}\right)

其中 $\Delta\varphi = \varphi_2 - \varphi_1$ ， $\Delta\lambda = \lambda_2 - \lambda_1$ ， $R = 6\,371\,000$ m（地球平均半径，比赤道半径略小）。

例子：北京 (116.40°E, 39.90°N) 到上海 (121.47°E, 31.23°N)：

$\Delta\varphi = -8.67° \approx -0.151\,3$ rad
$\Delta\lambda = 5.07° \approx 0.088\,5$ rad
代入 haversine： $d \approx 1\,067$ km ✓

5. 投影变形机制（Tissot’s indicatrix）

Tissot’s indicatrix（底索变形椭圆）是 19 世纪法国数学家 Tissot 发明的可视化工具：在球面上画一组等大的小圆，投影后看它们变成什么形状。

ASCII 示意：Web 墨卡托下的 Tissot 圆

   原始球面：纬度均匀分布的等大圆         Web 墨卡托投影后：圆变椭圆

   纬度 80°  ●  ●  ●  ●  ●              ◯ ▭ ◯ ▭ ◯ ▭ ◯ ▭ ◯    ← 极端纵向拉长
   纬度 60°  ●  ●  ●  ●  ●              ◯  ◯  ◯  ◯  ◯         ← 高纬拉长
   纬度 40°  ●  ●  ●  ●  ●              ◯  ◯  ◯  ◯  ◯         ← 中纬轻微
   纬度 20°  ●  ●  ●  ●  ●              ●  ●  ●  ●  ●         ← 接近圆形
   纬度  0°  ●  ●  ●  ●  ●              ●  ●  ●  ●  ●         ← 赤道完全圆形
   纬度-20°  ●  ●  ●  ●  ●              ●  ●  ●  ●  ●
   纬度-40°  ●  ●  ●  ●  ●              ◯  ◯  ◯  ◯  ◯
   纬度-60°  ●  ●  ●  ●  ●              ◯  ◯  ◯  ◯  ◯
   纬度-80°  ●  ●  ●  ●  ●              ◯ ▭ ◯ ▭ ◯ ▭ ◯ ▭ ◯

   圆都一样大                            越往两极椭圆越细长（无限趋向极点）

关键观察：

赤道附近圆保持圆形（投影是各向同性放大）
越往两极椭圆越在纵向拉长（这是”格陵兰看起来和非洲一样大”的根源——非洲在赤道带，格陵兰在北纬 60-80°）
角度（圆的形状原本对称）保持，所以是保角投影

三类投影对比

投影类型	变形方向	Tissot 椭圆形态	典型代表	适用场景
保角 (Conformal)	角度正确，面积失真	椭圆但局部形状保持	Mercator / Lambert Conformal	航海、导航（方向不能错）
等积 (Equal-area)	面积正确，角度失真	椭圆但面积守恒	Mollweide / Albers	统计、人口密度（面积不能错）
等距 (Equidistant)	某方向距离正确	椭圆各方向有差异	Equirectangular / Azimuthal	简化示意（教学、概览）

6. 变换链：从经纬度到屏幕像素

图 1：从经纬度到屏幕像素的完整变换链（虚线是 2D 地图的瓦片旁路）

两条路径并存：

3D 地球渲染：走 ① → ② → ③ → ④ → ⑤（完整 ECEF + 视图投影管线）
2D 地图浏览：走 ① → ⑥ → ⑦（直接 Web 墨卡托投影米 → 瓦片编号，03 篇详讲）

IMPORTANT

关键不变量：无论走哪条路径，原始输入永远是经纬度 $(\lambda, \varphi, h)$ 。任何中间坐标（ECEF、投影米、瓦片编号、屏幕像素）都可逆推回到经纬度——这是 GIS 引擎能在 3D/2D 视图间无缝切换的数学基础。

可视化对比与动手实验

对比表一：常见 EPSG 坐标系

EPSG	名称	单位	椭球	投影	失真	Web 地图
4326	WGS84 经纬度	度	WGS84	无（球面坐标）	—	OGC 标准交换格式
3857	Web Mercator	米	球面近似	Mercator	高纬面积放大	Google / OSM / Bing
3395	标准墨卡托	米	WGS84 椭球	Mercator	高纬面积放大（略小）	专业航海
32650	UTM Zone 50N	米	WGS84 椭球	横轴墨卡托（分带）	分带内极小	北京区域工程测量

取北京 (116.40°E, 39.90°N) 的数值对比：

EPSG:4326 → (116.40, 39.90)
EPSG:3857 → (12_958_039.79, 4_849_678.58)（米）
EPSG:32650 → (447,351.6, 4_419_379.7)（米，分带内坐标）

为什么差这么多：4326 是角度（球面），3857 是 Web 墨卡托投影后的米（平面），单位完全不同，差 5 个数量级。两者绝不能混用。

对比表二：4 种投影类型（按几何面分类）

类型	投影面	失真特点	典型投影	适用场景
圆柱投影 (Cylindrical)	圆柱包球	赤道附近准确，两极失真大	Mercator / Equirectangular	全球地图、热带区域
圆锥投影 (Conic)	圆锥包球	中纬度准确	Albers / Lambert Conformal	中纬度国家地图（美、中、欧）
方位投影 (Azimuthal)	平面切球	中心准确，边缘失真	Orthographic / Stereographic	极地、半球视图
椭圆/多圆锥	数学曲面	自定义平衡	Robinson / Winkel Tripel	世界地图可视化（无主导失真）

对比表三：保角 / 等积 / 等距

类型	保什么	失什么	典型	用途
保角	角度（局部形状）	面积（高纬放大）	Mercator	航海、Web 地图
等积	面积	角度（形状扭曲）	Mollweide	统计、密度图
等距	某方向距离	两者都不完美	Equirectangular	教学示意

TIP

选投影的速记口诀：要方向用 Mercator（保角），要面积用 Mollweide（等积），要简单用 Equirectangular（等距），要美观用 Robinson（妥协）。

动手实验一：lonLatToCartesian 伪代码

按上面的公式，自己用任意语言实现经纬度 → ECEF：

输入: (λ, φ, h)  -- 经度(度)、纬度(度)、椭球高(米)
输出: (X, Y, Z)  -- ECEF 笛卡尔(米)

常量:
  a = 6378137.0
  e² = 6.69437999014e-3

步骤:
  1. 把 λ, φ 从度转弧度
  2. N = a / sqrt(1 - e² · sin²(φ))
  3. X = (N + h) · cos(φ) · cos(λ)
  4. Y = (N + h) · cos(φ) · sin(λ)
  5. Z = (N · (1 - e²) + h) · sin(φ)

验证：北京 (116.40°E, 39.90°N, h=44m) 应得到 ECEF 约 (-2_175_059, 4_389_724, 4_080_472)（米）。注意 $X$ 是负的（北京在东半球本应 $X$ 正？不对—— $X$ 轴指向本初子午线，东经 116° 在 $X$ 负、 $Y$ 正方向）。

动手实验二：epsg.io 查询

打开 epsg.io，搜索北京坐标 (116.40, 39.90)：

在右上角输入坐标 + 选 EPSG:4326
切换到 EPSG:3857，看坐标值变化
切换到 EPSG:32650（UTM 50N），看分带后坐标
注意 Web Mercator 的 $y \approx 4.85 \times 10^6$ m，远大于 UTM 的 $y \approx 4.42 \times 10^6$ m

动手实验三：Tissot 互动 demo

打开 Jason Davies - Tissot’s Indicatrix（如可访问）：

选 Mercator 投影，观察高纬椭圆纵向拉长
切换 Mollweide，观察椭圆变扁但面积守恒
切换 Equirectangular，观察椭圆各方向都变形

常见误区

WARNING

误区 1：把 GPS 给的 height 当海拔用。 ✅ 正确：height 是椭球面以上的高度，海拔是大地水准面以上的高度。两者差值最大 ±100 m（看 EGM2008 模型）。详见 05 地形与 Worker。

WARNING

误区 2：经纬度直接当屏幕像素用。 ✅ 正确：屏幕是 2D 投影后空间。必须先经纬度 → 投影米（或 ECEF）→ 视图矩阵 → 屏幕像素。

WARNING

误区 3：用 EPSG:4326 经纬度算两点距离直接勾股。 ✅ 正确：经纬度是球面坐标，必须用 haversine 公式（或更精确的 Vincenty 公式）。北京 (116, 40) 与 (117, 40) 不是”经度差 1° = 111 km”——只有赤道附近才近似成立。

WARNING

误区 4：把 Web 墨卡托（EPSG:3857）等同于标准墨卡托（EPSG:3395）。 ✅ 正确：前者用球面近似（假设地球完美球体），后者用 WGS84 椭球。两者数值会差几米到几十米，工程上不能互换。

WARNING

误区 5：以为 WGS84 是唯一参考椭球。 ✅ 正确：还有 GRS80、CGCS2000（中国国家标准）、Krassovsky（俄罗斯）等。它们与 WGS84 的差异 < 1 m，Web GIS 通常统一用 WGS84。

WARNING

误区 6：以为 Mercator 投影是”错的”（因为格陵兰看着和非洲一样大）。 ✅ 正确：Mercator 是保角投影，方向正确——这对航海导航至关重要。面积失真是保角的代价，不是设计失误。

延伸阅读与自测

权威参考

EPSG Registry — 国际坐标参考系统权威注册库
epsg.io — 最方便的 EPSG 在线查询工具
Wikipedia — Map projection — 投影学综述
Wikipedia — Web Mercator projection — Web 墨卡托专题
Wikipedia — Geodesy — 大地测量学综述
Wikipedia — Tissot’s indicatrix — 投影变形可视化原理
Wikipedia — Haversine formula — 球面距离公式
Map Projections - A Working Manual (Snyder, USGS) — 投影学经典教材（免费 PDF）

自测题

ECEF 手算：给定北京天安门 (116.397°E, 39.908°N, h=44m)，手算 ECEF 坐标。需要哪些 WGS84 参数？结果应该是 $X \approx -2.175 \times 10^6$ m， $Y \approx 4.389 \times 10^6$ m， $Z \approx 4.081 \times 10^6$ m。
截断推导：为什么 Web 墨卡托在纬度 ±85.0511° 截断？从 $y = R\ln\tan(\pi/4 + \varphi/2)$ 出发，让 $y_{max} = R\pi$ 反推 $\varphi_{max}$ 。
球面距离：用 haversine 公式算北京到上海的球面距离（结果应约 1067 km）。如果用墨卡托投影后勾股定理算，误差是多少？为什么？
投影选择：要做”全球人口密度可视化”，会选保角还是等积投影？为什么？做”导航 App 方向指示”呢？
EPSG 数值差异：北京同一点在 EPSG:4326 是 (116.40, 39.90)，在 EPSG:3857 是 (12_958_039, 4_849_678)。差 5 个数量级，为什么？两者各自适合什么用途？

下一篇导引：03 四叉树与 LOD 将打开瓦片系统的黑盒——把本篇讲的”投影后的米坐标”切成 z/x/y 编号的瓦片，按相机距离动态调度。Web 地图的”无限缩放”和”流畅响应”，都靠四叉树 LOD 算法支撑。