QGIS e lo spatial join condizionato

Lo spatial join è una comunissima operazione in ambito GIS che permette di trasferire attributi da un layer ad un altro in base alle loro relazioni spaziali, questa operazione viene usata spesso per rispondere alla domanda: quali geometrie sono più vicine ad altre? negli strumenti di QGIS, questa operazione lavora prendendo come parametro solo la geometria, non ci sono algoritmi o funzioni che permettono di prendere in considerazioni altri parametri oltre alla geometria, come ad esempio un attributo alfanumerico, quindi è impossibile rispondere alla domanda: quali geometrie sono più vicine ad altre ma che hanno un determinato valore di attributo?

Per risolvere questo problema occorre, quindi, scrivere una espressione complessa o usare le query spaziali.

Vediamo come risolvere lo spatial join condizionato usando le espressioni di QGIS.

Il quesito è preso da una domanda posta su un gruppo Facebook:

Sintetizzando: occorre acchiappare il punto (del layer2) più vicino (al layer1) ricadenti nello stesso poligono (grid) o che abbia lo stesso valore di un attributo.

l punto 762 è collegato correttamente al punto 6173, che è il punto più vicino rimanendo dentro il poligono; il punto 4219 è il punto più vicino geometricamente, ma non è nello stesso poligono
Ogni poligono può contenere uno o più punti del layer1 e del layer2, quindi possono verificarsi questi casi descritti sopra: due punti layer1 collegati allo stesso punto del layer2.

Un membro del gruppo, Andrea Giudiceandrea (grazie), ha realizzato un dataset di prova che condivido anche qui:

i layer (layer1 e layer2), per facilitare i test, contengono già ID univoco (IDp) del poligono che li contiene, quindi possiamo concentrarci solo nei due layer (questo IDp è l’attributo da prendere in considerazione assieme alla geometria).

Prima soluzione di Andrea ( >=QGIS 3.16 )

attribute( 
array_filter( 
overlay_nearest('layer2',$currentfeature,
                           limit:=-1,max_distance:=2000), 
attribute( @element, 'IDp' ) = "IDp" )[0], 'IDl2' )

Seconda soluzione di Totò ( >=QGIS 3.18 )

relation_aggregate( 
	relation:='rel',
	aggregate:='array_agg',
	expression:="IDl2")
	[
		with_variable('fufu',
		array_foreach(
			relation_aggregate( 
				relation:='rel',
				aggregate:='array_agg',
				expression:=$geometry),
		distance($geometry,@element)),
		array_find(@fufu,array_min(@fufu)))
	]

entrambe le soluzioni risolvono il problema, ma la complessità viene pagata tramite un tempo di elaborazione lungo su grossi dataset (per il dataset allegato e nel mio PC [Dell G5]: circa 7 e 3 sec, rispettivamente la prima e seconda soluzione); sotto un esempio della tabella degli attributi del layer1:

Osservazione sulla prima soluzione:

La funzione overlay_nearest() dovrebbe essere utilizzata senza il parametro max_distance:=2000 altrimenti risulta non generalizzabile e dipendente dalle dimensione dei poligoni, ma in questo caso impiegherebbe troppo tempo (circa 100 sec), quindi occorre valutare il valore in funzione delle dimensioni dei poligoni utilizzati, nel nostro caso 2000 metri vanno bene, ma solo perché abbiamo volutamente creato esagoni regolari.

Osservazione sulla seconda soluzione:

Occorre creare una relazione di progetto (rel) e, inoltre, nessun punto deve cadere fuori dai poligoni altrimenti l’espressione segnala un errore; per evitare l’errore occorre modificare l’espressione in questo modo:

if ("IDp" is NULL, NULL,
     with_variable('related_features',
           relation_aggregate(
           relation:='rel',
           aggregate:='array_agg',
           expression:="IDl2"),
           if(array_length(@related_features)<1,
               NULL, @related_features
  [
            with_variable('fufu',
                 array_foreach(
                 relation_aggregate(
                 relation:='rel',
                 aggregate:='array_agg',
                 expression:=$geometry),
           distance($geometry,@element)),
           array_find(@fufu,array_min(@fufu)))
  ])))

Per far tracciare il segmento che unisce i punti dei due layer possiamo creare un nuovo layer (come Generatore Geometria) nel tema del layer1:

e usare l’espressione:

make_line($geometry,
geometry(get_feature(
		'layer2',
		'IDl2', 
		attribute($currentfeature,'toto')))
		)

Soluzione tramite query spaziale ( >=SQLite 3.33)

--
-- traccia un segmento che collega i punti 
-- del layer1 con quelli del layer2
--
SELECT l1."IDl1" AS la1, 
               shortestline(l1."geom", l2."geom") as geom
FROM "layer1" l1 , "layer2" l2 
WHERE l1."IDp"=l2."IDp"
GROUP BY l1."IDp"
HAVING min(st_distance(l1."geom", l2."geom"))
--
-- per aggiornare il campo IDl2 nel layer1
--
UPDATE layer1 SET IDl2="IDl2" 
FROM 
(SELECT l2."IDl2", l2."IDp"
FROM"layer1" l1 , "layer2" l2 
WHERE l1."IDp"=l2."IDp"
GROUP BY l1."IDp"
HAVING min(st_distance(l1."geom", l2."geom"))) as t
WHERE layer1."IDp"=t."IDp"

l’output delle query è velocissimo anche su grossi dataset.


NOTE FINALI: Le query spaziali SQL risolvono brillantemente questi quesiti, sia in termini di facilità di scrittura della query che di tempo. Per fortuna QGIS ha la possibilità di creare layer a partire da query SQL usando vari dialetti: SQLite/SpatiaLite, PostgrSQL/PostGIS e altri.

[30-10-2021] Andrea Giudiceandrea ha proposto una patch per contemplare il filtro dinamico nella funzioni overlay_nearest per risolvere questi casi.
[31-10-2021] La patch è stata accettata e unita alla master e verrà fatto anche un backport nella 3.22 e forse anche nella 3.16.


RIFERIMENTI

RINGRAZIAMENTI


Approfondimenti

Per evitare di creare una relazione di progetto (per utilizzare la mia soluzione), potremmo creare una relazione al volo e usarla direttamente nella espressione, eccola:

attribute(
with_variable('relation',
aggregate( 
layer:='layer2',
aggregate:='array_agg',
expression:= get_feature_by_id('layer2',$id),
filter:= "IDp" = attribute(@parent,'IDp')),
@relation
	[
		with_variable('fufu',
		array_foreach(@relation,
		distance($geometry,geometry(@element))),
		array_find(@fufu,array_min(@fufu)))
	]),
'IDl2')

Questa espressione è molto più lenta (da 3 sec passa a 45) rispetto al caso della creazione di una vera relazione di progetto, ma la reputo altamente didattica.


Per utilizzare l’espressione con la relazione a partire da QGIS 3.16:

relation_aggregate( 
	relation:='rel',
	aggregate:='array_agg',
	expression:="IDl2")
	[
		with_variable('fufu',
		array_foreach(
			relation_aggregate( 
				relation:='rel',
				aggregate:='array_agg',
				expression:=$geometry),
		distance($geometry,@element)),
		array_find(@fufu,array_sort(@fufu)[0]))
	]

Osservazione: la funzione array_sort permette di sostituire array_min introdotta in QGIS 3.18.


Nuovissima patch appena introdotta nella master:

eval('overlay_nearest( \'layer2\', 
           "IDl2", 
            filter:=IDp=' || "IDp" || ')')[0]

questa espressione impiega meno di tre secondi per gli stessi dati.

In breve: la nuova patch permette di utilizzare un filtro (filter) dinamico che cambia per ogni caratteristica valutata del livello sorgente, cioè: eval('overlay_nearest(\'target_layer\',ID_target_layer,filter:=field_target_layer='||"field_source_layer"||')').

Un pensiero su “QGIS e lo spatial join condizionato

Rispondi

Inserisci i tuoi dati qui sotto o clicca su un'icona per effettuare l'accesso:

Logo di WordPress.com

Stai commentando usando il tuo account WordPress.com. Chiudi sessione /  Modifica )

Google photo

Stai commentando usando il tuo account Google. Chiudi sessione /  Modifica )

Foto Twitter

Stai commentando usando il tuo account Twitter. Chiudi sessione /  Modifica )

Foto di Facebook

Stai commentando usando il tuo account Facebook. Chiudi sessione /  Modifica )

Connessione a %s...

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