On Sun, 12 Nov 2017 11:14:00 +0100, Massimiliano Moraca wrote:
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 250m.
Il tutto voglio farlo confluire in una view chiamata *spatial_join*.


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.


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_CONTAINS (target.geom, buffer.geom)


- 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;


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;


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 milis, 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));


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));

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);

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

_______________________________________________
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e' una lista di discussione pubblica aperta a tutti.
I messaggi di questa lista non hanno relazione diretta con le posizioni 
dell'Associazione GFOSS.it.
801 iscritti al 19/07/2017

Rispondere a