L’elaborazione dei dati per la mappa del decreto “di Natale”, in chilometri

La versione “a riga di comando”

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 🤣