SpatiaLite: Update of "liblwgeom-4.0"

Many hyperlinks are disabled.
Use anonymous login to enable hyperlinks.

Overview

Artifact ID:	66eb25ca81a89c781b1cabf827ef03811c79cbcf
Page Name:	liblwgeom-4.0
Date:	2012-08-31 10:50:08
Original User:	sandro
Parent:	ee415825b4dcbdb2bb73a971e23173fd7a7d703a (diff)
Next	974b02eb4d86705ac227f81a0d0274309fb3df95

Content

SQL functions based on liblwgeom support in version 4.0.0

Back to main page

Making invalid Geometries to become perfectly valid ones

As you already surely know, not all geometries are valid ones: most notably in the specific case of Polygonal geometries there are many formal rules to be carefully respected. Infringing one of such rules directly leads to some invalid Geometry: and an invalid Geometry could cause invalid results to be returned, or could eventually cause some unexpected nasty crash in the worst case.

GEOS (and thus SpatiaLite) already supports the ST_IsValid() SQL function; by invoking this function you can easily identify all offending Geometries eventually contained within your tables; anyway you cannot attempt to sanitize them.
SpatiaLite already supported a ST_SanitizeGeometry() SQL function; but this was simply capable to effectively resolve just few invalidity causes, and wasn't at all a general solution for this problem.

Now, thanks to liblwgeom, SpatiaLite can support the same identical ST_MakeValid() already supported by PostGIS; few small implementation details differ (due to the huge architectural differences distinguishing PostGIS and SpatiaLite), but the underlying code is exactly one and the same for both.

a first basic example

We'll start loading the Local Councils administrative boundaries supplied by ISTAT (the Italian National Statics Agency): this dataset is freely available for download under a CC-BY license. Just a quick check, and we'll soon discover that this dataset actually contains several malformed Polygons:

SELECT Count(*)
FROM com2011
WHERE ST_IsValid(geometry) = 0;

Recovering all malformed Geometries is now absolutely simple and easy:

UPDATE com2011 SET geometry = ST_MakeValid(geometry)
WHERE ST_IsValid(geometry) = 0;

SELECT Count(*)
FROM com2011
WHERE ST_IsValid(geometry) = 0;

Why the Bronte Local Council boundary was malformed ? a quick analysis
malformed / invalid	valid

As you can easily notice, the invalid polygon was simply represent by the exterior ring: but there is a huge internal hole in this Polygon. This odd condition is reputed perfectly valid by some mainstream proprietary software; anyway, is actually invalid accordingly to standard rules. So the correct representation for this Polygon requires an exterior ring and a separate interior ring; ST_MakeValid() does the magic, thus recovering a full valid Polygon.

a second more elaborate example

This time we'll purposely create a severely malformed Polygon:

SELECT g, ST_MakeValid(g), ST_MakeValidDiscarded(g)
FROM (
SELECT ST_GeomFromText('POLYGON((0 0, 0 10, 11 10, 10 10, 10 1, 5 1, 5 9, 5 1, 0 1, 0 0))') AS g
);

This figure represents the invalid Polygon returned by ST_GeomFromText(). There are three spikes in the exterior ring, and this one surely is a severe invalidity cause.
This figure represents the valid Polygon returned by ST_MakeValid(). Now we have a nice regular rectangle, all spikes have been removed.
Anyway the removed spikes aren't simply vanished into nothing. You could eventually retrieve (and may be, saving somewhere for further processing / editing) all offending elements being discarded during the validation process. You simply have to invoke ST_MakeValidDiscarded() Please note: this is strongly different from the PostGIS own implementation.

Splitting geometries in two halves

The ST_Split() Spatial SQL function is intended to cut a Geometry.
This function always requires to pass two different Geometries:

the first Geometry is assumed to represent the input aka target to be split.
the second Geometry is assumed to represent the blade.

Only the following configurations are assumed to be valid:

Linestring target / Point blade.
MultiLinestring target / Point blade.
GeometryCollection (containing at least one Linestring) / Point blade.
Linestring target / Linestring blade.
MultiLinestring target / Linestring blade.
Polygon target / Linestring blade.
MultiPolygon target / Linestring blade.
GeometryCollection / Linestring blade.
Please notice: the target Geometry should never contain a Point.

Example #1 - splitting a Linestring by a Point blade

SELECT ST_Split(input.g, blade.g), ST_SplitLeft(input.g, blade.g), ST_SplitRight(input.g, blade.g)
FROM
(SELECT GeomFromText('LINESTRING(0 10, 2 0, 4 4, 6 0, 10 10)') AS g) AS input,
(SELECT GeomFromText('POINT(3 2)') AS g) AS blade;

The ST_Split() SQL function will simply return a collection aggregating all fragments deriving from the cut. This isn't really useful on many cases. Please note: collections in SpatiaLite behave quite differently from PostGis.
The ST_SplitLeft() SQL function will return instead a collection aggregating all fragments laying on the left side of the cut. For Linestrings you cannot intend left in the very literal sense; this really means the side where the start-point lay. Please note: if the blade doesn't intercepts the target at all, than no cut would be obviously possible. In this special case the original target Geometry will be always returned (absolutely unchanged) on the left side collection.
The ST_SplitRight() SQL function will return a collection aggregating all fragments laying on the right side of the cut. For Linestrings you cannot intend right in the very literal sense; this really means the side where the end-point lay. Please note: if the blade doesn't intercepts the target at all, than no cut would be obviously possible. In this special case the right side collection will be NULL.

Example #2 - splitting a Polygon by a Linestring blade

SELECT ST_Split(input.g, blade.g), ST_SplitLeft(input.g, blade.g), ST_SplitRight(input.g, blade.g)
FROM
(SELECT GeomFromText('POLYGON((0 1, 10 1, 10 9, 0 9, 0 1), (2 2, 7 2, 7 6, 2 6, 2 2))') AS g) AS input,
(SELECT GeomFromText('LINESTRING(2 0, 6 10)') AS g) AS blade;

The ST_Split() SQL function will simply return a collection aggregating all fragments deriving from the cut. Exactly as we have already previously seen on the Linestring case.
The ST_SplitLeft() SQL function will return instead a collection aggregating all fragments laying on the left side of the cut. For Polygons left really means left side (at least, this is true when the blade is almost vertically oriented). Please note: if the blade doesn't intercepts the target at all, than no cut would be obviously possible. In this special case the original target Geometry will be always returned (absolutely unchanged) on the left side collection.
The ST_SplitRight() SQL function will return a collection aggregating all fragments laying on the right side of the cut. Please note: if the blade doesn't intercepts the target at all, than no cut would be obviously possible. In this special case the right side collection will be NULL.

Segmentization

Sometimes it could be useful interpolating many further vertices into some very long Linestring or Polygon's Ring, so to obtain many smaller segments of the same identical length.

SELECT geometry, ST_Segmentize(geometry, 10.0)
FROM com2011
WHERE nome_com = 'Bronte';

original	segmentized

As you can easily notice, there are many more vertices in the segmentized Geometry. We imposed the constraint that no segment could be longer than 10m, thus causing many more vertices to be interpolated. The overall shape is absolutely unchanged.

Azimuth

Probably not the most interesting function supported by liblwgeom; anyway it's there, and is worth enough to be supported by SpatiaLite as well.
Probably highly appreciated on GPS devices so to determine the current bearing.

SELECT Degrees(ST_Azimuth(MakePoint(0, 0), MakePoint(0, 1)));

0.0

SELECT Degrees(ST_Azimuth(MakePoint(0, 0), MakePoint(1, 1)));

45.0

SELECT Degrees(ST_Azimuth(MakePoint(0, 0), MakePoint(1, 0)));

90.0

SELECT Degrees(ST_Azimuth(MakePoint(0, 0), MakePoint(0, -1)));

180.0

Snapping Geometries to a predefined grid

This feature really is the mother of topological consistency; imposing that all points / vertices should be strictly aligned to the same grid will surely make easier ensuring a strong consistency between different Geometries.

Please note: due to technical reasons the current version of liblwgeom doesn't exposes to the exterior the corresponding link symbols.
So the current implementation available on SpatiaLite is an independent one (although being directly inspired by the liblwgeom one).

SELECT ST_SnapToGrid(geometry, 250.0)
FROM com2011
WHERE nome_com = 'Bronte';

Yet again another time the Bronte Local Council boundary.
By invoking ST_SnapToGrid() this time the Geometry has been strongly simplified (please notice; many vertices have been suppressed).

The imposed grid size (250m) is absolutely unrealistic and exaggerated for any practical purpose.
Anyway a such strong factor is useful in order to clearly show what really happens:

all vertices are now strictly aligned to the grid nodes.
all vertices falling within the same grid cell have been suppressed, still preserving only one of them.
so the overall effect is the one leading to a rigorously normalized and (may be) simplified Geometry; and this will surely make easier ensuring the required topological consistency with any other adjacent Local Council.

Credits
This work (exposing liblwgeom APIs as SpatiaLite own SQL functions) has been entirely funded by: Tuscany Region - Territorial and Environmental Information System Regione Toscana - Settore Sistema Informativo Territoriale ed Ambientale.

Back to main page