ML GFOSS.it : una lezione su spatialite

In questo articolo voglio far vedere la bellezza dell’Open Source e della comunità che ci sta dietro.

Ecco una thread in ML GFOSS.it e la relativa risposta:

OGGETTO: Spatial Join con SpatiaLite

Buongiorno e buona domenica a tutti.
Sto provando uno spatial join sfruttando un buffer a 250 m da un punto, ma non riesco a capire perchè non funziona.

Come dato di base sto sfruttando le celle censuarie ISTAT che potete scaricare da qui
http://www.istat.it/storage/cartografia/basi_territoriali/WGS_84_UTM/2011/R15_11_WGS84.zip,
le ho chiamate test_celle.  Poi c’è un punto che ha come attributo il solo ID e di cui faccio (vorrei fare) un buffer a 250 m. Il tutto voglio farlo confluire in una view chiamata spatial_join.
Inserisco il codice che ho usato:


CREATE VIEW spatial_join AS
SELECT *
FROM test_celle AS target, (SELECT ST_BUFFER(geom, 250) AS geom FROM input) AS buffer
WHERE ST_INTERSECTS (target.geom, buffer.geom);

view raw

query_mm_1.sql

hosted with ❤ by GitHub

Dove sbaglio?

RISPOSTA:

(Alessandro Furieri papà di SpatiaLite)

ciao Massimiliano,

scusami per la risposta poco tempestiva ma nei giorni scorsi ero impegnato su tutt’altri versanti. Ho provato a replicare il tuo esempio a partire dalle Sezioni di Censimento ISTAT 2011 della Campania. Come punto di riferimento (visto che tu non hai fornito esempi) ho usato la posizione della Cappella San Severo nel cuore del centro storico di Napoli. Piu’ sotto troverai i miei commenti.


CREATE VIEW spatial_join AS
SELECT *
FROM test_celle AS target, (SELECT ST_BUFFER(geom, 250) AS geom FROM input) AS buffer
WHERE ST_INTERSECTS (target.geom, buffer.geom);

view raw

query_mm_1.sql

hosted with ❤ by GitHub

primo errore: per arrivare a creare una View non si deve mai   partire creando direttamente la View stessa, perche’ in questo   modo se qualcosa va storto non capirai piu’ nulla e non sarai   in grado di identificare e correggere i tuoi errori. Si parte sempre definendo (e testando accuratamente) lo  statement SELECT, che una volta messo opportunamente a punto sara’ poi facilissimo trasformare in una View.

secondo errore: mai usare “SELECT *” in una View; tu pensi in  questo modo di risparmiare tempo e fatica, ma finirari invece per perdere un sacco di tempo e per faticare un sacco rincorrendo errori “misteriosi”. Buona regola da usare _SEMPRE_ con le View: tutte le colonne vanno identificate esplicitamente con il proprio nome qualificato  (specificando cioe’ a quale tavola appartiene ciascuna colonna), ed e’ sempre opportuno specificare esplicitamente tutti gli alias names di colonna per evitare poi spiacevoli sorprese al momento dell’esecuzione.

terzo errore: (vedo che lo segnali tu stesso nella tua mail successiva) ST_Contains() non e’ l’operatore spaziale che devi usare, perche’ richiede che la prima geometria copra interamente la seconda. Ma nel nostro caso e’ molto piu’ verosimile che le Sezioni di Censimento ricadano solo parzialmente all’interno del Buffer, ragion per cui e’ molto  piu’ opportuno usare la ST_Intersects()

quarto errore (veniale): quando usi un alias-name per le tavole, cerca di usare una singola lettera, altrimenti poi ti troverai a dovere gestire nomi completi esageratamente  lunghi ed inutilmente complicati.

ed ecco qua la prima versione della nostra SELECT riscritta in modo piu’ ragionevole:


SELECT t.pk_uid AS rowid, t.cod_reg AS cod_reg, t.cod_istat AS cod_istat,
       t.pro_com AS pro_com, t.sez2011 AS sez201, t.sez AS sez,
       t.loc2011 AS loc2011, t.cod_loc AS cod_loc, t.geometry AS geom
FROM test_celle AS t,
     (SELECT ST_Buffer(geom, 250) AS geom FROM input) AS b
WHERE ST_Intersects (t.geometry, b.geom) = 1;

view raw

query_mm_2.sql

hosted with ❤ by GitHub

questa prima query identifica 46 sezioni di censimento in circa 1 secondo e 460 millisecondi; possiamo migliorare la velocita’ di esecuzione ?


SELECT t.pk_uid AS rowid, t.cod_reg AS cod_reg, t.cod_istat AS cod_istat,
       t.pro_com AS pro_com, t.sez2011 AS sez201, t.sez AS sez,
       t.loc2011 AS loc2011, t.cod_loc AS cod_loc, t.geometry AS geom
FROM test_celle AS t,
     (SELECT geom FROM input) AS b
WHERE ST_Distance (t.geometry, b.geom) <= 250;

view raw

query_mm_3.sql

hosted with ❤ by GitHub

basta evitare di usare la ST_Buffer + ST_Intersects, che nel caso di un singolo punto rappresenta un’inutile complicazione che impone un sacco di calcoli complessi di cui possiamo fare tranquillamente a meno. Basta usare semplicemente la ST_Distance per ottenere il medesimo risulato con tempi di esecuzione decisamente migliori. Questa seconda query impiega soli 600 millis, cioe’ meno della meta’ del tempo richiesto dalla precedente. Comunque possiamo arrivare ad una ottimizzazione ancora piu’ spinta.

SELECT CreateSpatialIndex(‘test_celle’, ‘geometry’);


SELECT t.pk_uid AS rowid, t.cod_reg AS cod_reg, t.cod_istat AS cod_istat,
       t.pro_com AS pro_com, t.sez2011 AS sez201, t.sez AS sez,
       t.loc2011 AS loc2011, t.cod_loc AS cod_loc, t.geometry AS geom
FROM test_celle AS t,
     (SELECT geom FROM input) AS b
WHERE ST_Distance (t.geometry, b.geom) <= 250 AND t.ROWID IN
     (SELECT rowid FROM SpatialIndex
      WHERE f_table_name = 'test_celle' AND search_frame =
            BuildCircleMbr(ST_X(b.geom), ST_Y(b.geom), 250));

view raw

query_mm_4.sql

hosted with ❤ by GitHub

per prima cosa andiamo a creare uno Spatial Index che supporti le sezioni di censimento. e poi riscriviamo la query precedente in modo tale che ora prenda in considerazione lo Spatial Index. Ora siamo scesi ad un tempo di esecuzione di soli 250 millis, che e’ decisamente soddisfacente.

Nota bene: il dataset sezioni di censimento campania contiene circa 25mila features, cioe’ un dataset MEDIO-PICCOLO; ed infatti i tempi di esecuzione sono abbastanza soddisfacenti anche quando non si usa lo Spatial Index. Ma se si fosse trattato di un dataset GRANDE (con milioni di features) l’uso dello Spatial Index o meno avrebbe
prodotto effetti decisamente critici.

a questo punto siamo finalmente pronti per andare a creare la nostra View:


CREATE VIEW spatial_join AS
SELECT t.pk_uid AS rowid, t.cod_reg AS cod_reg, t.cod_istat AS cod_istat,
       t.pro_com AS pro_com, t.sez2011 AS sez201, t.sez AS sez,
       t.loc2011 AS loc2011, t.cod_loc AS cod_loc, t.geometry AS geom
FROM test_celle AS t,
     (SELECT geom FROM input) AS b
WHERE ST_Distance (t.geometry, b.geom) <= 250 AND t.ROWID IN
     (SELECT rowid FROM SpatialIndex
      WHERE f_table_name = 'test_celle' AND search_frame =
            BuildCircleMbr(ST_X(b.geom), ST_Y(b.geom), 250));

view raw

query_mm_5.sql

hosted with ❤ by GitHub

SELECT * FROM spatial_join;

ok, ora abbiamo creato la View ed abbiamo verificato che funzioni correttamente.

Attenzione: non abbiamo ancora finito, perche’ per ora abbiamo semplicemente una View “non qualificata”, che il data provider di QGIS non sara’ in grado di visualizzare. Quello che dobbiamo fare a questo punto e’ registrare correttamente una Spatial View, che anche QGIS sapra’ riconoscere correttamente come un layer di tipo vettoriale.


INSERT INTO views_geometry_columns VALUES
('spatial_join', 'geom', 'rowid', 'test_celle', 'geometry', 1);

view raw

query_mm_6.sql

hosted with ❤ by GitHub

in cui:
– ‘spatial_join‘ e’ il nome della View che intendi registrare per promuoverla a Spatial View;

– ‘geom‘ e’ l’alias-name della geometria presente nella View

– ‘rowid‘ identifica la Primary Key corrispondente alla geometria (a questo punto ti dovrebbe risultare assolutamente chiaro perche’ abbiamo dovuto usare un sacco di noiosi alias-names nella definizione della View);

– ‘test_celle‘ e’ il nome della tavola che fornisce le geometrie

– ‘geometry‘ e’ il nome della colonna geometria dentro a test_celle (specificare queste ultime tre informazioni e’ indispensabile per consentire al data provider di QGIS di accedere correttamente agli Spatial Index anche nel caso delle Spatial Views).

– infine 1 significa semplicemente che questa Spatial View e’ di tipo read-only (puoi consultare le features, ma non le puoi modificare in nessun caso).

ora abbiamo veramente finito; fai partire QGIS e vedrai che funziona tutto 😉

ciao Sandro


Note finali: questa non è una semplice risposta ad un quesito ma è una piccola lezione sull’uso dei DBMS e, più in generale, nell’ottimizzazione di uno script SQL. Non è la prima volta che Furieri, nel rispondere a dei quesiti, fa una bella lezione.

W Open Source, GFOSS.it e A. Furieri.

Iscrivetevi alla mailing list GFOSS.it  e all’associazione GFOSS.it


 

8 pensieri su “ML GFOSS.it : una lezione su spatialite

  1. Complimenti per gli articoli e video su spatialite!
    molto pratici ed accessibili anche per coloro che sono alle prime armi.
    Un argomento interessante sono anche i Triggers!

    "Mi piace"

  2. Ciao Salvatore, sto cercando di fare una intersezione di una polilinea con dei poligoni. Utilizzando le views vorrei aggiornare il campo delle lunghezze dei vari tratti che ricadono dentro i diversi poligoni dei comuni dinamicamente. Quando faccio la select di intersezione con st_intersects la tabella risultante presenta tante righe quante sono le intersezioni: come posso farle diventare features distinte?
    Grazie in anticipo,
    Un caro saluto

    Piace a 1 persona

    1. Ciao, sperando di aver capito il problema:
      le view (con geometria) affinché possano essere viste come feature (geometrie con tutte le caratteristiche) occorre che siano registrate e questo è spiegato alla fine dell’articolo., a partire da: ATTENZIONE!!!

      "Mi piace"

      1. Il problema è che il layer generato con st_intersect contiene una feature unica anche se in tabella attributi ci sono tre righe, quante sono le intersezioni

        "Mi piace"

Lascia un commento

Questo sito utilizza Akismet per ridurre lo spam. Scopri come vengono elaborati i dati derivati dai commenti.