D 2012-08-30T17:15:23.443
L liblwgeom-4.0
U sandro
W 11913
<h2>SQL functions based on liblwgeom support in version 4.0.0</h2>
Back to <a href="https://www.gaia-gis.it/fossil/libspatialite/wiki?name=switching-to-4.0#liblwgeom">main page</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 detailed rules to be 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 support 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 exactely 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 datasat 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) 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 - splltting 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.<br>
This isn't really useful.<br><br>
<u>Please note</u>: <i>collections</i> in SpatiaLite behave quite diffently 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 - splltting 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>
Exectly 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>cutting</i> a very long Linestring or Polygon's Ring into many small 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>
Back to <a href="https://www.gaia-gis.it/fossil/libspatialite/wiki?name=switching-to-4.0#liblwgeom">main page</a>
Z c6bee7476e3d0032dced7d41efcf5e98