alsperGIS DATI
SOFTWARE


Accedere ai servizi online ArcGIS REST con QGIS (dati vettoriali)



Alcuni enti, come la Regione Lombardia, offrono l'accesso ad alcuni dati geografici mediante i servizi web di ArcGIS con possibilità di interfaccia REST. Questi dati possono essere visualizzati con un comune web browser o con software ESRI.
Rimanendo sull'esempio della Regione Lombardia ecco alcuni indirizzi web:
- http://www.cartografia.regione.lombardia.it/ArcGIS10/rest/services
- http://www.cartografia.regione.lombardia.it/ArcGIS10P/rest/services
- http://www.cartografia.regione.lombardia.it/ArcGIS93/rest/services
- http://www.cartografia.regione.lombardia.it/ArcGISIIT1/rest/services
da queste si può accedere attraverso sottogruppi fino al dato desiderato (terminante con "/MapServer") che nel caso di dati vettoriali presenta un elenco di layers (con i rispettivi numeri di identificazione tra parentesi) mentre per i dati raster c'è in genere un solo layer tiled (ovvero suddiviso in mosaici con differenti risoluzioni) e vengono indicati i vari livelli di dettaglio.


Al momento in cui scrivo (dicembre 2013), in QGIS manca ancora un'opzione integrata (come puo essere quella per i servizi WMS/WMTS) per accedere a questo tipo di servizio. Alla pagina http://hub.qgis.org/wiki/quantum-gis/Arcgis_rest viene spiegato come acquisire questi dati con QGIS o GDAL; funziona per alcuni servers, ma in altri casi possono esserci dei problemi.
Per i dati raster entrambi i metodi spiegati funzionano bene solo quando ogni livello di dettaglio ha una risoluzione/scala che è esattamente la metà rispetto al livello superiore altrimenti si verificano problemi (vedi http://gis.stackexchange.com/questions/72742/problem-with-gdal-wms-qgis-connecting-to-an-arcgis-server e http://osgeo-org.1560.x6.nabble.com/Problem-connecting-to-ArcGIS-rest-services-td5090795.html). Inoltre capita che dopo un po' il server impedisca l'accesso ai dati secondo criteri che non ho ancora capito. Per tali problemi con i dati raster non sono riuscito a trovare una soluzione.
Per i dati vettoriali le cose vanno meglio. Seguendo le indicazioni alla pagina http://hub.qgis.org/wiki/quantum-gis/Arcgis_rest si otterrebbe un file con tutti gli elementi (features) contenuti nel layer richiesto; in realtà il server può avere un limite di per cui l'utente può acquisire soli i primi 500 o 1000 elementi, ottenendo quindi un file incompleto. Per chi conosce un po' di Python è facile creare script per ottenere quegli elementi che cadono in un'area ristretta in modo da non superare il numero massimo e quindi ottenere un file completo per l'area di proprio interesse.

Lo script che riporto in questa pagina considera l'area della vista (map canvas) di QGIS, chiede al server il solo conteggio degli elementi che cadono in tale area, quindi scarica gli elementi e li conta. Se i due conteggi coincidono allora abbiamo tutte le geometrie, se invece manca qualcosa è meglio selezionare un'area minore.
Per usare questo script non serve conoscer Phyton; è sufficiente modificare le prime quattro righe specificando:
- url: l'indirizzo http principale (terminante con "/MapServer"),
- layers: i numeri identificativi dei layers che si vogliono acquisire (in forma di lista, ovvero tra parentesi quadre e separati da virgole),
- name: il nome base dei layers/files che verranno creati (a cui verrà aggiunto l'ID di ciascun layer),
- temp_folder: la cartella dove salvare i files.
Fatto questo è sufficiente copiare il tutto ed incollare nella console python di QGIS, quindi premere invio due volte.
Questo è lo script (provato su QGIS 2.0.1):


url = "http://www.cartografia.regione.lombardia.it/ArcGISIIT1/rest/services/DBTR/dbtr_dynamic/MapServer"
layers = [21,22]
name = "prova_"
temp_folder = "D:/Dati/_prove_rest-service/dbtr_dynamic"
# from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
from osgeo import ogr
import osgeo.ogr
import urllib
import json
#
# Get map extent
area = iface.mapCanvas().extent()
query_area = str(area.xMinimum())+','+str(area.yMinimum())+','+str(area.xMaximum())+','+str(area.yMaximum())
#
for numlayer in layers:
   #
   # Get count of features in map extent
   query_count_url = url + '/' + str(numlayer) + '/query?geometry=' + query_area + '&geometryType=esriGeometryEnvelope&returnCountOnly=true&f=json'
   query_count = urllib.urlopen (query_count_url)
   count_json = json.loads(query_count.read())
   count = count_json['count']
   if count == 0:
      print  "No features for layer %s." %(numlayer)
   else:
      print "Layer %s has %d features. Acquiring..." %(numlayer,count)
      #
      # Get features in map extent
      input = url + '/' + str(numlayer) + '/query?geometry=' + query_area + '&geometryType=esriGeometryEnvelope&outfields=*&f=json'
      layer_name = name + str(numlayer)
      layer = QgsVectorLayer(input, layer_name, "ogr")
      acqFeature = layer.featureCount()
      if acqFeature == 0:
         print "No features acquired for layer %d." %(numlayer)
      else:
         file = temp_folder + "/" + layer_name + ".shp"
         error = QgsVectorFileWriter.writeAsVectorFormat(layer, file, "CP1250", None, "ESRI Shapefile")
         if error == QgsVectorFileWriter.NoError:
            print "Layer %s acquired: %d features." %(numlayer,acqFeature)
            if acqFeature < count:
               print "WARNING: some features were not acquired."
            print "%s has been created." %(file)
            vlayer = QgsVectorLayer(file, layer_name, "ogr")
            addlayer = QgsMapLayerRegistry.instance().addMapLayer(vlayer)


Per richiedere gli elementi di un layer si interroga il server aggiungendo "/query?" all'url del layer, seguito poi dai criteri di ricerca (url/query?criteri).
Nel mio script ho utilizzato:
- geometry=xmin,ymin,xmax,ymax&geometryType=esriGeometryEnvelope per richiedere solo le geometrie che cadono nell'area rettangolare (envelope) definita dalle coordinate xmin, ymin, xmax, ymax (prese dal map canvas di QGIS)
- outfields=* per acquisire tutti gli attributi degli elementi
- f=json per ottenere la risposta in formato json (che può essere letto da OGR)
- per avere il solo conteggio degli elementi: returnCountOnly=true.



 
Argomenti correlati
>>> QGIS




Collegamenti esterni
>>> ArcGIS_rest via QGIS e GDAL
>>> Query su ArcGIS rest









Dicembre 2013
Alessandro Perego