PostGIS: query su dati georiferiti

PostGIS è l'estensione spaziale del server PostgreSQL che introduce i tipi di dato geometrico e le funzioni per lavorare con essi. Inoltre fornisce le definizioni dei sistemi di coordinate (derivati dalla classificazione EPSG) per eseguire trasformazioni tra geometrie materializzate in sistemi di riferimento diversi.

Quali tipi di oggetti geometrici possono essere utilizzati?

Attenendosi ai formati descritti dall'Open GIS Consortium, per garantire trasparenza ed interoperabilità, vengono usati i tipi: POINT, MULTIPOINT, LINESTRING, MULTILINESTRING, POLYGON, MULTIPOLYGON, GEOMETRYCOLLECTIONS (più le estensioni SQL-MM CIRCULARSTRING, COMPOUNDCURVE, CURVEPOLYGON, MULTICURVE, MULTISURFACE) con estensione XYZ,XYM, XYZM.

Come funziona un database spaziale?

Per ogni strato informativo occorre definire una colonna che contenga le coordinate, ed un sistema che gestisca in maniera differenziata le tabelle che contengono dati geometrici di tipo differente. Occore inoltre introdurre un insieme di vincoli per verificare che i dati inseriti siano congruenti: ad esempio dovrebbe essere sempre verificata la bi- o tridimensionalità nello stesso dataset. Interrogazioni veloci alle tabelle sono la ragione d'essere dei database (assieme al supporto per le transazioni). I database si differenziano proprio per la maggior robustezza e prestanza degli indici. Le elaborazioni GIS avvengono usando la sintassi SQL sui costrutti "spaziali". Esistono comunque funzioni che costruiscono immediatamente tutta la struttura necessaria alla gestione dei dati territoriali. Analogamente le interrogazioni al database avvengono utilizzando ``SQL query'' utilizzando combinazioni di dati e funzioni e di test vero/falso. Per comprendere come funzioni una query spaziale occorre tenere presenti due cose: Gli indici spaziali esistono per migliorare l'efficienza delle query. I dati geografici vengono "aggregati" e "amministrati" in gruppi spaziali distinti. In questo senso gli indici sono "lossy": sono definiti per semplificare e nella semplificazione perdono informazioni. Ad esempio: spesso si vuole utilizzare l'operatore intersezione && che testa se il rettangolo che circoscrive le entità geometriche ne intersechi altre. Questa funzione e' ottimizzata per l'utilizzo degli indici spaziali: serve a scremare i dati per esegure una ricerca più raffinata usando le funzioni Distance, Intersects, Contains, Within... Le funzioni che utilizzano l'indice spaziale servono per identificare le geometrie che possono soddisfare ad una condizione. Le funzioni spaziali testano la condizione esattamente.
PostgreSQL utilizza 3 tipi di indici: B-Tree, R-Tree e GiST (Generalized Search Tree). GRASS, per i dati vettoriali, usa R-Tree, PostGIS usa GiST. Per PostgreSQL, R-Tree ha l'inconveniente di NON essere "null safe": non possono essere immagazzinate geometrie NULLE, degeneri o tipi diversi all'interno della stessa colonna. Sinteticamente, PostGIS manipola dati "congruenti".

Prerequisiti

Si descrive la procedura per importare dati geometrici in formato shapefile nel database. Esiste un programma di importazione/esportazione chiamato loader/dumper: dapprima si esporti il vettoriale di GRASS "geology" in formato shapefile, tenendo presente, nell'esportazione in formato shapefile che si stanno esportando delle aree (opzione type=line,boundary,area), di esportare le linee come poligoni (flag -p) ed, ovviamente, di esportare solamente le geometrie alle quali sono associate delle informazioni (-c); ATTENZIONE! L'operazione può essere piuttosto lunga se i file sono di grandi dimensioni. Il formato shapefile risulta essere "de facto" lo standard di interscambio di dati georeferenziati: PostGIS fornisce quindi un programma dedicato per trasformare le informazioni contenute da formato ESRI binario ad un file di testo contenente gli SQL-statement che popoleranno la tabella nel database. Esistono altri metodi per importare (o scambiarsi) informazioni tra GIS: è possibile trasferire files direttamente da GRASS, oppure utilizzare la libreria OGR. I comandi da eseguire in una shell, una volta posizionatisi nella directory dove risiede lo shapefile sono:

shp2pgsql -s 32119 -I geology geology2postgis > geology.sql

La flag "-I" crea automaticamente un indice di tipo GiST (Generalized Search Tree) sulla colonna contenente le geometrie: ciò aumenta notevolmente la performance delle query su dati geografici in quanto, attraverso l'uso delle bounding box permette di effettuare le ricerce su subsets di dati. v.in.ogr non crea in automatico l'indice sulla colonna delle geometrie. La flag "-s" con il relativo numero associa ai dati geografici ad un preciso datum (quello del Carolina), streets è il nome dello shapefile, streets2postgis è il nome della tabella che verrà creata, streets.sql è il nome del file di testo contenente gli SQL-statements. Si procede quindi al caricamento dei dati nel database:

psql -d gisdemo_ncspm -U grass -f geology.sql

Adesso se andiamo ad osservare la tabella geometry_columns con PgAdminIII secondo le procedure descritte (collegandosi al database gisdemo_ncspm), si nota come le colonne abbiano nomi diversi a seconda del meccanismo di esportazione dei dati: "the_geom" nel caso del loader, "wkb_geometry" nel caso dell'OGR.
Si dovrebbe ripetere l'operazione con il file "streets" ma per comodità, visto che l'esportazione in shape è lenta a causa delle dimensioni del file lo shape file è scaricabile direttamente qui: streets.zip una volta scaricato e scompattato streets eseguire in una shell, una volta posizionatisi nella directory dove risiede lo shapefile:
shp2pgsql -s 32119 -I streets streets2postgis > streets.sql
e
psql -d gisdemo_ncspm -U grass -f streets.sql

Geometrie e datum

Qual'è la definizione del datum associato ai dati appena importati? Interrogando la tabella "spatial_ref_sys" è possibile risalire a questa informazione. Si deve usare un ambiente in cui sia possibile effettuare delle query: si può usare PgAdminIII oppure l'ambiente a riga di comando psql con:
psql -d gisdemo_ncspm -U grass
NOTA: per uscire dall'ambiente psql, digitare \q Si esegue quindi la query con:
SELECT srid,proj4text FROM spatial_ref_sys WHERE srid=32119;
che restituisce
 srid  |                                                                       proj4text
-------+--------------------------------------------------------------------------------------------------------------------------------------------------------
 32119 | +proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs

HWKB e coordinate

Per garantire l'efficienza e lo standard di immagazzinamento dei dati, OGC propone di trasformare le coordinate in Well-Known-Binary formato esadecimale. Infatti, volendo consocere le caratteristiche del primo tratto di autostrada si digita il comando SQL:

SELECT stype AS descrizione, the_geom AS coordinate FROM streets2postgis WHERE stype='HWY' LIMIT 1;

ottenendo:
 descrizione |                                                                  coordinate
-------------+----------------------------------------------------------------------------------------------------------------------------------------------
 HWY         | 0105000020777D0000010000000102000000030000001D384F205F312441379AB876B1DD0B411C2EB428C231244175D592737FDE0B4157BA70D7433224415587DB7D90DF0B41
Il formato risulta inintelleggibile, aiuta, quindi la funzione di trasformazione delle geometrie AsText():

SELECT stype AS descrizione, AsText(the_geom) AS coordinate FROM streets2postgis WHERE stype='HWY' LIMIT 1;

che restituisce:
 descrizione |                                                        coordinate
-------------+--------------------------------------------------------------------------------------------------------------------------
 HWY         | MULTILINESTRING((661679.563104394 228278.182969289,661729.079499665 228303.931432407,661793.920781921 228338.061453874))

Aggregazione dei dati

Si vuole rispondere a questa domanda: qual'è la lunghezza totale, espressa in chilometri, dei boulevard? Si utilizza allora la query:

SELECT sum(st_length(the_geom))/1000 AS km_roads FROM streets2postgis WHERE stype='BLVD';

In questo caso si crea un dato aggregato partendo filtrando tutti e 799 boulevards immagazzinati nel database nel database. Per ogni tratto estratto viene calcolata la lunghezza che viene poi sommata.
        km_roads
------------------
 183.098038103509 

Aggregazione 1:molti

Si vuole produrre l'elenco delle 5 aree geologiche più grandi con substrato di roccia granitica (classificazione USGS) ordinate per grandezza (in ettari):

SELECT cat AS n_area, ST_Area(the_geom)/10000 AS superficie from geology2postgis WHERE geo_name='DOg' ORDER BY superficie desc LIMIT 5;

Questa query combina una condizione sull'attributo ed una di calcolo spaziale:
 n_area |    superficie
--------+------------------
   1321 | 18276.4251066893
   1081 | 4632.82475102043
   1193 | 4513.41760449962
   1112 | 1477.15172200672
   1394 | 1288.61212637985

Aggregazione 1:1

Dai dati precedenti siamo incuriositi quanto sia esteso, in totale, il substrato geologico di roccia granitica (in ettari):

SELECT geo_name AS tipologia, sum(ST_Area(the_geom))/10000 AS superficie FROM geology2postgis GROUP BY geo_name HAVING geo_name='DOg';

Per rispondere a questa domanda si e' calcolata, con la funzione ST_Area, l'estensione di ogni singolo poligono.
 tipologia |    superficie
-----------+------------------
 DOg       | 32404.5197312962

Spatial join

Qual'è la lunghezza dei boulevard contenuti all'interno di ogni unità geologica?

SELECT g.geo_name AS tipologia, sum(st_length(s.the_geom))/1000 AS km_boulevard FROM geology2postgis AS g, streets2postgis AS s WHERE s.stype='BLVD' AND s.the_geom && g.the_geom AND st_contains(g.the_geom, s.The_geom) GROUP BY g.geo_name ORDER BY km_boulevard;

Questo è un esempio di "spatial join" poichè stiamo unendo i dati da due tavole (join) usando una operazione spaziale (st_contains) come condizione di join piuttosto che il tradizionale approccio di compiere un join tra tavole usando una chiave comune:
 tipologia |   km_boulevard                                                                                                                                               
-----------+------------------
 Km        | 1.13581818201085
 Tt        | 1.74814631019625
 CZbg      | 3.96241813789034
 CZve      |  6.0151866702079
 CZfg      |  7.1171721781027
 CZlg      | 12.1324471963735
 TRc       |  25.751978510117
 PPmg      | 48.1958958470183
 CZig      | 67.6738970427193

Miglioramento della viabilità

Il piano regionale di miglioramento della viabilità prevede un allargamento dei boulevard. L'entità dello sbancamento per l'allargamento delle strade (compresa l'area di rispetto) è di 6m. Nel progetto preliminare si evidenzia la presenza di solida roccia granitica (PPmg) difficilmente sbancabile. Considerando che i lavori di sbancamento con quelle condizioni aumentano di 150$ ogni mq di strada si vuole stimare a quanto lieviteranno i costi:
SELECT 
  g.geo_name AS tipologia, 
  round(sum(st_length(st_intersection(s.the_geom,g.the_geom)))) AS lunghezza_strada, 
  round(sum(st_area(st_buffer(st_intersection(s.the_geom,g.the_geom),3)))) AS area_allargamento, 
  round(sum(st_area(st_buffer(st_intersection(s.the_geom,g.the_geom),3)))*150) AS costo_aggiuntivo 
FROM 
  streets2postgis AS s, 
  geology2postgis AS g 
WHERE 
  s.the_geom && g.the_geom AND 
  intersects(s.the_geom, g.the_geom) AND 
  g.geo_name = 'PPmg' AND 
  s.stype = 'BLVD' 
GROUP BY 
  g.geo_name;
Il risultato della query è:
 tipologia | lunghezza_strada | area_allargamento | costo_aggiuntivo
-----------+------------------+-------------------+------------------
 PPmg      |            50223 |            307233 |         46084960
46 milioni di dollari! Come esercizio: si sarebbe potuto scrivere la query in modo più efficiente riducendo il numero di operazioni ed ottenendo un risultato più velocemente?


La stessa query Miglioramento della viabilità eseguita via interfaccia grafica di PgadminIII