JTS空间几何实战：关系判断、缓冲区与STRtree索引避坑

1. 先划清边界JTS 算的到底是什么几何接手过一个园区资产管理的数据核对需求业务方的诉求很朴素给一批设备坐标判断每台设备落在哪个管理区块里顺带算一下设备到区块边界的最近距离。数据库里坐标是经纬度我图省事直接用JTS把点的 WKT 和区块的 WKT 拼出来做contains判断跑完发现 80% 的设备都被判定为不在任何区块内距离算出来是零点几。当时第一反应是数据脏查了半天坐标最后发现问题在我自己身上——JTS 是个纯平面欧氏几何库它不知道什么是经度纬度也不知道什么是米。你把经纬度直接丢进去它就当成长度单位是度的笛卡尔坐标来算那个0.0003的距离单位是度不是你想要的米。这件事是我用 JTS 这些年踩过最大的一个坑也是我写这篇东西的直接动机。JTSJava Topology Suite是一套面向 Java 的二维空间几何算法库核心能力可以拆成四块几何对象的构建与序列化、空间关系谓词判断、空间叠加与缓冲区运算、空间索引与批量查询。它支撑的场景非常广——地图服务里的面面相交判断、轨迹数据处理中的缓冲区碰撞检测、后台的空间数据清洗、栅格化前的几何预处理甚至一些非 GIS 领域的排布碰撞检测都能用。适合谁看我觉得分两类一类是刚接触空间数据的后端开发被TopologyException和坐标系问题折磨得怀疑人生另一类是已经用了几年但对精度、有效性、性能这些底层机制一直模模糊糊想系统理一遍的同行。2. 几何对象构建从 WKT 到能参与运算的实体2.1 依赖引入与包名这个历史包袱先说个新手容易卡半天的点。JTS 经历过一次坐标迁移老版本在com.vividsolutions.jts这个包下1.15 之后迁到了org.locationtech.jts。你在网上搜到的示例代码一半是老包名直接抄会编译不过。现在正常的做法是引jts-core版本用 1.18 或 1.19 都可以1.19 在几何修复和索引查询 API 上做了不少易用性改进我个人更推荐。dependency groupIdorg.locationtech.jts/groupId artifactIdjts-core/artifactId version1.19.0/version /dependency依赖本身很轻没有传递依赖这点在服务端项目里挺友好。如果你还要读写 GeoJSON 或者做坐标投影转换那是另外的库的事JTS 自己不干这个别指望它。2.2 六种基本类型与它们的维度约定JTS 的几何类型族谱不复杂理解清楚继承关系能省很多instanceof的功夫。Geometry是抽象基类底下分Point维度 0、LineString维度 1、Polygon维度 2以及对应的多部件版本MultiPoint、MultiLineString、MultiPolygon再加一个万能容器GeometryCollection。每种几何都有一个getDimension()返回值就是 0、1、2写通用逻辑时用它做分支比判断具体类型稳。有几个细节值得单独拎出来。Polygon由外壳ExteriorRing和零到多个内洞InteriorRing组成构造函数传的是LinearRing数组第一个元素必须是外壳顺序错了直接抛异常。LinearRing是一个首尾坐标必须闭合的LineString构造时如果没闭合GeometryFactory会自动帮你补一个首点但如果你手写坐标数组又自己拼错了顺序后面算面积可能出负数。至于MultiPolygon它要求里面的每个Polygon互不相交不能有共享边这个约束很多人不知道拼错了isValid()会返回 false。2.3 GeometryFactory 与 PrecisionModel精度先定再谈几何GeometryFactory是造几何的入口几乎所有构建方法都挂在它上面。它有两个可以定制的构造参数一个是PrecisionModel一个是 SRID。这两个东西决定了后面所有运算的地基我建议在任何项目里都显式指定不要用默认的new GeometryFactory()图省事。PrecisionModel有三种模式模式构造方式适用场景副作用FLOATINGnew PrecisionModel()默认不做取舍保留 double 全精度浮点误差会累积长链条运算后可能出现极小缝隙FLOATING_SINGLEnew PrecisionModel(PrecisionModel.FLOATING_SINGLE)数据源本身是 float或需要更强的容差精度损失明显坐标量级大时误差可观FIXEDnew PrecisionModel(1000)固定小数位比如毫米级、0.001 度级坐标会被强制 snap 到网格可能引入新的拓扑问题选哪个我的经验是数据本身干净、运算链条短用默认 FLOATING 就行如果要做大量叠加运算、结果还要入库比对用 FIXED 配合一个合理的网格尺度能显著减少差一点点就对不上的问题。那个网格尺度怎么定看你数据的实际精度比如业务上坐标只到 0.000001 度约 0.1 米那就别把模型设成 1e9 那么细设成 1e6 反而稳。至于 SRID这里必须强调一句JTS 里的 SRID 只是一个整型标签不参与任何计算也不触发任何坐标转换。你给一个经纬度的几何标上 4326再给一个投影坐标的几何标上 3857两者做intersectsJTS 照样直接拿数字算结果是垃圾但不会报错。这个静默失败是最阴险的因为没有任何异常提示。PrecisionModel pm new PrecisionModel(1000000); // 1e-6 网格 GeometryFactory gf new GeometryFactory(pm, 4326); WKTReader reader new WKTReader(gf); Polygon block (Polygon) reader.read( POLYGON ((116.30 39.90, 116.32 39.90, 116.32 39.92, 116.30 39.92, 116.30 39.90)));3. 空间关系判断九种谓词里藏着的边界语义3.1 DE-9IM所有关系判断背后的统一模型JTS 的关系判断不是各写各的它们全部建立在DE-9IM这个模型上。简单说任意两个几何 A 和 B各自有内部Interior、边界Boundary、外部Exterior三个部分两两组合成 3×3 共九个格子每个格子记录这两个部分的交集维度-1 表示不相交0/1/2 表示点/线/面。relate(A, B)返回的就是这九个数字拼成的字符串比如212101212。intersects、contains、within、touches、crosses、overlaps、covers、coveredBy、equals、disjoint这些谓词本质上都是对 DE-9IM 矩阵的模式匹配。理解这一点你在遇到为什么这个点明明在面上contains却返回 false的时候就知道该去查矩阵而不是怀疑数据。调试时有个很实用的小技巧把relate的结果打出来看。比如判断点在不在面内point.relate(polygon)打出0FFFFF212说明点在面内部打出F0FFFF212说明点落在外壳线上也就是边界上。3.2 contains 和 covers 的差别害过太多人这是我最想强调的一节。contains的定义是B 的内部与 A 的内部有交集且 B 的内部不与 A 的外部相交且 B 不完全落在 A 的边界上。最后这半句是关键——如果一个点正好落在多边形的边线上polygon.contains(point)返回false。covers则宽松一些只要 B 的任意部分含边界和 A 的内部或边界有交集且 B 不在 A 外部就返回 true。所以点落在边线上时polygon.covers(point)返回true。实战里怎么选设备点位判断、围栏判断这类需求我建议默认用covers或者intersects因为业务方认为在边上就是在里面。只有当业务明确要求必须在内部、不能压线时才用contains。这个坑我见过的最惨案例是某次统计报表所有恰好落在区块边界线上因为坐标被取整过的设备全部被归到区块外报表数和现场对不上排查了两个礼拜。// 落在边线上的点 Point p gf.createPoint(new Coordinate(116.30, 39.91)); boolean c1 block.contains(p); // false —— 压线不算包含 boolean c2 block.covers(p); // true —— 压线算覆盖 boolean c3 block.intersects(p); // true —— 只要碰到就算3.3 关系判断的性能与批量场景的写法单个intersects的开销其实很小真正拖慢系统的是两两比较。一万个点对一万个面暴力双循环是一个亿次判断几秒钟跑不完数据量再上一个量级直接不可用。这时候的思路是先做包络矩形Envelope预筛把明显不相交的候选砍掉再对剩下的做精确判断。JTS 提供了STRtree就是干这个的后面单独一节说。另外提一个常被忽略的点Geometry对象不是线程安全的JTS 官方文档里写得很含蓄但实操中你如果在多线程里共享同一个Polygon实例去做运算偶发的诡异结果会很难复现。稳妥做法是每个线程各自构建或者加同步。还有JTS 的isValid()和getArea()这类方法内部有缓存字段虽然没有明显的并发问题但为了省心我还是倾向于不共享。4. 叠加运算与缓冲区参数选不对结果全白费4.1 buffer 的四个参数怎么定buffer是用得最多的运算之一做法是给几何向外扩一个距离生成新的几何。它最少要传一个distance但真正决定结果质量的是另外几个参数很多人根本不知道它们存在。quadrantSegments默认 8。它控制圆弧被切分成多少段每四分之一圆分几段。段数越多越接近真圆但顶点数暴涨。有个估算公式近似圆与真圆的最大径向误差约等于r × (1 - cos(π / (2 × quadrantSegments)))。代入默认值 8误差大约是半径的 0.48%。如果你缓冲半径是 100 米误差接近半米。要控到 1 厘米以内需要quadrantSegments到 40 左右。别盲目调大顶点数会跟着涨十几倍后续相交运算的开销也要涨。endCapStyle端点样式。CAP_ROUND圆头默认CAP_FLAT平头CAP_SQUARE方头。做路径覆盖范围时平头和方头能省掉大量顶点视觉上也更符合路段的业务预期。joinStyle拐角样式。JOIN_ROUND圆角默认JOIN_MITRE尖角JOIN_BEVEL斜切。尖角在某些锐角拐弯处会向外戳得很远mitre 长度有限制JTS 内部有mitreLimit控制做建筑轮廓时通常用尖角更贴合实际。simplifyFactor对结果做抽稀的容差默认不减。BufferParameters bp new BufferParameters(); bp.setQuadrantSegments(16); bp.setEndCapStyle(BufferParameters.CAP_FLAT); bp.setJoinStyle(BufferParameters.JOIN_MITRE); bp.setSimplifyFactor(0.01); Geometry buffered geom.buffer(50, bp);distance的正负也有讲究正数是外扩负数是内缩。对多边形内缩时如果距离超过多边形最窄处的半宽结果会变成空几何isEmpty()返回 true这不是 bug是数学上真实如此。做退界计算时务必先判断isEmpty()否则后面拿空几何去做相交会得到一堆空结果排查起来很费劲。4.2 叠加分析union、intersection、difference 的取舍叠加运算有四兄弟union并、intersection交、difference差、symDifference对称差。它们的语义很直白坑主要在输入数据质量上。最常见的异常是TopologyException: found non-noded intersection between LINESTRING ... and LINESTRING ...。这个报错的意思是运算过程中发现两条线在某个位置相交了但那个交点没有被显式地记录成一个顶点拓扑上有歧义。触发原因通常是输入几何本身无效自相交、重复点或者浮点误差导致理论上重合的边实际差了几微米。排查手段有两步第一步对两个输入都跑isValid()拿到 false 的那个用explain()看具体原因第二步如果输入都有效试着对参与运算的几何做一次GeometryPrecisionReducer把坐标统一到某个网格上很多时候浮点噪音吃掉就过了。还有一种粗暴但有效的办法对输入先buffer(0)这个操作本身就是在做一次清理自相交的拓扑重建。4.3 简化与抽稀两种 simplifier 不能随便换轨迹抽稀、轮廓压缩这类需求要用简化器。JTS 里有两个主要实现名字很像但行为差很多简化器策略优点缺点DouglasPeuckerSimplifier按点到弦的垂直距离做取舍快压缩率高形状保真度好可能产生自相交多边形可能变无效TopologyPreservingSimplifier在 DP 基础上加拓扑约束结果保证有效不会自相交慢一些压缩率略低我的建议很明确只要你的几何是多边形或者需要保持有效性一律用 TopologyPreservingSimplifier。用 DP 简化多边形然后拿去做面积统计你会得到一堆无效几何和错误面积。有个折中做法是先简化再用GeometryFixer.fix()修一遍但那就是绕远路了。5. 空间索引与批量处理STRtree 的正确打开方式5.1 构建与查询的基本套路回到前面点面判断的场景一万个点对一万个面。正确姿势是给面建一个STRtree然后每个点拿自己的包络矩形去查只对命中的候选做精确判断。STRtree index new STRtree(); for (Polygon block : blocks) { index.insert(block.getEnvelopeInternal(), block); } index.build(); // 显式构建之后查询才快 for (Point p : points) { List? candidates index.query(p.getEnvelopeInternal()); for (Object o : candidates) { Polygon block (Polygon) o; if (block.covers(p)) { // 命中记录对应关系 break; } } }这里有几个细节。insert的第一个参数是Envelope第二个是你要存的任意对象索引本身不关心它是不是几何你甚至可以存 ID 字符串。build()是可选的不调用的话第一次查询时会自动构建但显式调用能让构建时机可控也避免多线程环境下第一次查询时的竞态。query返回的是包络矩形相交的候选集不是精确相交的结果。拿到候选后必须自己做精确判断这一步千万不能省——省了就会得到一堆误判。我见过有人直接拿query的结果当最终答案因为测试数据里所有几何的包络恰好都不重叠上线后数据一乱套结果全错。5.2 一对一查找和批量查找的写法差异如果需求是A 集合里每个元素找到 B 集合里相交的那个除了逐个query还有一个更省事的思路把 A 和 B 合并成一个GeometryCollection用STRtree分别建两个索引或者干脆用 JTS 的OverlayNG加上索引参数做批量叠加效率能再提一截。不过这属于性能优化的深水区数据量没到十万以上其实没必要。真正需要提醒的是内存。STRtree会把所有几何的包络常驻在内存里一个Envelope对象大概几十字节十万个就是几兆问题不大。但如果你把完整的Geometry也一起挂在索引里通常会的几何本身的顶点数据才是内存大头。百万级顶点量的多边形集合堆内存直接吃掉几百兆服务端跑着跑着就 OOM。这时候要么分片处理要么先把几何做简化别硬扛。6. 常见问题与排查技巧实录6.1 异常和错误结果速查表这张表是我这些年攒下来的遇到问题先对号入座能省掉大量瞎猜时间现象最可能的原因排查手段处理方式距离算出来是零点几用了经纬度直接算平面距离检查坐标数值范围纬度应在 -90~90先做投影转换或改用球面距离公式自算contains返回 false 但视觉上点在面里点在边线上打relate矩阵确认改用coversTopologyException: found non-noded intersection输入无效或浮点误差对输入跑isValid()和explain()buffer(0)清理或统一精度模型isValid()返回 falseexplain说 Self-intersection多边形边自交定位到具体坐标GeometryFixer.fix()修复面积为负环方向反了外环顺时针检查isCCW调normalize()或重排坐标相交结果为空但两者明显重叠SRID 不一致一个投影一个经纬度打印两者 SRID 和坐标范围统一坐标系这是脏数据问题多线程下运算结果偶发异常共享了 Geometry 实例检查是否有静态复用每线程独立构建或用不可变副本简化后几何无效用了 DP 简化器换实现改用TopologyPreservingSimplifier点标注位置跑到多边形外面用了getCentroid()打标签凹多边形的重心可能在外部改用getInteriorPoint()getCentroid和getInteriorPoint这个坑值得单独说一句。centroid是几何的质心数学上正确但它可能落在多边形外部——一个 U 形多边形重心正好在凹口里。做地图标注、POI 落位这类需求必须用getInteriorPoint()它保证返回的点一定在多边形内部。这是个一行代码的改动但不知道的话会一直觉得是渲染库有问题。6.2 数值精度那些不报错但结果错的地方浮点误差在空间运算里是系统性存在的不是偶发。最典型的表现是两个理论上应该完全重合的面相交后边缘出现了一圈极细的碎面。这种碎面面积接近 0但如果你的逻辑是相交面积大于 0 就算重叠这些碎面会被全部判成重叠报表数字直接翻倍。处理办法有三个层次。第一层判断面积时加一个极小阈值比如if (area 1e-9)再过而不是 0。这个阈值要根据你的坐标单位定投影坐标下取 0.01 平方米合理经纬度下就没法直接用了。第二层在构建几何时就用 FIXED 精度模型把坐标 snap 到网格上从源头减少噪音。第三层对结果做一次GeometryPrecisionReducer.reduce()把运算产生的细微偏差量化掉。还有一类精度问题比较隐蔽坐标数值量级太大。如果坐标系单位是米原点又离你的数据区很远X 坐标可能是1.2e7这种量级double 在这个量级上做减法可用的有效位数会少好几位。这时候的稳妥做法是先把坐标整体平移到原点附近做运算算完再平移回去。这个技巧在大比例尺工程数据里很常见值不值得一提但救过命。7. WKT/WKB 序列化与跨系统交换7.1 WKT 适合调试WKB 适合传输WKTReader和WKTWriter是最常用的序列化组合优点是肉眼可读日志里打出来一眼就能看出几何形状对不对。但 WKT 文本冗长一个稍微复杂点的多边形字符串几十 KB高频接口里传这个不合适。跨进程传输或者存数据库用WKBWriter更合理二进制体积能小一个量级。WKBWriter writer new WKBWriter(2, ByteOrderValues.LITTLE_ENDIAN); byte[] bytes writer.write(geom); WKBReader reader new WKBReader(); Geometry restored reader.read(bytes);构造WKBWriter时第一个参数是输出维度传 2 只写 XY传 3 会带上 Z。这里有个坑如果你的几何里 Z 值是Coordinate.NULL_ORDINATE也就是 NaN写成 3 维会出问题读回来的时候 Z 可能变成垃圾值。所以除非确实需要高程一律用 2 维。另外WKTWriter有个setTrim(true)能去掉坐标后面那些没用的零和小数尾巴日志体积能小不少。还有个setPrecisionModel()可以直接在输出时做精度取舍比先reduce再输出方便。7.2 解析容错别让一行脏数据搞崩整个任务WKTReader.read()遇到格式不对的字符串会直接抛ParseException。批量处理时如果不用 try-catch 包住一千条数据里有一条格式错误整个批次全废。我的习惯是写一个包装方法把异常转成业务上的解析失败结果记录下来但继续往下走。还有一种情况是read成功了但得到的几何是空的POINT EMPTY或者维度维度不对这些不会抛异常需要主动判isEmpty()和getNumPoints() 0。别指望库帮你兜底它只负责语法正确。8. 一个完整的落地示例设备落区判断把前面的点串起来给一个能直接抄的完整流程。需求是读取区块 WKT 列表和设备坐标列表输出每台设备所属区块和到区块边界的最近距离。第一步坐标统一。如果原始数据是经纬度而距离需要以米为单位必须先转投影这一步 JTS 不管得用单独的投影库转完再进 JTS。第二步构建几何时统一用一个GeometryFactory和同一个PrecisionModel。第三步建STRtree索引。第四步逐点查询精确判断用covers。第五步距离用geom.distance(point)注意这个距离是点到几何边界的最短距离如果点在多边形内部返回的是到最近边的距离这个语义正好符合到边界距离的需求。double dist block.distance(p); // 到边界的最近距离 boolean inside block.covers(p); // 是否在区块内含边界 Point label block.getInteriorPoint(); // 标注点保证在内部整个流程里最容易出错的仍然是第一步和第二步。我现在的习惯是在数据进 JTS 之前先做一个前置校验抽样打印坐标范围确认 X 落在合理的经度区间、Y 落在合理的纬度区间或者投影坐标落在预期量级。这一步花两分钟能省掉后面两天的排查。最后分享一个我个人坚持的做法任何涉及 JTS 的服务日志里必须能打出原始 WKT 和运算结果 WKT。空间问题一旦出错没有原始几何基本没法复现尤其是那种只在特定数据组合下才触发的拓扑异常。日志多写两行排查时就是救命稻草。