Leggere il tema “ECOSTRESS” in particolare nel contesto di cambiamento climatico e con una mappa associata fa scattare subito nell’immaginario collettivo qualcosa che è in grado di calcolare lo “stress ecologico” del pianeta.

In realtà si tratta dell’azzeccato acronimo di Ecosystem Spaceborne Thermal Radiometer Experiment on Space Station ( = esperimento del radiometro termico spaziale per l’osservazione degli ecosistemi dalla stazione spaziale): un esperimento scientifico della NASA, partito nel 2018 e che si concluderà quest’anno, dove un radiometro avanzato montato sulla International Space Station (ISS) raccoglie misure dalle piante che rispondono a stress dovuto al riscaldamento globale.

I dati raccolti da ECOSTRESS servono a rispondere principalmente a queste domande:

1. Come sta rispondendo la biosfera terrestre ai cambiamenti nella disponibilità di acqua?
2. Come influiscono i cambiamenti nello stress idrico vegetale diurno sul ciclo globale del carbonio?
3. È possibile ridurre la vulnerabilità agricola mediante il monitoraggio avanzato del consumo di acqua agricola e un miglioramento della stima delle siccità?

E lo fa misurando accuratamente la temperatura delle piante. Le piante regolano la loro temperatura rilasciando acqua attraverso piccoli pori sulle loro foglie chiamati stomi. Se dispongono di acqua sufficiente, possono mantenere la loro temperatura, ma se l’acqua è insufficiente, le loro temperature aumentano e questo aumento di temperatura può essere misurato con ECOSTRESS. Le immagini acquisite da ECOSTRESS sono le immagini di temperatura più dettagliate mai acquisite dalla superficie dello spazio e possono essere utilizzate per misurare la temperatura di un singolo campo agricolo.

Uno dei principali prodotti della missione ECOSTRESS è l’Indice di Stress Evaporativo (ESI): un indicatore leader di siccità in grado di indicare quanto le piante sono sotto stress è pertanto con quale probabilità si può verificare una siccità.
Maggiori informazioni si trovano al sito ufficiale.

Mappe delle isole di calore

Nel 2022 l’Agenzia Spaziale Europea ha pubblicato una serie di immagini a partire da quei dati raccolti nel”articolo City heat extremes — all’interno del progetto — mostrando così le isole di calore di diverse città del mondo ed evidenziando ancora di più che la presenza di alberi è fondamentale.
Un esempio è la comparazione fra Praga e Milano dove le immagini confermano da sole l’importanza della presenza degli alberi.

La curiosità quindi degli appassionati di mappe diventa immediata nel momento in cui si domanda se ECOSTRESS raccoglie anche dati della zona in cui vive.

Scaricare i dati

Per scoprirlo è sufficiente navigare i dati che NASA espone sul suo portale Earthdata.
Cercando la parola “ecostress” e scegliendo il prodotto “ECOSTRESS Land Surface Temperature and Emissivity Daily L2 Global 70m V001″ ed infine disegnando un rettangolo sull’Italia (eventualmente inserendo anche il periodo storico che si vuole esaminare), l’interfaccia mostra tutte le zone che ECOSTRESS ha raccolto.

Per conoscere a fondo i dati è importante andare a leggere le informazioni sui dettagli che si raggiungono con un clic sulla barra verticale tratteggiata a sinistra del nome del file.

Questo apre una finestra con una importante serie di metadati che si consiglia di guardare con molta attenzione in particolare a risorse in PDF come questa guida ECO2_User_Guide_V1.pdf

Processare i dati per elaborali

questo anche perché la prima azione della persona desiderosa di vedere i dati è quella di scaricare e aprire con un software come QGIS

I file sono griglie di 70m (raster) in formato HDF5, pertanto ricchi di informazioni.
Infatti, l’apertura di uno di questi in QGIS, mostra diverse sorgenti.

Il file con le informazioni delle temperature di superficie è LST — Land Surface Temperature.

La sola visualizzazione anche come scale di grigi però non aiuta a comprendere dove ci si trova e sembra molto lontana da quella che è l’immagine di anteprima che il portale offre.

Inoltre non è presente nemmeno l’informazione sul sistema di riferimento e quindi questo rende il tutto ancora più difficile.

Questo perché il dato è grezzo ed ha bisogno di essere elaborato al fine di ottenere tutti i prodotti in formato geoTIFF.
La serie di operazioni da svolgere è contenuta nello script python ECOSTRESS_swath2grid.py

Leggendo le istruzioni nel repository che contiene il sorgente si scopre che è necessario avere anche i dati del prodotto ECO1BGEO che contiene il layer di dati per valori di latitudine e longitudine, informazioni sulla geometria solare e della vista, altezza della superficie e frazione di pixel sulla terra rispetto all’acqua necessari quindi per aver un raster completo e georiferito.

Questo prodotto si può scaricare sempre dal sito EarthData di NASA ma diventa estremamente più facile scaricarlo direttamente dal sito web che esporta i file verso l’esterno raggiungibile a questo indirizzo.

L’interfaccia è quella tipica delle directory esposte via web.

per raggiungere il file giusto si tratta quindi di dare una occhiata al nome del file .h5 scaricato precedentemente.
Nel mio caso il file si chiama ECOSTRESS_L2_LSTE_28491_008_20230716T002924_0601_02.h5

Il primo passaggio è quindi individuare la data che è 20230716 quindi 2023.07.16 e da qui cercare il file che contiene la stringa _28491_008_20230716T002924_0601 con estensione .h5 e scaricarlo nella stessa directory del file precedente

Ora l’operazione è più facile di quello che sembra.
Una volta scaricato lo script python e verificato di avere installato tutte le librerie dipendenti (h5py pyproj pyresample gdal pandas scipy) si tratta solo di eseguirlo specificando la directory dove si trovano i file da elaborare e se si vuole avere i dati proiettati usando unità in metri (UTM) oppure in gradi (GEO).
Ad esempio se i file sono nella directory download e si vorrà ottenere una proiezione UTM allora il comando sarà:

python ECOSTRESS_swath2grid.py --proj UTM --dir download/

Il computer comincerà a calcolare e, a fine calcolo, nella sottodirectory “output” si troveranno tutti i geoTIFF convertiti.

Il file che contiene il suffisso “LST_UTM.tif” è quello con le temperature di superficie che si presenterà — nel mio caso — con la proiezione WGS84 — UTM32N (epsg 32632)
Infatti, sovrapponendolo ad una mappa di OpenStreetMap compare sopra l’Italia

la cui colorazione poi può essere fatta in con una scala dal blu al rosso passando per il verde e il giallo.

ed evidenziare così quali sono le zone dove le temperature sono più o meno alte in una cella di 70m, giocando poi con zoom e trasparenze individuare le aree di una città.

Qui però arriva il momento più importante: studiare meglio i dati dalla documentazione fornita da ECOSTRESS verificando anche date e orari del giorno e proseguire con le proprie analisi e considerazioni.

Il consiglio è DIVERTIRSI e STUDIARE!

Esplorare ECOSTRESS: scaricare e organizzare i dati termici della Terra offerti da NASA was originally published in tantotanto on Medium, where people are continuing the conversation by highlighting and responding to this story.

A person I respect very much, has written to me an email in which he asks to me how to create an animated gif of cartographic data, using the simple and convenient ways typical of the command line.

Several ways to do it, here is the MapServer one.

The goal

What I did was build an animation GIF with my collecting sites in Tasmania. The data are from here: https://doi.org/10.3897/BDJ.9.e70463 but I extracted the Tasmanian collecting events and reduced them to a 3-column table — year, lat, lon — with unique entries (sort | uniq).

The data are points, and in the desired output the red points are those of the specified year and the blue ones, those of the previous years. This is the GIF created by my friend, but his problem was how slow he was doing it:

The slow part was making each of the year categories visible, one after the other, changing the marker stylings and exporting as a PNG with year as filename — 45 minutes.

The MapServer way

MapServer is an Open Source platform for publishing spatial data and interactive mapping applications to the web. Originally developed in the mid-1990’s at the University of Minnesota, MapServer is released under an MIT-style license, and runs on all major platforms (Windows, Linux, Mac OS X).

It has a cli — mapserv - that you can use to create outputs. Moreover it can use run-time substitution to change variables values at run-time.

Then I can use it to create a PNG output file for every year, and then to process these and create an animated GIF.

To use MapServer it’s necessary to create a map file, a configuration file in which set the geographical bounding box, the input layers you want to use, the graphic styles to apply, etc..
I will not explain here what the syntax is, it is not the purpose of this post, I'll just make a few comments.

Some points for the example map file I have created:

• I have set a decimal degrees bounding box (EPSG:4326);
• the coordinate reference system (CRS) of the output is again “latitude, longitude”, but it is possible to set the CRS you want;
• I have used 3 layers: 1) one to have a background map, the polygons of Natural Earh Admin 0 – Countries layer, low resolution, used to create an example; 2) the data layer, a 3 fields CSV (decimalLatitude, decimalLongitude and year); 3) the layer used to add the year as a label;
• I have set 2 styles for data layer, the red dots for the specified year and the blue dots the previous years. The definition is by variable, passed at run-time via cli (EXPRESSION ("[year]" < '%year%')).

The data layer, is defined as GDAL/OGR virtual layer, that is one of MapServer input formats. It’s based on a XML definition file, to map CSV coordinates fields (and some other metadata).

Then to create a PNG, just run something like:

mapserv -nh "QUERY_STRING=map=data.map&mode=map&year=1977" >1977.png

data.map is the Mapserver configuratione file and year is a variable that you can change at run-time. Then in bash it could be:

for i in {1975..2000}; do
  mapserv -nh "QUERY_STRING=map=data.map&mode=map&year=${i}" >${i}.png
done

To add for each year, the year as label, I have added this in the query (labelyear is the name that I use to add labels to the PNG output):

&map_layer[labelyear]=FEATURE+POINTS+300+-560+END+TEXT+"1977"+END

It’s possible to do it, because in MapServer you can change map file parameters via URL.

The bash script

I have created a bash script to run all:

• first of all it downloads the data layer and extract decimalLatitude, decimalLongitude and year;
• then it creates a PNG for every year, using MapServer cli;
• then it creates the below animated GIF using convert ImageMagick utility.

The tools I use in it are MapServer, Miller and ImageMagick.

The output

Create an animated gif of cartographic data via command line was originally published in tantotanto on Medium, where people are continuing the conversation by highlighting and responding to this story.

I cittadini di Ottawa sembrano avere una lunga tradizione con l’uso della bicicletta e non stupisce quindi che dal 1984 esista un’associazione senza fini di lucro dal nome “Bike Ottawa” che crede fortemente che muoversi in bicicletta ad Ottawa dovrebbe essere una scelta facile e sicura per chiunque e pertanto si impegnano a fare in modo che le strade della loro città siano sicure.

Negli anni il collettivo di Bike Ottawa così è cresciuto in maniera importante mantenendo sempre l’attenzione alta e costruendo strumenti utili alla discussione sul tema.

Sul sito ufficiale si trovano quindi tantissimi articoli, raccolta di documentazione sul tema della ciclabilità, report, analisi ed altro ancora che raccontano e stimolano a ragionare sul valore dell’uso della bicicletta.

Fra gli strumenti anche mappe con cui pianificare percorsi sicuri, raccogliere informazioni su punti di interesse, immaginare percorsi per il futuro, mostrare una mappa di come la ciclabilità cambia d’inverno, dove sono avvenuti incidenti negli anni ed altro ancora.

In questo “altro ancora” una mappa molto affascinante è quella del livello da stress che si subisce percorrendo una singola strada in bicicletta.

Uno stress che è calcolato sulla base di come si presenta il percorso ciclabile in relazione alla strada che lo contiene e nel rapporto con le strade che incrocia. Ad esempio si tratta di un percorso condiviso con auto o con pedoni o è vicino a strade alto scorrimento ecc…

Un modello che il data working group di Bike Ottawa ha costruito sulla base della conoscenza maturata in questi anni, fatto da una serie di regole basate sui dati con cui sono descritti ogni singoli tratti di strade.

L’aspetto più affascinante di tutto questo è poi la riproducibilità con cui il modello può essere applicato.

Questo perché, Bike Ottawa non solo rilascia tutto il codice che produce su GitHub, ma utilizza anche i dati di OpenStreetMap.

Questo fa sì che, dando una sbirciatina al codice e leggendo le poche e semplici istruzioni, chiunque è in grado di applicare il modello del livello di stress di ciclabilità su una qualsiasi parte del mondo.

E così, curiosando sul singolo repository dello stressmodel, ecco che noi due (Maurizio Napolitano e Matteo Fortini) ci siamo messi a sviluppare le mappe da “stress in bicicletta” delle città dove viviamo (Trento e Cento).

Da lì in poi, in brevissimo tempo, è nata l’idea di fare una mappa per tutta Italia dando vita a “Bici stressate dal traffico”.

Appena si visita la mappa il primo pensiero è quello di fare zoom sulla propria città per farsi subito un’idea se, quanto elaborato corrisponde a realtà e se ci sono caratteristiche da scoprire. E così ci si da da subito risposta a domande come: “la ciclabilità nella mia città è rilassante oppure bisogna sempre stare attenti?” oppure “Il percorso da casa al lavoro ha delle strade che non mi mettono in difficoltà?”.

C’è chi rimane soddisfatto, c’è invece che pensa che qualcosa proprio non torna.

In effetti, se si guarda la mappa di tutta Italia, ci si incuriosisce subito nel vedere che, sulla montagna, sembrano esserci percorsi a basso livello di stress se non addirittura percorribili da bambini.

Le motivazioni sono diverse, e il fatto che ci sia un errore e che questo venga evidenziato, a noi piace molto, perché scopo della mappa è anche quello di stimolare il dibattito, replicare, modificare il modello (è costruito su Ottawa e quindi qualcosa potrebbe essere diverso — le informazioni si trovano nel lavoro di Peter Furth delle Northeastern University)

Sicuramente uno dei problemi principali è anche la mancanza di dettaglio nei dati (o la totale mancanza).

La comunità di Bike Ottawa è molto attenta al dettaglio ed ha creato anche una guida su come mappare i dati della ciclabilità dentro OpenStreetMap.

Forse anche questa andrebbe rivista su un modello italiano, rimane però il fatto che la quantità di dettagli necessaria fa capire da subito quanto le persone di Bike Ottawa abbiano amore per il loro territorio e quanto credano fortemente in questo progetto.

Un suggerimento che ci piace dare a chi ci sta leggendo e a chi vorrebbe partecipare per migliorare il progetto è quello di contribuire a OpenStreetMap, magari cominciando a “giocare” con l’app per android “StreetComplete”.

Una volta installata l’applicazione ci si potrà divertire a gironzolare per la propria città raccogliendo dati di dettaglio su percorsi pedonali, ciclabili, altezze degli edifici, conferma della presenza di panchine, cestini, fontanelle dell’acqua, descrizione del manto stradale, numero di scalini … ecc.. ecc…

Piccoli dati di dettaglio che sono poi in grado di fare la differenza e di aiutare non solo ciclisti, ma chiunque abbia la necessità di raggiungere un luogo con mezzi diversi.

Il nostro desiderio poi è quello di aggiornare di volta in volta questa mappa che abbiamo creato, sperando che il dibattito sulla necessità di creare spazi per una mobilità sostenibile parta anche da considerazioni oggettive che guardano anche ai dati.

Grazie!
Maurizio e Matteo

Lo stress dei percorsi ciclabili was originally published in tantotanto on Medium, where people are continuing the conversation by highlighting and responding to this story.

Con Maurizio e Salvatore abbiamo raccontato Il decreto “di Natale”, in chilometri, ovvero come calcolare le aree in cui sarà possibile spostarsi nei giorni 28, 29, 30 dicembre 2020 e 4 gennaio 2021, secondo quanto indicato nel Decreto Legge numero 172 del 18 dicembre 2020.

Maurizio le ha calcolate usando Python, Salvatore con QGIS e SpatiaLite; io l’ho fatto “a riga di comando”, sfruttando alcune utility e uno script bash.

Utility usate

Lo script sfrutta principalmente queste 4 utility:

• GDAL/OGR, la più importante libreria open source per leggere e scrivere file geografici vettoriali e raster https://gdal.org/
• Mapshaper (usare versione >= 0.5.22), una straordinaria applicazione open source per modificare file in formato Shapefile, GeoJSON, TopoJSON, CSV, ecc. https://github.com/mbloch/mapshaper
• Miller, una straordinaria applicazione open source per elaborare file in formato CSV, TSV, ecc. https://miller.readthedocs.io/en/latest/features.html
• SpatiaLite, uno straordinario RDBMS open source con funzioni spaziali https://www.gaia-gis.it/fossil/libspatialite/index.

Per usare lo script, è necessario installarle.

Lo script

È diviso in tre parti:

• l’elaborazione dei dati sulla popolazione;
• l’elaborazione dei dati geografici;
• l’elaborazione per produrre i dati a supporto della mappa online che abbiamo realizzato.

Elaborazione dati sulla popolazione

Abbiamo scelto come fonte questa di ISTAT: http://demo.istat.it/pop2020/dati/comuni.zip.

Sono dati aggiornati a gennaio 2020 con queste caratteristiche:

• formato CSV, encoding UTF-8 e la , come separatore di campi;
• 807.233 righe;
• 2 righe di intestazione, una descrittiva e una con i nomi dei campi;
• 2 righe al piede, con note, a partire da riga 807.232;
• 19 campi. Tra questi: il codice ISTAT in formato numerico del comune, il nome del comune, la classe di età della popolazione (da 0 a 100 a passo di un anno, con uno speciale valore di 999 per il totale di popolazione per comune), il totale di popolazione dei maschi, il totale di popolazione delle femmine;
• per ogni comune un record per ogni classe di età.

Un estratto:

"Popolazione residente al 1° Gennaio 2020 per sesso età (b) e stato civile post censimento (n)"
Codice comune,Denominazione,Età,Maschi celibi,Maschi coniugati,Maschi divorziati,Maschi vedovi,Maschi uniti civilmente,Maschi già in unione civile (per scioglimento),Maschi già in unione civile (per decesso del partner),Totale Maschi,Femmine nubili,Femmine coniugate,Femmine divorziate,Femmine vedove,Femmine unite civilmente,Femmine già in unione civile (per scioglimento),Femmine già in unione civile (per decesso del partner),Totale Femmine
1001,Agliè,0,5,0,0,0,,,,5,11,0,0,0,,,,11
1001,Agliè,1,9,0,0,0,,,,9,15,0,0,0,,,,15
1001,Agliè,2,6,0,0,0,,,,6,8,0,0,0,,,,8

Le informazioni di sintesi sullo schema:

Per poterlo utilizzare è necessario rimuovere la prima riga di intestazione. Nello script è stata rimossa con tail, con l’opzione -n e l’argomento +2, per avere in output il file a partire dalla seconda riga:

tail <./comuni.csv -n +2

Successivamente vengono rimosse le righe al piede, che contengono le note. Viene fatto in due passaggi: prima estraendo la riga che contiene i nomi dei campi, e poi aggiungendo a questa tutte le righe che iniziano con un numero (quindi tutto, tranne le righe con le note, perché il primo campo contiene i codici ISTAT dei comuni).
Gli strumenti sono i due classici head, per estrarre soltanto la prima riga, e grep per filtrare tramite espressioni regolari le righe che iniziano con un numero:

# estri riga con nomi campo
head <./comuni.csv -n 1 >./tmp_comuni.csv
# aggiungi corpo, rimuovendo ciò che non inizia per codice comune (in modo da rimuovere footer)
grep <./comuni.csv -P '^[0-9]+' >>./tmp_comuni.csv
# rinomina file
mv ./tmp_comuni.csv ./comuni.csv

Il passo successivo è quello di estrarre il totale di popolazione per comune. Per come è fatto il file, è necessario estrarre tutte le righe in cuiEtà=="999" e fare poi la somma del totale di maschi e femmine.
A questo è stata aggiunta l’estrazione dei campi “Codice comune” e “Abitanti” (e la rimozione dei restanti) e il cambio del loro nome.
È stato usato Miller:

mlr -I --csv clean-whitespace \
  then filter -S '${Età}=="999"' \
  then put '$Abitanti=${Totale Maschi}+${Totale Femmine}' \
  then cut -f "Codice comune",Abitanti \
  then rename "Codice comune",PRO_COM_T ./comuni.csv

La sintassi è molto leggibile. Qualche nota:

• ai nomi di campo si fa riferimento con $ seguito da nome campo. Qui sono state aggiunte le parentesi graffe, perché siamo in presenza di nomi “speciali”, con spazi e accentate;
• viene applicato il comando clean-whitespace (in Miller si chiamano “verbi”), per rimuovere eventuali spazi bianchi “in più”, ovvero 1 o più a inizio e fine cella, e più di 1 tra caratteri in un cella;
• l’opzione -I fa in modo che il comando lavori in sovrascrittura sul file.

In ultimo si è scelto di “standardizzare” il codice ISTAT dei comuni, da campo numerico a campo testuale a 6 caratteri: ad esempio trasformare il codice 1001 del comune di Agliè in “001001”.
Con Miller sfruttano il “verbo” put e la funzione fmtnum:

mlr -I --csv put '$PRO_COM_T=fmtnum($PRO_COM_T,"%06d")' ./comuni.csv

Elaborazione dati geografici

Ancora una volta la fonte è ISTAT, in particolare la versione generalizzata dei “confini delle unità amministrative a fini statistici”:
https://www.istat.it/storage/cartografia/confini_amministrativi/generalizzati/Limiti01012020_g.zip.

Sono dati aggiornati a gennaio 2020 con queste caratteristiche:

• formato Shapefile
• encoding UTF-8;
• sistema di coordinate EPSG:32632;
• tipo geometrico Polygon;

Queste le informazioni di sintesi sullo schema:

Come scritto nell’articolo che presenta il progetto, dato un comune con non più di 5.000 abitanti, e il suo confine, per calcolare l’area in cui da questo è possibile spostarsi (nei giorni 28, 29, 30 dicembre 2020 e 4 gennaio 2021), è necessario per ogni comune:

• calcolare l’area di buffer attorno al confine, di 30.000 metri;
• rimuovere da questa l’eventuale area dei comuni capoluogo che ricadono all’interno;
• rimuovere la parte che va al di fuori dei confini nazionali.

Il buffer, il luogo dei punti distanti 30.000 metri dal confine comunale, è stato calcolato tramite ogr2ogr:

ogr2ogr -t_srs EPSG:4326 ./comuni_30cappa_5mila.shp ./comuni.shp \
-dialect sqlite \
-sql "SELECT PRO_COM_T,COD_REG,COMUNE,Abitanti,st_buffer(comuni.geometry,30000) AS geom FROM comuni where Abitanti <= 5000"

Alcune note:

• come sistema di coordinate di output è stato scelto EPSG:4326, perché compatibile in modo nativo con le librerie per la pubblicazione di mappe sul web;
• il buffer viene calcolato tramite una interrogazione SQL di tipo spaziale, sfruttando il dialetto “sqlite” e la funzione ST_BUFFER. Questa vuole come argomenti la colonna geometrica e la distanza (di default usando l’unità di misura nativa, che qui sono metri);
• vengono estratte soltanto alcune colonne (codice comunale a 6 caratteri, codice regionale numerico e numero di abitanti).

I poligoni dei capoluoghi sono tutti quelli che nel file di input hanno CC_UTS==1. Vengono estratti con Mapshaper:

mapshaper ./comuni.shp \
-filter 'CC_UTS==1' \
-filter-fields PRO_COM_T,COD_REG,COMUNE,Abitanti \
-proj wgs84 \
-o ./capoluoghi_4326.shp

La sintassi è ancora una volta molto leggibile. Vengono estratti soltanto alcuni campi (come nella creazione dei buffer) e per le ragioni precedenti viene scelto EPSG:4326 come sistema di coordinare di output.

Per creare il limite poligonale dello stato italiano, basta unire — fare il dissolve — dei poligoni dei comuni del file di input.
Con Mapshaper:

mapshaper ./comuni.shp -dissolve \
-proj wgs84 \
-o precision=0.000001 ./italia.shp

Il comando dissolve è il cuore del processo e ancora una volta in output EPSG:4326.
Questo file viene usato per ritagliare le aree di buffer dei comuni. Si utilizza il comando clip di Mapshaper:

mapshaper ./comuni_30cappa_5mila.shp -clip ./italia.shp -o ./output.shp

Elaborazione dati per il sito web

La pagina web per presentare le aree in cui — per ogni comune — è possibile spostarsi ha una struttura di URL di questo tipo: https://ondata.github.io/30cappa/mappa.html?id=067042.

Al cambio di id, che qui è il codice ISTAT del comune con la rappresentazione a 6 caratteri , viene aperta una mappa centrata sul buffer di quel comune, con evidenziati l’area in cui è possibile spostarsi e il limite comunale.

Il poligono del comune viene restituita in risposta a un’interrogazione alle API SQL di CARTO, mentre quello del buffer è in risposta a una chiamata diretta a un file GeoJSON statico.
È quindi necessario generare circa 5.500 file, perché tanti sono i comuni con abitanti <=5.000.

Il comando essenziale “geografico” mancante è quello che sottrae all’aerea in cui ci si può spostare, quella occupata dai comuni capoluogo. Nell’immagine di sotto, ad esempio, quella dei capoluoghi di Teramo e Ascoli Piceno.

Per rimuovere queste aree, è stato utilizzato SpatiaLite e una sua funzione specializzata, ST_CUTTER (attiva dalla versione 4.4 e finanziata dalla Regione Toscana). Questo per la sua semplicità di utilizzo e rapidità di esecuzione.

Si inizia allora dall’importare i file Shapefile in un nuovo database in formato SpatiaLite. La struttura del comando è:

ogr2ogr -f SQLite -dsco SPATIALITE=YES ./db.sqlite ./input.shp -nln nomeTabellaOutput -lco SPATIAL_INDEX=YES -nlt PROMOTE_TO_MULTI

Il layer dei buffer ha una certa complessità. È opportuno in questi casi imporre una validazione geometrica e conseguente correzione (una geometria non validata potrebbe produrre errori e bloccare il processo):

ogrinfo ./db.sqlite -sql 'UPDATE nomeTabella SET GEOMETRY = MakeValid(GEOMETRY) WHERE ST_IsValid(GEOMETRY) <> 1;'

Si tratta di una query SQL di aggiornamento, in cui viene applicata la funzione “MakeValid” alla colonna geometrica (qui l’elenco delle funzioni). L’utility ogrinfo viene usata come “intermediaria”: non è necessario infatti usare un client nativo SQLite o SpatiaLite.

Il passo successivo è quello di “bucare” i poligoni dei buffer con i poligoni dei limiti amministrativi dei capoluogo. Viene fatto appunto a partire dalla funzione ST_CUTTER di SpatiaLite.
Ha bisogno di due oggetti: quello di input e quello che farà da “lama” di taglio, ovvero i nostri poligoni di buffer e quelli dei limiti dei capoluoghi.

A partire ad esempio dai tre poligoni di buffer di sotto in grigio, il poligono del capoluogo in bianco sarà usato per ritagliare questi tre.

Il risultato è quello che si vede nell’immagine seguente (in cui i poligoni sono stati volutamente spostati e ridimensionati):

• ognuno dei poligoni di input viene “tagliato” dal poligono che fa da “lama”;
• ognuno dei poligoni di input viene, in questo caso, diviso in due parti.

Dati due layer di un database SpatiaLite, la query per realizzare quanto descritto è (per eseguire il “taglio”):

SELECT ST_Cutter('input-db', 'input-layer', 'input-geometry-field', 'lama-db', 'lama-layer', 'lama-geometry-field', 'output-layer', 1, 1);

Gli argomenti principali della funzione sono il nome del layer di output, quello del db che contiene il layer di input, il nome del layer di input e il nome del campo geometrico del layer di input. Lo stesso per la “lama”.
Se si importano nel db SpatiaLite i due layer e gli viene assegnato rispettivamente il nome di “input” e “lama”, la query sarà:

SELECT ST_Cutter(NULL, 'input', NULL, NULL, 'lama', NULL, 'output', 1, 1);

NOTA BENE: se il database corrente è quello che contiene sia il layer di input che la “lama”, il nome del db si può impostare a NULL. Se il campo geometrico dei due layer è uno, non è necessario specificarne il nome e si può impostare a NULL.

Il layer di output avrà associata una tabella come quella sottostante. Questa è generata a partire dalle geometrie delle immagini soprastanti.

Alcune note:

• “PK_UID” è la colonna con gli ID (gli identificativi numerici distinti) dei 6 poligoni di output;
• “input_input_pk_uid” è la colonna con gli ID dei 3 poligoni di input. È evidente che ognuno viene diviso in due parti;
• “blade_lama_pk_uid” è la colonna con gli ID dei poligoni che fanno da lama. Qui è uno solo. Contiene valori nulli per tutti gli oggetti che ricadono fuori la “lama”.

Quest’ultimo è un punto chiave: basterà rimuovere i poligoni con valori NULL per il campo “blade_lama_pk_uid” e si otterrà il layer dei poligoni di buffer “bucato”, in corrispondenza dei poligoni dei comuni capoluogo (che hanno fatto da “lama”).

Nello script il passo successivo è proprio questo. Da notare che i layer che fanno da input e da “lama” sono stati denominati “a” e “b”, mentre quello di output è denominato “out”.

ogrinfo ./db.sqlite -sql 'create table buffer_clipped AS
SELECT
  PRO_COM_T,
  COMUNE,
  geometry
FROM
(SELECT
  "input_a_ogc_fid",
  "blade_b_ogc_fid",
  a.pro_com_t PRO_COM_T,
  a.comune COMUNE,
  out."geometry"
FROM
  "out"
  LEFT JOIN a ON out.input_a_ogc_fid = a.ogc_fid
WHERE
  "blade_b_ogc_fid" IS NULL);'

Una volta creato e esportato questo layer, vengono uniti nuovamente quei poligoni appartenenti allo stesso comune, eventualmente separati in precedenza.

mapshaper ./input.shp -dissolve PRO_COM_T copy-fields=COMUNE -o precision=0.000001 ./output.shp

Viene usato ancora una volta il comando dissolve di Mapshaper, in cui è possibile impostare quale campo usare come criterio di “unione” (qui è PRO_COM_T) e quale/i campi copiare in output (qui è COMUNE).

Infine, solo ad uso delle mappa online che è stata prodotta, è stato generato un file GeoJSON per ognuno dei poligoni di buffer creato.

mapshaper ./buffer.shp -split PRO_COM_T -o format=geojson ./output_folder/

Qui il comando è split, a cui si passa come argomento il nome di campo in cui cercare valori distinti, a partire dai quali creare file con quel nome. In output circa 5.500 file.

Note conclusive

Lo script bash non è perfettamente coincidente con quanto descritto in questo articolo.
Qui alcuni comandi sono stati leggermente semplificati a vantaggio di una maggiore leggibilità.

Si sarebbe potuta usare anche una sola di queste utility per fare tutto (SQLite/SpatiaLite la candidata principe). Ma alcune sono più easy e rapide in certe operazioni e questa era anche una buona occasione didattica per “toccare” un po’ 4 straordinari esempi di applicazione open source a riga di comando.
Anche Mapshaper sarebbe stato adatto per processare la gran parte del flusso, ma sembra esserci un problema nella procedura di ritaglio di poligoni sovrapposti e con il comando erase.

L’occasione è stata ottima per avere confermato ancora una volta come il SQL geografico e SpatiaLite siano degli strumenti di grande efficienza e comodità (grazie a Sandro Furieri e a tutti gli sviluppatori di SpatiaLite).
Mi ha colpito molto scoprire come ci siano pochissime tracce web della funzione ST_CUTTER. Peccato, è un gioiellino.

Infine un grosso grazie a Maurizio e Salvatore, per la disponibilità costante al confronto, per le note, i rilanci e per essere come sono 🤣

L’elaborazione dei dati per la mappa del decreto “di Natale”, in chilometri was originally published in tantotanto on Medium, where people are continuing the conversation by highlighting and responding to this story.

Il decreto legge numero 172 del 18 dicembre 2020 — detto in maniera un po’ fastidiosa “di Natale” — definisce delle misure urgenti per le festività natalizie e di inizio anno nuovo.

In questo si legge:

nei giorni 28, 29, 30 dicembre 2020 e 4 gennaio 2021 si applicano le misure di cui all’articolo 2 del medesimo decreto del Presidente del Consiglio dei ministri 3 dicembre 2020, ma sono altresì consentiti gli spostamenti dai comuni con popolazione non superiore a 5.000 abitanti e per una distanza non superiore a 30 chilometri dai relativi confini, con esclusione in ogni caso degli spostamenti verso i capoluoghi di provincia.

Questo articolo a cura di aborruso, napo e pigreco, scritto per restituire un quadro sul contesto numerico e geografico di questa parte del decreto.

Se vuoi subito vedere in quale aree ti potrai spostare, a partire dal tuo “piccolo” comune — senza leggere tutto il post — fai click qui 🗺.

I comuni con popolazione non superiore ai 5.000 abitanti

Sono tantissimi in Italia, quelli in rosso nell’immagine sottostante.

In numeri (al giorno 27 settembre 2020):

• sono 5508 comuni su un totale di 7904 (circa il 70%);
• questi coprono circa il 55% della superficie del nostro territorio;
• sono abitati da circa 10.000.000 di persone, circa il 16% della popolazione complessiva.

La fonte è ISTAT sia per i dati sulla popolazione, che per i dati spaziali/geografici.

Quanto sono “grandi” 30 chilometri

La risposta è solo apparentemente tautologica. Nelle distanze c’è molto di percettivo, di sentito, di legato al contesto.

30 chilometri sono poco meno della distanza in linea retta dall’aeroporto di Fiumicino a Villa Ada (a Roma). Letta così a molti probabilmente farà pensare a una misura “piccola”.
In automobile quanti km sarebbero e soprattutto in quanti minuti si percorrerebbero? E in linea retta (quindi ad esempio a piedi) che pendenze si dovrebbero affrontare, che barriere fisiche sarebbe necessario superare?

Il decreto fissa la cosa in modo semplice: “una distanza non superiore a 30 chilometri dai relativi confini comunali”, ovvero quello che nei sistemi informativi geografici è un buffer di 30.000 metri attorno al limite comunale del comune di cui si vuole calcolare.

Prendiamo come esempio il comune di Torano Nuovo, in provincia di Teramo, in Abruzzo, con codice ISTAT 067042. È abitato da 1552 persone (dati ISTAT gennaio 2020) e copre un’area di circa 1050 ettari (il piccolo poligono di sotto).

L’area di 30 chilometri attorno al confine comunale — il buffer — è quella in arancione della successiva immagine.

Nel decreto è riportato inoltre “con esclusione in ogni caso degli spostamenti verso i capoluoghi di provincia”.

E all’interno di quest’area di 30 km ricadono i comuni capoluogo di Ascoli Piceno e Teramo (sotto in verde), che saranno quindi da rimuovere dall’area.

Nel precedente DPCM del 3 dicembre a cui il “decreto di Natale” fa riferimento si legge:

Ai sensi dell’articolo 1, comma 2, del decreto-legge 2 dicembre 2020, n. 158, dal 21 dicembre 2020 al 6 gennaio 2021 è vietato, nell’ambito del territorio nazionale, ogni spostamento in entrata e in uscita tra i territori di diverse regioni o province autonome

Avevamo quindi pensato che fosse necessario verificare anche se un comune è limitrofo a dei confini regionali. Il comune di Torano Nuovo — in Abruzzo — è vicino al confine regionale tra Abruzzo e Marche e pensavamo fosse necessario rimuovere la parte dell’area che ricade nelle Marche.
Sembra che non sia così 🤔

Quindi, dato un comune con non più di 5.000 abitanti, e il suo confine, per calcolare l’area in cui da questo è possibile spostarsi (nei giorni 28, 29, 30 dicembre 2020 e 4 gennaio 2021):

• si calcola area di buffer attorno al confine, di 30.000 metri;
• si rimuove da questa l’eventuale area dei comuni capoluogo che ricadono all’interno;
• si rimuove la parte che va fuori dai confini nazionali.

Gli abitanti del comune di Torano Nuovo potranno quindi nei giorni 28, 29, 30 dicembre 2020 e 4 gennaio 2021 — e nel rispetto di quanto scritto nel decreto legge e nei DPCM in questo citati — spostarsi nell’area arancione sottostante.

Una pagina per sfogliare il catalogo dei comuni coinvolti e visualizzare la mappa

Abbiamo realizzato una pagina in cui è possibile cercare il comune di interesse (con abitanti non superiori a 5.000), inserendo il nome del comune o anche il suo codice ISTAT. Nell’immagine di sotto è evidenziato in 1.
Per visualizzare la mappa delle aree in cui sono consentiti gli spostamenti, si potrà fare click su “vedi”. Sotto è evidenziato in 2.

Come è stato fatto

Tutto questo, partendo dai dati di origine, può essere realizzato attraverso vari strumenti/processi. Il “nostro” trio, ha deciso di farlo ognuno con lo strumento preferito, arrivando ad ottenere sempre lo stesso risultato.

In ordine di pubblicazione

• napo, ha creato un python notebook, basato principalmente su geopandas pronto per essere usato in Google Colab e MyBinder.
• pigreco ha fatto tutto con SpatiaLite e QGIS;
• aborruso “a riga di comando” con uno script shell.

“Se sono un cittadino di uno di questi Comuni e mi muovo nell’area calcolata, sto rispettando la legge?”

La domanda è più che lecita ed è quella che ci siamo posti anche noi, che da buoni “impallinati sul tema geospaziale” sappiamo benissimo quanti tranelli ci sono dietro ai confini comunali e la percezione delle distanze quando si parla di quelle aeree (“nel raggio di ..”) o di quelle su un grafo stradale.

Non possiamo dire che i nostri calcoli hanno valore legale, ma non possiamo nemmeno dire di avere sbagliato visto che le indicazioni fornite “sono […] consentiti gli spostamenti dai comuni con popolazione non superiore a 5.000 abitanti e per una distanza non superiore a 30 chilometri dai relativi confini, con esclusione in ogni caso degli spostamenti verso i capoluoghi di provincia.”

Rimane comunque una distanza — che una volta calcolata e rappresentata su mappa — non pare indifferente. Sarebbe stato meglio definirla in relazione al grafo stradale?

I tranelli dei confini comunali

Indicare le distanze senza spiegare come si percorrono è un problema classico: un luogo può essere a 200 metri in linea d’aria rispetto ad un altro, ma se nel mezzo c’è un fiume o altro ostacolo allora per raggiungerlo bisognerà trovare il primo passaggio di collegamento che probabilmente farà percorrere più della distanza scelta.

Quando poi la scala aumenta, come nel nostro caso, allora, oltre a mettersi di mezzo confini naturali come specchi d’acqua, montagne, costa ed altro ancora ci si aggiungono i problemi delle enclave ed exclave.
Su questo l’Italia offre un sacco di situazioni curiose e spesso sottovalutate.
Uno fra i più popolari è il caso di Campione d’Italia i cui comuni confinanti sono tutti della Svizzera italiana in quanto si trova proprio al di fuori del confine italiano. Un comune con 2003 abitanti ma che, per “usufruire’ dei 30km di distanza, comunque devono passare il confine svizzero.

Altro caso interessante è quello del territorio dell’Alta Valmarecchia dove, nel 2009, 10 comuni sono passati dalle Marche all’Emilia-Romagna, creando così una enclave.
La lista potrebbe proseguire a lungo. Di certo, in virtù dei 30 km, molte di queste si “fondono” e danno il ‘vantaggio” di allungare gli spazi come il caso del Comune di Fardella che con i suoi 577 abitanti e le sue enclave, riesce già a confinare con nove comuni.

Ragionamenti analoghi possiamo farli poi per quei comuni circondati dall’acqua come le isole (ma non solo).

Ulteriori problemi possono venire poi dai dati: il calcolo è fatto usando gli ultimi dati ISTAT disponibili per la popolazione il cui aggiornamento è di gennaio 2020 e qualcosa potrebbe essere cambiato in particolare nei valori prossimi al 5000 (nota: il comune di Castelnuovo Scrivia ha 5001 abitanti), così come il fatto che, per questioni di performance, abbiamo preferito usare i confini amministrativi generalizzati (dove i confini sono leggermente meno precisi).

Nota sull’output prodotto

Il processamento dei dati che abbiamo descritto, può produrre dei casi come quello dell’immagine sottostante, relativa all’area di 30 chilometri attorno al comune di Nissoria, in Sicilia.
In corrispondenza del punto 1 è evidenziata una piccola aerea arancione, isolata. Ci si può recare lì?
Nelle FAQ ufficiali sul decreto si legge “esclusione in ogni caso degli spostamenti verso i capoluoghi di Provincia”, quindi ad Enna non si può andare.
La nostra interpretazione è che ci si possa spostare verso queste aree isolate attraversando aree di capoluogo, purché all’interno dei 30 km e al di fuori del limite amministrativo di un capoluogo.

Un bot telegram dedicato

Salvatore Larosa ha creato un bot Telegram che — a partire dal nome di uno dei comuni con non più di 5000 abitanti — restituisce l’elenco dei comuni limitrofi in cui è possibile spostarsi e il link alle nostre mappe.

È bello che, prima di andare online e in poche ore, si è creato già un servizio correlato al tema e che dialoga con il nostro (un grazie a Salvatore).

In conclusione …

Per noi questo è stato un momento di confronto, di divertimento, di imparare l’uno dall’altro e di condividere questa storia con chi ci sta leggendo.

Di certo non è nostra intenzione quello di invitare gli abitanti dei “piccoli Comuni’ a muoversi liberamente: crediamo sia fondamentale il buon senso.

Forse siamo un po’ troppo pignoli, ma pensiamo che quando le parole non sono sufficienti per spiegare i concetti, forse è opportuno fornire strumenti come le mappe per far capire cosa si vuole ottenere.

Buon Natale ;)

Il decreto “di Natale”, in chilometri was originally published in tantotanto on Medium, where people are continuing the conversation by highlighting and responding to this story.

E ha senso finalmente dedicargli un mini tutorial

Alla fine del 2017 le mappe del catasto italiano venivano messe online dall'Agenzia delle Entrate, purtroppo con una licenza che le rendeva pressoché inutilizzabili: nessun uso commerciale e nessuna opera derivata, ovvero la licenza CC BY-NC-ND 4.0.

Ai tempi scrissi una frase ironica, ma non campata in aria:

Posso essenzialmente stamparmi qualche cartina, e appenderla nel bagno di casa. Già in un bagno pubblico vado fuori licenza.

Ma da ieri (24 settembre 2020) non è più così e sono pubblicate con una licenza che ne consente un utilizzo pieno, con il solo obbligo di citazione della fonte (la CC BY 4.0).

È un’ottima notizia!

La vecchia licenza per me era un blocco, anche creativo. Che senso aveva ad esempio fare un tutorial su come estrarre — a partire da una lista di coppie di coordinate — i numeri delle particelle corrispondenti?
La nuova licenza, mi “titilla il velopendulo” (cit.), mi rende libero ed ecco una piccola guida d’esempio.

Come estrarre le particelle catastali a partire da un elenco di coppie di coordinate

I dati catastali sono pubblicati al momento soltanto in WMS, un servizio pensato per la "rappresentazione" del dato: di base produce immagini.
L'URL del WMS del catasto è https://wms.cartografia.agenziaentrate.gov.it/inspire/wms/ows01.php.
Tutte le caratteriste (sistemi di coordinate, elenco e nomi dei layer, ecc.) del servizio sono disponibili facendo una richiesta GetCapabilities.

Nel WMS può essere reso disponibile il metodo GetFeatureInfo, che consente di estrarre gli attributi di un layer in una determinata posizione: al click su un pixel dell'immagine - al click su un punto - vengono restituite le informazioni sull'oggetto "toccato".

Per questi dati il metodo è reso disponibile ed è un modo per portare a termine l’obiettivo di questo tutorial.

Il layer delle particelle è interrogabile, difatti nelle capabilities si può leggere queryable="1" (quando non lo è, è 0):

<Layer queryable="1" opaque="0" cascaded="0">
  <Name>CP.CadastralParcel</Name>
  <Title>Particelle</Title>
  <Abstract>Il layer contiene le particelle catastali.Costituisce distinta particella catastale ogni porzione continua di terreno che sia della medesima qualità e classe o abbia la stessa destinazione d’uso, che sia situata in un medesimo comune ed appartenga allo stesso possessore.Le particelle, i fabbricati, le strade e le acque, in questo servizio presentati su layer separati, coprono l’intero territorio nazionale, ad eccezione dei territori nei quali il catasto è gestito, per delega dello Stato, dalle Province Autonome di Trento e di Bolzano. Nel layer Particelle sono rappresentati anche i bordi dei fabbricati eventualmente presenti all’interno delle particelle.</Abstract>
  <LatLonBoundingBox minx="2" miny="33" maxx="19" maxy="48"/>
  <BoundingBox SRS="EPSG:6706" minx="2" miny="33" maxx="19" maxy="48"/>
  <MetadataURL type="ISO19115:2003">
    <Format>application/vnd.ogc.csw.capabilities.response_xml</Format>
    <OnlineResource xlink:type="simple" xlink:href="https://csw.cartografia.agenziaentrate.gov.it/age-inspire/srv/ita/csw?SERVICE=CSW&amp;VERSION=2.0.2&amp;REQUEST=GetRecordById&amp;ID=age:D_0000_ITALIA&amp;ElementSetName=full&amp;outputSchema=csw:IsoRecord"
      xmlns:xlink="http://www.w3.org/1999/xlink"/>
  </MetadataURL>
  <Style>
    <Name>default</Name>
    <Title>default</Title>
    <LegendURL width="89" height="20">
      <Format>image/png</Format>
      <OnlineResource xlink:type="simple" xlink:href="https://wms.cartografia.agenziaentrate.gov.it/inspire/wms/ows01.php?language=ita&amp;version=1.1.1&amp;service=WMS&amp;request=GetLegendGraphic&amp;layer=CP.CadastralParcel&amp;format=image/png&amp;STYLE=default"
        xmlns:xlink="http://www.w3.org/1999/xlink"/>
    </LegendURL>
  </Style>
  <ScaleHint min="0.00997805696859274" max="2.49451424214819"/>
</Layer>

Una richiesta GetFeatureInfo per questo servizio ha questo schema di URL (gli a capo sono per renderlo leggibile):

https://wms.cartografia.agenziaentrate.gov.it/inspire/wms/ows01.php?
REQUEST=GetFeatureInfo&
SERVICE=WMS&
SRS=EPSG:3045&
STYLES=&
VERSION=1.1&
FORMAT=image/png&
BBOX=355714.5,4222470.1,355716.5,4222472.1&
HEIGHT=9&
WIDTH=9&
LAYERS=CP.CadastralParcel&
QUERY_LAYERS=CP.CadastralParcel&
INFO_FORMAT=text/html&
X=5&
Y=5

Per la descrizione di dettaglio di ogni parametro, suggerisco di leggere questa pagina.
Qualche nota qui utile all'obiettivo del tutorial.

SRS è per impostare il sistema di coordinate, che è da scegliere tra quelli previsti dal servizio. Qui verrà utilizzato il 3045 (ETRS89 / UTM zone 33N (N-E)).

BBOX (bounding box) è per definire l'area di cui si vuole generare l'immagine su cui fare la query. Perché il metodo WMS, vale la pena ribadirlo, espone una rappresentazione in pixel dei dati, su cui può essere possibile fare un interrogazione per singolo pixel. Le coordinate vanno inserite in questo ordine: minimo x, minimo y, massimo x, massimo y.

HEIGHT e WIDTH per definire le dimensioni in pixel dell'immagine su cui fare la query.

LAYERS e QUERY_LAYERS per selezionare il layer da interrogare, che qui è CP.CadastralParcel.

INFO_FORMAT, per impostare il formato di output della risposta alla query.

E infine X e Y. Sono le coordinate in pixel, e lo 0 è in alto a sinistra. Se ho quindi un immagine in pixel di 9 x 9, il pixel centrale sarà quello di coordinate 5,5.

Se apro nel browser la query di sopra avrò le informazioni utili sulla particella che ricade sotto il pixel di mio interesse:

Con i parametri che ho inserito, ho generato da un’area di un metro di lato un’immagine di 9 x 9 pixel, su cui è stata eseguita l’interrogazione sul pixel di coordinate 5,5. Quindi sto lavorando a una risoluzione di circa 10 centimetri.

L’output è in HTML ed è qualcosa come:

<table class="wmstable" border="0">
<tbody>
<tr><th class="wmstable" colspan=7>Strato CP.CadastralParcel 'Particelle'</th></tr>
<tr><th class="wmstable" scope="col">InspireId localId</th><td>IT.AGE.PLA.G273_003400.1298</td></tr>
<tr><th class="wmstable" scope="col">InspireId_namespace</th><td>IT.AGE.PLA.</td></tr>
<tr><th class="wmstable" scope="col">Label</th><td>1298</td></tr>
<tr><th class="wmstable" scope="col">NationalCadastralReference</th><td>G273_003400.1298</td></tr>
</tbody>
</table>

Per estrarre quindi da una lista di coordinate, per ciascuna coppia, i dati della particella sottostante, basterà scrivere uno script (in qualsiasi linguaggio) che faccia questo:

• per ogni riga di un file di input legga la coppia di coordinate;
• per ogni coppia costruisca il bounding box (per il parametroBBOX) che lo "circonda", ad esempio sottraendo e sommando 1 metro;
• per ogni bounding box faccia l’interrogazione sul punto centrale;
• per ogni output HTML estragga l’informazione utile (che nel caso di sopra è IT.AGE.PLA.G273_003400.1298).

A partire da una lista così fatta (colonne separate da tabulazione, è un TSV):

355461.5	4222490.7
356257.6	4222742.7
356259.7	4222236.7

Si può usare ad esempio uno script bash, che sfrutti curl per fare le richieste HTTP e grep per estrarre dall'HTML di output i dati di interesse. Qualcosa come:

#!/bin/bash

rm ./outputCatasto.txt

while IFS=$'\t' read -r X Y; do
  minX=$(bc <<< "$X - 1")
  minY=$(bc <<< "$Y - 1")
  maxX=$(bc <<< "$X + 1")
  maxY=$(bc <<< "$Y + 1")
  curl -sKL "https://wms.cartografia.agenziaentrate.gov.it/inspire/wms/ows01.php?REQUEST=GetFeatureInfo&SERVICE=WMS&SRS=EPSG:3045&STYLES=&VERSION=1.1&FORMAT=image/png&BBOX=$minX,$minY,$maxX,$maxY&HEIGHT=9&WIDTH=9&LAYERS=CP.CadastralParcel&QUERY_LAYERS=CP.CadastralParcel&INFO_FORMAT=text/html&X=5&Y=5" | grep -oE 'IT.+[0-9]' >>./outputCatasto.txt
done <./input.tsv

Restituirebbe:

IT.AGE.PLA.G273_0033A0.673
IT.AGE.PLA.G273_003500.3077
IT.AGE.PLA.G273_003500.2419

In cui ci sono tutti i dati di interesse (come particella e comune, che qui è Palermo, che ha codice G273).

Note conclusive

Mi ripeto: il tipo di licenza scelto non definisce soltanto cosa sia legalmente consentito, ma ha un impatto molto più ampio. Non mi sembra banale ad esempio il fatto che la precedente mi avesse tolto qualsiasi voglia di fare un po’ di didattica su questo servizio esposto dall'Agenzia delle Entrate e sulle richieste GetFeatureInfo dei servizi WMS.

Questa nuova licenza apre scenari creativi e di business, prima inimmaginabili e irrealizzabili. Penso che ne vedremo delle belle, è un’ottima cosa.

Per informazioni di questa natura, bisognerebbe esporre i dati anche con il servizio WFS, che consente di accedere ai dati e non a una loro rappresentazione. Se l’Agenzia è in ascolto, come sembra esserlo, questa è l’occasione per richiedere di farlo.

L’Agenzia delle Entrate dovrebbe pure offrire un accesso programmatico anche al servizio inverso: a partire dai dati di foglio e particella, restituire le coordinate dell’area corrispondente. È un servizio che è già attivo qui, ma è usabile soltanto a schermo e non anche come servizio.

Infine un grazie a Salvatore Fiandaca, che tanto tempo fa aveva provato a stimolarmi a scrivere questo post, ma insieme concordammo come — data la vecchia licenza — fosse inutile.

Altri articoli recenti e correlati a tema

• Visualizzare le mappe del catasto terreni con Google Maps, di Paolo Mistrangelo;
• WMS catasto ADE: estrarre i dati con il field calc di QGIS, di Salvatore Fiandaca;
• WMS CATASTO ADE: vettorializzare le particelle, di Salvatore Fiandaca.

Post Scriptum 1 (c’è un altro servizio attivo di query per coppia di coordinate)

Ho scoperto — qualche giorno dopo la pubblicazione del post — che c’è un altro servizio attivo, che risponde in JSON, che consente di avere i dati catastali a partire da una coppia di coordinate geografiche.

L’URL ha questa struttura (gli a capo per renderlo leggibile):

https://wms.cartografia.agenziaentrate.gov.it/inspire/ajax/ajax.php?
op=getDatiOggetto&
lon=13.62607730254396&
lat=37.327193489720905

E si può provare facendo click qui. In output si avrà:

{
  "SIGLA_PROV": "AG",
  "COD_COMUNE": "A089",
  "DENOM": "AGRIGENTO (AG)",
  "SEZIONE": " ",
  "FOGLIO": "131",
  "ALLEGATO": "null",
  "SVILUPPO": "null",
  "NUM_PART": "660",
  "TIPOLOGIA": "PARTICELLA"
}

Post Scriptum 2 (come usare GDAL per scaricare i dati come GeoTIFF)

In questo vecchio post abbiamo già descritto come fare. A partire da un file XML (qui specifiche per GDAL) in cui è definita la sorgente, come questo di sotto (qui in copia)

<GDAL_WMS>
  <Service name="WMS">
    <Version>1.1.1</Version>
    <ServerUrl>https://wms.cartografia.agenziaentrate.gov.it/inspire/wms/ows01.php?language=ita&amp;SERVICE=WMS</ServerUrl>
    <Layers>CP.CadastralParcel</Layers>
    <SRS>EPSG:6706</SRS>
    <ImageFormat>image/png</ImageFormat>
    <Transparent>FALSE</Transparent>
    <BBoxOrder>xyXY</BBoxOrder>
  </Service>
  <DataWindow>
    <UpperLeftX>2</UpperLeftX>
    <UpperLeftY>48</UpperLeftY>
    <LowerRightX>19</LowerRightX>
    <LowerRightY>33</LowerRightY>
    <SizeX>1073741824</SizeX>
    <SizeY>947419256</SizeY>
  </DataWindow>
  <BandsCount>3</BandsCount>
  <BlockSizeX>1024</BlockSizeX>
  <BlockSizeY>1024</BlockSizeY>
  <OverviewCount>20</OverviewCount>
</GDAL_WMS>

si lancia il comando gdal_translate con una sintassi come quella sottostante:

gdal_translate -projwin 13.34513 38.12016 13.35073 38.11676 -outsize 2000 1214 particelle.xml out.tif

In output avrò un file GeoTIFF che apparirà come la PNG sottostostante

Alcune note sul comando:

• -projwin per impostare il bounding box dell’area di cui voglio l’output;
• -outsize per impostare le dimensioni in pixel dell’immagine di output (e quindi la risoluzione.

Nota bene: alcuni layer (come quello dell’esempio soprastante) sono renderizzati e quindi esportabili soltanto in un certo intervallo di scala (fare riferimento al GetCapabilities del WMS).

Note sulla codifica dell’output della query al WMS dell’Agenzia delle Entrate

Nell’esempio di sopra, la query dà in output il codice G273_003400.1298. Qual è il significato?

La struttura della prima parte, come riportato nella documentazione ufficiale a pagina 4 (grazie Stefano Campus per avermela suggerita), è

CCCCZFFFFAS

ovvero:

• “CCCC”, rappresenta il codice nazionale del comune (es.: H282 per il comune di Rieti);
• “Z”, rappresenta il codice della sezione censuaria (es. A, oppure B). Se la sezione è assente si utilizza il carattere ‘_’
• “FFFF”, rappresenta il numero del foglio, riempito eventualmente con caratteri ‘0’ a sinistra se il numero ha meno di 4 cifre (es. 0001 per il foglio numero 1). Se la mappa rappresenta un quadro d’unione dei bordi di più mappe allora FFFF rappresenta il numero identificativo della richiesta (modulo 10000);
• “A”, rappresenta il codice allegato. Assume il valore 0 se la mappa non è un allegato. Se la mappa rappresenta un quadro d’unione dei bordi di più mappe allora “A” ha il valore ‘Q’;
• “S”, rappresenta il codice dello sviluppo. Assume il valore 0 se la mappa non è uno sviluppo. Se la mappa rappresenta un quadro d’unione dei bordi di più mappe allora S ha il valore ‘U’.

Mentre l’ultima parte, dopo il . è il numero di particella.

Le mappe catastali dell’Agenzia delle Entrate diventano finalmente utilizzabili was originally published in tantotanto on Medium, where people are continuing the conversation by highlighting and responding to this story.

L’idea di questo post si è inserita nel filone delle esperienze positive -nonostante tutto- di questo tempo, come questo risultato: pareva legalmente e amministrativamente impraticabile svolgere le riunioni e le assemblee di organismi istituzionali per via telematica. Quanti ostacoli, solo pochi mesi fa. Poi, puff! Tutti superati in un mese.

Una base di dati di interesse nazionale -leggiamo dal Codice dell’Amministrazione Digitale, CAD (art. 60) — è: “L’insieme delle informazioni raccolte e gestite digitalmente dalle pubbliche amministrazioni, omogenee per tipologia e contenuto e la cui conoscenza è rilevante per lo svolgimento delle funzioni istituzionali delle altre pubbliche amministrazioni, anche solo per fini statistici, nel rispetto delle competenze e delle normative vigenti”.

Già il CAD ne indica alcune, come il Repertorio nazionale dei dati territoriali, RNDT, l’Anagrafe nazionale della popolazione residente, il Casellario giudiziario e il Registro delle Imprese. Altre sono state aggiunte a seguito di altri provvedimenti legislativi. Questo insieme di banche dati viene continuamente ampliato dall’AgID: una lista aggiornata è disponibile qui.

Con la transizione dalla società industriale a quella -della nostra epoca- informazionale, i dati devono stare al centro dell’azione pubblica, anche in quanto valore per l’economia. Molte delle basi di dati di interesse nazionale catalogati da AgID rientrano anche tra le categorie tematiche di serie di dati di elevato valore stabilite dalla Direttiva (UE) 2019/1024 del 20 giugno 2019, relativa all’apertura dei dati e al riutilizzo dell’informazione del settore pubblico — Direttiva “Public Sector Information”. Sono classificate di elevato valore, ad esempio, basi di dati geospaziali, i database statistici, quelli relativi alle imprese e alla proprietà delle imprese e, ancora, alla mobilità.

Quindi, la Pubblica Amministrazione riveste per l’intera società un ruolo di catalizzatrice rispetto al tema “dati”. Questi ormai devono essere intesi come “materia prima” essenziale per il funzionamento dell’economia, e occorre che il nostro Paese si doti di una infrastruttura ad essi dedicata, così come in passato si è munito di una rete di trasporto e di telecomunicazione.

Tra le basi di dati geografici di elevato valore è risaputo che ha un notevole potenziale di riuso la banca dati degli indirizzi o, più formalmente, degli stradari comunali.

Qual è lo stato di questa importante banca dati in Italia?

Con Decreto del Presidente del Consiglio dei Ministri del 12 maggio 2016 “Censimento della popolazione e archivio nazionale dei numeri civici e delle strade urbane” è stato istituito l’Archivio Nazionale dei Numeri Civici e delle strade urbane, ANNCSU. Lo stesso DPCM ne ha assegnato la realizzazione all’ISTAT e all’Agenzia delle Entrate, mentre il suo aggiornamento è di competenza dei Comuni che, per assolvere questo compito, si possono avvalere della Regione quale intermediario infrastrutturale tra il livello centrale e locale “previa sottoscrizione di specifici accordi di servizio tra Regione, ISTAT, Agenzia delle Entrate e Comuni” (art. 4).

Dal 2016 ad oggi sono trascorsi 4 anni, ma questo Archivio non è ancora disponibile (provare qui).

Ovviamente vi sono state, a più riprese, richieste di notizie su questo importante archivio digitale, come si può constatare, per esempio, consultando lo Spazio di discussione sui servizi pubblici digitali, tramite la ricerca dell’acronimo “ANNCSU”.

Ne scriviamo oggi perché dell’ANNCSU si occupa il recentissimo Piano d’investimenti per lo sviluppo sostenibile delle città elaborato dall'Alleanza Italiana per lo Sviluppo Sostenibile, ASviS, e presentato lo scorso 28 maggio. Il documento dell’Alleanza “propone un programma da oltre 200 miliardi in 10 anni su ambiente, mobilità, trasformazione digitale, sanità e lotta alla povertà, elaborato da un gruppo di esperti”, tra cui anche di Stati Generali dell’Innovazione.

Nel raccomandare che “occorre completare nel più breve tempo possibile tutte le piattaforme centrali e locali previste dal Piano triennale per l’informatica della pubblica amministrazione 2019–2021 al fine di renderle interoperabili per l’accesso e lo scambio dei dati” (pag. 18), ASviS ricorda esplicitamente, tra gli altri, l’ANNCSU.

Il Piano ASviS, per esigenze di sintesi, si limita ad enfatizzare che gli interventi sollecitati “possono finalmente consentire lo switch off dei servizi di sportello della PA in favore della modalità digitale
anche attraverso l’utilizzo dell’applicazione IO”. Ma l’Archivio nazionale degli Stradari, lo ribadiamo, ha un potenziale di riuso affatto strategico e amplificabile (Google insegna), tale da contribuire concretamente a quell’uscita dalla crisi “rimbalzando in avanti”, auspicata dagli estensori del Piano ASviS, motivo che ne ha ispirato l’elaborazione.

Ecco perché ne scriviamo: vorremmo che grazie a questo autorevole richiamo e con nuova consapevolezza acquisita a seguito delle lezioni apprese in questi ultimi mesi, si avverasse un altro… puff! La pubblicazione dell’ANNCSU:

§ superando un parere del Garante per la protezione dei dati personali previsto da un normatore, forse eccessivamente prudente (per uno stradario?!?!);

§ creando tra ISTAT, l’Agenzia delle Entrate, i Comuni e tutti i portatori d’interesse -ad iniziare da ASviS- un forte legame e un impegno coeso volto a completare rapidamente la georeferenziazione dei numeri civici, l’integrazione con i CAP (Codici di Avviamento Postale) e altre features adeguate per fare quel “balzo in avanti” cui ambiamo;

§ facendo sentire fieri tutti gli attori che contribuiranno a raggiungere questo obiettivo, ma sempre con umiltà, quindi avendo atteggiamenti inclusivi, di apertura verso nuove idee e di curiosità per le soluzioni sviluppate da altri Paesi leader su questa tematica.

L’appartenenza ad una categoria di dati di elevato valore, secondo la Direttiva “Public Sector Information”, attesta l’ANNCSU come Bene Comune, qualità che prospetta -aderendo al paradigma dell’Open Government- la sperimentazione di modelli di produzione e di governance dei dati fortemente votati alla collaborazione e alla condivisione.

Si possono considerare soluzioni adottate da altri Paesi orientate in questa cornice (ad esempio: GB, Francia, Belgio, Paesi Bassi, Danimarca ed Estonia). È molto apprezzata in ambito internazionale, ad esempio, la Base Adresse Nationale (BAN), progetto iniziato nel 2015 dall’Istituto Nazionale per l’Informazione Geografica e Forestale, IGN, il Gruppo “La Poste”, l’Associazione OpenStreetMap France, il Dipartimento Etalab e con il supporto dell’Amministratore Generale dei Dati, AGD[1]. Cogestita attualmente anche insieme all’Autorità per la Regolazione delle Comunicazioni elettroniche, Poste e della distribuzione della Stampa, ARCEP, all’Istituto di Statistica e Studi Economici, INSEE, e alla Direzione Generale delle Finanze Pubbliche, DGFIP.

BAN è il riferimento per tutti gli indirizzi nel territorio francese: contiene la posizione geografica di oltre 25 milioni di indirizzi. Completano questo progetto altri servizi come, ad esempio, un’API di geocodifica liberamente accessibile. Il volume di utilizzo di questi strumenti dimostra la loro pertinenza e utilità: nel 2017, l’API di geocodifica della BAN ha registrato un forte aumento del suo pubblico con oltre un miliardo di richieste e 11,8 milioni di visitatori unici nei primi dieci mesi dell’anno contro 443 milioni di richieste e 5,5 milioni di visitatori unici nel 2016.

Questo archivio è originale e unico non solo per il suo contenuto -che lo rende ad oggi il database più completo per quanto riguarda gli indirizzi in Francia- ma anche per la sua governance inedita, perseguita da una compagine che riunisce amministrazioni, società pubbliche e collaboratori della società civile.

Si prendano ispirazioni da questo progetto francese, come da progetti di altre nazioni. Il nostro Paese possiede tutte le conoscenze e tecnologie che servono; certamente non mancheranno idee; si adotti una governance inclusiva — a proposito, perché non approcciare anche un dialogo con i giganti della geolocalizzazione multinazionali?

Il successo dell’ANNCSU dopo il CONV-19 è assolutamente a portata di mano in tempi brevi.

[1] Riferimenti al progetto DAB sono tratti dal “Rapport au Premier Ministre sur la Donnée dans les Administrations”, 2016–2017 e dal Rapporto Etalab 2019.

Banche dati di interesse nazionale: quanto ci serviranno i numeri civici per la ripresa! was originally published in tantotanto on Medium, where people are continuing the conversation by highlighting and responding to this story.

Datawrapper è uno dei migliori servizi online per raccontare storie basate su grafici, mappe e tabelle. È cosa nota e non ha bisogno di presentazioni.

Però leggo spesso in chat private, in email, ecc. “peccato non si possano creare mappe a partire da fonti live, che si aggiornano automaticamente …”.

Su Datawrapper invece si può fare. Questi gli step:

• crea una nuova mappa, scegli la tipologia (nell'esempio di sotto è una “Choropleth map”) e quale layer di base vuoi usare (ho scelto “Italy » Regions”);
• aggiungi i dati e guarda l’URL della mappa. Sarà qualcosa come

https://app.datawrapper.de/map/CHART-ID/data.

• copia il CHART-ID della mappa e usalo per costruire questo URL con il tuo CHART-ID

https://app.datawrapper.de/chart/Chart-ID/upload

• lancia l’URL così costruito. Si aprirà un’interfaccia per fare l’upload dei dati e qui potrai impostare la fonte live. Da qui procedi a selezionare una fonte esterna (nell'immagine di sotto ad esempio un CSV esposto da Google Drive), assicurandoti che la prima colonna contenga i geo-code o i nomi delle regioni (questo ovviamente nel caso di una mappa per regioni);

• quando farai click su “Proceed” potrai completare la personalizzazione della tua mappa.

Il risultato finale sarà una mappa che si aggiornerà automaticamente all’aggiornarsi della fonte.

NOTA BENE:

• nella tua tabella deve essere presente una colonna con i nomi dei luoghi (nel mio esempio le regioni) o i codici, coincidenti a quelli dei dataset di datawrapper (puoi copiarli e incollarli dalla tabella di preview). In questo modo verrà effettuato il JOIN tra i poligoni del layer di datawrapper e i tuoi dati;
• questa colonna deve essere la prima (è bene ripeterlo).

La mappa di sotto ad esempio, basata sul dataset ISTAT “Aspetti della vita quotidiana : Stato di salute — regioni e tipo di comune”, è generata in modalità live a partire da questa tabella.

PER FAVORE: NON FARE SOLO MAPPE SUL COVID-19

Nota bene: qualche settimana dopo la pubblicazione di questo post, è arrivata la guida ufficiale per farlo:
https://academy.datawrapper.de/article/236-how-to-create-a-live-updating-symbol-map-or-choropleth-map

Datawrapper: creare una mappa a partire da dati “live” was originally published in tantotanto on Medium, where people are continuing the conversation by highlighting and responding to this story.

L’atlante di QGIS consente di creare in modo automatico delle pubblicazioni cartografiche multipagina. Proprio come i meravigliosi volumi del tempo della scuola (in cui ci si perdeva letteralmente dentro, specie “noi di quando non c’era l’internèt”).

Esistono tanti tutorial sul come farne uno e non replicherò.
Documenterò qui una modalità atipica e di grande utilità. Lo spunto nasce da questo messaggio del mitico @puntofisso (Giuseppe Sollazzo), ricevuto a fine dicembre 2019:

Borruso, sei la mia ultima speranza, stile Obi Wan. Posso chiederti una domanda QGIS? :) Sto cercando di capire se mi stia sfuggendo qualcosa di ovvio sulla generazione di un Atlas […]. Ho un dataset che mostra le precipitazione storiche in UK in termini di decili, per mostrare quando le precipitazioni sono state fuori dal comune.[…]
Posso trasporlo, splittare le coordinate (in easting, northing), etc… ma faccio fatica a pensare a una struttura che sia compatibile con Atlas per cui possa alla fine ottenere una mappa per anno

A partire da un dataset come quello soprastante Giuseppe vuole realizzare un atlante con una pagina per anno, senza cambiare l’area geografica. Al variare di pagina, fare riferimento a un diverso campo, a una diversa fonte dati.
Di solito avviene proprio il contrario: la variazione di pagina è in base alle righe di un layer (ad esempio al cambio di valore nella colonna “nome regione”), con taglio geografico diverso per ognuna.

L’obiettivo è quello di ottenere qualcosa come quella di sotto (le n pagine corrispondenti agli anni dei dati sorgente).

Come creare l’atlante

Il punto di partenza è quello di caricare i layer cartografici di base. Nell'esempio creato per descrivere la procedura faccio uso soltanto di un file CSV (denominato decilesNE.csv), con una struttura come quella della tabella in alto, più due colonne con le coordinate (non sono riuscito a risalire al sistema di coordinate, ho usato EPSG:32632).
Una volta aggiunto, si formatta il layer sfruttando i valori di una delle colonne year..., che contengono valori da 0 a 9.
Sono dieci categorie a cui ad esempio associare dei toni dal bianco al blu, al passaggio da 0 a 9. Si tratta di punti che rappresentano quadrati di 5 chilometri e a cui verrà associato pertanto un simbolo quadrato di 5 km di lato.

Per poter “paginare” per anno verrà aggiunto un layer con i valori distinti per anno, ma disposti in righe. Basta quindi trasporre la riga con il nome delle colonne del file di input e ottenere un CSV come questo (che ho chiamato at.csv).
Una volta creato, si dovrà aggiunge al progetto.

Si potrà iniziare dal creare un nuovo layout e attivare in questo la generazione dell’atlante.
Due i parametri importanti in questa fase:

1. quale layer farà da coverage, ovvero quale "comanderà" il cambio pagina. Per questo progetto sarà at;
2. quale campo utilizzare per ordinare le pagine. Sarà item, l'unico del file at.csv, che contiene l'elenco dei nomi delle colonne (che corrispondono agli anni).

In questo modo QGIS saprà come “paginare” (vedi sotto).

Rimane però un PROBLEMA: come fare in modo che al cambio pagina la formattazione della mappa sia basata sul campo corrispondente all'anno selezionato?

Al momento la formattazione è basata su una delle colonne del layer. Quindi al cambio pagina non varierà.
QGIS espone alcune variabili legate agli atlanti, una è denominata @atlas_pagename e contiene il valore del nome della pagina corrente.
Quindi basterebbe usare questa variabile al posto del nome della colonna, e fare in modo che il layer spaziale vari al cambio di pagina.

Per farlo bisogna aprire il pannello di formattazione del layer e sostituire il nome del campo con l’espressione eval(@atlas_pagename) (vedi sotto).
Si usa la funzione eval per fare in modo che il valore della variabile (ad esempio year1958) sia mappato come nome campo e non come stringa.

A questo punto si potrà esportare l’atlante in PDF, in SVG e in PNG.

Ho creato una copia del file progetto con le risorse necessarie (creato con QGIS 3.4), in modo che sia possibile testare il tutto in autonomia. All'interno il layout con attivo l’atlante, denominato miSollazzo.

NOTA BENE: all'apertura del progetto il layer non sarà visualizzato, perché l'atlante non è attivo e quindi QGIS non sa quale sia la pagina e di conseguenza quale il campo da usare per formattare il layer. Allora, una volta aperto il progetto, si dovrà aprire il layout presente e attivare l'anteprima dell'atlante.
Il sistema di coordinate del layer è impostato in modo casuale, perché non sono riuscito a trovare informazioni sullo Spatial Reference System dei dati sorgente. Allo scopo dell’articolo, non è per fortuna qualcosa di bloccante.

Come farlo con pandas e matplotlib

Giuseppe Sollazzo ha lavorato su questi dati utilizzando un Jupyter Notebook Python, basato su Pandas e pyplot — Matplotlib.
La sua curiosità sul come farlo in QGIS nasce proprio da questi suoi “lavori in corso”.

L’ha ovviamente ultimato e lo trovate qui https://github.com/puntofisso/Notebooks/Rainfall.ipynb

Voglio bene a questo notebook, perché è stata l’occasione per “to play” con puntofisso 🥰

Divertissement

Durante la preparazione di questo testo mi sono anche divertito a utilizzare l’output dell’atlante per creare due cose:

• la griglia di immagini presente all'inizio;
• una gif animata che mostri le variazioni nel tempo.

Per farlo, ho esportato le pagine dell’atlante in PNG e le ho elaborate con ImageMagick.

La griglia l’ho creata con l’utility montage, con il comando:

montage output_*.png -tile 10x6 -geometry +5+5 out.png

In cui con -tile 10x6 si imposta la dimensione della griglia (che qui è 6 righe per 10 colonne).

La gif animata invece sfruttando l’utility convert, con il comando:

convert -delay 100x100  output_*.png -loop 0 movie.gif

QGIS: generare un atlante basato su campi was originally published in tantotanto on Medium, where people are continuing the conversation by highlighting and responding to this story.

Quando si lavora con dati spaziali uno degli obiettivi può essere quello di creare nuove geometrie in base a come si sovrappongono fra loro. Queste procedure sono spesso indicate usando il linguaggio degli insiemi: intersezione, unione, differenza, ecc..

Qui sotto un’immagine con alcune delle operazioni “classiche” (questa e il testo introduttivo dalla guida di GeoPandas).

In questo articolo verrà descritta una di queste e in particolare l’unione. È un’operazione realizzabile in diverse modalità e con diversi strumenti, ma diventa abbastanza complessa e/o non realizzabile in modo diretto se l’obiettivo è:

• realizzare l’unione tra oggetti di uno stesso layer di input;
• produrre in output non soltanto l’unione delle geometrie, ma anche quella degli attributi, sempre in base alla sovrapposizione delle geometrie.

Nell'immagine di sotto un esempio: a sinistra le geometrie di input e a destra l’output. I perimetri delle aree fanno da “lama” di ritaglio e laddove le aree si sovrappongono devono essere trasferiti gli attributi di input (in questo esempio i valori A, B e C); eventuali poligoni sovrapposti e coincidenti sono aggregati in un’unica geometria.

Un’applicazione con cui è molto semplice realizzare quanto descritto è il fantastico mapshaper. È uno strumento scritto in JavaScript, per modificare file Shapefile, GeoJSON, TopoJSON, CSV e altri formati, molto noto per le sue eccellenti doti nella semplificazione di geometrie.
Fa molto molto di più e questa è l’occasione per mostrare un esempio.

A Totò Fiandaca (grazie mille) si deve la redazione del testo che descrive come fare la cosa con QGIS e tramite SQL geografico (in particolare con SpatiaLite).

Unione dei poligoni di un layer con mapshaper

Una volta installato è utilizzabile o a riga di comando o tramite un’interfaccia web (questa la versione pubblica ufficiale).

Per replicare l’operazione descritta a seguire è stato predisposto questo file GeoJSON di input d’esempio, che ha un solo campo (id) ed è composto da tre poligoni (con valori di id pari ad A, B e C).

Questo è il comando tipo da lanciare da shell:

mapshaper -i ./inputLayer.geojson -mosaic calc='output=collect(id).toString()' -o ./output.geojson

Per punti:

• con -i si imposta il file di input;
• poi il comando da eseguire — -mosaic — che in questo caso è appunto l'unione dei poligoni;
• al comando mosaic si associa l'opzione calc, che consente di eseguire calcoli nelle aggregazioni molti-a-uno, utilizzando espressioni JavaScript. In questo caso viene generato un nuovo campo denominato output, in cui vengono aggregati i valori del campo id nei poligoni di intersezione risultanti;
• infine si definisce con -o il file di output.

Nell'immagine di sotto è rappresentato il processo. Si veda ad esempio come al poligono centrale, frutto dell’intersezione geometrica di tutti e tre i poligoni di input, venga associato il valore A,B,C.

In termini di geometrie si avranno 7 record di output (corrispondenti alle 7 intersezioni geometriche possibili) con questi attributi:

• A,B,C
• B,C
• A,C
• C
• A
• A,B
• B

Tantissimi i campi di applicazione di un processo come questo. Questi tre poligoni potrebbero essere gli areali “colpiti” dalla diffusione di inquinanti da tre punti sorgente, oppure le aree che distano una determinata misura da una fontanella d’acqua potabile, ecc..
L’unione, con intersezione geometrica e aggregazione di attributi, consente di produrre un output in cui è possibile leggere per ogni area il contributo di tutti i poligoni di input.

Altre modalità

Come dicevo in apertura, tutto questo è fattibile in altre modalità, ma non mi risulta ce ne sia una così semplice e diretta in altri ambienti. A seguire qualche esempio.

GeoPandas

GeoPandas è un famoso progetto open source per semplificare il lavoro con i dati geospaziali in Python. È un gran bel modulo!

Consente di eseguire il processo di unione (e altri), ma di base è pensato per essere eseguito tra due layer distinti di input.
A seguito di una issue in merito, aperta sul repository di progetto, è stata fornita questa risposta, in cui si evidenzia come non ci sia un modo diretto per farlo:

[…]there is right now no direct way to do this in geopandas. You could do a “cumulative” intersection manually rather easily in a loop, but that would need some more complexity to keep track of the values.[…]

Per l’occasione è stato creato un bel notebook in cui è mostrato come sia possibile eseguire l’intersezione tra poligoni sovrapposti; da notare la maggiore complessità e sopratutto come comunque non sia gestita l’aggregrazione molti-a-uno degli attributi (nel senso che c’è da scrivere questa parte di codice).

QGIS

A partire dalla versione QGIS 3.2 Bonn è disponibile il geoalgoritmo Union che svolge la parte geometrica di ritaglio del processo descritto sopra, ma non raggruppa le geometrie omologhe e di conseguenza gli attributi.
Una singola area condivisa tra n poligoni sarà ritagliata e produrrà in output n elementi geometrici. Nell'immagine di sotto ad esempio, in corrispondenza dell'area selezionata in giallo, tre record e tre geometrie; complessivamente così verranno quindi prodotti 12 record.

Per raggruppare le geometrie e gli attributi occorre un altro passaggio, oppure creare un modello grafico, come quello che trovate qui; in questo si utilizzano anche lo snapping e la riparazione delle geometrie, che servono in casi molto più complessi di quello attuale (vedi issue correlata).

La query utilizzata nel modello è la seguente:

SELECT group_concat ("id") AS ID, st_union (geometry) AS geometry
FROM input1
GROUP BY geometry

Fatto il raggruppamento si avranno in output ancora una volta 7 record.

SQL

Non esiste una funzione di SQL spaziale (SpatiaLite o PostGIS) che risolva in modo diretto la cosa; occorre analizzare step by step il da farsi:

1. estrarre i perimetri dai poligoni di input;
2. estrarre i punti di intersezione dei perimetri delle varie feature;
3. aggiungere i punti alla geometria dei perimetri;
4. “splittare” la geometria dei perimetri;
5. “poligonalizzare” a partire dai perimetri “splittati”;
6. associare gli attributi alle varie geometrie “poligonalizzate”.

Ad esempio in ambiente SpatiaLite si potrebbe partire con l’importazione del file di esempio:

# importa file inputLayer e crea file spatialite
ogr2ogr -append -t_srs EPSG:4326 -f SQLite ./dbOverlayUnion.sqlite ./inputLayer.geojson -nln "inputLayer" -dsco SPATIALITE=YES

E infine lanciare la seguente query:

-- algoritmo proiezioni, the SpatiaLite way
-- di Salvatore Fiandaca
-- e-mail: pigrecoinfinito@gmail.com
-- è stato di aiuto: http://blog.cleverelephant.ca/2019/07/postgis-overlays.html

-- crea geotabella polygonize
CREATE TABLE polygonize AS
SELECT St_Polygonize(t.geom) as geom
FROM
(SELECT id, St_Union(st_boundary(geometry)) as geom
FROM inputLayer) t;
SELECT RecoverGeometryColumn('polygonize','geom',4326,'MULTIPOLYGON','XY');

-- crea geotabella dalle componenti elementari della geotabella polygonize
SELECT DropGeoTable('elementi');
SELECT ElementaryGeometries( 'elementi' ,'geom' , 'polygonize' ,'out_pk' , 'out_multi_id', 1 ) as num, 'polygon splitted' as label;

-- crea poligoni di output con attributi
SELECT DropGeoTable( "OUTPUT");
CREATE TABLE OUTPUT AS
SELECT Group_Concat(id) as id, e.geom
FROM inputLayer p, elementi e
where st_intersects (ST_PointOnSurface(e.geom), p.geometry) = 1
GROUP BY e.geom;
SELECT RecoverGeometryColumn('OUTPUT','geom',4326,'POLYGON','XY');

-- aggiorna statistiche e VACUUM
UPDATE geometry_columns_statistics set last_verified = 0;
SELECT UpdateLayerStatistics('geometry_table_name');
VACUUM;

In output ancora una volta 7 record.

Nota bene: Totò ha poi dettagliato tutto in questo post https://pigrecoinfinito.com/2020/01/02/spatialite-overlays/

SAGA

Sezione a cura di Ludovico Frate (grazie).

In SAGA GIS è possibile utilizzare il comando Polygon Self-Intersection per ottenere 7 geometrie (record) derivanti dall'unione del vettore originale.
I passaggi da eseguire sono:

• Importazione del file geojson utilizzando il comando Import Shapes (Tools -> Import/Export->GDAL/OGR->Import Shapes);

• Eseguire l’unione eseguendo il comando Polygon Self-Intersection (Tools -> Shapes->Polygons->Polygons Self-Intersection) ed inserire in Polygons il geojson appena caricato, in identifier la colonna id.

In output i soliti 7 record.

R

Sezione a cura di Andrea Zedda (grazie). Qui il notebook.

Carico la libreria sf, una libreria che manipola geometrie spaziali utilizzando lo standard OGC “Simple feature”.
Sf è una libreria relativamente nuova in ambiente R ma data la semplicità d’uso, la versatilità e la piena compatibilità con la sintassi “tidy” è già divenuto lo standard incontrastato.

library(sf)

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

library(tidyverse)

## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0

## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

s <- st_read("inputLayer.geojson")
plot(s)

In in sf il concetto di unione è intersezione sembrano essere stati interpretati diversamente da geopandas e altri software. Per ottenere l’input lancio quindi semplicemente la funziona st_intersection che genera le features derivanti dalle sovrapposizioni e produce di default 2 campi: n.overlap che calcola quante volte la geometria si sovrappone con altre geometrie, e origins, una lista che contiene la combinazione degli indici delle features che overlappano.

i <- s %>% st_intersection()
i

## Simple feature collection with 7 features and 3 fields
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: 0 ymin: 0 xmax: 3 ymax: 3
## epsg (SRID):    NA
## proj4string:    NA
##     id n.overlaps origins                       geometry
## 1    A          1       1 POLYGON ((1.5 0.1339746, 1....
## 1.1  A          2    1, 2 POLYGON ((0.1339746 1.5, 0....
## 2    B          1       2 POLYGON ((2 2, 2 2, 2 2, 1....
## 1.2  A          2    1, 3 POLYGON ((2 1, 1.99863 0.94...
## 1.3  A          3 1, 2, 3 GEOMETRYCOLLECTION (LINESTR...
## 2.1  B          2    2, 3 POLYGON ((2 2, 1.99863 1.94...
## 3    C          1       3 MULTIPOLYGON (((3 1, 2.9986...

L’output generato sembra corretto, sono state prodotte 7 features tante quante sono le aree che si sovrappongono. Osservando la colonna delle gemotrie notiamo che non tutte sono dei poligoni ma abbiamo un multypolygon e addirittura una geometrycollection.

i %>% 
  ggplot()+
  geom_sf(aes(fill=row.names(i)), show.legend = F)+
  facet_wrap(.~row.names(i), )+
  theme_bw()

Questo è dovuto ad una delle caratteristiche dello standard simple features che per velocizzare alcune operazione elabora delle coordinate arrotondate. È un caso abbastanza raro ma quando si fanno delle operazioni unarie sulle geometrie è meglio settare la precisione.

i <-  st_intersection(st_set_precision(s, 10000))
i

## Simple feature collection with 7 features and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 0 ymin: 0 xmax: 3 ymax: 3
## epsg (SRID):    NA
## proj4string:    NA
##     id n.overlaps origins                       geometry
## 1    A          1       1 POLYGON ((1.5 0.134, 1.454 ...
## 1.1  A          2    1, 2 POLYGON ((0.134 1.5, 0.1613...
## 2    B          1       2 POLYGON ((1.5 1.866, 1.454 ...
## 1.2  A          2    1, 3 POLYGON ((2 1, 1.9986 0.947...
## 1.3  A          3 1, 2, 3 POLYGON ((1.5 1.866, 1.5446...
## 2.1  B          2    2, 3 POLYGON ((2 2, 1.9986 1.947...
## 3    C          1       3 POLYGON ((3 1, 2.9986 0.947...

Risolto il problema possiamo infine eseguire una seconda operazione che ci consente di portarci nell’output le elaborazioni sugli attributi. In questo caso inserisco nel pipe una funzione di mutate che crea una colonna output contenente una stringa degli attributi concatenati. Per far ciò si è fatto riferimento al campo origins che contiene come detto la posizione delle geometrie nell’input iniziale. Essendo origins una lista si estraggono i valori utilizzando la funzione sapply.

out <-  i %>% mutate(output=sapply(origins,
                     function(x) paste0(as.character(s$id)[x], collapse = ",")
))

Si ottiene quindi l’output desiderato.

out

## Simple feature collection with 7 features and 4 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 0 ymin: 0 xmax: 3 ymax: 3
## epsg (SRID):    NA
## proj4string:    NA
##   id n.overlaps origins                       geometry output
## 1  A          1       1 POLYGON ((1.5 0.134, 1.454 ...      A
## 2  A          2    1, 2 POLYGON ((0.134 1.5, 0.1613...    A,B
## 3  B          1       2 POLYGON ((1.5 1.866, 1.454 ...      B
## 4  A          2    1, 3 POLYGON ((2 1, 1.9986 0.947...    A,C
## 5  A          3 1, 2, 3 POLYGON ((1.5 1.866, 1.5446...  A,B,C
## 6  B          2    2, 3 POLYGON ((2 2, 1.9986 1.947...    B,C
## 7  C          1       3 POLYGON ((3 1, 2.9986 0.947...      C

out["output"] %>% plot

out %>% 
  ggplot()+
  geom_sf(aes(fill=output), show.legend = F)+
  facet_wrap(.~output)+
  theme_bw()

QGIS (seconda modalità)

Sezione a cura di Ivano Giuliano (grazie).

L’esempio grafico da risolvere mostrato in quest’articolo, è una rappresentazione di quello che viene definito nella teoria degli insiemi un diagramma di Venn, ovvero un diagramma che mostra tutte le possibili relazioni logiche tra un numero finito di insiemi (le tre circonferenze), che vengono a costituire le regioni di interesse. Il metodo si basa proprio su questo concetto, individuazione delle parti di sovrapposizione di regioni e quindi l’insieme di punti in esse ricadenti, nonchè le regioni di provenienza, che andranno a costituire 12 elementi nella prima fase e 7 nella seconda.

Algoritmo Union —> A-B-B,A-C-C,A-C,B-C,B,A (12 elementi)

Algoritmo Aggregate —> {A} - {B} - {B,A} - {C} - {C,A} - {C,B} - {C,B,A} ( 7 insiemi)

Il metodo che si andrà ad esporre qui di seguito, è eseguibile completamente dall’interno di QGIS che prevede due fasi sequenziali con i seguenti geoprocessi.

Geoalgoritmo Union

Riprendendo la definizione : Questo algoritmo controlla le sovrapposizioni tra le features all'interno del livello di Input (nell'esempio : polyover) e crea features separate per quelle parti sovrapposte e non sovrapposte. Ovvero per le parti di aree sovrapposte, verranno create tante features quante sono le effettive parti di sovrapposizione. La tabella degli attributi del livello Union che verrà generato, conterrà tutti i valori di attributo delle features originali per le parti non sovrapposte e valori di attributo per le parti sovrapposte, in totale 12 elementi come indicato qui sopra. L’ algoritmo prevede inoltre la possibilità di utilizzare un livello di sovrapposizione e quindi tutte le sue derivazioni, ma che in questo caso non vengono prese in considerazione.

Geoalgoritmo Aggregate

Riprendendo la definizione : Questo algoritmo a partire da un livello vettoriale in ingresso o tabella, permette l’ aggregazione di features tramite una espressione di raggruppamento Group By Expression . Le features che secondo l’espressione restituiranno lo stesso valore, saranno aggregate. Questo rappresenta la chiave di questa risoluzione al quesito proposto, tramite QGIS. L’algoritmo è molto ben strutturato e permette di effettuare tutta una serie di raggruppamenti mediante modalità parametriche diverse, in questo esempio, viene utilizzata la seguente espressione di raggruppamento:

geom_to_wkt($geometry)

classica funzione utilizzabile nel field calculator che permette in modalità Well-Known-Text (WKT), di rappresentare la geometria . La risultante del geoprocesso, sarà che le geometrie riconosciute, saranno combinate in una geometria multipart per ogni gruppo creato. I relativi attributi di output, verranno cosí calcolati sulla base di ciascuna aggregazione definita.

NOTA BENE: Salvatore Fiandaca ha realizzato (grazie) questo modello grafico QGIS, per questa modalità

Per concludere

L’autore di mapshaper è Matthew Bloch. È uno sviluppatore con doti al di fuori del comune, sia in termini professionali, che in termini di attitudine all'ascolto e allo scambio con gli altri.
Questo comando di unione è stato introdotto da poco (il 18 novembre 2019), in risposta a una richiesta ricevuta nel repository del progetto (da sottolineare che anche precedentemente fosse semplice realizzare la cosa).

È evidente come mapshaper sia (tra quelle di questo articolo) la soluzione più diretta e semplice. Lo è anche perché è stata costruita una funzione specializzata per svolgere questo compito; vista l’utilità dovrebbe/potrebbe essere disponibile anche in altri ambienti.

Non vi resta che installarlo e apprezzarne anche tutte le altre ottime caratteristiche.

Un grosso grazie a Matthew per quello che realizza e mette a disposizione e per il modo in cui lo fa!

“Unire” i poligoni di un layer con grande semplicità: è un lavoro (soltanto?) per mapshaper was originally published in tantotanto on Medium, where people are continuing the conversation by highlighting and responding to this story.