Come prendere ID di un layer sovrapposto con area maggiore

Dati due layer poligonali (per esempio le regioni italiane e il quadro di unione del DEM Tinitaly); assegnare, per ogni riquadro, il codice univoco (cod_reg) della regione che ha maggior estensione all’interno del riquadro stesso.

Per capire meglio, ecco uno screenshot di esempio:

all’interno del riquadro (rosso) ricade una porzione di territorio siciliano (cod_reg=19) e una porzione del territorio calabrese (cod_reg=18): assegnare al riquadro il cod_reg corrispondente alla porzione con area maggiore, in questo caso, la porzione di territorio siciliano è maggiore di quello calabrese.

Per versione di QGIS < 3.24 Tisler:

Il quesito è risolto usando le potenti espressioni di QGIS:

array_reverse(
        aggregate(
              layer:='Reg01012022_G_Wgs84', 
              aggregate:='array_agg',
              expression:="cod_reg", 
              filter:=intersects( $geometry,geometry(@parent)),
              order_by:=area(intersection(geometry(@parent),$geometry )))
                       )[0]

Aggrego in un Array tutti i cod_reg intersecantesi per ogni riquadro, li ordino in funzione dell’area e prendo quella con estensione maggiore.

Note aggiuntive:

  • l’opzione order_by ordina, by default, in senso ascendente e quindi occorre invertire il senso usando array_reverse, oppure, senza usare array_reverse ma estraendo l’ultimo valore dell’array sostituendo [0] con [-1].

Per versione di QGIS >= 3.24 Tisler:

overlay_intersects(
           'Reg01012022_G_Wgs84',
           "cod_reg",
            sort_by_intersection_size:='des')[0] -- sort_by_intersection_size Nuovo parametro

Utilizzando le funzioni di Overlay_* e il nuovo parametro sort_by_intersection_size che semplifica enormemente la ricerca dei valori (circa tre volte più veloce!!!).

Oppure con una query SQL in un database SpatiaLite/PostGIS:

SELECT
  t.pk AS pk,
  t.geom AS geom,
  t.max_area_regione AS COD_REG
FROM
  (
    SELECT
      re.pk,
      re.geom,
      FIRST_VALUE(COD_REG) OVER (
        PARTITION BY re.pk
        ORDER BY
          ST_area(ST_intersection(re.geom, ri.geom)) DESC
      ) AS max_area_regione
    FROM
      regioni_istat ri
      JOIN reticolo re ON st_intersects(ri.geom, re.geom)
  ) t
GROUP BY
  1,
  2

Nella query si fa uso delle Windows Function, studiate appositamente per risolvere questi casi.

NOTE FINALI: Le espressioni di QGIS sono sempre più potenti e permettono di risolvere molti quesiti sia alfanumerici che spaziali, quest’ultimo parametro (sort_by_intersection_size) introdotto di recente nella funzione di Overlay_intersects semplifica notevolmente i tempi di elaborazione in quanto permette di ordinare i risultati dell’Array in funzione della dimensione dell’intersezione.

RIFERIMENTI

RINGRAZIAMENTI

  • Andrea Borruso per aver dedicato il suo tempo alla risoluzione di un problema legato agli Array.

Dati e Progetto


I MIEI CANALI – ISCRIVITI


Se il blog post Ti è piaciuto cliccate su ‘Mi piace’, grazie!!!
if you liked the blog post click on ‘Like’, thank you !!!

SE IL POST/BLOG TI È STATO UTILE CONTRIBUISCI A MANTENERLO AGGIORNATO PAYPAL


Pubblicità

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 )

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.