sabato 13 dicembre 2008

ogr2ogr, codici epsg e Piani Forestali del Piemonte

In un precedente post ho usato ogr2ogr per convertire gli shapefiles dei Piani Forestali della Regione Piemonte da ED50 UTM 32N a WGS84 UTM 32N.
Quando si parla di conversione tra datum diversi divento diffidente: forse a causa della scarsa conoscenza della materia e delle ripetute fregature avute in passato...e nel presente!
In ogni caso voglio semplicemente provare ad usare ogr2ogr con diversi settaggi, tanto per chiarirmi le idee.
La sintassi (semplificata) di ogr2ogr per la conversione tra sistemi di riferimento diversi (di shapefiles) è la seguente:
ogr2ogr -s_srs (descrizione del sistema di riferimento dell'input)
-t_srs (descrizione del sistema di riferimento dell'output)
(nome shapefile di output) (nome shapefile di input)

Per quanto riguarda la descrizione del sistema di riferimento dell'imput (nel mio caso: ED50 UTM 32N) le possibilità che ho preso in considerazione sono:

  1. una variabile che immagazzina una lunga stringa wkt, come quella riportata qui:

    wkt_ed50='PROJCS["UTM Zone 32, Northern Hemisphere", GEOGCS["international",
    DATUM["European_Datum_1950",
    SPHEROID["International_1924",6378388,297],
    TOWGS84[-87.000,-98.000,-121.000]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]],
    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]]'

    Per ottenere questo wkt ho semplicemente creato una Location di GRASSS in ED50 UTM32N e da shell ho digitato g.proj -w. Lo stesso risultato potrebbe essere ottenuto da qui (fare attenzione ai parametri di traslazione TOWGS84...);
  2. il codice epsg:
    '+init=epsg:23032'
  3. il codice epsg (23032) e i parametri di traslazione:
    '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000'
  4. le stesse informazioni di cui sopra ma in formato proj:
    '+proj=utm +zone=32 +ellps=intl +units=m +towgs84=-87.000,-98.000,-121.000'
Per quanto riguarda l'output in WGS84 UTM 32N ho semplicemente settato al codice epsg 32632.
Ho quindi scritto questo semplice script per testare il tutto su uno shapefile dei Piani Forestali della Regione Piemonte:

#!/bin/sh
dir_work="$HOME/pft"
shpfile="004017"
suffix="_wkt"
ogr2ogr -s_srs "$wkt_ed50" \
    -t_srs '+init=epsg:32632' \
    "$dir_work/${shpfile}${suffix}.shp" \
    "$dir_work/${shpfile}.SHP" suffix="_epsg"

ogr2ogr -s_srs '+init=epsg:23032'\
    -t_srs '+init=epsg:32632' \
    "$dir_work/${shpfile}${suffix}.shp" \
    "$dir_work/${shpfile}.SHP" suffix="_epsgtowg"

ogr2ogr -s_srs '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000' \
    -t_srs '+init=epsg:32632' \
    "$dir_work/${shpfile}${suffix}.shp" \
    "$dir_work/${shpfile}.SHP" suffix="_defepsg"

ogr2ogr -s_srs '+proj=utm +zone=32 +ellps=intl +units=m +towgs84=-87.000,-98.000,-121.000' \
    -t_srs '+init=epsg:32632' \
    "$dir_work/${shpfile}${suffix}.shp" \
    "$dir_work/${shpfile}.SHP" suffix="traspunto"

wine /usr/local/traspunto/traspunto.exe shp pia pia \
    "${dir_work}/${shpfile}" \
    "${dir_work}/${shpfile}${suffix}" ED50 WGS84 32 32

dove nella variabile $wkt_ed50 è contenuto il file wkt visto sopra.

Da notare che uso anche traspunto, disponibile qui, con wine.

Questa è una prova (empirica) nella quale apprezzo visivamente i risultati; non ho qui l'interesse di fare alcuna analisi statistica del tipo: esiste una differenza statistica (e se si, quale?)nel posizionamento spaziale di una serie di punti trasformati con traspunto rispetto a ogr2ogr?

Una volta eseguito lo script visualizzo il tutto in OpenJUMP (ma, ovviamente, vanno bene anche Quantum Gis, gvSIG, udig, ecc).

Il primo confronto interessante è quello tra traspunto (linee in rosso) e ogr2ogr con '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000' (linee in verde).

Per apprezzare la differenza, nella figura ho settato una griglia di sfondo con una maglia di 10 metri.

Il secondo confronto riguarda i due settagggi di ogr2ogr: il primo con il solo codice epsg (in blu) (-s_srs '+init=epsg:23032'); il secondo (in verde) con il codice e i parametri di traslazione ('+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000'). Ho mantenuto la griglia a 10 m.




Com'era prevedibile non esistono differenze tra il risultato della conversione effettuata con il wkt, con '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000' e con
-s_srs '+proj=utm +zone=32 +ellps=intl +units=m +towgs84=-87.000,-98.000,-121.000'.

Conclusioni

Una prova rigorosa su modalità diverse di conversione può solo essere condotta con una serie di punti dei quali si conoscono a priori con esattezza le coordinate nei diversi sistemi di interesse.
Questa è una semplice prova per fissare alcune cose:
  • non fermarsi al codice epsg: cercare, se vi sono, i parametri di rototraslazione relativi alla zona di interesse;
  • nella definizione di una nuova Location porsi il problema di settare i parametri di rototraslazione (vedere g.proj);
  • l'uso di traspunto può essere un'alternativa;
  • l'uso del software Verto unitamente ai grigliati (in vendita presso l'Istituto Geografico Militare ) potrebbe costituire la migliore tra le alternative (?) ...

Nessun commento:

Posta un commento