Script x DTM5 Lombardia

Script per scaricare/ricostruire i dati DTM 5 m della Lombardia

Il presente script in linguaggio Python consente di interrogare automaticamente il servizion WMS relativo DTM 5x5m della regione Lombardia, estrarre la quota e ricostruire il modello digitale di terreno come file ASCII (praticamente un file di testo). Lo script scarica e crea il DTM con sistema di riferimento UTM 32 - WGS 84 (EPSG:32632).
Il file di testo creato può essere caricato in QGIS che lo vede come un normale raster.

L'utente deve cambiare le prime righe dello script per specificare:
- percorso e nome del file di testo e del relativo file prj che verranno creati,
- coordinate dell'angolo in basso a sinistra dell'area di proprio interesse (inserire multipli di 10),
- estensione in pixel (righe e colonne) dell'area di proprio interesse (1 pixel = 5x5 m);
poi è sufficiente copiare lo script, incollarlo nella console python di QGIS e premere due volte invio.
(Alla prima esecuzione dello script è meglio fare una prova con estensioni limitate, es. 10x10 pixel, per valutare la velocità.)

Questo è lo script:

# percorso e nome del file di testo da creare:
filepathname = "D:/Dati/Italia/prova/dtm1.txt"
fileprj = "D:/Dati/Italia/prova/dtm1.prj"

# coordinate dell'angolo in basso a sinistra (inserire multipli di 10):
x1 = 536500
y1 = 5078000

# numero di righe (height) e colonne (width) dell'area in pixel:
hgt = 10
wdt = 10

# risoluzione
res = 5

# valore da assegnare in assenza di dato
valNoData = 0

# parti fisse della richiesta GetFeatureInfo
request1 = 'http://www.cartografia.servizirl.it/arcgis/services/wms/DTM5_RL_wms/MapServer/WmsServer?service=wms&version=1.3.0&request=GetFeatureInfo&query_layers=DTM%205X5&info_format=text/xml&crs=EPSG:32632&bbox='

# importazione librerie
import xml.etree.ElementTree as ET
import urllib2
import time

# calcolo bbox
x2 = x1 + (wdt * res)
y2 = y1 + (hgt * res)
bbox = str(x1) + ',' + str(y1) + ',' + str(x2) + ',' + str(y2)
request2 = request1 + bbox + '&width=' + str(wdt) + '&height=' + str(hgt)

# verifica che l'area indicata sia interna all'estensione del DTM
if x1 > 692666 or x2 < 460001 or y1 > 5166901 or y2 <4944901:
   print 'Area esterna alla copertura del DTM'
else:
   #
   print 'Ok, script avviato'
   start = time.time()
   #
   # crea e inizia il file di output
   out_file = open(filepathname ,"w")
   out_file.write( 'ncols '+ str(wdt) + '\n' )
   out_file.write( 'nrows '+ str(hgt) + '\n' )
   xll = x1 - 3.1306
   yll = y1 + 1.2481
   out_file.write( 'xllcorner ' + str(xll) + '\n' )
   out_file.write( 'yllcorner ' + str(yll) + '\n' )
   out_file.write( 'cellsize ' + str(res) + '\n' )
   out_file.write( 'NODATA_value ' + str( valNoData ) + '\n' )
   #
   # azzeramento contatori valori nulli
   NoDataCount = 0
   NoDataCountTemp = 0
   NoRisp = 0
   NoRispTemp = 0
   NoPixelValue = 0
   NoPixelValueTemp = 0
   #
   # interrogazione del servizio
   #
   iy = 0
   while iy < hgt:
      ix = 0
      row = ''
      while ix < wdt:
         value = "NoData"
         iNoData = 0
         while value == "NoData" and iNoData < 10:
            iNoData = iNoData + 1
            #
            # interrogazione GetFeatureInfo del servizio wms
            request = request2 + '&x=' +str(ix) + '&y=' + str(iy)
            risposta = urllib2.urlopen(request)
            rispostaxml = risposta.read()
            risposta.close()
            #
            # estrazione degli attributi dalla risposta xml
            FeatureInfoResponse = ET.fromstring(rispostaxml)
            #
            # se la risposta è vuota viene incrementato il contatore relativo
            if len( FeatureInfoResponse ) == 0:
               NoPixelValueTemp = 0
               NoRispTemp = NoRispTemp + 1
            # altrimenti viene cercato l'attributo PixelValue
            else:
               NoRispTemp = 0
               fieldxml = FeatureInfoResponse[0]
               if 'PixelValue' in fieldxml.attrib:
                  NoPixelValueTemp = 0
                  value = fieldxml.attrib['PixelValue']
                  #
                  # se la risposta è NoData viene incrementato il contatore relativo
                  if value == "NoData":
                     NoDataCountTemp = NoDataCountTemp + 1
               else:
                  # se non è presente l'attributo PixelValue viene incrementato il contatore relativo
                  NoPixelValueTemp = NoPixelValueTemp +1
            #
         # conteggio dei NoData, risposte assenti, risposte senza PixelValue
         if value == "NoData":
            NoDataCount = NoDataCount + 1
            value = str( valNoData )
         if NoRispTemp > 0:
            NoRisp = NoRisp + 1
            # Togliere il cancelletto dalle righe sottostanti per visualizzare un avviso ad ogni risposta vuota definitiva
            #print 'Attenzione! risposta vuota al punto x=' + str(ix) + 'y=' + str(iy)
            #print rispostaxml
            NoRispTemp = 0
         if NoPixelValueTemp > 0:
            NoPixelValue = NoPixelValue + 1
            # Togliere il cancelletto dalle righe sottostanti per visualizzare un avviso ad ogni risposta definitiva senza PixelValue
            #print 'Attenzione! Pixel Value assente al punto x=' + str(ix) + 'y=' + str(iy)
            #print rispostaxml
            NoPixelValueTemp = 0
         #
         row = row + ' ' + value
         ix = ix +1
      #
      # scrittura della riga nel file di output
      out_file.write( row + '\n' )
      iy = iy + 1
      end = time.time()
      print 'riga ' + str(iy) + '; tempo impiegato: ' + str( end - start ) + ' secondi'
   #
   # chiusura del file di output
   out_file.close()
   #
   # creazione e scrittura del file prj per il sistema di riferimento UTM32-WGS84
   prj_file = open(fileprj ,"w")
   prj_file.write( 'PROJCS["WGS_1984_UTM_Zone_32N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]' )
   prj_file.close()
   #
   print ''
   print 'Creato file ' + filepathname
   print 'Valori NoData temporanei: ' + str(NoDataCountTemp)
   print 'Valori NoData definitivi: ' + str(NoDataCount)
   print 'Risposte vuote (probabilmente al di fuori della copertura geografica): ' + str(NoRisp)
   print 'Risposte senza PixelValue: ' + str(NoPixelValue)
   #

DESCRIZIONE

Parametri di base:

questi parametri possono essere cambiati dall'utente.
Le coordinate sono riferite al sistema UTM 32 - WGS 84 ed indicano l'angolo in basso a sinistra dell'area da scaricare/ricostruire. Lo script allinea automaticamente il DTM ricreato a quello effettivo presente sul server aggiungendo 1,2481 m in latitudine e sottraendo 3,1306 m in longitudine ma perché sia corretto le coordinate inserite devono essere multipli di 10 (ovvero numeri interi terminanti con 0).
Come numeri di righe e colonne è consigliabile iniziare con valori bassi (es. 10) per valutare il tempo necessario.
La risoluzione è impostata di default su 5 m per avere il massimo dettaglio.
valNoData indica il valore che viene assegnato alle celle se l'interrogazione non h arisposta o se la risposta contiene un valore nullo; poichè le quote in Lombardia sono tutte al di sopra del livello del mare si può usare un valore pari a 0.

Indirizzo internet del servizio wms e funzione GetFeatureInfo:

# parti fisse della richiesta GetFeatureInfo
request1 = 'http://www.cartografia.servizirl.it/arcgis/services/wms/DTM5_RL_wms/MapServer/WmsServer?service=wms&version=1.3.0&request=GetFeatureInfo&query_layers=DTM%205X5&info_format=text/xml&crs=EPSG:32632&bbox='

Il DTM viene interrogato grazie alla funzione GetFeatureInfo del servizio WMS (> DTM 5m Lombardia - interrogazione).
L'indirizzo di base è http://www.cartografia.servizirl.it/arcgis/services/wms/DTM5_RL_wms/MapServer/WmsServer al quale sono aggiunti i seguenti parametri:
- service=wms
- version=1.3.0
- request=GetFeatureInfo
- query_layers=DTM%205X5   per interrogare il livello DTM 5X5 (lo spazio è sostituito da %20)
- info_format=text/xml   per avere la risposta in formato xml da cui si possono estrarre i valori chiave
- crs=EPSG:32632
i parametri bbox, width, height, x, y vengono aggiunti automaticamente durante l'esecuzione dello script.

   Librerie:

# importazione librerie
import xml.etree.ElementTree as ET
import urllib2
import time

Per il funzionamento dello script sono necessarie queste librerie:
- xml.etree.ElementTree   per leggere la risposta in formato xml ed estrarre i valori
- urllib2 per accedere agli indirizi internet
- time   per misurare il tempo di esecuzione dello script

   Area:

# calcolo bbox
x2 = x1 + (wdt * res)
y2 = y1 + (hgt * res)
bbox = str(x1) + ',' + str(y1) + ',' + str(x2) + ',' + str(y2)
request2 = request1 + bbox + '&width=' + str(wdt) + '&height=' + str(hgt)

Vengono calcolate le coordinate dei limiti destro e superiore dell'area di interesse (bbox) e all'indirizzo di interrogazione sono aggiunti i parametri bbox, width e heigth.

Verifica che l'area indicata (bbox) sia coperta dal DTM:

...
# verifica che l'area indicata sia interna all'estensione del DTM
if x1 > 692666 or x2 < 460001 or y1 > 5166901 or y2 <4944901:
   print 'Area esterna alla copertura del DTM'
else:
   #
   print 'Ok, script avviato'
   start = time.time()
   ...

Lo script procede solo se l'area indicata ricade almeno parzialmente nell'estensione del DTM. In tal caso compare il messaggio 'Ok, script avviato' e viene registrato l'orario di inizio (start) con time.time().

Creazione e scrittura dei file di testo:

# percorso e nome del file di testo da creare:
filepathname = "D:/Dati/Italia/prova/dtm1.txt"
fileprj = "D:/Dati/Italia/prova/dtm1.prj"
...
   ...
   # crea e inizia il file di output
   out_file = open(filepathname ,"w")
   out_file.write( 'ncols '+ str(wdt) + '\n' )
   out_file.write( 'nrows '+ str(hgt) + '\n' )
   xll = x1 - 3.1306
   yll = y1 + 1.2481
   out_file.write( 'xllcorner ' + str(x1) + '\n' )
   out_file.write( 'yllcorner ' + str(y1) + '\n' )
   out_file.write( 'cellsize ' + str(res) + '\n' )
   out_file.write( 'NODATA_value ' + str( valNoData ) + '\n' )
   ...
      ...
      # scrittura della riga nel file di output
      out_file.write( row + '\n' )
      ...
   ...
   # chiusura del file di output
   out_file.close()
   #
   # creazione e scrittura del file prj per il sistema di riferimento UTM32-WGS84
   prj_file = open(fileprj ,"w")
   prj_file.write( 'PROJCS["WGS_1984_UTM_Zone_32N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]' )
   prj_file.close()
   #

I file di testo sono creati con il comando variabile_file = open( percorso/nome, "w") dove "w" indica la possibilità di scrittura.
Il testo viene inserito con variabile_file.write( testo ). \n indica il passaggio alla riga sottostante (a capo). Alla fine il file deve essere chiuso con .close().
In out_file vengono prima inseriti i parametri di base del file:
   ncols
   nrows
   xllcorner
   yllcorner
   cellsize
   NODATA_value
(Per allineare il DTM ricreato con quello effettivo lo script aggiunge 1,2481 m in latitudine e sottrae 3,1306 m in longitudine ma ciò risulta solo corretto solo se le coordinate inserite sono multipli di 10, ovvero numeri interi terminanti con 0).
Poi vengono aggiunti i valori delle celle di ogni riga una volta che tutta la riga è stata calcolata. I valori vanno da sinistra verso dastra e le righe dall'alto verso il basso.
Quindi risulta un file di testo tipo:

ncols 10
nrows 10
xllcorner 536496.8694
yllcorner 5078001.2481
cellsize 5
NODATA_value 0
1500.668579 1504.863159 1509.191895 1512.725342 1516.002197 ...
1502.546143 1506.083008 1509.546387 1513.095703 1516.451904 ...
1505.094482 1508.475220 1511.342407 1514.268188 1517.391846 ...
1507.908691 1510.892212 1513.254639 1515.883057 1518.879883 ...
...

Contatori:
contano le volte in cui non si riceve un valore di quota valido dal servizio WMS.

...
   # azzeramento contatori valori nulli
   NoDataCount = 0
   NoDataCountTemp = 0
   NoRisp = 0
   NoRispTemp = 0
   NoPixelValue = 0
   NoPixelValueTemp = 0
   ...
      ...
         value = "NoData"
         iNoData = 0
         while value == "NoData" and iNoData < 10:
            iNoData = iNoData + 1
            ...
            # se la risposta è vuota viene incrementato il contatore relativo
            if len( FeatureInfoResponse ) == 0:
               NoPixelValueTemp = 0
               NoRispTemp = NoRispTemp + 1
            # altrimenti viene cercato l'attributo PixelValue
            else:
               NoRispTemp = 0
               fieldxml = FeatureInfoResponse[0]
               if 'PixelValue' in fieldxml.attrib:
                  NoPixelValueTemp = 0
                  value = fieldxml.attrib['PixelValue']
                  #
                  # se la risposta è NoData viene incrementato il contatore relativo
                  if value == "NoData":
                     NoDataCountTemp = NoDataCountTemp + 1
               else:
                  # se non è presente l'attributo PixelValue viene incrementato il contatore relativo
                  NoPixelValueTemp = NoPixelValueTemp +1
            #
         # conteggio dei NoData, risposte assenti, risposte senza PixelValue
         if value == "NoData":
            NoDataCount = NoDataCount + 1
            value = str( valNoData )
         if NoRispTemp > 0:
            NoRisp = NoRisp + 1
            # Togliere il cancelletto dalle righe sottostanti per visualizzare un avviso ad ogni risposta vuota definitiva
            #print 'Attenzione! risposta vuota al punto x=' + str(ix) + 'y=' + str(iy)
            #print rispostaxml
            NoRispTemp = 0
         if NoPixelValueTemp > 0:
            NoPixelValue = NoPixelValue + 1
            # Togliere il cancelletto dalle righe sottostanti per visualizzare un avviso ad ogni risposta definitiva senza PixelValue
            #print 'Attenzione! Pixel Value assente al punto x=' + str(ix) + 'y=' + str(iy)
            #print rispostaxml
            NoPixelValueTemp = 0
      ...
   ...
   print 'Valori NoData temporanei: ' + str(NoDataCountTemp)
   print 'Valori NoData definitivi: ' + str(NoDataCount)
   print 'Risposte vuote (probabilmente al di fuori della copertura geografica): ' + str(NoRisp)
   print 'Risposte senza PixelValue: ' + str(NoPixelValue)

NoDataCount: numero di volte in cui il parametro PixelValue è uguale a NoData (succede quando viene interrogato un punto all'interno dell'area rettangolare del DTM ma fuori dai confini regionali).
NoRisp: numero di risposte vuote (succede in genere quando viene interrogoto un punto al di fuori della copertura rettangolare del DTM).
NoPixelValue: numero di risposte con schema xml ma senza PixelValue (non dovrebbe succedere).
A volte il servizio WMS sembra distrarsi e restituisce risposte non valide o valori nulli anche quando il punto interrogato ha un valore di quota valido. Per questo lo script insiste nel chiedere la quota finchè non ottiene un valore valido o finchè non si raggiunge il numero massimo di ripetizioni della domanda impostato a 10 (while value == "NoData" and iNoData < 10:). NoRispTemp e NoPixelValueTemp contano le risposte non valide temporanee ricevute in questi momenti di "distrazione". Se dopo 10 ripetizioni della damanda uno di questi contatori temporanei è maggiore di 0 allora viene incrementato il rispettivo contatore definitivo. NoDataCountTemp conte le risposte nulle temporanee a volte molto frequenti; non viene mai azzerato e il suo valore finale viene visualizzato al termine dell'esecuzione per dare un'idea di quanto affidabile è stato il servizio WMS.

Iterazione su righe e colonne:

   iy = 0
   while iy < hgt:
      ix = 0
      row = ''
      while ix < wdt:
         ...
         row = row + ' ' + value
         ix = ix +1
      #
      # scrittura della riga nel file di output
      out_file.write( row + '\n' )
      iy = iy + 1
      end = time.time()
      print 'riga ' + str(iy) + '; tempo impiegato: ' + str( end - start ) + ' secondi'
   ...

L'iterazione viene eseguita con cicli while prima sulle righe poi, al loro interno, sui pixel di ogni riga. Determinato il valore per un pixel esso viene aggiunto alla variabile row (una stringa di testo in cui i valori sono separati da uno spazio). Al termine di ogni riga la stringa row viene scritta nel file di output.

Interrogazione via internet:

            # interrogazione GetFeatureInfo del servizio wms
            request = request2 + '&x=' +str(ix) + '&y=' + str(iy)
            risposta = urllib2.urlopen(request)
            rispostaxml = risposta.read()
            risposta.close()

Una volta ricostruito l'indirizzo per la richiesta GetFeatureInfo (http://www.cartografia.servizirl.it/arcgis/services/wms/DTM5_RL_wms/MapServer/WmsServer?service=wms&version=1.3.0&request=GetFeatureInfo&query_layers=DTM%205X5&info_format=text/xml&crs=EPSG:32632&bbox=xmin,ymin,xmax,ymax&width=width&height=height&x=ix&y=iy)
viene assimilata la risposta via internet grazie alla libreria urllib2 ed in particolare al comando urllib2.urlopen(url). La risposta deve poi essere letta con .read() e salvata in una variabile che in questo caso è rispostaxml; infine è buona norma chiuderla con .close().
urllib vs urllib2: l'interrogazione via internet funziona anche con la libreria urllib tuttavia la versione 2 ha prestazioni in velocità notevolmente migliori; personalmente con urllib per scaricare un'area di 100x100 pixel ho impiagato quasi 2 ore mentre con urllib2 sono bastati 15 minuti.

Estrazione degli attributi dalla risposta xml:
Con la funzione GetFeatureInfo è stato chiesta la risposta in formato xml (info_format=text/xml) da cui ora possiamo estrarre gli attributi grazie alla libreria xml.etree.ElementTree.

            # estrazione degli attributi dalla risposta xml
            FeatureInfoResponse = ET.fromstring(rispostaxml)
            #
            # se la risposta è vuota viene incrementato il contatore relativo
            if len( FeatureInfoResponse ) == 0:
               NoPixelValueTemp = 0
               NoRispTemp = NoRispTemp + 1
            # altrimenti viene cercato l'attributo PixelValue
            else:
               NoRispTemp = 0
               fieldxml = FeatureInfoResponse[0]
               if 'PixelValue' in fieldxml.attrib:
                  NoPixelValueTemp = 0
                  value = fieldxml.attrib['PixelValue']

import xml.etree.ElementTree as ET   (vedi sopra: librerie) permette di abbreviare i riferimenti alla libreria con la sigla ET.
ET.fromstring(rispostaxml)   converte la rispostaxml in un formato da cui è possibile estrarre i componenti (_children). Visto che la nosrta risposta ha un solo componente possiamo indicarlo con il suo numero d'ordine che è [0] ed estrarne gli attributi con   .attrib['attributo'].

Se abbiamo una risposta xml tipo:
<FeatureInfoResponse xmlns:esri_wms="http://www.esri.com/wms" xmlns="http://www.esri.com/wms">
   <FIELDS Stretchedvalue="187" PixelValue="1568.901855"></FIELDS>
</FeatureInfoResponse>
ecco alcuni possibili comandi e i relativi risultati:

rispostaxml
'<?xml version="1.0" encoding="UTF-8"?>\n\r\n<FeatureInfoResponse xmlns:esri_wms="http://www.esri.com/wms" xmlns="http://www.esri.com/wms">\r\n<FIELDS Stretchedvalue="187" PixelValue="1568.901855"></FIELDS>\r\n</FeatureInfoResponse>\r\n'

ET.fromstring(rispostaxml)
<Element '{http://www.esri.com/wms}FeatureInfoResponse' at 0x20da7a20>

ET.fromstring(rispostaxml).__dict__
{'text': '\n', 'attrib': {}, 'tag': '{http://www.esri.com/wms}FeatureInfoResponse', '_children': [<Element '{http://www.esri.com/wms}FIELDS' at 0x20da7e48>]}

ET.fromstring(rispostaxml)._children
[<Element '{http://www.esri.com/wms}FIELDS' at 0x20da7e48>]

ET.fromstring(rispostaxml).attrib
{}

ET.fromstring(rispostaxml)[0]
<Element '{http://www.esri.com/wms}FIELDS' at 0x20da7e48>

ET.fromstring(rispostaxml)[0].__dict__
{'attrib': {'Stretchedvalue': '187', 'PixelValue': '1568.901855'}, 'tag': '{http://www.esri.com/wms}FIELDS', 'tail': '\n', '_children': []}

ET.fromstring(rispostaxml)[0].attrib
{'Stretchedvalue': '187', 'PixelValue': '1568.901855'}

ET.fromstring(rispostaxml)[0].attrib['PixelValue']
'1568.901855'

Pagina visitata volte

Argomenti correlati

> DTM 5x5 m Lombardia

> QGIS

> DEM/DTM

> WMS in QGIS

> DEMs globali gratuiti

Maggio 2016
Alessandro Perego