Rust 几何基础算法实现：从 geo crate 到 geo-toolbox

项目地址：https://github.com/Miku196/geo-toolbox

最近在看 geo-toolbox 这个项目，它是个用 Rust 写的地理空间工具箱，45 个包、15 个核心 crate 的规模，在 Rust GIS 圈子里不算小。我比较好奇它底层的几何算法是怎么组织的，翻了一下 core/geo-core 的源码，记录一些发现。

先说一句：下面涉及投影和距离选择的部分，有些是我一开始理解错了的，文中会标出来，避免误导。

geo-core 干了什么

先看它的 core/geo-core/src/lib.rs，模块声明就摆在那里：

//! geo-core: Shared types, geometry operations, and CRS registry.
//!
//! This crate is the foundation of geo-toolbox. All other crates
//! depend on it. It provides:
//!
//! - Unified error types ([`GeoError`])
//! - CRS registry with built-in common coordinate systems ([`crs::CrsRegistry`])
//! - Geometry type aliases and validation ([`types`])
pub mod config;
pub mod crs;
pub mod errors;
pub mod guard;
pub mod health;
pub mod observability;
pub mod plugin;
pub mod traits;
pub mod types;

关键点是：它没有自己的算法模块。几何基元直接从 geo-types 重导出（core/geo-core/src/types.rs 第 6 行就是 pub use geo_types;），空间谓词、布尔运算、距离计算、凸包、简化——这些全部交给 geo crate 处理。geo-core 自己写的，是类型系统、校验和 CRS 注册表。

这个选择挺务实。Rust 生态里 geo crate 已经很成熟了，DE-9IM 谓词、i_overlay 驱动的布尔运算、多种距离度量、QuickHull 凸包、RDP 简化，都有现成实现。geo-toolbox 的做法是在外面包一层业务约束，而不是重新造轮子。

空间谓词与 BBox 粗筛

geo crate 的空间谓词在概念上遵循 DE-9IM（九交矩阵）模型，contains、crosses、touches、within、intersects 等都基于它定义。实际实现里不一定每次都完整算九交矩阵——有些谓词会有专门的快速路径，但语义上以 DE-9IM 为准。

对于多边形，contains 遵循 OGC 规范，边界上的点不算包含，只有 intersects 才返回 true。这个坑我踩过。

geo-toolbox 没有重新实现这些谓词，而是在 core/geo-core/src/types.rs 里自己写了个 BBox 结构做快速相交检测（第 29–31 行）：

/// Check intersection with another BBox.
pub fn intersects(&self, other: &BBox) -> bool {
    self.min_x <= other.max_x
        && self.max_x >= other.min_x
        && self.min_y <= other.max_y
        && self.max_y >= other.min_y
}

就是最简单的 AABB 检测，O(1)。还有一个 union 方法做包围盒合并（第 32–34 行）：

/// Compute the union of two BBoxes.
pub fn union(&self, other: &BBox) -> BBox {
    BBox {
        min_x: self.min_x.min(other.min_x),
        min_y: self.min_y.min(other.min_y),
        max_x: self.max_x.max(other.max_x),
        max_y: self.max_y.max(other.max_y),
    }
}

思路很清楚：批量处理 GeoJSON 的时候，先用 BBox 过滤掉明显不相交的要素，再让 geo crate 做精确谓词判断。精确谓词的开销不小，粗筛能省很多时间。

不过这里有个适用范围的问题。BBox 相交检测只适合 intersects 这类谓词的粗筛。如果是 contains，粗筛逻辑要分情况看：

• 判断“点是否可能被几何体包含”，可以用 BBox::contains(x, y)（core/geo-core/src/types.rs 第 24 行），它判断的是点是否落在 BBox 范围内；
• 判断“几何体 B 是否可能被几何体 A 包含”，源码里没有现成的 BBox 对 BBox 包含方法，需要自己写 a.min_x <= b.min_x && a.max_x >= b.max_x && a.min_y <= b.min_y && a.max_y >= b.max_y 这类比较。

而且要注意：BBox 包含只是 contains 的必要条件，不是充分条件——B 的 BBox 落在 A 的 BBox 里，不代表 B 就真被 A 包含。粗筛之后仍然要走精确判断。

坐标校验：看起来简单但很有必要

core/geo-core/src/types.rs 里有个函数我一开始觉得没什么，后来想想挺关键：

/// Validate a single (lon, lat) coordinate pair.
///
/// Returns `Ok(())` if lon ∈ [-180, 180] and lat ∈ [-90, 90].
pub fn validate_coord(lon: f64, lat: f64) -> Result<(), crate::GeoError> {
    if !(-180.0..=180.0).contains(&lon) || !(-90.0..=90.0).contains(&lat) {
        Err(crate::GeoError::Validation(format!(
            "coordinate out of range: lon={lon}, lat={lat}"
        )))
    } else {
        Ok(())
    }
}

geo crate 本身不管你传什么坐标进去，都照单全收。但在实际 GIS 流水线里，经纬度超出范围基本意味着数据有问题——要么是 CRS 搞混了，把投影坐标当经纬度塞进来了，要么是数据本身损坏。与其让错误坐标一路流到后面的计算里产生莫名其妙的异常结果，不如在入口直接拒掉。

布尔运算：i_overlay 的整数量化

geo 的布尔运算现在底层用的是 i_overlay，不是早期的 Martinez-Rueda。切换的原因主要是浮点坐标在退化情况（比如边重合、顶点交叠）下容易崩，i_overlay 用整数坐标来提升鲁棒性。

具体做法是：浮点坐标先量化到整数网格，在整数域内做精确运算，算完再反量化回去。这里要区分两个概念——整数域内是精确的，但浮点转整数的量化本身有损，可能把邻近顶点合并、改变拓扑。所以它的定位是“鲁棒”，不是“无损精确”。做高精度场景时要留意这个区别。

另外 geo 提供了 unary_union 用来批量合并大量重叠多边形，比逐个调 union 快很多。

需要说明的是，上面这些是 geo crate 本身的能力，不是 geo-toolbox 自己实现的。geo-toolbox 的插件层（plugins/ 目录）在碳汇聚合这类场景中会调用到这些函数，但 geo-core 本身只是重导出类型，没有对布尔运算做额外封装。

距离计算与 CRS 的关系

这一节是我一开始理解错的地方，重新写一下。

core/geo-core/src/crs.rs 里有个枚举，把 CRS 按用途分了类（第 8–12 行）：

pub enum CrsCategory {
    /// Default storage CRS: EPSG:4326 (WGS84 lat/lon).
    Storage,
    /// Web map display: EPSG:3857 (Web Mercator).
    Display,
    /// Area-sensitive computations: EPSG:3405 (World Equal Area) or local UTM.
    Carbon,
    /// Local engineering / CAD coordinate system.
    CadLocal,
}

源码里这个枚举上面还有 #[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]，我上面为了简洁省略了。

顺便说一句，Carbon 这个注释本身写得也不够严谨——把等面积投影（EPSG:3405）和等角投影（UTM）并列放在一个类别里，容易让人误以为两者是一类东西。下面展开说。

先说投影类型。常见的几类：

• 等角投影：保持局部角度和形状，UTM、Web Mercator（EPSG:3857）都属于这类。面积有变形，且变形随纬度增大。
• 等面积投影：保持面积比例，EPSG:3405（World Equal Area）属于这类。角度和距离都有变形。
• 等距投影：保持某个方向或某个点的距离，其他方向会变形。

Carbon 类别里的 EPSG:3405 是等面积投影，用于面积计算（比如碳储量按面积折算），这个没问题。但本地 UTM 并不是等面积投影，它是横轴墨卡托（等角投影）。UTM 在小范围内面积变形很小，可以近似用于面积计算，但严格来说它不是等面积投影。之前我看到注释里把 UTM 和等面积并列，就想当然地写成“都是等面积投影”，这是错的。

再说距离。等面积投影只保证面积有意义，距离通常没有物理意义——在等面积投影上算平面欧几里得距离是错的方向。距离计算应该：

• 用测地线公式（GeodesicDistance 基于 Karney 2013 椭球算法，HaversineDistance 是球面近似），直接在地理坐标上算；或者
• 转换到局部投影（比如小范围 UTM，配合中央经线的比例因子 0.9996 修正）再用平面距离近似。

UTM 之所以能用于距离近似，不是因为它“等距”，而是因为小范围内变形小、比例因子可以修正。归类上它仍然是等角投影。

所以之前那句“输入是 4326 就用 Haversine，或者先转到投影坐标系再用欧几里得”要加个限定：转到投影坐标系之后能不能用欧几里得，取决于投影是不是等角/等距、范围多大、比例因子怎么修正。等面积投影不适合这个用途。

面积和距离是两个不同的需求，对应不同的投影选择——面积敏感的场景用等面积投影算面积；距离计算则单独用测地线或局部投影近似，不要指望等面积投影能同时解决距离问题。

geo crate 提供的距离 trait 有好几种：

• EuclideanDistance：平面勾股定理，只在合适的投影坐标系或小范围内有意义
• HaversineDistance：球面大圆距离，平均地球半径近似
• GeodesicDistance：Karney (2013) 椭球算法，精度最高
• VincentyDistance：经典椭球公式，精度高但某些情况（近对跖点）不收敛

geo-core 自己没有实现这些，它是通过 transform_point 在坐标变换层面保证单位一致。调用方根据 CrsCategory 自己选距离算法。

凸包与简化

geo 的 ConvexHull trait 实现了 QuickHull 算法，凸包顶点按逆时针排列。不同版本可能提供不同的后端选择（比如 qhull 和 graham），但这不是稳定保证，具体看版本和 feature 配置。

Simplify trait 是 RDP（Ramer–Douglas–Peucker）算法，还有个 SimplifyIdx 变体返回保留顶点的索引，方便回溯原始数据。

geo-toolbox 的野外 PWA 场景把这两个组合起来用：用户离线采集一批 GPS 点，先凸包生成区域范围，再用 RDP 简化边界。主要目的是降低存储量——IndexedDB 存东西也是要空间的。

看完源码的几点想法

翻完 core/geo-core，有几个地方印象比较深：

不重造轮子，但也不是无脑包一层。geo-core 没有试图“改进” geo 的算法，它的价值在于类型校验和 CRS 分类这些业务约束。这样 geo 自身升级（比如从 Martinez-Rueda 换到 i_overlay）不会波及上层代码。

校验前置。validate_coord 和 BBox 快速相交检测构成入口的第一道防线。GIS 流水线里脏数据越早拒绝越好，等到复杂计算中途才暴露问题，追溯起来很痛苦。

CRS 分类把隐式约定显式化了。CrsCategory 枚举把“该用哪种投影、该用哪种距离算法”这个决策变成类型系统的一部分，避免了业务代码里到处散落 if crs == 4326 这种判断。不过要记住的是：面积和距离是两个不同的需求，对应不同的投影选择——等面积投影适合算面积但不适合算距离，UTM 是等角投影不是等面积投影。

geo-toolbox 的 15 个核心 crate 里，geo-core 从模块数量看职责边界是最清晰的——9 个源文件，覆盖配置、CRS、错误、守卫、健康检查、可观测性、插件、trait 和类型。它定义了整个工作区的几何数据契约：哪些类型可用、什么坐标合法、什么 CRS 适合什么计算。上面的碳核算、水文分析、遥感辐射校正插件都基于这套语义构建。我觉得这个分层方式挺值得借鉴的。