Every so often I produce an image in the JTS TestBuilder which strikes me as worthy of capture. Here's one that seems pretty seasonal:
A (long) while ago I posted about "soon-to-be-released" JTS code for polygon triangulation using Ear Clipping. It turned out it was actually in the category of "never-to-be-released". However, later I worked with a student, Dan Tong, on a coding exercise sponsored by Facebook. We decided to tackle polygon triangulation. By the end of the project he implemented a functional algorithm, including handling holes (via the hole-joining technique described by Eberly). To make it release-ready the code needed performance and structural improvements. This has finally happened and the code is out!
The new codebase has been substantially rewritten to improve modularity and the API. A new data structure to represent triangulations has been introduced, called a Tri It is simpler, more memory efficient, and easier to work with than the QuadEdge data structure previously used in JTS. It provides an excellent basis for developing further algorithms based on triangulations (of which more below)
Most of the reworking involved implementing key performance optimizations, utilizing existing and new JTS spatial index structures:
Performance must be benchmarked against an existing implementation. Geometry maestro Vladimir Agafonkin has developed an Ear Clipping algorithm in Javascript called Earcut. It's claimed to be the fastest implementation in Javascript. Indeed the performance comparison shows Earcut to be over 4 times faster than some other Javascript triangulation implementations. Earcut has been ported to Java as earcut4j, which allows comparing JTS to it.
| Dataset | Size | JTS Time | Earcut4j Time |
|---|---|---|---|
| Lake Superior | 3,478 vertices / 28 holes | 18 ms | 13 ms |
| Canada | 10,522 vertices | 35 ms | 26 ms |
| Complex polygon | 44,051 vertices / 125 holes | 1.1 s | 0.64 s |
| North America with Lakes | 115,206 vertices / 1,660 holes | 13.3 s | 5.4 s |
The triangulation algorithm is exposed as the PolygonTriangulator class (deliberately named after the effect, rather than the algorithm producing it). Here's an example of the output:
The output of PolygonTriangulator is guaranteed to be a valid triangulation, but in the interests of performance it does not attempt to produce a high-quality one. Triangulation quality is measured in terms of minimizing the number of "skinny" triangles (or equivalently, maximizing the sum of the interior angles). Optimal quality is provided by Constrained Delaunay Triangulation (CDT). In addition to being optimal in terms of triangle shape, the CDT has the nice properties of being (essentially) unique, and of accurately representing the structure of the input polygon. It also removes the influence that Hole Joining holes has on the raw triangulation provided by Ear-Clipping.
I originally proposed that a CDT can be achieved by a Delaunay Improvement algorithm based on iterated triangle flipping. It turned out that this technique is not only effective, but also pleasingly performant. This is implemented in the ConstrainedDelaunayTriangulator class. Running it on the example above shows the improvement obtained. Note how there are fewer narrow triangles, and their arrangement more closely represents the structure of the polygon.
This code will appear in the next JTS version (1.18.3 or 1.19), and has already been released in GEOS 3.10.
It seems possible to adapt the machinery for Ear Clipping to provide two other useful constructions:
Recently a GEOS patch was contributed to change the KdTree query implementation to use an explicit stack rather than recursion. This has been ported to JTS as PR #779 (along with some refactoring).
The change was motivated by a QGIS issue in which a union of some large polygons caused a stack overflow during a KdTree query. The reason is that the poor vertex alignment of the input polygons causes the union overlay process to invoke the SnappingNoder. (For background about why this occurs see the post OverlayNG - Noding Strategies). The snapping noder snaps vertices using a KdTree with a tolerance. The polygon boundaries contain runs of coherent (nearly-monotonic) vertices (which is typical for high-resolution polygons delineating natural areas). When these are loaded directly into a KdTree the tree can become unbalanced.
An unbalanced tree has a relatively large depth for its size. The diagrams below shows the graph of the KdTree for a polygon containing 10,552 vertices. The unbalanced nature of the tree is evident from the few subtrees extending much deeper than the rest of the tree. The tree depth is 282, but most of that occurs in a single subtree.
An unbalanced KD-tree (size = 10,552, depth = 282)
Left: full graph. Right: close-up of graph structure
For large inputs querying such a deep subtree causes the recursive tree traversal to exceed the available call stack. Switching to an iterative implementation eliminates the stack overflow error, since the memory-based stack can grow to whatever size is needed.
However, stack overflow is not the only problem! Unbalanced KD-trees also exhibit poor query performance. This is because deep tree traversals make the search runtime more linear than logarithmic. Increasing the size of the query stack increases robustness, but does not improve the performance problem. The best solution is to build the KD-tree in a balanced way in the first place.
One way to do this is to to break up monotonic runs by randomizing the order of vertex insertion. This is actually implemented in the OverlayNG SnapRoundingNoder, as discussed in the post Randomization to the Rescue. However, using this approach in the SnappingNoder would effectively query the tree twice for each vertex. And in fact it's not necessary to randomize all the inserted points. A better-balanced tree can be produced by "seeding" it with a small set of well-chosen points.
But how can the points be chosen? They could be selected randomly, but true randomness can be quite "clumpy". Also, a non-deterministic algorithm is undesirable for repeatability. Both of these issues can be avoided by using a low-discrepancy quasi-random sequence. This concept has a fascinating theory, with many approaches available. The simplest to implement is an additive recurrence sequence with a constant α (the notation {...} means "fractional value of"):
R(α) : tn = { t0 + n * α }, n = 1,2,3,...
If α is irrational the sequence never repeats (within the limits of finite-precision floating point). The excellent article The Unreasonable Effectiveness of Quasirandom Sequences suggests that the lowest possible discrepancy occurs when α is the reciprocal of the Golden Ratio φ.
1 / φ = 0.6180339887498949...
Implementing this is straightforward, and it is remarkably effective. For the example above, the resulting KD-tree graph is significantly better-balanced, with lower depth (85 versus 282):
Adding KD-tree seeding to the SnappingNoder gives a huge performance boost to the QGIS test case. For the 526,466 input vertices the tree depth is reduced from 17,776 to 183 and the snapping time is improved by about 100 times, from ~80 s to ~800 ms! (In this case the overall time for the union operation is reduced by only about 50%, since the remainder of the overlay process is time-consuming due to the large number of holes in the output.)
In GEOS the performance difference is smaller, but still quite significant. The snapping time improves 14x, from 12.2 s to 0.82 s, with the union operation time decreasing by 40%, from 31 to 18 secs.
This improvement will be merged into JTS and GEOS soon.
In a previous post I described how the JTS Topology Suite operation IsSimpleOp has been completely rewritten to reduce code dependencies, improve performance, and provide a simpler, more understandable implementation. The post points out that the IsValidOp implementation would benefit from the same treatment. This work has now been carried out, with similar benefits achieved.
The original IsValidOp code used the GeometryGraph framework to represent the topology of geometries. This code reuse reduced development effort, but it carried some serious drawbacks:
The first check determines if the rings in the geometry have any invalid intersections. This includes the cases of a ring crossing itself or another ring, a ring intersecting itself or another ring in a line segment (a collinear intersection), and a ring self-touching (which is invalid in the OGC polygon model). This check can use the robust and performant machinery in the JTS noding package. If an invalid intersection is found, the algorithm can return that result immediately.
This check is an essential precursor to all the remaining checks. If a geometry has no invalid intersections this confirms that its rings do not self-cross or partially overlap. This means that the results of orientation and point-in-polygon computations are reliable. It also indicates that the rings are properly nested (although it remains to be determined if the nesting of shells and holes is valid).
Intersection search is the most computationally expensive phase of the algorithm. Because of this, information about the locations where rings touch is saved for use in the final phase of checking connected interiors.
The final check is that the interior of each polygon is connected. In a polygon with no self-touching rings, there are two ways that interior disconnection can happen:
As expected, the new codebase provides better performance, due to simpler code and the ability to short-circuit when an invalidity is found. Here's some performance comparisons for various datasets:
| Data | Size | New (ms) | Old (ms) | Improvement |
| world | 244 Polygons, 366K vertices | 340 | 1250 | 3.7 x |
| invalid-polys | 640 Polygons, 455K vertices | 126 | 334 | 2.6 x |
| valid-polys | 640 Polygons, 455K vertices | 244 | 487 | 2 x |
| australia | 1 MultiPolygon, 1,222K vertices | 1321 | 69026 | 52 x |
The improved IsValidOp will appear in JTS version 1.19. It has already been ported to GEOS, providing similar performance improvements. And since GEOS backs the PostGIS ST_IsValid function, that will become faster as well.
There is one remaining use of GeometryGraph in JTS for a non-constructive operation: the RelateOp class, used to compute all topological spatial predicates. Converting it should provide the same benefits of simpler code and improved performance. It should also improve robustness, and allow more lenient handling of invalid inputs. Watch this space!
Hard to believe that the JTS Topology Suite is almost 20 years old. That's 140 in dog years! Despite what they say about old dogs, one of the benefits of longevity is that you have the opportunity to learn a trick or two along the way. One of the key lessons learned after the initial release of JTS is that intersection (node) detection is a fundamental part of many spatial algorithms, and critical in terms of performance. This resulted in the development of the noding package to provide an API supporting many different kinds of intersection detection and insertion.
Prior to this, intersection detection was performed as part of the GeometryGraph framework, which combined it with topology graph formation and analysis. At the time this seemed like an elegant way to maximize code reuse across many JTS operations, including overlay, buffering, spatial predicates and validation. But as often the case, there are significant costs to such general-purpose code:
Another operation that relies on GeometryGraph is the IsSimpleOp class, which implements the OGC Simple Features isSimple predicate. The algorithm for isSimple essentially involves determining if the geometry linework contains a self-intersection. GeometryGraph is unnecessarily complex for this particular task, since there is no need to compute the entire topology graph in order to find a single self-intersection. Reworking the code to use the the MCIndexNoder class in the noding API produces a much simpler and more performant implementation. I also took the opportunity to move the code to the operation.valid package, since the operations of isSimple and isValid are somewhat complementary.
Now, isSimple is probably the least-used OGC operation. Its only real use is to test for self-intersections in lines or collections of lines, and that is not a critical issue for many workflows. However, there is one situation where it is quite useful: testing that linear network datasets are "vector-clean" - i.e. contain LineStrings which touch only at their endpoints.
A linear network containing non-simple intersections (isSimple == false)
To demonstrate the performance improvement, I'll use a dataset for Rivers of the US maintained by the US National Weather Service. It supplies two datasets: a full dataset of all rivers, and a subset of major rivers only. You might expect a hydrographic network to be "vector-clean", but in fact both of these datasets contain numerous instances of self-intersections and coincident linework.
Here's the results of running the isSimple predicate on the datasets. On the larger dataset the new implementation provides a 20x performance boost!
| Dataset | New time | Old time |
| Subset (909,865 pts) | 0.25 s | 1 s |
| Full (5,212,102 pts) | 2 s | 30 s |
The new codebase made it easy to add a functionality enhancement that computes the locations of all places where lines self-intersect. This can be used to for visual confirmation that the operation is working as expected, and to indicate places where data quality needs to be improved. Here's the non-simple intersection points found in the river network subset:
Closeups of some non-simple intersection locations:
IsSimpleOp is the easiest algorithm to convert over from using GeometryGraph. As such it serves as a good proof-of-viability, and establishes useful code patterns for further conversions.
Next up is to give IsValidOp the same treatment. This should provide similar benefits of simplicity and performance. And as always, porting the improved code to GEOS.
TLDR: JTS can now fix invalid geometry!
The JTS Topology Suite implements the Geometry model defined in the OGC Simple Features specification. An important part of the specification is the definition of what constitutes valid geometry. These are defined by rules about the structural and geometric characteristics of geometry objects. Some validity rules apply to all geometry; e.g. vertices must be defined by coordinates with finite numeric values (so that NaN and Inf ordinates are not valid). In addition, each geometric subtype (Point, LineString, LinearRing, Polygon, and Multi-geometry collections) has its own specific rules for validity.
The rules for Polygons and MultiPolygons are by far the most restrictive. They include the following constraints:
These rules guarantee that:
Given the highly-constrained definition of polygonal validity, it is not uncommon that real-world datasets contain polygons which do not satisfy all the rules, and hence are invalid. This occurs for various reasons:
Because of this, JTS does not enforce validity on geometry creation, apart from a few simple structural constraints (such as rings having identical first and last points). This allows invalid geometry to be represented as JTS geometry objects, and processed using JTS code. Some kinds of spatial algorithms can execute correctly on invalid geometry (e.g. determining the convex hull). But most algorithms require valid input in order to ensure correct results (e.g. the spatial predicates) or to avoid throwing exceptions (e.g. overlay operations). So the main reason for representing invalid geometry is to allow validity to be tested, to take appropriate action on failure.
Often users would like "appropriate action" to be Just Make It Work. This requires converting invalid geometry to be valid. Many spatial systems provide a way to do this:
But this has a been a conspicuous gap in the JTS API. While it is possible to test for validity, there has never been a way to fix an invalid geometry. To be fair, JTS has always had an unofficial way to make polygonal geometry valid. This is the well-known trick of computing geometry.buffer(0), which creates a valid output which often is a good match to the input. This has worked as a stop-gap for years (in spite of an issue which caused some problems, now fixed - see the post Fixing Buffer for fixing Polygons). However, using buffer(0) on self-intersecting "figure-8" polygons produces a "lossy" result. Specifically, it retains only the largest lobes of the input linework. This is undesirable for some uses (although it is advantageous in other situations, such as trimming off small self-intersections after polygon simplification).
So, it's about time that JTS stepped up to provide a supported, guaranteed way of fixing invalid geometry. This should handle all geometry, although polygonal geometry repair is the most critical requirement.
The other day this badge popped up on my GitHub profile:
Sure enough, there is JTS in the list of dependencies provided by the Jet Propulsion Laboratory (JPL).
So take that, Elon - JTS got there first! With a lot of help from NASA, of course.
A recent blog post by Elephant Tamer gives a critical appraisal of the improvements to overlay processing shipped in PostGIS 3.1 with GEOS 3.9. The author is disappointed that PostGIS still reports errors when overlay is used on invalid geometry. However, this is based on a misunderstanding of the technology.
GEOS 3.9 includes OverlayNG, ported from the JTS Topology Suite). It brings a major advance in overlay robustness, along with other improvements (described here and here). Previously, robustness limitations in the overlay algorithm could sometimes cause errors even for inputs which were topologically valid. This was doubly problematic because there was no fully effective way to process the input geometries to avoid the errors. Now, OverlayNG solves this problem completely. Valid inputs will always produce a valid and essentially correct(*) output.
(*) "Essentially" correct, because in order to achieve full robustness a snapping heuristic may be applied to the input geometry. However, this is done with a very fine tolerance, so should not appreciably alter the output from the theoretically correct value.
But for invalid inputs, OverlayNG will still report errors. The reason is that there is a wide variety of gruesome ways in which geometry can be invalid. Automated handling of invalidity would involve expensive extra processing, and also require making assumptions about what area a corrupt geometry is intended to represent. Rather than silently repairing invalid geometry and returning potentially incorrect results, the design decision is to report this situation as an error.
In fact, OverlayNG is able to handle "mildly" invalid polygons, as described in this post. This covers situations which are technically invalid according to the OGC SFS specification, but which still have well-defined topology. This includes self-touching rings (sometimes called "inverted polygons" or "exverted holes"), and zero-width gores and spikes.
Taking a detailed look at the data used in the blog post, we can see these improvements at work. The dataset is the ecology polygons obtained from the GDOS WFS server. This contains 7662 geometries, of which 10 are invalid. Using the old overlay algorithm, 9 of these invalid polygons cause TopologyException errors. Using OverlayNG, only 4 of them cause errors.
The polygons that can now be processed successfully are typical "OGC-invalid" situations, which do not materially affect the polygonal topology. These include self-touching rings with pinch points:
and zero-width gores:
Of the cases that still cause errors, two are classic small bow-tie errors:
And two are wildly invalid self-crossing rings:
The last two are good examples of geometry which is so invalid that it is impossible to unambiguously decide what area is represented (although ST_MakeValid will happily grind them into something that is technically valid).
Ultimately it is the user's responsibility to ensure that geometries to be processed by overlay (and many other PostGIS functions) have valid topology (as reported by ST_IsValid). Ideally this is done by correcting the data at source. But it can also be done a posteriori in the database itself, by either the ST_MakeValid function, or the well-known buffer(0) trick. (Which to use is a topic for another blog post...)
One improvement that could be made is to check for input validity when OverlayNG throws an error. Then PostGIS can report definitively that an overlay error is caused by invalid input. If there is an overlay error that is not caused by invalidity, the PostGIS team wants to hear about it!
And perhaps there is a case to be made for repairing invalid geometry automatically, even if the repair is suspect. Possibly this could be invoked via a flag parameter on the overlay functions. More research is required - feedback is welcome!
In a previous post I announced the new geosop command-line interface (CLI) for the GEOS geometry API. This post provides examples of using geosop for various tasks. (Refer to the README for more information about the various options used.)
The GEOS geometry API is used by many, many projects to do their heavy geometric lifting. But GEOS has always had a bit of a PR problem. Most of those projects provide a more accessible interface to perform GEOS operations. Some offer a high-level language like Python, R, or SQL (and these typically come with a REPL to make things even easier). Or there are GUIs like QGIS, or a command-line interface (CLI) like GDAL/OGR.
But you can't do much with GEOS on its own. It is a C/C++ library, and to use it you need to break out the compiler and start cutting code. It's essentially "headless". Even for GEOS developers, writing an entire C program just to try out a geometry operation on a dataset is painful, to say the least.
There is the GEOS XMLTester utility, of course. It processes carefully structured XML files, but that is hardly convenient. (And in case this brings to mind a snide comment like "2001 called and wants its file format back", XML actually works very well in JTS and GEOS as a portable and readable format for geometry tests. But I digress.)
JTS (on which GEOS is based) has the TestBuilder GUI, which works well for testing out and visualizing the results of JTS operations. JTS also has a CLI called JtsOp. Writing a GUI for GEOS would be a tall order. But a command-line interface (CLI) is much simpler to code, and has significant utility. In fact there is an interesting project called geos-cli that provides a simple CLI for GEOS. But it's ideal to have the CLI code as part of the GEOS project, since it ensures being up-to-date with the library code, and makes it easy to add operations to test new functionality.
This need has led to the development of geosop. It is a CLI for GEOS which performs a range of useful tasks:
geosop -a world.wkt -f wkt interiorPoint
geosop -a world.wkt --time buffer 1
geosop -a europe.wkb --collect -f wkb unaryUnion
The README gives many more examples of how to use the various command-line options. In a subsequent post I'll give some demonstrations of using geosop for various tasks including GEOS testing, performance tuning, and geoprocessing.
