Use the left/right arrow keys to navigate, 's' to enable/disable scrolling.

Introduzione ai GIS Liberi


CLI tools (PROJ4 - Gdal/OGR)

Indice


PROJ4

GDAL/OGR

Il formato HDF

PROJ4 utils


La libreria proj4 è alla base di tutti i processi di proiezione, conversione tra datum e passaggio tra sistemi di riferimento nei software geografici liberi


La distribuzione standard è corredata di alcuni utili programmi a linea di comando:


  • proj: effettua una proiezione diretta (da coordinate geografiche a coordinate cartesiane)
  • invproj: effettua una proiezione inversa (da coordinate cartesiane a coordinate geografiche)
  • cs2cs: effettua conversione di coordinate da un crs ad un altro
  • geod/invgeod: effettua calcoli geodetici (diretti e inversi)

Conversione diretta (proj)


Il comando proj è un tool che permette la trasformazione diretta di coordinate (da coordinate geografiche a coordinate cartografiche), all’interno dello stesso datum


È implementato come un classico filtro bash (legge dallo standard input e scrive sullo standard output)


Prende in input una coppia di coordinate geografiche (f,l), e le proietta sul sistema di riferimento cartografico specificato sulla linea di comando

EPSG


La distribuzione standard di proj comprende una serie di file che definiscono i principali sistemi di riferimento (geografici e cartografici) esistenti.


Ai fini del corso faremo riferimento unicamente al file /usr/share/proj/epsg, che contiente il database prodotto dall’EPSG, con le definizioni nella sintassi di proj


Per utilizzare un sistema di riferimento definito in tale file, è sufficiente utilizzare l’opzione +init=”epsg:srsid”

Conversione diretta (proj)



Es.: conversione da WGS84 a UTM33N (mediante file di inizializzazione):


$> echo "16E 45N"| proj -v +init="epsg:32633"

Es.: conversione da WGS84 a UTM33N (mediante definizione della proiezione):


$> echo "16E 45N"| proj -v +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs

Conversione inversa (invproj)


Il comando invproj effettua la trasformazione inversa di coordinate, all’interno dello stesso datum (da coordinate cartografiche a coordinate geografiche)


La sintassi richiesta è analoga a quella di proj


Nota: è possibile sostituire il comando invproj con il seguente: proj -i


Esempi:


$> echo "578815.30 4983436.77"| invproj +init="epsg:32633"
$> echo "578815.30 4983436.77"| proj -I +init="epsg:32633"

Formato dell’output


È possibile modificare il formato dell’output di proj e invproj, mediante l’opzione -f, che utilizza la stessa sintassi della classica printf del C:




$> echo "16E 45N"| proj -v -f "%.1f" +init="epsg:32633"

$> echo "578814 4983432"| proj -I -f "%.5f"  +init="epsg:32633"

Ellissoidi supportati


La libreria proj implementa nativamente un elevato numero di ellissoidi, grazie ai quali è possibile definire SRS


Per conoscere la lista degli ellissoidi supportati, lanciare il comando proj -le (list ellipsoid)


Tra gli ellissoidi supportati, sono di interesse ai fini del corso i seguenti:


  • intl a=6378388.0 rf=297 International 1909 (Hayford)
  • WGS84 a=6378137.0 rf=298.257223563 WGS84
  • bessel a=6377397.155 rf=299.1528128 Bessel 1841

Proiezioni supportate


Come per gli ellissoidi, proj implementa internamente un insieme di proiezioni utilizzabili in fase di definizione degli SRS


Una descrizione approfondita delle proiezioni disponibili è reperibile nella documentazione ufficiale del progetto:

proj3 Main User Manual


È possibile ottenere la lista delle proiezioni con il comando proj -lp (list projections):


  • utm : Universal Transverse Mercator (UTM)
  • tmerc : Transverse Mercator
  • sinu : Sinusoidal (Sanson-Flamsteed)

List datum supportati


È disponibile anche un esiguo numero di datum interni, la cui lista è ottenibile con il comando proj -ld (list datum)


La definizione di altri datum è possibile mediante alcuni parametri del comando proj (in particolare +towgs e +nadgrids)


Principali datum interni:


  • WGS84 WGS84 towgs84=0,0,0
  • NAD83 GRS80 towgs84=0,0,0 (North American Datum 1983)
  • NAD27 clrk66 nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat (North American Datum 1927)

Fattori di conversione metrica


Nel passaggio da unità cartesiane a unità metriche, proj può utilizzare una lunga lista di fattori di conversione, che è possibile visualizzare con il comando proj -lu (list units):



km 1000. Kilometer
m 1. Meter
kmi 1852.0 International Nautical Mile
in 0.0254 International Inch
ft 0.3048 International Foot
us-in 1./39.37 U.S. Surveyor's Inch
us-ft 0.304800609601219 U.S. Surveyor's Foot
[---]

Definizione sistemi di rif.to (1/2)


Segue una lista dei principali parametri di input per la definizione dei sistemi di riferimento:


OpzioneSignificato
+adefinizione del semiasse maggiore
+bdefinizione del semiasse minore
+datumnome del datum (proj -ld)
+ellpsnome dell'ellissoide (proj -le)
+k / +k_0fattore di scala
+lat_0latitudine dell'origine
+lon_0meridiano centrale

Definizione sistemi di rif.to (2/2)


OpzioneSignificato
+towgs84parametri (3 o 7) per il datum shifting (trasformazione di Helmert)
+nadgridsfile in formato ntv1/ntv2, per la trasformazione mediante reticolo
+no_defsinibisce l'uso del file /usr/share/proj/proj_def.dat per i parametri di default
+projnome della proiezione da utilizzare (proj -lp)
+zonedefinizione della zona UTM
+x_0falso est
+y_0falso nord
+unitsfattore di conversione da unità cartesiane a unità metriche

Esempi di definizioni SRS


WGS84 (coordinate geografiche)

+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs


Roma40 (con datum shifting)

+proj=longlat +ellps=intl +towgs=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +no_defs


UTM33N

+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs


GaussBoaga Est

+proj=tmerc +towgs=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m +no_defs

cs2cs


Il programma cs2cs effettua trasformazioni di coordinate tra un SRS origine ed uno destinazione (anche tra datum diversi)


La sintassi base del comando richiede la definizione di due sistemi di riferimento, separati dall’opzione +to (tutti i parametri successivi a tale opzione sono considerati relativi all’SRS target)


Nel caso in cui non venga fornito il secondo SRS viene implicitamente utilizzato il sistema di riferimento geografico relativo allo stesso datum del primo


Nota: i sistemi di riferimento specificati possono essere entrambi geografici, entrambi cartografici o misti

Esempi di utilizzo di cs2cs


Trasformazione da GBEst a Roma40

echo “2520000 456789” | cs2cs +init=”epsg:3004”


Trasformazione da Roma40 a UTM33N (senza datum shifting)

echo “16 45” | cs2cs +init=”epsg:4230” +to +proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs


echo “16 45” | cs2cs +init=”epsg:4230” +to +init=”epsg:32633”


Trasformazione da Roma40 a WGS84 (con datum shifting)

echo “16 45” | cs2cs +proj=latlong +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +no_defs +to +init=”epsg:32633”

GDAL/OGR utils


GDAL è una libreria di gestione di formati raster georeferenziati. OGR è l’analogo di GDAL per i formati vettoriali


La distribuzione standard della libreria raccoglie una serie di utili programmi a linea di comando per la traduzione interformato ed il goeprocessing:


  • gdalinfo: reperimento informazioni su raster georiferiti
  • gdal_translate: conversione tra formati
  • gdalwarp: riproiezione raster

  • ogrinfo: reperimento informazioni su dati vettoriali
  • ogr2ogr: conversione tra formati vettoriali

gdalinfo (1/2)


gdalinfo fornisce informazioni sulla sorgente raster in input


Le informazioni mostrate si dividono in tre categorie:

  • informazioni sull’immagine (dimensione dei pixel, formato…)
  • informazioni di carattere geografico (reperite dai metadati interni, se il formato lo permette, o da file specifici - es. world file)
  • informazioni di carattere statistico (numero di bande, distribuzione dei valori, …)

gdalinfo (2/2)


Lista dei formati supportati:

gdalinfo –formats

Informazioni su un formato:

gdalinfo –format GTiff


Uso base:

gdalinfo file_raster


Opzioni utili:

-noct: sopprime l’outpput delle color tables
-stats: calcola e mostra le statistiche

gdal_translate (1/2)


Il programma gdal_translate permette la conversione di dati raster tra formati differenti


All’interno dello stesso task è possibile istruire il processo per effettuare operazioni di rescaling, resampling e subsetting


Sintassi base:

gdal_translate [opzioni] src_dataset dst_dataset

gdal_translate (2/2)


OpzioneSignificato
-offormato di output
-otdata type della banda esportata
-bbanda/e da esportare
-a_srssovrascrive le informazioni sul datum
-outsizerescaling dell'output
-scaleresampling dell'output
-srcwinsubsetting
-moaggiunta di metatag all'output
-coopzioni di creazione specifiche del formato (v.di gdalinfo -format format_name)

Esempi:

gdal_translate -of JPEG data.tiff data.jpg


gdal_translate -of JPEG -co “PROGRESSIVE=TRUE” -co “QUALITY=90” -mo “PIPPO=DATA” -a_srs epsg:3004 data.tif data.jpg

gdalwarp (1/2)


Il programma gdalwarp permette mosaicatura, riproiezione e warping di dati raster tra due sistemi di riferimento differenti


Sintassi base:

gdalwarp [opzioni] src_dataset* dst_dataset

gdalwarp (2/2)


Opzioni comuni:

-s_srssource SRS
-t_srsdestination SRS
-orderordine del polinomio interpolatore (da 1 a 3)
-teestensione geografica dell'output (extension)
-trrisoluzione dell'output (in unità di misura del SRS target)
-rmetodo di ricampionamento (near, bilinear, cubic, cubicspline, lanczos)
-multimodalità di processamento multicore (parallela)
-offormato dell'output
-coopzioni di creazione dell'output (specifiche del formato)

Esempio:

gdalwarp -t_srs epsg:3004 -multi -order cubic srtm_40_04.tif srtm_40_05.tif dem.tif

Piramidi (1/2)


Generalmente i dati raster (e in particolar modo le immagini georiferite) sono caratterizzati da dimensioni notevoli ed elevate risoluzioni


Per velocizzare le operazioni di zoom (in e out) e migliorare la resa visiva nell’utilizzo all’interno di SW GIS, è possibile precalcolare viste della stessa immagine a risoluzioni inferiori, dette piramidi o overviews


Le piramidi, a seconda del formato utilizzato, possono essere memorizzate sia internamente al file di origine (ad es. nei raster GeoTiff), che in file esterni (estensione .ovr)

Piramidi (2/2)


La libreria gdal fornisce il comando gdaladdo per il calcolo delle piramidi


Sintassi base:

gdaladdo [opzioni] src_data livelli


Esempi

gdaladdo -r cubic data.tif 2 4 6 8 16
gdaladdo -r cubic -ro –config COMPRESS_OVERVIEW JPEG data.tif 2 4 6 8 16

gdaldem (1/2)


Il comando gdaldem offre una serie di utili modalità per la visualizzazione e l’analisi spaziale di file DEM (Digital Elevation Model)


Questo tipo di raster è molto utilizzato nello studio del suolo (pendenze, ombre, rischio idrogeologico, ecc.), ed è caratterizzato da celle che memorizzano i valori altimetrici del corrispondente punto della superficie Terrestre

gdaldem (2/2)


Il comando offre sette modalità di funzionamento:

  • hillshade: mappa dei rilievi ombreggiati (simil-3D)
  • slope: mappa delle pendenze
  • aspect: mappa dell’orientamento delle pendenze (azimut)
  • color-relief: mappa dei rilievi a colori
  • TRI (Terrain Ruggedness Index): carta dell’indice di asprezza dei pendii (differenza media tra un pixel centrale e gli otto che lo circondano)
  • TPI (Topographic Position Index): indeice definito come differenza tra un punto centrale e i pixel intorno
  • roughness: differenza massima tra un punto centrale e il suo intorno

Nota: ogni modalità presenta una serie di opzioni specifiche per modificare il comportamento del programma

Carta dei rilievi ombreggiati


Sintassi generale del comando:

gdaldem hillshade srcfile dstfile


Opzioni:

-zfattore di scala per le altitudini (asse z)
-srapporto tra unità verticali e orizzontali (utile se le latitudini sono espresse in gradi e le longitudini in metri)
-azazimut della fonte luce (in gradi). Il valore 0 corrisponde al lato superiore del raster, 90 alla direzione EST
-altaltitudine della fonte luce, in gradi

Carta dei rilievi a colori


Sintassi generale del comando:

gdaldem color-relief srcfile colorsfile dstfile


Color file:

La sintassi del file di configurazione dei colori è estremamente semplice, e mutuata dal comando grass r.colors


Si tratta di un file di testo in cui ciascuna riga è costituita da quattro valori separati da spazio: (altitudine R G B). Un quinto valore può essere aggiunto per caratterizzare la componente di opacità del pixel in uscita


Il valore speciale nv matcha le celle del raster con valore NULL

gdal_contour


L’utility gdal_contour genera un file vettoriale con le curve di livello (isoipse) a partire da un DEM

Nelle ultime versioni di gdal, le curve di livello calcolate sono orientate, per garantire una consistenza topologica


Sintassi base:

gdal_contour [opzioni] srcfile dstfile


Opzioni comuni:

-aimposta il nome della colonna dei valori altimetrici
-fformato di output (di default shp)
-iinposta l'intervallo delle isoipse
-flestrae solo le curve all'altitudine selezionata

gdal_rasterize


Il programma gdal_rasterize permette di rasterizzare un file vettoriale. È possibile unire più sorgenti dati vettoriali, purché siano espresse nello stesso sistema di riferimento (non è fornita la funzionalità di riproiezione al volo)


Sintassi base:

gdal_rasterize [opzioni] srcfile* dstfile


Opzioni comuni:

-bdefinizione delle bande su cui salvare i valori numerici
-burnvalori fissi da utilizzare nel processo di rasterizzazione
-anome dell'attributo del file vettoriale da cui prendere i valori numerici
-whereclausula WHERE per la selezione delle feature
-sqlstringa SQL per la selezione delle feature
-offormato del file di output
-a_srssovrascrive l'SRS del file in output

ogrinfo (1/2)


L’utility ogrinfo offre funzionalità analoghe a quelle fornite da gdalinfo, ma relative a dati in vari formati vettoriali


Sintassi generale:

ogrinfo [opzioni] srcfile
ogrinfo [opzioni] srcfile layername


In base alla sintassi utilizzata, il comando funziona in due modalità differenti:

  • lanciato con un solo parametro, mostra informazioni sui layer memorizzati nel file vettoriale
  • lanciato con il secondo parametro, mostra i dettagliate del layer specificato

ogrinfo (2/2)


Analogamente a gdalinfo, per conoscere la lista dei formati vettoriali supportati si utilizza l’opzione –formats


Opzioni comuni:

-alnella prima modalità di utilizzo, mostra tutte le geometrie di tutti i layers
-sonella seconda modalità di utilizzo, mostra solo i metadati e non le geometrie del layer
-spatseleziona solo le geometrie che intersecano il rettangolo specificato
-whereseleziona le geometrie da visualizzare in base alla clausola SQL WHERE specificata
-sqlseleziona le geometrie da visualizzare in base ad un'espressione SQL

ogr2ogr (1/2)


Il programma ogr2ogr fornisce le seguenti funzionalità:

  • conversione tra i formati vettoriali supportati
  • riproiezione tra differenti SRS
  • selezione di sottoinsiemi di geometrie e clipping
  • merge di geometrie (append mode)
  • segmentazione (su geometrie non puntiformi)

Sintassi generale:

ogr2ogr [opzioni] dstfile srcfile [opzioni layer] [layer [layer …]]

ogr2ogr (2/2)


Segue una lista delle opzioni di più comune utilizzo:


-fspecifica il formato di output
-a_srssovrascrive l'SRS del file in output
-s_srsspecifica l'SRS del file in input
-t_srsspecifica l'SRS del file in output
-whereseleziona le geometrie in base ad una clausola WHERE SQL
-sqlseleziona le geometrie in base ad una espressione SQL
-clipsrceffettua il clipping delle geometrie da selezionare in base a diversi parametri (ad es. un rettangolo)
-appendaggiunge geometrie ad un layer/datasource già esistente
-segmentizeeffettua la segmentazione delle geometrie rispetto ad una distanza massima

Il formato HDF4 (1/2)


HDF4 è un formato file multiobject per la memorizzazione di dati scientifici


Caratteristiche principali:

  • supporto a dati e metadati comunemente utilizzati in ambiti scientifici
  • modalità di memorizzazione efficiente, orientata all’accesso a grossi dataset
  • indipendenza dalla piattaforma ed estendibilità

I file HDF sono autodescrittivi, ovvero contengono le informazioni necessarie per comprendere la natura dei dati memorizzati e localizzarli all’interno dell’archivio


HDF è attualmente sviluppato e supportato da HDF Group, che rilascia librerie e tool per la manipolazione di tali file (Es. HDFView)

Il formato HDF4 (2/2)


Un file HDF (dataset) può contenere un numero illimitato di data structures, di 5 differenti tipologie:

  • Immagini raster
  • Palette
  • Annotations
  • Vdata (tabelle alfanumeriche)
  • Vgroups (gruppi di data structures)

Negli ultimi anni è stata sviluppata una nuova versione del formato (HDF5), concettualmente simile alla precedente, ma con essa incompatibile


Gdal fornisce un supporto trasparente per gestione dei file hdf (v4/v5), attraverso le librerie sviluppate dall’HDF Group

Reperire dati satellitari


È possibile scaricare a titolo gratuito (previa registrazione e con obbligo di citazione) una serie di dati provenienti dalle reti satellitari NASA, attraverso il portale del centro di elaborazione dell’EROS (Earth Resources Observation and Science) http://eros.usgs.gov



Alcuni dei dati scaricabili:

  • Modelli digitali del terreno (DEMs, SRTM, NED)
  • Dati provenienti dai sensori MODIS (satelliti Terra e Acqua)
  • Dati provenienti dai sensori ASTER (piattaforma Terra)
  • Mappe tematiche (Landsat ETM+)

Gestire dati MODIS con GDAL (1/4)


Nel seguito viene presentato un caso di studio, relativo all’acquisizione di dati raster da un dataset contenente gli LST (Land Surface Temperature), rilevati dai sensori MODIS a bordo dei satelliti Terra e Acqua


La NASA distribuisce dati LST preelaborati (a diversi livelli), in differenti risoluzioni temporali (giornalieri, settimanali..) e spaziali (da 250m a 1Km)


Il formato per la distribuzione è l’HDF4; ciascun dataset contiene numerose data structures, che necessitano di un processo di estrazione mediante i tool gdal


Maggiori informazioni sono reperibili all’URL: https://lpdaac.usgs.gov/lpdaac/products

Gestire dati MODIS con GDAL (2/4)


Convenzione sui nomi:

MOD11A1.A2008061.h19v04.005.2008062195705.hdf


  • MOD11A1: nome breve del tipo di prodotto
  • A2008061: data (gregoriana, A-YYYYDDD) di rilevazione
  • h19v04: identificativo della tile (nel sistema di proiezione utilizzato)
  • 005: versione della collezione dati
  • 2008062195705: data (YYYYDDD-HHMMSS) di produzione

Gestire dati MODIS con GDAL (3/4)


MODIS tiling system (proiezione sinusoidale)


Gestire dati MODIS con GDAL (4/4)


Processo di estrazione con gdal:


  • Ottenere informazioni sui dataset e le data structures: gdalinfo
  • Estrarre la data structure di interesse: gdal_translate
  • Riproiettare il raster estratto nell’SRS utilizzato: gdalwarp
  • Ottenere informazioni sui dati contenuti nel raster: gdalinfo
  • Postelaborazione (GRASS map calculator, R+perl/python, …)

Copyright 2010 - Alca Soc. Coop.


http://learn.alcacoop.it - learn@alcacoop.it



released under CreativeCommons 2.5 by-nc-sa