Per lavoro utilizzo abbastanza frequentemente alcuni script (procedure automatizzate e personalizzate) per velocizzare le attività che svolgo quotidianamente, i quali si rivelano molto utili quando il tempo a disposizione è veramente poco. Tra questi ce n’è uno in particolare che mi piacerebbe condividere e riguarda uno script per scaricare il DTM del Geoportale Nazionale.
Il Geoportale Nazionale mette a disposizione dei servizi web (i famosi servizi OGC) che consentono un accesso libero e completo ai dati geografici disponibili nel catalogo. Tra questi esistono due servizi in particolare che permettono il download di dataset di interesse, il WCS (Web Coverage Service) per i dati di tipo coverage ed il WFS (Web Feature Service) per i dati vettoriali.
Ormai, quasi tutti i software GIS permettono facilmente di caricare e visualizzare queste tipologie di servizi, fungendo da client per gli stessi servizi, i quali solitamente vengono forniti attraverso un indirizzo URL. E’ banale, per esempio con QGIS, accedere a tali servizi ed esistono centinaia, se non più, di tutorial su internet che illustrano passo dopo passo come eseguire una simile procedura.
In questo breve articolo vedremo come scaricare il DTM con risoluzione a 20 metri accessibile tramite il servizio WCS. La descrizione del dataset può essere consultata dall’XML prodotto nella richiesta di GetCapabilities:
Il lavoro e’ partito agli inizi degli anni 80. Il modello numerico del terreno e’ stato ottenuto interpolando le curve di livello dell’archivio IGM. Oltre a questi sono state utilizzate le linee di costa dei laghi, aventi quota pari al livello delle acque dei laghi stessi, e del mare (quest’ultima proveniente da cartografia scala 1:100.000) per arrestarne l’interpolazione. La quota, espressa in metri, e’ il valore altimetrico medio incontrato nell’area della cella elementare (pixel) riferito al baricentro della cella stessa. Le curve di livello sono state digitalizzate con varie metodologie quindi l’accuratezza del dato dipende dal tipo di terreno (pianura, collina, montagne). L’accuratezza planimetrica mediamente e’ di circa 3-4 decimi di millimetro grafico.
Un modo alternativo per scaricare il DTM, che a me è tornato utile, è quello illustrato di seguito in questo articolo, facendo uso di uno script python da eseguire comodamente nella Python Console di QGIS.
Come funziona?
Il funzionamento è veramente banale, prima di tutto bisogna copiare ed incollare lo script nella console python di QGIS e salvarlo successivamente nel filesystem locale. Una volta salvato possiamo eseguire lo script. Per scaricare il DTM è sufficiente trascinare il mouse sull’area della mappa (Map Canvas) per definire l’estensione spaziale dell’area di interesse (boundary box). That’s all!
…affinché il dato venga scaricato e caricato correttamente è necessario impostare il sistema di riferimento WGS84 (EPSG:4326) per il progetto QGIS.
Di seguito una GIF esemplificativa ed il codice dello script:
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| from processing.gui.RectangleMapTool import RectangleMapTool | |
| from owslib.wcs import WebCoverageService | |
| import tempfile | |
| import os | |
| ## risoluzione DEM, i possibili valori sono: 20, 40 e 75 | |
| RES = '20' | |
| def run(): | |
| ms = canvas.mapSettings() | |
| crs = ms.destinationCrs() | |
| msg_bar = iface.messageBar() | |
| r = tool.rectangle() | |
| s = (r.xMinimum(), r.xMaximum(), r.yMinimum(), r.yMaximum()) | |
| try: | |
| if crs.authid() != 'EPSG:4326': | |
| msg_bar.pushMessage("downdem script", | |
| "Il sistema di riferimento del progetto deve essere WGS84.", | |
| 0, 5) | |
| return | |
| dem = get_dem_from_extent(s) | |
| lyr = QgsRasterLayer(dem, "DEM") | |
| QgsMapLayerRegistry.instance().addMapLayers([lyr]) | |
| msg_bar.pushMessage("downdem script", | |
| "DEM caricato correttamente.", | |
| 3, 5) | |
| finally: | |
| tool.reset() | |
| canvas.setMapTool(map_tool) | |
| def get_dem_from_extent(s): | |
| dirpath = tempfile.mkdtemp() | |
| bbox = (s[0], s[3], s[1], s[2]) | |
| wcs_url = 'http://wms.pcn.minambiente.it/wcs/dtm_' + RES + 'm' | |
| w = WebCoverageService(wcs_url, version='1.0.0') | |
| output = w.getCoverage(identifier='EL.DTM.' + RES + 'M', | |
| format='GeoTIFF', | |
| bbox=bbox, | |
| crs='EPSG:4326', | |
| width=400, height=300) | |
| with open(os.path.join(dirpath, 'dem.tif'), 'wb') as f: | |
| f.write(output.read()) | |
| return os.path.join(dirpath, 'dem.tif') | |
| canvas = iface.mapCanvas() | |
| map_tool = canvas.mapTool() | |
| tool = RectangleMapTool(canvas) | |
| tool.rectangleCreated.connect(run) | |
| canvas.setMapTool(tool) |
.salvatorelarosa 
Grazie per la condivisione.
sono un neofita del linguaggio Python e mi piacerebbe sapere se, con lo stesso script (con opportune modifiche), fosse possibile scaricare anche il dtm40 e dtm80; basta cambiare solo la riga 36? cioè, è sufficiente questo http://wms.pcn.minambiente.it/wcs/dtm_40m per scaricare il dtm40?
grazie!!!
LikeLike
Ciao Totò, ho modificato lo script, adesso puoi scegliere la risoluzione del DTM, i valori possibili, da inserire nella variabile RES (riga 7), sono: 20, 40 e 75. Grazie per il suggerimento!
LikeLiked by 1 person
Grazie per la rapida risposta e scusami se continuo…
ho sostituito nella variabile RES sia 40 che 75, l’output è sempre di 20 m/px!!!
dove sbaglio?
LikeLike
Totò, strano. Ho provato e sembra funzioni. La variabile RES serve a modificare la URL del WCS ed il nome del Layer da recuperare (riga 39 e 42 rispettivamente). Sinceramente al momento non mi viene in mente niente, l’unica cosa si potrebbe fare un test direttamente dal client WCS con i tre servizi. Provo a capire dove sta il problema.
LikeLiked by 1 person
Scaricando uno stralcio del DTM 20m del Geoportale, zona Italia Centrale, ho visto che dà pixel 0.000845° che non corrispondono a 20m, come mai?
LikeLike
Ciao Italo, è evidente che qualcosa non torna. Quello che faccio non è altro che esporre la URL del WCS all’interno dello script. La URL è quella disponibile tra i servizi OGC del sito del fornitore. Ho utilizzato la URL per altre cose e devo dire che ho riscontrato anche io qualche anomalia. Investigando più a fondo sul problema mi sono accorto che il Geoportale Nazionale oltre alla URL presente nello script, alla quale si accede dalla pagina dei servizi di download, fornisce altre URL per il servizio medesimo. Basta controllare nel catalogo CSW dove viene riportato un indirizzo differente. Ho provato anche con quest’ultimo ma i risultati sembrano identici. Può essere comunque che entrambe le URL puntino alla stessa risorsa trattandosi semplicemente solo di differenti versioni del servizio WCS. Sinceramente ho provato a contattare il GN direttamente tempo fa ma ad oggi non ho ricevuto risposta. C’è da dire inoltre che lo script utilizza la libreria OWSLib per accedere al servizio e non mi sento di escludere che ci potrebbe essere qualche problema nella libreria.
Mi dispiace se qualcosa non va, ma noi sviluppatori solitamente abbiamo un difetto molto brutto che è la pigrizia :), per cui è veramente facile che qualcosa ci possa sfuggire, per questo confidiamo negli utilizzatori. Allo stesso tempo sono contento che in molti stanno testando anche perché tutto questa rete sociale che si forma in queste occasioni dimostra che quando le cose vengono utilizzate e testate si riesce bene o male a dare un segno costruttivo, sia a fornitori che utilizzatori, per verificare meglio e capire eventualmente dove sta il problema con la preghiera che qualcuno riesce a risolverlo, sempre se quel qualcuno a voglia e tempo per farlo. Grazie.
LikeLike