../../_images/logo_gdal.png ../../_images/OSGeo_project.png

Guide de démarrage pour GDAL/OGR

Ce guide de démarrage rapide est divisé en deux parties : GDAL (données raster) et OGR (données vectorielles).

Ce guide décrit comment :

GDAL

  • explorer vos données raster avec gdalinfo

  • convertir dans d’autres formats avec gdal_translate

  • reprojeter vos données avec gdalarp

  • créer une mosaïque à partir de vos données avec gdal_warp ou gdal_merge.py

  • créer un Shapefile d’index de tuiles raster avec gdalindex

OGR

  • obtenir des informations sur vos données avec ogrinfo

  • utiliser ogr2ogr pour convertir vos données dans d’autres formats

La seule chose dont vous aurez besoin pour ce guide de démarrage rapide est un terminal. Si vous souhaitez visualiser les résultats, vous pouvez utiliser l’une des applications SIG de bureau sur OSGeoLive comme QGIS.

Faire connaissance avec GDAL

Vous trouverez les données de démonstration dans /usr/local/share/data. Nous allons explorer les données Natural Earth dans ce guide de démarrage. Nous allons travailler avec une copie des données donc la première étape est de les copier dans votre répertoire utilisateur.

cd /home/user
cp -R /usr/local/share/data/natural_earth2/ ./gdal_natural_earth

Vous trouverez un fichier raster Natural Earth ainsi qu’un World-file .tfw ici

ls /home/user/gdal_natural_earth/HYP_50M_SR_W.*

Astuce

Ouvrez le fichier avec un logiciel SIG tel que QGIS et jetez-y un œil.

Obtenir des informations sur le raster avec gdalinfo

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
Notez :

  • le format est « GTiff/GeoTIFF »

  • la taille est de 10800x5400

  • les 3 bandes de type Byte

  • les coordonnées

  • l’absence de système de coordonnées

Changement de format de fichier simple

Premièrement, prenez connaissance des formats. L’option –formats de la commande gdal_translate peut être utilisée pour voir la liste des formats disponibles.

Chaque format indique s’il est :

  • en lecture seule (ro),

  • en lecture/écriture (rw),

  • en lecture/écriture/mise à jour (rw+).

gdal_translate --formats

L’option –format permet d’obtenir des détails sur un format en particulier dont les options de création et les types de données autorisés.

gdalinfo --format jpeg
gdal_translate --format png

Conversion

Les conversions sont réalisées avec la commande gdal_translate. Le format de sortie par défaut est le GeoTIFF. L’option -of est utilisée pour choisir un autre format de sortie et -co pour spécifier une option de création :

gdal_translate -of JPEG -co QUALITY=40 HYP_50M_SR_W.tif HYP_50M_SR_W.jpg

L’option -ot peut être utilisée pour modifier le type de données en sortie.

gdal_translate -ot Int16 HYP_50M_SR_W.tif HYP_50M_SR_W_Int16.tif

Utilisez gdalinfo pour vérifier le type de données.

Ré-échantillonner

L’option -outsize peut être utiliser pour spécifier la taille du raster en sortie.

gdal_translate -outsize 50% 50% HYP_50M_SR_W.tif  HYP_50M_SR_W_small.tif

Utilisez gdalinfo pour vérifier la taille.

L’option -scale peut être utilisée pour ré-échantillonner un raster. Un contrôle explicite des valeurs d’entrée et de sortie est également possible. L’option -mm de gdalinfo permet d’obtenir les valeurs min/max des pixels.

Découpons notre raster en deux avec l’option -srcwin qui permet de faire une copie en se basant sur la localisation des pixels/lignes (décalage en x, décalage en y, nombre de pixels en x, nombre de pixels en y). Vous pouvez également utiliser -projwin et définir les coordonnées des coins dans le système de coordonnées du raster (coordonnée en x du pixel supérieur gauche puis en y, coordonnée en x du pixel inférieur droit puis en y).

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

Index de tuiles raster avec gdalindex

Vous pouvez générer un Shapefile d’index des tuiles raster. Pour chaque raster, un polygone est créé selon son emprise ainsi que le chemin vers le fichier.

gdaltindex index_natural_earth.shp *st.tif

Jetez un œil au Shapefile en sortie avec QGIS ainsi qu’avec ogrinfo (plus d’informations sur ogrinfo suivent dans ce tutoriel).

../../_images/gdal_gdaltindex.png 
ogrinfo index_natural_earth.shp index_natural_earth
INFO: Open of `index_natural_earth.shp'
    using driver `ESRI Shapefile' successful.

Layer name: index_natural_earth
Geometry: Polygon
Feature Count: 2
Extent: (-180.000000, -90.000000) - (180.000000, 90.000000)
Layer SRS WKT: (unknown)
location: String (255.0)
OGRFeature(index_natural_earth):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_natural_earth):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))

Reprojection

Pour ce traitement, nous supposons que le raster HYP_50M_SR_W.tif a été correctement créé. Comme nous l’avons vu précédemment avec gdalinfo, aucun système de coordonnées n’a été spécifié. Nous assignons donc premièrement le WGS84 à notre raster.

gdal_translate -a_srs WGS84 HYP_50M_SR_W.tif HYP_50M_SR_W_4326.tif

La commande gdalwarp permet de reprojeter un raster. Ici, nous reprojetons le raster WGS84 en projection Mercator :

gdalwarp -t_srs '+proj=merc +datum=WGS84' HYP_50M_SR_W_4326.tif mercator.tif

Utilisez gdalinfo pour vérifier le changement et jeter un œil au raster.

../../_images/gdal_mercator.png

Ici nous reprojetons dans la projection Ortho.

gdalwarp -t_srs '+proj=ortho +datum=WGS84' HYP_50M_SR_W_4326.tif ortho.tif
../../_images/gdal_ortho.png

Notez comment les pôles sont coupés? C’est parce que les bords aux pôles ne peut pas être reprojetés car gdalwarp ne lit pas toutes les données. Nous pouvons forcer gdalwarp à lire un supplément de données autour des morceaux comme un moyen de résoudre ce problème. En savoir plus à ce sujet dans le tutoriel Raster: https://trac.osgeo.org/gdal/wiki/UserDocs/RasterProcTutorial.

Créer une mosaïque

gdal_merge.py est un script Python qui peut être utiliser pour créer des mosaïques simples. Mosaïquez les raster east.tif et west.tif dans un fichier unique :

gdal_merge.py  east.tif west.tif -o merged.tif

La même tâche peut être accomplie avec gdalwarp. gdlawarp présente de nombreux avantages par rapport à gdal_merge mais peut être lent pour fusionner de multiples fichiers :

gdalwarp east.tif west.tif warpmerged.tif

Faire connaissance avec OGR

cd /home/user/gdal_natural_earth/

Astuce

Ouvrez le Shapefile dans un logiciel SIG tel QGIS et jetez-y un œil.

Obtenir des informations sur la couche vecteur avec ogrinfo

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)

Obtenez un récapitulatif de vos données avec ogrinfo et l’option -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)

Si vous utilisez ogrinfo sans option, vous obtiendrez un récapitulatif de vos données suivi d’une section pour chaque donnée.

ogrinfo -ro ne_10m_admin_0_countries.shp ne_10m_admin_0_countries

Vous pouvez envoyer le résultat de la commande ogrinfo à grep pour le filtrer et n’obtenir que les informations sur le champ 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.

Vous pouvez convertir vos données dans d’autres formats. Obtenez la liste de tous les formats gérés avec –formats.

Utiliser ogr2ogr pour convertir la couche dans d’autres formats

Vous pouvez utiliser ogr2ogr pour convertir vos couches vecteur dans d’autres formats. –formats vous donne la liste de tous les formats gérés avec les possibilités de lecture/écriture.

Convertissez les pays en GML.

ogr2ogr --formats
ogr2ogr -f GML countries.xml ne_10m_admin_0_countries.shp

Choses à essayer

Voici quelques défis supplémentaires que vous pouvez relever :

  • Essayer gdalwarp ou gdal_merge.py pour mosaïquer vos données

  • Essayer gdaladdo pour créer des vignettes intégrées

  • QGIS utilise également GDAL/OGR pour la gestion de nombreux formats. Il propose aussi les outils GdalTools pour traiter les données raster.

  • Essayer ogr2ogr pour importer/exporter vos données vecteur dans d’autres formats comme PostGIS. Jeter un œil aux options proposées par ogr2ogr.

  • Essayer l’extension QGIS OGR-Layer-Konverter.

Ensuite ?

C’est seulement la première étape sur la route de l’utilisation de GDAL et OGR. Il y a beaucoup plus de fonctionnalités que vous pouvez essayer.

Page d’accueil de GDAL

Tutoriel pour GDAL raster

Workshop sur GDAL