Skip navigation.
Home
Istituto Scienze della Terra - SUPSI

Conversione di geodati con FwTools v2.1.0

Indice

--> Contesto generale

--> Trasformazioni più frequenti DI DATI VETTORIALI con linee di comando nel Bash di FWTools

--> Trasformazioni DI DATI RASTER con linee di comando nel Bash di FWTools

--> Indicatori di Tissot ed utilizzo passo passo di FwTools


Contesto generale

Vista la crescente facilità nel disporre di dati geografici (acquisizione GPS con cellulare ad esempio) e l'utilizzo sempre più frequente e famigliare della cartografia (si pensi a Google Earth), la conversione di dati tra sistemi di riferimento é un operazione assolutamente indispensabile per permettere all'utente di poter utilizzare queste risorse nel miglior modo possibile. Gli Open Source con FwTools e le librerie ogr2ogr consentono d'effettuare gratuitamente e facilmente questo tipo di lavoro. Viene proposto un esempio pratico che consente al lettore d'imparare ad utilizzare tali strumenti.

Trasformazioni più frequenti DI DATI VETTORIALI con linee di comando nel Bash di FwTools

Link per scaricare FWTools (nel caso degli esempi FWTools 2.1.0 per Windows)

http://fwtools.maptools.org/

Tappa 1 - Nel Bash di Fwtools entrare nella directory nella quale ci sono i dati da trasformare (esempio)

 


Trasformazione di formato shp da sistema svizzero CH1903 (EPSG:21781) a sistema WGS84 (EPSG:4326)

ogr2ogr.exe -f "ESRI Shapefile" -s_srs EPSG:21781 -t_srs EPSG:4326 output.shp input.shp

 


Trasformazione di formato shp da sistema WGS84 (EPSG:4326) a sistema svizzero CH1903 (EPSG:21781)

ogr2ogr.exe -f "ESRI Shapefile" -t_srs EPSG:21781 -s_srs EPSG:4326 output.shp input.shp

da notare che rispetto alla linea di comando precedente cambiano unicamente -t_srs (trasforma in ...) e -s_srs (associa il datum EPSG prima di trasformare in ...)
 


Trasformazione di formato shp da sistema UTM32/WGS84 (EPSG:32632) a sistema svizzero CH1903 (EPSG:21781)

ogr2ogr.exe -f "ESRI Shapefile" -t_srs EPSG:21781 -s_srs EPSG:32632 output.shp input.shp

da notare : la proiezione UTM32 utilizza l'elissoide globale WGS84 ed é una proiezione UTM del fuso 32 utilizzato spesso per rappresentare dati italiani (nord Italia in particolare). Proponiamo questo esempio da UTM32 a CH1903 visto che per molti progetti GIS transfrontalieri (tra Ticino e Lombardia per esempio) si ha la necessità di trasformare dati da un sistema di riferimento ad un altro e viceversa. Come esempio un punto sul Monte Bisbino (tra Mendrisiotto e provincia di Como) con coordinate UTM32 : E 505123.83275, N 5080063.1649  --> CH1903 : Y 726361.171435854 X 81587.205339726


 

Trasformazione di formato shp da sistema di riferimento WGS84 a sistema rumeno stereo70 (EPSG specificato con il Datum e non tramite codice)

ogr2ogr.exe -f "ESRI Shapefile" -s_srs EPSG:4326 -t_srs "+proj=stere +lat_0=46 +lon_0=25  +k=0.999750 +x_0=500000 +y_0=500000 \ 
+ellps=krass +towgs84=28,-121,-77,0,0,0,0 +units=m +no_defs no_defs <>" output.shp imput.shp

da notare che anche se non esiste un codice EPSG che identifica il sistema di riferimento (in questo caso Stereo70) c'é la possibilità di mettere direttamente a comando il Datum (ellissoide di riferimento e proiezione)

 

Codice che permette la conversione di coordinate in modo testuale senza passare via un formato geografico come lo shp. Esempio inserisco lon e lat (25.769943 45.583957)
e mi restituisce le coordinate rumente s70

cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=stere +lat_0=46 +lon_0=25  +k=0.999750 +x_0=500000 \ 
+y_0=500000 +ellps=krass +towgs84=28,-121,-77,0,0,0,0 +units=m +no_defs
[enter]
25.769943 45.583957
[enter]


Per trasformare più coordinate in blocco e stamparle su di un file fare come segue :

1) inserire il comando di cs2cs (ad esempio come sopra) mettendo il nome del file di testo da creare con i risultati
con > nome_file_output (esempio > stampa_risultato.txt) e schiacciare enter

2) inserire le coordinate (ogni riga rappresenta una coppia di coordinate), enter

3) CRTL+Z per uscire dall'applicazione ed aggiornare il file output


 


Trasformazione da formato kml (Google Earth già in WGS84) a formato shp in sistema svizzero CH1903 (EPSG:21781)

ogr2ogr.exe -f "ESRI Shapefile" output.shp input.kml -t_srs "EPSG:21781"

attenzione : il kml di input deve avere geometria unica senò errore. Se ci sono geometrie miste nel kml il comando non funziona
 


Trasformazione da formato shp in sistema svizzero CH1903 (EPSG:21781) a formato kml (Google Earth già in WGS84)

ogr2ogr.exe –s_srs "EPSG:21781" –f "KML" output.kml input.shp  -t_srs "EPSG:4326"

 


Trasformazione nella proiezione di Mercator di uno shp

ogr2ogr.exe -t_srs "+proj=merc" output_merc.shp input.shp

da notare che come per altre proiezioni già l'abbreviazione proj=merc é riconosciuta come proiezione di Mercator


Trasformazioni DI DATI RASTER con linee di comando nel Bash di FwTools

Abbiamo visto precedentemente che per la trasformazione di dati VETTORIALI nei vari sistemi di riferimento in FWTools si utilizzano le librerie ogr2ogr. Invece, quando si tratta di dati RASTER, magari un DTM o un Ortofoto che vogliamo convertire (per esempio dal WGS 84 --> al CH1903 o un altro sistema) utilizziamo le librerie gdal ed in particolare gdalwarp

con il comando gdalwarp --formats posso vedere tutti i formati supportati in input da gdal. C'é da dire che il GeoTiff é il formato di output per default, ma eventualmente per definire un altro formato output c'é da inserire

-of format:
    Select the output format. The default is GeoTIFF (GTiff). Use the short format name

Esempio di schermata FWTools e utilizzo di gdalwarp

 

Tappa 1 - Nel Bash di Fwtools entrare nella directory nella quale ci sono i dati da trasformare

 


Trasformazione di un DTM nel formato GeoTiff a GeoTiff  da sistema WGS84 (EPSG:4326) a sistema svizzero CH1903 (EPSG:21781)

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:21781 input.tif output.tif

da notare che rispetto alla linea di comando precedente cambiano unicamente -t_srs (trasforma in ...) e -s_srs (associa il datum EPSG prima di trasformare in ...)


Trasformazione di un DTM nel formato GeoTiff a GeoTiff da sistema di riferimento WGS84 a sistema rumeno stereo70 (EPSG specificato con il Datum e non tramite codice)

gdalwarp -s_srs EPSG:4326 -t_srs "+proj=stere +lat_0=46 +lon_0=25 +k=0.999750 +x_0=500000 +y_0=500000  +ellps=krass \
+towgs84=28,-121,-77,0,0,0,0 +units=m +no_defs  no_defs <>" input.tif output.tif

da notare che anche se non esiste un codice EPSG che identifica il sistema di riferimento (in questo caso Stereo70) c'é la possibilità di mettere direttamente a comando il Datum (ellissoide di riferimento e proiezione)


 

Link utili per questo argomento :

http://casoilresource.lawr.ucdavis.edu/drupal/node/98

http://www.gdal.org/gdalwarp.html

http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp

(nell'utlimo sito possibilità di scaricare dati DTM 90m gratuiti in tutto il mondo anche nel formato GeoTiff che é riconosciuto e consigliato quando utilizzo gdal. I DTM scaricati sono in formato WGS84 ma hanno l'informazione di quota non nell'ellissoide WGS84 ma già corretta a quote usuali o s.l.m.


 

Utilizzo di gdal_translate : all'interno delle librerie gdal, oltre a gdalwarp appena visto, abbiamo la possibilità di convertire dati raster in formati diversi, comprimerli e trattarli in modo tale da poter essere utilizzati nel modo più adeguato possibile magari all'interno di un geoservizio. Esempio: in Geoserver con Openlayers ho la priorità di vedere un raster con un buon dettaglio ma al tempo stesso questo dev'essere il più veloce possibile quando faccio le operazioni di navigazione della mappa quali zoom, pan, ecc. Ecco che utilizzerò FWTools con gdal_translate per ottenere il raster più adatto.

 

In FWTools

- comando gdal_translate : mi da semplicemente l'elenco dei formati raster supportati

- gdalinfo merano.jpg (dove merano.jpg é il nome del raster) : mi da tutte le informazioni sull'header del raster, se questo é compresso, se ci sono tile o overview, ecc.

- tradurre un .jpg in .geotiff (formato di default per l'output)

gdal_translate merano.jpg merano.tiff

- tradurre un .jpg in .geotiff dandogli diverse caratteristiche come assegnare il sistema di coordinate che prima non esisteva (attenzione non fa una trasformazione), dividere il raster in tile e queste tile le comprime in jpg con una certa qualità d'immagine. Da notare che già questa operazione permette d'ottenere un'imagine che con Geoserver é bella veloce.

gdal_translate a_srs EPSG 4326 -co "TILED=YES" -co "COMPRESS=JPEG" -co "JPEG_QUALITY=90" merano.jpg merano.tiff

 

Utili (lista del formati gdal + specifiche possibili per ogni formato con gdal_translate) :

http://www.gdal.org/formats_list.html

http://www.gdal.org/frmt_gtiff.html 

nota: i file non compressi sono i più veloci in Geoserver ma occupano molta memori nel disco duro. Inoltre non posso caricare immagini maggiori di 2 GB. per questa ragione cerco di comprimere fin quando trovo un giusto compromesso tra velocità e memoria occupata dal file. In ogni caso in Geoserver vengono usati anche degli acceleratori detti *.j e *.jjo, ma questa é un'altra storia. Vedi piuttosto corso di Andrea Aime su Geoserver (IST-SUPSI, dicembre 2008)


 

Indicatori di Tissot ed utilizzo passo passo di FwTools

Problematica : tramite uno script in Python (file generate_tissot.py) creato da PerryGeo (vedi http://www.perrygeo.net/wordpress/?p=4 oppure http://yukongis.ca/bin/view/Main/ProjectionDistortion) creo gli indicatori di Tissot in formato .shp all’interno di una region specifica con intervalli e grandezze definite (invece utilizzando lo script senza modificarlo gli indicatori sono creati per l’intero mondo).

Questi indicatori, come pure altri shape, per esempio quelli del mondo, posso riproiettarli tramite le librerie di ogr2ogr utilizzando l’FWTools.

Gli indicatori di Tissot riproiettati nel sistema voluto sono molto preziosi perché consentono di capire esattamente l’influenza della proiezione sulla deformazione delle figure. Salta per esempio subito all’occhio se si tratta di una proiezione conforme (conserva gli angoli) o equivalente (conserva le superfici). Essendo questi indicatori degli shape possiamo direttamente misurarli tirando fuori delle quantità per eventuali apprezzamenti, comparazioni (esempio sapere che per la proiezione svizzera abbiamo una deformazione di ca. 30 cm / km nel sud del Ticino).

Indicatori di Tissot della Proiezione di Mercator

 

Indicatori di Tissot della Proiezione di Lambert

 

Compreso lo scopo dell’esercizio, qui spiegherò principalmente come utilizzare “a livello informatico” lo script di OGR/Python e l’ogr2ogr per le trasformazioni tra sistemi di riferimento. Il tutto con l’ausilio di FWTools.

 

Tappa 1) Creo lo shp degli indicatori di Tissot (tissot.shp) con lo script OGR/Python

Premessa: il file generate_tissot.py (scarica) si trova nella cartella
D:\Alessio_Spataro\Proiezioni_cartografiche

1.    Apro FwTools
2.    entro in D: (enter) --> cd Alessio_Spataro (enter) --> cd Proiezioni_cartografiche (enter) e poi chiamo l’applicazione ogr2ogr.exe (enter)

Questa operazione serve per “settare” l’FWTools in modo tale che d’ora in avanti posso lavorare con file all’interno della cartella D:\Alessio_Spataro\Proiezioni_cartografiche senza dover ogni volta specificare il percorso del file input ed output. Nel caso in cui volessi specificare il percorso del file dovrei farlo come segue, vedi esempio:

C:\Program Files\FWTools2.1.0\bin>ogr2ogr.exe –s_srs “EPSG:21781” –f “KML” D:\Alessio_Spataro\Proiezioni_cartografiche\output_prova.kml D:\Alessio_Spataro\Proiezioni_cartografiche\input_prova.shp  -t_srs “EPSG:4326”

Nell’esempio trasformo sempre tramite ogr2ogr un file .shp con proj CH1903 in .kml con proj WGS84 per utilizzarlo poi in google map.

3.    Chiamo nella finesta FWTools il mio script in Python che mi crea gli .shp di Tissot in questo modo: python generate_tissot.py (enter)

L’FWTools esegue lo script e crea tissot.shp con le sue caratteristiche come definito nello script. Attenzione: lo script si può modificare ad esempio settando una region diversa, con intervalli diversi o dimensione dei cerchi di Tissot più grandi o più piccoli. Vedi a titolo d’esempio generate_tissot_ch.py e comparalo con generate_tissot.py

 

Esempio generate_tissot_ch.py (crea i tissot per la region Svizzera)

#!/usr/bin/env python
# Tissot Circles
# Represent perfect circles of equal area on a globe
# but will appear distorted in ANY 2d projection.
# Used to show the size, shape and directional distortion
# by Matthew T. Perry
# 12/10/2005
   
import ogr
import os
import osr
   
output = 'tissot_ch.shp'
debug = False
   
# Create the Shapefile
driver = ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists(output):
        driver.DeleteDataSource(output)
ds = driver.CreateDataSource(output)
layer = ds.CreateLayer(output, geom_type=ogr.wkbPolygon)
   
# Set up spatial reference systems
latlong = osr.SpatialReference()
ortho = osr.SpatialReference()
latlong.ImportFromProj4('+proj=latlong')
   
# For each grid point, reproject to ortho centered on itself,
# buffer by 640,000 meters, reproject back to latlong,
# and output the latlong ellipse to shapefile
for x in range(50,150,1):
    for y in range (450,550,1):
        f= ogr.Feature(feature_def=layer.GetLayerDefn())
        wkt = 'POINT(%f %f)' % (float(x)/10, float(y)/10)
        p = ogr.CreateGeometryFromWkt(wkt)
        p.AssignSpatialReference(latlong)
        proj = '+proj=ortho +lon_0=%f +lat_0=%f' % (float(x)/10,float(y)/10)
        ortho.ImportFromProj4(proj)
        p.TransformTo(ortho)
        b = p.Buffer(2500)
        b.AssignSpatialReference(ortho)
        b.TransformTo(latlong)
        f.SetGeometryDirectly(b)
        layer.CreateFeature(f)
        f.Destroy()
   
ds.Destroy()

Tappa 2) Dopo aver creato il tissot.sh posso divertirmi, magari insieme allo shp delle nazioni Mondo.shp già in dotazione, ad effettuare le trasformazioni cartografiche classiche oppure le più diverse. Ecco qualche esempio per capire la sintassi per la trasformazione tra sistemi con FWTools

Proiezione di Mercator , prima delle nazioni Mondo.shp (input), poi tissot.shp (input). Il risultato (output) è Mondo_merc.shp e rispettivamente tissot_merc.shp

ogr2ogr -t_srs “+proj=merc” Mondo_merc.shp Mondo.shp (enter)

ogr2ogr -t_srs “+proj=merc” tissot_merc.shp tissot.shp (enter)

Proiezione di Lambert (corrisponde ad epsg:2163)

ogr2ogr -t_srs “epsg:2163″ Mondo_lambert.shp Mondo.shp
(enter)

ogr2ogr -t_srs “epsg:2163″ tissot_lambert.shp tissot.shp (enter)

 

Si consiglia per comodità di copia/incolla di lavorare con un editor di testo (Notepad o altro) per elaborare il testo che poi verrà copiato nella finestra FWTools (Copy e Paste con il mouse)

 

Esempio (sotto) di visualizzazione in ArcGis di Mondo_lamb.shp e tissot_shp.shp

L’indicatore di Tissot è molto importante. Con questo si può apprezzare l’entità della deformazione cartografica in termini quantitativi. Per farlo si consiglia di utilizzare le funzioni di Calcolate Values … di ArcGis comparando il tissot.shp originale con quello della trasformazione studiata.

 ********************************************************************************************************************************************************

altri comandi utili frequenti in Dos (funzionali per lavorare per chi se li dimentica facilmente)

cd .. (torna indietro di cartella o folder) es. D:\Alessio\prova , con cd .. torna a D:\Alessio\

mkdir /Prova es. crea la cartella Prova

mv *.aux/Prova es. Sposta (move) tutti i file di estensione .aux nella cartella Prova

ls >> lista.txt Crea un file lista.txt con la lista dei nomi di tutti i file contenuti in una cartella

MORE file.txt es. MORE lasersc.txt (nel prompt di dos in windows) : legge le prime righe di un file anche molto pesante. Particolarmente utile per file di rilevo es. laserscanner dei quali non si conosce la struttura. Per copiare poi le prime righe visibili nel prompt in un altro file fare:

Edit (destro mouse sulla cornice Prompt) --> Mark (selezionare a schermo le linee volute) --> Enter (funziona come un copia) --> aprire un file ed incollare (paste).

help : presenta tutti i comandi possibili