Elenco ordinato di Comuni attraversati da un itinerario Domodossola-Aosta.

Pubblicato il da Totò

Ieri nel gruppo Telegram di QGIS Italia, il vice-presidente di GFOSS.it (Stefano Campus) pone questa domanda molto interessante:

qualcuno risponde:

Altri membri del gruppo rispondono con varie idee, in alcuni casi buone idee ma nulla di concreto. Dopo qualche ora presento la mia soluzione usando solo le espressioni e il calcolatore di campi.

Come nasce una espressione complessa

Per risolvere qualsiasi tipo di problema occorre una idea di base da tradurre, poi, in espressione utilizzando le varie funzioni messe a disposizione dal field calc.

Nel gruppo ci sono due suggerimenti:

L’idea è quella di fare una intersezione tra la linea che rappresenta l’itinerario e i poligoni comunali.

Step 1 : ottenere le linee di intersezione tra itinerario e poligoni

intersection($geometry, geometry(get_feature('percorso','fid',1)))

utilizzo l’espressione di sopra e genero una nuova geometria:

Step 2: per ogni linea ricadente nel poligono comunale, calcolo il baricentro creando una nuova geometria puntuale a partire dall’espressione di sopra (che rappresenta la geometria dell’ itinerari spezzato per ogni comune):

centroid(intersection($geometry, geometry(get_feature('percorso','fid',1))))

ma i centroid NON cadono lungo le linee e questo potrebbe creare problemi, per ovviare uso altra funzione :

line_interpolate_point(intersection($geometry, geometry(get_feature('percorso','fid',1))),
length(intersection($geometry, geometry(get_feature('percorso','fid',1)))))

Step 3: per ottenere la lunghezza del percorso, rispetto l’origine dell’itinerario e per ogni punto baricentro dell’area poligonale, occorre utilizzare altra funzione per popolare il campo lunghezza (sempre sul layer poligonale):

line_locate_point(geometry(get_feature('percorso','fid',1)),
line_interpolate_point(intersection($geometry, geometry(get_feature('percorso','fid',1))),
length(intersection($geometry, geometry(get_feature('percorso','fid',1))))/2))

e qui trovo un primo problema (segnalato da Stefano), cioè se l’itinerario parte da un Comune, passa in altro Comune e rientra nel Comune di partenza (vedi screenshot di sotto), si manifesta un problema: il secondo Comune viene conteggiato come se fosse il primo:

questo perché il centroide della linea del secondo Comune ha una lunghezza inferiore, dall’origine dell’itinerario, rispetto al primo Comune.

Per risolvere il problema, non posso utilizzare il centroide ma un altro punto lungo la linea. La funzione line_locate_point permette di esprimere una lunghezza in termini percentuali e quindi utilizzo il 10% o il 90% a seconda dell’origine dell’itinerario.

L’espressione finale è:

with_variable('toto',  reverse( geometry(get_feature('percorso','fid',1))),
    array_find( 
        array_filter(
            array_sort(
                array_agg(
                    line_locate_point(@toto, line_interpolate_point( intersection($geometry,@toto ), 
                    length(intersection($geometry,@toto))*0.1 ))
                         )
                      ), @element >0),
                    line_locate_point(@toto, line_interpolate_point( intersection($geometry,@toto ), 
                    length(intersection($geometry,@toto)) *0.1 ))
                )
)

dove ho creato una variabile @toto che equivale a un pezzo di espressione: geometry(get_feature('percorso','fid',1)).

Un piccolo cenno agli Array utilizzati nell’espressione:

  • array_agg permette di aggregare tutti i valori in un unico array;
  • array_sort permette di ordinare i valori;
  • array_filter filtra i valori in modo da lavorare solo con lunghezze > 0;
  • array_find permette di usare i numeri da 0 a n al posto delle lunghezze che possono essere numero grossi.

I Comuni sono stati ordinati a partire da Domodossola, ma se volessi farlo al contrario, basterebbe sostituire due solo valori nell’espressione:

  1. usare la funzione reverse che inverte il senso di marcia dell’itinerario;
  2. utilizzare il 90% della lunghezza al posto del 10%.(length(intersection($geometry,@toto)) *0.9)

NOTE FINALI: questa è l’ennesima prova che dimostra la potenza delle espressioni e del calcolatore di campi di QGIS.

Download dati e progetto in formato gpkg (per gentile concessione di Stefano)

Riferiementi

Ringraziamenti

Per chi volesse approfondire l’uso del Field Calc:
http://hfcqgis.opendatasicilia.it/it/latest/corso_formazione/corso_di_formazione.html

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

SE IL POST VI È STATO UTILE CONTRIBUITE A MANTENERLO AGGIORNATO PAYPAL

14 pensieri su “Elenco ordinato di Comuni attraversati da un itinerario Domodossola-Aosta.”

  2. Ciao Salvatore,
    la prima riga dell’espressione finale è indicata come:
    with_variable(‘toto’, reverse( geometry(get_feature(‘percorso’,’fid’,1))),

    ma non dovrebbe essere:
    with_variable(‘toto’, geometry(get_feature(‘percorso’,’fid’,1)),
    ?

    Piace a 1 persona

    Rispondi

  3. Ciao Salvatore,
    nel codice indicato c’è
    with_variable(‘toto’, reverse( geometry(get_feature(‘percorso’,’fid’,1))),
    e poi
    length(intersection($geometry,@toto))*0.1 ))

    Più sotto nel post c’è scritto che:
    se volessi farlo al contrario, basterebbe sostituire due solo valori nell’espressione:
    – usare la funzione reverse che inverte il senso di marcia dell’itinerario;
    – utilizzare il 90% della lunghezza al posto del 10%.(length(intersection($geometry,@toto)) *0.9)

    Quindi resto confuso… con la funzione reverse tu suggerisci di usare “*0.1” come nel codice o “*0.9” come è scritto in fondo al post?

    Piace a 1 persona

    Rispondi

    1. Rieccomi,
      0.1 o 0.9 sono due fattori che servono a posizionare il punto lungo la linea altrimenti il risultato finale sarebbe non corretto.
      Il reverse va usato assieme e a 0.1 o 0.9 a seconda della configurazione che hai tra i due layer.
      Ora non ricordo in dettaglio, ma fai delle prove e vedi se il risultato è corretto.

      "Mi piace"

      Rispondi

  4. Grazie per l’interessamento.

    Ho fatto alcune prove, ma non sono riuscito ad ottenere il risultato atteso…

    Nel GeoPackage allegato al tuo articolo, ci sono: un progetto QGIS, un layer lineare con tre percorsi e un layer poligonale con i territori comunali e i campi fid, comune_nom, sort_A-D e sort_D-A.

    Nel layer poligonale ci sono, tra gli altri, 64 record con valori da 0 a 63 nei due campi sort_A-D (da 0 Domodossola a 63 Carema) e sort_D-A (da 0 Carema a 63 Domodossola).
    Questi 64 record sono effettivamente i 64 territori comunali attraversati del percorso 1.

    Ho utilizzato la tua “espressione finale” (che contiene la funzione “reverse” e la moltiplicazione per 0.1) con il calcolatore di campi per creare un nuovo campo e verificare che il risultato fosse uguale ai valori già presenti.
    Tuttavia il nuovo campo è diverso da NULL solo per 55 comuni, invece che 64, con valori da 0 a 54 (da 0 Domodossola a 54 Carema): ci sono quindi 9 comuni che, benché attraversati dal percorso 1, non vengono conteggiati quali comuni attraversati.

    Ho provato anche le altre tre combinazioni (sempre con “reverse”, ma con 0.9; senza reverse, con 0.1 e con 0.9) ottenendo risultati diversi, rispettivamente: 48 comuni, 55 comuni (ma non tutti uguali a quelli del primo caso), 49 comuni.

    Piace a 1 persona

    Rispondi

    1. Ciao,
      ho riascaricato il geopackage e ho popolato nuovamente un nuovo campo utilizzando l’espressione finale e tutto funziona, ottengo valori tra 0 e 63.
      Purteoppo non posso sapere cosa sbagli, se puoi registra un video e fammi vedere cosa fai.

      "Mi piace"

      Rispondi

  5. Trovato l’inghippo e scoperto l’arcano! Grazie per l’assistenza e grazie anche per aver scritto questo articolo!

    Stavo provando la tua espressione usando il computer di un amico con installato QGIS 3.10.3….

    Per fortuna mi sono ricordato che il comportamento della funzione QgsGeometry::interpolate(), su cui si basa line_interpolate_point(), è stato cambiato a partire da QGIS 3.4 e poi di nuovo a partire da QGIS 3.10.4: nelle versioni tra 3.4 e 3.10.3 poteva interpolare solo nella prima parte delle linee multipart; da 3.10.4 in poi l’interpolazione è stata estesa a tutte le parti (come era prima della versione 3.4).

    Quando il percorso interseca più volte lo stesso Comune, si generano linee multipart e se la prima parte ha una lunghezza inferiore al 10% (in base al parametro stabilito) della lunghezza totale del percorso in quel dato Comune, QGIS 3.10.3 non calcola correttamente l’interpolazione. Per avere più chance di non ricadere in tale casistica e ottenere una corretta interpolazione, ho ridotto il fattore da 0.1 a 0.00001.

    Dai test fatti poi con QGIS 3.10.7 quello che riscontro è che purtroppo non sembra comunque possibile gestire correttamente proprio il caso in cui il percorso attraversa più volte lo stesso Comune.

    Per esempio per il Comune di Roasio e Villa del Bosco l’ordinamento presente nella colonna sort_D-A è il seguente:
    …
    Brusnengo 32
    Villa del Bosco 33
    Roasio 34
    Lozzolo 35
    …

    In realtà, come si può vedere seguendo la linea del percorso, dopo Brusnengo esso attraversa prima Roasio, poi Villa del Bosco, poi di nuovo Roasio, per poi passare da Lozzolo:
    …
    Brusnengo
    Roasio
    Villa del Bosco
    Roasio
    Lozzolo
    …

    Altri comuni con questo problema sono Samone e Salerano Canavese, Anzola d’Ossola e Premosello-Chiovenda.

    Purtroppo, penso inevitabilmente con il solo calcolatore di campi e senza creare un nuovo layer, non si può gestire correttamente questa problematica (perché non si possono aggiungere ulteriori records in questo modo).

    L’uso del fattore 0.1 o 0.9 allora serve a “dirimere la questione” nel caso un cui percorso attraversi il territorio di un Comune più di una volta:
    utilizzando il fattore 0.1 quel Comune verrà posto in elenco quando è attraversato dal percorso per la prima volta, mentre il passaggio per le volte successive verrà ignorato
    …
    Brusnengo 32
    Roasio 33
    Villa del Bosco 34
    Lozzolo 35
    …
    utilizzando il fattore 0.9 quel Comune verrà posto in elenco quando è attraversato per l’ultima volta, ignorando le precedenti
    …
    Brusnengo 32
    Villa del Bosco 33
    Roasio 34
    Lozzolo 35
    …

    Spero che queste note un po’ prolisse possano essere d’aiuto a chi si è trovato con lo stesso mio problema e con gli stessi dubbi.

    Grazie ancora per il supporto.

    Piace a 1 persona

    Rispondi

    2. Ho provato a semplificare la tua espressione eliminando la necessità di utilizzare le funzioni “line_interpolate_point” e “length” e i due fattori arbitrari 0.1 e 0.9 (utilizzando al loro posto le funzioni “start_point e “end_point”) e di calcolare due volte la stessa intersezione per ogni “line_locate_point”.

      In sostanza ho sostituito le due ricorrenze di

      line_locate_point(@toto, line_interpolate_point( intersection($geometry,@toto ), length(intersection($geometry,@toto))*0.1 ))

      con

      line_locate_point(@toto, start_point(intersection($geometry, @toto)))

      Per cui l’espressione diventa:

      with_variable(‘toto’, reverse(geometry(get_feature(‘percorso’, ‘fid’, 1))),
      array_find(
      array_filter(
      array_sort(
      array_agg(line_locate_point(@toto, start_point(intersection($geometry, @toto))))
      ), @element>=0),
      line_locate_point(@toto, start_point(intersection($geometry, @toto)))
      )
      )

      (Nota anche “@element>=0” invece che “@element>0”)

      In questa maniera, oltre ad evitare il bug di “line_interpolate_point” per le versioni inferiori a 3.10.4, si riducono anche i tempi di esecuzione:
      con l’espressione originale il mio pc impiega 13 secondi, con l’espressione modificata 8 secondi.

      Ho verificato che è possibile ridurre ulteriormente il tempo d’esecuzione facendo in modo che l’espressione venga calcolata solo per i poligoni che effettivamente sono attraversati dal percorso, aggiungendo un’istruzione “if”:

      with_variable(‘toto’, reverse(geometry(get_feature(‘percorso’, ‘fid’, 1))),
      if(intersects($geometry, @toto),
      array_find(
      array_filter(
      array_sort(
      array_agg(
      line_locate_point(@toto, start_point(intersection($geometry, @toto))))
      ), @element>=0),
      line_locate_point(@toto, start_point(intersection($geometry, @toto)))
      ), NULL)
      )

      In questo modo il tempo d’esecuzione è di 5 secondi.

      Piace a 2 people

      Rispondi

  6. Banalmente, usando un cronometro 🙂
    Ma volendo si potrebbe usare il Field Calculator degli strumenti di processing tramite la Python Console e misurare il tempo di esecuzione con process_time() per esempio.
    Però il Field Calculator degli strumenti di processing è più lento perché scritto in Python ed è necessario generare sempre un nuovo file di output, quindi i tempi di esecuzione dovrebbero essere più lunghi.
    Comunque le differenze nei tempi di esecuzione di varie espressioni dovrebbero potersi apprezzare ugualmente.

    Piace a 1 persona

    Rispondi

Lascia un commento

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

Sono un ingegnere civile, nel 2008 seguii un corso di "Alta formazione per progettisti di Sistemi Informativi Territoriali" (SIT) dove iniziai a conoscere SO linux, SW GIS ( Autodesk, ESRI e FOSS) e mi 'innamorai' della 'Geomatica'.