PostGIS: query su dati georiferiti

PostGIS è l'estensione spaziale al server PostgreSQL. Essa introduce tipi di dati geometrici e funzioni per lavorare con essi. Inoltre fornisce un set di sistemi di coordinate EPSG per lavorare con questi dati geometrici ed eseguire trasformazioni tra sistemi di riferimento.

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, line, polygon, multipolygon e geometrycollections 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 sistema che verifichi che i dati inseriti siano congruenti e che venga sempre verificata la bi- o tridimensionalità. 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. Tutto avviene usando la sintassi SQL. 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, funzioni e di test vero/falso. Per comprendere come funzioni una query spaziale occorre tenere presenti due cose:

esistono gli indici spaziali;
le interrogazioni spaziali sono molto dispendiose in termini di calcolo se effettuate su un gran numero di entità geometriche.

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 interseca 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". PostgreSQL non vuole, insomma, geometrie NULLE o tipi diversi all'interno della stessa colonna. PostgreSQL è un database che, per come è stato costruito, manipola dati "congruenti".

Prerequisiti

Per svolgere questa esercitazione è necessario avere svolto l'esercitazione precedente (PostgreSQL e PostGIS) in cui si crea la tabella fields2postgis. Per provare ad usare il programma di importazione dati di PostgresSQL (loader/dumper) dapprima si esporti il vettoriale di GRASS "roads" in formato shapefile, creando uno shapefile di nome roads con

v.out.ogr input=roads type=line,boundary dsn=/home/ubuntu olayer=roads layer=1 format=ESRI_Shapefile

PostGIS fornisce un programma per trasformare le informazioni contenute da formato ESRI Shape ad un file di testo contenente gli SQL-statement che popoleranno la tabella nel database. I comandi da eseguire in una shell, una volta posizionatisi nella directory dove risiede lo shapefile sono:

shp2pgsql -I -s 32767 roads roads2postgis > roads.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 boundingbox, permette di effettuare le ricerce su subsets di dati (per contro v.out.ogr non crea in automatico l'indice sulla colonna delle geometrie). La flag "-s" con il relativo numero associa ai dati geografici un preciso datum (in questo caso quello della location Spearfish, NAD27 - Utm zona 13), roads è il nome dello shapefile, roads2postgis è il nome della tabella che verrà creata, roads.sql è il nome del file di testo contenente gli SQL-statements. Si procede quindi al caricamento dei dati nel database in un terminale (in GRASS o meno):

psql -d spearpg -U grass -f roads.sql

Osservando con pgAdmin la tabella geometry_columns si nota come le colonne abbiano nomi diversi a seconda del meccanismo di esportazione dei dati: "the_geom" nel caso di shape2pgsql (loader), "wkb_geometry" nel caso dell'OGR (cioè per fields2postgis).

Spearfish e datum

Qual è la definizione del datum associato ai dati di Spearfish? Interrogando la tabella "spatial_ref_sys" è possibile risalire a questa informazione, usando ad esempio il Query Builder di pgAdmin:

SELECT srid,proj4text FROM spatial_ref_sys WHERE srid=32767;

che restituisce

srid   |                            proj4text
-------+------------------------------------------------------------------
 32767 | +proj=utm +zone=13 +ellps=clrk66 +datum=NAD27 +units=m +no_defs
(1 row)

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 conoscere le caratteristiche del primo tratto di autostrada digitalizzata in Spearfish si digita il comando SQL:

SELECT label AS descrizione, the_geom AS coordinate FROM roads2postgis WHERE label ~ 'primary highway' LIMIT 1;

ottenendo:

          descrizione          |                                         coordinate
-------------------------------+--------------------------------------------------------------
 primary highway, hard surface | 0105000020FF7F000001000000010200000003000000D7E164E20C2322...

Il formato risulta inintelleggibile, aiuta quindi la funzione di trasformazione delle geometrie AsText():

SELECT label AS descrizione, AsText(the_geom) AS coordinate FROM roads2postgis WHERE label ~ 'primary highway' LIMIT 1;

che restituisce:

          descrizione          |                                         coordinate
-------------------------------+--------------------------------------------------------------
 primary highway, hard surface | MULTILINESTRING((594310.44217592 4925247.08617608,594308.8...

Aggregazione dei dati

Si vuole rispondere a questa domanda: qual'è la lunghezza totale, espressa in chilometri, dei sentieri nella contea di Spearfish? Si utilizza allora la query:

SELECT sum(st_length(the_geom))/1000 AS km_roads FROM roads2postgis WHERE label='unimproved road';

In questo caso si crea un dato aggregato selezionando tutti e 153 i sentieri immagazzinati (con label='unimproved road') nel database. Per ogni tratto estratto viene calcolata la lunghezza che viene poi sommata.

     km_roads
------------------
 136.313999118462

Estratto tavolare

Si vuole produrre l'elenco delle 10 più grandi particelle fondiarie di Spearfish ordinate per grandezza:

SELECT ogc_fid as particella, owner as proprietario, st_area(wkb_geometry)/10000 AS superficie FROM fields2postgis ORDER BY superficie DESC LIMIT 10;

Questa query combina una condizione sull'attributo ed una di calcolo spaziale:

 particella |       proprietario       |    superficie
------------+--------------------------+------------------
         24 | Black Hills Natl. Forest | 11651.4876800005
         38 | D. Portillo              | 503.438798999917
          1 | P. Biggam                | 305.164704499996
         25 | R. Reed                  | 292.013480499935
         62 | B. Cartwright            | 282.227241999912
         21 | G. Tandy                 | 278.218926000047
          9 | D. Liston                | 262.542519000065
          5 | P. Biggam                | 262.541994500053
         48 | C. Krenshaw              | 230.196791000044
         40 | D. Seagall               | 196.777930999994

Pur essendo Mr.Portillo il detentore del record "il fondo più grande", si vede che Mr.Biggam non se la cava poi così male, detenendo una superficie totale maggiore, almeno tra le prime 10.

Estratto di proprietà

Dai dati precedenti siamo incuriositi sull'effettiva disponibilità di terreno agricolo (in ettari) di Mr. Biggam.

SELECT owner AS proprietario, sum(st_area(wkb_geometry))/10000 AS superficie FROM fields2postgis GROUP BY owner HAVING owner~'Biggam';

Per rispondere a questa domanda si e' dovuto calcolare l'area di ogni singolo poligono. Se dovessimo compiere spesso questo calcolo avrebbe senso aggiungere una colonna AREA alla tabella ed un ulteriore indice per migliorare le prestazioni della query.

 proprietario |    superficie
--------------+------------------
 P. Biggam    | 891.545676000082

Spatial join

Qual e' la lunghezza dei sentieri contenuti all'interno di ogni proprietà?

SELECT f.owner AS proprietario, sum(st_length(r.the_geom))/1000 AS km_sentieri FROM fields2postgis AS f, roads2postgis AS r WHERE r.the_geom && f.wkb_geometry AND st_contains (f.wkb_geometry, r.the_geom) GROUP BY f.owner ORDER BY km_sentieri;

Questo e' un esempio di "spatial join" poiche' stiamo unendo i dati da due tavole (join) usando una condizione di interazione spaziale (st_contains) come condizione di join piuttosto che il tradizionale approccio di compiere un join tra tavole usando una chiave comune:

       proprietario       |    km_sentieri
--------------------------+--------------------
 P. Biggam                | 0.0730895522860136
 D. Seagall               |  0.463857456898311
 C. Mitchell              |   1.05953562101462
 Black Hills Natl. Forest |    75.190978824988
(4 rows)

Esproprio

Il piano regolatore di Spearfish prevede che i sentieri agricoli vengano trasformati in strade. La larghezza delle strade (compresa dell'area di rispetto) è di 6m. Si agirà per esproprio. Considerando che il valore del terreno agricolo a Spearfish è di 35$ al mq si vuole stimare a quanto ammonterà l'indennizzo per Mr. Mitchell.

SELECT
	f.owner AS proprietario, 
	round(sum(st_length(st_intersection(r.the_geom,f.wkb_geometry)))) AS lunghezza_strada,
	round(sum(st_area(st_buffer(st_intersection(r.the_geom,f.wkb_geometry),3)))) AS area_esproprio, 
	round(sum(st_area(st_buffer(st_intersection(r.the_geom,f.wkb_geometry),3)))*35) AS valore_terreno
FROM roads2postgis AS r, fields2postgis AS f
WHERE r.the_geom && f.wkb_geometry 
AND intersects(r.the_geom, f.wkb_geometry)
AND f.owner = 'C. Mitchell'
AND r.label = 'unimproved road'
GROUP BY f.owner
HAVING f.owner = 'C. Mitchell';

Il risultato della query è:

 proprietario | lunghezza_strada | area_esproprio | valore_terreno
--------------+------------------+----------------+----------------
 C. Mitchell  |             1420 |           8659 |         303081
(1 row)

La domanda ora è: mr.Mitchell si arrabbierà?


Previous: POSTGRESQL E GRASS	Up: ESERCITAZIONI	Next: NVIZ