Auteur : Paul Passy

Licence : cc_by_nc_sa

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.

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. 226).

alternate text

Fig. 226 Le raster à découper selon le carré rouge.

L’opération se fait simplement en allant dans le menu Raster ‣ Extraction ‣ Découper un raster selon une emprise… La fenêtre suivante s’affiche (Fig. 227).

alternate text

Fig. 227 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 icone_browse. 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

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 Raster ‣ Extraction ‣ Découper un raster selon une emprise… 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. 228).

alternate text

Fig. 228 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 icone_plus 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 icone_browse. 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. 229).

alternate text

Fig. 229 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 selon une couche vecteur

En plus de découper un raster selon une emprise, il est parfois utile de découper un raster selon un polygone issu d’une couche vecteur. Cette opération de découpe est aussi connu sous le nom de clip. Par exemple si vous travaillez sur un département précis, il peut être intéressant de découper vos rasters (de bandes spectrales, de MNT et autres) selon les contours du département étudié. Cette manipulation permet d’alléger le poids des données et de ne travailler que sur les données précisément contenues dans la zone d’étude.

Cette manipulation est très courante et peut être effectuée avec la plupart des logiciels de géomatique et des langages comme R ou Python.

Dans cet exemple nous allons découper une (puis plusieurs) bande spectrale Landsat 8 selon les contours du département du Val-de-Marne dans l’agglomération parisienne.

Découper un raster dans QGIS

Version de QGIS : 3.18.3

Nous commençons par charger dans QGIS le raster que nous souhaitons découper ainsi que la couche vecteur du polygone qui nous servira pour découper (Fig. 230). Notons que les deux couches doivent être dans le même SCR.

alternate text

Fig. 230 Le raster à découper selon les contours du Val-de-Marne.

Une fois ces couches chargées, nous allons dans le menu Raster ‣ Extraction ‣ Découper un raster selon une couche de masque… Le menu suivant apparaît (Fig. 231).

alternate text

Fig. 231 Paramétrage du découpage du raster.

Dans le champ Couche source nous sélectionnons le raster à découper. Dans le champ Couche de masquage, nous indiquons la couche vecteur de type polygone qui doit servir de masque. Les champs SCR peuvent normalement rester vides. Si la couche de découpage est plus petite que le raster à masquer nous laissons cocher le champ Faire coïncider l'emprise du raster découpé avec l'emprise du calque de masque. Grâce à cette option le raster sera redimensionné aux dimensions du masque. Par contre, si la couche de découpage est plus grande que le raster à découper nous décochons cette case sinon le raster sera « agrandi » par l’ajout de pixels en « no data » aux dimensions de la couche de masque. Il peut être judicieux de définir la valeur que prendront les pixels en « no data » c’est-à-dire ceux qui seront en dehors du masque. En effet, par défaut ils seront mis à 0 ce qui n’est pas toujours pertinent. Cette valeur peut, par exemple, être mise à -99 pour vraiment différencier les « no data » du reste. Enfin dans le champ Découpé nous indiquons un chemin pour sauver le raster masqué qui sera créé. Il n’y plus qu’à cliquer sur Exécuter. Le nouveau raster apparaît automatiquement.

Astuce

Si nous découpons un raster multi-bandes, toutes les bandes sont découpées et nous obtenons bien un raster multi-bandes masqué.

Nous pouvons noter que ce processus de découpage dans QGIS n’est en fait (comme souvent dans les prétraitements de rasters) qu’un interfaçage d’une commande GDAL. La commande sous-jacente s’affiche dans le bloc Console GDAL/OGR en bas de la fenêtre précédente (Fig. 231).

Note

Comme pour la fonction Découper, il est possible de découper plusieurs rasters en lot via la fonction Exécuter comme processus de lot se trouvant en bas de la fenêtre du menu.

Découper un raster avec SCP

Version de QGIS : 3.18.3

Version de SCP : 7.8.21

Le module SCP de QGIS permet également de découper des rasters. Comme à chaque fois que nous travaillons avec ce module, il est nécessaire de commencer par définir un Jeu de bandes de travail. Cela se fait via le menu SCP ‣ Jeu de bandes. Une fois le jeu de bandes défini, le découpage s’effectue en allant dans le menu SCP ‣ Pré-traitement ‣ Découper plusieurs rasters. Le menu suivant s’affiche (Fig. 232).

alternate text

Fig. 232 Découper un raster avec SCP.

Dans le champ Sélectionner un jeu de bandes nous sélectionnons le jeu de bandes à découper. Dans le champ Préfixe de sortie nous indiquons un préfixe qui sera donné aux rasters découpées qui seront créés. Nous cochons ensuite la case Utiliser un vecteur de découpe et nous sélectionnons le vecteur à utiliser. Si celui-ci n’apparaît pas dans la liste déroulante il peut être nécessaire de la rafraîchir avec l’icône icone_refresh. Nous cliquons ensuite sur Lancer. Il nous est demandé de sélectionner un répertoire de sortie et les rasters découpés apparaissent automatiquement.

Avertissement

Les rasters découpés seront créés sous forme de rasters individuels. Même si notre jeu de bandes a été créé avec un raster multi-bandes, les rasters découpés seront tout de même indépendants.

Nous pouvons noter que nous pourrions simplement découper nos rasters selon une emprise sans utiliser de couche de masque. Pour cela nous ne sélectionnons pas l’option Utiliser un vecteur de découpe (Fig. 232) mais nous utilisons l’icône icone_trace pour tracer une emprise sur la fenêtre principale de QGIS.

Découper un raster avec R

Version de R : 4.8.1

Il est possible de découper un raster chargé dans R à l’aide d’un vecteur lui aussi chargé dans R. Dans l’exemple suivant, nous découperons un raster chargé dans R sous forme d’une objet stars par un vecteur contenant un polygone (une limite de département ou de bassin-versant par exemple) chargé sous forme d’un objet sf. Dans cet exemple, le raster stars est nommé bande, le vecteur sf est nommé bassin et nous stockons le raster découpé dans une variable nommé bande_bassin. La commande à utiliser est une commande du package stars et se nomme st_crop. Elle prend simplement en entrée le raster à découper suivi du vecteur à utiliser pour le découpage.

library(stars)
bande_bassin <- st_crop(bande, bassin)

Note

Il est conseillé que les couches rasters et vecteurs soient géoréférencées dans le même système de coordonnées.

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. 233).

alternate text

Fig. 233 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 icone_identifier 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 Raster ‣ Divers ‣ Fusion… Le menu suivant apparaît (Fig. 234).

alternate text

Fig. 234 Paramétrage de la mosaïque.

Dans le champ Couches en entrée nous sélectionnons les rasters à mosaïquer via l’icône icone_browse. 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. 235).

alternate text

Fig. 235 Résultat de la mosaïque de deux images Landsat 8.

Astuce

Sur la figure Fig. 235, 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.

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

Dans QGIS, il existe deux possibilités de concaténation de rasters. La première consiste à créer un vrai raster multi-bandes en dur et la seconde consiste à créer un raster virtuel pointant vers les rasters initiaux.

Raster multi-bandes en dur

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 Raster ‣ Divers ‣ Fusionner …. Le menu de fusion s’ouvre (Fig. 236).

alternate text

Fig. 236 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 icone_browse (Fig. 237).

alternate text

Fig. 237 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 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 Raster ‣ Divers ‣ Construire un raster virtuel. La fenêtre Construire un raster virtuel s’affiche alors (Fig. 238).

alternate text

Fig. 238 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. 239).

alternate text

Fig. 239 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.

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 SCP ‣ Pré-traitement ‣ Séparer les bandes. Le menu suivant s’affiche (Fig. 240).

alternate text

Fig. 240 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 icone_refresh. 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.

Tableau 9 Classification visée

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 Boîte à outils de traitements ‣ Analyse raster ‣ Rapport sur les valeurs uniques de la couche raster. Le menu suivant apparaît (Fig. 241).

alternate text

Fig. 241 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. 242).

alternate text

Fig. 242 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 SCP ‣ Post-traitement ‣ Rapport de classification. 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. 243).

alternate text

Fig. 243 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. 244).

alternate text

Fig. 244 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.

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. 245). Sur cette figure, nous considérons le cas d’un raster composé de 9 pixels, en 3 lignes et 3 colonnes.

alternate text

Fig. 245 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 Raster ‣ Calculatrice Raster. La calculatrice s’ouvre (Fig. 246).

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. 246). Nous cliquons ensuite sur OK.

alternate text

Fig. 246 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. 247 B).

alternate text

Fig. 247 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. 247 B).

alternate text

Fig. 248 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) dans le panneau d’expression en faisant bien attention aux parenthèses (Fig. 249).

Avertissement

Les différents rasters qui entrent dans le calcul doivent être dans le même SCR.

alternate text

Fig. 249 Calcul de NDVI avec la calculatrice raster de QGIS.

Nous obtenons bien le raster de NDVI de la zone (Fig. 250).

alternate text

Fig. 250 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 SCP ‣ Calcul de bande et la calculatrice de bande de SCP s’ouvre (Fig. 251). 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 icone_refresh. 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.

alternate text

Fig. 251 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. 252).

alternate text

Fig. 252 Différence d’altitude entre chaque pixel et l’altitude moyenne de la région.

Sur le raster résultat (Fig. 252), 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. 253).

alternate text

Fig. 253 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 Affichage ‣ Navigateur d’OTB-Applications. 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. 254).

alternate text

Fig. 254 Calcul raster avec l’application BandMath de OTB.

Dans le panneau Name, grâce à l’icône en forme de « + » icone_otb_plus 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)

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 bien se placer dan 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 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 que les rasters ont été chargés sous forme d’objets raster ou stars. Nous verrons ici le cas pour des rasters chargés sous forme d’objets stars.

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 (stars) 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 de pixels 1000 fois plus grandes. Ce raster sera aussi un objet stars que nous stockerons dans une variable nommée pir_1000.

pir_1000 <- 1000 * proche_infrarouge

Nous constatons que l’opération est très simple.

Entre plusieurs rasters

Il est également possible d’effectuer des opérations entre deux, ou plus, 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 stars. Nous stockerons le résultat dans une variable nommée ndvi.

ndvi <- (pir- rouge) / (pir + rouge)

La formulation est très simple. Le raster résultat sera également un objet stars.

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 un raster de NDVI (préalablement chargé au format stars) afin d’extraire les zones les plus végétalisées. Ces zones correspondront aux pixels présentant un NDVI supérieur à 0.8. Le calcul à entrer est alors le suivant.

zones_vegetalisees <- ndvi > 0.8

Aucune difficulté dans la formulation, et le raster correspondant est bien un raster binaire dans lequel les pixels à 1 correspondent aux pixels ayant un NDVI supérieur à 0.8, soit les zones les plus végétalisées.

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. 255).

alternate text

Fig. 255 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 SCP ‣ Jeu de bandes 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 SCP ‣ Calcul de bandes. Le menu suivant s’affiche (Fig. 256).

alternate text

Fig. 256 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. 257).

alternate text

Fig. 257 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 icone_otb_plus. 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)

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.

Tableau 10 Occupation du sol initiale

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.

Tableau 11 Reclassification (regroupement) de l’occupation du sol

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 Boîte à outils de traitements ‣ Analyse raster ‣ Reclassification par table. Le menu s’ouvre (Fig. 258).

alternate text

Fig. 258 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 icone_browse en bout de ligne. Le menu de paramétrage de la table de reclassification s’ouvre (Fig. 259).

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. 259).

alternate text

Fig. 259 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. 260).

alternate text

Fig. 260 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. 261).

alternate text

Fig. 261 Occupation du sol initiale en 7 classes (A) et occupation du sol après regroupement en 4 classes (B).