Non sarà necessario nient’altro che un terminale per questo tutorial. Se volete visualizzare i risultati, potete usare uno dei Software Desktop GIS presenti in OSGeo-Live come Quantum GIS (QGIS).
Questa guida rapida è divisa in due parti GDAL (dati raster) e OGR (dati vettoriali). Si incomincerà con GDAL.
Questo tutorial descrive come:
Si troveranno i dati di demo in /usr/local/share/data. Si potrebbe volere vedere il documento dati Natural Earth per questo tutorial. Si lavorerà su una copia dei dati. Perciò la prima cosa è copiare i dati nella directory home.
cd /home/user
cp -R /usr/local/share/data/natural_earth2/ ./gdal_natural_earth
Quindi troverete un NaturalEarth Raster file e un .tfw World-file:
ls /home/user/gdal_natural_earth/HYP_50M_SR_W.*
Tip
Oprire il file con un Desktop GIS tipo QGIS. E date un’occhiata al file.
gdalinfo HYP_50M_SR_W.tif
Driver: GTiff/GeoTIFF
Files: HYP_50M_SR_W.tif
HYP_50M_SR_W.tfw
Size is 10800, 5400
Coordinate System is `'
Origin = (-179.999999999999972,90.000000000000000)
Pixel Size = (0.033333333333330,-0.033333333333330)
Metadata:
TIFFTAG_SOFTWARE=Adobe Photoshop CS3 Macintosh
TIFFTAG_DATETIME=2009:09:19 10:13:17
TIFFTAG_XRESOLUTION=342.85699
TIFFTAG_YRESOLUTION=342.85699
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (-180.0000000, 90.0000000)
Lower Left (-180.0000000, -90.0000000)
Upper Right ( 180.0000000, 90.0000000)
Lower Right ( 180.0000000, -90.0000000)
Center ( -0.0000000, 0.0000000)
Band 1 Block=10800x1 Type=Byte, ColorInterp=Red
Band 2 Block=10800x1 Type=Byte, ColorInterp=Green
Band 3 Block=10800x1 Type=Byte, ColorInterp=Blue
Per prima cosa si può ottenere i formati riconosciuti. L’opzione –formats di gdal_translate può essere usate per vedere una lista di formati disponibili.
gdal_translate --formats
L’opzione –format può essere usate per interrogare maggiori informazioni riguardo un particolare formato, incluso le opzioni per la creazione di nuovo dati, e le tipologie di dati permesse.
gdalinfo --format jpeg
gdal_translate --format png
Conversioni sono effettuare con il comando gdal_translate. Il formato di output di default è il GeoTIFF. L’opzione -of è usata per selezionare un formato di output format e l’opzione -co è usate per specificare un’opzione di creazione:
gdal_translate -of JPEG -co QUALITY=40 HYP_50M_SR_W.tif HYP_50M_SR_W.jpg
L’opzione -ot può essere usata per modificare la tipologia del dato di output.
gdal_translate -ot Int16 HYP_50M_SR_W.tif HYP_50M_SR_W_Int16.tif
Usare gdalinfo per verificare la tipologia del dato
L’opzione -outsize può essere usata per impostare la dimensione del file di output.
gdal_translate -outsize 50% 50% HYP_50M_SR_W.tif HYP_50M_SR_W_small.tif
Usare gdalinfo per verificare la dimensione.
L’opzione -scale può essere usata per riscalare i dati. È disponibile controllo esplicito dei range di input e output. L’opzione -mm di gdalinfo può essere usata per vedere i valori minimo e massimo dei pixel.
Si divida l’immagine in due con -srcwin che crea una copia basata sulla posizione pixel/linea (xoff yoff xsize ysize). È possibile anche usare -projwin e definire gli angoli in coordinate georiferite (ulx uly lrx lry).
gdalinfo -mm HYP_50M_SR_W.tif
gdal_translate -srcwin 0 0 5400 5400 HYP_50M_SR_W.tif west.tif
gdal_translate -srcwin 5400 0 5400 5400 HYP_50M_SR_W.tif east.tif
Si può creare uno shapefile come un tileindex di raster. Per ogni immagine un poligono è generato con il perimetro dell’estensione del raster e il percorso al file.
gdaltindex index_natural_earth.shp *st.tif
Controllate lo shapefile di output con QGIS e ogrinfo (imparerete di più su ogrinfo dopo in questo tutorial)
ogrinfo ../HYP_50M_SR_W/ index
INFO: Open of `../HYP_50M_SR_W/'
using driver `ESRI Shapefile' successful.
Layer name: index
Geometry: Polygon
Feature Count: 2
Extent: (-180.000000, -90.000000) - (180.000000, 90.000000)
Layer SRS WKT: (unknown)
location: String (255.0)
OGRFeature(index):0
location (String) = east.tif
POLYGON ((-0.00000000001796 90.0,179.999999999964047 90.0,179.999999999964047 -89.999999999982009,-0.00000000001796 -89.999999999982009,-0.00000000001796 90.0))
OGRFeature(index):1
location (String) = west.tif
POLYGON ((-179.999999999999972 90.0,-0.00000000001796 90.0,-0.00000000001796 -89.999999999982009,-179.999999999999972 -89.999999999982009,-179.999999999999972 90.0))
Per questo processo si assume che HYP_50M_SR_W.tif è stato correttamente creato con i confini. Come si è visto precedentemente con gdainfo nessun sistema di coordinate è stato impostato. Si assegna all’immagine WGS84 come sistema di coordinate nel primo passaggio.
gdal_translate -a_srs WGS84 HYP_50M_SR_W.tif HYP_50M_SR_W_4326.tif
Il comando gdalwarp può essere usato per riproiettare le immagini. Qui si riproietta l’immagine WGS84 in una proiezione di Mercatore:
gdalwarp -t_srs '+proj=merc +datum=WGS84' HYP_50M_SR_W_4326.tif mercator.tif
Usare gdalinfo per verificare il cambio e visualizzate l’immagine.
Qui si riproietta in una proiezione Ortho.
gdalwarp -t_srs '+proj=ortho +datum=WGS84' HYP_50M_SR_W_4326.tif ortho.tif
Notare come i poli sono tagliati? Questo perchè gli angoli al polo non possono essere riproiettati gdalwarp non legge tutti i dati. Come soluzione si può forzare gdalwarp a leggere un mucchio di dati eccedenti. Approfondire questo argomento nel RasterTutorial http://trac.osgeo.org/gdal/wiki/UserDocs/RasterProcTutorial.
gdal_merge.py è un script Python che può essere utile per semplici mosaicature. Mosaicare east.tif e west.tif in un singolo file:
gdal_merge.py east.tif west.tif -o merged.tif
La stessa operazione può essere fatta con gdalwarp. gdalwarp ha una varietà di vantaggi rispetto gdal_merge, ma può essere lento per unire molti file:
gdalwarp east.tif west.tif warpmerged.tif
cd /home/usr/gdal_natural_earth/
Tip
Aprire lo shapefile con un Desktop GIS tipo QGIS. E date un’occhiata al file.
ogrinfo -ro /home/user/gdal_natural_earth
INFO: Open of `/home/user/gdal_natural_earth'
using driver `ESRI Shapefile' successful.
1: ne_10m_populated_places (3D Point)
2: ne_10m_geography_regions_polys (3D Polygon)
3: ne_10m_admin_1_states_provinces_shp (3D Polygon)
4: ne_10m_urban_areas (3D Polygon)
5: ne_10m_geography_marine_polys (3D Polygon)
6: ne_10m_land (3D Polygon)
7: ne_10m_geography_regions_elevation_points (3D Point)
8: ne_10m_admin_0_countries (3D Polygon)
9: ne_10m_rivers_lake_centerlines (3D Line String)
10: ne_10m_lakes (3D Polygon)
11: ne_10m_geography_regions_points (3D Point)
12: ne_10m_ocean (3D Polygon)
Ottenere un riepilogo dei dati con ogrinfo insieme a -so.
ogrinfo -ro -so ne_10m_admin_0_countries.shp ne_10m_admin_0_countries
INFO: Open of `ne_10m_admin_0_countries.shp'
using driver `ESRI Shapefile' successful.
Layer name: ne_10m_admin_0_countries
Geometry: 3D Polygon
Feature Count: 254
Extent: (-180.000000, -90.000000) - (180.000000, 83.634101)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
DATUM["WGS_1984",
SPHEROID["WGS_84",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]]
scalerank: Integer (4.0)
featurecla: String (30.0)
labelrank: Real (16.6)
sovereignt: String (254.0)
sov_a3: String (254.0)
adm0_dif: Real (16.6)
level: Real (16.6)
type: String (254.0)
admin: String (254.0)
adm0_a3: String (254.0)
geou_dif: Real (16.6)
geounit: String (254.0)
gu_a3: String (254.0)
su_dif: Real (16.6)
subunit: String (254.0)
su_a3: String (254.0)
brk_diff: Real (16.6)
name: String (254.0)
name_long: String (254.0)
brk_a3: String (254.0)
brk_name: String (254.0)
brk_group: String (254.0)
abbrev: String (254.0)
postal: String (254.0)
formal_en: String (254.0)
formal_fr: String (254.0)
note_adm0: String (254.0)
note_brk: String (254.0)
name_sort: String (254.0)
name_alt: String (254.0)
mapcolor7: Real (16.6)
mapcolor8: Real (16.6)
mapcolor9: Real (16.6)
mapcolor13: Real (16.6)
pop_est: Real (16.6)
gdp_md_est: Real (16.6)
pop_year: Real (16.6)
lastcensus: Real (16.6)
gdp_year: Real (16.6)
economy: String (254.0)
income_grp: String (254.0)
wikipedia: Real (16.6)
fips_10: String (254.0)
iso_a2: String (254.0)
iso_a3: String (254.0)
iso_n3: String (254.0)
un_a3: String (254.0)
wb_a2: String (254.0)
wb_a3: String (254.0)
woe_id: Real (16.6)
adm0_a3_is: String (254.0)
adm0_a3_us: String (254.0)
adm0_a3_un: Real (16.6)
adm0_a3_wb: Real (16.6)
continent: String (254.0)
region_un: String (254.0)
subregion: String (254.0)
region_wb: String (254.0)
name_len: Real (16.6)
long_len: Real (16.6)
abbrev_len: Real (16.6)
tiny: Real (16.6)
homepart: Real (16.6)
Se si lancia ogrinfo senza un parametro otterrete un riepilogo dei dati e successivamente una sezione per ogni dataset.
ogrinfo -ro ne_10m_admin_0_countries.shp ne_10m_admin_0_countries
È possibile inoltrare il risultato di ogrinfo a grep per filtrare e ottenere solo l’attributo COUNTRY.
ogrinfo ne_10m_admin_0_countries.shp ne_10m_admin_0_countries | grep 'admin '
admin (String) = Aruba
admin (String) = Afghanistan
admin (String) = Angola
admin (String) = Anguilla
admin (String) = Albania
admin (String) = Aland
admin (String) = Andorra
etc.
È possibile convertire i dati in altri formati. Si può avere la lista dei formati supportati con –formats.
Si può usare ogr2ogr per convertire dati con elementi semplici tra diversi formati. Si può avere la lista dei formati supportati con –formats. con i permessi di lettura/scrittura
Convertire gli stati in GML.
ogr2ogr --formats
ogr2ogr -f GML countries.xml ne_10m_admin_0_countries.shp
Qui alcune ulteriori prove da provare:
Questo è solo il primo passaggio sulla strada per usare GDAL e OGR. Ci sono un sacco di ulteriori funzionalità che potete provare.
GDAL Project
Tutto su OGR
GDAL Tutorial