RelateNG Performance

A previous post introduced a new algorithm in the JTS Topology Suite called RelateNG.  It computes topological relationships between geometries using the Dimensionally-Extended 9 Intersection Model  (DE-9IM) model.  

This algorithm is fundamental to a large proportion of spatial queries executed in numerous geospatial environments.  It would not be surprising to learn that the Relate algorithm is executed billions of times per day across the world's data centres.  Given this, performance is a critical metric.  

The original JTS RelateOp algorithm was hamstrung by its design, which required computing a full topology graph for each relationship evaluation.  The subsequent development of PreparedGeometry provided a significant performance boost for common spatial predicates such as intersects and covers.  A useful side-effect is that it is less susceptible to geometric robustness problems.  However, it has some drawbacks:

• it only implements a small set of predicates
• the design does not easily extend to more complex predicates
• it does not support GeometryCollection inputs
• it is a completely separate codebase to the general-purpose RelateOp, which increases maintenance effort and increases the risk of bugs and behaviour discrepancies

RelateNG solves all of these problems.  In addition, it provides much better performance for non-prepared (stateless) calls.  This is important for systems which do not execute spatial queries in a stateful way (due to architectural limitations, or simply lack of attention to performance), and thus cannot use prepared mode.

Performance Tests

Here are some examples of performance metrics.  

Note that the data used in the tests is sometimes subsetted to avoid skewing the results with many false results due to non-interaction between the input geometries.  This reflects the typical deployment of these predicates in systems which perform a primary extent filter (i.e. via a spatial index) before executing the full predicate.

Some of datasets are in geodetic coordinates, which strictly speaking JTS does not handle. But this does not matter for the purposes of performance testing.

Point / Polygon

This test queries a MultiPolygon of the outline of Europe (including larger islands) against a synthetic dataset of random points.  The test evaluates the intersects predicate, which is really the only one of interest in the Point / Polygon case.  

The test data metrics are:

• Query MultiPolygon: 461 polygons, 22,675 vertices
• Target Points: 10,000 points

Timings

Op	Intersects	Intersects Prep
Relate	61.7 s	21 ms
RelateNG	0.1 s	19 ms

The result clearly shows the enormous improvement in the stateless case, since RelateNG does not build topology on the input polygon.  The results are very similar in the prepared case.  This is as expected, since both codebases run a simple Point-in-Polygon test using a cached spatial index.

Line / Line

This test uses a dataset of major rivers of the continental United States.  A single river mainstem is queried against a subset of river branches using the intersects and touches relationships, in stateless and prepared modes.

The test data metrics are:

• Query Line: 6,975 vertices
• Target Lines: 407 lines with 47,328 vertices

 

Timings (per 100 executions of the operation)

Op	Intersects	Intersects Prep	Touches	Touches Prep
Relate	38.2 s	133 ms	36 s	N/A
RelateNG	1.18 s	142 ms	2.05 s	2.03 s

This shows the much better performance of RelateNG.  RelateNG can evaluate touches in prepared mode, but the performance is similar to stateless mode, since currently indexes are not cached for Line/Line cases.  This could be improved in a future release.

Polygon / Polygon

This test uses data from two polygonal datasets:

• The Britsh Columbia Bedrock Geology polygonal coverage
• The GADM Level 2 boundaries for Canada

The tests query an administrative unit polygon against a subset of Bedrock polygons which intersect its extent, using intersects and covers in stateless and prepared modes. 

The test data metrics are:

• GADM unit: 4,017 vertices
• Bedrock polygons: 4,318 polygon with 337,650 vertices

Timings (per 100 executions of the operation)

Op	Intersects	Intersects Prep	Covers	Covers Prep
Relate	61.7 s	534 ms	54.9 s	842 ms
RelateNG	5.8 s	595 ms	6.4 s	943 ms

Polygon / Polygon - Custom Relate Pattern

The tests shows the ability of RelateNG to efficiently evaluate arbitrary Relate Intersection Matrix patterns.  In this case, the pattern used is F***0****, which corresponds to a relationship which could be called "point-touches": two geometries have boundaries which intersect in only one (or more) points (dimension 0), but whose interiors do not intersect.

This test uses data from The GADM Level 1 boundaries for Canada.  This has the advantage that Canada contains a rare example of 4 boundaries meeting at a single point (the jurisdictions of Saskatchewan, Manitoba, Northwest Territories, and Nunavut).

The test data metrics are:

• GADM Canada Level 1: 13 polygons, 4,005,926 vertices

Timings

Op	F***0****	F***0**** Prep
Relate	504 s	N/A
RelateNG	9.8 s	6.6 s

As expected, for this test RelateNG is far more performant than Relate.  This is due to its ability to analyze Intersection Matrix patterns and perform only topology tests that are necessary, as well as not building a full topological structure of the inputs.  The test shows the effect of caching spatial indexes in prepared mode, although stateless mode is very efficient as well.

Analysis of Results

Clearly RelateNG is far more performant than Relate when running in non-prepared mode.  The PreparedGeometry implementation is slightly faster (which confirms the efficiency of its original design), but not by very much.  The difference may be a consequence of the general-purpose and thus more complex code and data structures in RelateNG. Closing this gap may be an area for future research.

One thing is certain: the RelateNG algorithm design provides much more scope for adding case-specific optimizations.  If you have a significant use case which could be improved by further targeted optimization improvements, let me know in the comments!

It will be interesting to re-evaluate these tests once RelateNG is ported to GEOS.  Sometimes (but not always) a C++ implementation can be significantly faster than Java, due to more opportunities for code and compiler optimization.

JTS Topological Relationships - the Next Generation

The most fundamental and widely-used operations in the JTS Topology Suite are the ones that evaluate topological relationships between geometries.  JTS implements the Dimensionally-Extended 9 Intersection Model (DE-9IM), as defined in the OGC Simple Features specification, in the RelateOp API.  

DE-9IM matrix for overlapping polygons

The RelateOp algorithm was the very first one implemented during the initial JTS development, over 20 years ago.  At that time it was an appealing idea to implement a general-purpose topology framework (the GeometryGraph package), and use it to support topological predicates, overlay, and buffering.  However, some disadvantages of this approach have become evident over time:

• the need to create a topological graph structure limits the ability to improve performance. This has led to the implementation of PreparedGeometry - but that adds further complexity to the codebase, and supports only a limited set of predicates.
• a large number of code dependencies make it hard to fix problems and improve semantics
• constructing a full topology graph increases exposure to geometric robustness errors
• GeometryCollections are not supported (initially because the OGC did not define the semantics for this, and now because adding this capability is difficult)

There's an extensive list of issues relating (ha!) to RelateOp here.

The importance of this functionality is especially significant since the same algorithm is implemented in GEOS.  That codebase is used to evaluate spatial queries in popular spatial APIs such as Shapely and R-sf, and numerous systems such as PostGIS, DuckDB, SpatialLite, QGIS, and GDAL (to name just a few).  It would not be surprising to learn that the RelateOp algorithm is executed billions of times per day across the world's CPUs.

During the subsequent years of working on JTS I realized that there was a better way to evaluate topological relationships.  It would required a ground-up rewrite, but would avoid the shortcomings of RelateOp and provide better performance and a more tractable codebase.  Thanks to my employer Crunchy Data I have finally been able to make this idea a reality.  Soon JTS will provide a new algorithm for topological relationships called RelateNG.

Key Features of RelateNG

The RelateNG algorithm incorporates a broad spectrum of improvements over RelateOp in the areas of functionality, robustness, and performance.  It provides the following features:

• Efficient short-circuited evaluation of topological predicates (including matching custom DE-9IM matrix patterns)
• Optimized repeated evaluation of predicates against a single geometry via cached spatial indexes (AKA "prepared mode")
• Robust computation (only point-local geometry topology is computed, so invalid topology does not cause failures)
• GeometryCollection inputs containing mixed types and overlapping polygons are supported, using union semantics.
• Zero-length LineStrings are treated as being topologically identical to Points.
• Support for BoundaryNodeRules.

Using the RelateNG API

The main entry point is the RelateNG class.  It supports evaluating topological relationships in three different ways:

• Evaluating a standard OGC named boolean binary predicate, specified via a TopologyPredicate instance.  Standard predicates are obtained from the RelatePredicate factory functions intersects, contains, overlaps, etc.
• Testing an arbitrary DE-9IM relationship by matching an intersection matrix pattern (e.g. "T**FF*FF*", which is the pattern for a relation called Contains Properly).
• Computing the full value of a DE-9IM IntersectionMatrix.

RelateNG operations can be executed in two modes: stateless and prepared.  Stateless mode is provided by the relate static functions, which accept two input geometries and evaluate a predicate.  For prepared mode, a single geometry is provided to the prepare functions to create a RelateNG instance.  Then the evaluate instance methods can be used to evaluate spatial predicates on a sequence of geometries. The instance creates spatial indexes (lazily) to make topology evaluation much more efficient.  Note that even stateless mode can be significantly faster than the current implementation, due to short-circuiting and other internal heuristics.

Here is an example of matching an intersection matrix pattern, in stateless mode:

boolean isMatched = RelateNG.relate(geomA, geomB, "T**FF*FF*");

Here is an example of setting up a geometry in prepared mode, and evaluating a named predicate on it:

RelateNG rng = RelateNG.prepare(geomA);
for (Geometry geomB : geomSet) {
    boolean predValue = rng.evaluate(geomB, RelatePredicate.intersects());
}

Rolling It Out

It's exciting to launch a major improvement on such a core piece of spatial functionality.  The Crunchy spatial team will get busy on porting this algorithm to GEOS.  From there it should get extensive usage in downstream projects.  We're looking forward to hearing feedback from our own PostGIS clients as well as other users.  We're always happy to be able to reduce query times and equally importantly, carbon footprints.  

In further blog posts I'll describe the RelateNG algorithm design and provide some examples of performance metrics.

Future Ideas

The RelateNG implementation provides an excellent foundation to build out some interesting extensions to the fundamental DE-9IM concept.

Extended Patterns

The current DE-9IM pattern language is quite limited.  In fact, it's not even powerful enough to express the standard named predicates.  It could be improved by adding features like:

• disjunctive combinations of patterns.  For example, touches is defined by  "FT******* | F**T***** | F***T****"
• dimension guards to specify which dimensions a pattern applies to.  For example, overlaps is defined by "[0,2] T*T***T** | [1] 1*T***T**"
• while we're at it, might as well support dotted notation and spaces for readability; e.g. "FT*.***.***"

Approximate Spatial Relationships

A requirement that comes up occasionally is to compute approximate spatial relationships between imprecise data using a distance tolerance.  This is useful in datasets where features don't match exactly due to data inaccuracy. Because RelateNG creates topology only in the neighbourhood of specified points, it should be possible to specify the neighbourhood size using a distance tolerance.  Essentially, vertices and line intersections will "snap together" to determine the effective topology.  Currently within-distance queries are often used to compute "approximate intersects".  Adding a distance tolerance capability to RelateNG will allow other kinds of approximate relationships to be evaluated as well.

Further Performance Optimizations?

A challenge with implementing algorithms over a wide variety of spatial types and use cases is how to provide general-purpose code which matches (or exceeds) the efficiency of more targeted implementations.  RelateNG analyzes the input geometries and the predicate under evaluation to tune strategies to reduce the amount of work needed to evaluate the DE-9IM.  It may be that profiling specific use cases reveals further hotspots in the code which can be improved by additional optimizations.  

Curve Support

GEOS has recently added support for representing geometries with curves.  The RelateNG design is modular enough that it should be possible to extend it to allow evaluating relationships for geometries with curves.  

Relational Properties of DE-9IM spatial predicates

There is an elegant mathematical theory of binary relations.  Homogeneous relations are an important subclass of binary relations in which both domains are the same.  A homogeneous relation R is a subset of all ordered pairs (x,y) with x and y elements of the domain.  This can be thought of as a boolean-valued function R(x,y), which is true if the pair has the relationship and false if not.

The restricted structure of homogeneous relations allows describing them by various properties, including:

Reflexive - R(x, x) is always true

Irreflexive - R(x, x) is never true

Symmetric - if R(x, y), then R(y, x)

Antisymmetric - if R(x, y) and R(y, x), then x = y

Transitive - if R(x, y) and R(y, z), then R(x, z)

The Dimensionally Extended 9-Intersection Model (DE-9IM) represents the topological relationship between two geometries.  

Various useful subsets of spatial relationships are specified by named spatial predicates.  These are the basis for spatial querying in many spatial systems including the JTS Topology Suite, GEOS and PostGIS.

The spatial predicates are homogeneous binary relations over the Geometry domain  They can thus be categorized in terms of the relational properties above.  The following table shows the properties of the standard predicates:

Predicate	Reflexive / Irreflexive	Symmetric / Antisymmetric	Transitive
Equals	R	S	T
Intersects	R	S	-
Disjoint	R	S	-
Contains	R	A	-
Within	R	A	-
Covers	R	A	T
CoveredBy	R	A	T
Crosses	I	S	-
Overlaps	I	S	-
Touches	I	S	-

Notes

• Contains and Within are not Transitive because of the quirk that "Polygons do not contain their Boundary" (explained in this post).  A counterexample is a Polygon that contains a LineString lying in its boundary and interior, with the LineString containing a Point that lies in the Polygon boundary.  The Polygon does not contain the Point.  So Contains(poly, line) = true and Contains(line, pt) = true, but Contains(poly, pt) = false.   The predicates Covers and CoveredBy do not have this idiosyncrasy, and thus are transitive.

JTS IsValidOp Built Back Better

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 GeometryGraph structure was used for many JTS operations, including overlay, buffer, and spatial predicate evaluation.  This made the code complex, hard to understand, and difficult to maintain and enhance
• The GeometryGraph structure computes the full topology graph of the input geometry, in order to support constructive operations such as overlay and buffer.  This makes it slower and more subject to robustness problems.  Non-constructive operations such as IsValidOp and IsSimpleOp can compute topology "on-the-fly", which allows short-circuiting processing as soon as an invalidity is found.  

The new IsValidOp implementation is more self-contained, relying only on the noding framework and some basic spatial predicates.  Using the MCIndexNoder allows short-circuited detection of some kinds of invalidity, which improves performance.  And it will be easy to switch to a different noding strategy should a superior one arise.

However, dropping GeometryGraph means that the topology information required to confirm validity needs to be computed some other way.  The goal is to compute just enough topological structure to evaluate validity, and do that efficiently on an as-needed basis.  This required some deep thinking about validity to pare the topology information required down to a minimum.  This was assisted by the extensive set of unit tests for IsValidOp, which ensured that all possible situations were handled by the new logic (although there was an untested situation which uncovered a bug late in development.) The resulting algorithm uses some elegant logic and data structures, which are explained in detail below. 

OGC Validity

JTS implements the OGC Simple Features geometry model.  To understand the validity algorithm, it helps to review the rules of geometry validity in the OGC specification:

For non-polygonal geometry types, the validity rules are simple;

1. Coordinates must contain valid numbers
2. LineStrings must not be zero-length
3. LinearRings must be simple (have no self-intersections)

For polygonal geometry, on the other hand, the validity rules are quite stringent (sometimes considered too much so!).  The rules are:

1. Rings must be simple (have no self-intersections)
2. Rings must not cross, and can touch only at discrete points
3. Holes must lie inside their shell
4. Holes must not lie inside another hole of their polygon
5. The interiors of polygons must be connected
6. In a MultiPolygon, the element interiors must not intersect (so polygons cannot overlap)
7. In a MultiPolygon, the element boundaries may intersect only at a finite number of points

Validity Algorithm for Polygonal Geometry

The new IsValidOp uses an entirely new algorithm to validate polygonal geometry.  It has the following improvements:

• It supports short-circuiting when an invalid situation is found, to improve performance
• Spatial indexing is used in more places
• The code is simpler and more modular.  This makes it easier to understand, and easier to adapt to new requirements (e.g. such as the ability to validate self-touching rings described below)

The algorithm consists of a sequence of checks for different kinds of invalid situations.  The order of the checks is important, since the logic for some checks depends on previous checks passing. 

1. Check Ring Intersections

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.

Kinds of invalid intersections: 

(i) a self-touch; (ii) collinear; (iii) segment crossing; (iv) vertex crossing

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.  

2. Check Shell and Hole nesting

If no invalid intersections are found, then the rings are properly nested.  Now it is necessary to verify that shells and holes are correctly positioned relative to each other:

• A hole must lie inside its shell
• A hole may not lie inside another hole. 
• In MultiPolygons, a shell may not lie inside another shell, unless the inner shell lies in a hole of the outer one.  

Ring position can often be tested using a simple Point-In-Polygon test.  But it can happen that the ring point tested lies on the boundary of the other ring.  In fact, it is possible that all points of a ring lie on another ring.  In this case the Point-In-Polygon test is not sufficient to provide information about the relative position of the rings.  Instead, the topology of the incident segments at the intersection is used to determine the relative position.  Since the rings are known to not cross, it is sufficient to test whether a single segment of one ring lies in the interior or exterior of the other ring.

MultiPolygon with an element where all vertices lie on the boundary

Checking nesting of holes and shells requires comparing all interacting pairs of rings.  A simple looping algorithm has quadratic complexity, which is too slow for large geometries.  Instead, spatial indexes (using the STRtree) are used to make this process performant. 

3. Check Connected Interior 

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:

• a chain of one or more touching holes touches the shell at both ends, thus splitting the polygon in two
• a chain of of two or more touching holes forms a loop, thus enclosing a portion of the polygon interior  

To check that these conditions do not occur, the only information needed is the set of locations where two rings touch.  These induce the structure of an undirected graph, with the rings being the graph vertices, and the touches forming the edges.  The situations above correspond to the touch graph containing a cycle (the splitting situation is a cycle because the shell ring is also a graph vertex).  (Equivalently, the touch graph of a valid polygon is a tree, or more accurately a forest, since there may be sets of holes forming disconnected subgraphs.) So a polygon can be verified to have a connected interior by checking that its touch graph has no cycles.  This is done using a simple graph traversal, detecting vertices (rings) which are visited twice.

Situations causing a disconnected interior: 

(i) a chain of holes; (ii) a cycle of holes

If no invalid situations are discovered during the above checks, the polygonal geometry is valid according to the OGC rules.

Validity with Self-touching Rings

Some spatial systems allow a polygon model in which "inverted shells" are valid.  These are polygon shells which contain a self-touch in a way that encloses some exterior area.  They also allow "exverted holes", which contain self-touches that disconnect the hole into two or more lobes.  A key point is that they do not allow "exverted shells" or "inverted holes"; self-touches may disconnect the polygon exterior, but not the interior.  Visually this looks like:

Left: Inverted shell and exverted hole (valid)

Right: Exverted shell and inverted hole (invalid)

This kind of topology is invalid under the OGC polygon model, which prohibits all ring self-intersections.  However, the more lenient model still provides unambiguous topology, and most spatial algorithms can handle geometry of this form perfectly well.  So there is a school of thought that considers the OGC model to be overly restrictive.

To support these systems JTS provides the setSelfTouchingRingFormingHoleValid option for IsValidOp to allow this topology model to be validated.  Re-implementing this functionality in the new codebase initially presented a dilemma. It appeared that building the touch graph to check connectivity required computing the geometry of the additional rings formed by the inversions and exversions.  This needed a significant amount of additional code and data structures, eliminating the simplicity of the graph-based approach.  

However, it turns out that the solution is much simpler.  It relies on the fact that valid self-touches can only disconnect the exterior, not the interior.  This means that self-touches of inverted shells and exverted holes can be validated by the local topology at the self-intersection node.  This condition can be tested by the same logic already used to test for nested touching shells.  Even better, if the self-touches are valid, then the touch-graph algorithm to check connectivity still works.   (This is the same phenomenon that allows most spatial algorithms to work correctly on the inverted/exverted ring model.)  

Performance

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

• world is the standard dataset of world country outlines
• invalid-polys is a dataset of all invalid polygons.  valid-polys is the same dataset processed by GeometryFixer.  The timings show the effect of short-circuiting invalidity to improve performance. 
• australia is a geometry of the Australian coastline and near-shore islands, with 6,697 polygon elements.  The dramatic difference in performance reflects the addition of a spatial index for checking nested shells.

Australia and 6,696 islands

Next Steps

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!

JTS IsSimple gets simpler (and faster)

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:

• The overall codebase is substantially more complex
• A performance penalty is imposed on algorithms which don't require topology construction
• Algorithms are harder to read and understand.  
• The code is brittle, and so hard to modify
• Porting the code is more difficult  

Because of this, a focus of JTS development is to free operations from their dependency on GeometryGraph - with the ultimate goal of expunging it from the JTS codebase.  A major step along this road was the rewrite of the overlay operations.

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

Finding Non-Simple Locations

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.

The complexity of simple

Recently Sean Gillies and others have been looking at the concept of geometric simplicity as implemented in JTS/GEOS, specifically in the case of polygons.  The TL;DR is that JTS has a bug because it reports the following polygon as simple:

POLYGON ((39.14513303 23.75100816, 97.05372209 40.54550088, 105.26527014 48.13302213, 100.91752685 58.43336815, 71.56081448 83.55114801, 60.71189168 86.25316099, 62.00469808 75.1478176, 83.16310007 42.82071673, 92.82305862 37.19175582, 95.99401129 26.47051246, 106.22054482 15.51975192, 39.14513303 23.75100816))

This polygon is invalid, because the shell ring self-intersects.  This means the ring is non-simple - but it's not so clear-cut for the polygon itself.

In fact, reporting the polygon as simple isn't a bug - it was a deliberate design decision.  The original rationale was:

• the OGC Simple Features Specification in section 2.1.10 defining polygon geometries states that "Polygons are simple geometries".  (This is consistent with the requirement that valid polygons contain no self-intersections).  
• JTS usually follows the principle that computed results are only well-defined when the input geometries are valid.  This is because it can be expensive to check for validity in order to handle invalid inputs gracefully.  Validity testing is factored out into a separate method (isValid) in order to allow the client to decide when to incur the code of executing it.

For these reasons the initial JTS implementation of isSimple returned true for all polygonal inputs, thus avoiding the cost of checking for self-intersections.

However, perhaps it would be more useful to actually carry out some testing of simplicity.  This raises the question of what exactly are the semantics of isSimple for polygons?  Looking to the SFS for guidance, it has the following somewhat imprecise general definition of simplicity:

the Geometry has no anomalous geometric points, such as self intersection or self tangency. The description of each instantiable geometric class will include the specific conditions that cause an instance of that class to be classified as not simple.

For polygons the spec simply states:

Polygons are simple geometries

(In contrast, for linear geometry the specification has a rigorous mathematical definition of simplicity.)

Apart from the "always simple" semantics, there are at least two other options:

1. Polygons are simple if their component rings are simple (i.e. have no self-intersections)
2. Polygons are simple if they have no topological anomalies.  This is equivalent to checking whether the polygon is valid.

Because option #2 simply replicates the semantics of isValid, it seems more useful to adopt option #1 to provide more possibilities for analyzing polygon topology (but I'm definitely open to suggestions about why #2 would be better!)

This change to the semantics of isSimple has been implemented in JTS and will appear in the 1.13 release.  The code will also be extended to handle GeometryCollections, with the semantics that they are simple if all their components are simple.

PreparedGeometry - efficient batch spatial operations for JTS

A common usage pattern for JTS spatial operations (in particular, predicates such as intersects and contains, and overlay operations such as intersection and union) involves the operation being invoked between a single target geometry and a batch of test geometries. For instance, this occurs in queries in a spatial database such as PostGIS, where each geometry in the driving table interacts with a set of nearby geometries via the required operation. By taking advantage of the fact that the target geometry does not change, it is possible to significantly reduce the overall execution cost of the operation. Information such as an index on the line segments in the target geometry can be computed once and cached for reuse with all of the test geometries.

In order for this to work, the calling code must be able to indicate that the operation is being invoked in a batch setting. The standard JTS Geometry API doesn't allow differentiating this case, so the new concept of PreparedGeometry has been developed. (This name intentionally echoes the concept of prepared statement in SQL API's, which implements the same strategy.) Creating a PreparedGeometry on a target Geometry enables precomputation of various data structures, which then allows subsequent operations to be evaluated with maximum efficiency.

In addition to caching, it is also possible to optimize the evaluation of spatial predicates. Optimization takes the form of short-circuiting predicate evaluation when certain situations are detected. While not all OGC spatial predicates are amenable to optimization, fortunately the most commonly usd predicates, intersects and contains, can be optimized significantly. Examples of situations which can be optimized are:

• For intersects, any intersection between line segments indicates a true result. If there is no intersection between any line segments of the operands, and the target geometry is polygonal, then one or more fast point-in-polygon tests can be performed to compute the result value.
• Since it provides more information than intersects, the contains predicate is more complex to optimize. In the case where the target geometry is polygonal, evaluation can be short-circuited when the test geometry has no line segment intersection. In this case, containment can be tested by an appropriate set of point-in-polygon tests. It is also possible to short-circuit in some specific cases where a proper intersection is detected between two line segments. Both of these situations can be tested efficiently.

The performance gains from using PreparedGeometry in real-world situations can be dramatic - up to 40 times improvement in speed. This functionality will appear in the next release of JTS, and hopefully will make its way into GEOS and PostGIS as well.

Quirks of the "Contains" Spatial Predicate

One of the most useful OGC standards is their specification of the Dimensionally-Extended 9 Intersection Model for spatial relationships. (To give credit where is it due, this was originally developed by Egenhofer, Clementini and others). In its most general form this model is quite complex, so to make it usable for mortal programmers, a set of commonly-useful "named predicates" have been specified. These include obvious relationships such as intersects, disjoint, touches, equals, and contains.

But are they really so obvious? In fact not - several of them have subtle aspects to their definition which are contrary to intuition.

In particular, "contains" (and its converse "within") has an aspect of its definition which may produce unexpected behaviour. This quirk can be expressed as "Polygons do not contain their boundary". More precisely, the definition of contains is:

Geometry A contains Geometry B iff no points of B lie in the exterior of A, and at least one point of the interior of B lies in the interior of A

That last clause causes the trap - because of it, a LineString which is completely contained in the boundary of a Polygon is not considered to be contained in that Polygon!

This behaviour could easily trip up someone who is simply trying to find all LineStrings which have no points outside a given Polygon. In fact, this is probably the most common usage of contains. For this reason it's useful to define another predicate called covers, which has the intuitively expected semantics:

Geometry A covers Geometry B iff no points of B lie in the exterior of A

Its converse coveredBy is also useful as well. It's a bit of a mystery why OGC did not define this predicate. In any case, Oracle Spatial does provide this predicate. I have added it to JTS as well.

There's a further bonus to defining the covers predicate. It is much easier to optimize the common use case:

covers(Rectangle, Geometry)

All that is necessary to determine this condition is to perform a simple bounding box comparison. This is not possible with contains, because even if the bounding box of Geometry is covered by the Rectangle, a further expensive operation is required to test if the Geometry lies wholly in the boundary of the Rectangle (in which case the predicate fails).

Covers "simplifies" the defintion of contains by making it more general (inclusive). There's another possibility as well, which is to simplify in the direction of being less inclusive. This would produce a predicate which might be called containsProperly, with the definition:

Geometry A containsProperly Geometry B iff all points of B lie in the interior of A

Interestingly, some recent work I'm doing indicates that this predicate might play a very useful role - but that's a story for another post.


References

DE-9IM: Clementini, Eliseo, Di Felice, and van Osstrom; "A Small Set of Formal Topological Relationships Suitable for End-User Interaction," D. Abel and B.C. Ooi (Ed.), Advances in Spatial Database—Third International Symposium. SSD '93. LNCS 692. Pp. 277-295.

9 Intersection model: M. J. Egenhofer and J. Herring; "Categorizing binary topological relationships between regions, lines, and points in geographic databases," Tech. Report, Department of Surveying Engineering, University of Maine, Orono, ME 1991.