.quanti Sistemi di Riferimento conosci? Io almeno 6000!

Sono incappato in questo sito epsg.io nei giorni scorsi e l’ho trovato molto interessante oltre che utile. Si tratta di un progetto libero e a codice aperto disponibile anche su Github, che consente di recuperare molte informazioni sulla quasi totalità dei sitemi di riferimento spaziali (CRS – Coordinate Reference System) esistenti utilizzati in tutto il mondo, sia globali che locali. Il numero di CRS presenti attualmente è quello dell’archivio EPSG (European Petroleum Survey Group), divenuto ormai uno standard internazionale per la codifica dei sistemi geodetici.

Non sarò ripetitivo su cos’è e come è nato l’archivio EPSG, ma vi rimando alla lettura di alcune risorse che aiutano ad approfondire l’argomento:

  1. Nota Istituto Geografico Militare
  2. Registro ufficiale EPSG dell’IOGP (International Association of Oil & Gas Producers)

Diciamo che il mio sito di riferimento per i sistemi geodetici è sempre stato http://www.spatialreference.org insieme a http://www.epsg-registry.org, ma devo ammettere che EPSG.io è un’ottima alternativa. Il sito ha molte funzionalità che vanno dalla semplice ricerca testuale, inserendo per esempio il nome del paese, del codice epsg o direttamente il nome del sistema di riferimento, ad una ricerca più avanzata passando dei dati di input attraverso una query string, sotto forma di coppie &chiave=valore, direttamente nell’URL del servizio. Il risultato della ricerca può essere comodamente esportato in formato JSON. Forte!

Schermata 2016-04-17 alle 21.26.18

Anche l’utilizzo dell’interfaccia web non è male: per sapere tutto sul Sistema Monte Mario/Italy 2 basta andare su epsg.io/3004, wow! un URL così corto sarà difficile da dimenticare! E’ possibile anche scaricare la definizione del CRS per differenti formati, per esempio epsg.io/3004.wkt?download scaricherà il file 3004.prj, in un modo molto simile si possono scaricare tutti gli altri formati disponibili. Di seguito un estratto dalla pagina del sito dove sono elencate le principali funzionalità offerte da EPSG.io.

Schermata 2016-04-17 alle 08.56.39

Vediamo adesso come sfruttare questo servizio, recuperando informazioni sul sistema di riferimento desiderato. In particolare, facendo uso di poche righe di codice, vedremo come esportare le definizioni del sistema di riferimento scelto in due formati differenti, PROJ4 e WKT.

Iniziamo con il leggere il risultato dalla pagina del servizio. Supponiamo di voler ottenere le informazioni sulla definizione del sistema di coordinate EPSG:3004, Monte Mario/ Italy 2, per prima cosa dobbiamo recuperare il file JSON dal servizio EPSG.io attraverso la seguente query string:

?q=3004&format=json&trans=1

(più informazioni sulle possibili query si possono trovare consultando direttamente il README presente nella repository Github)

quindi copiamo ed incolliamo il seguente codice

import json
import urllib

EPSG = '3004'

url = 'http://epsg.io/?q=' + EPSG + '&format=json&trans=1'
epsgio_response = urllib.urlopen(url)
json_response = json.loads(epsgio_response.read())

Una volta eseguito il codice sopra, all’interno della variabile json_response troveremo tutto quello che ci serve in formato JSON. Alcuni CRS possono contenere più definizioni per le diverse trasformazioni (i famosi TOWGS84 o anche conosciuti come parametri di Bursa Wolf), il 3004 per esempio ne ha 11. La keyword da usare nel JSON per ottenere tutte le transformazioni disponibili è [“trans”]. Associamo le trasformazioni così recuperate alla variabile transformations e stampiamone il risultato che in questo caso sarà il codice EPSG del CRS, la sua stringa di definizione PROJ4 (nell’esempio sotto l’ho accorciata per motivi di spazio) e l’area di riferimento del CRS selezionato

transformations = json_response['results'][0]['trans']

if transformations:
    for i in range(len(transformations)):
        # transformation code
        print transformations[i]['code_trans']
        # proj4 definition
        print transformations[i]['proj4']
        # area of use of the transformation
        print transformations[i]['area']
1088
+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=in...
Italy - offshore - Adriatic Sea - North Ancona.
1089
+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=in...
Italy - offshore - Adriatic Sea - South Ancona and North Gargano.
1090
+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=in...
Italy - offshore - Adriatic Sea - South Gargano.
1091
+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=in...
Italy - offshore - Otranto channel.
1092
+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=in...
Italy - offshore - north Ionian Sea.
1093
+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=in...
Italy - offshore - Strait of Sicily - east of 13°E (of Greenwich).
1094
+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=in...
Italy - offshore - Strait of Sicily - west of 13°E (of Greenwich).
1169
+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=in...
Italy - Sardinia onshore.
1660
+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=in...
Italy - mainland including San Marino and Vatican City State.
1662
+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=in...
Italy - Sardinia onshore.
1664
+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=in...
Italy - Sicily onshore.

Adesso, dalla console python di QGIS, eseguiremo due operazioni, nella prima utilizzeremo la definizione PROJ4 per associare un sistema di riferimento ad un layer caricato in QGIS e nella seconda useremo la stringa in WKT per creare un file *.qpj.

Il file *.qpj è un file molto simile al *.prj che troviamo comunemente associato insieme ai file che compongono uno shapefile, consente al software di riconoscere il sistema di riferimento del dataset indispensabile per eseguire le corrette trasformazioni sul dataset stesso tra i diversi sistemi di riferimento.

Per eseguire le operazioni useremo la trasformazione di default del 3004, cioè la 1660 (in grassetto). La prima operazione consiste quindi nell’associare il CRS ad un layer caricato in QGIS attraverso la stringa PROJ4

crs = QgsCoordinateReferenceSystem()
crs.createFromProj4(json_response['results'][0]['proj4']) 

mapLayer = iface.activeLayer()
mapLayer.setCrs(crs)

Mentre, con la seconda operazione andiamo a scrivere la stringa WKT di definizione del CRS all’interno di un file che dovrà avere lo stesso nome dello shapefile a cui sarà associato {nome_shapefile}.qpj

with open('path/shapefile_name.qpj', 'w') as f:
    f.write(json_response['results'][0]['wkt'])

Entrambe le operazioni associano solamente il corretto sistema di coordinate al dataset, non eseguono alcuna riproiezione.

Finito, cotto e mangiato!

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s