Il decreto “di Natale”, in chilometri

Pubblicato il da Totò

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!!!

In grigio i Comuni con al massimo 5.000 abitanti, in blu quelli oltre i 5.000 abitanti – dati ISTAT 2020

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):

Modello grafico di QGIS

gli algoritmi usati sono:

  1. Estrai elementi per attributi: per filtrare i comuni con meno di 5.000 abitanti;
  2. Estrai elementi per attributi: per estrarre i comuni capoluogo di provincia;
  3. Dissolvi: per fondere tutti i comuni e ottenere i confini nazionali;
  4. Buffer: per creare area di influenza entro 30 km per ogni punto del confine amministrativo del comune;
  5. Differenza: per togliere dai buffer (30 km) i comuni capoluoghi;
  6. Intersezione: per sagomare tutti i buffer rispetto i confini nazionali;
  7. Riorganizzazione campi: per ordinare e selezionare solo gli attributi che ci servono;
  8. 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):

I Comuni in Blu hanno una popolazione oltre i 5.000 abitanti, grigio meno di 5.000
Buffer da 30 km senza comune capoluogo e ritagliato secondo il limite amministrativo siciliano

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:

spatiaLite_gui 2.1 con spatialite 5.0

alcuni dati statistici:

  1. il comune che ha la possibilità di raggiungere più comuni è Pontida (427);
  2. il comune che può raggiungere il minor numero di comuni è Ventotene (1).
Atlas di QGIS – Comune di Pontida BG
Atlas di QGIS – Comune di Ventotene LT

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 arancione l’area di tutti i comuni raggiungibili

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

  1. Avviare progetto QGIS (trovi tutto qui);
  2. dal pannello dei layer, selezionare il layer Aree30cappa;
  3. aprire la tabella attributi F6 (selezionate l’icona porta in cima);
  4. pigiare il tasto funzione F3, comparirà una maschera di ricerca;
  5. cercare un comune e pigiare su Seleziona Elementi;
  6. tasto destro mouse sulla riga selezionata nella tabella attributi e cliccare su Imposta come elemento...
  7. si avvierà il compositore di stampe con i dati sul Comune selezionato:
come usare Atlas

Esempio di ricerca:

come fare una ricerca per comune

ecco l’Atlas:

Pagina del Comune Campione d’Italia

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

RINGRAZIAMENTI

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”

  1. Pingback: Il decreto ‘di Natale’ sugli spostamenti piccoli comuni | pigrecoinfinito

Rispondi

Inserisci i tuoi dati qui sotto o clicca su un'icona per effettuare l'accesso:

Gravatar
Logo di WordPress.com

Stai commentando usando il tuo account WordPress.com. ( Chiudi sessione /  Modifica )

Google photo

Stai commentando usando il tuo account Google. ( Chiudi sessione /  Modifica )

Foto Twitter

Stai commentando usando il tuo account Twitter. ( Chiudi sessione /  Modifica )

Foto di Facebook

Stai commentando usando il tuo account Facebook. ( Chiudi sessione /  Modifica )

Annulla

Connessione a %s...

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

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