In questo blog post spiego gli strumenti usati per la mia versione di #30cappa del decreto “di Natale” e segue quanto descritto in questo blog post di lancio (invito a leggerlo prima di continuare) a cura di aborruso, napo e pigreco (lo scrivente).
EDIT: hanno aggiornato i dati sulla popolazione ISTAT il 21/12/2020, quindi i dati di questo blog non sono allineati (ci sono circa 25 comuni in più sotto i 5000 abitanti) – a breve aggiorno!!!

Ho realizzato le Aree30cappa
utilizzando due software GIS diversi: QGIS 3.16 Hannover e SpatiaLite 5.0 NextGen
QGIS
QGIS è un software GIS desktop Libero e Open Source che Crea, modifica, visualizza, analizza e pubblica informazioni geospaziali su Windows, Mac, Linux, BSD e dispositivi mobili.
QGIS permette di realizzare complesse mappe colorate e di fare analisi con potenti strumenti senza scrivere una linea di codice, un esempio sono gli Strumenti di Processing che hanno migliaia di geo-algoritmi pronti all’uso; io utilizzerò, per questo caso, il modellatore grafico per mettere assieme e in cascata gli algoritmi necessari per ottenere le aree30cappa
, ovvero i buffer di 30 km attorno ai confini dei comuni con al massimo 5.000 abitanti, togliendo i comuni capoluogo di provincia e i limiti nazionali.
Avviando il modello comparirà la seguente finestra:
che richiede solo un Vettore di input (caricato in QGIS o presente nel nostro PC) da selezionare o da ricercare nel PC; in output possiamo salvare il file in qualsiasi formato (gestito da QGIS: shapefile, KML, geojson ecc…) o lasciarlo temporaneo; dopo Esegui
comparirà la maschera:
il risultato è:
tante bolle (buffer da 30 km), i buchi sono i capoluogo di provincia, che non possono essere raggiunti secondo il decreto.
Vediamo cosa c’è dentro il modello (QGIS 3.16.2 Hannover):
gli algoritmi usati sono:
- Estrai elementi per attributi: per filtrare i comuni con meno di 5.000 abitanti;
- Estrai elementi per attributi: per estrarre i comuni capoluogo di provincia;
- Dissolvi: per fondere tutti i comuni e ottenere i confini nazionali;
- Buffer: per creare area di influenza entro 30 km per ogni punto del confine amministrativo del comune;
- Differenza: per togliere dai buffer (30 km) i comuni capoluoghi;
- Intersezione: per sagomare tutti i buffer rispetto i confini nazionali;
- Riorganizzazione campi: per ordinare e selezionare solo gli attributi che ci servono;
- vettore finale!
come visibile nello screenshot di sopra, QGIS impiega circa 5 minuti per ultimare il processo (il mio laptop è del 2015).
Ecco un esempio su un Comune molto allungato (Petralia Sottana – PA):
A questo punto mi sono chiesto, ma quanti comuni sono ‘toccati’ da questo Buffer da 30 km? e quanta popolazione ci rientra? Ecco che nasce l’Atlas!!!
Atlas di QGIS
Per poter rispondere alle domande poste sopra occorre, per ogni poligono buffer da 30 km (sono circa 5.500), calcolare il numero di comuni che intersecano il buffer e “catturare” i loro attributi, per esempio la popolazione.
È un processo molto oneroso e la tabella di output ha oltre 800.000 righe.
Ho scelto, questa volta, di usare SpatiaLite e dopo aver creato un database e importato i dati (Buffer da 30 km e i Comuni ISTAT) ho lanciato questo script:
-- -- the present SQL script is intended to be executed from SpatiaLite_gui -- SELECT InitSpatialMetadata(1); -- -- il vettore Com01012020_g_WGS84 contiene già il campo Abitanti demoistat2020 -- in OUTPUT crea le due tabelle di sotto -- -- crea tabella - con attributi su tutti i comuni per ogni aree30cappa CREATE TABLE comune_e_comuni AS SELECT b.cod_reg,b.cod_prov,b.pro_com_t,b.comune,b.Abitanti AS ABITANTI, p.cod_reg AS COD_REG_2,p.pro_com_t AS PRO_COM_T_2,p.comune AS COMUNE_2 FROM Com01012020_g_WGS84 b, Aree30cappa p WHERE ST_intersects (b.geom , p.geom) = 1 AND b.ROWID IN ( SELECT ROWID FROM SpatialIndex WHERE f_table_name = 'Com01012020_g_WGS84' AND search_frame = p.geom) AND b.CC_UTS != 1 AND b.comune != p.comune; -- -- crea tabella - per ogni aree30cappa lista dei comuni ed altro CREATE TABLE lista_nro_comuni_x_comune5k AS SELECT cc."COD_REG_2" AS COD_REG, cc."PRO_COM_T_2" AS PRO_COM_T, cc."COMUNE_2" AS COMUNI, group_concat(distinct cc."COD_REG_2") AS lista_COD_REG_2, group_concat(cc."PRO_COM_T") AS lista_PRO_COM_T, group_concat(cc."COMUNE") AS lista_comuni, sum(cc."Abitanti") AS tot_abitanti, count(*) AS nro_comuni, CastToMultiPolygon(ST_UNION (c.geom)) as geom FROM "comune_e_comuni" cc join "Com01012020_g_WGS84" c USING (PRO_COM_T) GROUP BY 3 SELECT RecoverGeometryColumn('lista_nro_comuni_x_comune5k','geom',32632,'MULTIPOLYGON','XY'); -- -- cancella geotabelle inutili DROP TABLE IF EXISTS comune_e_comuni; -- -- aggiorno statistiche e VACUUM (nel mio vecchio laptop impiega circa 7 minuti) -- UPDATE geometry_columns_statistics set last_verified = 0; SELECT UpdateLayerStatistics('geometry_table_name'); VACUUM;
che mi restituisce una tabella con circa 5.500 righe (i comuni con meno di 5.000 abitanti) e i dati cercati: numero comuni raggiungibili, lista dei comuni, e totale popolazione residente nei comuni raggiunti:
alcuni dati statistici:
- il comune che ha la possibilità di raggiungere più comuni è Pontida (427);
- il comune che può raggiungere il minor numero di comuni è Ventotene (1).

Questo ultimo processo (lungo e noioso) è necessario solo se volessimo esportare i dati, ma il mio obiettivo è quello di creare un Atlas e quindi posso far calcolare al volo i dati che mi occorrono, demandando il computo all’Atlas stesso e alle incredibili espressioni di QGIS.
Ho utilizzato come Layer di copertura
dell’Atlas il vettore Aree30cappa
(i buffer da 30km bucati) e come nome pagina il campo comune
, per calcolare il numero di comuni ‘toccati‘ basta utilizzare questa espressione:
array_length(overlay_intersects('comuni',$geometry, filter:="cc_uts"!=1))-1
elenco dei comuni, la seguente espressione:
array_to_string( array_filter( overlay_intersects('comuni',"comune", filter:="cc_uts"!=1), @element!=@atlas_pagename))
somma degli abitanti:
array_sum(overlay_intersects('comuni',"abitanti",filter:= "cc_uts" != 1 ))
inoltre posso creare nell’Atlas, sempre al volo e con le espressioni, una collezione di geometrie (i Comuni) intersecantesi con l’area di buffer da 30 km e ottenere l’effetto di sotto:

in questo caso ho aggiunto un simbolo (layer di stile) usando il generatore di geometrie e l’espressione:
with_variable('com30cappa', collect_geometries( overlay_intersects( 'comuni',$geometry, filter:="cc_uts"!=1)), intersection( bounds(@com30cappa),@com30cappa))
Come usare Atlas di QGIS
- Avviare progetto QGIS (trovi tutto qui);
- dal pannello dei layer, selezionare il layer
Aree30cappa
; - aprire la tabella attributi F6 (selezionate l’icona porta in cima);
- pigiare il tasto funzione F3, comparirà una maschera di ricerca;
- cercare un comune e pigiare su Seleziona Elementi;
- tasto destro mouse sulla riga selezionata nella tabella attributi e cliccare su
Imposta come elemento...
- si avvierà il compositore di stampe con i dati sul Comune selezionato:
Esempio di ricerca:
ecco l’Atlas:
l’output dell’Atlas è possibile stamparlo per singola pagina o in serie, per tutta Italia (5496 pagine!) o per una Provincia o per Regione, basta definire un filtro e poi ci pensa QGIS; le pagine possono essere esportate come immagini (jpg, png, SVG) o come PDF.
SpatiaLite
SpatiaLite è una libreria Open Source pensata per estendere il core SQLite per supportare le funzionalità Spatial SQL a tutti gli effetti (realizzata e curata da Alessandro Furieri).
Ho realizzato un unico script SQL che riproduce lo stesso workflow adottato in QGIS ottenendo gli stessi risultati, lo script è il seguente:
-- -- the present SQL script is intended to be executed from SpatiaLite_gui -- SELECT InitSpatialMetadata(1); -- -- il vettore Com01012020_g_WGS84 contiene già il campo Abitanti demoistat2020 -- in OUTPUT crea un vettore aree30cappa -- -- crea tabella - buffer da 30km su comuni da 5k SELECT DropGeoTable( "b30k_comuni5k"); CREATE TABLE b30k_comuni5k AS SELECT pro_com_t,cod_reg,comune,Abitanti, st_buffer (geom,30000) as geom FROM Com01012020_g_WGS84 WHERE abitanti <=5000; SELECT RecoverGeometryColumn('b30k_comuni5k','geom',32632,'POLYGON','XY'); SELECT CreateSpatialIndex('b30k_comuni5k', 'geom'); -- -- crea tabella - capoluoghi di provincia SELECT DropGeoTable( "capoluoghi_prov"); CREATE TABLE capoluoghi_prov AS SELECT pro_com_t,cod_rip,cod_reg,comune,geom FROM Com01012020_g_WGS84 WHERE cc_uts != 0; SELECT RecoverGeometryColumn('capoluoghi_prov','geom',32632,'MULTIPOLYGON','XY'); SELECT CreateSpatialIndex('capoluoghi_prov', 'geom'); -- -- crea tabella temporanea unendo i capoluoghi di provincia SELECT DropGeoTable( "tmp_cap_prov"); CREATE TABLE tmp_cap_prov AS SELECT CastToMultiPolygon(ST_UNION(geom)) AS geom FROM capoluoghi_prov; SELECT RecoverGeometryColumn('tmp_cap_prov','geom',32632,'MULTIPOLYGON','XY'); -- -- crea tabella temporanea unendo tutti i comuni SELECT DropGeoTable( "tmp_italia"); CREATE TABLE tmp_italia AS SELECT CastToMultiPolygon(ST_UNION(geom)) AS geom FROM Com01012020_g_WGS84; SELECT RecoverGeometryColumn('tmp_italia','geom',32632,'MULTIPOLYGON','XY'); -- -- clippa italia con i buffer da 30km e comuni 5k SELECT DropGeoTable( "b30k_comuni5k_italia"); CREATE TABLE b30k_comuni5k_italia AS SELECT b.pro_com_t,b.cod_reg,b.comune,CastToMultiPolygon(ST_Intersection(b.geom ,p.geom)) as geom FROM b30k_comuni5k b, tmp_italia p WHERE ST_intersects (b.geom , p.geom) = 1; SELECT RecoverGeometryColumn('b30k_comuni5k_italia','geom',32632,'MULTIPOLYGON','XY'); -- -- clippa i capoluoghi con i buffer da 30km e comuni 5k clippati con italia SELECT DropGeoTable( "aree30cappa"); CREATE TABLE aree30cappa AS SELECT b.pro_com_t,b.cod_reg,b.comune,CastToMultiPolygon(ST_DIFFERENCE(b.geom , p.geom)) as geom FROM b30k_comuni5k_italia b, tmp_cap_prov p WHERE ST_intersects (b.geom , p.geom) = 1 union SELECT b.pro_com_t,b.cod_reg,b.comune,CastToMultiPolygon(b.geom) as geom FROM b30k_comuni5k_italia b,tmp_cap_prov p WHERE ST_intersects (b.geom , p.geom) = 0; SELECT RecoverGeometryColumn('aree30cappa','geom',32632,'MULTIPOLYGON','XY'); SELECT CreateSpatialIndex('aree30cappa', 'geom'); -- -- cancella geotabelle inutili DROP TABLE IF EXISTS tmp_cap_prov; DROP TABLE IF EXISTS tmp_italia; DROP TABLE IF EXISTS b30k_comuni5k; DROP TABLE IF EXISTS capoluoghi_prov; DROP TABLE IF EXISTS tmp_cap_prov; DROP TABLE IF EXISTS tmp_italia; DROP TABLE IF EXISTS b30k_comuni5k_italia; -- -- aggiorno statistiche e VACUUM (nel mio vecchio laptop impiega circa 7 minuti) -- UPDATE geometry_columns_statistics set last_verified = 0; SELECT UpdateLayerStatistics('geometry_table_name'); VACUUM;
basta lanciarlo e dopo qualche minuto otterremo il vettore Aree30cappa
.
EDIT: (29/12/2020) sotto una nuova versione dello scriptSQL che usa due funzioni poco note e con poche tracce nel web: ST_CUTTER e ST_SUBDIVIDE:
-- -- the present SQL script is intended to be executed from SpatiaLite_gui 2.10 -- -- il vettore Com01012020_g_WGS84 contiene già il campo Abitanti demoistat2020 -- in OUTPUT crea un vettore aree30cappa -- -- crea tabella - buffer da 30 km su comuni da 5k DROP TABLE IF EXISTS "b30k_comuni5k"; CREATE TABLE "b30k_comuni5k" ("pk_uid" integer PRIMARY KEY autoincrement NOT NULL, "pro_com_t" text); -- aggiunge campo geom SELECT AddGeometryColumn ('b30k_comuni5k','geom',32632,'MULTIPOLYGON','XY'); -- INSERT INTO "b30k_comuni5k" (pk_uid, pro_com_t, geom) SELECT pk_uid, pro_com_t, CastToMultiPolygon(ST_Buffer (geom,30000)) AS geom FROM (SELECT pk_uid,pro_com_t, geom FROM "Com01012020_g_WGS84" WHERE abitanti <=5000); -- SELECT CreateSpatialIndex('b30k_comuni5k', 'geom'); -- -- crea tabella - capoluoghi di provincia DROP TABLE IF EXISTS "capoluoghi_prov"; CREATE TABLE "capoluoghi_prov" ("pk_uid" integer PRIMARY KEY autoincrement NOT NULL, "pro_com_t" text); -- aggiunge campo geom SELECT AddGeometryColumn ('capoluoghi_prov','geom',32632,'MULTIPOLYGON','XY'); -- popola la tabella INSERT INTO "capoluoghi_prov" (pk_uid, pro_com_t, geom) SELECT pk_uid, pro_com_t, geom FROM (SELECT pk_uid,pro_com_t,geom FROM Com01012020_g_WGS84 WHERE cc_uts != 0); -- -- crea tabella temporanea unendo tutti i comuni DROP TABLE IF EXISTS "tmp_italia"; CREATE TABLE "tmp_italia" ("pk_uid" INTEGER PRIMARY KEY); -- aggiunge campo geom SELECT AddGeometryColumn ('tmp_italia','geom',32632,'MULTIPOLYGON','XY'); -- popola la tabella INSERT INTO "tmp_italia" VALUES (1,NULL); -- aggiorna tabella UPDATE "tmp_italia" SET geom = (SELECT CastToMultiPolygon(ST_UNION(geom)) AS geom FROM "Com01012020_g_WGS84") WHERE pk_uid = 1; -- -- suddivide l'intero poligono nazionale DROP TABLE IF EXISTS "italia_subd2048"; CREATE TABLE "italia_subd2048" AS SELECT ST_Subdivide(geom,2048) AS geom FROM "tmp_italia"; SELECT RecoverGeometryColumn('italia_subd2048','geom',32632,'MULTIPOLYGON','XY'); -- -- estrae i poligoni elementari SELECT ElementaryGeometries( 'italia_subd2048', 'geom', 'italia_subd2048_elem', 'pk_elem', 'out_multi_id', 1 ) as num; -- -- taglia i buffer secondo il vettore dato SELECT ST_Cutter(NULL, 'b30k_comuni5k', NULL, NULL, 'italia_subd2048_elem', NULL, 'b30k_com5k_italy_subd_elem', 1, 1); -- -- ricompone i pezzi - ricorda che pk_id si riferisce alla tabella comuni 'Com01012020_g_WGS84' DROP TABLE IF EXISTS "b30k_com5k_italy"; CREATE TABLE "b30k_com5k_italy" ("pk_uid" INTEGER PRIMARY KEY); -- aggiunge campo geom SELECT AddGeometryColumn ('b30k_com5k_italy','geom',32632,'MULTIPOLYGON','XY'); INSERT INTO "b30k_com5k_italy" (pk_uid,geom) SELECT k."input_b30k_comuni5k_pk_uid" AS pk_uid, k.geom FROM (SELECT "input_b30k_comuni5k_pk_uid", "blade_italia_subd2048_elem_pk_elem", CastToMultiPolygon(st_union("geom")) AS geom FROM "b30k_com5k_italy_subd_elem" WHERE "blade_italia_subd2048_elem_pk_elem" is NOT NULL GROUP BY "input_b30k_comuni5k_pk_uid") k; SELECT RecoverGeometryColumn('b30k_com5k_italy','geom',32632,'MULTIPOLYGON','XY'); -- -- taglia i buffer sagomati con lo stivale secondo il vettore dei capoluoghi SELECT ST_Cutter(NULL, 'b30k_com5k_italy', NULL, NULL, 'capoluoghi_prov', NULL, 'b30k_com5k_italy_subd_elem_capolp', 1, 1); -- -- ricompone i pezzi - ricorda che pk_uid si riferisce alla tabella comuni 'Com01012020_g_WGS84' DROP TABLE IF EXISTS "finale"; CREATE TABLE "finale" ("pk_uid" INTEGER PRIMARY KEY); -- aggiunge campo geom SELECT AddGeometryColumn ('finale','geom',32632,'MULTIPOLYGON','XY'); INSERT INTO "finale" (pk_uid,geom) SELECT k."input_b30k_com5k_italy_pk_uid" AS pk_uid, k.geom FROM (SELECT "input_b30k_com5k_italy_pk_uid", CastToMultiPolygon(st_union("geom")) AS geom FROM "b30k_com5k_italy_subd_elem_capolp" WHERE "blade_capoluoghi_prov_pk_uid" IS NULL AND ST_Area(geom) > 1 group by "input_b30k_com5k_italy_pk_uid") k; SELECT RecoverGeometryColumn('finale','geom',32632,'MULTIPOLYGON','XY'); -- -- associa campo pro_com_t ed altri DROP TABLE IF EXISTS "aree30cappa"; CREATE TABLE "aree30cappa" AS SELECT f."pk_uid",c."pro_com_t",c."cod_reg",c."comune",c."abitanti",f."geom" FROM "finale" f left join "Com01012020_g_WGS84" c USING ("pk_uid"); SELECT RecoverGeometryColumn('aree30cappa','geom',32632,'MULTIPOLYGON','XY'); -- -- cancella geotabelle inutili DROP TABLE IF EXISTS "b30k_comuni5k"; DROP TABLE IF EXISTS "capoluoghi_prov"; DROP TABLE IF EXISTS "tmp_italia"; DROP TABLE IF EXISTS "italia_subd2048"; DROP TABLE IF EXISTS "italia_subd2048_elem"; DROP TABLE IF EXISTS "b30k_com5k_italy_subd_elem"; DROP TABLE IF EXISTS "b30k_com5k_italy"; DROP TABLE IF EXISTS "b30k_com5k_italy_subd_elem_capolp"; DROP TABLE IF EXISTS "finale"; -- -- aggiorno statistiche e VACUUM (nel mio vecchio laptop impiega circa 4 minuti) -- SELECT UpdateLayerStatistics('aree30cappa'); VACUUM;
VIDEO DEMO
in lavorazione …
NOTE FINALI: per chi volesse replicare le due procedure (con QGIS o/e SpatiaLite) trova tutti i dati e progetto nel repository messo a disposizione dall’Associazione OnData, nella cartella script/QGIS
(oppure scarica la cartella da qui).
Per risolvere il problema e semplificare un po’ la procedura, trovi nella cartella script/QGIS/dati
lo shapefile comuni.shp
(Com01012020_g_WGS84) che contiene già il campo Abitanti
(ISTAT 2020).
Nella cartella script/QGIS/output
trovi i risultati delle elaborazioni più il modello 30cappaAtlas.model
utilizzato per l’elaborazione dei dati (il modello è presente nel progetto, lo trovi negli strumenti di processing).
Per quanto riguarda la procedura SpatiaLite, crea una db vuoto, importa solo lo shapefile comuni.shp
rinominandolo in Com01012020_g_WGS84
; successivamente avvia lo script.sql
.
RIFERIMENTI
- Post di lancio su Tanto: https://medium.com/tantotanto/il-decreto-di-natale-in-chilometri-8af38744a7d5
- Sito 30cappa: https://ondata.github.io/30cappa/
- Popolazione ISTAT 2020: http://demo.istat.it/
- Limiti amministrativi ISTAT 2020 (generalizzati): https://www.istat.it/it/archivio/222527
- FAQ: https://ondata.github.io/30cappa/#faq
- QGIS: https://qgis.org/it/site/
- SpatiaLite: https://www.gaia-gis.it/fossil/libspatialite/index
- Decreto legge numero 172 del 18 dicembre 2020: https://www.gazzettaufficiale.it/eli/id/2020/12/18/20G00196/s
- 👉 Prova tu, scarica i dati e il progetto: https://drive.google.com/file/d/1AxZw8jxYokaDFRTKQ30qgm6hq_btwl1_/view?usp=sharing
RINGRAZIAMENTI
- Andrea Borruso – @aborruso
- Maurizio Napolitano – @napo
Licenza WTFPL (da una idea di @napo)
Vieni a trovarmi nel canale Telegram Pigrecoinfinito, ISCRIVITI!!! https://t.me/pigrecoinfinito
Un pensiero su “Il decreto “di Natale”, in chilometri”