MySQL è uno dei più diffusi database server, e non deve stupire che abbia un’estensione spaziale (peraltro attiva di default). Mi stupisce però che le funzioni per testare le relazioni spaziali tra oggetti siano limitate.

Mi spiego con un esempio classico: selezionare tutti i punti di un layer che cadono all’interno del perimetro di un poligono presente in un altro layer. E’ un’operazione tipica dei GIS, e presente in tutti database con estensione spaziale.

In MySQL basta lanciare una query spaziale di questo tipo:

SELECT * FROM poligoni as g1, punti as g2 WHERE Contains(g1.geometry,g2.geometry) = 1

Purtroppo però l’output non è costituito da tutti i punti contenuti nel perimetro di mio interesse, ma da tutti quelli contenuti nel rettangolo che lo include. Quest’ultimo è il classico Bounding Box, o come si definisce in ambiente MySQL Minimal Bounding Rectangles (MBR). Nel manuale online è riportato:

The OpenGIS specification defines the following functions. They test the relationship between two geometry values g1 and g2.

Questo dovrebbe garantire che l’output della funzione Contains() sia relativo al perimetro dell’oggetto “target”, ma poco più avanti si legge:

Currently, MySQL does not implement these functions according to the specification. Those that are implemented return the same result as the corresponding MBR-based functions.

Immagino che i GeoSpatial Developers abbiano diverse frecce al loro arco per superare questo problema, così come ne esisteranno diverse lato utente. Una potrebbe essere quella di non usare questo database server e fare tutto con PostGIS, Oracle Spatial, Spatialite, ecc., ma è troppo facile e non sempre si può scegliere.

Io ho pensato a GDAL/OGR ed alle sue utility per oggetti vettoriali, ed alla possibilità (che esiste dalla versione 1.10 della libreria) di usare il dialetto SQLite/Spatialite. Questo dialetto estende di molto quello che queste utility fanno (egregiamente) di default, ovvero eseguire delle query sql all’interno di un comando; ad esempio:

ogrinfo province.shp province -sql "SELECT nome FROM province WHERE ID_PRO = 2"

Con il dialetto SQLite/Spatialite ho a disposizione anche le funzioni per verificare relazioni spaziali tra oggetti, anche per quelle basi dati che non prevedono intrinsecamente la possibilità di farlo, proprio come MySQL spatial.

Andando nel concreto dovrei scrivere una cosa di questo tipo:

ogrinfo MYSQL:"mydb,user=myuse,password=mypwd,port=3306" -dialect sqlite -sql "SELECT *FROM poligoni as g1, punti as g2 WHERE Contains(g1.geometry,g2.geometry) = 1"

La sorpresa è che, nonostante si dichiari nella stringa il “dialetto” sqlite, si ottiene sempre come output quello relativo ai punti contenuti nel bounding box del perimetro di interesse. E’ un baco? Nella documentazione di OGR si legge:

All OGR drivers for database systems: MySQL, PostgreSQL and PostGIS (PG), Oracle (OCI), SQLite, ODBC, ESRI Personal Geodatabase (PGeo) and MS SQL Spatial (MSSQLSpatial), override the OGRDataSource::ExecuteSQL() function with dedicated implementation and, by default, pass the SQL statements directly to the underlying RDBMS. In these cases the SQL syntax varies in some particulars from OGR SQL. Also, anything possible in SQL can then be accomplished for these particular databases. Only the result of SQL WHERE statements will be returned as layers.

Ma si legge anche:

The SQLite dialect may be used with any OGR datasource, like the OGR SQL dialect. It is available through the OGRDataSource::ExecuteSQL() method by specifying the pszDialect to “SQLITE”. For the ogrinfo or ogr2ogr utility, you must specify the “-dialect SQLITE” option.

La soluzione l’ho trovata in un test stupido che ho voluto fare: accedere alla fonte MySQL non direttamente, ma tramite il Virtual Format di OGR, che in qualche modo astrae l’accesso al formato di input.

A partire quindi da una fonte MySQL costruita secondo le specifiche del Virtual Format e salvata come “test.vrt”:

<ogrvrtdatasource>
<ogrvrtlayer name="poligoni">
<srcdatasource>MYSQL:mydb,user=myuser,
password=mypwd,port=3306,host=127.0.0.1
</srcdatasource>
<srclayer>poligoni</srclayer>
<geometrytype>wkbPolygon</geometrytype>
<layersrs>epsg:4326</layersrs>
</ogrvrtlayer>
<ogrvrtlayer name="punti">
<srcdatasource>MYSQL:mydb,user=myuser,
password=mypwd,port=3306,host=127.0.0.1
</srcdatasource>
<srclayer>punti</srclayer>
<geometrytype>wkbPoint</geometrytype>
<layersrs>epsg:4326</layersrs>
</ogrvrtlayer>
</ogrvrtdatasource>

Posso ad esempio lanciare:

ogrinfo test.vrt -dialect sqlite -sql "SELECT *FROM poligoni as g1, punti as g2 WHERE Contains(g1.geometry,g2.geometry) = 1"

Il risultato ottenuto sarà così quello desiderato e l’output corrisponderà a tutti i punti che ricadono nel perimetro di interesse. Nell’immagine di sotto, ho raccontato il tutto in modo visuale.

Concludo sottolineando  quanto siano interessanti, utili, belle e di alto livello professionale alcune delle dinamiche delle comunità open source. Ho segnalato quello che a me sembrava un baco alla lista GDAL – DEV. Non solo ho ottenuto subito delle risposte che mi hanno consentito di capirne molto di più, ma si sono messi già al lavoro sia in termini di codice che di documentazione. Un semplice “grazie” sembra veramente poco.

Questo post è dedicato a Flaviano, un mio amico che grazie anche ad alcuni tips & tricks su GDAL/OGR ha trovato un bel posto di lavoro in Qatar

Pensato e scritto da Andrea Borruso

Sono un geologo marino e nel 1999 ho iniziato con Antonio D'Argenio la vita da libero professionista nel mondo dei Sistemi Informativi Geografici. Nei primi anni - di tanto studio e di poche fatture - è stato fatale l’incontro con MapServer: il web-mapping e le applicazioni web sono diventati due pezzi importanti della mia professione. Negli ultimi anni mi sono occupato anche di formazione e di recente sono entrato nello strano mondo del telerilevamento di prossimità. Mi piace molto leggere di GIS, informazione geografica, e neogeography, e parlarne con gli amici che incontro spesso al bar dietro al router.