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.

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

• QGIS : https://www.qgis.org/it/site/
• DATI : https://drive.google.com/file/d/1CJaX-eZ2sAqtUxw445S3-7-lxFlAXKFg/view
• Guida HfcQGIS : https://hfcqgis.opendatasicilia.it/
• overlay_nearest: https://hfcqgis.opendatasicilia.it/gr_funzioni/geometria/geometria_unico_con_troppi_colori/#overlay_nearest
• SQLite : https://sqlite.org/lang_update.html#upfrom
• Caso simile su Stack Exchange: https://gis.stackexchange.com/questions/391120/qgis-expression-with-overlay-fuction-filter-condition-based-on-comparison-of-at

RINGRAZIAMENTI

• Andrea Giudiceandrea
• Gruppo Facebook

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"||')').