D 2014-08-27T13:00:34.949
L liblwgeom-4.0
P 974b02eb4d86705ac227f81a0d0274309fb3df95
U sandro
W 14346
<h2>SQL functions based on liblwgeom support in version 4.0.0</h2>
<a href="https://www.gaia-gis.it/fossil/libspatialite/wiki?name=4.2.0-doc">back</a>
<h3>Making invalid Geometries to become perfectly valid ones</h3>
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.<br><br>
GEOS (and thus SpatiaLite) already supports the <b>ST_IsValid()</b> SQL function; by invoking this function you can easily identify all offending Geometries eventually contained within your tables; anyway you cannot attempt to <i>sanitize</i> them.<br>
SpatiaLite already supported a <b>ST_SanitizeGeometry()</b> 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.<br><br>
Now, thanks to <i>liblwgeom</i>, SpatiaLite can support the same identical <b>ST_MakeValid()</b> 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.
<h4>a first basic example</h4>
We'll start loading the Local Councils administrative boundaries supplied by ISTAT (the Italian National Statics Agency): this dataset is freely available for <a href="http://www.istat.it/uploads/com2011.zip">download</a> under a CC-BY license.
Just a quick check, and we'll soon discover that this dataset actually contains several malformed Polygons:
<table bgcolor="#e0e0e0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT Count(*)<br>
FROM com2011<br>
WHERE ST_IsValid(geometry) = 0;</b>
<hr>
19
</td></tr>
</table><br>
Recovering all malformed Geometries is now absolutely simple and easy:
<table bgcolor="#e0e0e0" cellspacing="4" cellpadding="8"><tr><td>
<b>UPDATE com2011 SET geometry = ST_MakeValid(geometry)<br>
WHERE ST_IsValid(geometry) = 0;</b>
</td></tr>
<tr><td>
<b>SELECT Count(*)<br>
FROM com2011<br>
WHERE ST_IsValid(geometry) = 0;</b>
<hr>
0
</td></tr>
</table><br>
<table border="1" cellspacing="4" cellpadding="4">
<tr><th colspan="2">Why the Bronte Local Council boundary was malformed ?<br>
a quick analysis</th></tr>
<tr><th>malformed / invalid</th><th>valid</th></tr>
<tr><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/bronte-invalid.png" alt="bronte invalid">
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/bronte-valid.png" alt="bronte valid">
</td></tr>
<tr><td colspan="2">
As you can easily notice, the <i>invalid</i> polygon was simply represent by the <i>exterior ring</i>: 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 <i>invalid</i> accordingly to standard rules.<br>
So the correct representation for this Polygon requires an <i>exterior ring</i> and a separate <i>interior ring</i>; ST_MakeValid() does the magic, thus recovering a full valid Polygon.
</td></tr>
</table>
<h4>a second more elaborate example</h4>
This time we'll purposely create a severely malformed Polygon:
<table bgcolor="#e0e0e0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT g, ST_MakeValid(g), ST_MakeValidDiscarded(g)<br>
FROM (<br>
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<br>
);</b>
</td></tr>
</table><br>
<table cellspacing="4" cellpadding="4">
<tr><td>
This figure represents the <i>invalid</i> Polygon returned by <b>ST_GeomFromText()</b>.<br><br>
There are three <i>spikes</i> in the <i>exterior ring</i>, and this one surely is a severe invalidity cause.
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/spike-invalid.png" alt="spike invalid">
</td></tr>
<tr><td>
This figure represents the <i>valid</i> Polygon returned by <b>ST_MakeValid()</b>.<br><br>
Now we have a nice regular rectangle, all <i>spikes</i> have been removed.
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/spike-valid.png" alt="spike valid">
</td></tr>
<tr><td>
Anyway the removed <i>spikes</i> 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.<br><br>
You simply have to invoke <b>ST_MakeValidDiscarded()</b><br><br>
<u>Please note</u>: this is strongly different from the PostGIS own implementation.
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/spike-discarded.png" alt="spike discarded">
</td></tr>
</table>
<br><hr>
<h3>Splitting geometries in two halves</h3>
The <b>ST_Split()</b> Spatial SQL function is intended to <i>cut</i> a Geometry.<br>
This function always requires to pass two different Geometries:<ul>
<li>the first Geometry is assumed to represent the <b>input</b> aka <b>target</b> to be split.</li>
<li>the second Geometry is assumed to represent the <b>blade</b>.</li>
</ul>
Only the following configurations are assumed to be valid:<ul>
<li>Linestring target / Point blade.</li>
<li>MultiLinestring target / Point blade.</li>
<li>GeometryCollection (containing at least one Linestring) / Point blade.</li>
<li>Linestring target / Linestring blade.</li>
<li>MultiLinestring target / Linestring blade.</li>
<li>Polygon target / Linestring blade.</li>
<li>MultiPolygon target / Linestring blade.</li>
<li>GeometryCollection / Linestring blade.</li>
<li>Please notice: the target Geometry should never contain a Point.</li>
</ul>
<h4>Example #1 - splitting a Linestring by a Point blade</h4>
<table bgcolor="#e0e0e0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT ST_Split(input.g, blade.g), ST_SplitLeft(input.g, blade.g), ST_SplitRight(input.g, blade.g)<br>
FROM<br>
(SELECT GeomFromText('LINESTRING(0 10, 2 0, 4 4, 6 0, 10 10)') AS g) AS input,<br>
(SELECT GeomFromText('POINT(3 2)') AS g) AS blade;</b>
</td></tr></table>
<table cellspacing="4" cellpadding="4">
<tr><td>
The <b>ST_Split()</b> SQL function will simply return a <i>collection</i> aggregating all fragments deriving from the cut.
This isn't really useful on many cases.<br><br>
<u>Please note</u>: <i>collections</i> in SpatiaLite behave quite differently from PostGis.
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/split-ln.png" alt="split line">
</td></tr>
<tr><td>
The <b>ST_SplitLeft()</b> SQL function will return instead a <i>collection</i> aggregating all fragments laying on the <i>left</i> side of the cut.<br>
For Linestrings you cannot intend <i>left</i> in the very literal sense; this really means <i>the side where the start-point lay</i>.<br><br>
<u>Please note</u>: if the <i>blade</i> 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 <i>left</i> side collection.
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/split-ln-L.png" alt="split-left line">
</td></tr>
<tr><td>
The <b>ST_SplitRight()</b> SQL function will return a <i>collection</i> aggregating all fragments laying on the <i>right</i> side of the cut.<br>
For Linestrings you cannot intend <i>right</i> in the very literal sense; this really means <i>the side where the end-point lay</i>.<br><br>
<u>Please note</u>: if the <i>blade</i> doesn't intercepts the target at all, than no cut would be obviously possible. In this special case the <i>right</i> side collection will be NULL.
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/split-ln-R.png" alt="split-right line">
</td></tr>
</table>
<h4>Example #2 - splitting a Polygon by a Linestring blade</h4>
<table bgcolor="#e0e0e0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT ST_Split(input.g, blade.g), ST_SplitLeft(input.g, blade.g), ST_SplitRight(input.g, blade.g)<br>
FROM<br>
(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,<br>
(SELECT GeomFromText('LINESTRING(2 0, 6 10)') AS g) AS blade;</b>
</td></tr>
</table><br>
<table cellspacing="4" cellpadding="4">
<tr><td>
The <b>ST_Split()</b> SQL function will simply return a <i>collection</i> aggregating all fragments deriving from the cut.<br>
Exactly as we have already previously seen on the Linestring case.<br><br>
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/split-pg.png" alt="split polygon">
</td></tr>
<tr><td>
The <b>ST_SplitLeft()</b> SQL function will return instead a <i>collection</i> aggregating all fragments laying on the <i>left</i> side of the cut.<br>
For Polygons <i>left</i> really means <i>left side</i> (at least, this is true when the blade is almost vertically oriented).<br><br>
<u>Please note</u>: if the <i>blade</i> 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 <i>left</i> side collection.
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/split-pg-L.png" alt="split-left polygon">
</td></tr>
<tr><td>
The <b>ST_SplitRight()</b> SQL function will return a <i>collection</i> aggregating all fragments laying on the <i>right</i> side of the cut.<br>
<u>Please note</u>: if the <i>blade</i> doesn't intercepts the target at all, than no cut would be obviously possible. In this special case the <i>right</i> side collection will be NULL.
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/split-pg-R.png" alt="split-right polygon">
</td></tr>
</table>
<br><hr>
<h3>Segmentization</h3>
Sometimes it could be useful <i>interpolating</i> many further vertices into some very long Linestring or Polygon's Ring, so to obtain many smaller segments of the same identical length.
<table bgcolor="#e0e0e0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT geometry, ST_Segmentize(geometry, 10.0)<br>
FROM com2011<br>
WHERE nome_com = 'Bronte';</b>
</td></tr></table><br>
<table border="1" cellspacing="4" cellpadding="4">
<tr><th>original</th><th>segmentized</th></tr>
<tr><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/bronte-original.png" alt="bronte original">
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/bronte-segmentized.png" alt="bronte segmentized">
</td></tr>
<tr><td colspan="2">
As you can easily notice, there are many more vertices in the segmentized Geometry.<br>
We imposed the constraint that no segment could be longer than <b>10m</b>, thus causing many more vertices to be interpolated.<br>
The overall shape is absolutely unchanged.
</td></tr>
</table>
<br><hr>
<h3>Azimuth</h3>
Probably not the most interesting function supported by <i>liblwgeom</i>; anyway it's there, and is worth enough to be supported by SpatiaLite as well.<br>
Probably highly appreciated on <b>GPS devices</b> so to determine the current bearing.
<table bgcolor="#e0e0e0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT Degrees(ST_Azimuth(MakePoint(0, 0), MakePoint(0, 1)));</b>
<hr>
0.0
</td></tr>
<tr><td>
<b>SELECT Degrees(ST_Azimuth(MakePoint(0, 0), MakePoint(1, 1)));</b>
<hr>
45.0
</td></tr>
<tr><td>
<b>SELECT Degrees(ST_Azimuth(MakePoint(0, 0), MakePoint(1, 0)));</b>
<hr>
90.0
</td></tr>
<tr><td>
<b>SELECT Degrees(ST_Azimuth(MakePoint(0, 0), MakePoint(0, -1)));</b>
<hr>
180.0
</td></tr>
</table><br>
<br><hr>
<h3>Snapping Geometries to a predefined grid</h3>
This feature really is <i><u>the mother of topological consistency</u></i>; imposing that all points / vertices should be strictly aligned to the same grid will surely make easier ensuring a strong consistency between different Geometries.<br><br>
<u>Please note</u>: due to technical reasons the current version of <i>liblwgeom</i> doesn't exposes to the exterior the corresponding link symbols.<br>
So the current implementation available on SpatiaLite is an independent one (although being directly inspired by the <i>liblwgeom</i> one).<br>
<table bgcolor="#e0e0e0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT ST_SnapToGrid(geometry, 250.0)<br>
FROM com2011<br>
WHERE nome_com = 'Bronte';</b>
</td></tr></table><br>
<table cellspacing="4" cellpadding="4">
<tr><td>
Yet again another time the Bronte Local Council boundary.<br>
By invoking <b>ST_SnapToGrid()</b> this time the Geometry has been strongly simplified (please notice; many vertices have been suppressed).<br><br>
The imposed <i>grid size</i> (<b>250m</b>) is absolutely unrealistic and exaggerated for any practical purpose.<br>
Anyway a such strong factor is useful in order to clearly show what really happens:<ul>
<li>all vertices are now strictly aligned to the grid nodes.</li>
<li>all vertices falling within the same grid cell have been suppressed, still preserving only one of them.</li>
<li>so the overall effect is the one leading to a rigorously <i>normalized</i> and (may be) <i>simplified</i> Geometry; and this will surely make easier ensuring the required topological consistency with any other adjacent Local Council.</li>
</ul>
</td><td>
<img src="http://www.gaia-gis.it/gaia-sins/write-view-pics/snap-to-grid.png" alt="snap to grid">
</td></tr>
</table>
<br><hr>
<table bgcolor="#d7ffff" cellspacing="6" cellpadding="10" width="100%">
<tr><th>Credits</th></tr>
<tr><td align="center">
This work (exposing <i>liblwgeom</i> APIs as SpatiaLite own SQL functions) has been entirely funded by:<br>
<a href="http://en.wikipedia.org/wiki/Tuscany">Tuscany Region</a> - Territorial and Environmental Information System<br>
<a href="http://www.regione.toscana.it/ambienteeterritorio/geografiageologia/index.html">Regione Toscana</a> - Settore Sistema Informativo Territoriale ed Ambientale.
</td></tr>
</table>
<br><hr>
<a href="https://www.gaia-gis.it/fossil/libspatialite/wiki?name=4.2.0-doc">main page</a>
Z bb2f74fefe4472c606bd1d189ef3b5d4