TOC > Aperçu de l'utilisation > Tutoriel GDAL en ligne de commande
Tutoriel GDAL en ligne de commande avec des données météorologiques¶
Introduction¶
Les données de GeoMet du SMC et du Datamart du SMC peuvent être manipulées via la ligne de commande utilisant GDAL, un bibliothèque logicielle très connue utilisée pour lire et écrire des données géospatiales matricielles et vectorielles. Dans les exemples qui suivent, vous utiliserez un fichier GeoTIFF récupéré via une requête Web Coverage Service (WCS) sur GeoMet du SMC. Ce tutoriel vous montrera comment:
- Afficher la version de GDAL installée sur votre système
- Enregistrer la sortie d'une requête WCS sur votre ordinateur
- Lister les informations/metadonnées reliées au fichier matriciel
- Reprojeter un fichier matriciel
- Convertir un fichier GeoTIFF au format de fichier NetCDF
- Obtenir la valeur d'un point spécifique basé sur un emplacement en longitude/latitude
Il existe de multiples façons d'installer GDAL, veuillez vous référez au site web https://gdal.org/ pour plus d'informations.
Pour exécuter les exemples GDAL de ligne de commande, vous devez avoir une connaissance de base de l'utilisation du terminal de ligne de commande. Ces exemples fonctionnent dans un terminal bash.
La version interactive de ce notebook Jupyter est disponible.
Afficher la version de GDAL¶
GDAL est une suite de nombreux outils de ligne de commande. Quand vous installez GDAL, vous obtenez tous les différents outils de ligne de commande. L'outil de base est gdalinfo, qui peut être utilisé pour obtenir des informations relatives à votre installation de GDAL et afficher des informations sur une fichier matriciel.
%%bash
gdalinfo --version
GDAL 3.1.1, released 2020/06/22
Sauvegarder une sortie de requête WCS sur le disque¶
Les requêtes de l'OGC Web Coverage Service permettent au client de récupérer les informations de couverture d'un fichier matriciel pour une zone d'intérêt donnée. Les requêtes WCS sont effectuées sur Internet (HTTPS) et donnent à l'utilisateur plus de flexibilité lorsqu'il demande des informations sur la couverture d'une couche par rapport à la méthode plus traditionnelle de télécharger des fichiers statiques. Le Web Coverage Service permet de nombreux types de requêtes, chacune d'entre elle étant décrite en plus ample détails ci-dessous. Pour plus d'informations sur les paramètres des requêtes WCS, veuillez vous référer à la page WCS GetCoverage de GeoMet du SMC.
Nous allons utiliser une commande curl pour enregistrer le résultat de la requête WCS sur le disque, le fichier sera nommé CMC_glb_TMP.tif. Le résultat est un fichier GeoTIFF, montrant la température (°C) d'un sous-ensemble du Système Global de Prévision Déterministe (SGPD) du SMC.
%%bash
curl "https://geo.weather.gc.ca/geomet?SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=GDPS.ETA_TT&SUBSETTINGCRS=EPSG:4326&SUBSET=x(-120,-85)&SUBSET=y(48,66)&RESOLUTION=x(0.24)&RESOLUTION=y(0.24)&FORMAT=image/tiff" > CMC_glb_TMP.tif
% Total % Received % Xferd Average Speed Time Time Time Current
Dload Upload Total Spent Left Speed
100 44266 100 44266 0 0 78208 0 --:--:-- --:--:-- --:--:-- 78208
Lister les informations à propos d'un fichier matriciel¶
L'outil gdalinfo peut être utilisé pour récupérer les métadonnées du fichier matriciel téléchargé. La sortie de la commande va lister certaines informations sur le fichier, telles que:
- Format du fichier
- Taille du fichier
- Système de coordonnées
- Métadonnées
- Informations sur la bande
%%bash
gdalinfo CMC_glb_TMP.tif
Driver: GTiff/GeoTIFF
Files: CMC_glb_TMP.tif
Size is 146, 75
Coordinate System is:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-119.999862068965513,66.000000000000000)
Pixel Size = (0.239724137931034,-0.240000000000000)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
TIFFTAG_XRESOLUTION=72
TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-119.9998621, 66.0000000) (119d59'59.50"W, 66d 0' 0.00"N)
Lower Left (-119.9998621, 48.0000000) (119d59'59.50"W, 48d 0' 0.00"N)
Upper Right ( -85.0001379, 66.0000000) ( 85d 0' 0.50"W, 66d 0' 0.00"N)
Lower Right ( -85.0001379, 48.0000000) ( 85d 0' 0.50"W, 48d 0' 0.00"N)
Center (-102.5000000, 57.0000000) (102d30' 0.00"W, 57d 0' 0.00"N)
Band 1 Block=146x14 Type=Float32, ColorInterp=Gray
Il est aussi possible d'utiliser gdalinfo pour obtenir quelques statistiques de base sur le fichier matriciel, comme la valeur minimale et maximale en ajoutant l'option -mm. Notez que les valeurs résultantes sont en °C.
%%bash
gdalinfo -mm CMC_glb_TMP.tif | grep Min/Max
Computed Min/Max=0.728,27.828
L'ajout de l'option -proj4 à gdalinfo va produire la définition de la projection en tant qu'un chaîne de caractères proj4:
%%bash
gdalinfo -proj4 CMC_glb_TMP.tif | grep PROJ.4 -A 1
PROJ.4 string is: '+proj=longlat +datum=WGS84 +no_defs'
Reprojeter un fichier raster¶
En utilisant la commande gdalwarp, nous pouvons reprojeter un fichier matriciel. En utilisant les données de GeoMet du SMC et du Datamart du SMC, vous devez seulement fournir une définition de projection de sortie correspondant à un code EPSG, ou vous pouvez utiliser une chaîne de caractères proj4.
L'exemple suivant reprojette le fichier GeoTIFF dans une projection EPSG:3857. Le fichier de sortie est nommé CMC_glb_TMP_epsg3857.tif
%%bash
gdalwarp -t_srs EPSG:3857 CMC_glb_TMP.tif CMC_glb_TMP_epsg3857.tif
Creating output file that is 118P x 114L. Processing CMC_glb_TMP.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.
Ensuite, nous pouvons utiliser gdalinfo pour regarder les coordonnées et la chaîne de caractères proj4 pour s'assurer que la projection de CMC_glb_TMP_epsg3857.tif est réellement différente du fichier original.
%%bash
gdalinfo -proj4 epsg3857.tif | grep -E '(PROJ.4|Corner Coordinates:)' -A 5
PROJ.4 string is: '+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs' Origin = (-20037500.189506221562624,17439592.350557137280703) Pixel Size = (19796.259162699658191,-19796.259162699658191) Metadata: AREA_OR_POINT=Area -- Corner Coordinates: Upper Left (-20037500.190,17439592.351) (179d59'59.74"W, 82d34' 7.50"N) Lower Left (-20037500.190,-17441416.294) (179d59'59.74"W, 82d34'15.13"S) Upper Right (20030128.356,17439592.351) (179d56' 1.34"E, 82d34' 7.50"N) Lower Right (20030128.356,-17441416.294) (179d56' 1.34"E, 82d34'15.13"S) Center ( -3685.917, -911.972) ( 0d 1'59.20"W, 0d 0'29.49"S)
Convertir un fichier GeoTIFF au format de fichier NetCDF¶
En utilisant la commande gdal_translate, nous pouvons convertir un fichier matriciel de n'importe quel format supporté (gdalinfo --formats) en un autre format de fichier matriciel
Dans cet exemple, nous convertissons notre fichier GeoTIFF en un fichier NetCDF. Le paramètre -of NetCDF indique à gdal_translate dans quel format faire la projection. Le fichier de sortie sera nommé CMC_glb_TMP.nc
%%bash
gdal_translate -of NetCDF CMC_glb_TMP.tif CMC_glb_TMP.nc
Input file size is 146, 75 0...10...20...30...40...50...60...70...80...90...100 - done.
Ensuite, en utilisant gdalinfo nous pouvons nous assurer que le fichier NetCDF en sortie est un fichier matriciel valide.
%%bash
gdalinfo CMC_glb_TMP.nc
Driver: netCDF/Network Common Data Format
Files: CMC_glb_TMP.nc
Size is 146, 75
Coordinate System is:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-119.999862068965513,66.000000000000000)
Pixel Size = (0.239724137931034,-0.240000000000000)
Metadata:
Band1#grid_mapping=crs
Band1#long_name=GDAL Band Number 1
Band1#_FillValue=9.96921e+36
crs#GeoTransform=-119.9998620689655 0.2397241379310344 0 66 0 -0.24
crs#grid_mapping_name=latitude_longitude
crs#inverse_flattening=298.257223563
crs#longitude_of_prime_meridian=0
crs#long_name=CRS definition
crs#semi_major_axis=6378137
crs#spatial_ref=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
lat#long_name=latitude
lat#standard_name=latitude
lat#units=degrees_north
lon#long_name=longitude
lon#standard_name=longitude
lon#units=degrees_east
NC_GLOBAL#Conventions=CF-1.5
NC_GLOBAL#GDAL=GDAL 3.1.1, released 2020/06/22
NC_GLOBAL#GDAL_AREA_OR_POINT=Area
NC_GLOBAL#GDAL_TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
NC_GLOBAL#GDAL_TIFFTAG_XRESOLUTION=72
NC_GLOBAL#GDAL_TIFFTAG_YRESOLUTION=72
NC_GLOBAL#history=Tue Jul 21 18:51:21 2020: GDAL CreateCopy( CMC_glb_TMP.nc, ... )
Corner Coordinates:
Upper Left (-119.9998621, 66.0000000) (119d59'59.50"W, 66d 0' 0.00"N)
Lower Left (-119.9998621, 48.0000000) (119d59'59.50"W, 48d 0' 0.00"N)
Upper Right ( -85.0001379, 66.0000000) ( 85d 0' 0.50"W, 66d 0' 0.00"N)
Lower Right ( -85.0001379, 48.0000000) ( 85d 0' 0.50"W, 48d 0' 0.00"N)
Center (-102.5000000, 57.0000000) (102d30' 0.00"W, 57d 0' 0.00"N)
Band 1 Block=146x1 Type=Float32, ColorInterp=Undefined
NoData Value=9.96920996838686905e+36
Metadata:
grid_mapping=crs
long_name=GDAL Band Number 1
NETCDF_VARNAME=Band1
_FillValue=9.96921e+36
Obtenir la valeur pour un point spécifique basé sur un lieu en longitude/latitude¶
En utilisant la commande gdallocationinfo nous pouvons obtenir la valeur brute d'un pixel en spécifiant un emplacement en longitude/latitude ou en spécifiant une position de pixel.
Dans l'exemple suivant, nous utilisons la longitude/latitude. La valeur résultante est en °C.
%%bash
gdallocationinfo -wgs84 CMC_glb_TMP.tif -100 50
Report:
Location: (83P,66L)
Band 1:
Value: 16.4780216217041