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.  

Query KD-trees 100x faster with this one weird trick!

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(α) :  t_n = { t₀ + 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):

Better-balanced KD-tree (size = 10,552,  depth = 85)

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.

Introducing geosop - a CLI for GEOS

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:

• Run GEOS operations to confirm their semantics
• Test the behaviour of GEOS on specific geometric data
• Time the performance of operation execution
• Profile GEOS code to find hotspots
• Check memory usage characteristics of GEOS code
• Generate spatial data for use in visualization or testing
• Convert datasets between WKT and WKB

geosop has the following capabilities:

• Read WKT and WKB from files, standard input, or command-line literals
• Execute GEOS operations on the list(s) of input geometries.  Binary operations are executed on every pair of input geometries (i.e. the cross join aka Cartesian product)
• Output geometry results in WKT or WKB (or text, for non-geometric results)
• Display the execution time of data input and operations
• Display a full log of the command processing

Here's a look at how it works.  

geosop -h gives a list of the options and operations available:

geosop - GEOS v. 3.10.0dev

Executes GEOS geometry operations

Usage:

  geosop [OPTION...] opName opArg

  -a arg               source for A geometries (WKT, WKB, file, stdin,

                       stdin.wkb)

  -b arg               source for B geometries (WKT, WKB, file, stdin,

                       stdin.wkb)

      --alimit arg     Limit number of A geometries read

  -c, --collect        Collect input into single geometry

  -e, --explode        Explode result

  -f, --format arg     Output format

  -h, --help           Print help

  -p, --precision arg  Sets number of decimal places in WKT output

  -r, --repeat arg     Repeat operation N times

  -t, --time           Print execution time

  -v, --verbose        Verbose output

Operations:

  area A - computes area for geometry A

  boundary A - computes boundary for geometry A

  buffer A N - cmputes the buffer of geometry A

  centroid A - computes centroid for geometry A

  contains A B - tests if geometry A contains geometry B

  containsPrep A B - tests if geometry A contains geometry B, using PreparedGeometry

  containsProperlyPrep A B - tests if geometry A properly contains geometry B using PreparedGeometry

  convexHull A - computes convexHull for geometry A

  copy A - computes copy for geometry A

  covers A B - tests if geometry A covers geometry B

  coversPrep A B - tests if geometry A covers geometry B using PreparedGeometry

  difference A B - computes difference of geometry A from B

  differenceSR A B - computes difference of geometry A from B rounding to a precision scale factor

  distance A B - computes distance between geometry A and B

  distancePrep A B - computes distance between geometry A and B using PreparedGeometry

  envelope A - computes envelope for geometry A

  interiorPoint A - computes interiorPoint for geometry A

  intersection A B - computes intersection of geometry A and B

  intersectionSR A B - computes intersection of geometry A and B

  intersects A B - tests if geometry A and B intersect

  intersectsPrep A B - tests if geometry A intersects B using PreparedGeometry

  isValid A - tests if geometry A is valid

  length A - computes length for geometry A

  makeValid A - computes makeValid for geometry A

  nearestPoints A B - computes nearest points of geometry A and B

  nearestPointsPrep A B - computes nearest points of geometry A and B using PreparedGeometry

  polygonize A - computes polygonize for geometry A

  reducePrecision A N - reduces precision of geometry to a precision scale factor

  relate A B - computes DE-9IM matrix for geometry A and B

  symDifference A B - computes symmetric difference of geometry A and B

  symDifferenceSR A B - computes symmetric difference of geometry A and B

  unaryUnion A - computes unaryUnion for geometry A

  union A B - computes union of geometry A and B

  unionSR A B - computes union of geometry A and B

Most GEOS operations are provided, and the list will be completed soon.

Some examples of using geosop are below.

• Compute the interior point for each country in a world polygons dataset, and output them as WKT:

geosop -a world.wkt -f wkt interiorPoint

• Determine the time required to compute buffers of distance 1 for each country in the world:

geosop -a world.wkt --time buffer 1

• Compute the union of all countries in Europe:

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.

Future Work

There's potential to make geosop even more useful: 

• GeoJSON is a popular format for use in spatial toolchains.  Adding GeoJSON reading and writing would allow geosop to be more widely used for geo-processing.  
• Adding SVG output would provide a way to visualize the results of GEOS operations.   
• Improve support for performance testing by adding operations to generate various kinds of standard test datasets (such as point grids, polygon grids, and random point fields).  

And of course, work will be ongoing to keep geosop up-to-date as new operations and functionality are added to GEOS.

Variable-distance buffering in JTS

The operation of buffering a geometry is a core geospatial concept.  Standard buffers are computed using a fixed distance around the input geometry. JTS has provided a buffer implementation since its inception, and this is used in GEOS and all the other downstream projects as well. 

One way to generalize the buffer concept is to allow the buffer distance to vary along the geometry.  As often the case with geospatial concepts this construction has a few different names, including tapered buffer, cone buffer, varying buffer, and variable-distance buffer.  The classic use case for variable-distance buffers is generating polygons for the "cone of uncertainty" along predicted hurricane tracks:

Another use case is depicting rivers showing the width of reaches along the river course. 

I took a crack at prototyping "variable width buffers" a few years ago, in the JTS Lab.  Recently there were a couple of GIS StackExchange posts looking for this functionality in PostGIS and Shapely   They were good motivation to buff up the prototype and move it into JTS core.  But I was dismayed to realize that the output had some serious deficiencies, due to an overly-simplistic algorithm.  The idea in the Lab code was appealingly simple: compute buffer circles around each line vertex at the specified distance, and union (merge) them with trapezoids computed around the connecting line segment.  But this produced ugly discontinuities for large deltas in buffer distances. 

For such a seemingly simple concept there is surprisingly little prior art to be found.  Even the GIS That Shall Not Be Named does not seem to provide a native variable buffer function (although recently I found a couple of user contribs which do, to some extent).   There's an implementation in QGIS - but since it seems to be based on the original problematic JTS code that didn't really help.   So I had the fun of coding it up from scratch.  

The problem with the original code is that it should have used the outer tangent lines to the buffer circles at each vertex.  Wikipedia has a good discussion of this construction.  Even better, it provides an elegant mathematical algorithm, which worked perfectly when coded up.

The construction computes a single tangent segment.  The opposite segment is simply the reflection in the line between the circle centres.  The geometric math to handle this is now provided for reuse as  LineSegment.reflect().  Finally, it is notoriously tricky to produce buffer curves with high quality in the fine details.  In this case, generating the line segment buffer caps required care to avoid creating out-of-phase vertices, which would produce a "bumpy" result when merged.

This is now available in JTS as the VariableBuffer class.  In addition to the general-purpose API which allows specifying a distance at each vertex, it provides a couple of simpler functions which accept a start and end distance; and a start, middle and end distance.  The first is an easy way to produce the classic cone of uncertainty.  The latter is useful for buffering rings, or just creating interesting shapes.

The next step is to port this to GEOS. Then it can be exposed in PostGIS and all the other downstream projects like Shapely and R-SF, and maybe even QGIS.

Future Work

There's some enhancements that could be made:

• Some or all of the buffer end cap and join styles provided by the classic buffer operation could be supported, as well as the ability to set the quadrant segment count
• Variable buffering of polygons and multigeometries is straightforward to add.  
• It might be possible to allow different distances for each side of the input line (although it may be tricky to get the joins right)
• The approach of using union to merge the individual segment buffers works well. But it would improve performance to feed the entire generated variable buffer curve directly into the existing buffer generation code

More experimentally, variable buffering might provide a way to construct a decent approximation to a true geodetic buffer.  The vertex distances would be computed using true geodetic distance.  There might be some complications around long line segments, but perhaps these can be finessed (e.g. by densification). 

Fast Geometry Distance in JTS

The second-most important criteria for a spatial algorithm is that it be fast.  (The most important is that it's correct!)  Many spatial algorithms have a simple implementation available, but with performance of O(n²) (or worse).  This is unacceptably slow for production usage, since it results in long runtimes for data of any significant size. In JTS a lot of effort has gone into identifying O(n²) performance hotspots and engineering efficient replacements for them.

One long-standing hotspot is the algorithm for computing Euclidean distance between geometries.  The obvious distance algorithm is a brute-force O(MxN) comparison between the vertices and edges (facets) of the input geometries.  This is simple to implement, but very slow for large inputs.  Surprisingly, there seems to be little in the computational geometry literature about more efficient distance algorithms.  Perhaps because of this, many geometric libraries provide only the slow brute force algorithm - including JTS (until now). 

Happily, it turns out there is a faster approach to distance computation.  It uses data structures and algorithms which are already provided in JTS, so it's relatively easy to implement. The basic idea is to build a spatial index on each of the input geometries, and then use a Branch-and-Bound search algorithm to efficiently traverse the index trees to find for the minimum distance between geometry facets.  This is a generalization of the R-tree Nearest Neighbour algorithm described in the classic paper by Rousssopoulos et al.  [1]. 

JTS has the STRtree R-tree index implementation (a packed R-tree using the Sort-Tile-Recursive algorithm). This has recently been enhanced with several kinds of nearest-neighbour searches.  In particular, it now supports a method to find the nearest neighbours between two different trees.  The IndexedFacetDistance class uses this capability to implement fast distance searching on the facets of two geometries.

Another benefit of this approach is that it allows caching the index of one geometry.  This further increases performance in the common case of repeated distance calculations against a fixed geometry.

The performance improvement is impressive.  Here's the timings for computing the distance from Antarctica to other world countries:

Source Data size	Target Data size	Time Indexed	Time Brute-Force	Improvement
1 polygon (19,487 vertices)	244 polygons (366,951 vertices)	164 ms	136 s	x 830

Branch-and-bound search also speeds up isWithinDistance queries.  Here's a within-distance selection query between another antipodean continent and a large set of small rectangles:

Source Data size	Target Data size	Time	Time Brute-Force	Improvement
1 polygon (7,316 vertices)	100,000 polygons (500,000 vertices)	53 ms	10.03 s	x 19

A small fly in the algorithmic ointment is that Indexed Distance is not always better than the brute-force approach.  For small geometries (such as points or rectangles) a simple scan is actually faster, since it avoids the overhead of building indexes.  It may be possible to determine a tuning parameter that allows automatically choosing the fastest option.  Or the client can choose the faster approach, using knowledge of the use case.

Future Work

A few further ideas to build or investigate:

• Implement a caching FastDistanceOp using IndexedFacetDistance and indexed Point-In-Polygon.  This can be used to add a fast distance() method to PreparedGeometry 
• Investigate improving isWithinDistance by using the MINMAXDISTANCE metric for envelopes.  This allows earlier detection of index nodes satisfying the distance constraint.
• Investigate alternative R-Tree packing algorithms (such as Hilbert packing or sequence packing) to see if they improve performance

[1] Roussopoulos, Nick, Stephen Kelley, and Frédéric Vincent. "Nearest neighbor queries."  ACM SIGMOD record. Vol. 24. No. 2. ACM, 1995.

Joining Crunchy Data to work on PostGIS

I'm happy to announce that I am taking a position with Crunchy Data as a Senior Geospatial Engineer.  I'm working alongside fellow Victorian geospatial maven extraordinaire Paul Ramsey, as the core of a proposed Geospatial Data Centre of Excellence.

Our mission statement is simple:  make PostGIS bigger, better, faster!

• Bigger - more spatial algorithms and functions
• Better - enhance existing functionality to make it easier to use, more powerful, and more robust
• Faster - keep looking for algorithmic optimizations and ways to use the power of Postgres to make spatial processing faster

A lot of this work will involve enhancements to the core GEOS geometry library.  Part of the goal is to keep JTS and GEOS aligned, so this should produce a nice boost to JTS as well.

Having been lurking in the background for many years now, I'm stoked to be (finally) able to work directly on PostGIS.  And I'm excited to be part of the Crunchy team. They have some of the leading Postgres experts in-house, so I'm expecting that it will be a great learning experience.  And their client list promises to expose us to some fascinating large-scale use cases for spatial data processing, which can only be good for the power and robustness of PostGIS.

I'm also looking forward to re-engaging with the geospatial open source community, and learning more about the (even bigger) open source Postgres community.  Great things lie ahead!

The world needs a new flavour of SOSS!

Yes, you don't what the acronym SOSS means, because I just made it up. SOSS stands for Standard Open Simple Spatial format.

It's crazy that in the 21st century the most common de facto standard spatial format is based on 30-year old technology, is proprietary, and has silly limitations such as 11-character uppercase attribute names.

I'm talking, of course, about shapefiles.

Surely we can do better than this?!

Now there are actually a few things that shapefiles get right. For instance, the shapefile's simplistic tabular data model gets two full marks for being - simple and tabular! Hierarchical data models are very cool and highly expressive, but overkill and too complex for 80% of the use cases out there.

Another useful feature of shapefiles is that they store floating point data with full precision - i.e. in binary. Representing binary floating point numbers as textual decimal values is inherently lossy, and causes all kinds of subtle and annoying problems. (I'm always surprised that this doesn't crop up more often as a serious limitation of GML.)

So what are the current leading contenders for a SOSS format? Here's an opinionated list, with pros and cons

Format	Pro	Con
Shapefile	tabular, lossless numerics	proprietary, antiquated, limited
GML	complex to model and parse, lossy numerics, poor schema handling	open, flexible
KML	proprietary, lossy, limited attribute handling, designed for presentation	relatively simple, well documented
GeoRSS	not appropriate as a full-featured SOSS, lossy
GeoJSON	too tied to Javascript, lossy, no schema standard
YAML	needs a spatial profile

Conspicuous by its absence on this list is XML. In fact XML is a meta-format, not a format. To utilize XML would require defining an appropriate profile (which would need to be highly restricted to meet the criteria of simple). The major drawback of XML is that specifying the profile almost inevitably drags one into the mind-bending hell of XML Schema. (There are other schema languages, such as RelaxNG, but they involve similar complexity and have even less traction).

There's also more esoteric formats such as NetCDF, but it fails the simplicity test, and it's unclear how well it supports Geometry types.

SliceGraphs in JEQL

In the interests of increasing blog output, I'm going to experiment with using a terser style in posts where the content is mostly self-explanatory.  Sort of an embedded microblog...

Inspiration: original Population Lines print , this plot from the Line Graphs in R blog post:

Data: SEDAC World Population grid, subsampled

JEQL script

Plot:

Lead, the criminal element

I've heard before about the postulated link between atmospheric lead levels (courtesy of the leaded gasoline used through the middle decades of the 20th century) and crime levels.  This Mother Jones article America's Real Criminal Element: Lead is the best explanation I've seen so far (and has links to the original papers).  It really sounds like this hypothesis is fully confirmed - and the best thing about this story is that it has a happy ending.  (Unless you're trying to get elected as mayor - or Prime Minister - on a tough-on-crime platform).

There is a nice geospatial connection here.  As with many epidemiological issues, spatial locality is an important aspect of the analyses that lead (ahem) to the conclusion.  The article is chock-full of references to the spatial nature of the problem, such as:

We now have studies at the international level, the national level, the state level, the city level, and even the individual level

and my favourite:

a good rule of thumb for categorizing epidemics: If it spreads along lines of communication, he says, the cause is information. Think Bieber Fever. If it travels along major transportation routes, the cause is microbial. Think influenza. If it spreads out like a fan, the cause is an insect. Think malaria. But if it's everywhere, all at once—as both the rise of crime in the '60s and '70s and the fall of crime in the '90s seemed to be—the cause is a molecule.

YAUSEM (Yet Another US Election Map)

The usual US election map is a starkly simplistic collection of red and blue blobs.  This does reflect the esoteric (to a Canadian) Electoral College first-past-the-post arrangement.  But after seeing how close the actual vote counts were in most states, it seemed to me like this doesn't really reflect the actual political reality of the US.  Really it's pretty much a purely purple country.  John Nelson has a nice map that elegantly visualizes this reality at a county level (and Brian Timoney explains why it's the only map that isn't a lie).

As another attempt at map truthiness, I used JEQL to produce the following map.  It shows actual vote numbers at a state level, color themed along two dimensions:

• The hue shows the relative proportion of Democrat VS Republican votes (using the now-canonical blue and red).  For reference, Florida is almost exactly 50-50.  This nicely shows that really the US is just varying shades of purple.
• The saturation is proportional to the relative population of the state.  California is fully saturated, since it's the most populous state.  The the inland Western and the far Northeastern states are pretty pale, since they have fairly low populations.  This is roughly proportional to the weight of the state's Electoral College votes, although there are amusing anomalies.

I make no claim that this map represents any valid statistics - it's just a fun exercise in using JEQL to do spatial visualization.  For reference, the script is:

CSVReader t colSep: "\t" file: "us_vote_raw.txt";

t = select String.trim(col1) name,
    Val.toInt(String.keepChars(col2, "0123456789")) ecVote,
    Val.toDouble(String.keepChars(col4, "0123456789")) demvote,
    Val.toDouble(String.keepChars(col5, "0123456789")) repvote from t;

Print t;

maxVote = val(select max(demvote+repvote) from t);

ShapefileReader tus file: "us_state.shp";

tvote = select name, demvote, repvote, GEOMETRY,
        demvote / (demvote+repvote) demfrac,
        demvote+repvote totvote,        
       (demvote+repvote)/maxVote density
    from t
    join tus on t.name == String.trim(tus.STATE_NAME);

Mem tvote;
 
tplot = select GEOMETRY,
     styleFill,
    #ffffff styleStroke,
    1 styleStrokeWidth
    with {
        clr = Color.interpolate("ff0000", "0000ff", demfrac);
        h = Color.getH(clr);
        s = density;
        v = 1;
        styleFill = Color.toRGBfromHSV(h,s,v);    }   
    from tvote;
extent = BOX(-128 20, -65 50);

Plot    data: tplot
    extent: val(extent)
    width: 800
    height: 400
    background: "0000aa"
    file: "us_vote.png";

The raw data came from Wikipedia via simple cut-and-paste to a text file.

Battle of the Spatial Scripting Languages!

The gloves are down in this exchange on gis.stackexchange, over LOC counts for Fiona, JEQL, and GeoScript-JS.

Fiona takes one on the chin in the first round!  But Sean recovers quickly, and makes a hit by pointing out that Fiona allows defining new functions right in the scripting language.  This would be feasible in JEQL as well, by providing hooks to into JVM-based scripting languages.  And it's also very easy to develop extensions as simple Java functions or classes.

Now, this is really just a schoolyard scuffle over a very simple use case for spatial scripting.  What the fans really want to see is a face-off over a more complex task...

Barnes Analysis for Surface interpolation

Recently I've been working on generating surfaces from irregular sets of data observations, for thematic rendering purposes.  The data I'm working with is meteorological measurements (e..g max and min temperature, average rainfall, etc.).  Here's a sample:

Maximum Daily Temperatures for November 30, 2010

There are many, many different ways of interpolating surfaces - Wikipedia has a long list.  It turns out that a classic approach for interpolation of meteorological data is Barnes Interpolation (also called Barnes Analysis).  Barnes Interpolation is a surface estimating technique which uses an initial pass to provide a smoothed estimate of the surface, and optionally further refinement passes to improve the surface to more closely match the observation values.  The initial pass produces an averaged (smoothed) surface, using a summation of exponential (Gaussian) decay functions around each observation point.  Subsequent refinement passes compute an error surface using the delta between the previous estimated surface and the observations.  The error surface is added to the previous estimate surface to refine the estimate (by reducing the delta to the observations).

I'm more of a vector guy than a raster guy, but I figured that Google and Wikipedia would make short work of this - after all, how hard can raster be?  Unfortunately, the Wikipedia article is not very helpful - the description is unclear, and it doesn't even define all the symbols used in the presented equations.  I found some better references here and here, but they still contained some confusing discrepancies in the equations.  Barnes' original paper is online as well, but as is often the case with primary sources, it is more concerned with the physical properties of the analysis, and doesn't have a very crisp description of the actual algorithm.

Triangulating between the various sources, I eventually arrived at the following algorithm definition (provided herewith without warranty!):

Barnes Interpolation Algorithm For the first pass, the estimated value Eg at each grid point gx,y is:    Eg = ∑( wi * oi) / ∑( wi ) where: oi is the value of the i'th observation point wi is the weight (decay function) value for the i'th observation point, defined as       wi = exp(- di2   /   L2C  ) where: di is the distance from the grid point to the i'th observation point L is the length scale, which is determined by the observation spacing and the natural scale of the phenomena being measured.   C is the convergence factor, which controls how much refinement takes place during each refinement step.  In the first pass the convergence factor is 1. For subsequent passes a value in the range 0.2 - 0.3 is effective. During refinement passes the estimate at each grid point is re-computed using:   E'g = Eg + ∑( wi * (oi - Ei) ) / ∑( wi ) where: Ei is the estimated value at the grid cell containing the i'th observation point

The target development platform is Java (of course).  For prototyping purposes I used JEQL to read sample datasets, run the algorithm with various parameter choices, and produce plots to visualize the results.  Here's the surface plot for the dataset above:

This looks pretty reasonable, and provides a nice visualization of the large-scale data trend. (Surprise!  It's hotter down south...)  The ultimate goal is much more exciting.  The algorithm will be provided as a Rendering Transformation in GeoServer, allowing surfaces to be generated dynamically from any supported datastore.  This will open the way to providing all kinds of further surface visualizations, including other interpolation techniques - as well as everyone's favourite, Heat Maps (aka Multivariate Kernel Density Estimation).

Now, the pure mathematical description of the algorithm is sufficient for a proof-of-concept implementation producing single images.  But to make it usable in a dynamic way in a high-demand environment, there are more details which need to be worked out to improve performance and functionality.  For instance:

• How large an extent of data points is necessary to produce a stable surface for a given area?  This is important for dynamic presentation, since zooming, panning and tiling should not change the appearance of the surface.  Of course, all data points could always be used, but more selectivity improves performance.
• How far should the surface be extrapolated into areas of few or no data observations?  Apart from the validity implications, limiting the extent of the computed surface provides a boost in performance. 
• A further performance increase may be gained by computing the surface on a coarse grid, then up-sampling in an appropriate way to the display resolution.
• For data which is global in scope, how does the algorithm need to be changed to accommodate computation on the spheroid?  It seems likely that it's necessary to use a geodetic distance calculation, and there may be other impacts as well.

Let's do it again at FOSS4G-NA!

FOSS4G 2011 in Denver last September was such a good time, we're going to do it all over again at FOSS4G North America 2012! 

It will be interesting to see who and how many attend.  Some conference budgets might have been blown last year - but note that this event is cleverly timed to fall into a new fiscal!  And being in D.C. there should be lots of suits and spooks in attendance.

I'm definitely looking forward to building on the momentum from the last conference, especially given my new gig with OpenGeo (ok, I'm definitely not booth-babe material - but being a booth-geek is probably more fun). 

Your Favourite Map Projection

...and what it says about you.

XKCD goes geospatial!  (Happens to everything eventually...)

Data model diagrams for GTFS

Recently I've been doing some work with the General Transit Feed Specification (GTFS).  This is a simple data model initiated by Google, in order to allow transit organizations to publish their schedules for routing and visualization purposes.

As usual, Google has nice, succinct documentation for the GTFS format (which is a set of CSV files embodying an implicit data model).  But I was surprised to find that they don't provide any kind of data model diagram showing the relationships between the various entities (files) in the model.  I have a feeling that they think this would be "too technical" for users - or maybe they just don't believe in using a long-established body of theory?

So, to help my understanding of GTFS I put together a simple diagram of the data model.  I used Google Diagrams, which is why this looks a little primitive.   (Hmmm... maybe another reason they don't do diagrams is that they are religious about dog-fooding?)

After doing this I found the paper A Transmodel based XML schema for the Google Transit Feed Specification, which has a much more detailed UML model of GTFS:

The end goal of doing all this was to be able to transform a GTFS data feed into shapefile format.  For this I used a simple JEQL script.  But that's another blog post.