INTRODUZIONE AI DATABASE GEOSPAZIALI

https://enricofer.github.io/spatial_sql_workshop/

FOSS4G Italia 2021

Esercitazione guidata

SINTASSI SQL

In linguaggio SQL sono possibili parecchi tipi di Query:

query di creazione tabella:

CREATE TABLE tab(id INT PRIMARY KEY NOT NULL, campo1 CHAR(50), campo2 INT);

query di inserimento dati:

INSERT INTO tab (id, campo1, campo2) VALUES (1, 'valore 1', 120);

query di cancellazione tabella:

DROP TABLE tab;

query di aggiornamento:

UPDATE tab SET campo1 = 'valore2' WHERE id = 1;

SQL SELECT

ma quello che ci interessa oggi è la query di ricerca

        SELECT id, campo1   -- prendi id e campo1
        FROM tabella        -- da tabella
        WHERE campo2 > 100; -- dalle righe che soddisfano una condizione

        SELECT * 
        FROM tabella -- elenca tutti i campi e tutte le righe di una tabella

I dialetti SQL sono molti e ci sono grandi differenze di dettaglio. Ma la struttura della sintassi principale per l'estrazione dei dati rimane praticamente la stessa per tutti i RDBMS.
Oggi lavoreremo con SQlite, nella sua versione con estensioni spaziali Spatialite il cui sviluppo è open source ed è curato da Alessandro Furieri

L'ambiente di lavoro del Workshop

Per facilitare e velocizzare il workshop le esercitazioni sono completamente basate su un una versione speciale di Spatialite compilata per l'utilizzo da browser. Il database è residente nella memoria del browser, ma è facilmente scaricabile per l'utilizzo con QGIS per chi lo desidera:dati esercitazione

E' possibile una visualizzazione tabellare:

      SELECT * FROM comuni LIMIT 10

o una visualizzazione su mappa: Le tabelle contengono un campo *geom* con le informazioni spaziali memorizzate in latitudine e longitudine (WGS84 epsg:4326) e visualizzabili producendo una geometria geojson per mezzo della funzione AsGeoJSON()

      SELECT *, AsGeoJSON(geom) AS geojson
      FROM comuni

Il formato geojson è un formato testuale facilmente interpretabile da javascript come un'oggetto (dati serializzati) contentente dati geografici e relativi attributi. Il formato è uno standard codificato di internet e costituisce una specie di lingua universale per l'interscambio di dati geografici nel web

I dati dell’esercitazione

Dati ISTAT (http://www.istat.it/it/archivio/104317)

Limiti provinciali del Piemonte (province poligonale)
Località abitate (localita poligonale)
Tabella TIPO_LOC Tipologia di località 2001/2011. (tipo_localita.csv)
Il campo può assumere i seguenti valori:
1. centro abitato
2. nucleo abitato
3. località produttiva
4. case sparse

Da Geoportale Regione Piemonte (http://www.geoportale.piemonte.it/cms/)
Limiti amministrativi comunali (comuni poligonale)
Una selezione della viabilità principale (autostrade lineare)

SELECT ... FROM ... WHERE ...

query di ricerca SELECT

Disponiamo delle seguenti tabelle, visualizzabili eseguendo le query:

SELECT * FROM comuni     -- Confini comunali del Piemonte (geometrie poligonali)

SELECT * FROM viab_princ     -- Viabilità principale (geometrie lineari)

SELECT * FROM localita   -- Località ISTAT (geometrie poligonali)

SELECT * FROM province      -- Province (geometrie poligonali)

SELECT * FROM tipo_localita     -- Tabella di decodifica tipo di località (senza geometria)

L'elenco delle tabelle e dei relativi campi disponibili può essere visualizzato in ogni momento premendo il bottone "elenco tabelle"

Esercitazione SELECT

--ES 3
SELECT 
  t1.descrizione, localita.id, 
  localita.geom, localita.loc2011, 
  localita.pro_com, localita.tipo_loc, 
  localita.popres as residenti, localita.denominazi as nome
FROM 
  localita, 
  tipo_localita as t1
WHERE 
  localita.tipo_loc = t1.tipo_loc
  --and popres > 20000
LIMIT 10

Indicare i campi delle tavole che si desidera selezionare, è possibile indicare anche un ALIAS. Elenco delle tavole, anche in questo caso è possibile indicare un ALIAS. Indicare la clausola di JOIN tra le due tabelle.

Durante la fase di test è consigliato impostare un limite massimo di risultati con la clausola: Limit
ad esempio: Limit 10

FUNZIONI GEOSPAZIALI

Funzioni geometriche di Spatialite

L'estensione spaziale di SQlite è Spatialite, e dispone di una serie molto completa di funzioni per la manipolazione di informazioni geografiche: Reference List

Le funzioni basilari che utilizzeremo nel workshop sono le seguenti

condizioni topologiche (condizioni spaziali): dentro, fuori, contenuto da, coincidente etc...
manipolazione geometrica (aggregazione): unione, differenza, intersezione, buffer, generalizzazione
misurazione geometrica: lunghezza, area, distanza
estrazione geometrica: centroidi, linee di unione
trasformazioni: riproiezioni, conversioni di formati geometrici

Il campo geometria

I tipi geometrici possono essere i seguenti:

geometria puntuale (POINT e MULTIPOINT)
geometria lineare (LINESTRING e MULTILINESTRING)
geometria poligonale (POLYGON e MULTIPOLYGON)

L'informazione spaziale è contenuta in almeno un campo geometrico. Nelle tabelle del workshop il campo geografico è per convenzione il campo "geom", ma nulla vieta che possa avere un altro nome.
Il formato di memorizzazione è binario, e quindi non visualizzabile "per gli umani" ma disponiamo di funzioni che permettono di trasformarlo in un formato più amichevole:

SELECT
	geom,
  ST_AsText(geom) AS wkt,
  Hex(ST_AsBinary(geom)) AS wkb
FROM localita

Esercitazione 4

--ES 4
SELECT 
  l.id, 
  ST_Centroid(l.geom) as geom,  asgeojson(ST_Centroid(l.geom)) as geojson,
  l.loc2011 as codice_localita,   l.pro_com as codice_comune,   l.tipo_loc, 
  tl.descrizione as desc_tipo_localita,     l.denominazi as nome_localita,
  l.popres as popolazione_residente,  l.maschi as maschi,
  l.popres - l.maschi as femmine
FROM 
  localita l,   tipo_localita tl
WHERE 
  --localita.tipo_loc = t1.tipo_loc;
  cast(l.tipo_loc as integer) = cast(tl.tipo_loc as integer)
LIMIT 10;

In questa query utilizziamo la funzione ST_Centroid per trasformare un elemento poligonale in puntuale. Esiste anche la funzione ST_Pointonsurface che restituisce un punto interno alla geometria.

VISUALIZZAZIONE

campi calcolati ed alias

select 
 id, 
 geom,
 st_astext(geom) as wkt,
 asgeojson(geom) as geojson
from 
  localita

ESERCITAZIONE 6 - Le subquery

select 
 c.id, c.comune, l.nome_localita, l.numero_famiglie, l.numero_abitazioni
from 
 comuni c, 
 (
	SELECT 
	  l.id, 
	  ST_Centroid(l.geom) as geom, asgeojson(ST_Centroid(l.geom)) as geojson,
	  l.loc2011 as codice_localita, l.pro_com as codice_comune, l.tipo_loc, tl.descrizione as desc_tipo_localita, 
	  l.denominazi as nome_localita,
	  l.popres as popolazione_residente, l.maschi as maschi, l.popres - l.maschi as femmine,
	  l.famiglie as numero_famiglie, l.abitazioni as numero_abitazioni, l.edifici as numero_edifici, l.shape_area as superficie_in_mq
	FROM 
	  localita l, tipo_localita tl
	WHERE 
	  --localita.tipo_loc = t1.tipo_loc;
	  cast(l.tipo_loc as integer) = cast(tl.tipo_loc as integer)
 ) as l
where 
 st_intersects(l.geom,c.geom)
 and c.comune like 'Torino';

In questa query facciamo dei calcoli: deriviamo il numero di femmine come sottrazione tra la popolazione totale e i maschi.

ESERCITAZIONE 7 - test geospaziali

--ES 7
SELECT 
 c.comune, count(*), AsGeoJSON(c.geom) as geojson
FROM 
  comuni AS c, localita AS l
WHERE 
c.cod_pro = 5 and
st_intersects(st_envelope(c.geom),st_envelope(l.geom)) and
st_intersects(c.geom,st_centroid(l.geom)) 
GROUP BY c.comune;

In questa query viene calcolato il numero di località presente in ogni comune. Dato che l'operazione viene fatta sulla combinazione dei dati di due tabelle (record tabella1 x record tabella2), il numero di dati su cui si va ad operare è enorme, quindi è opportuno restringere il campo di ricerca inserendo delle condizioni in WHERE

In tutti i database la valutazione delle condizioni WHERE si dice che è "lazy" ovvero alla prima condizione false si tralascia la valutazione delle restanti condizioni

Quindi prima restringiamo la ricerca ad un codica provincia, poi valutiamo se le estensioni delle geometrie coincidono con st_envelope e solo allora andremo ad effettuare il test sull'intersezione tra geometrie con st_intersects che è una funzione molto pesante dal punto di vista computazionale

Tipi geometrici

Visualizzazione in mappa

Il risultato della query può essere visualizzato in mappa quando la tabella di output contiene un campo "geojson" contente le geometrie da visualizzare in formato geojson:

Il formato geojson è l'espressione testuale di un'oggetto javascript prodotta secondo lo standard "json" ed è facilmente "digeribile" dalle applicazioni web

L'altra condizione perchè la geometria possa essere visualizzata correttamente nel visualizzatore è che il dato geografico deve essere espresso in gradi di latitudine e longitudine. Le tabelle predefinite sono nel formato ETRS89-epsg4258, sono quindi compatibili con il visualizzatore

SELECT
  AsGeoJSON(
    ST_GeomFromText('LINESTRING(7 45,8 45,8 44,7 44, 7 45)')
  ) AS geojson

Condizioni topologiche

SELECT 
	comuni.comune,
	AsGeoJSON(geom) as geojson
FROM comuni
WHERE ST_Intersects(geom,ST_GeomFromText('POLYGON((7 45,8 45,8 44,7 44, 7 45))'))

AGGREGAZIONE SPAZIALE

Le funzioni spaziali possono essere usate per collezionare le geometrie aggregate da una query di raggruppamento

SELECT
	cod_pro,
  COUNT(comune) as num_comuni,
  SUM(shape_area)/10000 as area_ha,
  AVG(shape_area)/10000 as area_media,
  AsGeoJSON(ST_union(geom)) as geojson
FROM comuni
GROUP BY cod_pro

L'intersezione tra geometrie lineari genera geometrie puntuali:

SELECT
  asgeojson(st_intersection(viab_princ.geom,padanasuperiore.geom)) as geojson,
  viab_princ.id,
  padanasuperiore.id
FROM viab_princ, padanasuperiore
WHERE st_intersects(viab_princ.geom,padanasuperiore.geom)

MISURE, TRASFORMAZIONI E BUFFER

E' possibile acquisire dati dimensionali delle geometrie inserite lunghezze, aree, distanze

SELECT comune,St_area(geom) AS superficie, St_length(geom) AS perimetro FROM comuni

SELECT id,St_length(geom) AS lunghezza FROM padanasuperiore

SELECT id,St_distance(geom,ST_GeomFromText('POINT(7 45)')) AS dist_pd FROM centri

Essendo le geometrie espresse in latitudine e longitudine, I dati geometrici sono espressi in gradi. E' quindi necessario trasformare le coordinate in un sistema di riferimento proiettato per ottenere le corrette dimensioni in metri

SELECT id,name,St_length(st_transform(geom,3003)) AS lunghezza FROM viab_princ
limit 10

La funzione di buffer genera un poligono i cui vertici distano ad una distanza costante da una geometria generatrice

SELECT  asgeojson(st_buffer(viab_princ.geom,0.03)) as geojson FROM viab_princ
limit 10

ELABORAZIONI

Quali sono i centri serviti dalla ferrovia?

Individuazione di un buffer di 10 km dalla ferrovia

SELECT AsGeoJson(st_buffer(st_union(geom), 0.03)) AS geojson
FROM padanasuperiore

individuazione dei comuni ricadenti nel buffer appena individuato

SELECT *,ASgeoJSON(geom) as geojson
FROM comuni
WHERE st_intersects(geom, (
    SELECT st_buffer(st_union(geom), 0.03)
    FROM padanasuperiore
))

metodo alternativo meno performante ma che recupera i dati

SELECT comuni.comune,st_distance(comuni.geom,padanasuperiore.geom),ASgeoJSON(comuni.geom) as geojson
FROM comuni JOIN padanasuperiore 
	ON st_intersects(comuni.geom, st_buffer(padanasuperiore.geom, 0.03))
GROUP BY comuni.comune

NB: i dati geografici sono espressi in latitudine e longitudine quindi la distanza di buffer va espressa in gradi (0.01 gradi corrispondono circa a 950 m)

Qual'è la popolazione dell'area metropolitana di Torino?

Individuazione dei comuni limitrofi

SELECT comune, asgeojson(geom) as geojson
FROM comuni 
WHERE st_touches(geom, (
  SELECT geom
  FROM comuni
  WHERE comune = 'Torino'
))

Usando st_intersects al posto di st_touches si include anche il comune stesso di Torino

SELECT comune, asgeojson(geom) as geojson
FROM comuni 
WHERE st_intersects(geom, (
  SELECT geom
  FROM comuni
  WHERE comune = 'Torino'
))

Quanto è lungo Corso Francia?

Individuazione del fiume: Trattandosi di un grafo idrologico, il fiume è individuato da una serie di archi:

SELECT  name,asgeojson(geom) as geojson 
FROM viab_princ 
WHERE name = 'Corso Francia'

Metodo 1: quindi si può raggruppare sommando la lunghezza di ugni arco

SELECT  
	name, 		
    st_length(st_transform(st_union(geom),3003)),
    asgeojson(st_union(geom)) as geojson 
FROM viab_princ 
WHERE name = 'Corso Francia'
GROUP BY name

Nota bene: Questa query non funzionerebbe in postgis poichè nel caso di raggruppamento Postgresql richiede di specificare una funzione di aggregazione per ogni campo

QGIS

Spatialite in QGIS

E' possibile replicare in ambiente QGIS le esercitazioni del workshop. In tal caso è necessario:

Scaricare in locale il database dell'esercitazione: veneto.sqlite
Stabilire una connessione con il database
Aprire DB MANAGER
Incollare le query nella finestra SQL ed aggiungerle in mappa

Virtual layers in QGIS

Un'altro modo per usare in modo produttivo le query geospaziali in QGIS sono i Virtual Layers:

I virtual Layers sono una caratteristica poco nota ma molto potente di QGIS, e permettono di creare relazioni SQL tra i layer correntemente caricati da qualunque datasource. Il virtual layer viene rigenerato in modo automatico e trasparente ogni qualvolta vengono aggiornati i layers collegati al layer virtuale:

La funzionalità è resa possibile replicando i layers in un database Spatialite locale gestito direttamente da QGIS, dunque la sintassi da utilizzare nelle query è quella di Spatialite. Su questa funzionalità è disponibile un breve ma chiaro tutorial di Totò Fiandaca