Archivio degli autori ing.Massimiliano Moraca

Calcolo della superficie reale: confronto tra DTM

Qualche settimana fa nel gruppo GIS Italia fu posta questa domanda:

Ciao a tutti! Qualcuno sa come faccio a calcolare l’area reale (cioè l’area corretta in base all’altitudine) delimitata da un poligono in QGIS?

NB: Ho provato a calcolare l’area con il modulo di GRASS “r.surf.area”, dove viene inserito anche il dem, ma ottengo lo stesso valore che ottengo calcolando l’area senza tenere in conto l’altitudine.

Grazie

La risposta che ho dato a questa problematica è che il calcolo della superficie reale di un DEM è funzione sia della morfologia dell’area di interesse che della risoluzione del DEM stesso. In base a come variano questi due dati si può avere o non si può avere uno scostamento consistente tra il calcolo della superficie in piano e quella proiettata. La ragazza che ha posto la domanda in seguito alle risposte ricevute ha verificato che l’aria in esame era essenzialmente in piano per cui lo scostamento tra le superfici calcolate era molto esiguo.

Solo ora, con questo articolo, ho avuto la possibilità di dare una piccola dimostrazione di ciò che avevo scritto all’epoca.

Di seguito andrò a fare un confronto tra il DTM Ispra che ha la risoluzione di 20m/px ed il DTM da LiDAR della Città Metropolitana di Napoli che ha risoluzione 1m/px.

Il DTM dell’ISPRA copre tutta l’Italia ed è disponibile qui, mentre quello della Città Metropolitana di Napoli copra la sola area metropolitana del capoluogo campano ed è disponibile qui.

Il primo passo che ho effettuato è stato ricercare un’area essenzialmente pianeggiante; l’ho trovata nella area a nord di Napoli.

superficie dtm
Nota l’area ho eseguito un clip su entrambi i DTM che avevo a disposizione, poi usando il geoalgoritmo r.reclass ho definito tre fasce altimetriche. Preventivamente avevo analizzato i metadati dei raster scoprendo che il massimo ed il minimo di quota si aggirava, rispettivamente, su 35 ed i 336 metri. In base a questi dati ho diviso l’area nelle seguenti tre fasce altimetriche:

  • 35 – 100
  • 101 – 200
  • 201 – 336

Una volta fatto ciò ho preferito usare il risultato proveniente dal DTM LiDAR per gli studi successivi poichè più accurato. Dopo l’estrazione delle fasce quindi ho convertito il raster in vettore usando il geoalgoritmo Poligonize (Raster to Vector).

superficie dtm

Fatto ciò ho usato il geoalgoritmo r.surface.area per l’estrazione del raster delle superfici reali. Vi consiglio di guardare il video di Salvatore Fiandaca per vederne il funzionamento.

L’ultimo passaggio è stato usare la statistiche zonali per inserire nel vettore delle fasce altimentriche l’estensione reale delle tre fasce in base ad i due DTM.

Superficie (ha)
Fascia altimetrica DTM 20m/px DTM 1m/px Vettore
35 – 100 3.506 3.523 3.515
101 – 200 2.381 2.414 2.384
201 – 336 234 254 233
Totale 6.121 6.191 6.132

La precedente tabella comparativa va interpretata così: con il DTM 20m/px si tende a sottostimare la reale superficie delle aree, cosa che non avviene con il DTM 1m/px in cui si ha un risultato più preciso. Ricordo che l’area in esame è per lo più in piano ma nonostante ciò si ha una variazione totale di 59 ha se si confronta l’area totale del vettore e quella del DTM 1m/px; 1 ha è 10.000 m2, quasi 1 campo da calcio a 11 regolamentare.

A seguire il link da cui poter scaricare il progetto QGIS che contiene i file usati per questo articolo.

SS268, un problema non ancora risolto

La Strada Statale 268 è una superstrada che attraversa la zona est della Città Metropolitana di Napoli, come indica Wikipedia. E’ un lungo serpentone che abbraccia da nord-est il Vesuvio, è infatti detta anche Strada Statale del Vesuvio, ha una lunghezza di circa 27 Km. E’ detta anche strada della morte visti i continui incidenti mortali che si susseguono da decenni.

E’ anche e soprattutto una delle arterie individuate nel piano di evacuazione della zona rossa a seguito di una imminente eruzione del Vesuvio.

Oltre alle infrastrutture autostradali di cui si è detto, all’interno od in prossimità della Zona Rossa
vengono impegnati dai flussi di esodo i seguenti archi della rete stradale primaria regionale:
1) SS268 “del Vesuvio”
2) SS268 Raccordo
3) ex SS162dir “del Centro Direzionale”
4) SP1 “Circumvallazione di Napoli”
5) SS7bis variante “Asse di Supporto”
6) SS7IV “Domitiana”

Da Pianificazione di allontanamento – Delibera della Regione Campania n.8 del 17 gennaio 2017 pagina 19.

Fatta questa premessa, la finalità di questo articolo è andare a stimare la popolazione che risiede nel raggio di 1.000 m e 2.500 m da questa importante arteria stradale e cercarne di valutare l’efficienza e l’utilizzabilità in caso di eruzione del Vesuvio.

Per farlo ho estratto da OpenStreetMap la sola SS268, dal portale dell’ISTAT ho invece scaricato le celle censuarie e le tabelle della popolazione residente al 2011(ultimo anno di censimento ufficiale).

Ottenuti questi dati ho effettuato un buffer dalla SS268 con le distanze indicate in precedenza, quindi 1 e 2,5 KM. Successivamente con una operazione di selezione spaziale, tra i due buffer ottenuti verso le celle censuarie, ho estratto le sole celle censuarie che ricadono nelle aree di buffer in esame.

Da queste operazioni ne viene fuori che:

  • entro 1.000 m dalla SS268 risiedono 112.319 persone
  • entro 2.500 m dalla SS268 risiedono 325.874 persone

Considerando che, secondo i dati ISTAT, le famiglie italiane sono composte in media quasi da 3 persone, è stimabile in 37.400 e 108.625 il numero di famiglie che ricadono rispettivamente in un’area di raggio 1.000 m e 2.500 m dalla Strada Statale del Vesuvio.

Considerando che in ogni famiglia italiana c’è almeno un’automobile e considerando che la lunghezza media di una familiare è 4 m, è stimabile in 149.760 m e 434.500 m la lunghezza della coda che si genererebbe rispettivamente dalle automobili di coloro che risiedono in un’area di raggio 1.000 m e 2.500 m dalla Strada Statale del Vesuvio.

Dalle immagini che avete visto in precedenza, liberamente consultabili su Google Street View, è evidente che la superstrada in questione è ad una corsia per senso di marcia a parte un tratto di pochi km in cui ci sono i lavori per il suo ampliamento(lavori che vanno a singhiozzo da anni). Ciò significa questa importante arteria di circa 27 Km di lunghezza(54 Km considerando ambo i sensi di marcia), che in caso di calamità vedrà anche il passaggio dei mezzi di soccorso, è inadeguata per accogliere coloro che saranno in fuga dall’eruzione.

Certo non è detto che tutte le famiglie residenti in queste aree individuate scelgano di usare la SS268 ma è molto probabile che lo facciano visto che è la superstrada più vicina.

E’ attiva una convenzione tra ANAS e Regione Campania per effettuare il raddoppio della SS268, convenzione che dal 2016 al 2020 prevede un investimento di 304 milioni di euro. Dalle slide presenti al link precedente sembrerebbe che oramai i lavori sono agli sgoccioli ma la realtà è ben diversa ed i lavori sono tutt’altro che in fase di ultimazione. Auguriamoci che i lavori di ammodernamento e raddoppio finiscano prima della prossima eruzione del Vesuvio e non durino quanto quelli per l’ultimazione della Salerno-Reggio Calabria!

L’uso del suolo in Campania. Dati dal progetto Corine Land Cover

Guardando al susseguirsi dei dati, dal 2000 al 2012, del progetto Corine Land Cover potrebbe sembrare che tutto sommato il modo in cui viene sfruttato il suolo in Campania (per fini antropici e non) è rimasto invariato. Ma è andando a studiare il terzo livello del DB geografico del progetto che si scoprono grandi e, a mio avviso, preoccupanti verità.

Nella tabella che segue ho raggruppato le singole categorie della Corine Land Cover nei tre censimenti, sommandone le aree; poi ho effettuato una differenza tra le stesse categorie a differenti anni. In rosso sono evidenziate aree in decrescita rispetto al censimento precedente. Un ettaro (ha) è circa un campo da calcio 11 regolamentare.

Come si può notare il “Tessuto urbano continuo” è andato via via decrescendo in maniera abbastanza constante, potrebbe sembrare una cosa positiva ma andando a guardare il “Tessuto urbano discontinuo” si nota come questo sia aumentato in maniera altrettanto costante. Per avere un paragone, in Campania, i comuni di Lacedonia e San Bartolomeo in Galdo hanno una estensione molto simile all’espansione del tessuto urbano discontinuo. Tra il 2006 ed il 2012 la Campania ha “acquisito” quindi altri due comuni in termini di antropizzazione e di consumo del suolo. Non facciamoci ingannare dalla dicitura “Tessuto urbano discontinuo” poichè anch’esso toglie spazio al suolo naturale avendo necessità di infrastrutture stradali, sottoservizi, opere che rendono impermeabile il suolo. Ed infatti si nota un aumento anche delle reti stradali. La domanda che mi pongo è:

A fronte di un aumento delle reti stradali c’è stato un parallelo miglioramento del tessuto stradale preesistente?

Sempre restando nella macro area delle “Aree edificate” si nota un aumento significativo anche dei siti industriali, aree di estrazione mineraria, discariche. Volendo non vedere tutto grigio, fa ben sperare il dato dell’aumento del “Verde urbano“, si spera che lo stesso aumento venga rilevato anche nella prossima versione della Corine Land Cover.

Le perdite nella macro area “Aree agricole” potrebbero essere lette come una lenta ma costante perdita di suolo agricolo e più in generale di quella cura che mettevano i contadini di una volta nel manutenere i campi e le aree non boschive come i pascoli, anch’essi in declino. Non avendo una formazione agraria non so dare altre letture.

Questo slideshow richiede JavaScript.

Nemmeno i boschi se la passano bene in Campania. Quasi tutte le categorie del terzo livello della CLC della macro area “Aree boschive e naturali” sono in rosso e non basta il solo fenomeno degli incendi estivi a spiegare la lenta ma inesorabile scomparsa dei boschi campani. La sommatoria dei minus delle aree boschive si avvicina quasi all’estensione delle aree antropizzate ed in particolare del tessuto urbano discontinuo. Le aree boschive giocano un ruolo importante sia nella regolazione della temperatura delle aree a loro vicine sia nei riguardi degli eventi piovosi in quanto trattengono l’acqua meteorica facendo in modo che sia poca quella che arriva al suolo, evitando così fenomeni di dilavamento degli strati superficiali del terreno o fenomeni peggiori come le frane; inoltre ospitano diverse specie faunistiche. Un’altra riflessione che mi viene in mente quindi è:

Cosa si è fatto in questi anni a livello regionale e non per mitigare i fenomeni di dissesto idrogeologico in una Regione già di per se con grossi problemi di “instabilità”?

In ultimo, la riduzione delle aree paludose interne e dei corsi d’acqua a favore dell’aumento degli specchi d’acqua fa pensare ad un aumento dell’uso industriale/civile della potenza dell’acqua in movimento o comunque ad suo uno stoccaggio in bacini ad uso potabile.

Quali azioni sono state messe in campo dalle Autorità di Bacino e dagli enti gestori e fornitori di servizi affinchè venga ridotto, se non annullato, il quantitativo d’acqua che viene disperso per condotte ammalorate? Flora e fauna locale come sono state tutelate dalla riduzione dell’apporto idrico?

Questa analisi è stata fatta lavorando su dati pubblici e liberamente scaricabili da chiunque. Chi volesse effettuare ulteriori analisi potrà scaricare i dati già rielaborati della CLC dai link che seguono:

I dati sono in formato GeoPackage.

Cartografia elaborata in formato immagine:

Cosa sono le proiezioni cartografiche?

Spesso nelle richieste di supporto che mi arrivano, la problematica riscontrata è dovuta ad una non conoscenza dei concetti alla base dei sistemi di riferimento. Ho notato che uno dei problemi più ricorrenti sono le proiezioni cartografiche. In questo articolo cercherò di fare un po’ chiarezza su cosa sono.

Dal punto di vista geometrico-matematico una proiezione è una “azione”, un insieme di funzioni che sono strutturate in modo da “inviare” ogni punto di una superficie su un’altra. Certo questa che ho scritto è una semplificazione estrema ma, visto che il pubblico che si interfaccia con l’ambiente GIS non è fatto solo di professionisti che hanno sostenuto gli esami di analisi 1 e 2 ed algebra lineare, ce la faremo bastare!

Andando un po’ più nello specifico del mondo GIS una proiezione cartografica è una funzione che ha la finalità di convertire le coordinate geografiche(sferiche) in coordinate cartesiane(metriche), quindi ha il compito di portare su un piano ciò che nella realtà è tridimensionale per meglio rappresentarlo. Una proiezione cartografica se mantiene costante i rapporti tra le superfici si dice equivalente(es. la proiezione sinusoidale), se mantiene costante i rapporti tra le distanze di dice equidistante(es. la proiezione cilindrica equidistante), mentre se mantiene costante gli angoli si dice conforme o equiangola(es. la proiezione conforme di Gauss).

proiezioni cartograficheIl tipo di funzione ci consente di conoscere la tipologia di superficie geometrica su cui proiettare la superficie terreste. Questa potrebbe essere infatti una superficie cilindrica, una superficie conica o una superficie piana.

L’orientamento ci fa capire come interagiscono le superfici tra loro. Nello specifico, se consideriamo un parallelo, la superficie di proiezione può essere normale, si ha quando le due entità sono parallele, trasversa, quando la superficie è ortogonale al parallelo, oppure obliqua, quando c’è una certa inclinazione tra parallelo e superficie.

 

Il tipo di incidenza definisce come le due superfici vengono in contatto. Ipotizziamo che la superficie di proiezione sia un piano, esso può essere tangente oppure secante alla superficie terrestre.

proiezioni cartograficheIl centro di emanazione(P) definisce il tipo di proiezione prospettica. Se esso è coincidente con il centro della superficie terrestre abbiamo una proiezione centrografica o gnomonica(a), se esso è posizionato dalla parte opposta all’intersezione tra superficie terrestre e superficie di proiezione abbiamo una proiezione stereografica(b), se esso è posizionato dalla parte opposta all’intersezione tra superficie terrestre e superficie di proiezione ed è esterno alla superficie terrestre abbiamo una proiezione scenografica(c) mentre se è posizionato sul piano ortogonale all’intersezione tra superficie terrestre e piano di proiezione abbiamo una proiezione ortografica(d).

 

La scelta del tipo di proiezione ci restituirà una rappresentazione della superficie terrestre che avrà delle alterazioni in alcuni punti. Prendiamo il caso della proiezione di Mercatore, nella webmap che segue si può notare come aumentano le distorsioni man mano che ci allontaniamo dall’equatore.

Dopo quanto scritto sarà senz’altro chiaro che non bisogna trattare con leggerezza questo tema e che il passaggio da un sistema proiettivo ad un altro, se fatto senza criterio, può comportare grossi errori.

UPDATE 7 maggio 2018

In QUESTA pagina sono raccolte quasi tutte le proiezioni cartografiche esistenti, dagli uno sguardo!

Idraulica&GIS: quota serbatoio e tracciato condotta di avvicinamento

Nella progettazione di una rete idrica è importante andare ad individuare la posizione più idonea per il/i serbatoio/i a servizio della rete interna ed il tracciato della condotta di avvicinamento.

Questo articolo nasce per puro caso. Qualche giorno fa sistemando la libreria ho ritrovato i miei appunti del corso di Infrastrutture Idrauliche seguito quando ero studente di ingegneria. Insieme a quegli appunti c’era anche la relazione, i calcoli, i grafici ed i disegni del progetto di un acquedotto da discutere all’esame. I primissimi passaggi erano proprio la ricerca del punto migliore in cui posizionare il serbatoio a servizio della rete interna e l’individuazione del tracciato per la condotta di avvicinamento. Di seguito spiegherò come è possibile con QGIS definire queste due variabili progettuali alla base della progettazione di un acquedotto.

NB: l’intento di questo articolo non è quello di essere un tutorial passo passo poichè il suo unico scopo è divulgativo.

All’università ricordo che ci fu data una carta IGM 1:25.000 su cui era individuato il Comune da servire e la sorgente da cui prelevare l’acqua. Con solo questi dati individuammo(parlo al plurale perchè ero in un gruppo di lavoro) la quota del serbatoio ed aiutandoci con le isoipse definimmo il tracciato con il relativo profilo altimetrico. In realtà la parte più scocciante e lunga fu proprio individuare il tracciato e definire il profilo altimetrico perchè andammo ad individuare 3-4 soluzioni. Alla fine scegliemmo quella che ci sembrava meno gravosa dal punto di vista dei costi sia di scavo che di esproprio oltre che per la lunghezza del tracciato in se. Ricordo che per calcolare il profilo altimetrico usammo un applicativo in riga di comando, sembrava il DOS, in cui inserimmo i dati che avevamo rilevato sulla carta IGM e lui alla fine ci disegnò il profilo altimetrico inserendo la distanza progressiva ed assoluta sezione per sezione oltre che la pendenza.

Finito l’excursus romantico andiamo al sodo!

La situazione in avvio di progetto è questa:

  • sorgente a quota 444 metri;
  • quota del serbatoio, calcolata tenendo conto della quota massima del centro abitato da servire, dell’altezza dell’edificio più alto, delle perdite di carico, pari a 430 metri.

In aggiunta a questi dati ed alla carta IGM ho usato il DTM dell’ISPRA con risoluzione di 20 metri per pixel. Ho riclassificato un clip dell’area di studio secondo tre fasce altimetriche:

  • Classe 0: quote comprese tra il minimo altimetrico del clip e quota 427 metri;
  • Classe 1: quote comprese tra 428 e 432 metri;
  • Classe 2: quote comprese tra 433 metri ed il massimo altimetrico del clip.

La classe di mio interesse è la 1 che comprende la quota del serbatoio +/- 2 metri. Sovrapponendo il raster riclassificato con un buffer a 250 metri dal centro abitato sarà possibile individuare l’area migliore per posizionare il serbatoio. I pixel bianchi sono quelli in corrispondenza di aree utili allo scopo.

Individuata l’area migliore per posizionare il serbatoio sono passato alle step successivo: il tracciamento del percorso della condotta di avvicinamento.

Nel video spiego tutti i passaggi che ho effettuato fino all’individuazione del tracciato in se.

Noto il tracciato mi sono servito del plugin qProf per tracciare il profilo altimetrico con distanza progressiva tra cuspidi positive e negative, distanza totale e pendenza dei tratti. E’ importante sapere che il profilo ottenuto è in funzione della risoluzione del DEM da cui ha origine. Se si ha un DEM a 1m/px il risultato sarà sicuramente migliore del mio a 20m/px.

Nel video che segue spiego come ricavare il profilo altimetrico ed esportarlo come csv per essere poi elaborato in un foglio di calcolo.

L’esercizio è finito, alla prossima 🙂

SCARICA QUI I DATI USATI NELL’ARTICOLO

Il GeoPackage: una valida alternativa al formato shape!

In questo articolo vi spiegherò perchè con il GeoPackage potete liberarvi degli shapefile

Il GeoPackage è un formato dati open che fu pubblicato per la prima volta nel 2014. E’ un formato che segue gli standard OGC (Open Geospatial Consortium) ed ha tutta una serie di qualità che lo rendono migliore del vecchio ed oramai superato formato shape.

Prima di conoscere in maniera quanto più approfondita possibile il GeoPackage, andremo a vedere i motivi per cui il formato shape DEVE essere abbandonato. In rete c’è un articolo dal titolo abbastanza significativo “Switch from Shapefile” che vi invito a leggere. In questo articolo sono elencate tutte le problematiche legate all’utilizzo di uno shapefile; di seguito mi limiterò ad elencarne alcune tra quelle più comuni che ho personalmente incontrato:

  • Il nome di un field della tabella attributi non può contenere più di 10 caratteri. Personalmente trovo fastidiosissimo dover troncare i nomi delle colonne! Se devo chiamare un field

    La struttura di uno shapefile

    Superfici_in_evoluzione“, con il formato shp sono costretto a scrivere cose tipo “sup_in_evo“. Brutto da vedere, brutto da leggere, brutto da gestire per la costruzione di una legenda sia in TOC che in layout di stampa.

  • Non si può andare oltre i 255 caratteri per il tipo di dato testuale. Sono davvero tantissimi i casi in cui 255 caratteri sono davvero pochi. Per la costruzione di una carta geologica, ad esempio, quando si ha la necessità di descrivere il tipo di unità geologica, 255 caratteri stanno stretti. Una (non)soluzione al problema che a volte ho adottato in questi casi è creare una tabella esterna con i testi lunghi che ho poi associato allo shp. Questo comporta la necessità di aggiungere un altro file alla già lunga lista di file associata allo shp con l’aumento della possibilità di copia parziale dei file se si devono inviare ad altre persone i dati. Sarà indispensabile inviare anche il file progetto per non perdere l’unione tra tabella esterna ed interna(un altro file da aggiungere!).
  • Non è possibile impostare un set di codifica caratteri. Vi siete mai trovati ad usare dati in una lingua diversa? Non per forza cinese o arabo. Mal di pancia vengono già con l’utilizzo di testi scritti in francese! Caratteri illeggibili, un disastro!
  • Lo shapefile non è un formato dati topologico! Questo significa che le geometrie che vengono generate non sono topologiche e quindi è probabile che incapperete in errori, ad esempio, in geoprocessing. Vi riporto una mia esperienza recente.
    Avevo uno shp pieno di geometrie non valide. Lo importai in PostGIS e con la funzione ST_MakeValid lo resi topologico. Quando lo esportai in shp di nuovo non lo era. Per prova lo esportai da PostGIS in SpatiaLite e risultava corretto anche in SpatiaLite. Come ulteriore prova lo esportai da SpatiaLite in shp ed anche in questo caso mi ritrovai un vettore con geometrie invalide.

Per un professionista del settore GIS queste già dovrebbero essere argomentazioni sufficienti a dire addio allo shapefile; nell’articolo che vi ho linkato ce ne sono molte altre!

Veniamo al GeoPackage. E’ un formato di dati geografici molto duttile, può infatti contenere sia vettori che raster. In realtà questo formato è un data container basato su SQLite. Questo ci consente di avere tanti vantaggi, vediamone alcuni.

ALL IN ONE. Il Geopackage può contenere al suo interno le geometrie, la tabella attributi associata alle geometrie, la topologia delle geometrie, il sistema di riferimento del vettore ed i tematismi!! Il formato shp per fare quasi la stessa cosa ha bisogno di almeno 6 file.

Essendo basato su SQLite possiamo senza problemi scrivere per esteso il titolo dei nostri field. Se il field in questione è di tipo testo non abbiamo nessuna limitazione ai caratteri e quindi possiamo inserire senza difficoltà testi molto lunghi. Essendo un formato che ricalca in pieno gli standard OGC avremmo massima libertà di uso ed interoperabilità tra piattaforme software, inoltre può immagazzinare più di 2GB di dati cosa che uno shapefile non può fare.

Un po’ di mesi fa fu lanciato su Twitter un sondaggio sullo switch dal formato shp, ed il GeoPackage è stato scelto come miglior sostituto surclassando il formato GeoJson.

Un’altra cosa interessante del GeoPackage è il fatto che può contenere raster. Con QGIS è semplicissimo passare da un .tiff ad un .gpkg(l’estensione del GeoPackage), basta effettuare una conversione di formato.

Visti i tanti pro perchè non passate al formato .gpkg?

Fotovoltaico e GIS: come individuare rapidamente l’utilizzabilità di un tetto!

Nella mia precedente vita da commerciale per un’azienda che vende impianti fotovoltaici, mi sono sempre chiesto se e come fosse possibile effettuare un rilievo speditivo su un’area di interesse per individuare quali tetti fossero più idonei per l’installazione di un impianto fotovoltaico e stimarne il potenziale in Kwp.

Oggi, a distanza di un po’ di anni da quell’esperienza, ho finalmente trovato un po’ di tempo per scrivere questo articolo in cui, con la metodologia GIS, spero di dare una buona mano ai miei ex colleghi commerciali 😉

Quello di cui abbiamo bisogno è:

  • un vettore poligonale che riproduce l’edificato dell’area di interesse;
  • un DSM da LiDAR o comunque con risoluzione non superiore ad 1 m;
  • competenze in spatial analisys.

Iniziamo effettuando un clip del nostro raster usando, come poligono di clip, l’edificato. Tematizzando il DEM appena ottenuto dovreste avere qualcosa di simile a quello che vedete di seguito. La scelta dei colori è arbitraria; al variare dell’intensità del blu aumenta l’altezza dell’edificio.

fotovoltaico gis

Dando un’occhiata alle pendenze potremmo già farci una prima idea su quale tetto è idoneo e quale non lo è. Le pendenze le estraiamo con slope da Processing Toolbox di QGIS; anche in questo caso la tematizzazione è a vostro gusto.

fotovoltaico gis

NB: le pendenze fanno riferimento al piano cartografico! Quindi 0° fanno riferimento ad un tetto piano, con 85.8° di pendenza abbiamo un tetto quasi verticale(e ortogonale al piano cartografico)!

Sicuramente un tetto a quasi 90°, alle nostre latitudini, non è il più idoneo all’installazione di un impianto fotovoltaico!

Iniziamo ad estrarre i soli tetti piani. Per farlo dobbiamo eseguire queste tre operazioni in sequenza:

  1. Riclassificazione del raster delle pendenze;
  2. Generazione di un vettore dal raster appena riclassificato;
  3. Selezione e salvataggio dei soli poligoni che rientrano nella classe di interesse.

Per riclassificare il nostro raster delle pendenze useremo il geoalgoritmo r.reclass di GRASS GIS ricercabile dal Processing Toolbox. Importante è la corretta definizione delle classi!

Ho suddiviso le pendenze in due classi: maggiore di 3° e minore di 3°. Tutti i tetti con pendenza massima di 3° li ho considerati piani ed è questa la classe che ci interessa in questa fase. Al termine della riclassificazione ho convertito in poligoni le classi raster usando Polygonize sempre dal Processing Toolbox. Sul vettore poligonale così ottenuto ho effettuato una query per estrarmi i soli tetti che rientrano nella categoria di nostro interesse.

Ecco il risultato! Fatelo anche voi, sembra difficile ma non lo è!

fotovoltaico gis

Torniamo ora a ragionare sul nostro DEM clippato ad inizio articolo. Abbiamo le pendenze si ma a noi interessa anche l’esposizione. Un tetto inclinato ed esposto a SUD è la soluzione che più auspichiamo ma come lo individuiamo? Ora andremo a vedere come fare!

Usiamo, sul DEM clippato ad inizio articolo, il geoalgoritmo Aspect ricercabile da Processing Toolbox.

fotovoltaico gis

Riclassifichiamo il risultato del precedente geoalgoritmo in modo da ottenere le classi di esposizione di nostro interesse.

NB: vi consiglio di leggere attentamente il contenuto nel link che ho inserito precedentemente in Aspect!

Per semplicità e rapidità ho creato tre classi di esposizione: NORD-EST, EST-SUD-OVEST e OVEST-NORD. Quella che interessa a noi è EST-SUD-OVEST. Convertiamo in poligoni il raster ed effettuando una selezione dei soli poligoni di interesse, andremo a salvare la selezione così come fatto per l’individuazione dei tetti piani.

Il risultato è quello che segue e c’è anche una piccola chicca che avvalora il metodo usato 😉

fotovoltaico gis

Verificate la curiosità qui!

Effettuiamo un Union tra i due vettori appena creati(post selezione!) ed eliminiamo eventuali anomalie come tetti che possono essere contemporaneamente piani ed inclinati. Queste anomalie possono esserci perchè le operazione che abbiamo effettuato dipendono dalla risoluzione del DSM su cui abbiamo operato fin dall’inizio. Più andiamo al di sotto del metro di risoluzione meno errori avremmo!

Ora attiviamo il Field Calculator, nella tabella attributi del vettore di output dell’Union, e facciamo pulizia di tutte le aree inferiori a 7 mq. Ho utilizzato questo parametro perchè indicativamente 1 Kwp copre 7 mq. Ovviamente per fare questa pulizia dobbiamo calcolare preventivamente le aree dei singoli poligoni!

fotovoltaico gis

Siamo alla fine della nostra ricerca! Nell’immagine precedente possiamo vedere il vettore dell’edificato e quello dei tetti potenzialmente utili. Man mano che si intensifica il viola abbiamo una maggiore produzione di elettricità. Quindi, dei tanti edifici presenti in questa area di 25 Ha, possiamo vedere che solo 3-4 hanno un elevata produzione di energia (siamo tra i 70 ed i 90 Kwp) mentre gli altri hanno una produzione bassa.

Anche la mappa di calore ci evidenzia dove sono i potenziali target!

fotovoltaico gisE’ possibile affinare ancora di più la ricerca? Si! Ad esempio facendo un passo indietro rispetto all’Aspect ed estraendo delle pendenze di tetto utili scartando le altre, per poi determinare su queste ultime l’esposizione.

Provate voi! 😉

Eliminare i duplicati da un vettore

Lavorando sui dati di OSM mi sono trovato a dover affrontare un piccolo grande problema di sovrapposizione di vettori lineari.

Ho scaricato una grossa mole di dati relativi all’idrografia della Francia orientale usando il plugin QuickOSM e, siccome l’area era davvero vasta, l’ho dovuta suddividere in quadranti di 100km di lato. Solo in questo modo, sfruttando ogni singolo quadrante, ho potuto effettuare il download in maniera corretta e senza superare il tempo massimo di risposta del server.

Questa procedura ha però fatto nascere la problematica che mi ha portato a fare questo articolo ed il video tutorial che troverai in fondo. In pratica, ogni qualvolta scaricavo da un quadrante l’idrografia ed esso era attraversato da fiumi che a loro volta attraversavano altri quadranti, questi venivano caricati anche nel quadrante attivo in quel momenti e, visto che tutto confluiva in un unico vettore lineare presente in un GeoDB contenuto in PostGIS, mi sono ritrovato tantissimi duplicati; nello specifico 11.919! I duplicati oltre a rendere il vettore pesante lo rendevano affetto da errori topologici.

Mi sono accorto di questi duplicati perchè ho fatto una verifica topologica con il plugin Topology Checker impostando la verifica sui duplicati, così come visualizzato nell’immagine che segue.

duplicati

Effettuata la verifica il risultato è stato quello che vedi di seguito.

duplicati

Indagando la tabella attributi ho notato la presenza della colonna full_id e, guardando il feature ID delle geometrie duplicate presenti nel report di Topology Checker, ho notato che in quella colonna i duplicati erano facilmente individuabili perchè avevano lo stesso full_id.

duplicati

Istintivamente mi sono così concentrato su quella colonna e visto che volevo risolvere il tutto usando l’SQL mi sono documentato in giro sul web su come poter eliminare il problema.

Ho prima lanciato questa query:

SELECT
full_id,
COUNT(full_id) AS counter
FROM waterway_new
GROUP BY full_id
ORDER BY counter;

Individuando così un elenco di duplicati. Alcuni fiumi erano duplicati anche 7 volte!

Poi ho usato la query che segue per creare un nuovo vettore depurato dai duplicati:

CREATE TABLE waterway_new_noduplicate AS
SELECT DISTINCT ON (full_id) *
FROM waterway_new;

Tutta la procedura l’ho racchiusa nel video tutorial che trovi di seguito. Un’altra strada sarebbe potuta essere quella di usare il plugin MMQGIS senza quindi usare l’SQL.

 

Corso di Webmapping GRATUITO!

Descrizione del corso

Le webmap, mappe interattive, si stanno diffondendo come strumenti di informazione semplici ma di grande effetto. Possiamo vederle anche sui siti delle testate giornalistiche ad esempio o come strumento di promozione territoriale e più in generale sono uno strumento di rapida divulgazione di informazioni di qualunque genere, rapidità garantita anche dalla loro utilizzabilità su smartphone e tablet. Inoltre imparare i rudimenti sulla crazione di webmap è alla base dello sviluppo di sistemi WebGIS.


Esempio di webmap

Il corso è mirato ad impartire le conoscenze basilari per la produzione di webmap. Al termine del corso i partecipanti saranno in grado di sviluppare in proprio una webmap. Si partirà dal foglio bianco spiegando il minimo indispensabile di HTML, CSS e JavaScript per poter creare semplici webmap ma di grande effetto.

Il corso dura 4 ore, è suddiviso in 4 moduli più una esercitazione di riepilogo:

  • MODULO_1 Dal foglio bianco alla WebMap. In questo modulo si tratteranno le sintassi minime di HTML, CSS e JavaScript, necessarie al caricamento del tile layer in una WebMap basata sulle librerie Leafletjs;
  • MODULO_2 Choropleth map. In questo modulo si impartiranno le informazioni necessarie alla realizzazione di una mappa coropletica; si insegnerà quindi ad inserire un vettoriale di tipo poligonale in una WebMap basata sul tile layer visto nel modulo precedente;
  • MODULO_3 Mappa degli infopoint. Dopo aver visto, nel modulo precedente, come inserire e tematizzare un vettore poligonale, con questo modulo si andrà ad inserire in una WebMap un vettori di tipo puntuale;
  • MODULO_4 Mappa dei percorsi turistici. In questo modulo si vedrà come inserire un vettore di tipo lineare in una WebMap andando a tematizzare i percorsi in funzione del dato contenuto in ogni singola feauture lineare;
  • Esercitazione: Crea la tua prima webmap!

Temi trattati

  • Elementi di HTML, CSS e JavaScript
  • Libreria LeafletJS
  • Plugin qgis2web

Prerequisiti

Conoscenze elementari sui GIS, familiarità minima con QGIS.

Requisiti software

GUARDA LE VIDEO LEZIONI

La playlist comprende tutte le lezioni del corso, puoi visionarla qui o direttamente su Youtube, puoi visionare la singola lezione o solo quella che ti interessa.

Clicca qui e scarica le WebMap realizzate nei moduli

Analisi dei dati Copernicus relativi agli incendi sul Vesuvio di luglio 2017

In luglio un incendio doloso di vastissime proporzioni ha interessato l’area del Parco Nazionale del Vesuvio generando ingenti danni. Con questo articolo, in cui ho rielaborato i dati Copernicus, ho provato a quantificare l’estensione delle tipologie di suolo distrutte avvalendomi anche dei dati Corine Land Cover e del DTM LiDAR della Città Metropolitana di Napoli.

In un mio allenamento in bici di quei giorni mi sono avvicinato al Monte Somma e lo scenario era quello che vedete nel video che segue:

Il Parco ha bruciato per quasi una settimana sia lato Monte Somma che Vesuvio.

L’immagine che segue mostra l’uso del suolo prima degli incendi, i dati sono una estrapolazione della CLC2000 presente nel vector package scaricabile da qui.

vesuvio

L’estensione totale dei vettori che compongono i diversi tipi di uso di suolo è pari a circa 4.400 ha, risulta essere più piccola rispetto ai dati di Copernicus solo perchè ho preferito focalizzarmi sulle aree di maggior danno.

Descrizione suolo CLC2000 Estensione tipo
di suolo(ha)
Estensione tipo
di suolo(%)
Aree a vegetazione boschiva ed arbustiva in evoluzione 932,96 21
Aree con vegetazione rada 89,8 2
Aree prevalentemente occupate da colture agrarie con presenze di spazi naturali importanti 255,31 6
Boschi di conifere 677,04 15
Boschi di latifoglie 811,5 18
Boschi misti conifere e latifogli 732,78 16
Frutteti e frutti minori 59,6 1
Rocce nude, falesie, rupi, affioramenti 257,17 6
Sistemi colturali e particellari complessi 653,18 15
Estensione totale 4.469,34

 

Le estensioni maggiori riguardano le aree boschive e purtroppo sono loro ad aver avuto i maggiori danni. L’immagine che segue rappresenta lo scenario post incendi.

vesuvio

Descrizione suolo CLC2000 Estensione
tipo di suolo
pre incendi(ha)
Aree completamente
distrutte (ha)
Aree fortemente
danneggiate (ha)
Aree debolmente
danneggiate (ha)
Estensione
tipo di suolo
post incendi(ha)
Aree a vegetazione boschiva ed arbustiva in evoluzione 932,96 447,29 114,31 47,46 323,90
Aree con vegetazione rada 89,8 15,88 7,72 3,93 62,7
Aree prevalentemente occupate da colture agrarie con presenze di spazi naturali importanti 255,31 7,78 1,81 0,65 245,07
Boschi di conifere 677,04 334,90 88,02 50,68 203,44
Boschi di latifoglie 811,5 50,09 30,25 47,90 683,26
Boschi misti conifere e latifogli 732,78 108,42 35,94 18,89 569,44
Frutteti e frutti minori 59,6 9,40 5,37 8,11 36,72
Rocce nude, falesie, rupi, affioramenti 257,17 38,62 7,96 1,49 209,10
Sistemi colturali e particellari complessi 653,18 32,79 11,84 5,63 602,92
Estensione totale 4.469,34 1.045,17 303,22 184,83 2.936,12

Analizzando i dati nella tabella precedente è facile valutare che è andato in fumo:

  • il 70% dei boschi di conifere, passando da 677,04 ha a 203, 44 ha;
  • il 65% della aree a vegetazione boschiva ed arbustiva in evoluzione, passando da 932,96 ha a 323,90 ha;
  • il 38% delle aree a frutteti e frutti minori, passando da 59,6 ha a 36,72 ha.

Questo solo per citare le prime tre tipologie di suole per danni subiti. Le aree totalmente distrutte, ricadenti nel Parco Nazionale del Vesuvio, ricoprono una superficie di 1.045,17 ha; per fare un confronto, un campo da calcio a 11 regolamentare è circa 1 ha.

Di seguito sono presenti alcuni diagrammi di confronto.

vesuvio vesuvio vesuvio

vesuvio

Dati riferiti all’estensione dell’area di studio

vesuvio

Dati riferiti all’estensione dell’area di studio

 

Cosa dobbiamo aspettarci con le prossime piogge, tipicamente abbondanti soprattutto tra settembre ed ottobre? Staremo a vedere, purtroppo.

Intanto mi rincuora il fatto che la natura piano piano si sta riprendendo i suoi spazi; per ripristinare i boschi però ci vorranno anni se non secoli!

Di seguito è possibile scaricare:


Ricostruzione in 3D dell’area interessata dagli incendi

copernicus

Clicca sull’immagine per andare alle vista 3D

Software usati per l’elaborazione:

Si ringraziano i seguenti enti: