Tutoriel: Introduction à GDAL

Site:	IHE DELFT OPENCOURSEWARE
Cours:	FAO CB4WA: Python pour les applications hydrologiques géospatiales
Livre:	Tutoriel: Introduction à GDAL

Imprimé par:	Guest user
Date:	samedi 27 juillet 2024, 14:45

Table des matières

1. Introduction
2. Récupération d'informations de données SIG
- 2.1. Récupérer les informations de données raster
- 2.2. Récupérer les informations de données vectorielles
3. Reprojection de données SIG
- 3.1. Reprojecter des données raster
- 3.2. Reprojecter des données vectorielles
4. Convertir des formats SIG
5. Conversion par lots
6. Conclusion

1. Introduction

Qu'est-ce que GDAL

GDAL est une bibliothèque de traduction pour les formats de données géospatiales raster. Elle est publiée sous une licence Open Source de type X/MIT par l'Open Source Geospatial Foundation. En tant que bibliothèque, elle présente un modèle d'abstraction unique de données pour l'application concernée, et pour tous les formats pris en charge. Elle est également fournie avec une multitude d'utilitaires de ligne de commande utiles pour la traduction et le traitement des données. Pour plus d'informations, voir: http://www.gdal.org

Objectifs d'apprentissage

Après ce cours, vous serez en mesure de:

Récupérer des informations de données SIG
Reprojeter des données SIG
Modifier les propriétés d'un raster
Convertir des formats raster
Convertir des formats vectoriels
Appliquer des requêtes spatiales sur des vecteurs
Convertir des fichiers de valeurs séparées par des virgules
Effectuer des conversions par lots

Logiciel

Pour ces exercices, GDAL doit être installé, de préférence en utilisant le paquet de distribution OSGEO4W. Si vous avez installé QGIS, la distribution OSGEO4W est déjà présente. Si vous devez télécharger QGIS, accédez à la page de téléchargement de QGIS et sélectionnez la version à long terme (LTR) adaptée à votre ordinateur.

Données d'exercice

Les données de l'exercice peuvent être téléchargées à partir de la page principale. Elles contiennent:

roads.shp: carte routière provenant de OpenStreetMap (http://openstreetmap.org)

srtm_37_02.tif: feuille d'un Modèle Numérique d'Elévation (DEM) de Shuttle Radar Topography Mission (SRTM) (http://www2.jpl.nasa.gov/srtm/)

gem_2011_gn1.shp: frontières des communes néerlandaises, gratuitement disponibles à CBS (Statistics Netherlands) et Kadaster (http://www.cbs.nl).

Locations.csv: tableau avec des emplacements d'objets.

landuse.zip: contient des séries temporelles sur l'utilisation des terres au format IDRISI.

2. Récupération d'informations de données SIG

Dans ce chapitre, vous apprendrez à:

récupérer les informations de données raster
récupérer les informations de données vectorielles

Nous présenterons les commandes gdalinfo et ogrinfo.

2.1. Récupérer les informations de données raster

Une des commandes les plus simples et les plus utiles de GDAL est gdalinfo. Lorsqu'on lui donne une image comme argument, elle récupère et imprime toutes les informations pertinentes connues sur le fichier. Ceci est particulièrement utile si l'image contient des données supplémentaires, comme c'est le cas avec les fichiers TIF. Lorsque vous travaillez avec des images satellites, c'est un moyen extrêmement utile de garder la trace de l'emplacement des images en coordonnées long/lat ainsi que la projection de l'image.

1. Ouvrez l'invite où vous avez configuré GDAL. Si vous avez installé QGIS, vous pouvez utiliser OSGeo4W Shell

2. Avec l'invite de commande, allez dans le répertoire où vous avez sauvegardé les données d'exercice, par exemple Z:\gdal_exercises

3. Exécutez la commande suivante:

gdalinfo srtm_37_02.tif <ENTER>

4. Essayez de répondre aux questions suivantes en vous basant sur les informations affichées à l'écran:

Quelle est la taille de l'image?
Quel est le système de coordonnées?
Quel est le code EPSG?
Quelles sont les valeurs minimales et maximales?

Vous trouverez plus d'informations sur l'utilisation de la commande gdalinfo à l'adresse https://gdal.org/programs/gdalinfo.html

Les codes EPSG sont utilisés pour définir les projections. Plus d'informations dans la vidéo ci-dessous.

2.2. Récupérer les informations de données vectorielles

Pour récupérer les informations de données vectorielles, nous utilisons la commande ogrinfo .

1. Exécutez les commandes suivantes:
ogrinfo -al roads.shp | more <ENTER>
ogrinfo -al gem_2011_gn1.shp | more <ENTER>

2. Répondez aux questions suivantes:

Quelles sont les projections de ces deux fichiersde forme (shapefile)?
Quels sont les codes EPSG?

Comme vous pouvez le constater, plusieurs codes EPSG sont indiqués. Cela est dû au fait que différents éléments de la projection ont leur propre code EPSG. La règle générale est que vous pouvez prendre le dernier code EPSG. Ainsi, la projection de roads.shp est EPSG: 4326, ce qui signifie qu'il est en système de coordonnées géographiques (GCS) avec des coordonnées en latitude/longitude.

Pour gem_2011_gn1.shp la projection est EPSG: 28992, qui est la projection néerlandaise (Amersfoort RD/New).

Pour plus d'informations sur la commande ogrinfo, cliquez ici.

3. Reprojection de données SIG

Dans ce chapitre, vous apprendrez à:

Reprojeter des données ratser
Reprojeter des données vectorielles

Vous utiliserez les commandes gdalwarp et ogr2ogr.

Dans cet exercice, nous voulons faire une carte de la commune de Delft avec les principales routes et le relief. Comme les données sont dans des formats différents, nous devons les reprojeter dans un système de coordonnées commun. Ici, nous reprojetons tous les ensembles de données à la projection néerlandaise Amersfoort/RD New.

3.1. Reprojecter des données raster

GDAL permet de changer le système de coordonnées d'un raster en utilisant la syntaxe suivante:

gdalwarp -t_srs EPSG:... <entrée> <sortie>

L'argument -t_srs spécifie le système de coordonnées cible. Si le système de coordonnées source est inconnu, il doit être spécifié avec l'argument -s_srs . EPSG:... spécifie le code EPSG de la projection. <entrée> et <sortie> sont respectivement les données d'entrée et de sortie.

Nous allons maintenant reprojeter un Modèle Numérique d'Elévation (DEM) acquis par Shuttle Radar Topography Mission (SRTM). Vous pouvez télécharger des DEM pour votre propre zone d'intérêt à partir de USGS Earth Explorer. Ici, nous utiliserons les données du cours.

Afin de reprojeter le MNE de WGS-84 lat/lon à Amersfoort/RD New nous utilisons cette commande:

gdalwarp -t_srs EPSG:XXXXX srtm_37_02.tif dem_rd.tif

1. Remplacez XXXXX par le code EPSG approprié pour Amersfoort/RD New (voir une de vos réponses précédentes en utilisant ogrinfo).

2. Exécutez la commande:

gdalwarp -t_srs EPSG:28992 srtm_37_02.tif dem_rd.tif <ENTER>

3. Visualisez le résultat dans QGIS.

Pour plus d'informations sur la commande gdalwarp, cliquez ici.

3.2. Reprojecter des données vectorielles

Pour des données vectorielles, on utilise la commande ogr2ogr. La syntaxe est la suivante:

ogr2ogr -t_srs EPSG:XXXXX <sortie> <entrée>

Remarquez que <sortie> et <entrée> sont dans un ordre différent comparé à gdalwarp!

Ici, nous allons convertir les données routières d'OpenStreetMaps en projection Amersfoort/RD New.

1. Exécutez:
ogr2ogr -t_srs EPSG:28992 roadsreprojected.shp roads.shp

La commande s’ exécute (elle prend un certain temps), mais donne un avertissement.

2. Que signifie cet avertissement?

Vous trouverez plus d'informations sur la commande ogr2ogr ici.

4. Convertir des formats SIG

Dans ce chapitre, vous apprendrez à:

Convertir différents formats raster
Convertir différents formats vectoriels

Vous utiliserez les commandes gdal_translate et ogr2ogr.

4.1. Convertir des formats raster

La fonction principale de gdal_translate est de changer de format d'image. La syntaxe de base est la suivante:

gdal_translate -of FORMAT fichierentrée fichiersortie

Tous les formats pris en charge peuvent être trouvés ici.

Maintenant nous allons convertir le DEM de geoTiff au format SAGA. SAGA est un SIG qui a son propre format binaire supporté par gdal avec l'extension .sdat.

1. Exécutez la commande suivante:
gdal_translate -of SAGA dem_rd.tif dem_rd.sdat

Certains formats nécessitent plus d'arguments. PCRaster, par exemple, nécessite une spécification du type de données (booléen, nominal, scalaire, etc.). Convertissons les mêmes données au format PCRaster.

2. Exécutez la commande suivante:
gdal_translate -of PCRaster -ot Float32 -mo "VS_SCALAR" dem_rd.tif dem_rd.map

Pour convertir au format PCRaster, vous devez connaître le type de données du raster. Dans l'exemple ci-dessus, le DEM (MNE) est une donnée continue, le type de données est donc scalaire. Si la couche est discrète (classes), le type de données est nominal.

Type de données	-of	-mo
Booléen	Byte	"VS_BOOLEAN"
Nominal	Int32	"VS_NOMINAL"
Ordinal	Int32	"VS_ORDINAL"
Scalaire	Float32 or Float64	"VS_SCALAR"
Direction	Float32 or Float64	"VS_DIRECTION"
LDD	Int32	"VS_LDD"

4.2. Convertir des formats vectoriels

GDAL permet des conversions entre formats vectoriels en utilisant la commande ogr2ogr.

Trouvez ici la liste des formats pris en charge.

La syntaxe générale est:

ogr2ogr -f <"Format"> <sortie> <entrée>

Le "Format" est le nom abrégé dans le tableau des pilotes des vecteurs et doit être donné entre guillemets.

Nous allons convertir gem_2011_gn1.shp en un fichier Google KML qui peut être ouvert dans Google Earth.

1. Exécutez la commande suivante:
ogr2ogr -f "KML" gem.kml gem_2011_gn1.shp

Ignorez les avertissements. Ils sont liés à des formats et des caractères qui ne sont pas pris en charge.

Vous trouverez ici toutes les options de la commande ogr2ogr. Notez que vous pouvez également reprojeter et convertir en même temps en utilisant -t_srs. Il y a beaucoup d'autres choses que vous pouvez faire, mais nous ne les abordons pas dans ce tutoriel de base.

2. Ouvrez le fichier dans Google Earth pour vérifier le résultat.

4.3. Requêtes spatiales sur des données vectorielles

Pour notre carte de Delft, nous voulons faire l'analyse SIG suivante:

Sélectionner la commune de Delft sur la carte des communes et la sauvegarder dans un nouveau fichier de forme (shapefile);
Intersecter les limites de la commune de Delft avec la carte routière des Pays-Bas.

Nous pouvons utiliser une requête spatiale pour sélectionner une caractéristique dans une carte vectorielle.

1. Quel est l'attribut de la carte des communes contenant les noms des communes ? Vous pouvez utiliser soit ogrinfo ou QGIS pour répondre à cette question.

2. Exécutez la commande suivante:

ogr2ogr
-f "ESRI Shapefile" -where GM_NAAM='Delft' 
-a_srs EPSG:28992 delft.shp gem_2011_gn1.shp

Ceci sauvegardera l'entité avec GM_NAAM Delft dans un nouveau fichier appelé delft.shp. L'argument -a_srs EPSG:28992 est utilisé pour affecter la projection Amersfoort/RD New au fichier de sortie. L'argument -f définit le format de sortie.

3. Nous pouvons découper le vecteur roadsreprojected.shp afin qu'il ne couvre que la municipalité de Delft. Nous utilisons pour cela l'argument -clipsrc.

ogr2ogr -clipsrc delft.shp roadsdelft.shp roadsreprojected.shp

4. Ouvrez maintenant dans QGIS le DEM (MNE) reprojeté (dem_rd.tif), la carte routière découpée (roadsdelft.shp) et la commune de Delft (delft.shp).

4.4. Convertir des valeurs séparées par des virgules (CSV)

Parfois, vous souhaitez reprojeter des coordonnées se trouvant dans un fichier ASCII, par exemple, d’un fichier qui a été enregistré en texte dans un tableur. Ici, nous allons convertir les coordonnées d'un fichier ASCII séparé par des virgules, locations.csv,en un nouveau fichier ASCII locations_reprojected.csv.

Il est recommandé de visualiser d'abord le contenu d'un fichier CSV et de vérifier (1) les coordonnées, et (2) le séparateur de colonnes. Le séparateur n'est pas toujours une virgule et dépend parfois des paramètres de langues utilisés lors de l'exportation à partir du tableur.

1. Visualisez le contenu du fichier locations.csv . Vous pouvez utiliser la commande type à partir de l'invite de commande comme vous l'avez appris dans le tutoriel Ligne de commande.

Nous pouvons voir que les coordonnées sont en lat/lon, ce qui signifie que nous pouvons utiliser EPSG:4326. Le séparateur de colonne est une virgule. La dernière colonne donne les objets avec une chaîne de caractères entre guillemets.

Pour changer la projection du CSV, nous devons d'abord créer une source de données virtuelle en créant un fichier de contrôle XML.

2. Sur votre ordinateur Windows, ouvrez Notepad et tapez/copiez le code XML indiqué ci-dessous. Utilisez des marges de trois espaces.

<OGRVRTDataSource>
   <OGRVRTLayer name="locations">
      <SrcDataSource>locations.csv</SrcDataSource>
      <GeometryType>wkbPoint</GeometryType>
      <LayerSRS>EPSG:4326</LayerSRS>
      <GeometryField encoding="PointFromColumns" x="lon" y="lat"/>
   </OGRVRTLayer>
</OGRVRTDataSource>

3. Enregistrez le fichier sous le nom de locations.vrt dans le dossier gdal_exercises.

Quelques explications à propos du fichier XML:

<OGRVRTLayer name="locations"> doit correspondre à <SrcDataSource>locations.csv</SrcDataSource>
<LayerSRS>EPSG:4326</LayerSRS> doit correspondre au code EPSG des colonnes de coordonnées
<GeometryField encoding="PointFromColumns" x="lon" y="lat"/> indique les colonnes avec les coordonnées que vous voulez convertir.

4. Exécutez la commande suivante:

ogr2ogr
-t_srs EPSG:28992 -f "CSV" locations_reprojected.csv locations.vrt -lco
GEOMETRY=AS_XY

Dans cet exemple locations.csv avec des coordonnées lat/lon WGS-84 est converti en locations_projected.csv avec la projection Amersfoort/RD New.

5. Utilisez Notepad pour vérifier locations_reprojected.csv. Qu'est-ce qui est enregistré dans chaque colonne?

De la même manière, nous pouvons convertir le fichier séparé par des virgules en un fichier shapefile.

6. Exécutez la commande suivante:

ogr2ogr
-f "ESRI Shapefile" -t_srs
EPSG:28992 locations.shp locations.vrt

7. Vérifiez le résultat avec ogrinfo.

8. Visualiser le shapefile dans QGIS en traçant les emplacements sur le DEM (MNE), la carte routière et la frontière de la commune de Delft. Faites une belle carte.

9. Convertissez le fichier de localisation en Google KML, ouvrez-le dans Google Earth et découvrez ce à quoi correspond l'emplacement des objets.

5. Conversion par lots

Les programmes SIG de bureau sont très utiles pour des opérations SIG, mais sont difficiles à utiliser si nous devons répéter la même tâche pour de nombreuses couches SIG. Le scriptage peut alors être une solution.

Nous avons ici un exemple d’ensemble de données provenant d'un modèle d'occupation des terres de Dublin. Les données sont au format raster IDRISI (.rst), avec une couche pour chaque année entre 1990 et 2030. Notre tâche consiste à convertir toutes les couches au format GeoTiff (.tif ).

1. Dézippez le fichier landuse.zip fourni avec les données du cours dans un nouveau dossier et vérifiez le contenu.

Rappelez-vous que pour la conversion d'un fichier raster, nous utilisons ceci:
gdal_translate -of GTiff 01_State19900101.rst 01_State19900101.tif

Maintenant nous allons créer un fichier batch (voir le tutoriel Ligne de Commande) qui inclut une boucle pour convertir tous les fichiers du dossier.

2. Ouvrez un éditeur de texte, par exemple Notepad

3. Ajoutez le code:

for %%f in (*.rst) do (


       echo %%~nf
       gdal_translate -of GTiff %%f %%~nf.tif

)

4. Enregistrez le fichier batch sous le nom rst2tif.bat dans le dossier contenant les rasters d'utilisation des terres (n'oubliez pas de changer l'extension si vous utilisez Notepad, erreur classique !).

Essayez de comprendre le code. C'est une boucle qui passe en revue tous les fichiers *.rst se trouvant dans le dossier. %%f est la variable qui contient le nom de chaque fichier. Avec echo nous pouvons afficher quelque chose à l'écran. Ici, nous affichons %%~nf ,qui est la partie du nom de fichier avant le point qui le sépare de l'extension. Ensuite nous utilisons la commande gdal_translate avec le format de sortie GeoTiff. A la fin de la ligne, nous ajoutons l'extension .tif au nom de fichier.

5.Exécutez le fichier batch. Tapez
rst2tif <ENTER>

6. Vérifiez les résultats

6. Conclusion

Dans cette leçon, vous avez appris à:

Récupérer des informations à partir de données SIG
Reprojeter des données SIG
Changer les propriétés d’un raster
Convertir des formats de raster
Convertir des formats vectoriels
Appliquer des requêtes spatiales à des vecteurs
Convertir des fichiers de valeurs séparées par des virgules
Effectuer une conversion par lots

Cette vidéo montre les résultats du tutoriel GDAL: