Auteur : Paul Passy

Licence :

Combinaison de rasters et vecteurs

Il est possible de croiser des couches rasters et vecteurs afin d’en tirer une nouvelle information. Les deux formats dialoguent très bien et de nombreux outils permettent de croiser ces deux types de données. Nous en verrons ici quelques un comme les statistiques zonales ou la conversion d’un raster vers un vecteur et vice versa.

Avertissement

Pour tous les outils qui font intervenir une superposition de rasters et de vecteurs, il est vivement conseillé de travailler avec des couches partageant le même système de coordonnées de référence.

Table des matières

Découper un raster selon une couche vecteur
Statistiques zonales
Histogramme zonal
- Histogramme zonal avec QGIS
Échantillonnage de raster selon des points
- Échantillonnage avec QGIS
- Échantillonnage avec R
Raster vers vecteur
Découper un raster selon une grille
- Découper un raster selon une grille dans QGIS

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 avec un vecteur 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. 281). Notons que les deux couches doivent être dans le même SCR.

alternate text — Fig. 281 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. 282).

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

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 un vecteur dans 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. 283).

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 . 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. 283) mais nous utilisons l’icône 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 selon une couche vecteur. Ici, nous prendrons l’exemple d’une bande spectrale Landsat (LC08_L2SP_196030_20190613_20200828_02_T1_SR_B5.TIF) chargé sous forme d’un objet SpatRaster issu du package terra. Cette bande correspond à une scène Landsat prise au-dessus du sud de la France. Nous allons la déocuper selon les limites du département des Bouches-du-Rhône (dep_13_buffer_32631.gpkg). Ces limites sont un polygone que nous chargerons sous forme d’un objet SpatVector. La manipulation s’effectue en deux temps. D’abord nous découpons le raster selon l’emprise de la couche vecteur et ensuite nous masquons ce sous raster selon les limites du polygone présent sur la couche vecteur. Notons que le raster et le vecteur de découpage doivent être dans le même système de coordonnées de référence. La mnipulation est présentée dans les lignes qui suivent.

Avertissement

La fonction mask() de terra prend en entrée un objet SpatRaster et un objet SpatVector. Si le vecteur de découpage est chargé au format sf, il sera nécessaire de le transformer en un objet SpatVector à l’aide de la fonction vect() de terra.

# chargement du raster à découper
b5 <- terra::rast('Landsat_13/LC08_L2SP_196030_20190613_20200828_02_T1_SR_B5.TIF')

# chargement du polygone vecteur à utiliser pour le découpage
dep13 <- terra::vect('dep_13_buffer_32631.gpkg')

# découpage du raster selon l'emprise du vecteur
b5_emprise_dep13 <- terra::crop(b5, dep13)

# masque de ce sous raster selon le polygone du département 13
b5_dep13 <- terra::mask(b5_emprise_dep13, dep13)

# affichage du résultat
plot(b5_dep13, col=gray.colors(100, start = 0.1, end = 0.9))

Nous obtenons bien la bande spectrale initiale mais masquée selon les contours du département des Bouches-du-Rhône (Fig. 284).

Note

Dans le cas d’un raster multi-bandes, la procédure est exactement la même.

Statistiques zonales 

Cette partie présente la manière de procéder pour calculer des statistiques zonales. Par statistiques zonales nous entendons le fait de calculer des statistiques comme une moyenne, une médiane, un écart-type … à partir d’une couche raster et pour chaque entité d’une couche vecteur de type polygones. Le principe est de récupérer les valeurs de tous les pixels se trouvant sous un polygone et d’en calculer une statistique descriptive comme la moyenne ou autre (Fig. 285).

Sur la figure précédente (Fig. 285), nous pouvons par exemple calculer une statistique zonale pour le polygone vert. Le principe est dans un premier temps de repérer tous les pixels recouverts par le polygone vert lorsque la couche vecteur de polygones est superposée à la couche raster. Puis dans un second temps de calculer la statistique souhaitée sur ces pixels. Cette statistique peut être le maximum, le minimum, la moyenne, la variance, un percentile… Il est également possible, en théorie, de calculer des statistiques zonales selon le même principe mais pour une couche vecteur de type lignes. Dans ce cas aussi, les pixels recouverts par la ligne seront repérés puis une statistique sera calculée sur ces pixels. Notons que les statistiques zonales sous des lignes ne sont pas calculables avec tous les outils.

Note

Les vecteurs de type polygones, et à fortiori de type lignes, ne recouvrent pas tous les pixels à 100 %. En bordure, des pixels se retrouvent recouverts par le polygone en partie seulement. Faut-il alors les compter dans les statistiques zonales ? La question reste entière et il n’est pas toujours simple de savoir ce que fait chaque algorithme selon les outils. Dans la plupart des cas, nous fermons les yeux sur cette question…

Cette manipulation est assez fréquente en SIG. Par exemple, si vous disposez d’un raster de modèle numérique de terrain (MNT) et d’une couche vectorielle présentant des polygones de bassins-versants, grâce à cette manipulation vous pourrez calculer l” altitude moyenne / médiane / maximum … pour chacun de vos bassins. Dans cet exemple, nous allons calculer l’altitude moyenne pour chaque bassin-versant unitaire du bassin de la Roya dans les Alpes-Maritimes. Pour effectuer ce calcul, nous disposons du MNT de la Roya (srtm_roya.tif) et d’une couche vecteur de polygones (hydro_bv_roya.gpkg).

Cette manipulation peut s’effectuer avec tous les logiciels de SIG ainsi qu’avec R ou Python.

Statistiques zonales avec QGIS 

Version de QGIS : 3.16.1

Avec des polygones 

Il faut tout d’abord charger dans QGIS votre MNT et votre couche des bassins-versants. La manipulation s’effectue alors simplement en cherchant le module nommé Statistiques de zone dans la Boîte à outils de traitements de QGIS. À l’ouverture de ce module, la fenêtre de configuration suivante apparaît (Fig. 286).

Dans le champ Couche source, il faut définir la couche vecteur sur laquelle nous voulons calculer les statistiques, à savoir hydro_bv_roya. Dans le champ Couche raster, nous spécifions le raster depuis lequel nous calculons les statistiques, ici srtm_roya. Notons que dans le cas d’un raster multi-bandes, comme une image satellite multispectrale par exemple, nous pouvons choisir la bande à étudier dans le champ Bande raster. Dans le champ Préfixe de la colonne en sortie, nous pouvons mettre alti pour altitude. Le underscore va permettre de nommer le champ résultant alti_moyenne.

Dans le champ Statistiques à calculer, nous cliquons sur l’icône et nous sélectionnons la statistique que nous désirons parmi toutes celles proposées, moyenne dans notre cas. Enfin, dans le champ Statistiques de zone, nous spécifions le chemin d’export de la couche vecteur qui contiendra les statistiques. En effet, une nouvelle couche sera créée, il n’est pas possible de mettre directement à jour la couche. La nouvelle couche apparaît automatiquement et elle contient un nouveau champ nommé alti_mean qui contient l’altitude moyenne de chaque bassin-versant.

Avec des polylignes 

Il est possible de calculer des statistiques zonales le long de polylignes dans QGIS en utilisant un outil issu de GRASS GIS. Comme les outils GRASS sont interfacés nativement par QGIS, il n’y a aucune difficulté. Dans l’exemple, nous allons calculer quelques statistiques zonales d’altitude le long du réseau hydrographique du bassin de la Roya (Roya_hydro.gpkg). Nous commençons par charger la couche du réseau hydrographique sous forme de polylignes. Dans le fichier employé, chaque tronçon de rivière compris entre une source et la première confluence ou entre deux confluences est considéré comme une entité. Une fois les couches chargées, nous ouvrons le menu v.rast.stats qui se trouve dans la Boîte à outils de traitements ‣ GRASS ‣ Vector (v.*). La fenêtre suivante s’ouvre (Fig. 287).

À la ligne Name of vector polygon map nous indiquons le nom de la couche de points pour laquelle nous souhaitons les statistiques, à savoir Roya_hydro. Bien qu’il soit mentionné Polygon, il est possible de mettre des lignes en entrée. À la ligne Name of raster map to calculate statistics from, nous indiquons le raster à utiliser pour le calcul des statistiques, à savoir srtm_roya. La ligne Column prefix for new attribute columns est pratique. Il s’agît de définir le préfixe que prendront les colonnes de statistiques, ici hydro. Dans le panneau The methods to use, nous précisons les statistiques à calculer en cliquant sur l’icône . Pour l’exemple nous choisissons la moyenne (average) et l’écart-type (stdev). Enfin, à la ligne Rast stats nous indiquons un chemin et un nom pour la nouvelle couche de type lignes qui contiendra deux nouvelles colonnes dans la table attributaire, une pour les altitudes moyennes et une autre pour l’écart-type des altitudes de chaque tronçon de rivière. Nous pouvons la nommer Roya_hydro_stats_alti.gpkg. Puis nous cliquons sur Exécuter. La nouvelle couche, en tout point identique à la couche de lignes initiale, apparaît. Mais nous y retrouvons bien deux nouvelles colonnes dans la table attributaire.

Statistiques zonales sur raster multi-bandes avec QGIS 

Version de QGIS : 3.34.0

Il est possible de calculer des statistiques zonales sur plusieurs rasters à la fois en calculant ces statistiques sur un raster multi-bandes. Ce raster multi-bandes peut être un « vrai » multi-bandes au format .tif ou un bien un raster virtuel au format .vrt.

Note

Ce module ne prend en entrée que des couches vecteurs de type polygones.

Dans cet exemple, nous allons calculer les réflectances moyennes sur 7 bandes spectrales Landsat 8 pour quelques communes du Val-d’Oise.

Cette fonctionnalité nécessite d’installer au préalable une extension supplémentaire nommée Zonal Statistics for Multiband Rasters. Cette extension n’est pas disponible dans les dépôts officiels de QGIS, mais dans un dépôt maintenu par la société Dymaxion labs. Comme indiqué sur la page de documentation du plugin, il suffit d’aller dans le menu Extensions ‣ Installer/Gérer les extensions ‣ Paramètres. La fenêtre des paramètres des extensions apparaît (Fig. 288).

Dans ce menu, nous commençons par cocher la case Afficher les extensions expérimentales puis nous cliquons sur Ajouter…. Dans la fenêtre qui apparaît, dans le champ Nom, nous mettons un nom libre, par exemple Dymaxion. Dans le champ URL, nous renseignons l’URL du dépôt https://dymaxionlabs.github.io/qgis-repository/plugins.xml puis nous cliquons sur OK. Sur la fenêtre précédente, il n’y a plus qu’à cliquer sur Recharger tous les dépôts, à retourner dans l’onglet Toutes et à rechercher l’extension Zonal Statistics for Multiband Rasters.

Une fois l’extension installée, elle apparaît comme un module dans la Boîte à outils de traitements` sous l’entrée Dymaxion Labs. Nous chargeons ensuite le raster multi-bandes à analyser. Il peut être en .tif ou bien sous la forme d’un raster virtuel. Nous chargeons également la couche des entités polygonales sur lesquelles nous souhaitons calculer les statistiques.

Nous ouvrons cette extension en allant dans la Boîte à outils de traitements ‣ Dymaxion Labs ‣ Raster Analysis ‣ Zonal statistics (multiband). Le menu suivant apparaît (Fig. 290).

Dans le champ Raster layer nous indiquons le raster multi-bandes à analyser, ici L8_stack. Dans le champ Vector layer containing zones, nous sélectionnons la couche vecteur de polygones, ici communes_95_str. Dans le champ Output column prefix nous pouvons spécifier un préfixe pour le nouveau champ qui sera calculé. Comme nous allons calculer la moyenne des réflectances, nous le nommons moy_ . Enfin, dans le champ Statistics to calculate nous sélectionnons une ou plusieurs statistiques à calculer. Ici, nous ne calculons que la moyenne. Nous cliquons ensuite sur Exécuter.

Contrairement à l’outil interne à QGIS, les statistiques sont directement ajoutées à la couche vecteur d’entrée. Si nous ouvrons sa table attributaire, nous remarquons les nouveaux champs calculés à la fin. Ils se nomment moy_b1_mean, moy_b2_mean, etc

Astuce

Ce module peut s’avérer très pratique. Il peut même être pertinent de créer un raster multi-bandes spécialement pour utiliser cette extension. Par exemple, lorsque nous souhaitons calculer des primitives pour classifier des segments suite à une segmentation d’images, cette extension peut nous faire gagner du temps. Nous pouvons également l’utiliser sur un raster mono-bande si nous souhaitons « mettre à jour » une couche vecteur avec les statistiques zonales sans devoir en créer une nouvelle.

Statistiques zonales avec OTB 

Version de OTB : 7.2.0

Orfeo ToolBox propose un module de calcul de statistiques zonales. Comme tous les modules OTB, il est utilisable via l’interface Monteveri ou via la ligne de commende otbcli. Dans l’exemple nous allons calculer l’altitude moyenne de quelques communes localisées dans le bassin-versant de la Roya.

Avertissement

OTB ne permet pas de calculer des statistiques zonales sur autre chose que des polygones.

Via Monteverdi 

Dans l’interface, le module dédié se nomme ZonalStatistics et se trouve dans le Navigateur d’OTB-Applications ‣ Image Analysis ‣ ZonalStatistics. Une fois sélectionné, le menu suivant s’ouvre (Fig. 291).

À la ligne Input Image, nous sélectionnons le raster sur lequel calculer les statistiques, à savoir srtm_roya.tif. Dans le panneau Type of input for the zone definitions, nous spécifions si nos zones sont sous forme de fichier vecteur ou raster. Ici, nous travaillons avec une couche vecteur de type polygones représentant les communes, nous sélectionnons donc Input objects from vector data ; et à la ligne Input vector data nous indiquons où se trouve le fichier vecteur à utiliser. Si jamais le fichier vecteur n’est pas dans le même système de coordonnées de référence que la couche raster d’entrée, nous pouvons la reprojeter à la volée en mettant à ON l’option Reproject the input vector. Nous devons ensuite choisir un type de sortie dans le panneau Format of the output stats. Si nous choisissons une sortie en vecteur, nous sélectionnons Output vector data et nous le nommons à la ligne Filename for the output vector data. La sortie peut être en geopackage. Puis nous cliquons sur Execute. Un fichier vecteur en tout point identique au fichier d’entrée est créé. Mais il y a maintenant cinq nouvelles colonnes dans la table attributaire :

count : le nombre de pixels sous chaque polygone

mean_0 : la moyenne des pixels sous chaque polygone et dans la première bande du raster d’entrée (qu’une bande ici)

stdev_0 : l’écart-type des pixels

min_0 : le minimum des pixels

max_0 : le maximum des pixels

Il n’y a pas le choix des statistiques calculées avec OTB.

Astuce

Si le raster en entrée est un Raster multi-bandes, les statistiques zonales sont calculées sur toutes les bandes dans le même ordre que dans le raster. Les colonnes de statistiques ont alors un suffixe _0, _1, …

Ce module propose également de sortir le résultat sous forme d’un raster plutôt que d’un vecteur. Le résultat est alors un raster multi-bandes où chaque bande est une statistique et dans le même ordre que présenté dans le cas d’une sortie vecteur. Cette sortie sous forme de raster peut être utile dans le cas où cette sortie sert d’entrée à une classification.

Via la ligne de commande 

Le calcul des statistiques zonales est également possible via la ligne de commande otbcli. Pour une documentation détaillée, voire la référence sur le site de OTB. La ligne de commande idoine prend la forme suivante.

otbcli_ZonalStatistics -in srtm_roya.tif -inzone.vector.in roya_communes.gpkg -out.vector.filename roya_communes_alti.gpkg

Avec :

otbcli_ZonalStatistics : l’appel à la fonction dédiée
-in srtm_roya.tif : le raster duquel calculer les statistiques à utiliser en entrée
-inzone.vector.in roya_communes.gpkg : les zones à utiliser, ici sous forme de vecteur
-out.vector.filename roya_communes_alti.gpkg : le fichier de sortie, ici sous forme de vecteur

Statistiques zonales avec R 

Il est possible de calculer des statistiques zonales avec R. Il est nécessaire d’avoir au préalable chargé une couche raster et une couche vecteur avec des librairies dédiées. Dans les exemples ci-dessous, nous allons calculer l’altitude médiane pour quelques communes du bassin-versant de la Roya. Nous utiliserons le MNT fourni par le SRTM (srtm_roya_L93.tif) et une couche des communes au format geopackage sous forme de polygones (Roya_communes.gpkg).

Avec Terra 

Version de R : 4.3.1

Version de terra : 1.7.29

Nous commençons par charger le raster dans un objet SpatRaster et le vecteur dans un objet SpatVector. Si jamais le vecteur est chargé en sf, il suffit de transformer cet objet sf en SpatVector (Format SpatVector de terra). Le calcul proprement dit se fait à l’aide de la fonction extract de terra.

# chargement du raster et du vecteur
mnt_roya <- terra::rast('srtm_roya_L93.tif')
communes_roya <- terra::vect('communes_bv_Roya.gpkg')

# calcul de l'altitude médiane par commune
communes_alti <- terra::extract(mnt_roya, communes_roya, fun=median, bind=TRUE)

# on renomme le champ créé avec un nom explicite
names(communes_alti)[names(communes_alti) == names(mnt_roya)] <- 'Alti_mediane'

Il faut faire attention à deux choses. Tout d’abord, il ne faut pas oublier le mot clef bind=TRUE de la fonction extract, il permet de récupérer le résultat dans un objet SpatVector avec tous les champs initiaux. Sans ça, nous récupérons un simple dataframe duquel les champs initiaux ont disparu. Ensuite, le champ de statistique créé porte par défaut le nom du raster, ici srtm_roya_L93, ce qui n’est pas explicite. C’est pour cette raison, qu’il est bien de le renommer avec la fonction names.

Histogramme zonal 

Cette partie présente comment calculer ce que nous appelons un histogramme zonal. Contrairement aux statistiques zonales, il ne s’agira pas de calculer une moyenne ou un maximum d’une couche raster sous un polygone mais de compter le nombre d’occurrences de chaque valeur différente d’une couche raster sous un polygone. Par exemple, sur la figure suivante, nous avons superposé une couche vecteur de polygones de départements à une couche d’occupation du sol de l’Île-de-France (Corine Land Cover). Le calcul d’un histogramme zonal va nous permettre de compter le nombre de pixels de chaque valeur d’occupation du sol sous chaque département. Nous connaîtrons ainsi, par exemple, le nombre de pixels ayant le code 313 (forêt mixte) sous chaque département francilien. Nous pourrons ensuite, si nous le souhaitons, connaître la superficie de chaque occupation du sol en multipliant le nombre de pixels trouvé par la taille du pixel.

Cette manipulation est très utile, notamment lorsque nous souhaitons rester en mode raster et éviter de convertir des données en vecteur pour calculer des superficies.

Avertissement

Le calcul d’un histogramme zonal doit se faire sur un raster discret, i.e. qui ne contient qu’un nombre limité de valeurs possibles (comme un raster d’occupation du sol par exemple). Si un histogramme zonal est calculé depuis un raster continu comme une bande spectrale, le calcul plantera certainement du fait du trop grand nombre de valeurs différentes rencontrées.

Histogramme zonal avec QGIS 

Version de R : 3.24.3

Le calcul d’un histogramme zonal dans QGIS est très simple. Il se fait à l’aide du menu Histogramme zonal qui se trouve dans la Boîte à outils de traitements ‣ Analyse raster ‣ Histogramme zonal. Le menu suivant s’ouvre (Fig. 293).

Dans le menu Couche raster, nous spécifions la couche raster à partir de laquelle nous souhaitons calculer l’histogramme, à savoir CLC_2018_IDF.gpkg. Si nous avions affaire à un raster multi-bandes, nous pourrions choisir la bande à traiter via le menu Numéro de bande. Puis, dans le menu Couche vecteur contenant les zones, nous indiquons la couche de polygones à utiliser, ici departements_IDF. Nous pouvons changer le Préfixe de la colonne en sortie par un préfixe explicite, comme CLC_ (pour Corine Land Cover) dans notre cas. Enfin, à la ligne Zones de sortie, nous indiquons un chemin et un nom pour la couche vecteur de polygones qui sera créée.

À la fin de ce calcul, la nouvelle couche s’affiche. Elle est géographiquement identique à la couche de polygones de départ, mais sa table attributaire a été amendée avec le compte de chaque occurrence de pixel (Fig. 294). Dans cette table, nous retrouvons les champs attributaires de la couche initiale plus les nouveaux champs, ici commençant par CLC_, le préfixe que nous avions choisi précédemment. Ainsi, par exemple, pour le département 78, il y a 205 pixels classés en 205.

Échantillonnage de raster selon des points 

Il est possible de superposer une couche vecteur ponctuelle sur une couche raster et de relever les valeurs des pixels recouverts par les points. Ainsi, nous ajoutons dans la table attributaire de la couche de points les valeurs des pixels du raster sous-jacents à chacun des points. Cette manipulation peut se faire selon différents outils. Cette manipulation est aussi connue sous son appellation en anglais Point sampling tool.

Dans les exemples suivants, nous releverons l’altitude à partir d’un MNT pour quelques points superposés à ce MNT (Fig. 295).

Échantillonnage avec QGIS 

Version de QGIS : 3.22.0

Nous commençons par charger dans QGIS le raster à échantillonner et la couche de points pour lesquels nous souhaitons relever les échantillons. Le menu à utiliser se trouve dans la Boîte à outils de traitements ‣ Analyse raster ‣ Prélever des valeurs rasters. Le menu suivant s’ouvre (Fig. 296).

Dans le menu Couche source nous spécifions la couche de points à utiliser. Nous pouvons, si nous le souhaitons, cocher la case Entité(s) séléectionnée(s) uniquement si nous souhaitons simplement échantillonner certains points précédemment sélectionnés. Dans le menu Couche raster nous indiquons le raster à échantillonner. Nous pouvons personnaliser le Préfixe de la colonne en sortie, ici alti_ par exemple. Enfin, dans le champ Échantillonné nous indiquons un chemin et un nom pour la couche de points qui contiendra les valeurs du raster échantillonné.

Nous obtenons une couche de points identiques à la couche de points initial mais avec un nouveau champ intitulé alti_1 où nous trouvons la valeur du pixel du MNT sous-jacent à chaque point.

Note

Si notre MNT avait été multi-couches, nous aurions obtenu autant de nouveaux champs que de couches dans le raster, nommés alti_1, alti_2, …

Échantillonnage avec R 

Version de R : 4.3.1

Version de terra : 1.7.29

Il est possible d’échantillonner un raster selon une couche de points dans R en utilisant la librairie terra et sa fonction extract. Cela se fait simplement comme présenté ci-dessous.

# chargement du raster et du vecteur de points
mnt_roya <- terra::rast('srtm_roya_L93.tif')
points_bv <- terra::vect('points_bv_Roya.gpkg')
# extraction des valeurs
mnt_points <- terra::extract(mnt_roya, points_bv, bind=TRUE)

Nous ajoutons le paramètre bind=TRUE afin de récupérer le résultat sous forme d’un objet SpatVector. Si nous mettons ce paramètre à False, nous obtenons le résultat sous forme d’un simple tableau.

Raster vers vecteur 

Il est possible de convertir un raster en un vecteur. C’est une manipulation notamment utilisé pour transformer un raster d’occupation du sol en une couche vecteur sur laquelle il sera possible de calculer des superficies ou de faire des croisements avec d’autres couches géospatiales. C’est également utilisé à la suite de la construction d’un raster binaire mettant en avant un aspect particulier du terrain. La manipulation peut se faire selon différents outils qui peuvent proposer plus ou moins d’options. Dans les exemples suivants nous convertirons un (petit) raster d’occupation du sol.

Seuls certains outils proposent de convertir autre chose que des surfaces, comme des lignes par exemple. C’est le cas de l’outil proposé par GRASS GIS (Cas pour des lignes).

Avertissement

Après une conversion de couche raster en couche vecteur, même si le résultat est visuellement propre, il peut y avoir des erreurs de topologie dans le vecteur résultat. Ces erreurs empêcheront d’effectuer des géotraitements à partir de ce vecteur. Il sera nécessaire de corriger les erreurs générées en port-traitement. Ces corrections peuvent se faire dans QGIS avec l’outil Réparer les géométries.

Raster vers vecteur dans QGIS 

Version de QGIS : 3.20.1

Une fois le raster chargé dans QGIS, il est possible de le transformer une couche vecteur. Pour ce faire, nous allons dans le menu Raster ‣ Conversion ‣ Polygoniser (raster vers vecteur)… Le menu dédié s’ouvre (Fig. 297).

Dans la ligne Couche source, nous indiquons le raster à convertir, ici xtr_usol_2011. Nous pouvons choisir la bande à convertir dans le menu déroulant Numéro de bande, mais ici nous n’avons qu’une bande dans notre raster. Nous indiquons le nom du champ qui sera créé et qui contiendra les valeurs des pixels à la ligne Nom du champ à créer, nous pouvons indiquer usol. La case Use 8-connectedness permet de regrouper dans un même polygone deux entités qui ne se toucheraient que par un coin. Si cette case est décochée, deux entités se touchant par un coin seront considérées comme distinctes. Enfin, dans la ligne Vectorisé nous indiquons un chemin et un nom pour le fichier vecteur qui sera créé. Notons que ce menu est en fait un simple interface d’une commande GDAL comme nous pouvons le voir dans le panneau Console GDAL/OGR. Nous cliquons sur Exécuter pour lancer la conversion. Nous obtenons bien un fichier vecteur de type polygones (Fig. 298 B) possédant une table attributaire associée.

Note

Nous constatons qu’avec ce menu il n’est possible d’obtenir que des vecteurs de type polygones (i.e. surfacique), il n’est pas possible de vectoriser des lignes comme un réseau hydrographique par exemple.

Raster vers vecteur avec SCP 

Version de QGIS : 3.20.1

Version de SCP : 7.9.5

Le module SCP de QGIS (Module SCP dans QGIS) propose lui aussi un menu de conversion de raster vers vecteur. Une fois le raster chargé dans QGIS, nous allons dans le menu SCP ‣ Post-traitement ‣ Classification vers vecteur. Le menu s’ouvre (Fig. 299). Dans le panneau Clasification vers vecteur, nous sélectionnons le raster à convertir. Si besoin nous rafraîchissons la liste à l’aide de l’icône dédiée . Si le raster à convertir est issu d’une classification réalisé avec SCP nous pouvons récupérer directement son style en cochant la case Utiliser les codes de la liste des signatures. Nous cliquons ensuite sur Lancer.

Nous constatons que le fichier vecteur créé est découpé en bandes horizontales. C’est un artefact utilisé par ce module pour accélérer le traitement. Il suffit de regrouper les polygones par classe d’occupation du sol pour s’affranchir à postériori de cet artefact.

Raster vers vecteur avec GRASS GIS dans QGIS 

Version de QGIS : 3.20.1

Le logiciel GRASS GIS (GRASS GIS) propose un module de transformation de raster vers vecteur. Comme GRASS GIS est interfacé en natif par QGIS, ce module est utilisable directement depuis QGIS. Une fois le raster chargé, nous sélectionnons le menu r.to.vect dans la Boîte à outils de traitements ‣ GRASS ‣ Raster ‣ r.to.vect. Le menu GRASS apparaît (Fig. 300).

Cas pour des polygones 

À la ligne Couche raster en entrée nous indiquons le raster à convertir, ici xtr_usol_2011. Dans le menu Type d'entités nous spécifions que nous allons convertir des surfaces (i.e. des polygones) en sélectionnant area. Dans le champ Nom de la colonne d'attribut nous pouvons indiquer le nom que prendra le champ qui contiendra les valeurs prises du raster, usol par exemple. Nous pouvons laisser les cases suivantes décochées. Enfin, dans Vectorisé, nous indiquons un chemin et un nom pour le raster qui sera créé. Puis nous cliquons sur Exécuter. Le fichier converti apparaît automatiquement. Plusieurs champs sont créés, mais seul le champ que nous avons baptisé usol nous importe. Les autres champs peuvent être supprimés.

Cas pour des lignes 

Un aspect intéressant de cet outil de GRASS GIS est qu’il permet de vectoriser des lignes. C’est tout à fait intéressant dans le cas où nous avons un raster présentant un réseau hydrographique, par exemple issu d’une extraction automatique de réseau hydrographique à partir d’un modèle numérique de terrain. Le point délicat d’une telle conversion de lignes est que le raster de départ doit rigoureusement représenter des lignes. C’est-à-dire que chaque pixel doit toucher au plus deux pixels (et pas un de plus). Heureusement, GRASS propose également un outil qui permet d’amincir un raster de lignes nommé r.thin. Ainsi, la conversion d’un raster de lignes va se faire en deux temps. D’abord, l’amincissement avec r.thin puis la conversions proprement dite avec r.to.vect. Dans l’exemple, nous allons convertir en vecteur le réseau hydrographique de la région de la Roya que nous avons en raster. Nous commençons par charger ce raster dans QGIS puis nous ouvrons le menu r.thin via la Boîte à outils de traitements ‣ GRASS ‣ Raster ‣ r.to.thin. Le menu suivant apparaît (Fig. 301).

Dans le menu Input raster layer to thin nous sélectionnons le raster à amincir hydro_Roya. Le Nombre maximum d'itérations doit être d’autant plus élevé que le raster a besoin d’être aminci, la valeur de 200 peut être laissée par défaut. Enfin, à la ligne Aminci, nous indiquons le chemin et le nom du raster aminci qui sera créé. Il ne reste plus qu’à cliquer sur Exécuter. Le raster résultat est très similaire au raster initial, seuls quelques pixels ont été enlevés (Fig. 302).

Nous constatons sur la figure précédente (Fig. 302) que le pixel en rouge du raster initial touchait 3 autres pixels. Il a donc été enlevé par la procédure d’amincissement, ainsi chaque pixel ne touche au plus que 2 pixels. Nous pouvons maintenant lancer la procédure de conversion en vecteur en lançant le menu r.to.vect (Fig. 303).

Dans le menu Couche raster en entrée nous indiquons le raster aminci, ici hydro_thin. Dans Type d'entités nous précisons qu’il s’agît de lignes (line). Nous donnons éventuellement un nom pour le champ qui contiendra les valeurs du raster dans le menu Nom de la colonne d'attribut et enfin dans Vectorisé, nous précisons un chemin et un nom pour le fichier vecteur de type lignes qui sera créé. Pour finir nous cliquons sur Exécuter. Nous avons bien un vecteur de type lignes qui apparaît.

Raster vers vecteur avec GDAL 

Il est possible de faire une conversion de raster en vecteur avec GDAL en ligne de commande. Cette procédure est rapide et peut être intéressante dans une optique d’automatisation. Il est juste nécessaire de faire attention à la syntaxe qui est somme toute assez simple. Il faut, dans une console, se placer dans le répertoire qui contient le raster à convertir et entrer la commande suivante :

gdal_polygonize.py xtr_usol_2011.tif -f "GPKG"  xtr_usol_2011_gdal.gpkg usol_2011 classe

avec les éléments suivants :

gdal_polygonize.py : le nom de la commande GDAL à utiliser
xtr_usol_2011.tif : le nom du raster à convertir
-f « GPKG » : le format de la couche vecteur à créer. Mettre « SHP » pour créer un shapefile.
xtr_usol_2011_gdal.gpkg : le nom du fichier vecteur qui sera créé
usol_2011 : le nom de la couche qui sera créé dans le geopackage. Cette option n’est possible que pour les formats vecteur proposant de stocker plusieurs couches (comme le geopackage).
classe : le nom du champ qui contiendra les valeurs du raster

Note

Si nous ne nous plaçons pas dans le répertoire qui contient le raster à convertir, nous devons pointer explicitement sur ce raster en inscrivant son chemin dans la ligne de commande, par exemple C:\Test\Truc\mon_raster.tif sous Windows ou /home/SIG/truc/mon_raster.tif sous Linux. Faire de même pour l’export du vecteur si désiré.

Raster vers vecteur avec R 

Version de R : 4.8.1

Il est possible de convertir un raster vers un vecteur en utilisant R. La méthode à employer sera différente selon le package raster employé. Dans cet exemple, nous travaillerons sur un objet SpatRaster chargé via terra. Le raster chargé est un MNT du bassin de la Roya issu de la mission SRTM. À partir de ce MNT, nous allons extraire les pixels dont l’altitude est supérieure à 1500 m en appliquant un seuil, puis nous transformerons le résultat de ce raster binaire en une couche vecteur. Cette transformation d’un raster en un vecteur se fait via la fonction terra nommée as.polygons(). Il suffit de changer le mot clé polygons par lines ou points si il s’agît de vectoriser des lignes ou des points respectivement. Le code suivant présente la marche à suivre.

# chargement du raster de MNT
mnt_roya <- terra::rast("srtm_roya_L93.tif")

# application du seuil pour extraire les pixels d'altitude supérieures à 1500 m
mnt_roya_1500 <- mnt_roya > 1500

# vectorisation du raster binaire obtenu
seuil_1500_vect <- as.polygons(mnt_roya_1500, extent=FALSE)

Nous obtenons un vecteur de type polygone stocké en format SpatVector.

Découper un raster selon une grille 

Lorsque nous disposons d’un raster assez lourd à manipuler il peut être intéressant de le découper en dalles plus petites qui seront plus faciles à manipuler par la suite. Dans les exemples ci-après, nous allons découper une scène Landsat 8 prise au dessus de Shangai selon une grille composée de 25 carreaux (Fig. 304). À l’issue du traitement, nous disposerons ainsi de 25 dalles rasters correspondant exactement à notre image Landsat initiale par assemblage.

Découper un raster selon une grille dans QGIS 

Version de QGIS : 3.24.3

Dans QGIS il n’existe pas (à ma connaissance) d’outils permettant de faire cette manipulation en une seule fois de façon directe. L’astuce est de procéder en trois étapes :

création d’une grille de polygones
séparation de chacun des carreaux de la grille dans une couche indépendante
découpage du raster initial par chacun des carreaux

La première étape se fait en suivant les explications fournies dans la partie dédiée à la construction de grilles (Création de grilles dans QGIS). La deuxième étape se fait également en suivant les explications fournies dans la partie qui traite de la séparation d’une couche vecteur (Séparer une couche vecteur dans QGIS). La troisième partie est expliquée dans les lignes qui suivent.

L’idée est maintenant de découper notre raster selon les différents carreaux constitutifs de notre grille qui sont maintenant individualisés (Fig. 305).

Il serait fastidieux de répéter manuellement la manipulation de découpage du raster carreau par carreau. Nous utiliserons donc la possibilité de lancer des processus par lot offerte par QGIS. Pour cela nous ouvrons le menu de découpage d’un raster par un vecteur en allant dans le menu Raster ‣ Extraction ‣ Découper un raster selon une couche de masque… Et au lieu de directement paramétrer la fenêtre qui s’ouvre, nous cliquons sur Exécuter comme processus de lot… en bas de la fenêtre. Une grande fenêtre s’ouvre alors. Dans cet interface, chaque ligne correspondra à un processus à traiter ligne par ligne.

Seulement trois colonnes devront être renseignées : la première nommée Couche source, la deuxième nommée Couche de masque et la dernière nommée Découpé (masque). Nous allons commencer par renseigner la deuxième colonne qui renseignera tous nos carreaux à utiliser. Pour tout renseigner en une seule fois, dans cette colonne Couche de masque, nous cliquons sur Auto-remplissage... et dans le menu qui apparaît nous choisissons l’entrée Sélectionner à partir des couches chargées.... Dans le menu qui apparaît nous sélectionnons les couches correspondant à tous les carreaux (Fig. 306).

Une fois les couches choisies il suffit de cliquer sur OK. Autant de lignes que de couches sélectionnées sont créées dans le menu du processus de lot. Nous pouvons maintenant spécifier quel raster nous souhaitons découper, dans la colonne Couche source. Le même raster apparaîtra à chaque ligne puisque c’est bien un seul raster que nous souhaitons découper selon plusieurs carreaux. Dans le menu déroulant de la première ligne, dans la colonne Couche source nous choisissons le raster à découper, à savoir ici L8_Shangai_B5. Afin de le répéter sur toutes les lignes, il suffit de cliquer sur la zone Auto-remplissage... de cette colonne et de choisir l’entrée Remplir. Le même raster d’entrée est alors spécifié pour toutes les lignes (Fig. 307).

Maintenant nous devons régler la dernière colonne Découpé (masque) qui spécifiera le nom et le chemin des dalles créées. Dans cette colonne, à la première ligne, nous cliquons sur l’icône et nous indiquons un répertoire dans lequel nous sauverons les dalles ainsi que le nom générique de chaque dalle, par exemple dalle_ (avec un tiret bas à la fin). Lorsque nous cliquons sur Enregistrer, une fenêtre nous demande si nous souhaitons compléter le nom générique de chaque dalle avec un paramètre permettant son identification individuelle (Fig. 308).

Dans le menu Mode de remplissage automatique nous pouvons sélectionner Remplir avec les valeurs du paramètre et dans le menu Paramètre à utiliser nous pouvons choisir Couche de masque. Ainsi, le morceau de raster découpé avec le carreau nommé id_01 sera nommée dalle_id_01 et ainsi de suite. En cliquant sur OK, chaque ligne de cette colonne est maintenant renseignée comme il faut (Fig. 309).

Il n’y a plus qu’à cliquer sur Exécuter pour lancer le processus. Si nous avons coché la case Charger les couches dans le menu de processus de lot, toutes les dalles sont chargées. Sinon il suffit de les charger depuis le répertoire où elles ont été créées. Nous avons bien maintenant autant de dalles rasters que de carreaux dans la grille initiale (Fig. 310).