D 2014-08-27T13:02:24.187
L tesselations-4.0
P abe370e2a435d8480c00917a9d29a310aa7ab72a
U sandro
W 20062
<h2>Tessellation-related SQL functions supported in version 4.0.0</h2>
<a href="https://www.gaia-gis.it/fossil/libspatialite/wiki?name=4.2.0-doc">back</a>
<h3>Generalities about tessellations</h3>
<i>Tessellation is the process of creating a two-dimensional plane using the repetition of a geometric shape with no overlaps and no gaps.</i>
<a href="http://en.wikipedia.org/wiki/Tessellation">read more</a><br><br>
A tessellation could be eventually based on identical <i>cells</i>, all of exactly the same shape and size.
In this case we'll have a <i>regular tessellation</i>.<br>
Just a quick recall of elementary geometry; there are simply <b>3</b> regular polygonal shapes we can use in order to get a regular tessellation: the <b>equilateral triangle</b>, the <b>square</b> and the <b>regular hexagon</b>.<br><br>
On the other way many tessellations aren't regular at all, because each single cell has an individual size and shape of its own.<br>
It's really interesting to note that many <i>natural shapes</i> closely resemble a tessellation: going from biology to crystallography since geology and landscapes it's not at all difficult to identify many natural tessellation examples <i>on the wild</i>.
So it's not at all surprising to discover that tessellations are often useful in geography as well.
<h3>Setting up a testbed DB</h3>
In this short tutorial we'll use a very simple SpatiaLite DB, just containing the following Geometry tables:<ul>
<li>Administrative boundaries for Italy regions: you can <a href="http://www.istat.it/it/files/2011/04/reg2011.zip">download</a> this dataset from ISTAT (released on CC-BY license terms).<br></li>
<li>Main Italian populated places (namely, Local Councils): you can <a href="http://download.geonames.org/export/dump/allCountries.zip">download</a> this dataset from GeoNames (released on CC-BY license terms).<br>
Please note: this one is a worldwide dataset; italian populated places have then been extracted imposing the SQL clause<br>
<b>WHERE county_code = 'IT' AND population > 0</b></li>
</ul><br>
Just to keep any example as simple as possible, both datasets have been referenced into the <b>SRID=23032 - ED50 / UTM zone 32 N</b>
<br><br><hr>
<h3>Creating regular tessellations in Spatial SQL</h3>
<table cellspacing="4" cellpadding="4">
<tr><td colspan="2">
<table bgcolor="#f0f0f0" cellspacing="4" cellpadding="8">
<tr><td>
<i>generic prototypes</i>:<hr>
ST_xxxxGrid( input <i>Geometry</i>, size <i>double precision</i> ) : grid <i>Geometry</i><br>
ST_xxxxGrid( input <i>Geometry</i>, size <i>double precision</i>, edges_only <i>boolean</i> ) : grid <i>Geometry</i><br>
ST_xxxxGrid( input <i>Geometry</i>, size <i>double precision</i>, edges_only <i>boolean</i>, origin <i>Geometry</i> ) : grid <i>Geometry</i>
</td></tr></table>
<ul>
<li>the <b>input</b> Geometry is always expected to be a Polygon or a MultiPolygon, and will be exactly covered by the returned grid.</li>
<li>the <b>size</b> argument identifies the edge length of the grid cell.</li>
<li>the optional <b>edges_only</b> argument will be interpreted as follows:<ul>
<li>if <b>FALSE</b> (<i>default value</i>) a MultiPolygon will be returned.</li>
<li>if <b>TRUE</b> a MultiLinestring will be returned (simply representing the cells edges).</li>
</ul></li>
<li>the optional <b>origin</b> Geometry is always assumed to be a Point, and will identify the grid's origin.
By default a <b>(0, 0)</b> origin will be assumed.</li>
</ul></td></tr>
<tr><td>
<h4>using Square cells</h4>
<table bgcolor="#f0f0f0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT ST_SquareGrid(geometry, 10000)<br>
FROM regions<br>WHERE cod_reg = 9;</b>
</td></tr>
</table><br>
This SQL query will return a regular grid (square cells) covering Tuscany (<b>cod_reg=9</b>).<br>
Each grid's cell will have an edge length of exactly <b>10 Km</b>
</td>
<td><img border="1" src="http://www.gaia-gis.it/gaia-sins/write-view-pics/square-grid.png" alt="square grid"></td></tr>
<tr><td>
<h4>using Triangular cells</h4>
<table bgcolor="#f0f0f0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT ST_TriangularGrid(geometry, 10000)<br>
FROM regions<br>WHERE cod_reg = 9;</b>
</td></tr>
</table></td>
<td><img border="1" src="http://www.gaia-gis.it/gaia-sins/write-view-pics/tri-grid.png" alt="triangular grid"></td></tr>
<tr><td>
<h4>using Hexagonal cells</h4>
<table bgcolor="#f0f0f0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT ST_HexagonalGrid(geometry, 10000)<br>
FROM regions<br>WHERE cod_reg = 9;</b>
</td></tr>
</table></td>
<td><img border="1" src="http://www.gaia-gis.it/gaia-sins/write-view-pics/hex-grid.png" alt="hexagonal grid"></td></tr>
</table>
<br><hr>
<h3>Delaunay Triangulations</h3>
Triangulations simply represent a special case of tessellations: in this case all cells are represented by generic <b>triangles</b>, not necessarily of the equilateral kind.<br>
A <b>Delaunay Triangulation</b> is very peculiar because a specific constraint is imposed.<br>
All triangles in a Delaunay triangulation must satisfy the <u><i>empty circle</i></u> property: i.e. for each edge we can find a circle containing the edge's endpoints but not containing any other points.
<a href="http://en.wikipedia.org/wiki/Delaunay_triangulation">read more</a><br><br>
Computing a Delaunay Triangulation representing many points is a very complex operation, and one possibly imposing a huge computational load and may be requiring a long time.
Happily enough many highly efficient algorithms have been already developed for the Delaunay problem.<br>
The next-to-come <a href="http://trac.osgeo.org/geos/">GEOS 3.4.0</a> will support Delaunay Triangulations; this GEOS version is currently still under active development, but the <i>experimental</i> base-code (<i>trunk</i>) seems to be stable enough to be safely tested.
SpatiaLite already relies on GEOS for many tasks, so integrating a smooth support for Delaunay Triangulation as well wasn't at all difficult.<br><br>
<table bgcolor="#f0f0f0" cellspacing="4" cellpadding="8">
<tr><td>
<i>prototype</i>:<hr>
ST_DelaunayTriangulation( input <i>Geometry</i> ) : delaunay <i>Geometry</i><br>
ST_DelaunayTriangulation( input <i>Geometry</i>, edges_only <i>boolean</i> ) : delaunay <i>Geometry</i><br>
ST_DelaunayTriangulation( input <i>Geometry</i>, edges_only <i>boolean</i>, tolerance <i>double precision</i> ) : delaunay <i>Geometry</i>
</td></tr></table>
<ul>
<li>the <b>input</b> Geometry can be of absolutely arbitrary type; all Linestrings and / or Polygowns will eventually be dissolved into Points corresponding to vertices.<br>
So after all ST_DelaunayTriangulation() will always process a MultiPoint.</li>
<li>the optional <b>edges_only</b> argument will be interpreted as follows:<ul>
<li>if <b>FALSE</b> (<i>default value</i>) a MultiPolygon will be returned.</li>
<li>if <b>TRUE</b> a MultiLinestring will be returned (simply representing the triangles edges).</li>
</ul></li>
<li>the optional <b>tolerance</b> argument is intended to normalize the input point-set, suppressing repeated (or too much close) points.
By default a <b>0.0</b> tolerance will be assumed.</li>
</ul></td></tr></table>
<table cellspacing="4" cellpadding="4">
<tr><td>
<table bgcolor="#f0f0f0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT ST_DelaunayTriangulation(ST_Collect(geometry))<br>
FROM italy_populated_places;</b>
</td></tr>
</table><br>
This SQL query will return a Delaunay Triangulation based on Italy's populated places (about 8,000+ Points).<br>
The visual example simply covers Tuscany, so the ensure an easy readability.
</td>
<td><img border="1" src="http://www.gaia-gis.it/gaia-sins/write-view-pics/delaunay.png" alt="delaunay triangulation"></td></tr>
</table>
<br><hr>
<h3>Voronoj Diagrams</h3>
<table cellspacing="4" cellpadding="4">
<tr><td align="center"><img border="1" src="http://www.gaia-gis.it/gaia-sins/write-view-pics/voronoj-1.png" alt="voronoj-delaunay relationship">
</td><td>
The <b>Voronoj Diagram</b> simply is the <i><u>dual graph</u></i> of the <b>Delaunay Triangulation</b>.
A Voronoj Diagram still is a tessellation, and cells in a Voronoj can have any arbitrary polygonal (irregular) shape.<br><br>
This figure clearly shows the relation joining the Delaunay Triangulation and the corresponding Voronoj Diagram.<br><br>
Each single Voronoj's cell is obtained by connecting all the <i>circumcenters</i> of adjacent Delaunay's triangles; consequently each Voronoj cell surely contains one (and only one) Delaunay's node.<br><br>
The cell in the Voronoj Diagram presents an interesting property: all points falling within the same cell are ensured to be nearest to the Delaunay's node placed on the cell itself than to any other Delaunay's node placed in a different cell.<br><br>
So the Voronoj Diagram is a very effective conceptual tool allowing to divide an arbitrary space region into many <i>rationally defined</i> cells, and is thus widely used on many applicative fields based on Computational Geometry.
This obviously including Geography itself.<br><br>
<a href="http://en.wikipedia.org/wiki/Voronoi_tessellation">read more</a>
</td></tr>
<tr><td align="center">
<img border="1" src="http://www.gaia-gis.it/gaia-sins/write-view-pics/voronoj-open.png" alt="voronoj open">
</td><td>
<u>Please note</u>: any Voronoj Diagram will always present an unpleasant complication; when reaching the exterior boundary of the input point-set, there aren't enough points allowing to properly close a cell.<br>
And consequently many cells laying near to the exterior boundary actually are <i>open cells</i>, presenting edges going directly to the <i>infinite</i>.<br><br>
The directions of such edges are absolutely well determined; but assigning them a finite length isn't allowed at all. And consequently isn't possible defining a polygon representing these cells.<br><br>
This figure shows a limited Voronoj Diagram always based on the same dataset as above, but this time including only Populated Places laying on the Island of Elba.
</td></tr>
<tr><td align="center">
<img border="1" src="http://www.gaia-gis.it/gaia-sins/write-view-pics/voronoj-frame.png" alt="voronoj frame">
</td><td>
Anyway we can usefully apply a <i>trick</i> effectively allowing to ensure that all cells in a Voronoj will be properly <i>closed</i>.<br><br>
We simply have to constrain the Voronoj Diagram within a rectangular frame of arbitrary size, then computing the intersections between infinite lines and the frame itself.<br><br>
This figure shows the same Voronoj as before, this time opportunely framed; as you can easily notice, now all cells are surely closed and can be safely be represented by a well defined polygon.<br><br>
The size (aka extent) of the frame is absolutely arbitrary, and can be freely set as you wish better accordingly to your specific requirements.<br>
<u>Please notice</u>: the blue and the red frame have a strongly different extent, but both them have the same identical overall shape.
</td></tr>
</table>
<a href="http://trac.osgeo.org/geos/">GEOS 3.4.0</a> will not directly support Voronoj; so the SpatiaLite own support is based on an original implementation (based in turn on the Delaunay support made available by GEOS).<br><br>
<table bgcolor="#f0f0f0" cellspacing="4" cellpadding="8">
<tr><td>
<i>prototype</i>:<hr>
ST_VoronojDiagram( input <i>Geometry</i> ) : voronoj <i>Geometry</i><br>
ST_VoronojDiagram( input <i>Geometry</i>, edges_only <i>boolean</i> ) : voronoj <i>Geometry</i><br>
ST_VoronojDiagram( input <i>Geometry</i>, edges_only <i>boolean</i>, extra_frame_size <i>double precision</i> ) : voronoj <i>Geometry</i><br>
ST_VoronojDiagram( input <i>Geometry</i>, edges_only <i>boolean</i>, extra_frame_size <i>double precision</i>, tolerance <i>double precision</i> ) : voronoj <i>Geometry</i>
</td></tr></table>
<ul>
<li>the <b>input</b> Geometry can be of absolutely arbitrary type; all Linestrings and / or Polygowns will eventually be dissolved into Points corresponding to vertices.<br>
So after all ST_VoronojDiagram() (exactly as ST_DelaunayTriangulation) will always process a MultiPoint.</li>
<li>the optional <b>edges_only</b> argument will be interpreted as follows:<ul>
<li>if <b>FALSE</b> (<i>default value</i>) a MultiPolygon will be returned.</li>
<li>if <b>TRUE</b> a MultiLinestring will be returned (simply representing the triangles edges).</li>
</ul></li>
<li>the optional <b>extra_frame_size</b> allows to freely set the extent of the frame.<br>
This argument is intended as a <i>percent</i> increase of the <i>natural extent</i> of Voronoj as determined by evaluating all Delaunay's nodes and circumcenters.<br>
By default a <b>5% extra_frame_size</b> is assumed.</li>
<li>the optional <b>tolerance</b> argument is intended to normalize the input point-set, suppressing repeated points (simply used when internally computing the Delaunay Triangulation).
By default a <b>0.0</b> tolerance will be assumed.</li>
</ul></td></tr></table>
<table cellspacing="4" cellpadding="4">
<tr><td>
<table bgcolor="#f0f0f0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT ST_VoronojDiagram(ST_Collect(geometry))<br>
FROM italy_populated_places;</b>
</td></tr>
</table><br>
This SQL query will return a Voronoj Diagram based on Italy's populated places (about 8,000+ Points).<br>
The visual example simply covers Tuscany, so the ensure an easy readability.<br>
All <i>populated places</i> (aka cell <i>seeds</i>) are explicitly represented.
</td>
<td><img border="1" src="http://www.gaia-gis.it/gaia-sins/write-view-pics/voronoj-2.png" alt="voronoj diagram"></td></tr>
</table>
<br><hr>
<h3>Delaunay Triangulation, Convex Hull and Concave Hull</h3>
<table cellspacing="4" cellpadding="4">
<tr><td><img border="1" src="http://www.gaia-gis.it/gaia-sins/write-view-pics/convex.png" alt="delaunay-convexhull relationship">
</td><td>
An interesting point absolutely worth to be explicitly noted: the boundary of any <b>Delaunay Triangulation</b> exactly corresponds to the <b>Convex Hull</b> for the same input Geometry.
<a href="http://en.wikipedia.org/wiki/Convex_hull">read more</a><br><br>
And this in turn opens the way to a further consideration: we could purposely simplify a Delaunay Triangulation so to get a <b>concave hull</b>.<br>
You can get more extensive informations about this approach from <a href="http://grass.osgeo.org/wiki/Create_concave_hull">here</a><br><br>
<u>Please note well</u>: the <b>convex hull</b> concept corresponds to a robust and formal mathematical definition. On the other side the <b>concave hull</b> is a much more vague and undetermined notion.<br>
There is one and only one ConvexHull for a given Geometry; but many ConcaveHulls are possible.
Choosing the one or the other is much more a matter of personal taste than a mathematical operation formally defined.<br>
Computing a ConcaveHull always is an inherently arbitrary and heuristic process.
</td></tr>
<tr><td>
<h4>The SpatiaLite's own approach to ConcaveHull</h4>
<ol>
<li>the Delaunay Triangulation corresponding to the input Geometry will be computed first.</li>
<li>then the statistical distribution of all triangle edge's lengths will be evaluated, so to determine <b>σ</b> (aka the <i>standard deviation</i>)</li>
<li>a second pass will now examine yet again all Delaunay's triangles; any triangle presenting at least one edge longer than <b>σ * factor</b> will be discarded.</li>
<li>and finally all filtered triangles will be dissolved so to form the ConcaveHull to be returned.</li>
</ol>
<u>Please note</u>: by setting an appropriate value to <b>factor</b> you can freely influence how much <i>aggressive</i> the filtering pass will be.<br>
Adopting a very high <b>factor</b> value practically means applying a very bland filtering (very few triangles will be discarded, and you'll consequently get a rather <i>convex</i> shape).
On the other side adopting a very low <b>factor</b> values will apply a very strong filtering (many triangles will be now discarded, and you'll get a very <i>concave</i> shape).<br><br>
<u>Useful constants</u>: assuming a perfectly normal distribution of edge lengths (a by far unrealistic assumption for real world datasets), <b>3σ</b> will imply suppressing about <b>0.1%</b> triangles (only the few ones presenting abnormally lengthy edges), <b>2σ</b> corresponds to about <b>2.1%</b>, and <b>1σ</b> roughly corresponds to <b>15.8%</b>.<br>
Usually values ranging between <b>3σ</b> and <b>2σ</b> are the most appropriate to be used.<br><br>
The figure shows what you can actually get by applying a <b>3σ</b> filtering to Italian Populated Places; for clarity all Italian Regions are represented in yellow, and the ConvexHull is represented in blue. The ConcaveHull itself is represented in red.
</td><td>
<img border="1" src="http://www.gaia-gis.it/gaia-sins/write-view-pics/concave3.png" alt="concavehull-1">
</td></tr>
</table>
<br>
<b><u>Please note well</u></b>: PostGIS supports its own implementation of ST_ConcaveHull(), being based on a completely different approach. You can get more details from <a href="http://www.bostongis.com/postgis_concavehull.snippet">here</a><br>
Please, don't be confused: the SpatiaLite and PostGIS implementations are completely unrelated. Anyway, it could be interesting noting that the SpatiaLite's own implementation seems to be much more efficient and fast, and can be realistically applied even to very huge point-sets.<br><br>
<table bgcolor="#f0f0f0" cellspacing="4" cellpadding="8">
<tr><td>
<i>prototype</i>:<hr>
ST_ConcaveHull( input <i>Geometry</i> ) : concave <i>Geometry</i><br>
ST_ConcaveHull( input <i>Geometry</i>, factor <i>double precision</i> ) : concave <i>Geometry</i><br>
ST_ConcaveHull( input <i>Geometry</i>, factor <i>double precision</i>, allow_holes <i>boolean</i> ) : concave <i>Geometry</i><br>
ST_ConcaveHull( input <i>Geometry</i>, factor <i>double precision</i>, allow_holes <i>boolean</i>, tolerance <i>double precision</i> ) : concave <i>Geometry</i>
</td></tr></table>
<ul>
<li>the <b>input</b> Geometry can be of absolutely arbitrary type; all Linestrings and / or Polygowns will eventually be dissolved into Points corresponding to vertices.<br>
So after all ST_ConcaveHull() (exactly as ST_DelaunayTriangulation) will always process a MultiPoint.</li>
<li>the optional <b>factor</b> is intended to be a <b>σ</b> multiplier, thus determining the filtering threshold.<br>
By default a <b>3σ factor</b> is assumed.</li>
<li>the optional <b>allow_holes</b> argument will be interpreted as follows:<ul>
<li>if <b>FALSE</b> (<i>default value</i>) any eventual interior hole will be suppressed.</li>
<li>if <b>TRUE</b> any eventual interior holes will be preserved.</li>
</ul></li>
<li>the optional <b>tolerance</b> argument is intended to normalize the input point-set, suppressing repeated points (simply used when internally computing the Delaunay Triangulation).
By default a <b>0.0</b> tolerance will be assumed.</li>
</ul></td></tr></table>
<table cellspacing="4" cellpadding="4">
<tr><td>
<table bgcolor="#f0f0f0" cellspacing="4" cellpadding="8">
<tr><td>
<b>SELECT ST_ConcaveHull(ST_Collect(geometry), 2.5)<br>
FROM italy_populated_places;</b>
</td></tr>
</table><br>
This SQL query will return a ConcaveHull based on Italy's populated places (about 8,000+ Points) by applying a <b>2.5σ</b> filtering factor.<br><br>
The figure also shows the ConcaveHull returned by the previous example by applying a <b>3σ</b> factor (red boundary), so to evidentiate the relative differences.
</td>
<td><img border="1" src="http://www.gaia-gis.it/gaia-sins/write-view-pics/concave25.png" alt="concavehull-2"></td></tr>
</table>
<br><hr>
<a href="https://www.gaia-gis.it/fossil/libspatialite/wiki?name=4.2.0-doc">back</a>
Z 082c79084944472884d17a7a1e79458d