Fast detection of narrow polygons with JTS

Some spatial use cases require identifying "narrow" or "skinny" polygons.  A classic example is the process of cleaning polygonal coverages. Coverages (even clean ones) may contain narrow gaps between polygons which are unwanted.  But they may also contain larger gaps which are valid.  A geometric test is required to distinguish narrow polygons which are gaps (that can be merged) from wider ones (which should be retained). 

Narrowness Measures

There are various criteria which can be used to determine if a polygon is "narrow".  One approach is to use an estimate of circularity (also called compactness or roundness) such as the Polsby-Popper measure.  This is a measure based on the ratio of polygon area to the square of its perimeter.  Polygons with a measure closer to 1 are more circular than ones with lower values (which are thus more likely to be narrow).  This is easy to implement and fast to compute, but it has significant limitations:

• it is not obvious how to pick a measure value that determines "narrowness"
• it is based on shape-wide metrics (area and perimeter), so it may not accurately reflect the local properties of a shape (i.e. a polygon may have a low circularity measure and still contain a wide "bulge" in one place.)
• because it is a ratio it is size-independent, so very small gaps (which should be removed) may be compact according to the measure, and are thus retained. (Adding an area limit can mitigate this, but that introduces yet another difficult-to-choose parameter.)

Maximum Inscribed Circle isRadiusWithin

Really, the best way is to determine the narrowness of a polygon is to compute its actual width.  An effective way to do this is to compute the diameter of the largest circle which lies within the polygon.  This is known as the Maximum Inscribed Circle.  The JTS Topology Suite provides the MaximumInscribedCircle class to compute the centre point and radius of the required circle (see the blog post). (The algorithm is also available in downstream projects including GEOS, PostGIS and Shapely.)  The algorithm uses an iterative approximation, which currently requires specifying an approximation tolerance distance.

However, the real requirement is only to test whether the MIC radius is less than a given length value.  The actual value of the radius is not required.  Generally, an algorithm that tests for a value can be faster and more robust than one that computes the value explicitly.  This is a perfect example of this principle.  By formulating the function as a test against a provided radius length execution can be much faster, since computation can terminate as soon as it is known that the radius must be smaller or larger than the required length.  (Determining that the MIC radius is larger than the limit is easy: an internal point is found with distance to the polygon boundary which is greater than the limit.  The radius is known to be smaller than the limit once the grid size used to generate internal points is less than one-half the limit value.  If neither of those situations occurs the approximation tolerance is used to terminate the computation.)  This has recently been implemented as the MaximumInscribedCircle.isRadiusWithin predicate.  

A further benefit of the isRadiusWithin implementation is that it removes the need to specify the approximation tolerance.  The tolerance can be determined automatically, as a fraction of the specified radius limit.  This is significant, because selecting an appropriate distance tolerance for a Maximum Inscribed Circle computation is difficult, since it has no relationship to the size of the input geometry (i.e. a geometry may have a large extent, but still be very narrow.)  This is particularly challenging when a single tolerance must be chosen for a set of inputs of varying shape - precisely the use case for which this is being built.

Demonstration

As expected, evaluating the isRadiusWithin predicate is much faster than computing the actual value of the Maximum Inscribed Circle radius.  For a demonstration we use this dataset of world country borders (exploded to treat each island as a separate polygon).   This dataset has 8,393 polygons containing a total of 366,951 vertices.  

Computing the Maximum Inscribed Circle radius for each polygon, with a tolerance of 0.00001 degrees, takes 2.68 seconds.  

Evaluating the isRadiusWithin predicate, using a distance tolerance of 1 degree, executes in 0.034 seconds - almost 80 times faster.  

Usage in Coverage Cleaning

This function is used in the CoverageCleaner class to provide a fast test for narrow polygons (as specified by the user using the gapWidth parameter).  To see how this works on a real use case,  consider a polygonal coverage dataset of communes in Departement 71 in France (Saône-et-Loire).  The dataset has 564 polygons, with 612,238 vertices. As shown in a previous post on coverage cleaning, this dataset has highly non-aligned linework.  This creates lots of gaps.  The gaps can be found by cleaning the coverage with a gap tolerance distance of 0, which reveals 28,427 gaps, with 349,042 vertices.

Departement 71 communes showing the many gaps present

Close-up of gaps (note scale)

Computing the Maximum Inscribed Circle radius explicitly for each gap takes 25.4 seconds.  In contrast, using isRadiusWithin with a distance of 0.5 takes only 1.35 seconds - 18 times faster.  Since the cleaning process itself takes only 4.9 seconds, the reduction in time for the gap merging phase is significant.  

Future Work

It would be very nice to extend the automatic tolerance setting to the Maximum Inscribed Circle computation.  More research is required to see whether this can be done.

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.  

Polygonal Coverages and Operations in JTS

An important concept in spatial data modelling is that of a coverage.  A coverage models a two-dimensional region in which every point has a value out of a range (which may be defined over one or a set of attributes).  Coverages can be represented in both of the main physical spatial data models: raster and vector.  In the raster data model a coverage is represented by a grid of cells with varying values.  In the vector data model a coverage is a set of non-overlapping polygons (which usually, but not always, cover a contiguous area).  

This post is about the vector data coverage model, which is termed (naturally) a polygonal coverage. These are used to model regions which are occupied by discrete sub-regions with varying sets of attribute values.  The sub-regions are modelled by simple polygons.  The coverage may contain gaps between polygons, and may cover multiple disjoint areas.  The essential characteristics are:

• polygon interiors do not overlap
• the common boundary of adjacent polygons has the same set of vertices in both polygons.

There are many types of data which are naturally modelled by polygonal coverages.  Classic examples include:

• Man-made boundaries
• parcel fabrics
• political jurisdictions
• Natural boundaries
• vegetation cover
• land use

A polygonal coverage of regions of France

Topological and Discrete Polygonal Coverages

There are two ways to represent polygonal coverages: as a topological data structure, or as a set of discrete polygons.  

A coverage topology consists of linked faces, edges and nodes. The edges between two nodes form the shared boundary between two faces.  The coverage polygons can be reconstituted from the edges delimiting each face.  

The discrete polygon representation is simpler, and aligns better with the OGC Simple Features model.  It is simply a collection of polygons which satisfy the coverage validity criteria given above.

Most common spatial data formats support only a discrete polygon model, and many coverage datasets are provided in this form.  However, the lack of inherent topology means that datasets must be carefully constructed to ensure they have valid coverage topology.  In fact, many available datasets contain coverage invalidities.  A current focus of JTS development is to provide algorithms to detect this situation and provide the locations where the polygons fail to form a valid coverage.

Polygonal Coverage Operations

Operations which can be performed on polygonal coverages include:

• Validation - check that a set of discrete polygons forms a valid coverage
• Gap Detection - check if a polygonal coverage contains narrow gaps (using a given distance tolerance)
• Cleaning - fix errors such as gaps, overlaps and slivers in a polygonal dataset to ensure that it forms a clean, valid coverage
• Simplification - simplify (generalize) polygon boundary linework, ensuring coverage topology is preserved
• Precision Reduction - reduce precision of polygon coordinates, ensuring coverage topology is preserved
• Union - merge all or portions of the coverage polygons into a single polygon (or multipolygon, if the input contains disjoint regions)
• Overlay - compute the intersection of two coverages, producing a coverage of resultant polygons 

Implementing polygonal coverage operations is a current focus for development in the JTS Topology Suite.  Since most operations require a valid coverage as input, the first goal is to provide Coverage Validation.  Cleaning and Simplification are priority targets as well.  Coverage Union is already available, as is Overlay (in a slightly sub-optimal way).  In addition,  a Topology data structure will be provided to support the edge-node representation.  (Yes, the Topology Suite will provide topology at last!).  Stay tuned for further blog posts as functionality is rolled out.

As usual, the coverage algorithms developed in JTS will be ported to GEOS, and will thus be available to downstream projects like PostGIS.

Outer and Inner Polygon Hulls in JTS

The JTS Topology Suite recently gained the ability to compute concave hulls.  The Concave Hull algorithm computes a polygon enclosing a set of points using a parameter to determine the "tightness".  However, for polygonal inputs the computed concave hull is built only using the polygon vertices, and so does not always respect the polygon boundaries.  This means the concave hull may not contain the input polygon.

It would be useful to be able to compute the "outer hull" of a polygon.  This is a valid polygon formed by a subset of the vertices of the input polygon which fully contains the input polygon.  Vertices can be eliminated as long as the resulting boundary does not self-intersect, and does not cross into the original polygon.

An outer hull of a polygon representing Switzerland

As with point-set concave hulls, the vertex reduction is controlled by a numeric parameter. This creates a sequence of hulls of increasingly larger area with smaller vertex counts.  At an extreme value of the parameter, the outer hull is the same as the convex hull of the input.

A sequence of outer hulls of New Zealand's North Island

The outer hull concept extends to handle holes and MultiPolygons.  In all cases the hull boundaries are constructed so that they do not cross each other, thus ensuring the validity of the result.

An outer hull of a MultiPolygon for the coast of Denmark.  The hull polygons do not cross.

It's also possible to construct inner hulls of polygons, where the constructed hull is fully within the original polygon.

An inner hull of Switzerland

Inner hulls also support holes and MultiPolygons.  At an extreme value of the control parameter, holes become convex hulls, and a polygon shell reduces to a triangle (unless blocked by the presence of holes).

An inner hull of a lake with islands.  The island holes become convex hulls, and prevent the outer shell from reducing fully to a triangle

A hull can provide a significant reduction in the vertex size of a polygon for a minimal change in area. This could allow faster evaluation of spatial predicates, by pre-filtering with smaller hulls of polygons.

An outer hull of Brazil provides a 10x reduction in vertex size, with only ~1% change in area.

This has been on the JTS To-Do list for a while (I first proposed it back in 2009).  At that time it was presented as a way of simplifying polygonal geometry. Of course JTS has had the TopologyPreservingSimplifier for many years.  But it doesn't compute a strictly outer hull.  Also, it's based on Douglas-Peucker simplification, which isn't ideal for polygons.  

It seems there's quite a need for this functionality, as shown in these GIS-StackExchange posts (1, 2, 3, 4).  There's even existing implementations on Github: rdp-expansion-only and simplipy (both in Python) - but both of these sound like they have some significant issues. 

Recent JTS R&D (on concave hulls and polygon triangulation) has provided the basis for an effective, performant polygonal concave hull algorithm.  This is now released as the PolygonHullSimplifier class in JTS.

The PolygonHullSimplifier API

Polygon hulls have the following characteristics:

• Hulls can be constructed for Polygons and MultiPolygons, including holes.
• Hull geometries have the same structure as the input.  There is a one-to-one correspondence for  elements, shells and holes.
• Hulls are valid polygonal geometries.
• The hull vertices are a subset of the input vertices.

The PolygonHullSimplifier algorithm supports computing both outer and inner hulls. 

• Outer hulls contain the input geometry.  Vertices forming concave corners (convex for holes) are removed.  The maximum outer hull is the convex hull(s) of the input polygon(s), with holes reduced to triangles.
• Inner hulls are contained within the input geometry.  Vertices forming convex corners (concave for holes) are removed.   The minimum inner hull is a triangle contained in (each) polygon, with holes expanded to their convex hulls.  

The number of vertices removed is controlled by a numeric parameter.  Two different parameters are provided:

• the Vertex Number Fraction specifies the desired result vertex count as a fraction of the number of input vertices.  The value 1 produces the original geometry.  Smaller values produce simpler hulls.  The value 0 produces the maximum outer or minimum inner hull.
• the Area Delta Ratio specifies the desired maximum change in the ratio of the result area to the input area.  The value 0 produces the original geometry.  Larger values produce simpler hulls. 

Defining the parameters as ratios means they are independent of the size of the input geometry, and thus easier to specify for a range of inputs.  Both parameters are targets rather than absolutes; the validity constraint means the result hull may not attain the specified value in some cases. 

Algorithm Description

The algorithm removes vertices via "corner clipping".  Corners are triangles formed by three consecutive vertices in a (current) boundary ring of a polygon.  Corners are removed when they meet certain criteria.  For an outer hull, a corner can be removed if it is concave (for shell rings) or convex (for hole rings).  For an inner hull the removable corner orientations are reversed.  

In both variants, corners are removed only if the triangle they form does not contain other vertices of the (current) boundary rings.  This condition prevents self-intersections from occurring within or between rings. This ensures the resulting hull geometry is topologically valid.  Detecting triangle-vertex intersections is made performant by maintaining a spatial index on the vertices in the rings.  This is supported by an index structure called a VertexSequencePackedRtree.  This is a semi-static R-tree built on the list of vertices of each polygon boundary ring.  Vertex lists typically have a high degree of spatial coherency, so the constructed R-tree generally provides good space utilization.  It provides fast bounding-box search, and supports item removal (allowing the index to stay consistent as ring vertices are removed).

Corners that are candidates for removal are kept in a priority queue ordered by area.  Corners are removed in order of smallest area first.  This minimizes the amount of change for a given vertex count, and produces a better quality result.  Removing a corner may create new corners, which are inserted in the priority queue for processing.  Corners in the queue may be invalidated if one of the corner side vertices has previously been removed; invalid corners are discarded. 

This algorithm uses techniques originated for the Ear-Clipping approach used in the JTS PolgyonTriangulator implementation. It also has a similarity to the Visvalingham-Whyatt simplification algorithm.  But as far as I know using this approach for computing outer and inner hulls is novel. (After the fact I found a recent paper about a similar construction called a Shortcut Hull [Bonerath et al 2020], but it uses a different approach).

Further Work

It should be straightforward to use this same approach to implement a variant of Topology-Preserving Simplifier using the corner-area-removal approach (as in Visvalingham-Whyatt simplification).  The result would be a simplified, topologically-valid polygonal geometry.  The simplification parameter  limits the number of result vertices, or the net change in area.  The resulting shape would be a good approximation of the input, but is not necessarily be either wholly inside or outside.

JTS Offset Curves

Offset curves (also known as parallel curves) are an oft-requested feature in JTS.  They are a natural extension to the concept of buffering, and are useful for things like placing labels along rivers.  As far as I know there is no hard-and-fast definition for how an offset curve should be constructed, but a reasonable semantic would seem to be "a line lying on one side of another line at a given offset distance". 

Here's an image of offset curves on both sides of a river reach from a US rivers dataset:

GEOS Implementation

GEOS acquired an implementation to produce offset curves a few years ago.  However, it has some known bugs (such as producing disconnected curves, and including extraneous segments).  It also has the feature (or quirk?) of potentially producing a result with multiple disjoint curves for a single input line. 

The GEOS implementation is based on the concept that an offset curve is just a portion of the corresponding buffer boundary.  So to generate an offset curve the algorithm extracts a portion of the buffer polygon boundary.  The trick is deciding which portion!  

The GEOS implementation generates a raw offset curve (potentially full of self-intersections and unwanted linework) and then determines the intersection of that curve with the buffer boundary.  However, the use of intersection on linework is always tricky.  The buffer geometry is necessarily only a (close) approximation, and the buffer algorithm takes advantage of this to use various heuristics to improve the quality and robustness of the generated buffer linework.  This can cause the buffer linework to diverge from the raw offset curve.  The divergence makes the intersection result susceptible to errors caused by slight differences between the generated curves. The two issues above are caused by this limitation. 

JTS Implementation

Instead of using intersection, an effective technique is to match geometry linework using a distance tolerance.  This is the approach taken in the new JTS Offset Curve algorithm. The high-level design is

1. The buffer is generated at the offset distance
2. The raw offset curve is generated at the same distance
3. The raw curve segments are matched to the buffer boundary using a distance tolerance
4. The offset curve is the section of the buffer boundary between the first and last matching points.

To make the matching efficient a spatial index is created on the buffer curve segments.  

This algorithm provides the following semantics:

• The offset curve of a line is a single LineString
• The offset curve lies on one side of the input line (to the left if the offset distance is positive, and to the right if negative)
• The offset curve has the same direction as the input line
• The distance between the input line and the offset curve equals the offset distance (up to the limits of curve quantization and numerical precision)

This algorithm does a fine job at generating offset curves for typical simple linework.  The image below shows a family of offset curves on both sides of a convoluted line.  

This resolves both of the GEOS code issues.  It also supports parameters for join style (round, bevel, and mitre), number of quadrant segments (for round joins) and mitre limit:

Join Style = MITRE

Join Style = ROUND; Quadrant Segments = 2

There are also a few nuances required to handle some tricky situations; in particular, cases when the input line string curves back on itself and when it self-intersects.  These do not have a clear definition for what the result should be.  Moreover, the algorithm itself imposes some constraints on how these cases can be handled.  The images below show how the algorithm behaves in these cases.

"Loopbacks"produce offset curves representing only the "exposed" sections of the input linework: 

Offset Curves of a Line with "loop-backs"

For a self-intersecting input line, the offset curve starts at the beginning of the input line on the specified side, and continues only up to where the line side changes due to the crossing.  The length of the offset curve is also reduced by the requirement that it be no closer than specified offset distance to the input line: 

Offset Curves of a Line with a self-intersection

This algorithm is now in JTS, and has been ported to GEOS.

JTS Polygon Triangulation, at last

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!  

Polygon Triangulation

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: 

• the Hole Joining phase now uses a spatial index to optimize the test for whether a join line is internal to a polygon
• the Ear-Clipping phase has to check whether a candidate ear contains any vertices of the polygon shell being clipped.  This has been made performant by using a new spatial index structure called the VertexSequencePackedRtree

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

So JTS ear clipping is not quite as fast as Earcut4j.  But it's still plenty fast enough for production use!

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:

Ear-Clipping Triangulation of Lake Superior

Constrained Delaunay Triangulation

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.

Constrained Delaunay Triangulation of Lake Superior

This code will appear in the next JTS version (1.18.3 or 1.19), and has already been released in GEOS 3.10.

Future Work

It seems possible to adapt the machinery for Ear Clipping to provide two other useful constructions:

• Polygon Concave Hull ("Outer Hull") - this requires a triangulation of the area between the polygon and its convex hull.  That can be obtained by running Ear Clipping outwards, rather than inwards. A Concave Hull can be constructed by limiting the number of triangles generated based on some criteria such as area, longest edge, or number of vertices.  The result is a simplified polygon with fewer vertices which is guaranteed to contain the original polygon. 
• Polygon Inner Hull - the triangulation constructed by Ear Clipping can be "eroded" to provide a simpler polygon which is guaranteed to be contained in the original.

The property of the Constrained Delaunay Triangulation of representing the structure of the input polygon could provide a basis for computing the following constructions: 

• Approximate Skeleton construction (also known as the Medial Axis) by creating edges through the CDT triangles
• Equal-Area Polygon Subdivision, by partitioning the graph induced by the CDT

The Tri data structure provides a simple representation of a triangulation.  I expect that it will facilitate the development of further JTS triangulation algorithms, such as:

• Concave Hull of Point Set
• Constrained Delaunay Triangulation of Lines and Points

Using the GEOS geosop CLI

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.)

API Testing

geosop makes testing the GEOS library much easier.  Previously, testing the behaviour of the API usually required breaking out the C compiler (and updating autotools and cmake build files, and deciding whether to commit the new code for later use or throw it away, etc, etc).  Now testing is often just a matter of invoking a geosop operation on appropriate data, or at worst adding a few lines of code to the exiting framework.

For example, there is a long-standing issue with how GEOS handles number formatting in WKT output.  There are recent bug reports about this in GeoSwift and the Julia LibGEOS.  geosop makes it easy to run the test cases and see the less-than-desirable output:

geosop -a "POINT (654321.12 0.12)" -f wkt

POINT (654321.1199999999953434 0.1200000000000000)

geosop -a "POINT (-0.4225977234 46.3406448)" -f wkt

POINT (-0.4225977234000000 46.3406447999999997)

There's also an issue with precision handling.  To test this we added a --precision parameter to geosop. (This is the kind of rapid development enabled by having the CLI codebase co-resident with the API.)

geosop -a "POINT (654321.126 0.126)" --precision 2 -f wkt

POINT (6.5e+05 0.13)

Again we see undesirable behaviour. Using scientific notation for small numbers is unnecessary and difficult to read. And the precision value is determining the number of significant digits, not the number of decimal places as intended by the GEOS WKTWriter.setRoundingPrecision API.

These were all caused by using standard C/C++ numeric formatting, which is surprisingly limited and non-useful. After some fine work by Paul Ramsey to integrate the much better Ryu library, GEOS now has WKT output that is sensible and handles precision in a useful way.

By default, WKTWriter now nicely round-trips WKT text:

geosop -a "POINT (654321.12 0.12)" -f wkt

POINT (654321.12 0.12)

geosop -a "POINT (-0.4225977234 46.3406448)" -f wkt

POINT (-0.4225977234 46.3406448)

If WKTWriter.setRoundingPrecision or GEOSWKTWriter_setRoundingPrecision is called, the precision value applies to the decimal part of the number: 

geosop -a "POINT (654321.1234567 0.126)" --precision 0 -f wkt

POINT (654321 0)

geosop -a "POINT (654321.1234567 0.126)" --precision 2 -f wkt

POINT (654321.12 0.13)

geosop -a "POINT (654321.1234567 0.126)" --precision 4 -f wkt

POINT (654321.1235 0.126)

geosop -a "POINT (654321.1234567 0.126)" --precision 6 -f wkt

POINT (654321.123457 0.126)

Performance Testing

A key use case for geosop is to provide easy performance testing.  Performance of geometric operations is highly data-dependent.  It's useful to be able to run operations over different datasets and measure performance.  This allows detecting performance hotspots, and confirming the efficiency of algorithms.

As a simple example of performance testing, for many years GEOS has provided optimized spatial predicates using the concept of a prepared geometry.  Prepared geometry uses cached spatial indexes to dramatically improve performance for repeated spatial operations against a geometry.  Here is a performance comparison of the intersects spatial predicate in its basic and prepared form. 

geosop -a world.wkt -b world.wkt -t intersects

Ran 59,536 intersects ops ( 179,072,088 vertices)

  -- 16,726,188 usec    (GEOS 3.10.0dev)

geosop -a world.wkt -b world.wkt -t intersectsPrep

Ran 59,536 intersectsPrep ops ( 179,072,088 vertices)

  -- 1,278,348 usec    (GEOS 3.10.0dev)

The example of testing the world countries dataset against itself is artificial, but it shows off the dramatic 16x performance boost provided by using a prepared operation. 

Another interesting use is longitudinal testing of GEOS performance across different versions of the library.  For instance, here's a comparison of the performance of interiorPoint between GEOS 3.7 and 3.8.  The interiorPoint algorithm in GEOS 3.7 relied on overlay operations, which made it slow and sensitive to geometry invalidity. GEOS 3.8 included a major improvement which greatly improved the performance, and made the algorithm more robust.  Note that the dataset being used contains some (mildly) invalid geometries towards its end, which produces an error in GEOS 3.7 if the entire dataset is run.  The  --alimit 3800 option limits the number of geometries processed to avoid this issue.

GEOS 3.7

geosop -a ne_10m_admin_1_states_provinces.wkb --alimit 3800 -t interiorPoint

Ran 3,800 operations ( 1,154,703 vertices)  -- 2,452,540 usec  (GEOS 3.7)

GEOS 3.8

geosop -a ne_10m_admin_1_states_provinces.wkb --alimit 3800 -t interiorPoint

Ran 3,800 operations ( 1,154,703 vertices)  -- 35,665 usec  (GEOS 3.8)

The dramatic improvement in interiorPoint performance is clearly visible.

Geoprocessing

The raison d'etre of GEOS is to carry out geoprocessing.  Most users will likely do this using one of the numerous languages, applications and databases that include GEOS.  But it is still useful, convenient, and perhaps more performant to be able to process geometric data natively in GEOS.  

The design of geosop enables more complex geoprocessing via chaining operations together using shell pipes.  For example, here is a process which creates a Voronoi diagram of some points located in the British Isles, and then clips the Voronoi polygons to the outline of the islands.  This also shows a few more capabilities of geosop:

• input can be supplied as WKT (or WKB) geometry literals on the command-line
• input can be read from the standard input (here as WKB)
• the output data is sent to the standard output, so can be directed into a file

geosop 

  -a "MULTIPOINT ((1342 1227.5), (1312 1246.5), (1330 1270), (1316.5 1306.5), (1301 1323), (1298.5 1356), (1247.5 1288.5), (1237 1260))" 

  -f wkb voronoi 

| geosop -a stdin.wkb -b uk.wkt -f wkt intersection 

> uk-vor.wkt

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.