第09章 - 叠加分析详解
叠加分析(overlay analysis)是 GIS 里出现频率最高的一类几何运算:把两层矢量数据按空间位置叠在一起,求交集、并集、差集,从而回答“两块地重合的部分在哪”“这条管线穿过哪些行政区”“这个园区去掉水域之后剩多少面积”。在 JTS 中,这一族操作由 Geometry.intersection、Geometry.union、Geometry.difference、Geometry.symDifference 四个方法对外暴露,内部则牵出两套完全不同的算法实现:2016 年前后的旧引擎 OverlayOp,以及 1.18 引入、持续演进的 OverlayNG。
本章的目标有三个层次。第一层是会用:搞清四个方法的语义、返回类型、空几何与几何集合的处理规则,知道什么调用会抛异常、什么调用会静默给出错误结果。第二层是懂原理:拆开 OverlayNG 的流水线,看 noding、成边、建图、标注、构面这五步各自解决什么问题,理解 OverlayLabel 为什么要把边分成边界边、塌缩边、线边、非本几何边四种状态。第三层是能做工程决策:在浮点精度与固定精度之间取舍,知道 OverlayNGRobust 的重试阶梯在什么情况下救得了你,知道 CoverageUnion 什么时候能把并集从秒级压到毫秒级、又会在什么输入上悄悄算错。
阅读本章需要的前置知识:第 04 章的几何对象模型(Geometry、Polygon、MultiPolygon、GeometryFactory)、第 05 章的几何构造与 PrecisionModel、第 08 章的空间关系谓词与 DE-9IM。如果你只想先把功能跑起来,可以重点看 9.2、9.9、9.10 三节;如果你正被 TopologyException 折磨,直接跳到 9.7 与 9.11。
版本说明:本章以 JTS 1.21.0-SNAPSHOT 源码为准,示例代码在 1.18 以上均可编译。文中标注“1.18 新增”“1.20 新增”的特性,在更早版本中不存在。
9.1 叠加分析的数学语义
9.1.1 从集合论看四种运算
JTS 把几何对象看作平面上的点集(point set)。一个 Polygon 是它闭包内的所有点,一条 LineString 是其上所有点构成的连续点集,一个 Point 是单点集。在这个视角下,四个叠加运算就是集合的四种基本运算:
- 交(intersection):属于 A 且属于 B 的点,记作 A ∩ B。
- 并(union):属于 A 或属于 B 的点,记作 A ∪ B。
- 差(difference):属于 A 但不属于 B 的点,记作 A − B。注意它不对称,A − B 与 B − A 是两个不同结果。
- 对称差(symDifference):属于 A 或属于 B,但不同时属于两者的点,记作 A △ B,等价于 (A − B) ∪ (B − A)。
下面这张图用两个边长为 10 的正方形说明四种运算的几何形态。A 的坐标范围是 [0,10]×[0,10],B 是 [5,15]×[5,15],两者重叠一个 5×5 的方块:
┌───────────────┐
│ A │
│ ┌────────┼────────┐
│ │ A ∩ B │ B │
│ │ 25 │ │
└──────┼────────┘ │
│ │
└─────────────────┘
A = 100 B = 100
A ∩ B = 25 A ∪ B = 175
A − B = 75 B − A = 75
A △ B = 150(两个互不相连的 75 方块)
| 这里有个容易忽略的细节:面积满足容斥关系 | A ∪ B | = | A | + | B | − | A ∩ B | ,即 100 + 100 − 25 = 175;但对称差不是简单的 | A | + | B | − 2 | A ∩ B | 之外还要考虑边界重复计算——对区域几何而言这个式子是成立的(150),对线状几何则要小心重合线段被“去重”的影响。 |
9.1.2 四种运算的语义对照
下表把数学定义、JTS 方法、结果维度语义与常见用途放在一起,方便对照查阅:
| 运算 | 数学定义 | JTS 方法 | 结果维度 | 典型用途 |
|---|---|---|---|---|
| 交 | A ∩ B | intersection |
≤ min(dimA, dimB) | 提取重叠地块、裁剪、面与线求穿行段 |
| 并 | A ∪ B | union(other) |
≤ max(dimA, dimB) | 合并相邻宗地、溶解行政边界 |
| 差 | A − B | difference |
≤ dimA | 挖洞、扣除保护区、计算可用地 |
| 对称差 | A △ B | symDifference |
≤ max(dimA, dimB) | 变更检测、两期数据差异 |
注意“结果维度”一列用的是“小于等于”。原因在于结果可能发生降维:两个多边形只在一条边上相接,其交集是一条 LineString;两个多边形只在一个角点相接,其交集是一个 Point。这就是所谓拓扑塌缩(topology collapse),是叠加分析里大量“意外结果”的根源。
9.1.3 代数不变量与它们为何会失效
在理想数学世界里,四个运算满足一组漂亮的恒等式:
- 交换律:A ∩ B = B ∩ A,A ∪ B = B ∪ A,A △ B = B △ A。
- 结合律:(A ∪ B) ∪ C = A ∪ (B ∪ C),(A ∩ B) ∩ C = A ∩ (B ∩ C)。
- 分配律:A ∩ (B ∪ C) = (A ∩ B) ∪ (A ∩ C)。
- 对偶:A − (B ∪ C) = (A − B) ∩ (A − C);A − (B ∩ C) = (A − B) ∪ (A − C)。
但在浮点坐标上,这些恒等式只是“近似成立”。原因是每做一次叠加,输出顶点都会被重新计算——要么保留原始坐标(浮点精度),要么四舍五入到网格(固定精度)。一旦顶点坐标变化,下一次运算看到的就是另一个几何,于是出现一些工程上很恼人的现象:
- (A ∪ B) − B 不一定等于 A − B。并集可能把 B 的某个顶点吸附到 A 的边上,使得后续求差时多切掉一小块。
- 交换律在数值层面成立,在坐标输出顺序上不成立。
A.union(B)与B.union(A)的顶点集合通常相同,但环的起点、方向、以及MultiPolygon内各元素的顺序可能不同。做数据比对时不能直接比较 WKT 字符串,应先norm()再比较,或比较面积、周长、顶点集合等不变量。 - 面积不严格可加。对浮点结果,
(A ∪ B).getArea()与A.getArea() + B.getArea() - A.intersection(B).getArea()只在小数点后若干位上相等。
因此,工程上评价叠加结果时,应优先检查下面这组更稳健的不变量,而不是逐字符比较几何:
| 不变量 | 含义 | 检查方式 |
|---|---|---|
| 有效性 | 结果必须是合法几何 | result.isValid(),必要时 IsValidOp 定位问题 |
| 维数 | 结果维度不超过语义上界 | result.getDimension() |
| 面积/长度守恒 | 并、交、差之间的容斥关系 | 用相对误差 1e-9 量级比较,不要用 == |
| 包含关系 | A ∩ B 应被 A 与 B cover | A.covers(inter)、B.covers(inter) |
| 不重叠 | 结果各多边形内部不重叠 | result.isValid() 已隐含;覆盖数据可用 CoverageValidator |
9.1.4 边界归属约定
JTS 遵循 OGC Simple Features 的边界约定:一个几何由内部(interior)、边界(boundary)、外部(exterior)三部分组成,三者互不相交。对叠加运算是封闭点集运算,边界点的归属由运算本身决定:
- 交集包含 A 与 B 的公共边界(例如两个相邻多边形的公共边属于交集,哪怕交集面积为 0,结果会是一条线)。
- 并集只包含边界一次,公共边在结果中不会重复出现。
- 差集保留 A 的完整闭包,再减去 B 的闭包;A 与 B 的公共边界会从结果中移除。
这解释了为什么 A.difference(B) 在 A、B 相邻时得到的是“去掉公共边之后”的多边形——公共边既不属于 A − B 的内部,也不属于它的边界。若需要保留这条边(例如做地块合并时保留历史界址线),应改用 LineMerger 或在线状图层上另行处理,而不能指望差集。
9.2 顶层 API:Geometry 的四个方法
9.2.1 方法签名与返回约定
四个叠加方法都定义在 Geometry 抽象基类上,签名一致——接受一个 Geometry 参数,返回一个新的 Geometry,不修改任何一方:
import org.locationtech.jts.geom.Geometry;
// 交集
public Geometry intersection(Geometry other);
// 并集(二元)
public Geometry union(Geometry other);
// 差集:this − other
public Geometry difference(Geometry other);
// 对称差
public Geometry symDifference(Geometry other);
另有三个一元/便捷方法容易混淆,先在此厘清:
union()无参版本:对本几何的所有组成元素求并,允许输入是GeometryCollection,返回结果一定是最高的输入维度(多边形输入一定得到Polygonal)。详见 9.9。buffer(double distance):缓冲区,见第 10 章。它常被当成“修复工具”使用(geom.buffer(0)),但不是叠加运算。convexHull():凸包,见第 12 章。
返回值类型遵循“最小充分类型”原则:结果是单个多边形就返回 Polygon,多个不相连的多边形返回 MultiPolygon,混合维度(例如一个面加一条线)返回 GeometryCollection。写代码时不要假设返回类型,用 getGeometryType() 或 instanceof 判断。
9.2.2 第一个完整示例
下面这段代码演示四种运算,并打印结果的类型、面积、周长与元素个数。输出值来自实际运行,可直接对照:
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;
public class OverlayBasic {
public static void main(String[] args) throws Exception {
WKTReader reader = new WKTReader();
Geometry a = reader.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
Geometry b = reader.read("POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))");
report("A ∩ B", a.intersection(b));
report("A ∪ B", a.union(b));
report("A − B", a.difference(b));
report("B − A", b.difference(a));
report("A △ B", a.symDifference(b));
}
static void report(String label, Geometry g) {
System.out.println(label + " => " + g);
System.out.printf(" type=%s dim=%d area=%.1f length=%.1f num=%d%n",
g.getGeometryType(), g.getDimension(), g.getArea(),
g.getLength(), g.getNumGeometries());
}
}
运行结果:
A ∩ B => POLYGON ((5 10, 10 10, 10 5, 5 5, 5 10))
type=Polygon dim=2 area=25.0 length=20.0 num=1
A ∪ B => POLYGON ((10 5, 10 0, 0 0, 0 10, 5 10, 5 15, 15 15, 15 5, 10 5))
type=Polygon dim=2 area=175.0 length=60.0 num=1
A − B => POLYGON ((10 5, 10 0, 0 0, 0 10, 5 10, 5 5, 10 5))
type=Polygon dim=2 area=75.0 length=40.0 num=1
B − A => POLYGON ((5 10, 5 15, 15 15, 15 5, 10 5, 10 10, 5 10))
type=Polygon dim=2 area=75.0 length=40.0 num=1
A △ B => MULTIPOLYGON (((10 5, 10 0, 0 0, 0 10, 5 10, 5 5, 10 5)), ((10 5, 10 10, 5 10, 5 15, 15 15, 15 5, 10 5)))
type=MultiPolygon dim=2 area=150.0 length=80.0 num=2
可以验证两个细节:并集的周长 60 而不是 40+40=80,因为重叠部分的边被溶解了;对称差的周长 80 正好是 A、B 两个“L 形”的边界之和,因为两者在 (10,5)–(5,10) 处仅以点相接,不共享边。
9.2.3 空几何的短路规则
GeometryOverlay 在进入算法前先处理空几何,避免无谓的图构建。规则可以总结为:
| 场景 | 交集 | 并集 | 差集 | 对称差 |
|---|---|---|---|---|
| 两者都空 | 空(维度取输入较低者) | 空(维度取输入较高者) | 空 | 空 |
| A 空、B 非空 | 空 | B 的副本 | 空 | B 的副本 |
| A 非空、B 空 | 空 | A 的副本 | A 的副本 | A 的副本 |
“副本”意味着返回的是 copy() 而非原对象引用,但坐标与结构完全相同。空结果的维度规则由 OverlayOp.createEmptyResult / OverlayUtil.resultDimension 决定,例如 POLYGON EMPTY 与 POINT EMPTY 求交,得到的是 LINESTRING EMPTY 之外的哪一种,取决于输入维度而非运算。实际写代码时建议统一用 result.isEmpty() 判断,不要依赖空结果的类型。
9.2.4 GeometryCollection 的限制
四个方法对 GeometryCollection 的态度并不一致,这是实践中很容易踩的坑。用第 04 章的思路理解:JTS 的叠加引擎要求输入“同质”(homogeneous),即一个 MultiPolygon 里的元素必须都是多边形,一个 MultiLineString 里必须都是线。
- 二元并、差、对称差:遇到非空
GeometryCollection直接抛IllegalArgumentException,消息为Operation does not support GeometryCollection arguments。 - 交集:
GeometryOverlay.intersection对GeometryCollection做了特殊处理——把集合拆开,逐元素与另一个几何求交,再用GeometryCollectionMapper合并。所以它“能跑”,但结果是逐元素相交的拼接,可能包含大量低维碎片。 - 一元并集
union():完整支持异构GeometryCollection,是清理多类型数据的推荐入口(见 9.9)。
下面这段代码演示三种行为。第二个 catch 分支的异常消息是实际输出,可以据此在日志里快速识别问题类型:
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;
public class OverlayWithCollection {
public static void main(String[] args) throws Exception {
WKTReader reader = new WKTReader();
Geometry poly = reader.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
Geometry gc = reader.read(
"GEOMETRYCOLLECTION (POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0)), LINESTRING (20 20, 30 30))");
// 交集:逐元素映射,能算
System.out.println("intersection => " + poly.intersection(gc));
// 一元并集:支持异构集合
System.out.println("unary union => " + gc.union());
// 二元并集:拒绝 GeometryCollection
try {
poly.union(gc);
} catch (RuntimeException ex) {
System.out.println("union failed => " + ex.getClass().getSimpleName()
+ ": " + ex.getMessage());
}
}
}
输出:
intersection => POLYGON ((10 0, 0 0, 0 10, 10 10, 10 0))
unary union => GEOMETRYCOLLECTION (LINESTRING (20 20, 30 30), POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0)))
union failed => IllegalArgumentException: Operation does not support GeometryCollection arguments
一个实用的规避手法是:先按维度把集合拆成面层、线层、点层(PolygonExtracter、LineStringExtracter、PointExtracter 或 GeometryExtracter),分别在各自图层上做叠加,最后按需合并。这既符合业务语义,也避开了异构输入的限制。
9.3 两套引擎与全局开关
9.3.1 为什么会有两套引擎
JTS 的叠加分析最初只有一套实现,就是今天的“旧引擎”OverlayOp。它建立在 geomgraph 包的平面图结构上:把所有输入的边插入一个带标签的平面图,做 noding(求交打断),再按标签规则遍历子图、收集结果环。这套代码从 JTS 1.x 一直用到 1.17,支撑了 PostGIS、GeoTools、Hibernate Spatial 等一大票下游项目,稳定性经过了长期验证。
但它有两个绕不过去的短板。第一,精度问题:旧引擎的 noding 在双精度浮点上做,两条线段“数学上应该相交在一点”,浮点计算却可能给出两个略有差异的交点,于是 noding 不完整(non-noded intersection),最后在构面时抛 TopologyException。第二,性能与可维护性:平面图结构推理困难,遇到塌缩边、重合边时要靠大量启发式补丁。
1.18 引入的 OverlayNG(NG = Next Generation)从设计上就是把“精度”摆在第一位的重写。它的核心决策是:noding 必须在一个明确的精度模型下进行。若精度模型是固定的(PrecisionModel.FIXED),就用 Snap-Rounding Noder,把交点吸附到网格上,从根上保证 noding 完整;若精度模型是浮点的,就退回到普通的 MCIndexNoder,用校验器兜底,仍可能抛异常。
9.3.2 GeometryOverlay:隐藏的分发器
Geometry 的四个方法并不直接调用算法,而是委托给包级私有类 org.locationtech.jts.geom.GeometryOverlay。这个类做三件事:空几何短路、GeometryCollection 特殊处理、按全局开关选择引擎。它的核心分支如下:
Geometry.union(a, b)
│
▼
┌───────────────────────┐
│ GeometryOverlay │
│ (geom 包,包级私有) │
└───────────┬───────────┘
│ isOverlayNG ?
┌─────────────┴──────────────┐
false │ │ true
▼ ▼
┌──────────────────────┐ ┌───────────────────────┐
│ SnapIfNeededOverlayOp│ │ OverlayNGRobust │
│ ├─ OverlayOp(旧) │ │ └─ OverlayNG + 重试 │
│ └─ SnapOverlayOp │ └───────────────────────┘
└──────────────────────┘
开关由系统属性 jts.overlay 控制,取值只有两个:old(默认)与 ng。源码中对应 GeometryOverlay.OVERLAY_NG_DEFAULT = false,也就是说在 1.21.0-SNAPSHOT 中,高层 union/intersection/difference/symDifference 默认仍然走旧引擎。这一点与很多教程中“JTS 已全面切换到 OverlayNG”的说法不一致,写代码时必须自己确认。
可以用一行代码在启动时切换:
# 让所有 Geometry 高层叠加方法使用 OverlayNG(带鲁棒重试)
java -Djts.overlay=ng -cp jts-core-1.20.0.jar:app.jar com.example.Main
# 显式回退到旧引擎(也是默认值)
java -Djts.overlay=old -cp jts-core-1.20.0.jar:app.jar com.example.Main
如果你不想依赖全局状态(在多租户服务、并发测试里尤其危险),就不要碰这个属性,而是在需要的地方显式调用 OverlayNG 或 OverlayNGRobust 的静态方法。这样每个调用点的行为都是确定的。
9.3.3 引擎对照表
| 维度 | 旧引擎 OverlayOp |
新引擎 OverlayNG |
|---|---|---|
| 引入版本 | JTS 1.0 时代 | 1.18(2020-12) |
| 默认启用 | 是(jts.overlay=old) |
否,需 jts.overlay=ng 或显式调用 |
| 精度控制 | 跟随输入的 PrecisionModel,无独立开关 |
可独立传入 PrecisionModel,与输入无关 |
| noding 实现 | MCIndexNoder + RobustLineIntersector |
固定精度用 SnapRoundingNoder,浮点用校验过的 MCIndexNoder |
| 鲁棒性 | 依靠 SnapIfNeededOverlayOp 重试一次 |
依靠 OverlayNGRobust 五级重试阶梯 |
| 结果语义 | 非严格(可能混合维度) | 默认非严格,可 setStrictMode(true) |
| GeometryCollection 输入 | 二元运算不支持 | 支持“简单”GC(同质、无重叠多边形) |
| Z 值处理 | 不插值,保留首个来源坐标 | ElevationModel 沿边线性插值(1.18 起) |
| 结构 | geomgraph 平面图 |
轻量 Edge/OverlayEdge 图 |
| 可扩展点 | 几乎没有 | 可注入 Noder(CoverageUnion 即如此) |
9.3.4 该选哪个引擎
实践中的建议可以归纳成三条。第一,新代码优先显式使用 OverlayNGRobust:它把重试逻辑内置好了,输入再脏也有很大概率算出一个合法结果,代价是失败时比一次 OverlayNG 慢(因为要重试)。第二,批量数据处理优先考虑固定精度:把坐标按业务精度(例如 0.001 米)建立 PrecisionModel,一次性舍入后再算,既快又可复现。第三,不要混用:一个数据集里如果一半走旧引擎一半走新引擎,两者对塌缩边的处理规则不同,合并结果时会出现难以解释的缝隙与重叠。
至于结果差异,后面 9.12 节会用同一组输入给出对照。这里先记住一个结论:对“正常”的输入(本身合法、无近重合边),两套引擎的结果在拓扑上等价,差异主要体现在顶点坐标的舍入与环的起点方向;对“病态”输入,差异可能直接表现为一个能算、一个抛异常。
9.4 旧引擎 OverlayOp 的流程
9.4.1 流水线总览
虽然新项目不建议主动使用旧引擎,但读懂它的流程有两个好处:一是看到 TopologyException 时能判断异常从哪一步冒出来;二是理解 SnapIfNeededOverlayOp 为什么那样重试。旧引擎的主体在 org.locationtech.jts.operation.overlay.OverlayOp,流程大致如下:
输入 geom0, geom1, opCode
│
▼
① OverlayNodeFactory + PlanarGraph
│ 把两侧的边与节点加入图
▼
② EdgeNodingBuilder / EdgeSetNoder
│ 求交打断(noding),生成带 Label 的边
▼
③ OverlayOp.computeOverlay()
│ 按 opCode 计算每个节点的标签
▼
④ 遍历 MaximalEdgeRing → MinimalEdgeRing
│ 收集属于结果的面
▼
⑤ PointBuilder / LineBuilder / PolygonBuilder
│ 分别构造结果中的点、线、面
▼
buildGeometry(...) → 结果
其中第 ② 步是异常高发区。旧引擎用 RobustLineIntersector 计算线段交点,再把交点插入到边链中;如果同一对线段被算出两个“几乎相同但不相等”的交点,插入后图里就会出现两条本该相交却错开的边,第 ④ 步构环时便无法闭合,于是抛出:
org.locationtech.jts.geom.TopologyException: found non-noded intersection between
LINESTRING ( 0 0, 10 10 ) and LINESTRING ( 10 0, 0 10 ) [ (5.0, 5.0, NaN) ]
这条消息是经典款,几乎每个用过 JTS 的人都见过。它告诉你:在 (5,5) 附近,两条线段本应相交但 noding 没成功。注意方括号里给出的是引擎认为“应该存在交点”的位置估计,实际计算出的交点可能差了几个 ulp。
9.4.2 Self-intersecting 输入触发的异常
最容易被复现的异常来源是自相交多边形(俗称蝴蝶结,bow-tie)。下面这段代码直接对蝴蝶结多边形做交运算,异常消息就是上一节展示的那条:
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;
public class BowtieOverlay {
public static void main(String[] args) throws Exception {
WKTReader reader = new WKTReader();
Geometry bowtie = reader.read("POLYGON ((0 0, 10 10, 10 0, 0 10, 0 0))");
Geometry square = reader.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
System.out.println("bowtie isValid = " + bowtie.isValid());
try {
System.out.println(bowtie.intersection(square));
} catch (RuntimeException ex) {
System.out.println("intersection failed: " + ex.getMessage());
}
// buffer(0) 是常用的“合法化”手段
Geometry fixed = bowtie.buffer(0);
System.out.println("buffer(0) => " + fixed);
System.out.println("fixed isValid = " + fixed.isValid());
}
}
输出:
bowtie isValid = false
intersection failed: found non-noded intersection between LINESTRING ( 0 0, 10 10 ) and LINESTRING ( 10 0, 0 10 ) [ (5.0, 5.0, NaN) ]
buffer(0) => POLYGON ((5 5, 10 10, 10 0, 5 5))
fixed isValid = true
这里有两个值得记住的点。其一,buffer(0) 只保留“最大面积的那一瓣”,所以蝴蝶结被削成了一个三角形(面积 25),另一瓣被丢弃——这不是通用的修复方案,第 17 章会给出更合适的 GeometryFixer。其二,SnapIfNeededOverlayOp 的重试在这里没能救场:它会用 GeometrySnapper 的容差重新吸附后再算一遍,但这个容差(由 GeometrySnapper.computeSizeBasedSnapTolerance 给出,等于输入包络短边 × 1e-9,本例为 1e-8)远小于蝴蝶结两瓣的交叠尺度(这里是 5 个单位),吸附后依然是同一条自交边,所以第二次尝试照样失败,最终把原始异常抛给调用者。
9.4.3 SnapIfNeededOverlayOp 的重试机制
SnapIfNeededOverlayOp 的逻辑非常直白,可以概括为“一次正式 + 一次补救”:
// 伪代码,对应 SnapIfNeededOverlayOp.getResultGeometry
try {
result = OverlayOp.overlayOp(g0, g1, opCode); // 直接用原始坐标算
isSuccess = true;
} catch (RuntimeException ex) {
savedException = ex; // 记住原始异常
}
if (!isSuccess) {
try {
result = SnapOverlayOp.overlayOp(g0, g1, opCode); // 吸附后再算
} catch (RuntimeException ex) {
throw savedException; // 还是失败就抛最初那个
}
}
SnapOverlayOp 内部做了三件事:先用 CommonBitsRemover 去掉两个几何的公共高位(把大坐标平移成小坐标,减少精度损失),再用 GeometrySnapper.snap 把两边的顶点互相吸附,然后交给 OverlayOp 计算,最后把公共位加回去。吸附容差由 GeometrySnapper.computeOverlaySnapTolerance(g0, g1) 给出:取两个几何各自的“基于尺寸容差”中较小者,而单个几何的基于尺寸容差是 min(width, height) × 1e-9;若几何使用固定精度模型,则容差至少为 (1/scale) × 2 / 1.415。
这解释了为什么旧引擎对“差一点点就重合”的输入很有效(例如两个测绘单位分别数字化同一条界址线,顶点偏差 1e-7 米),而对“尺度级别”的病态输入无能为力——容差是相对尺度的,不是相对错误的。
| 阶段 | 做什么 | 容差来源 | 失败表现 |
|---|---|---|---|
| 第一次尝试 | 原始坐标直接叠加 | 无(用输入 PrecisionModel) |
抛 TopologyException |
| 第二次尝试 | 去公共位 → 顶点吸附 → 叠加 → 加回公共位 | min(w,h) × 1e-9,固定精度另加下限 |
抛原始异常 |
这条重试链的能力上限,正好解释了 OverlayNGRobust 存在的意义:它把“吸附一次”扩展为“吸附多次、容差递增、最后退回 Snap-Rounding”,覆盖了旧引擎够不着的场景。
9.5 OverlayNG 架构
9.5.1 五阶段流水线
OverlayNG 把叠加拆成五个阶段,每个阶段职责单一、数据形态明确。理解这条流水线,是读懂 TopologyException、读懂“为什么固定精度更慢但更稳”的关键。
geom0, geom1, pm, opCode
│
▼
┌──────────────────────────────────────────────┐
│ 1. nodeEdges() —— EdgeNodingBuilder │
│ · 抽取边(RingClipper/LineLimiter 可裁剪) │
│ · 按 pm 选 Noder:SnapRounding / MCIndex │
│ · 丢弃完全塌缩的边,成 Edge 并按 EdgeKey 合并│
└───────────────────┬──────────────────────────┘
▼ List<Edge>(唯一、已 noding)
┌──────────────────────────────────────────────┐
│ 2. buildGraph() —— OverlayGraph │
│ 每条 Edge → 一对反向 OverlayEdge(共享 Label)│
└───────────────────┬──────────────────────────┘
▼
┌──────────────────────────────────────────────┐
│ 3. labelGraph() —— OverlayLabeller │
│ computeLabelling() 求各边两侧 Location │
│ markResultAreaEdges(opCode) 标记结果边 │
└───────────────────┬──────────────────────────┘
▼
┌──────────────────────────────────────────────┐
│ 4. extractResult() —— 按维度分层构结果 │
│ PolygonBuilder → LineBuilder → PointBuilder│
└───────────────────┬──────────────────────────┘
▼
┌──────────────────────────────────────────────┐
│ 5. ElevationModel.populateZ()(仅当有 Z) │
└───────────────────┬──────────────────────────┘
▼
result
阶段 1 的产物是一组唯一的边:两条完全重合的边(比如两个相邻多边形共享的公共边)会被 EdgeMerger 合并成一条,并在标签里记录“它同时来自几何 0 与几何 1”。这一步决定了结果的正确性下限——如果 noding 没把该断开的边断开,后面所有步骤都白搭。
阶段 3 的产物是每条边两侧的 Location(INTERIOR / BOUNDARY / EXTERIOR)。OverlayLabeller 从一个已知位置的边出发,沿着节点周围的边按角度顺序推进,逐步推导出整张图所有边的位置。这是叠加分析里最“几何”的一步,也是性能瓶颈之一。
阶段 4 按运算类型决定要构哪些维度的结果。源码中 extractResult 的规则值得记住:先构面,若没有面再构线,若连线都没有且运算是交集,才构点。也就是说,交集的面优先于线、线优先于点,这是 isAllowMixedIntResult = !isStrictMode 与“最高维度优先”语义共同作用的结果。
9.5.2 Edge 与 OverlayLabel:四种边的状态
Edge 是“可见的一段线”,OverlayLabel 是“这段线相对于两个输入几何的拓扑身份”。标签用 dim 字段记录边的状态,只有四种取值:
| 状态常量 | 含义 | 记录的信息 |
|---|---|---|
DIM_BOUNDARY |
输入面的边界边 | 左右两侧 locLeft/locRight、该边属于外环还是内环(isHole) |
DIM_COLLAPSE |
输入面因吸附/舍入而塌缩的边 | 相对于有效面的位置 locLine(尖刺为 EXTERIOR,缝隙为 INTERIOR) |
DIM_LINE |
输入线的边 | locLine,初始为未知,后续定位 |
NOT_PART |
不属于该输入几何 | 无额外信息 |
这四种状态之所以够用,靠两条前提:输入必须同质(一条边不可能在同一个几何里既是塌缩边又是线边);两个输入几何各自的维度可以不同,所以同一条边可以在几何 0 里是边界边、在几何 1 里是线边。
一条边在两个输入几何中的标签组合(示意)
┌──────────────┬──────────────┬─────────────────────────┐
│ 相对 geom0 │ 相对 geom1 │ 在结果中的角色 │
├──────────────┼──────────────┼─────────────────────────┤
│ BOUNDARY/壳 │ NOT_PART │ A 独有的外边 → 并集保留 │
│ BOUNDARY/壳 │ BOUNDARY/壳 │ 重合边 → 按 opCode 决定 │
│ BOUNDARY/壳 │ COLLAPSE │ A 的边+B 的塌缩 → 需判定 │
│ COLLAPSE │ COLLAPSE │ 双塌缩 → 通常丢弃 │
│ LINE │ NOT_PART │ B 输入是面时 → 按位置判 │
└──────────────┴──────────────┴─────────────────────────┘
OverlayEdge 则是一对方向相反的半边的容器:一条 Edge 生成两条 OverlayEdge(正向与反向),它们共享同一个 OverlayLabel。这样“边的左侧”与“边的右侧”就分别对应两条半边各自的视角,避免了左右信息在合并时丢失。
9.5.3 裁剪与短路优化
OverlayNG 默认开启优化(isOptimized = true)。优化分两类。
第一类是包络裁剪。交集的结果不可能超出两个输入包络的交集,差集的结果不可能超出 A 的包络,所以在这两种运算下,OverlayUtil.clippingEnvelope 会算出一个安全的裁剪包络,交给 RingClipper 把面裁剪进去、交给 LineLimiter 把线截断。安全外扩量是浮点精度下包络短边的 10%、固定精度下 3 倍网格边长,确保舍入后的顶点仍落在包络内。并集与对称差的结果范围就是两输入的全范围,无法裁剪,此时返回 null,不做优化。
第二类是结果面积启发式检查。浮点精度下,OverlayNG 在返回结果前会用 OverlayUtil.isResultAreaConsistent 校验一次:交集的面积不应大于任一输入,并集的面积不应小于任一输入,差集的面积应满足上下界关系,容差系数是 0.1。若明显违反,说明 noding 让图“翻转”了,直接抛 TopologyException("Result area inconsistent with overlay operation")。这条检查是 1.18 之后新增的,它把一部分原本静默的错误结果变成了显式异常,是好事。
FastOverlayFilter 是另一层短路:若两个几何的包络根本不相交,或其中一个完全包含另一个,则不必进入完整流水线。它只在“结果必然为空或必然等于某一输入”时才生效,属于纯加速手段。
9.5.4 用代码观察中间产物
OverlayNG 提供了三个调试用输出开关:setOutputNodedEdges(true) 输出 noding 之后的全部边,setOutputEdges(true) 输出进入图的边,setOutputResultEdges(true) 输出被标记为结果的边。它们返回的是 MultiLineString 而不是最终几何,非常适合排查“为什么这里多出一条线”。
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.overlayng.OverlayNG;
public class OverlayDebugEdges {
public static void main(String[] args) throws Exception {
WKTReader reader = new WKTReader();
Geometry a = reader.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
Geometry b = reader.read("POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))");
OverlayNG union = new OverlayNG(a, b, null, OverlayNG.UNION);
union.setOutputNodedEdges(true);
System.out.println("noded edges = " + union.getResult());
OverlayNG inter = new OverlayNG(a, b, null, OverlayNG.INTERSECTION);
inter.setOutputResultEdges(true);
System.out.println("result edges = " + inter.getResult());
// 关掉裁剪优化,对比结果是否变化
OverlayNG noOpt = new OverlayNG(a, b, null, OverlayNG.INTERSECTION);
noOpt.setOptimized(false);
System.out.println("no-opt result = " + noOpt.getResult());
}
}
输出(节选,noded edges 共 12 段):
noded edges = MULTILINESTRING ((0 0, 10 0, 10 5), (10 5, 10 0, 0 0), (10 5, 10 10, 5 10), ... , (5 5, 5 10))
result edges = MULTILINESTRING ((5 10, 10 10, 10 5), (10 5, 5 5), (5 5, 5 10))
no-opt result = POLYGON ((5 10, 10 10, 10 5, 5 5, 5 10))
注意 noded edges 里每条公共边都出现了两次、方向相反——这正是“一对反向 OverlayEdge”在调试输出里的直观体现。而 result edges 只剩三条边,正好围成那个 5×5 的交集方块。关掉优化后结果与开启优化完全一致,说明裁剪只影响性能、不影响语义(前提是安全包络计算正确)。
9.6 精度模型与 PrecisionReducer
9.6.1 浮点与固定精度
PrecisionModel 有两种类型:FLOATING(默认,直接用 double)与 FIXED(按固定比例尺舍入到网格)。两者的差异在叠加分析里会被放大,因为它决定了 noding 策略:
| 精度类型 | 构造方式 | noding 实现 | 鲁棒性 | 输出坐标 |
|---|---|---|---|---|
| FLOATING | new PrecisionModel() 或几何默认 |
MCIndexNoder + 校验 |
可能失败,抛异常 | 保留计算值 |
| FIXED | new PrecisionModel(scale) |
SnapRoundingNoder |
在 13 位有效数字内可靠 | 舍入到 1/scale 网格 |
scale 的含义是“每单位长度分成多少格”,例如 new PrecisionModel(1000) 表示保留到 0.001。源码注释里有一条重要限制:固定精度的可靠计算上限约 13 位十进制有效数字(PrecisionUtil.MAX_ROBUST_DP_DIGITS = 14,留了一位余量)。因为坐标量级越大,可用于小数的有效位越少——量级 1e6 的坐标再乘 1e9 的 scale,就已经逼近 double 的 15 到 16 位有效数字极限。
9.6.2 三个精度工具函数
PrecisionUtil 提供了一组静态方法,作用是把“经验”变成“可计算的数”。最容易混淆的是 safeScale 与 robustScale:
| 方法 | 计算依据 | 典型返回值 | 用途 |
|---|---|---|---|
safeScale(value) / safeScale(geom) / safeScale(a,b) |
坐标最大量级,保证舍入后不超过 14 位有效数字 | 1.0E12 |
OverlayNGRobust 最后一级 Snap-Rounding 使用 |
inherentScale(geom) |
扫描全部坐标,取能完整表示现有小数位的最小 scale | 1.0 |
判断数据“本身”的精度,文档标注 “very slow” |
robustScale(a,b) / robustPM(a,b) |
在安全前提下尽量保留输入精度 | 1.0 |
需要“既保留输入精度又不溢出”时使用 |
对本章那对 0 到 15 量级的正方形,safeScale 给出 1e12,意味着“在这个量级下,保留 12 位小数仍然安全”;而 inherentScale 只给出 1.0,因为所有坐标都是整数——数据本身没有小数位,没必要用更高精度。robustScale 取两者中更保守的一个,所以得到 1.0。
这三个函数的取舍逻辑值得展开。safeScale 只关心“会不会溢出精度”,所以对整数坐标也会给出 1e12;inherentScale 只关心“数据长什么样”,所以浮点噪声也会被当成有效精度;robustScale 则是“在安全的前提下尽量保留输入精度”。实际工程里,如果你明确知道业务精度是 0.001 米,直接用 new PrecisionModel(1000) 比任何启发式都可靠。
9.6.3 固定精度叠加的实际效果
下面这段代码对同一对“歪斜”多边形分别做浮点叠加与两种固定精度叠加,并打印面积。输入刻意带了 9 位小数,方便观察舍入效果:
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.overlayng.OverlayNG;
public class OverlayPrecision {
public static void main(String[] args) throws Exception {
WKTReader reader = new WKTReader();
Geometry c = reader.read(
"POLYGON ((0.123456789 0.1, 10.987654321 0.2, 10.5 10.4, 0.3 9.9, 0.123456789 0.1))");
Geometry d = reader.read(
"POLYGON ((5.111111111 5.2, 15.2 5.3, 15.4 15.6, 5.5 15.1, 5.111111111 5.2))");
show("float ", OverlayNG.overlay(c, d, OverlayNG.INTERSECTION));
show("fixed 1 ", OverlayNG.overlay(c, d, OverlayNG.INTERSECTION, new PrecisionModel(1.0)));
show("fixed 10 ", OverlayNG.overlay(c, d, OverlayNG.INTERSECTION, new PrecisionModel(10.0)));
}
static void show(String label, Geometry g) {
System.out.println(label + " => " + g);
System.out.printf(" area=%.15f numPts=%d%n", g.getArea(), g.getNumPoints());
}
}
输出:
float => POLYGON ((5.305373344890099 10.145361438475005, 10.5 10.4, 10.745937852375098 5.2558518069072875, 5.111111111 5.2, 5.305373344890099 10.145361438475005))
area=27.319979782567877 numPts=5
fixed 1 => POLYGON ((5 10, 11 10, 11 5, 5 5, 5 10))
area=30.000000000000000 numPts=5
fixed 10 => POLYGON ((5.3 10.1, 10.5 10.4, 10.7 5.3, 5.1 5.2, 5.3 10.1))
area=27.000000000000000 numPts=5
三条结论可以直接写进工程规范。其一,固定精度的输出坐标全部落在网格上(5.3、10.7),这正是“结果可复现”的基础:同样的输入、同样的 scale,任何机器上都得到同一份 WKT。其二,精度越低,面积误差越大:scale=1(整米网格)把面积从 27.32 拉到了 30.00,误差接近 10%;scale=10(0.1 米)得到 27.00,误差约 1.2%。其三,浮点结果虽然“最接近真值”,但顶点坐标带着 15 位小数的噪声,做后续比较、入库、切片都会引入不确定性。
9.6.4 PrecisionReducer:只降精度、不改变拓扑意图
PrecisionReducer.reducePrecision(geom, pm) 是“把几何舍入到指定精度”的专用工具。它内部并不是简单地遍历坐标四舍五入,而是用 OverlayNG 的单输入模式(new OverlayNG(geom, pm),opCode 固定为 UNION)在舍入后的精度下重新计算一遍——这样即使舍入导致自相交或相邻边粘连,输出也是合法几何。对面积维度的输入,它会设置 setAreaResultOnly(true),避免塌缩出线状碎片。
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.overlayng.PrecisionReducer;
public class ReducePrecisionDemo {
public static void main(String[] args) throws Exception {
WKTReader reader = new WKTReader();
Geometry g = reader.read(
"POLYGON ((0.2 0.3, 9.7 0.6, 9.9 9.8, 0.1 9.4, 0.2 0.3))");
Geometry p1 = PrecisionReducer.reducePrecision(g, new PrecisionModel(1.0));
Geometry p10 = PrecisionReducer.reducePrecision(g, new PrecisionModel(10.0));
System.out.println("pm=1 => " + p1 + " area=" + p1.getArea());
System.out.println("pm=10 => " + p10 + " area=" + p10.getArea());
System.out.println("origin area=" + g.getArea());
}
}
输出:
pm=1 => POLYGON ((0 0, 0 9, 10 10, 10 1, 0 0)) area=90.0
pm=10 => POLYGON ((0.2 0.3, 0.1 9.4, 9.9 9.8, 9.7 0.6, 0.2 0.3)) area=88.28
origin area=88.28
pm=10 时面积没变(88.28 恰好等于原值,因为原坐标本来就只有一位小数),只是环的方向与起点变了——又一次说明“面积、顶点集合”比“WKT 字符串”更适合做断言。pm=1 时面积从 88.28 变成 90,四个顶点被分别吸到最近的整数格点,这是可接受的(也是可预期的)精度损失。
注意:
org.locationtech.jts.precision.GeometryPrecisionReducer是另一个类,功能更全(支持点、线、面,支持按元素处理),1.18 起内部同样改用OverlayNG与 Snap-Rounding。两者不要混用:PrecisionReducer只处理线面输入,遇到非法输入抛IllegalArgumentException。
9.7 OverlayNGRobust 的重试阶梯
9.7.1 为什么需要“重试”
OverlayNG 的浮点模式并不鲁棒:它用的 MCIndexNoder 在 noding 出错时只会让 ValidatingNoder 抛异常,不会自动修复。而现实数据里,noding 出错是常态——不同来源的图层边界共线但顶点不重合、坐标经过多次投影转换后带着细微抖动、从 Shapefile 读进来的双精度值本身就不可靠。于是 JTS 提供了 OverlayNGRobust:一个把“多次尝试、逐级放宽”封装好的门面类。
它的算法思想非常朴素:先按最保真的方式算,算不出来就牺牲一点精度再算,直到算出来为止。源码的类注释把策略列为五条,对应到实现里就是三级阶梯:
输入 geom0, geom1, opCode
│
▼
┌─────────────────────────────────────────┐
│ 第 1 级:OverlayNG + 浮点 noding │ 快,坐标改动最小
│ 失败 → 捕获原始异常保存起来 │
└───────────────────┬─────────────────────┘
▼
┌─────────────────────────────────────────┐
│ 第 2 级:SnappingNoder,最多 5 轮 │
│ 容差 tol = max(|ordinate|) / 1e12 │
│ 每轮:① 两几何互相吸附后叠加 │
│ ② 各自先自吸附、再叠加 │
│ 失败 → tol *= 10,重试 │
└───────────────────┬─────────────────────┘
▼
┌─────────────────────────────────────────┐
│ 第 3 级:SnapRoundingNoder + 安全比例尺 │
│ scale = PrecisionUtil.safeScale(a,b) │
│ 失败 → 抛出第 1 级保存的原始异常 │
└───────────────────┬─────────────────────┘
▼
结果 / 异常
9.7.2 三级策略的细节
第 1 级直接调用 OverlayNG.overlay(geom0, geom1, opCode),用几何自身的精度模型(通常是浮点)。成功就返回,失败则把异常暂存起来——注意是“暂存”而不是“记录并继续”,因为如果后面所有策略都失败,抛给调用者的应该是这个最原始的异常,它的消息里带着最接近真实的交点位置,最有诊断价值。
第 2 级用 SnappingNoder,容差由 snapTolerance 计算:取两个几何包络中绝对值最大的 ordinate,除以 SNAP_TOL_FACTOR = 1e12。以本章的正方形为例,最大 ordinate 是 15,容差就是 1.5e-11。每一轮做两次尝试:先让两个几何互相吸附后叠加(overlaySnapping),再让每个几何先自吸附、然后叠加(overlaySnapBoth)。都失败就把容差乘 10 进入下一轮,最多 NUM_SNAP_TRIES = 5 轮。容差从 1e-11 一路放大到 1e-7,覆盖了“浮点噪声”到“数据录入误差”的大部分场景。
第 3 级退回 Snap-Rounding:用 PrecisionUtil.safeScale(a, b) 算出一个保证不超过 14 位有效数字的比例尺,建 PrecisionModel,再调一次固定精度的 OverlayNG。这一级会把坐标吸附到网格上,是最“暴力”也最可靠的一级。
值得强调的是:第 2 级的自吸附(snapSelf)内部用 OverlayNG(geom, null) 加 setStrictMode(true),也就是“把几何与空几何做并集”。这相当于一次“自我 noding 加去重”,能把几何内部那些“几乎重合的相邻边”先粘在一起,从而消除大部分自相交类问题。这个技巧在工程上也可以单独使用:对一个脏的多边形执行 new OverlayNG(geom, null).getResult(),往往就能得到合法几何。
9.7.3 什么情况下它救不了你
OverlayNGRobust 并非万能。它救不了的三类情况:
- 输入真正无效且错误尺度大于容差。两个多边形大面积重叠、自相交的蝴蝶结、环方向错误导致的“面积在外部”,这些不是精度问题而是数据问题,容差放大一万倍也没用。
- 输入维度混合。
OverlayNG明确拒绝异构GeometryCollection,异常消息是Overlay input is mixed-dimension;这类参数错误不会因为重试而消失。 - 拓扑语义本身不成立。例如对一对“一边完全在另一个内部”的多边形求差,结果为空是正确行为,不应期待异常;反之,如果你期待得到非空结果而得到空,那是业务假设错了,不是引擎错了。
9.7.4 用代码观察三个入口的差异
OverlayNGRobust 没有暴露“第几级成功”的回调,但可以从结果结构上间接观察。下面这段代码对一对“顶点有微小偏差”的相邻多边形做并集,对比三个入口的返回结构:
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.overlayng.OverlayNG;
import org.locationtech.jts.operation.overlayng.OverlayNGRobust;
public class RobustRetryDemo {
public static void main(String[] args) throws Exception {
WKTReader reader = new WKTReader();
// 两个本应相邻的多边形,接缝处顶点偏差 1e-9
Geometry a = reader.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
Geometry b = reader.read(
"POLYGON ((10.000000001 0, 20 0, 20 10, 10.000000001 9.999999999, 10.000000001 0))");
report("OverlayNG.union ", OverlayNG.overlay(a, b, OverlayNG.UNION));
report("OverlayNGRobust.union", OverlayNGRobust.overlay(a, b, OverlayNG.UNION));
report("Geometry.union ", a.union(b));
}
static void report(String label, Geometry g) {
System.out.printf("%s => num=%d area=%.9f%n",
label, g.getNumGeometries(), g.getArea());
}
}
输出:
OverlayNG.union => num=1 area=200.000000000
OverlayNGRobust.union => num=1 area=200.000000000
Geometry.union => num=2 area=200.000000001
三个入口都能给出结果,但 Geometry.union 走的是旧引擎(默认 jts.overlay=old),结果是两个互不相连的多边形(num=2,因为 1e-9 的缝隙没有被识别为重合边);而 OverlayNG 系列把缝隙正确地消除成了 num=1。面积差异只在小数点后第 9 位,肉眼几乎不可见,但几何结构完全不同——这正是“精度决策会改变拓扑结论”的典型例证。
如果想手工复现第 3 级的效果,可以直接用 PrecisionUtil.safeScale 构造安全比例尺:
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.overlayng.OverlayNG;
import org.locationtech.jts.operation.overlayng.PrecisionUtil;
public class SafeScaleOverlay {
public static void main(String[] args) throws Exception {
WKTReader reader = new WKTReader();
Geometry a = reader.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
Geometry b = reader.read("POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))");
double scale = PrecisionUtil.safeScale(a, b);
System.out.println("safeScale = " + scale);
System.out.println("robustScale = " + PrecisionUtil.robustScale(a, b));
Geometry g = OverlayNG.overlay(a, b, OverlayNG.INTERSECTION, new PrecisionModel(scale));
System.out.println("result = " + g);
}
}
输出:
safeScale = 1.0E12
robustScale = 1.0
result = POLYGON ((5 10, 10 10, 10 5, 5 5, 5 10))
9.8 严格模式与结果维度
9.8.1 非严格模式下的混合维度结果
默认情况下(isStrictMode = false),OverlayNG 允许结果混合维度。最典型的例子是交集:两个 MultiPolygon 相交,既有一块面积重合的区域,又有一对只在角点相接的多边形,那么结果就是一个 GeometryCollection,里面既有 Polygon 又有 Point:
A = MULTIPOLYGON( [0,10]×[0,10] , [20,30]×[0,10] )
B = MULTIPOLYGON( [5,15]×[5,15] , [30,40]×[10,20] )
A ∩ B :第一对相交于 5×5 的面(面积 25)
第二对仅在点 (30,10) 相接
→ 非严格模式:GEOMETRYCOLLECTION (POLYGON(25), POINT(30 10))
→ 严格模式 :POLYGON(25)
这个行为继承了旧引擎的语义,对某些分析场景有用(例如你想同时知道“面接触”和“点接触”的位置),但对数据处理管线是灾难:一个 GEOMETRYCOLLECTION 传进下一个算子,多数算子会抛异常或返回空。
9.8.2 严格模式的三条规则
调用 setStrictMode(true) 后,结果语义收紧为三条:
- 因拓扑塌缩产生的线、点不进入结果。例如两个多边形只共享一条边,交集在非严格模式下是
LINESTRING,严格模式下是空。 INTERSECTION与DIFFERENCE的结果一定是同质的(homogeneous)。- 若两个输入维度相同,
UNION与SYMDIFFERENCE的结果也同质。
严格模式带来的好处是“可链式调用”:结果永远是干净的面或干净的线,不需要每次都判断“这次会不会又有低维碎片”。代价是丢失信息,因此严格模式适合做数据加工流水线,不适合做空间关系诊断。
9.8.3 代码对比
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.overlayng.OverlayNG;
public class StrictModeDemo {
public static void main(String[] args) throws Exception {
WKTReader reader = new WKTReader();
Geometry mp = reader.read(
"MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0)), ((20 0, 30 0, 30 10, 20 10, 20 0)))");
Geometry mp2 = reader.read(
"MULTIPOLYGON (((5 5, 15 5, 15 15, 5 15, 5 5)), ((30 10, 40 10, 40 20, 30 20, 30 10)))");
Geometry loose = OverlayNG.overlay(mp, mp2, OverlayNG.INTERSECTION);
System.out.println("non-strict => " + loose);
System.out.println(" type=" + loose.getGeometryType()
+ " num=" + loose.getNumGeometries());
OverlayNG ng = new OverlayNG(mp, mp2, null, OverlayNG.INTERSECTION);
ng.setStrictMode(true);
Geometry strict = ng.getResult();
System.out.println("strict => " + strict);
System.out.println(" type=" + strict.getGeometryType()
+ " num=" + strict.getNumGeometries());
}
}
输出:
non-strict => GEOMETRYCOLLECTION (POLYGON ((5 10, 10 10, 10 5, 5 5, 5 10)), POINT (30 10))
type=GeometryCollection num=2
strict => POLYGON ((5 10, 10 10, 10 5, 5 5, 5 10))
type=Polygon num=1
注意严格模式并不是“扔掉所有低维元素”那么简单:如果交集中只有线、没有面(例如两条线交叉),严格模式下的交集结果会是空的 LineString,而不是保留那个交点。做“线网相交检测”时不要用叠加结果去判断,应该用第 08 章的 RelateNG 或 IntersectionMatrix。
9.8.4 与旧引擎语义的对照
| 场景 | 旧引擎 OverlayOp |
OverlayNG 默认 |
OverlayNG 严格模式 |
|---|---|---|---|
| 面-面只在边相接,求交 | LineString(保留) |
LineString(保留) |
空 |
| 面-面只在点相接,求交 | Point(保留) |
Point(保留) |
空 |
| 面-线求交 | LineString / MultiLineString |
同左 | 同左 |
| 多元素分别产生面与点 | GeometryCollection |
GeometryCollection |
只留面 |
| 塌缩产生的零面积多边形 | 可能保留为线 | 可能保留为线 | 丢弃 |
如果不方便构造 OverlayNG 实例,也可以全局切换:-Djts.overlay=ng 之后,Geometry 的高层方法会走 OverlayNGRobust,但它仍然是非严格模式。严格模式只能通过显式构造 OverlayNG 并调用 setStrictMode(true) 获得,没有全局开关。
9.9 一元并集家族
9.9.1 为什么需要专门的一元并集
把 1000 个相邻宗地合并成一个街道,最直接的做法是循环调用 a = a.union(b)。这个朴素实现有两个致命问题:第一,每次二元并集都会重建整张拓扑图,总代价是 O(n) 次图构建,远高于必要的量级;第二,中间结果越来越复杂,后续每次并集的输入规模都在膨胀,而实际上大部分内部边在最终结果里会被溶解掉。
JTS 为此提供了“一元并集”这一族 API,核心思路是先收集所有边界边,一次性做 noding,再统一成面,避免反复构建中间几何。几个类的分工如下:
| 类 | 包 | 输入 | 策略 | 适用场景 |
|---|---|---|---|---|
UnaryUnionOp |
operation.union |
集合或单个几何(可为异构 GC) | 按维度分层,面走级联并集,点走去重 | 通用入口,语义最完整 |
CascadedPolygonUnion |
operation.union |
仅多边形集合 | 按 STRtree 空间邻近性两两合并(级联) | 大量多边形的批量溶解 |
OverlayNGRobust.union |
operation.overlayng |
单个几何或集合 | UnaryUnionOp + OverlayNG 引擎 |
需要新引擎鲁棒性 |
UnaryUnionNG |
operation.overlayng |
集合 | 固定精度下一次算完 | 已知精度、追求可复现 |
CoverageUnion |
operation.overlayng |
合法覆盖(无重叠无缝隙) | 跳过 noding 与标注,直接拼接边界 | 覆盖面数据的高性能溶解 |
OverlapUnion |
operation.union |
两个几何 | 只在重叠区域做真正的叠加 | 大量“仅边界轻微重叠”的合并 |
9.9.2 UnaryUnionOp:默认入口
UnaryUnionOp.union(geom) 与 Geometry.union() 等价,是“把所有东西并起来”的默认选择。它的内部流程分三步:先用 InputExtracter 把输入按维度拆成面、线、点三组;面组交给 CascadedPolygonUnion,线组走 noding 加 LineMerger 合并,点组做去重;最后按维度从高到低拼装结果。
它的默认并集函数是 CascadedPolygonUnion.CLASSIC_UNION,也就是旧引擎。这一点在源码里写得很清楚:字段声明为 private UnionStrategy unionFunction = CascadedPolygonUnion.CLASSIC_UNION;,而 CLASSIC_UNION 的实现就是调 OverlayOp.overlayOp(g0, g1, OverlayOp.UNION)。所以,Geometry.union() 默认也走旧引擎。若要让一元并集走新引擎,有两个办法:显式调 OverlayNGRobust.union(geom),或者设置 -Djts.overlay=ng(此时 GeometryOverlay.union(a) 会转发到 OverlayNGRobust.union)。
9.9.3 级联并集的空间索引思想
CascadedPolygonUnion 的“级联”不是简单的两两迭代,而是把输入按空间邻近性组织:每次从集合里取一个多边形,在 STRtree 里查与它相交的候选,把候选取出来先合并成一个小结果,再放回集合继续。这样每次合并只处理“真正相邻”的几何,避免了大范围无效计算。
输入:N 个多边形
│
▼
┌────────────────────────────┐
│ STRtree(node capacity = 4)│
└─────────────┬──────────────┘
│ 取一个多边形 p
▼
┌────────────────────────────┐
│ 查询与 p 包络相交的候选集合 │
└─────────────┬──────────────┘
│ 候选逐个与 p 合并
▼
┌────────────────────────────┐
│ 合并结果放回集合,重复直到 1 │
└─────────────┬──────────────┘
▼
最终并集
这段逻辑在源码里的常量是 STRTREE_NODE_CAPACITY = 4,与第 13 章要讲的 STRtree 默认值一致。级联策略对小规模数据(几十个面)优势不明显,但在上万面的批量溶解中通常能带来数倍加速。
9.9.4 四者结果一致、性能不同
下面这段代码用同一组 400 个互相重叠的小方块,对比四个入口的结果面积与耗时。四个入口都给出面积 78400(因为方块按 1.4 倍间距铺排,重叠部分被正确溶解),差别只在耗时:
import java.util.ArrayList;
import java.util.List;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.operation.overlayng.OverlayNGRobust;
import org.locationtech.jts.operation.overlayng.UnaryUnionNG;
import org.locationtech.jts.operation.union.CascadedPolygonUnion;
import org.locationtech.jts.operation.union.UnaryUnionOp;
public class UnaryUnionCompare {
static final GeometryFactory GF = new GeometryFactory();
static List<Geometry> squares(int n, double step, double offset) {
List<Geometry> out = new ArrayList<>();
int side = (int) Math.ceil(Math.sqrt(n));
for (int i = 0; i < n; i++) {
double x = (i % side) * step + offset;
double y = (i / side) * step + offset;
out.add(GF.toGeometry(new Envelope(x, x + step * 1.4, y, y + step * 1.4)));
}
return out;
}
public static void main(String[] args) {
List<Geometry> geoms = squares(400, 10, 0);
time("UnaryUnionOp ", () -> UnaryUnionOp.union(geoms));
time("CascadedPolygonUnion", () -> CascadedPolygonUnion.union(geoms));
time("OverlayNGRobust ", () -> OverlayNGRobust.union(geoms));
time("UnaryUnionNG pm=1 ", () -> UnaryUnionNG.union(geoms, new PrecisionModel(1.0)));
}
static void time(String label, java.util.function.Supplier<Geometry> op) {
op.get(); // 预热,避免把 JIT 时间算进去
long t0 = System.nanoTime();
Geometry g = op.get();
long ms = (System.nanoTime() - t0) / 1_000_000;
System.out.printf("%s => area=%.0f num=%d %d ms%n",
label, g.getArea(), g.getNumGeometries(), ms);
}
}
一台普通笔记本上的实测(耗时为第二次运行,仅供量级参考):
UnaryUnionOp => area=78400 num=1 14 ms
CascadedPolygonUnion => area=78400 num=1 16 ms
OverlayNGRobust => area=78400 num=1 19 ms
UnaryUnionNG pm=1 => area=78400 num=1 14 ms
结论与直觉略有出入:在这个规模下四个入口相差不大,新引擎甚至略慢(因为它多了校验与重试框架的开销)。级联策略的收益要到大几千个面才明显。因此选型时不要“唯新是从”,而应按数据规模与脏净程度决定:小数据用默认入口,大数据优先试 CascadedPolygonUnion 与 UnaryUnionNG,脏数据用 OverlayNGRobust。
补充:1.18 的版本历史里有一条 “Improve performance of UnaryUnionOp by removing OverlayUnion optimization (#644)”,说明早期
UnaryUnionOp里曾有一条更激进的优化路径被移除,因为带来的边界情况比收益多。读老版本源码时不要被那部分代码误导。
9.10 CoverageUnion:覆盖面专用溶解
9.10.1 前提条件与收益
CoverageUnion(1.20 新增)针对的是“多边形覆盖”(polygonal coverage)这种特殊数据:一组多边形,彼此不重叠、不留缝隙,拼起来正好铺满一片区域。行政区划、地籍图、土地利用图都属于这一类。
对这类数据,溶解内部边其实不需要完整的叠加计算——只要把所有多边形的边界边收集起来,去掉那些“被两个多边形共享”的边,剩下的边围成的面就是并集。这正是 CoverageUnion 的做法:它用一个 BoundaryChainNoder 把边界链切开,按边的共享关系去重,再用 OverlayNG 的严格模式拼面。省掉了 noding 求交与标签推导两大开销,因此比通用并集快得多。
它的入口非常简洁,只有一个条带参数的静态方法:
// 签名:public static Geometry union(Geometry coverage)
Geometry dissolved = CoverageUnion.union(coverageMultiPolygon);
9.10.2 用错前提的后果
致命的一点是:CoverageUnion 不校验输入是否为合法覆盖。如果传进去的多边形互相重叠,它不会报错,而是按“边的共享次数”这个纯拓扑规则硬算,结果可能把重叠区域重复计入。下面这段代码就是反例:
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.overlayng.CoverageUnion;
public class CoverageUnionPitfall {
public static void main(String[] args) throws Exception {
WKTReader reader = new WKTReader();
// 两个重叠的多边形:不是合法覆盖
Geometry badCoverage = reader.read(
"MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0)), ((5 5, 15 5, 15 15, 5 15, 5 5)))");
Geometry result = CoverageUnion.union(badCoverage);
System.out.println("result = " + result);
System.out.println("num = " + result.getNumGeometries());
System.out.println("area = " + result.getArea());
System.out.println("expected union area = 175");
}
}
输出(注意这是错误结果,但没有任何异常):
result = MULTIPOLYGON (((0 0, 0 10, 10 10, 10 0, 0 0)), ((5 5, 5 15, 15 15, 15 5, 5 5)))
num = 2
area = 200.0
expected union area = 175
两个多边形只是被原样保留了下来,面积 200 而不是正确的 175(重叠的 25 被重复计算)。所以工程上必须遵守一条硬规则:调用 CoverageUnion 之前,先用 CoverageValidator.isValid(coverage) 校验。校验通过再溶解;校验失败则改用 UnaryUnionOp(慢但正确),或先走第 17 章的清洗流程。
9.10.3 与通用并集的对比
| 维度 | UnaryUnionOp.union |
CoverageUnion.union |
|---|---|---|
| 输入前提 | 任意(可重叠、可自相交) | 必须是合法覆盖(无重叠、无缝隙) |
| 是否做 noding | 是 | 否(用边界链去重替代) |
| 是否校验输入 | 否 | 否(必须调用方自己校验) |
| 结果正确性 | 对任意输入正确 | 仅对合法覆盖正确 |
| 典型耗时 | 1 倍 | 约 0.1 到 0.3 倍 |
| 适用数据 | 任意多边形集合 | 行政区、地籍、土地利用覆盖 |
9.10.4 组合使用的最佳实践
覆盖面的完整处理流程通常是:先用 CoverageValidator.validate 找出问题,用 CoverageCleaner(1.21 新增)或 GeometryFixer 修掉,再用 CoverageUnion 溶解,最后对结果做一次 PrecisionReducer 统一精度。第 17 章会展开质检细节,这里给一段骨架代码:
import java.util.Arrays;
import org.locationtech.jts.coverage.CoverageCleaner;
import org.locationtech.jts.coverage.CoverageValidator;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.MultiPolygon;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.operation.overlayng.CoverageUnion;
public class CoveragePipeline {
public static void main(String[] args) {
Geometry[] polygons = loadCoverage();
if (!CoverageValidator.isValid(polygons)) {
System.out.println("覆盖不合法,先清洗");
// 第二个参数为负表示由实现自动推导吸附容差,第三个参数为可填补的最大缝隙宽度
polygons = CoverageCleaner.clean(polygons, -1, 1.5);
}
MultiPolygon coverage = polygons[0].getFactory()
.createMultiPolygon(Arrays.stream(polygons)
.map(g -> (Polygon) g)
.toArray(Polygon[]::new));
Geometry dissolved = CoverageUnion.union(coverage);
System.out.println("dissolved area = " + dissolved.getArea());
System.out.println("dissolved num = " + dissolved.getNumGeometries());
}
static Geometry[] loadCoverage() {
throw new UnsupportedOperationException("替换为真实读取逻辑");
}
}
9.11 TopologyException 的成因与容差策略
9.11.1 三类异常消息
叠加分析抛出的异常几乎全部是 TopologyException(继承自 RuntimeException),消息可以归成三类,每一类指向不同的根因:
| 异常消息(节选) | 抛出位置 | 根因 | 处理方向 |
|---|---|---|---|
found non-noded intersection between ... |
旧引擎构环阶段 | noding 未把应相交的边打断 | 提高精度、吸附、改固定精度 |
side location conflict ... |
新旧引擎的标签推导 | 同一条边两侧的 Location 矛盾,通常是输入非法或有重叠 | 先修复输入几何 |
Result area inconsistent with overlay operation |
OverlayNG 面积启发式 |
浮点 noding 使图翻转 | 改用固定精度或 OverlayNGRobust |
第一类最常见,也最容易误判成“JTS 的 bug”。实际上它几乎总是“输入本身不可靠”或“浮点精度不足”导致的。第二类更严重:它通常意味着输入几何无效——例如两个多边形互相重叠(自相交的多边形集合本身就违反“面内部不相交”的定义),这时无论怎么重试都不会得到正确结果。第三类是保护性检查,它宁可抛异常也不返回错误几何,遇到时直接降精度重算即可。
9.11.2 一次完整的“踩坑到解决”过程
下面这段代码构造了一对“大量共享边且互相重叠”的方块集合,它们在旧引擎与浮点新引擎下都会抛出 side location conflict。这不是构造的假想,而是本章基准程序里真实发生的一次失败:
import java.util.ArrayList;
import java.util.List;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.operation.overlay.OverlayOp;
import org.locationtech.jts.operation.overlay.snap.SnapIfNeededOverlayOp;
import org.locationtech.jts.operation.overlayng.OverlayNG;
public class OverlayFailureDemo {
static final GeometryFactory GF = new GeometryFactory();
static Geometry grid(int n, double step, double size, double offset) {
List<Geometry> parts = new ArrayList<>();
int side = (int) Math.ceil(Math.sqrt(n));
for (int i = 0; i < n; i++) {
double x = (i % side) * step + offset;
double y = (i / side) * step + offset;
parts.add(GF.toGeometry(new Envelope(x, x + size, y, y + size)));
}
return GF.buildGeometry(parts);
}
public static void main(String[] args) {
Geometry a = grid(400, 10, 14, 0);
Geometry b = grid(400, 10, 14, 5);
attempt("old OverlayOp ", () -> OverlayOp.overlayOp(a, b, OverlayOp.INTERSECTION));
attempt("SnapIfNeeded ", () -> SnapIfNeededOverlayOp.intersection(a, b));
attempt("OverlayNG float ", () -> OverlayNG.overlay(a, b, OverlayNG.INTERSECTION));
attempt("OverlayNG fixed1", () ->
OverlayNG.overlay(a, b, OverlayNG.INTERSECTION, new PrecisionModel(1.0)));
}
static void attempt(String label, java.util.function.Supplier<Geometry> op) {
try {
Geometry g = op.get();
System.out.println(label + " => OK area=" + g.getArea()
+ " num=" + g.getNumGeometries());
} catch (RuntimeException ex) {
System.out.println(label + " => FAILED " + ex.getMessage());
}
}
}
输出(本机实测):
old OverlayOp => FAILED side location conflict [ (0.0, 10.0, NaN) ]
SnapIfNeeded => FAILED side location conflict [ (0.0, 10.0, NaN) ]
OverlayNG float => FAILED side location conflict: arg 0 [ (50.0, 114.0, NaN) ]
OverlayNG fixed1 => OK area=39601.0 num=1
只有固定精度活了下来。原因是这些方块大量共享边、且互相重叠,浮点 noding 会在共享边的端点处产生方向与位置矛盾,而 Snap-Rounding 把所有顶点先吸到整数网格上,让共享边真正重合,矛盾随之消失。这是前文“固定精度更稳”的最直接证据。
9.11.3 容差策略的四条经验
结合源码里的容差计算与上面的实验,可以总结出四条可操作的经验:
第一条,先判断是数据问题还是算法问题。用 geom.isValid() 与 IsValidOp 检查输入;如果输入本身重叠或自相交,先去修数据(第 17 章),不要指望换引擎或加容差。
第二条,对“本应重合但不完全重合”的边,主动吸附。可以用 GeometrySnapper.snapToSelf(geom, tol, true) 对每个输入自吸附,容差取业务精度的 1 到 10 倍;也可以直接用固定精度叠加,让引擎替你吸。
第三条,容差要与坐标量级匹配。GeometrySnapper 的默认容差是包络短边乘 1e-9,OverlayNGRobust 的初始吸附容差是最大 ordinate 乘 1e-12;这两个值在 1e2 量级的数据上分别是 1e-7 与 1e-10,都很小。如果数据是经纬度(量级 1e2),1e-10 度约等于 1e-5 米,勉强够用;如果是投影坐标(量级 1e7 米),1e-12 乘 1e7 等于 1e-5 米,也够用。但如果业务精度是厘米级而输入噪声是米级,这些容差完全不够,必须显式指定固定精度。
第四条,捕获异常后要有降级路径。生产代码建议按“OverlayNGRobust → 固定精度重算 → 记录并跳过”三级降级,而不是让一个要素的拓扑异常拖垮整个批处理:
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.operation.overlayng.OverlayNG;
import org.locationtech.jts.operation.overlayng.OverlayNGRobust;
import org.locationtech.jts.operation.overlayng.PrecisionUtil;
public class OverlayWithFallback {
public static Geometry safeIntersection(Geometry a, Geometry b) {
try {
return OverlayNGRobust.overlay(a, b, OverlayNG.INTERSECTION);
} catch (RuntimeException first) {
double scale = PrecisionUtil.safeScale(a, b);
try {
return OverlayNG.overlay(a, b, OverlayNG.INTERSECTION,
new PrecisionModel(scale));
} catch (RuntimeException second) {
throw new IllegalStateException(
"叠加失败: " + first.getMessage() + " / " + second.getMessage(), first);
}
}
}
}
9.12 新旧引擎的结果与性能对比
9.12.1 结果层面的差异
对合法、无近重合边的输入,两套引擎的结果在拓扑上等价,差异只体现在三处:顶点坐标(新引擎在浮点模式下通常保留更多有效位)、环的起点与方向、MultiPolygon 内部元素顺序。用 norm() 规范化后再比较,通常可以判定为“完全相同”。
对“接缝有微小偏差”的输入,两者会给出结构不同的结果。9.7.4 的例子已经展示过:旧引擎把 1e-9 的缝隙当成真实缝隙,返回两个多边形;新引擎把它当作重合边,返回一个多边形。这不是谁对谁错,而是精度假设不同——旧引擎的吸附只在抛异常之后才发生(SnapIfNeededOverlayOp),而新引擎在固定精度下的 Snap-Rounding 是默认行为。
对“病态”输入,9.11.2 给出了更极端的情况:旧引擎与新引擎的浮点模式都失败,只有固定精度成功。
| 输入类型 | 旧引擎(默认) | OverlayNG 浮点 | OverlayNG 固定精度 |
|---|---|---|---|
| 合法、无近重合边 | 正确 | 正确 | 正确(坐标舍入) |
| 接缝偏差 1e-9 | 保守:视为分离 | 视为重合 | 视为重合 |
| 大量共享边且重叠 | 抛 side location conflict |
抛 side location conflict |
正确 |
| 自相交蝴蝶结 | 抛 non-noded intersection |
抛异常 | 抛异常(输入问题) |
| 含 Z 值 | 不插值 | 线性插值 | 线性插值 |
9.12.2 性能层面的差异
性能是选型的另一个维度。9.9.4 的数据显示:在小规模、干净的输入上,旧引擎往往更快,因为它的数据结构更紧凑、没有校验与重试框架;在需要大量 noding 的大输入上,新引擎的裁剪优化与更轻的边结构开始占优。
用锯齿多边形做基准(两个各含 8005 个顶点的“梳齿”多边形,互相错位 0.5 重叠)可以得到一组有代表性的对比。构造方法如下:
static Geometry comb(int teeth, double toothW, double toothH, double offset) {
List<Coordinate> cs = new ArrayList<>();
cs.add(new Coordinate(0, 0));
for (int i = 0; i < teeth; i++) {
double x0 = offset + i * toothW * 2;
cs.add(new Coordinate(x0, 0));
cs.add(new Coordinate(x0, toothH));
cs.add(new Coordinate(x0 + toothW, toothH));
cs.add(new Coordinate(x0 + toothW, 0));
}
cs.add(new Coordinate(offset + teeth * toothW * 2, 0));
cs.add(new Coordinate(offset + teeth * toothW * 2, -5));
cs.add(new Coordinate(0, -5));
cs.add(new Coordinate(0, 0));
return new GeometryFactory().createPolygon(cs.toArray(new Coordinate[0]));
}
本机实测(每个操作先预热一次,取第二次耗时;绝对值随机器变化,只比较同一台机器上的相对关系):
| 运算 | 旧引擎 | OverlayNG 浮点 | OverlayNG 固定精度 | OverlayNGRobust |
|---|---|---|---|---|
| 交集 | 106 ms | 76 ms | 488 ms | 36 ms |
| 并集 | 54 ms | 37 ms | 未测 | 未测 |
| 差集 | 54 ms | 34 ms | 未测 | 未测 |
三点解读。第一,浮点新引擎在中大型输入上比旧引擎快 30% 到 40%,收益来自更轻的边结构与裁剪优化。第二,固定精度慢得多(交集 488 ms),因为 Snap-Rounding 要把每个交点吸附到网格并重建拓扑;但它换来的是“绝不失败”,这笔账是否划算取决于业务。第三,OverlayNGRobust 的交集只要 36 ms,比纯浮点路径更快——因为它内部第一级就用浮点算成功了,而它的调用路径上少了一层旧引擎的平面图结构开销。
重要提醒:固定精度下若 scale 选得过大(网格过粗),薄特征会塌缩。上例中 0.5 宽的梳齿在 scale=1.0(整米网格)下全部被吸成零面积,交集结果从 1 个面积 70000 的多边形变成 4000 个面积 20000 的碎片。选 scale 时必须保证“最窄业务特征远大于 1/scale”。
9.12.3 决策清单
| 你的情况 | 建议 |
| — | — |
| 数据来历不明、批量处理、不能中断 | OverlayNGRobust.overlay 加固定精度兜底 |
| 需要结果可复现、可入库比对 | 固定精度 OverlayNG,scale 取业务精度的 10 倍 |
| 单次小规模计算、追求低延迟 | 默认 Geometry.union(旧引擎)即可 |
| 合法覆盖面批量溶解 | CoverageUnion(务必先跑 CoverageValidator) |
| 大量多边形溶解 | CascadedPolygonUnion 或 UnaryUnionNG |
| 需要混合维度结果做诊断 | OverlayNG 非严格模式 |
| 下游算子要求同质输入 | OverlayNG 加 setStrictMode(true) |
9.13 常见问题与避坑
问题 1:两个多边形明明相交,intersection 却返回 POLYGON EMPTY。
先确认两个几何的包络是否相交(getEnvelopeInternal().intersects(other))。若包络相交但结果为空,多半是坐标系不一致(一个经纬度、一个投影坐标),导致两者在数值上相距极远。JTS 不含 CRS 概念,这类错误不会被引擎发现,必须由上层校验。
问题 2:结果面积比两个输入之和还大。
典型的重复计算,常见于三种情况:对互相重叠的多边形用 CoverageUnion(见 9.10.2);对同一个几何做了两次并集且中间结果未被正确去重;用了 GeometryCollection 输入而交集被逐元素映射后简单拼接。解决方法是改用 UnaryUnionOp.union,并对结果做 isValid() 检查。
问题 3:union 报 Operation does not support GeometryCollection arguments。
二元并、差、对称差不接受非空 GeometryCollection。按维度拆分输入(PolygonExtracter、LineStringExtracter、PointExtracter),或改用一元并集 geom.union()。
问题 4:同样的代码,本地结果和服务器结果不一致。
两种可能:一是全局属性不同(服务器上设了 -Djts.overlay=ng),二是浮点运算在不同 JIT 与 CPU 上的差异。要彻底消除,必须用固定精度;只靠浮点无法保证跨平台逐位一致。
问题 5:TopologyException 捕获后重试仍然失败。
先看消息类型。side location conflict 与 non-noded intersection 若在固定精度下依旧出现,说明输入几何本身无效(自相交、环方向错误、坐标含 NaN),应先跑 GeometryFixer.fix。IsValidOp 能定位到具体的问题点,比盲目重试高效得多。
问题 6:叠加之后顶点暴增,内存吃紧。
交集与差集在裁剪优化下只处理相关区域,但并集必须处理全部边。若结果顶点过多,可对结果做 PrecisionReducer.reducePrecision 降精度,或用第 12 章的 DouglasPeuckerSimplifier 化简(注意化简必须在叠加之后,且要考虑拓扑保持)。
问题 7:Z 值在结果里丢了。
WKTWriter 默认输出二维,看起来像丢了;实际上 OverlayNG 的 ElevationModel 会沿边线性插值 Z。确认方式是用 new WKTWriter(3) 输出,或在代码里读 getCoordinate().getZ()。旧引擎不插值,只保留来源坐标,所以来回切换引擎会看到 Z 值差异。
问题 8:difference 之后公共边消失了。
这是 OGC 语义决定的,不是 bug(见 9.1.4)。需要保留界址线时,应在线图层单独求交:boundaryA.intersection(boundaryB)。
问题 9:symDifference 的结果总是 MultiPolygon,无法直接入库。
对称差天然可能产生多个不相连的块,返回 MultiPolygon 是正确行为。若业务表只接受 Polygon,需要按面积排序取主块,或改用两次 difference 分别处理。
问题 10:大量小多边形求并集时 OOM。
改用 CascadedPolygonUnion 分批处理;或先把输入按网格分块(STRtree 查询),块内溶解后再块间溶解。不要一次性把上百万个多边形塞进 UnaryUnionOp。
9.14 本章小结
本章从数学语义一路讲到源码实现与工程决策,要点可以浓缩成十条:
- 四个叠加运算在数学上是集合的交、并、差、对称差,JTS 用
Geometry的四个方法对外暴露,返回值遵循“最小充分类型”,可能是Polygon、MultiPolygon或GeometryCollection。 - 结果维度可能低于输入维度(拓扑塌缩),这是语义的一部分而不是异常;需要同质结果时用
OverlayNG.setStrictMode(true)。 - 代数恒等式在浮点坐标下只是近似成立;断言结果时应比较面积、周长、顶点集合等不变量,而不是 WKT 字符串。
- 空几何与
GeometryCollection有专门的短路与限制规则:交集支持 GC(逐元素映射),二元并、差、对称差不支持。 - JTS 有两套引擎:旧
OverlayOp与新OverlayNG;高层方法默认仍走旧引擎,开关是系统属性jts.overlay,取值为old或ng。 OverlayNG的流水线是“noding → 建图 → 标注 → 按维度构结果 → 插值 Z”,其中 noding 的精度模型决定鲁棒性,标注阶段决定结果语义。OverlayLabel用边界边、塌缩边、线边、非本几何边四种状态描述一条边的身份,这个设计使“重合边合并”与“塌缩边丢弃”成为统一规则的特例。- 精度是可控的:
PrecisionUtil给出安全比例尺与固有比例尺,PrecisionReducer负责降精度且保证输出合法;固定精度更稳但更慢,且可能让薄特征塌缩。 OverlayNGRobust用“浮点 → 五轮递增容差吸附 → Snap-Rounding”三级阶梯兜底,仍失败则抛出最初的异常;它救不了非法输入与混合维度输入。- 一元并集是一族工具而非一个方法:
UnaryUnionOp通用、CascadedPolygonUnion适合大批量、CoverageUnion只对合法覆盖正确且最快;用错前提会静默给出错误结果。
下一章进入缓冲区与偏移曲线:叠加回答“交集在哪”,缓冲区回答“距离它多远的范围内”。两者在实现上共享 noding 与构面技术,但在参数化与误差控制上有完全不同的关注点。
相关资源
- JTS 官方 Javadoc:https://locationtech.github.io/jts/javadoc/ ,重点读
OverlayNG、OverlayNGRobust、PrecisionUtil、CoverageUnion的类注释,本章很多结论直接从这些注释与源码得来。 - 源码位置:
modules/core/src/main/java/org/locationtech/jts/operation/overlayng/(新引擎)、同目录下的operation/overlay/(旧引擎)、operation/union/(一元并集)。 - 版本历史:
doc/JTS_Version_History.md。OverlayNG 的引入见 1.18.0 条目(#599、#615、#645);CoverageUnion 见 1.20.0 条目(#900);OverlayNG 面积启发式见 1.18.1 条目(#812)。 - 本站第 08 章《空间关系谓词与 DE-9IM》:https://znlgis.github.io/gis/tutorial/jts/第08章-空间关系谓词与DE-9IM/
- 本站第 10 章《缓冲区与偏移曲线》:https://znlgis.github.io/gis/tutorial/jts/第10章-缓冲区与偏移曲线/
- 本站第 17 章《有效性与 Coverage 质检》:https://znlgis.github.io/gis/tutorial/jts/第17章-有效性与Coverage质检/
- 本站第 18 章《精度模型与数值鲁棒性》:https://znlgis.github.io/gis/tutorial/jts/第18章-精度模型与数值鲁棒性/