Auteur : Paul Passy
Licence : 
Outils raster
Il est très souvent nécessaire d’appliquer quelques traitements sur les données rasters avant de les interpréter ou pour en tirer des informations supplémentaires. Il existe de nombreux outils de géomatique dédiés à l’exploitation des données rasters. Nous verrons ici quelques outils très utiles et disponibles dans différents logiciels. Nous détaillerons par exemple comment découper ou masquer un raster ou comment appliquer des opérations algébriques sur une ou plusieurs couches rasters.
Table des matières
Histogramme d’un raster
Il est souvent utile de regarder la forme de la distribution des valeurs des pixels d’un raster donné en traçant ce que nous appelons l’histogramme du raster considéré. Sur ce type d’histogramme, il est possible de connaître le nombre de pixels se trouvant dans chaque tranche de valeurs. Ces tranches de valeurs sont en général définies par l’utilisateur au moment de la construction de l’histogramme. La construction de ce genre d’histogramme est possible dans tous les logiciels de géomatique et dans les différents langages de programmation.
Dans les différents exemples ci-après, nous tracerons l’histogramme du MNT du bassin de la Roya (Alpes-Maritime) issu du SRTM et nommé srtm_roya_L93.tif.
Histogramme d’un raster avec QGIS
Version de QGIS : 3.26.1
Tout d’abord nous chargeons le raster à étudier dans QGIS. Il est possible d’avoir un premier aperçu de l’histogramme du raster en allant dans ses Propriétés comme présenté dans la partie dédiée à l’import de raster avec QGIS (Exploration des propriétés d’un raster). Cependant, cet histogramme est assez sommaire et il existe une autre fonctionnalité de QGIS pour tracer un histogramme plus ergonomique. Nous choisissons le menu . La fenêtre suivante s’affiche (Fig. 234).
Fig. 234 Menu de création d’un histogramme raster dans QGIS.
À la ligne Couche source nous indiquons le raster à analyser, ici srtm_roya_L93.tif. Dans le menu Numéro de bande nous indiquons la bande à analyser, dans le cas d’un raster multi-bandes. Pour le menu Nombre de classes, nous indiquons le nombre d’intervalles à considérer dans l’histogramme, nous pouvons mettre 100 ici. Enfin, à la ligne Histogramme nous indiquons un chemin et un nom pour l’histogramme qui sera produit. Remarquons que l’histogramme produit est au format html. En effet, l’histogramme sera tracé en utilisant la technologie Plotly qui permet de créer des graphiques dynamiques. Enfin, nous cliquons sur Exécuter.
À l’exécution rien n’apparaît mais dans l’onglet Journal de la fenêtre précédente nous trouvons le chemin du fichier html produit (Fig. 235).
Fig. 235 Le chemin de l’histogramme html produit, à la ligne Output.
Il suffit alors de copier le chemin, d’ouvrir un navigateur web comme Firefox ou autre et de coller ce chemin dans la barre d’adresse. L’histogramme s’ouvre alors dans le navigateur. Grâce à l’utilisation de Plotly, nous pouvons survoler les différents intervalles de l’histogramme pour faire apparaître les informations quantitatives sous forme d’infobulles. Il est également possible de zoomer, de se déplacer et d’exporter le graphique en png.
Histogramme d’un raster avec R
Version de R : 4.2.1
Il est tout à fait possible de tracer l’histogramme d’un raster avec R. Ici, nous le tracerons sur un raster stocké au format SpatRaster en utilisant le package terra. Après avoir chargé le raster dans le code R, nous utilisons simplement la fonction hist du package terra. Nous pouvons y ajouter des options pour augmenter le nombre d’intervalles à prendre en compte et améliorer l’habillage. Le tout est résumé dans les lignes suivantes.
# chargement de la librairie terra
library(terra)
# chargement du raster à analyser
mnt_roya <- terra::rast("srtm_roya_L93.tif")
# construction de l'histogramme avec un habillage pertinent
hist(mnt_roya,
breaks = 100,
main='Histogramme du MNT de la Roya',
col='lightblue',
xlab='Altitudes (m)',
ylab='Fréquence')
Ici, nous avons utilisé les options suivantes pour améliorer le rendu de l’histogramme :
breaks : le nombre d’intervalles à prendre en compte
main : définition du titre de la figure
col : définition de la couleur des bâtons
xlab : le titre de l’axe des abscisses
ylab : le titre de l’axe des ordonnées
Nous obtenons au final le résultat suivant (Fig. 236).
Fig. 236 Histogramme d’un raster calculé avec R.
Découper un raster
Il est souvent utile de découper un raster afin d’appliquer des traitements seulement sur une zone en particulier. Par découper, nous entendons extraire une sous partie (rectangulaire) d’un raster « trop grand » par rapport à notre zone d’études. Par exemple, si nous disposons d’un modèle numérique de terrain (MNT) à l’échelle de la France et que nous souhaitons simplement étudier les Alpes, il sera pertinent d’extraire simplement la région des Alpes de notre MNT général. Faire ce type d’extraction nous permet de gagner du temps de calcul, d’alléger nos données d’entrée et de cibler spécifiquement notre méthode sur notre zone d’intérêt.
Il s’agît d’une opération extrêmement courante qu’il est possible de faire avec de nombreux outils différents.
Dans les exemples suivants, nous allons extraire d’une (ou plusieurs) bande spectrale de la région parisienne, un zoom sur la ville de Paris.
Découper un raster dans QGIS
Version de QGIS : 3.18.1
Nous commençons par charger dans QGIS le raster que nous souhaitons découper et nous repérons visuellement la zone à extraire (Fig. 237).
Fig. 237 Le raster à découper selon le carré rouge.
L’opération se fait simplement en allant dans le menu La fenêtre suivante s’affiche (Fig. 238).
Fig. 238 Découper un raster dans QGIS.
Dans le champ Couche source, nous indiquons le raster à découper. Dans le champ Étendue de découpage, soit nous écrivons à la main les coordonnées du rectangle à découper (ce qui est plutôt rare…), soit nous récupérons ces coordonnées selon trois façons possibles en cliquant sur l’icône 
. Lors du clic sur cette icône, il nous est proposé soit de calculer ces coordonnées selon l’étendue d’une autre couche (Calculer depuis la couche), soit de faire coïncider ces coordonnées avec la vue actuelle dans QGIS (Use Map Canvas Extent) soit de récupérer ces coordonnées en traçant à la main un rectangle dans la fenêtre principale de QGIS (Dessiner sur le canevas). C’est très souvent cette dernière option qu’on utilise. Lors de la sélection de cette option, la fenêtre se cache et le curseur se transforme en croix. Il suffit alors de tracer le rectangle correspondant à la zone à extraire. Une fois ce triangle tracé, les coordonnées s’affichent dans le champ correspondant.
Enfin, dans le champ Découpé, nous spécifions le chemin et le nom du raster découpé qui sera créé. Nous pouvons voir dans le panneau Console GDAL/OGR, qu’en fait ce module ne fait qu’interfacer une commande GDAL. Après avoir cliqué sur Exécuter, le raster découpé apparaît dans la fenêtre principale de QGIS.
Astuce
Il est tout à fait possible de découper de la même manière un raster multi-bandes, qu’il soit sous forme de raster virtuel (.vrt) ou de raster « en dur » (.tif). Dans ce cas, le raster découpé est également multi-bandes. Il est même possible de convertir un raster virtuel en tif multi-bandes et vice versa.
Découper des rasters par lot dans QGIS
Version de QGIS : 3.18.1
Il est tout à fait possible de découper plusieurs rasters se superposant en une seule fois dans QGIS. C’est un cas qui arrive fréquemment lorsque nous travaillons avec des images de télédétection multi-spectrales par exemple. Comme dit précédemment, si nous travaillons avec un raster virtuel ou un raster multi-bandes, le découpage par lot est implicite. Ici, nous verrons comment découper plusieurs rasters indépendants mais qui se superposent. Dans cet exemple, nous allons découper plusieurs bandes spectrales Landsat centrées sur la région parisienne afin de n’avoir que des zooms sur l’étendue de la ville de Paris.
Une fois les rasters chargés dans QGIS, comme pour découper un seul raster, nous allons dans le menu La fenêtre dédiée de paramétrage du découpage s’ouvre. À ce niveau, nous cliquons sur Exécuter comme processus de lot… dans le bas de la fenêtre. La fenêtre suivante apparaît (Fig. 239).
Fig. 239 Découper des rasters par lot dans QGIS, initialisation.
Chaque ligne correspondra à un raster associé aux paramètres de découpage. Il faut donc ajouter autant de lignes que de rasters à découper. Soit on ajoute les lignes une par une en cliquant sur l’icône 
soit on sélectionne les rasters à découper parmi ceux déjà chargés dans QGIS. Nous cliquons alors sur Auto-remplissage… dans la colonne Couche source, puis nous sélectionnons Sélectionner à partir des couches chargées. Il est également possible d’ajouter des rasters depuis un répertoire sans les avoir chargés au préalable dans QGIS. Une fois cette sélection faite, autant de lignes que de rasters sont apparues, et chaque ligne correspond à un raster d’entrée.
Ensuite il est nécessaire de définir l”Étendue de découpage. Comme dans le cas d’un seul raster, nous pouvons définir l’étendue à la main en traçant un rectangle via l’option Dessiner sur le canevas accessible via l’icône . Afin d’avoir rigoureusement la même étendue de découpage pour chacun des rasters nous faisons un copier-coller des étendues de découpage de ligne en ligne.
Les deux champs suivants peuvent être laissées par défaut et nous renseignons finalement le chemin et le nom des rasters découpés qui seront créés dans la colonne Découpé (étendue). Il est possible de faire des copier-coller des chemins en changeant simplement le nom du raster créé (Fig. 240).
Fig. 240 Découper des rasters par lot dans QGIS.
Si nous avons coché la case Charger les couches dans la fenêtre précédente, les rasters découpés apparaissent automatiquement, sinon il faut les charger manuellement.
Découper un raster dans R
Version de R : 3.18.1
Version de terra : 1.7.29
Il est possible de découper un raster selon une l’emprise d’une couche veceur (ou raster). Pour ce faire nous pouvons passer par la librairie terra et charger un raster au format SpatRaster et un vecteur au format SpatVector. Ici nous découperons une bande spectrale Landsat selon l’emprise du département des Bouches-du-Rhône.
# chargement du raster et du vectuer
b4 <- terra::rast('Landsat_13/LC08_L2SP_196030_20190613_20200828_02_T1_SR_B4.TIF')
dep_13 <- terra::vect('dep_13_buffer_32631.gpkg')
# découpage du raster selon l'emprise du vecteur
b4_emprise_dep13 <- terra::crop(b4, dep13)
Si maintenant le but est de vraiment masquer le raster selon les limites précises du vecteur, dans ce cas là il s’agît de découper un raster selon une couche vecteur, comme plus loin (Découper un raster avec R).
Créer une mosaïque de rasters
Avant de se lancer dans des traitements sur des rasters il peut être utile de créer des mosaïques de rasters. Par exemple, si vous travaillez sur un bassin-versant qui se trouve à cheval sur deux dalles de modèles numériques de terrain, il peut être pertinent de mosaïquer les deux dalles avant d’effectuer des traitements. C’est notamment indispensable si vous souhaitez extraire un réseau hydrographique à partir de ce MNT.
Un autre cas peut être l’étude par télédétection d’une région qui se trouve à cheval sur deux images satellites. Dans ce cas une mosaïque des images peut être intéressante avant les traitements. Il faudra tout de même faire attention à ce que les images à mosaïquer soit similaires (saison, couverture nuageuse, niveau de traitement…).
Pour illustrer nous allons créer une mosaïque de deux rasters multi-bandes contenant 6 bandes Landsat 8. Ces deux images Landsat 8 recouvrent respectivement le département du Val-de-Marne et la ville de Paris.
Notons qu’il est plus prudent que les rasters à mosaïquer soient tous dans le même système de coordonnées.
Mosaïque de rasters dans QGIS
Version de QGIS : 3.18.3
Nous commençons par charger dans QGIS les deux rasters correspondant à nos deux images Landsat 8 (Fig. 241).
Fig. 241 Les deux rasters à mosaïquer.
Nous voyons tout de suite que les cadres noirs vont gêner le processus de mosaïque. Ces cadres noirs sont bien sûr les pixels en no data. La première étape consiste à relever la valeur numérique de ces pixels. Cela nous permettra de les ignorer par la suite. Nous relevons les valeurs de ces pixels via l’outil Identifier des entités 
de la barre d’outils principale de QGIS. Dans notre cas nous constatons que ces no data ont une valeur de 0 sur toutes les bandes (ce qui n’est pas des plus judicieux mais bon…).
Pour créer la mosaïque, nous allons dans le menu Le menu suivant apparaît (Fig. 242).
Fig. 242 Paramétrage de la mosaïque.
Dans le champ Couches en entrée nous sélectionnons les rasters à mosaïquer via l’icône . Il est possible, mais pas forcément utile, de cocher la case Placer chaque fichier en entrée dans une bande séparée. Le point intéressant ici est de spécifier la valeur de no data à utiliser. Ce réglage se fait dans le champ Valeur de pixel à considérer comme NoData dans le bloc Paramètres avancés. Ici nous la paramétrons à 0. Enfin, dans le champ Fusionné nous indiquons le chemin et le nom du raster de mosaïque qui sera créé. Nous pouvons retrouver la commande GDAL sous-jacente à cette manipulation dans le bloc Console GDAL/OGR. Nous cliquons ensuite sur Exécuter. La mosaïque apparaît automatiquement (Fig. 243).
Fig. 243 Résultat de la mosaïque de deux images Landsat 8.
Astuce
Sur la figure Fig. 243, dans la symbologie nous avons réglé les no data à 0 qui apparaissent alors en transparent, donc en blanc.
Le point intéressant est qu’avec cet outil les rasters multi-bandes sont mosaïqués bande à bande. La mosaïque résultante est donc également multi-bande.
Mosaïque de rasters dans R
Version de R : 4.3.1
Version de terra : 1.7.9
Lorsque deux rasters adjacents ou se superposant sont chargés en tant qu’objets SpatRaster, il est possible de les mosaïquer grâce à la fonction merge de la librairie terra. Le mot clef first=TRUE signifie que là où les deux rasters se superposent, les valeurs du premier raster sont prises pour la mosaïque. Il suffit de mettre ce mot clef à False pour ce soient les valeurs du second raster.
# chargement des deux rasters
rast1 <- terra::rast('Landsat_Paris1.tif')
rast2 <- terra::rast('Landsat_Paris2.tif')
# mosaïque des deux rasters
mos <- terra::merge(rast1, rast2, first=TRUE)
# affichage de la mosaïque
terra::plot(mos)
Astuce
La mosaïque avec terra fonctionne aussi sur des rasters multi-bandes.
Raster multi-bandes
La création d’un raster multi-bandes est souvent nécessaire pour différents traitements géomatiques, notamment en télédétection. Un raster multi-bandes est un fichier raster qui contient plusieurs sous rasters superposés les uns aux autres. C’est typiquement utilisé pour la manipulation d’images de télédétection multi-spectrales comme les images Landsat ou Sentinel par exemple. Selon le logiciel utilisé, la création d’un tel raster est nommé concaténation ou fusion. En franglais il est également possible de parler de stack (pile).
Une telle concaténation est souvent nécessaire pour effectuer une classification d’image satellite, une composition colorée ou extraire des signatures spectrales. Cette manipulation est faisable avec la plupart des outils de géomatique.
Création d’un raster multi-bandes avec QGIS
Version de QGIS : 3.26.1
Dans QGIS, il existe trois possibilités de concaténation de rasters. Les deux premières consistent à créer un vrai raster multi-bandes en dur et la troisième consiste à créer un raster virtuel pointant vers les rasters initiaux. La différence entre les deux premières est que dans un cas nous laissons QGIS simplement numéroter les bandes et dans le second cas nous utilisons un module supplémentaire qui permet de nommer les bandes à l’intérieur du raster fusionner.
Raster multi-bandes en dur (natif)
Avertissement
Si vous souhaitez pouvoir nommer explicitement les bandes de votre raster fusionné, allez directement à la section Raster multi-bandes en dur (avec noms).
Dans cet exemple, nous allons concaténer six bandes spectrales Landsat en un seul raster multi-bandes. Nous commençons par charger dans QGIS nos rasters de bandes spectrales. Nous avons bien six rasters correspondant à nos six bandes spectrales. Pour concaténer ces rasters, nous allons dans le menu . Le menu de fusion s’ouvre (Fig. 244).
Fig. 244 Concaténation de rasters dans QGIS.
Dans le champ Couches en entrée nous spécifions les rasters à concaténer. Il est possible de les sélectionner en cliquant sur l’icône (Fig. 245).
Fig. 245 Sélection des rasters à fusionner.
Ici, il faut bien faire attention à l’ordre des rasters car ce même ordre sera repris dans le raster concaténé. Nous sélectionnons les rasters à fusionner, ici nous pouvons cliquer sur Sélectionner tout. La fenêtre précédente se met à jour et nous voyons bien que dans le champ Couches en entrée nous avons 6 inputs selected.
Il est ensuite nécessaire de bien sélectionner Placer chaque fichier en entrée dans une bande séparée, sinon le résultat ne sera pas celui attendu.
Dans le champ Fusionné nous spécifions un chemin et un nom pour le raster concaténé résultant. Nous pouvons noter que ce menu n’est finalement qu’un interface à une commande GDAL. Nous avons accès à cette commande dans le panneau Commande GDAL/OGR. En cliquant sur Exécuter, le processus se lance et le raster concaténé apparaît automatiquement dans QGIS. Nous pouvons remarquer qu’une composition colorée par défaut est automatiquement générée.
Avertissement
Le numéro des bandes est un point important à bien saisir et qui prête à confusion. Le numéro des bandes du raster concaténé est totalement indépendant du « numéro » des bandes en entrée. Par exemple, si nous souhaitons concaténer les bandes 2 (Bleu), 3 (Vert) et 4 (Rouge) d’une image Landsat 8, nous obtiendrons un raster concaténé avec des bandes numérotées de 1 à 3. Mais il faudra bien faire attention au fait que la bande 1 du raster concaténé ne correspondra pas à la bande 1 de Landsat, mais à la bande 2 et ainsi de suite.
Raster multi-bandes en dur (avec noms)
La solution native de QGIS (Raster multi-bandes en dur (natif)) fonctionne parfaitement mais les bandes du raster fusionné ne sont pas nommées, ce qui peut s’avérer fastidueux pour, par exemple, construire des compositions colorées avec QGIS (Composition colorée avec QGIS). Il faut se souvenir à quoi correspond la bande 1, la bande 2, etc du raster multi-bandes. Ça peut être source de confusion. Pour palier cette limitation, il est possible d’utiliser le module supplémentaire nommé Rename Bands. Il s’installe comme toute extension de QGIS (Installation d’extensions dans QGIS). La seule chose est que cette extension est considérée comme étant expérimentale (malgré qu’elle fonctionne très bien). Pour y avoir accès, il est donc nécessaire d’aller dans le menu . Dans cet onglet, il faut cocher Afficher les extensions expérimentales. Une fois cette option cochée, retourner dans Toutes et l’extension se trouve maintenant bien dans la liste.
Concrètement, ce module va créer un raster fusionné en dur classique mais il va l’accompagner d’un fichier texte au format .xml qui portera le même nom que le raster fusionné et dans lequel les bandes seront nommées. Il faudra bien faire attention à toujours conserver ces deux fichiers ensemble si nous souhaitons conserver les noms des bandes.
Une fois le module installé, nous commençons par charger les rasters à fusionner dans QGIS. Ensuite, le module est accessible dans la . Une fois le module lancé, la fenêtre suivante s’affiche (Fig. 246).
Fig. 246 Création d’un raster virtuel dans QGIS en nommant les bandes.
Dans le panneau Raster layers la liste des rasters chargés dans QGIS s’affiche. Nous devons sélectionner les rasters à fusionner dans ce panneau. Une fois sélectionnés ils apparaissent surlignés en bleu. Les rasters sélectionnés apparaissent dans le panneau suivant, à la colonne Band name. À ce niveau il est possible de changer l’ordre des bandes à l’aide des icônes 
et 
et d’ajouter ou enlever un raster avec les icônes et 
.
Dans ce même panneau, à la colonne Band future name nous indiquons le nom explicite à donner à chaque bande. Il suffit de cliquer sur la ligne et une fenêtre d’entrée de texte Band name surgit. Nous pouvons y inscrire ce que nous souhaitons. Ici, nous travaillons sur une image Landsat 8 (OLI/TIRS - Landsat 8 / OLI-2/TIRS-2 - Landsat 9), nous choisissons donc les noms suivants B1 côtier, B2 Bleu, B3 Vert, etc
Si nous le souhaitons, à la ligne No data value nous pouvons indiquer une valeur qui correspondra aux pixels sans données. Enfin, à la ligne Raster output, nous indiquons un chemin et le nom du raster multi-bandes qui sera créé. Nous terminons en cliquant sur Exécuter. Le raster multi-bandes apparaît dans QGIS et nous constatons que les bandes sont bien nommées comme souhaité (Fig. 247).
Fig. 247 Raster fusionné avec noms de bandes explicites.
Notons que dans le répertoire d’export, nous trouvons bien notre raster fusionné au format .tif ainsi que le fichier texte .aux.xml qui contient le nom des bandes. Par curiosité il est possible d’ouvrir ce fichier XML avec un éditeur de texte. Nous y retrouvons bien le nom de chaque bande à chaque balise Description (Fig. 248).
Fig. 248 Le fichier XML accompagnant le raster fusionné.
Raster virtuel
Avec QGIS il est également possible de créer un raster virtuel. Un tel raster est qualifié de virtuel car le fichier généré ne sera pas un fichier raster mais un simple fichier texte pointant vers les rasters initiaux. Il s’agît d’une sorte de raccourci. L’avantage d’un raster virtuel est de ne pas dupliquer les données comme c’est le cas avec une concaténation en dur.
Avertissement
Un raster virtuel ne sera souvent pas accepté comme format d’entrée pour une classification d’images. Il faudra alors préférer un vrai raster concaténé.
Dans cet exemple, nous allons de nouveau concaténer nos six bandes spectrales Landsat. Une fois ces bandes chargées, nous allons dans le menu . La fenêtre Construire un raster virtuel s’affiche alors (Fig. 249).
Fig. 249 Création d’un raster virtuel dans QGIS.
La procédure est très similaire à celle de la concaténation de rasters en dur. Dans le champ Inuput layers nous spécifions les rasters à concaténer de la même manière que précédemment. L’ordre des couches en entrée est important car l’ordre sera le même dans raster virtuel final.
Le champ Resolution est important. Le raster virtuel homogénéisera obligatoirement les résolutions spatiales des différents rasters en entrée. Il est possible de prendre la résolution moyenne (Average), la plus haute (Highest) ou la plus basse (Lowest). Lorsque toutes les images ont la même résolution, comme dans cet exemple, l’option choisie n’a pas d’importance. Par contre si un des rasters a une résolution plus fine, le choix sera important. De plus, lors de l’homogénéisation des résolutions, une interpolation spatiale sera effectuée de façon implicite. Il sera possible de choisir la nature de cette interpolation dans les Paramètres avancés.
Il ne faut pas oublier de cocher la case Place each input file into a separate band. Enfin, dans le champ Virtuel nous spécifions un chemin et un nom pour le raster virtuel résultat. Il s’agît encore d’une commande GDAL interfacée, comme nous pouvons le voir dans le panneau Console GDAL/OGR. À la fin de l’exécution du menu, le raster virtuel apparaît automatiquement dans QGIS associée à une composition colorée par défaut. Comme pour le raster concaténé, il faut être prudent sur les numéros de bandes, voir l’avertissement à ce sujet dans la section précédente.
Si nous regardons dans notre explorateur de fichiers, nous constatons bien que le fichier créé est en fait un simple fichier texte avec l’extension .vrt. Il est possible de l’ouvrir dans un éditeur de texte de type bloc-note (Fig. 250).
Fig. 250 Raster virtuel dans un éditeur de texte.
Ce fichier texte est un fichier XML. Nous y retrouvons les chemins vers les rasters initiaux. Par exemple, ici nous retrouvons le chemin du premier raster au niveau de la première balise SourceFilename. Nous retrouvons également d’autres informations relatives aux dimensions du raster, à son encodage, à la méthode de rééchantillonnage…
Avertissement
Il faut être prudent lorsque nous déplaçons ou faisons un copier-coller d’un raster virtuel via l’explorateur de fichiers. Comme nous l’avons vu, il s’agît d’un raccourci pointant vers les rasters initiaux. Lorsque nous déplaçons un fichier .vrt il faut également déplacer les rasters initiaux. Le plus simple est de créer le raster virtuel dans le même répertoire que les rasters initiaux et déplacer tous les fichiers en même temps lorsque c’est nécessaire.
Création d’un raster multi-bandes avec R
Version de R : 4.2.1
Il est très simple de créer un raster multi-bandes avec R. Ici nous présenterons le cas de la création d’un raster contenant quatre bandes spectrales Landsat avec le package terra. Après avoir chargé les rasters mono bandes, nous les combinons grâce à la fonction c() qui signifie concatenate. Le tout se fait avec les lignes suivantes.
# chargement du package terra
library(terra)
# chargement dess rasters mono-bandes à fusionner
b2 <- terra::rast("Landsat_13/LC08_L2SP_196030_20190613_20200828_02_T1_SR_B2.TIF")
b3 <- terra::rast("Landsat_13/LC08_L2SP_196030_20190613_20200828_02_T1_SR_B3.TIF")
b4 <- terra::rast("Landsat_13/LC08_L2SP_196030_20190613_20200828_02_T1_SR_B4.TIF")
b5 <- terra::rast("Landsat_13/LC08_L2SP_196030_20190613_20200828_02_T1_SR_B5.TIF")
# création du raster multi-bandes
multi_bandes <- c(b2, b3, b4, b5)
# affichage du résultat
plot(multi_bandes, col=gray.colors(100, start = 0.1, end = 0.9))
Nous obtenons bien un objet SpatRaster à 4 bandes (Fig. 251).
Fig. 251 Création d’un raster multi-bandes avec R.
Une fois ce raster multi-bandes créé, il est possible de changer le nom des bandes, grâce à la fonction names(), comme indiqué ci-après.
# renommage des bandes selon leur numéro
names(multi_bandes[[1]]) <- 'Bleu'
names(multi_bandes[[2]]) <- 'Vert'
names(multi_bandes[[3]]) <- 'Rouge'
names(multi_bandes[[4]]) <- 'Proche infrarouge'
# affichage du résultat
plot(multi_bandes, col=gray.colors(100, start = 0.1, end = 0.9))
Les bandes du raster multi-bandes ont bien été renommées comme souhaité (Fig. 252).
Fig. 252 Renommer les bandes d’un raster multi-bandes dans R.
Défusionner un raster multi-bandes avec SCP (QGIS)
Il est tout à fait possible de « défusionner » un raster multi-bandes afin de récupérer chaque bande dans un fichier .tif indépendant. Une solution consiste à passer par le module SCP.
Une fois le raster multi-bandes au format .tif chargé dans QGIS, il suffit d’aller dans le menu . Le menu suivant s’affiche (Fig. 253).
Fig. 253 Séparer un raster multi-bandes avec SCP.
Dans le champ Sélectionner un raster multibandes, nous sélectionnons le raster à séparer. Si nécessaire nous actualisons la liste des rasters en cliquant sur l’icône 
. Nous pouvons spécifier un préfixe pour les rasters indépendants qui seront créés, mais ce champ peut rester vide. Enfin nous cliquons sur Lancer. Il nous est alors demander un répertoire pour sauver les rasters séparés. Les rasters indépendants s’affichent automatiquement.
Rapport de valeurs uniques
Dans de nombreux cas il est utile de compter le nombre d’occurrences de chaque valeur de pixels sur un raster. C’est-à-dire de lister toutes les valeurs que prennent les pixels d’un raster donné et, pour chaque valeur rencontrée, compter le nombre de fois où cette valeur apparaît.
Cette manipulation est souvent appelée Rapport de valeurs uniques.
Évidemment, ce rapport de valeurs uniques se calcule sur des rasters discrets, qui ne prennent qu’une certaine gammes de valeurs. C’est typiquement employé sur les rasters d’usage du sol. Pour chaque classe d’usage du sol, nous pouvons obtenir le nombre de pixels. Et connaissant la résolution spatiale du raster (donc la superficie d’un pixel), il est possible de calculer la superficie occupée par chaque classe d’usage du sol.
Cette manipulation peut s’effectuer de différentes manières, avec différents outils. Dans les exemples présentés ci-dessous, nous calculerons le rapport de valeurs uniques sur un raster d’occupation du sol de la région du Caire. Ce raster a été obtenu par classification supervisée et présente 7 classes d’occupations du sol comme présentées dans le tableau suivant.
ID classe |
Label classe |
|---|---|
1 |
Eau |
2 |
Végétation inondée |
3 |
Végétation dense |
4 |
Sol nu sableux |
5 |
Sol nu rocheux |
6 |
Bâti dense |
7 |
Bâti diffus |
Rapport de valeurs uniques avec QGIS
Version de QGIS : 3.18.3
Nous commençons par charger dans QGIS le raster sur lequel nous souhaitons calculer le rapport de valeurs uniques. Ensuite, nous allons dans la . Le menu suivant apparaît (Fig. 254).
Fig. 254 Calcul du rapport de valeurs uniques dans QGIS.
Dans le champ Couche source nous indiquons le raster sur lequel nous travaillons. Dans le champ Numéro de bande nous sélectionnons la bande du raster à utiliser (seulement utile dans le cas d’un raster multi-bandes). Dans le champ Rapport de valeurs uniques nous indiquons un chemin et un nom pour le rapport qui sera créé. Ce rapport sera au format HTML et sera donc consultable avec un navigateur Web. Puis dans le champ Table de valeurs uniques nous indiquons un chemin et un nom pour la table contenant le rapport qui sera créé. Ces deux sorties sont optionnelles. Puis nous cliquons sur Exécuter.
La table de rapport créée apparaît automatiquement et est nommée par défaut report. Il s’agît d’un fichier GeoPackage sans attribut géographique. Nous pouvons ouvrir ses attributs en faisant un clic droit dessus dans le panneau des couches et en sélectionnant Ouvrir la Table d'Attributs (Fig. 255).
Fig. 255 Table de rapport des valeurs uniques.
La première colonne fid est juste une colonne technique. La colonne value contient toutes les valeurs uniques trouvées dans le raster, soit, ici, les identifiants des 7 classes d’occupations du sol. La valeur 0 correspond aux no data. La colonne count contient le nombre de pixels trouvés pour chaque valeur. La colonne m2 est la conversion du nombre de pixels en superficie exprimée en unités du raster. Cette superficie est calculée selon la résolution spatiale du raster.
Si nous l’avons exporté, nous trouvons également un fichier HTML de rapport nommé report.html contenant ces mêmes informations plus quelques autres comme le SCR du raster, son emprise…
Note
Le fait d’exporter le rapport en une table GeoPackage peut être utile si nous souhaitons ensuite faire une jointure attributaire avec ce rapport.
Rapport de valeurs uniques avec SCP
Version de QGIS : 3.18.3
Version de SCP : 7.8.21
Le module supplémentaire SCP de QGIS offre aussi une fonction pour calculer un rapport de valeurs uniques. Nous commençons par charger le raster à analyser dans QGIS. Ensuite nous allons dans le menu . Bien que la fonction soit ici nommée « Rapport de classification », elle est applicable à tout raster discret même non issu de classification. Le menu suivant s’affiche (Fig. 256).
Fig. 256 Rapport de valeurs uniques avec SCP.
Dans le champ Sélectionner la classification, nous sélectionnons la classification à analyser depuis la liste déroulante. Il peut être nécessaire de rafraîchir la liste en cliquant sur l’icône. Il est possible de définir une valeur de no data à ne pas compter dans le rapport. Ici ce sera la valeur 0. Il n’y a plus qu’à cliquer sur Lancer. Il est demandé un chemin et un nom pour le fichier de rapport au format csv qui sera créé. Puis le rapport apparaît sous forme de texte dans le panneau Sortie du menu précédent (Fig. 257).
Fig. 257 Sortie du rapport de valeurs uniques avec SCP.
La colonne class répertorie toutes les valeurs de pixels identifiées sur le raster. La colonne PixelSum présente le nombre d’occurrences de chaque valeur de pixel. La colonne Percentage % indique le pourcentage que représente chaque valeur de pixel par rapport au raster total. Enfin, la colonne Area traduit le nombre de pixels en superficie en mètres carrés. Ce rapport texte se retrouve exactement dans le fichier csv exporté.
Note
Le petit plus de cette façon de procéder est de sortir un fichier csv facilement chargeable dans un tableur de type Excel ou LibreOffice Calc ou dans un script R ou Python. Le fait d’avoir directement les pourcentages est également intéressant.
Rapport de valeurs uniques avec R
Version de R : 4.8.1
Il est possible de sortir un rapport de valeurs uniques dans R à partir d’un raster. Dans cet exemple, nous allons générer un rapport de valeurs uniques pour un raster d’occupation du sol d’une région de la ville du Caire. Ce raster d’occupation du sol sera chargé sous forme d’un objet SpatRaster via terra. La fonction à utiliser pour générer ce rapport est freq(). La procédure à suivre est décrite ci-après.
# chargement de la librairie terra
library(terra)
# chargement du raster d'occupation du sol
usol <- terra::rast("xtr_usol_2011.tif")
# rapport de valeurs uniques
val_uniques <- freq(usol)
# rapport de valeurs uniques pour une valeur donnée
val_uniques <- freq(usol, value=5)
Le rapport est ici stocké sous forme d’un dataframe dans la variable val_uniques. Comme présenté dans le code ci-dessus, il est possible de ne générer ce rapport que pour une valeur donnée grâce à l’ajout de value=ma valeur.
Calcul raster mono-bande
Fréquemment il est utile d’effectuer des calculs sur une ou plusieurs couches rasters. Un tel calcul peut par exemple être utile pour isoler les altitudes supérieures à un certain seuil sur un MNT ou sur tout autre raster ou pour calculer un indice radiométrique à partir de rasters de bandes spectrales comme le NDVI.
En général, il y a trois formes de calculs rasters possibles. La première consiste à combiner un raster et un scalaire par une opération arithmétique. La deuxième consiste à appliquer un opérateur logique sur un raster afin d’identifier les pixels respectant la condition. Enfin, la troisième consiste à combiner plusieurs rasters par un ou plusieurs opérateurs arithmétiques. Dans ce dernier cas, les rasters en entrée doivent, dans la plupart des cas, présenter la même résolution spatiale, la même étendue et être géoréférencés dans le même système de coordonnées. Bien sûr, ces trois formes peuvent faire intervenir plusieurs rasters ou plusieurs scalaires en même temps, ainsi que plusieurs opérateurs arithmétiques et logiques.
Dans tous les cas, la logique est que le calcul s’effectue pixel par pixel comme schématisé sur la figure suivante (Fig. 258). Sur cette figure, nous considérons le cas d’un raster composé de 9 pixels, en 3 lignes et 3 colonnes.
Fig. 258 Les trois principaux cas de figures du calcul raster.
Dans l’exemple précédent, pour le premier cas de figure, chaque pixel est multiplié par 10. Le raster résultat est géographiquement identique au raster d’entrée mais chaque pixel vaut maintenant 10 fois plus que le pixel d’entrée. Dans le deuxième cas, nous avons extrait les pixels dont la valeur est strictement supérieure à 3. Pour chaque pixel, nous avons comparé sa valeur au chiffre 3. Si le pixel était supérieur, alors on le recodait en 1 et dans le cas contraire en 0. En effet, en informatique il n’est pas possible de stocker directement les booléens Vrai et Faux sur les pixels. Les pixels sont obligatoirement de type numérique. Par défaut, Vrai se code en 1 et Faux en 0. Mais notons bien que ces valeurs de 1 et 0 n’ont aucune signification quantitative. Enfin, par défaut (souvent mais pas toujours) à l’affichage, les pixels à 1 apparaissent en blanc et les pixels à 0 en noir. Enfin, dans le troisième cas, nous avons additionné deux rasters géographiquement identiques. Le raster produit contient pour chaque pixel la somme des pixels des deux rasters d’entrée.
Calcul raster avec QGIS
Version de QGIS : 3.20.1
Le calcul raster peut se faire très simplement avec QGIS en utilisant la calculatrice raster intégrée. Elle permet de faire des calculs numériques et logiques.
Calcul logique
Dans cet exemple, nous allons extraire les zones du bassin de la Roya dont l’altitude est supérieure à 2000 m. Nous commençons par charger le MNT du bassin dans QGIS. Puis l’extraction va se faire grâce à l’opérateur logique » > « . Pour ouvrir la calculatrice raster nous allons dans le menu . La calculatrice s’ouvre (Fig. 259).
Dans le panneau Bandes raster nous trouvons la liste des rasters chargés dans QGIS. Nous remarquons qu’après le nom du raster nous trouvons un @ suivi d’un numéro (srtm_roya@1). Le numéro correspond au numéro de bande dans le raster. Dans le cas d’un raster multi-bandes, nous aurions mon_raster@1, mon_raster@2 et ainsi de suite. Dans le panneau Opérateur nous trouvons les opérateurs disponibles dans la calculatrice et dans le panneau Couche résultat les réglages pour le raster qui sera créé.
Nous allons maintenant entrer l’expression qui va nous permettre d’extraire les altitudes supérieures à 2000 m dans le panneau Expression de la calculatrice raster. Nous commençons par sélectionner le raster qui rentre dans le calcul en double-cliquant sur son nom dans le panneau Bandes raster. Le raster apparaît alors dans le panneau d’expression. Nous entrons ensuite la suite de l’expression soit au clavier, soit en nous servant des opérateurs du panneau dédié.
« srtm_roya@1 » > 2000
À la ligne Couche en sortie, nous spécifions un chemin et un nom vers le raster qui sera créé. Nous pouvons également régler le SCR du raster résultat à la ligne SCR en sortie, généralement nous laissons le même SCR que le raster en entrée. Les réglages de l’emprise du résultat peuvent garder leurs valeurs par défaut (Fig. 259). Nous cliquons ensuite sur OK.
Fig. 259 Extraction des altitudes supérieures à 2000 m du bassin de la Roya via la calculatrice raster de QGIS.
Avec ce type de calcul logique, pour chaque pixel du raster nous regardons si son altitude est supérieure ou non à 2000. La réponse est Oui ou Non, c’est binaire. Par convention, en logique numérique, Oui est transcrit en 1 et Non est transcrit en 0. Ainsi, sur le raster résultat, les pixels supérieurs à 2000 m seront codés en 1 et ceux inférieurs à 2000 m seront codés en 0. Par défaut, dans la plupart des logiciels, les pixels à 1 sont colorés en blanc et les pixels à 0 le sont en noir. Nous obtenons ainsi le raster binaire suivant (Fig. 260 B).
Fig. 260 MNT initial (A), raster binaire des altitudes supérieures à 2000 m en blanc et inférieures en noir (B), altitudes supérieures à 2000 m en vert et inférieures en transparent (C).
Il est tout à fait possible de complexifier ce genre d’extraction à l’aide d’opérateurs booléens comme le ET (AND) et le OU (OR) Ainsi, si nous souhaitons extraire à la fois les altitudes comprises entre 500 m et 700 m et celles comprises entre 1000 m et 1200 m, il nous suffira d’entrer l’expression suivante dans le panneau d’expression de la calculatrice raster.
(« srtm_roya@1 » > 500 AND « srtm_roya@1 »< 700) OR ( « srtm_roya@1 » > 1000 AND « srtm_roya@1 » < 1200)
Les AND servent à définir les bandes d’altitude et le OR sert à sélectionner les deux bandes. Nous obtenons bien un raster binaire correspondant à ces deux gammes d’altitudes (Fig. 260 B).
Fig. 261 Raster binaire des deux bandes d’altitudes extraites.
Calcul sur plusieurs rasters
Il est également possible d’effectuer des calculs combinant plusieurs rasters. Dans l’exemple suivant nous allons calculer un NDVI sur une image Landsat 8 prise au-dessus de Shangaï en août 2020. Après avoir chargé les deux rasters dans QGIS, nous ouvrons la calculatrice raster et nous entrons simplement la formule du NDVI avec les bonnes bandes (OLI/TIRS - Landsat 8 / OLI-2/TIRS-2 - Landsat 9) dans le panneau d’expression en faisant bien attention aux parenthèses (Fig. 262).
Avertissement
Les différents rasters qui entrent dans le calcul doivent être dans le même SCR.
Fig. 262 Calcul de NDVI avec la calculatrice raster de QGIS.
Nous obtenons bien le raster de NDVI de la zone (Fig. 263).
Fig. 263 NDVI calculé pour la zone via la calculatrice raster de QGIS.
Avertissement
Le numéro qui suit le @ dans le nom des rasters dans la calculatrice raster correspond à la position de la bande dans un raster multi-bandes. Il ne correspond en aucun cas au numéro de la bande spectrale d’une image satellite. Par exemple, si vous travaillez avec un raster multi-bandes contenant les bandes 4 et 5 d’une image Landsat 8, ces bandes seront connues comme étant les bandes 1 et 2 du raster multi-bandes.
Calcul raster avec SCP
Version de QGIS : 3.20.1
Version de SCP : 7.9.5
Le module SCP (Module SCP dans QGIS) possède une calculatrice raster similaire à celle de QGIS mais avec quelques fonctionnalités en plus. Il est par exemple possible d’utiliser des opérateurs statistiques et des clauses conditionnelles.
Opérateurs statistiques
Ici nous allons voir un exemple de calcul avec un opérateur statistique appliqué sur le MNT de la Roya. Pour l’exercice nous souhaitons savoir de combien de mètres s’écarte chaque pixel de l’altitude moyenne de la zone. Nous allons ainsi faire la différence entre la valeur de chaque pixel et la valeur moyenne du raster, ce qui est possible en une seule opération avec SCP. Nous allons dans le menu et la calculatrice de bande de SCP s’ouvre (Fig. 264). Il peut être nécessaire de rafraîchir la liste des bandes dans le panneau Liste de bandes en cliquant sur l’icône d’actualisation . Nous entrerons ensuite l’expression dans le panneau Expression en nous servant des fonctions proposées dans le panneau Fonctions. Pour répondre à notre question nous entrons ainsi la formule suivante.
« srtm_roya » - mean(« srtm_roya »)
Note
Nous pouvons récupérer le nom du raster sans l’entrer à la main en double-cliquant sur son nom dans le panneau Liste de bandes. De même, il suffit de double-cliquer sur la fonction mean dans le panneau Fonctions. Nous évitons ainsi des fautes de frappe.
Fig. 264 Calcul avec opérateur statistique dans SCP.
Il suffit alors de cliquer sur Lancer. Il nous est demandé d’indiquer un chemin et un nom pour le raster résultat, puis le raster calculé apparaît (Fig. 265).
Fig. 265 Différence d’altitude entre chaque pixel et l’altitude moyenne de la région.
Sur le raster résultat (Fig. 265), les pixels allant vers le bleu sont les pixels dont l’altitude est inférieure à l’altitude moyenne de la région, et ceux tirant vers le rouge ceux dont l’altitude est supérieure à la moyenne.
Opérateurs logiques
Dans ce module SCP, nous trouvons également des opérateurs de conditions de type Si telle condition est remplie, alors ceci, sinon cela. Cet opérateur est identifié dans SCP sous le nom de where. La syntaxe est where(condition, opération si Vrai, opération si Faux). Par exemple, si nous souhaitons recoder notre MNT en mettant à 10 les pixels dont l’altitude est supérieure à 1000 m et à 20 les pixels dont l’altitude est inférieure à 1000 m, nous entrons la formule suivante dans le panneau d’expression : where(« srtm_roya » > 1000 , 10, 20) (Fig. 266).
Fig. 266 Calcul avec opérateur conditionnel dans SCP.
Le raster résultat est similaire à un raster binaire sauf que, grâce à cet outil de condition, nous pouvons le coder en d’autres valeurs que 0 et 1. Il est même possible d’emboîter les conditions afin d’obtenir plus de deux valeurs.
Calcul raster avec OTB
Version de OTB : 7.2.0
Le logiciel OrfeoToolbox (Orfeo ToolBox) propose également un module de calcul raster. Ce module est basé sur une librairie écrite en C++ nommée muparser optimisé pour le calcul matriciel. Ainsi, la calculatrice raster de OTB est performante en terme de puissance. L’autre avantage de OTB est de pouvoir être utilisé via l’interface graphique Monteverdi ou via la console, ce qui est intéressant dans une optique d’automatisation. Dans les exemples suivants, nous calculerons un NDVI à partir d’une image Landsat 8.
Avec l’interface Monteverdi
Une fois l’interface Monteverdi lancée, nous affichons le Navigateur d'OTB-Applications en allant dans le menu . Le navigateur des applications s’affiche alors sur la partie droite de l’interface dans un onglet dédié. Le module permettant les calculs rasters se nomme BandMath. Une description détaillée est disponible sur la documentation en ligne. Nous cherchons ce module dans la barre de recherche des applications et nous le lançons (Fig. 267).
Fig. 267 Calcul raster avec l’application BandMath de OTB.
Dans le panneau Name, grâce à l’icône en forme de « + » 
nous ajoutons les deux bandes spectrales Landsat 8 qui vont nous servir pour le calcul du NDVI, les bandes Rouge (4) et proche infrarouge (5). Sur la ligne Ouput image, nous entrons le chemin et le nom du raster qui sera créé. Puis dans la ligne Expression, nous entrons l’expression à calculer en utilisant la syntaxe adéquate. Les rasters sont numérotés selon leur ordre dans la liste du panneau Name. Le premier de la liste sera référencé comme im1, le deuxième comme im2 et ainsi de suite. Ensuite, nous devons renseigner la bande à considérer dans chacun de ces rasters. Pour faire un calcul sur la première bande du premier raster nous indiquons im1b1 et ainsi de suite. Dans notre cas nous avons des rasters mono-bande, la formule est donc :
(im2b1 - im1b1) / (im2b1 + im1b1)
où im2b1 correspond à la première bande du deuxième raster chargé, soit bien le raster de la bande du proche infrarouge L8_Shangai_B5.tif
et im1b1 correspond à la première bande du premier raster chargé, soit bien le raster de la bande du rouge L8_Shangai_B4.tif
Nous cliquons sur Execute pour lancer le calcul.
Avec la commande otbcli
En ligne de commande, l’opération est simple, il suffit de bien respecter la syntaxe et de bien se placer dans le répertoire qui contient nos rasters. Nous commençons par appeler le module BandMath, nous définissons ensuite les rasters en entrée via le paramètre -il, puis nous définissons le nom du raster résultat via le paramètre -out et enfin nous entrons la même expression que précédemment mais placée entre guillemets via le paramètre -exp. Dans notre cas, la ligne de commande est la suivante :
otbcli_BandMath -il L8_Shangai_B4.tif L8_Shangai_B5.tif -out ndvi_otbcli.tif -exp "(im2b1 - im1b1) / (im2b1 + im1b1)"
Le point important est de bien faire attention à l’ordre des rasters. Le premier raster de la commande sera le im1, le deuxième im2 et ainsi de suite. Dans tous les cas, le raster calculé est exporté en .tif et peut être ouvert dans QGIS ou tout autre logiciel de géomatique.
Avec un opérateur de condition
Ce module BandMath propose aussi un opérateur de condition, comme vu dans la partie dédié au calcul raster avec SCP (Opérateurs logiques). La syntaxe est (condition ? valeur si Vraie : valeur si Faux). Ainsi, dans notre exemple de coder les altitudes supérieures à 1000 m en 10 et celles inférieures à 1000 m en 20 sur le MNT de la Roya, l’expression à entrer à la ligne Expression (ou après -exp) est :
(im1b1 > 1000 ? 10 : 20)
où im1b1 correspont au raster de MNT.
Calcul raster avec R
Version de R : 4.8.1
Lorsqu’un ou plusieur rasters sont chargés dans R, il est possible d’effectuer des opérations mathématiques dessus. Le calcul raster avec R est très simple et très utile. La façon de calculer est légèrement différente selon les packages utilisés. Ici, nous montrerons la marche à suivre pour un raster chargé sous forme d’un objet SpatRaster via terra. Nous travaillerons sur des bandes spectrales Landsat 8 prises au-dessus des Bouches-du-Rhône.
Raster et scalaire
Il est possible d’appliquer une opération sur un raster, i.e. sur chacun des pixels de ce raster. Dans notre exemple nous multiplierons un raster, de type SpatRaster stocké dans une variable nommé proche_infrarouge, par 1000. Concrètement, nous allons multiplier chacun des pixels constitutifs du raster proche_infrarouge par 1000. Nous obtiendrons ainsi en sortie un nouveau raster aux même propriétés géométriques que le raster initial mais avec des valeurs 1000 fois plus grandes pour chacun des pixels du raster. Ce raster sera aussi un objet SpatRaster que nous stockerons dans une variable nommée pir_1000.
# chargement de la librairie terra
library(terra)
# chargement du raster
pir <- terra::rast("Landsat_13/dep_13_L8_20190613_B5.tif")
# multiplication de ce raster par 1000
pir_1000 <- pir * 1000
# affichage du résultat
plot(pir_1000, col=gray.colors(100, start = 0.1, end = 0.9))
Nous constatons que l’opération est très simple, et nous obtenons bien notre raster initial où chacun des pixels ont été multipliés par 1000 (Fig. 268).
Fig. 268 Résultat de la bande spectrale initiale multipliée par 1000 (remarquer le facteur 1000 dans la légende).
Entre plusieurs rasters
Il est également possible d’effectuer des opérations entre deux, ou plus de deux, rasters. Il est conseillé que les rasters en entrée aient les mêmes dimensions et soient dans les mêmes systèmes de coordonnées. Ici, nous calculerons un NDVI à partir de deux rasters de bandes spectrales du proche infrarouge (pir) et du rouge (rouge). Ces rasters sont au format SpatRaster. Nous stockerons le résultat dans une variable nommée ndvi.
# chargement de la librairie terra
library(terra)
# chargement des rasters des bandes spectrales Rouge et PIR
rouge <- terra::rast("Landsat_13/dep_13_L8_20190613_B4.tif")
pir <- terra::rast("Landsat_13/dep_13_L8_20190613_B5.tif")
# calcul du NDVI
ndvi <- (pir - rouge) / (pir + rouge)
La formulation est très simple. Le raster résultat sera également un objet SpatRaster.
Opération logique
Dans certains cas, nous sommes amenés à appliquer des opérations logiques sur des rasters. À l’issue d’une telle opération logique nous obtenons un nouveau raster dans lequel les pixels répondant Vrai à la condition sont codés en 1 et les pixels répondant Faux sont codés à 0. Le raster obtenu est dit binaire car exprimé en 0 et 1. Ce type d’opération est courante lors de l’application de seuils. Dans l’exemple ici, nous appliquerons un seuil sur le MNT du bassin de la Roya afin d’extraire les pixels dont l’altitude est spérieure à 1500 m.
# chargement de la librairie terra
library(terra)
# chargement du MNT
mnt_roya <- terra::rast("srtm_roya_L93.tif")
# application du seuil sur les altitudes
mnt_roya_1500 <- mnt_roya > 1500
# affichage du résultat
plot(mnt_roya_1500)
Aucune difficulté dans la formulation, et le raster correspondant est bien un raster binaire dans lequel les pixels à 1 correspondent aux pixels ayant une altitude supérieure à 1500 m (Fig. 269).
Fig. 269 Raster binaire après un calcul raster de type seuil dans R.
Calcul raster avec GDAL
Version de GDAL : 3.0.4
GDAL a une fonctionnalité permettant de faire du calcul raster. Comme GDAL s’exécute en lignes de commandes, son usage peut être intéressant pour faire un calcul rapide ou bien un calcul raster dans un script. La commande dédiée est gdal_calc.py. C’est une extension Python de GDAl installée par défaut avec GDAL. Les détails de la commande et la liste des options possibles se trouvent sur la page du manuel dédiée. Nous verrons ici les fonctionnalités les plus fréquentes.
Raster et scalaire
Il est très simple d’associer un raster à un scalaire par un calcul. Par exemple, si nous souhaitons multiplier les altitudes d’un MNT, nommé srtm_Roya.tif, par 10, la commande à entrer est la suivante.
gdal_calc.py -A srtm_roya.tif --outfile=mnt_10.tif --calc="A*10"
Avec :
gdal_calc.py : appel à la commande GDAL
-A srtm_roya.tif : le raster à prendre en entrée qu’on stocke dans la variable A
–outfile=mnt_10.tif : le nom du raster résultat en sortie
–calc= »A*10 » : l’expression à calculer. Bien noter que cette expression fait appel à la variable A et pas directement au raster et que cette expression est entre guillemets et sans espaces
En sortie, nous obtenons bien un raster au format GeoTiff que nous pouvons charger dans QGIS. Sur ce raster, les pixels correspondent aux pixels du MNT de départ mais multipliés par 10.
Raster et calcul logique
Sur le même principe, il est facile d’appliquer une opération logique sur un raster. Par exemple, si nous souhaitons extraire les pixels dont l’altitude est comprise entre 1000 m et 2000 m à partir d’un MNT nommé srtm_Roya.tif, la commande à entrer est la suivante.
gdal_calc.py -A srtm_roya.tif --outfile=alti_1000-2000.tif --calc="(A>1000)&(A<2000)"
Notez que l’opérateur à employer pour le ET est l’esperluette &. Le résultat est un raster binaire où les pixels à 1 correspondent aux pixels respectant la condition et les pixels à 0 correspondent à tous les autres pixels. Cette expression peut également prendre la forme suivante.
gdal_calc.py -A srtm_roya.tif --outfile=alti_1000-2000.tif --calc="logical_and(A>1000, A<2000)"
Ici nous employons le mot clef logical_and() avec la condition à l’intérieur des parenthèses.
Si le calcul fait appel à une condition OU, il faut alors employer le mot clef logical_or(). Par exemple, si nous disposons d’un raster d’occupation du sol de type Corine Land Cover et que nous souhaitons extraire les occupations du sol de classe 211 et 511, la commande GDAL sera la suivante.
gdal_calc.py -A CLC.tif --outfile=CLC_extract.tif --calc="logical_or(A==211,A==511)"
Notez l’emploi de deux signes = pour l’égalité. Nous obtenons bien un nouveau raster binaire avec les pixels à 1 qui correspondent aux pixels satisfaisant la condition.
Calcul sur plusieurs rasters
Il est également tout à fait possible de combiner plusieurs rasters lors d’un calcul. Les rasters d’entrée doivent par contre être de mêmes dimensions et dans le même système de projection. Par exemple, pour calculer un NDVI à partir d’une bande de proche infrarouge nommée L8_B5.TIF et d’une bande rouge nommée L8_B4.TIF, nous entrons la commande suivante.
gdal_calc.py -A L8_B4.TIF -B L8_B5.TIF --outfile=ndvi.tif --calc="(B-A)/(B+A)"
La différence avec précédemment est que nous spécifions deux rasters en entrée stockés respectivement dans les variables A et B. Ce sont ces deux variables qui sont utilisées dans l’expression de calcul que nous retrouvons dans –calc= »(B-A)/(B+A) ». En sortie, nous obtenons bien notre raster de NDVI.
Avertissement
Le raster de sortie d’un calcul raster avec GDAL est du même type que les rasters en entrée. Ainsi, si les rasters en entrée sont de type Integer (Entiers), le résultat sera aussi de type entier, ce qui peut amener à des surprises inattendues et non souhaitées. Il est censé être possible de forcer le type souhaité en spécifiant par exemple A.astype(numpy.float64) dans le calcul mais ça n’a pas l’air de toujours fonctionner…
Si les rasters à combiner dans le calcul ne sont pas constitués de fichiers séparés mais proviennent de bandes différentes se trouvant au sein d’un raster multi-bandes, la syntaxe doit être modifiée. Par exemple, ici nous disposons d’un raster multi-bandes contenant des bandes spectrales Landsat 8 de 1 à 7. La bande du rouge est donc la bande 4 et celle du proche infrarouge la bande 5. Pour calculer le NDVI à partir de ce raster multi-bandes, la syntaxe est la suivante.
gdal_calc.py -A L8_stack.tif --A_band=5 -B L8_stack.tif --B_band=4 --outfile=ndvi_stack.tif --calc="(A-B)/(A+B)"
Notez que nous stockons le raster multi-bandes une première fois dans la variable A et que nous spécifions que nous y prenons simplement la bande 5 via le mot clef –A_band=5. Puis nous stockons le raster multi-bandes une seconde fois dans la variable B et nous spécifions que nous y prenons cette fois-ci la bande 4 via le mot clef –B_band=4. Ainsi les A et B appelés dans l’expression du calcul correspondent bien aux bandes désirées.
Calcul raster avec statistiques
En télédétection il est fréquent de manipuler des rasters multi-bandes. Dans la plupart des cas, ces rasters multi-bandes sont des regroupements de réflectances dans différentes bandes spectrales pour une date donnée ou bien des regroupements d’une même réflectance ou d’un indice pour plusieurs dates. Dans ce second cas nous parlons de séries temporelles. Pour effectuer des calculs, ou des statistiques pixel à pixel, sur ces rasters multi-bandes, il n’est pas possible d’utiliser les outils basiques de calculs rasters. Il est nécessaire d’employer des outils dédiés. Dans tous les cas, le raster résultat sera géographiquement identique aux rasters en entrée mais les valeurs des pixels auront pris la valeur de la statistique utilisée en fonction des rasters en entrée, comme présenté sur la figure suivante (Fig. 270).
Fig. 270 Calcul de statistiques sur un stack de rasters.
Dans cet exemple théorique nous avons en entrée 4 rasters de 9 pixels chacun. Dans un premier temps, nous calculons le raster issu de la statistique Maximum appliquée au stack de ces 4 rasters. Nous obtenons un nouveau raster dans lequel chaque pixel prend la valeur maximale rencontrée dans les 4 pixels d’entrée. Dans un second temps, nous calculons le raster issu de la statistique Moyenne appliquée au stack des 4 rasters. Nous obtenons un nouveau raster dans lequel chaque pixel prend la valeur moyenne des 4 pixels d’entrée.
Calcul raster et stats avec QGIS et SCP
Version de QGIS : 3.22.0
Version de SCP : 7.10.5
Le plugin SCP de QGIS propose une fonctionnalité intéressante pour calculer des statistiques pixel à pixel sur un raster multi-bandes. L’idée est de créer un nouveau raster sur lequel les valeurs des pixels correspondront à un résumé statistique pour chacun des pixels. Dans l’exemple ci-dessous, nous allons calculer la réflectance moyenne dans le proche-infrarouge de 4 images Sentinel-2 prises sur la même zone mais à 4 dates différentes. Nous commençons par charger dans QGIS ces 4 rasters. Puis nous allons dans le menu pour créer un jeu de bandes qui contiendra nos 4 rasters de réflectance.
Une fois le jeu de bandes (Band set) défini, nous allons dans le menu . Le menu suivant s’affiche (Fig. 271).
Fig. 271 Moyenne pixel à pixel sur un raster multi-bandes avec SCP dans QGIS.
Dans le panneau Expression nous entrons l’expression suivante : mean("bandset#b*"). Par défaut, le bandset est le Jeu de bandes actif. Cette syntaxe signifie que nous travaillons sur le bandset actif et sur toutes les bandes le composant #b*. Et sur toutes ces bandes nous demandons la moyenne mean(). Nous cliquons ensuite sur Exécuter et nous spécifions un chemin et un nom pour le raster résultat. Chaque pixel de ce raster résultat correspondra à la moyenne des réflectances sur les 4 dates.
Calcul raster et stats avec OTB
Version de OTB : 7.2.0
OrfeoToolbox propose une fonctionnalité qui permet de faire du calcul multibandes combinant des statistiques. Dans l’exemple ci-après, nous allons combiner 4 rasters de NDVI calculés sur une même région mais à 4 dates différentes au sein d’une année. Nous allons plus précisément construire un nouveau raster, qui sera géographiquement identique aux rasters de NDVI, mais dont les pixels correspondront au NDVI maximum des 4 dates.
Avec l’interface Monteverdi
Avec Monteverdi, ce module s’ouvre en cherchant BandMathX dans la barre de recherche du panneau Navigateur d'OTB-Applications. La fenêtre suivante s’ouvre (Fig. 272).
Fig. 272 Calcul raster avec statistiques sur plusieurs bandes.
Dans le panneau Input list nous renseignons les rasters sur lesquels nous souhaitons calculer des statistiques. Ces rasters doivent avoir exactement les mêmes dimensions et les mêmes systèmes de projection. Ces rasters sont à ajouter à l’aide de l’icône . Ici nous ajoutons nos 4 NDVI calculés pour des dates différentes de l’année 2016. À la ligne Output image, nous spécifions un chemin et un nom pour le raster de sortie. Enfin, le plus important est l’expression à entrer à la ligne Expressions. Dans le cas du calcul d’un maximum, l’expression est :
max(im1b1,im2b1,im3b1,im4b1)
où max est la statistique à calculer et im1b1 signifie que nous prenons la bande 1 (b1) de l’image 1 (im1) correspondant au raster chargé en premier, et ainsi de suite.
Note
Dans le cas où nous ne travaillons pas sur des rasters séparés mais sur des bandes au sein d’un raster multi-bandes, l’expression devient simplement : max(im1b1,im1b2,im1b3,im1b4). Nous ne travaillons qu’avec l’image 1 im1 mais nous y appelons les différentes bandes de b1 à b4.
Une fois cliqué sur Exécuter, le raster contenant la valeur maximums des NDVI est calculé. Il est possible de le charger dans QGIS.
Avec la commande otbcli
Dans une optique d’automatisation de tâches, il est possible d’appeler ce module BandMathX en ligne de commandes otbcli. Cette ligne s’écrit alors simplement comme ci-dessous.
otbcli_BandMathX -il ndvi_2016-01-25.tif ndvi_2016-03-15.tif ndvi_2016-05-04.tif ndvi_2016-08-22.tif -out max_ndvi2.tif -exp "max(im1b1,im2b1,im3b1,im4b1)"
- Où :
otbcli_BandMathX est l’appel au module BandMathX
-il est le mot clef pour spécifier les rasters en entrée
-out est le mot clef pour spécifier le nom du raster en sortie
-exp est le mot clef pour entrer l’expression au même format que vu dans le paragraphe précédent
Note
Il peut être laborieux d’entrer le nom à la main de tous les rasters en entrée. Le plus simple est de créer au préalable un Raster multi-bandes puis d’appeler les différentes bandes de ce raster dans l’expression du calcul.
Reclassifier un raster
Il est souvent utile de reclassifier un raster, c’est-à-dire d’assigner de nouvelles valeurs aux pixels en fonction de règles. Par exemple, il est possible de changer tous les pixels dont les valeurs sont comprises entre une valeur a et une valeur b par une valeur c. Il est également possible de changer tous les pixels qui ont exactement la valeur a par la valeur c. Dans les exemples qui suivent, nous reclassifierons un raster d’usage du sol du Caire en 7 classes, obtenu par classification supervisée, pour obtenir un nouveau raster présentant de nouveaux identifiants de classe. Les 7 classes initiales sont rappelées dans le tableau suivant.
ID classe |
Label classe |
|---|---|
1 |
Eau |
2 |
Végétation inondée |
3 |
Végétation dense |
4 |
Sol nu sableux |
5 |
Sol nu rocheux |
6 |
Bâti dense |
7 |
Bâti diffus |
La reclassification nous permettra d’assigner de nouveaux numéros de classes aux différentes classes ou de regrouper plusieurs classes en une seule classe. Pour l’exemple, dans un premier temps nous assignerons l’identifiant 11 à la classe 1, 22 à la classe 2 et ainsi de suite jusqu’à 77 à la classe 7. Dans un second temps, nous regrouperons les classes par grandes catégories comme indiqué dans le tableau suivant.
ID classe old |
Label old |
ID classe new |
Label new |
|---|---|---|---|
1 |
Eau |
1 |
Eau |
2 |
Végétation inondée |
2 |
Végétation |
3 |
Végétation dense |
2 |
Végétation |
4 |
Sol nu sableux |
3 |
Sol nu |
5 |
Sol nu rocheux |
3 |
Sol nu |
6 |
Bâti dense |
4 |
Bâti |
7 |
Bâti diffus |
4 |
Bâti |
Reclassifier un raster avec QGIS
Version de QGIS : 3.20.3
Une fois le raster à reclassifier chargé, nous ouvrons le menu Reclassification par table dans le panneau . Le menu s’ouvre (Fig. 273).
Fig. 273 Reclassification de raster avec QGIS.
Reclassification classe par classe
Nous commençons par reclassifier le raster classe par classe, sans regroupement. Nous changeons simplement l’identifiant de chaque classe par le nouvel identifiant (1 par 11, 2 par 22…).
À la ligne Couche raster, nous sélectionnons le raster à reclassifier, ici classif_7. Dans le cas d’un raster multi-bandes, nous précisons la bande à reclassifier à la ligne Numéro de bande. Ensuite, le travail de reclassification s’effectue via une table de reclassification à paramétrer à la ligne Table de reclassification en cliquant sur l’icône en bout de ligne. Le menu de paramétrage de la table de reclassification s’ouvre (Fig. 274).
Nous commençons par ajouter autant de lignes que de classes à reclassifier en cliquant autant de fois que nécessaire sur Ajouter une ligne. Ensuite, nous indiquons la borne inférieure et supérieure des gammes de valeurs à changer dans les colonnes Minimum et Maximum et nous renseignons la nouvelle valeur dans la colonne Valeur. Ici, comme nous changeons les valeurs de toutes les classes, classe par classe, les bornes minimum et maximum sont les mêmes et correspondent aux classes à changer (Fig. 274).
Fig. 274 Table de reclassification de raster classe par classe.
Une fois ce paramétrage effectué, nous cliquons sur OK pour revenir à la fenêtre principale du module. Il faut maintenant régler la règle de reclassification en dépliant le panneau Paramètres avancés. Dans ce panneau, à la ligne Limites de plages, nous spécifions la règle à suivre pour la reclassification. Ici, les valeurs maximum et minimum correspondent directement à nos classes. Nous choisissons donc la règle min <= valeur <= max. Nous renseignons ensuite le chemin d’export du raster qui sera créé à la ligne Raster reclassifié. Puis nous cliquons sur Exécuter. Le raster reclassifié apparaît automatiquement. Les classes ont bien été réattribuées.
Reclassification pour regroupement de classes
Dans ce cas de figure, la table de reclassification doit être paramétrée différemment. Les bornes minimum et maximum doivent être renseignées de telle sorte que ce soient des gammes de valeurs qui vont être réassignées. Ici, pour opérer un regroupement de classes tel que présenté dans l’introduction de cette section (Reclassifier un raster), nous paramétrons la table tel que présenté ci-dessous (Fig. 275).
Fig. 275 Table de reclassification de raster pour regroupement de classes.
Puis nous réglons la règle de reclassification dans le panneau Paramètres avancés. Tel que nous avons paramétré nos bornes de classes, la règle appropriée est alors min < valeur <= max. Nous renseignons le chemin du raster de regroupement de classes et nous cliquons sur Exécuter. Le raster aux classes regroupées apparaît (Fig. 276).
Fig. 276 Occupation du sol initiale en 7 classes (A) et occupation du sol après regroupement en 4 classes (B).
Reclassifier un raster avec R
Version de R : 4.3.1
Version de terra : 1.7.29
Nous allons reproduire les manipulations de reclassification présentées kdans l’introduction de cette partie (Reclassifier un raster) avec R en utilisant la librairie terra.
Reclassification classe par classe
Nous commencerons par reclassifier un raster d’occupation du sol classe par classe. À chaque classe initiale nous allons attribuer une nouveau numéro de classe. Cela se fait avec la méthode classify de terra, comme présenté ci-après.
# import de l'occupation du sol à reclasser
usol <- terra::rast('classif_C_2011.tif')
# reclass classe par classe
# on créé une matrice de reclassification
m <- c(1, 11,
2, 22,
3, 33,
4, 44,
5, 55,
6, 66,
7, 77)
rclmat <- matrix(m, ncol=2, byrow=TRUE)
# on applique cette matrice de reclassification
usol_reclass <- terra::classify(usol, rclmat)
Comme nous pouvons le voir, nous commençons par créer une matrice de reclassification de type ancienne valeur, nouvelle valeur puis nous appliquons cette matrice au raster à reclassifier.
Reclassification pour regroupement de classes
Il est également possible de reclassifier en regroupant des classes. Dans cet exemple, nous allons simplifier notre occupation du sol initialement en 7 classes en une occupation du sol en 4 classes (Fig. 275). Nous utilisons également la fonction classify.
# import de l'occupation du sol à reclasser
usol <- terra::rast('classif_C_2011.tif')
# reclass classe par classe
# on créé la matrice de reclassification
m <- c(1, 1, 1,
1, 3, 2,
3, 5, 3,
5, 7, 4)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
# on applique la matrice de reclassification
usol_reclass <- terra::classify(usol, rclmat)
Cette fois-ci la matrice s’écrit sur 3 colonnes de la façon de telle valeur, à telle valeur, nouvelle valeur. Nous pouvons toutefois noter que par défaut, la fonction classify ne prend pas la valeur minimum de chaque intervalle mais prend la valeur maximale. C’est pourquoi nous lui indiquons que les pixels allant de 1 (non inclus) à 3 (inclus) doivent prendre la nouvelle valeur de 2.
Il est également possible de reclassifier spécifiquement certaines valeurs et de mettre toutes les autres valeurs à telle autre valeur. Par exemple, si nous souhaitons reclassifier notre occupation du sol en deux classes Végéation et Pas végétation en mettant les deux classes de végétation à 1 et tout le reste à 0, nous pouvons utiliser le mot clef others, comme présenté ci-dessous.
# import de l'occupation du sol à reclasser
usol <- terra::rast('classif_C_2011.tif')
# reclass en deux classes
# on créé la matrice de reclassification simplement pour la classe d'intérêt
m <- c(1, 3, 1)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
# on applique la matrice de reclassification en spécifiant "tout le reste à 0"
usol_reclass <- terra::classify(usol, rclmat, others=0)
Note
La fonction classify de terra présente d’autres options intéressantes pour des usages plus précis, à regarder de plus près en cas de besoin.
Compresser un raster
Il arrive d’avoir à faire à des rasters qui deviennent trop volumineux pour un ordinateur donné ou pour un logiciel. Une façon de faire peut être de le découper en tuiles, comme expliqué dans la rubrique dédiée (Découper un raster selon une grille). Mais il est également possible de réduite le poids d’un raster en utilisant des algorithmes de compression. Une présentation détaillée de ces algorithmes peut être trouvée en ligne. Ici, nous nous contenterons d’en utiliser un sans regarder les processus sous-jacents. Néanmoins, le point clé à bien retenir est que cette compression n’altère pas les données. Ainsi, les valeurs des pixels sont bien conservées ce qui permet de garder l’intégrité de l’information du raster.
Dans les exemples ci après nous verrons comment compresser un raster d’une bande spectrale Landsat 8.
Compresser un raster avec QGIS
Version de QGIS : 3.26.0
Tout d’abord nous chargeons le raster à compresser dans QGIS. Une fois le raster chargé, il suffit de faire un clic droit sur la couche dans le panneau des couches. Puis nous sélectionnons Le menu suivant s’affiche (Fig. 277).
Fig. 277 Compression d’un raster avec QGIS.
Dans le champ Format nous laissons GeoTIFF. Dans le champ Nom de fichier, nous précisons un chemin et un nom pour le raster compressé qui sera créé. Nous pouvons laisser le champ SCR à sa valeur pour conserver le même SCR. Les options de compression s’obtiennent en dépliant le menu Options de création. Nous choisissons Compression élevée dans le menu Profil et nous pouvons laisser les valeurs proposées par défaut pour les variables Compress, Predictor et Zlevel qui produisent de bons résultats. Pour plus de détails techniques et des comparaisons de différentes compressions vous pouvez vous référer à cette page. Il ne reste qu’à cliquer sur OK.
Le raster passe ainsi d’un poids de 228.6 MiO à 85.2 MiO, soit une diminution de poids de près de 60 %, ce qui n’est pas négligeable du tout. De plus, comme dit en préambule, les valeurs des pixels sont conservées à l’identique.
Compresser un raster avec GDAL
Version de GDAL : 3.0.4
Il peut être intéressant de passer par GDAL en ligne de commandes pour compresser un raster. Ce type d’usage permet d’aller vite ou de chaîner des traitements facilement. Pour cela, il suffit que GDAL soit bien renseigné dans le PATH de votre système d’exploitation. La ligne de commande à utiliser est la suivante.
gdal_translate -of GTiff -co "COMPRESS=DEFLATE" -co "PREDICTOR=2" -co "ZLEVEL=9" LC08_L2SP_196030_20190613_20200828_02_T1_SR_B5.TIF gdal-deflate_2-9.tiff
Cette ligne de commande se comprend de la façon suivante :
gdal_translate : la fonction GDAL dédié à la conversion de format des rasters
-of GTiff : spécification du format de raster en sortie, ici GeoTIFF
-co « COMPRESS=DEFLATE » : l’option qui stipule qu’on va compresser le raster avec l’option DEFLATE. D’autres options sont possibles.
-co « PREDICTOR=2 » -co « ZLEVEL=9 » : les deux options qui précisent le niveau de compression DEFLATE. Il est possible de changer ces valeurs mais celles-ci donnent de bons résultats.
LC08_L2SP_196030_20190613_20200828_02_T1_SR_B5.TIF : le raster à compresser
gdal-deflate_2-9.tiff : le nom du raster compressé qui sera créé
Le raster passe ainsi d’un poids de 228.6 MiO à 85.2 MiO.
Aligner deux rasters
Pour certains traitements, comme le calcul entre plusieurs rasters, il est nécessaire d’aligner des rasters lorsque ceux-ci ne sont pas de la même résolution et donc pas alignés sur une même grille. C’est par exemple le cas lorsque nous calculons certains indices radiométriques faisant appel à des bandes spectrales de différentes résolutions spatiales. Il faudra alors choisir un raster qui servira de référence et aligner le, ou les autres rasters, sur ce raster de référence. Il existe différents outils pour aligner un raster selon un autre raster.
Avertissement
Certains logiciels, comme QGIS, permettent d’emblée de faire des calculs rasters entre rasters de différentes résolutions. Dans ce cas, le logiciel réalise par lui même, de façon cachée, un alignement entre les différents rasters. L’avantage est que l’utilisateur n’a pas à se soucier de cette étape, mais le désavantage est que l’utilisateur n’a pas la main sur l’alignement. Quel raster utilisé comme référence ? Quelle méthode de rééchantillonnage pour l’alignement ? Le logiciel décide sans rien dire…
Avertissement
Un alignement entre deux rasters entraîne obligatoirement un ré-échantillonnage de l’un des rasters. Or la méthode de ré-échantillonnage a son importance. Si nous souhaitons aligner un raster discret (une occupation du sol par exemple) il faut choisir le plus proche voisin, sinon une interpolation linéaire entre les valeurs de classes sera faite, ce qui n’a pas de sens. Par contre, sur un raster continu (comme un MNT ou une bande spectrale) il est possible de choisir une méthode qui interpole.
Aligner deux rasters avec QGIS
Version de QGIS : 3.22.1
Même si QGIS permet de faire des calculs entre rasters de différentes résolutions, il est possible d’aligner des rasters manuellement en choisissant le raster de référence et la méthode de ré-échantillonnage. Pour l’exemple nous allons aligner un raster correspondant à la bande 3 à 10 m d’une image Sentinel-2 sur un raster correspondant à la bande 11 à 20 m.
Une fois les deux rasters chargés dans QGIS, nous allons dans le menu . Dans le menu qui s’ouvre, nous commençons par choisir les rasters à aligner entre eux en cliquant sur le . La fenêtre suivante s’affiche (Fig. 278).
Fig. 278 Sélection du premier raster à aligner dans QGIS.
À Couche d'entrée nous choisissons au moins les deux rasters à aligner. Nous pouvons commencer par le raster à aligner, à savoir celui à 10 m nommé S2_B03_10m. À la ligne Nom de du fichier de raster de sortie nous indiquons un chemin et un nom pour le raster qui sera aligné, par exemple S2_B03_10m_resample.tif. Puis nous devons choisir une Méthode de ré-échantillonnage parmi celles proposées dans le menu déroulant. Nous pouvons choisir Bilinéaire dans notre cas. Enfin nous cliquons sur OK. Nous recliquons sur le pour ajouter le second raster, ici S2_B11_20m. Nous indiquons ce raster à la ligne Nom de du fichier de raster de sortie et nous choisissons aussi un nom de sortie et une méthode de ré-échantillonnage pour ce raster. La fenêtre précédente s’est mise à jour comme suit (Fig. 279).
Fig. 279 Paramétrage de l’alignement des rasters.
Nos deux rasters apparaissent dans le cadre Couches rasters à aligner. Il faut ensuite choisir une Couche de référence, ici S2_B11_20m. Nous pouvons laisser les autres options vides. Si nos deux rasters ne se superposaient pas totalement, nous pourrions régler l’emprise du résultat dans le cadre Découper selon l'emprise actuelle. Enfin nous cliquons sur OK. Bien que deux fichiers rasters ont été générés, celui qui a servi de référence n’est en fait pas modifié.
Aligner deux rasters avec OTB
Version de OTB : 8.1.1
Le logiciel OrfeoToolbox (Orfeo ToolBox) propose un module pour aligner un raster sur un autre. Ce module se nomme Superimpose. Comme tout module de OTB, il est possible de l’utiliser via l’interface graphique Monteverdi ou via la ligne de commande.
Avec l’interface Monteverdi
Une fois l’interface chargé et le panneau des applications affichés, nous cherchons l’application Superimpose dans la barre de recherches du Navigateur d'OTB Applications. Une fois l’application lancée, la fenêtre suivante apparaît (Fig. 280).
Fig. 280 Aligner deux rasters dans OTB.
À la ligne Reference input nous pointons vers le raster à utiliser comme référence, à savoir ici S2_B11_20m.tiff. À la ligne, The image to reproject nous pointons vers le raster à aligner, ici S2_B03_10m.tiff. Puis, nous renseignons le chemin et le nom du raster résultat à la ligne Output image. Une fois ces rasters renseignés, dans le bloc Interpolation nous choisissons une méthode d’interpolation parmi celles proposées dans le menu déroulant. Il ne reste plus qu’à cliquer sur Execute.
Avec la commande otbcli
Il est possible d’exécuter l’application Superimpose via la ligne de commande otbcli. La synataxe est la suivante :
otbcli_Superimpose -inr S2_B11_20m.tiff -inm S2_B03_10m.tiff -interpolator bco -out S2_B03_resample_otbcli.tiff
- Où :
otbcli_Superimpose est l’appel au module Superimpose
-inr est le mot clef pour spécifier le raster à utiliser comme référence
-inm est le mot clef pour spécifier le raster à aligner
-interpolator est le mot clef pour choisir la méthode d’interpolation (bco pour Bicubic interpolation et nn pour Nearest neighbor interpolation)
-out est le mot clef pour le raster de sortie aligné
Aligner deux rasters avec R
Version de R : 4.2.2
Version de Terra : 1.5.34
À l’aide du package terra de R, il est facile d’aligner un raster sur un autre. Il suffit d’utiliser la fonction resample, comme présenté ci-après, où la bande 3 à 10 m d’une image Sentinel-2 est alignée sur la bande 11 à 20 m.
# chargement de la librairie terra
library(terra)
# chargement du raster qui devra être aligné sur un autre
s2_b3_10m <- terra::rast('S2_B03_10m.tiff')
# chargement du raster sur lequel le premier sera aligné
s2_b11_20m <- terra::rast('S2_B11_20m.tiff')
# alignement de la bande à 10m sur la bande à 20m
b3_20m <- terra::resample(s2_b3_10m, s2_b11_20m)
À la fin de ce code, un nouveau raster de type SpatRaster est créé dans la variable b3_20m. Et ce nouveau raster est identique au premier mais qui a été aligné sur le second.