Auteur : Paul Passy

Licence : cc_by_nc_sa

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.

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

alternate text

Fig. 284 Principe du calcul de statistiques zonales par recouvrement d’une couche raster (en gris) par une couche vecteur de type polygone en vert et de type ligne en rouge.

Sur la figure précédente (Fig. 284), 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. 285).

alternate text

Fig. 285 Paramétrage des statistiques zonales dans QGIS.

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 icone_browse 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. 286).

alternate text

Fig. 286 Statistiques zonales le long de lignes dans QGIS.

À 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 icone_browse. 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.18.3

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 est expérimentale (mais fonctionne tout de même bien). Pour l’installer, il est nécessaire d’autoriser les extensions expérimentales en cochant la case Afficher les extensions expérimentales dans le menu Extensions ‣ Installer/Gérer les extensions ‣ Paramètres. Une fois cette option cochée, nous retournons dans l’onglet Toutes et nous recherchons cette extension dans la barre de recherche (Fig. 287).

alternate text

Fig. 287 Installation de l’extension de calcul de statistiques zonales multi-bandes.

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

alternate text

Fig. 288 Calcul de statistiques zonales multi-bandes.

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

alternate text

Fig. 289 Calcul de statistiques zonales avec OTB via Monteverdi.

À 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

Version de R : 4.1.2

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 et sf

Pour ce calcul de statistiques zonales, il est nécessaire de charger au préalable la couche raster au format SpatRaster via la librairie Terra et la couche vecteur au format sf. Ensuite, le processus se fait en plusieurs temps. Il est tout d’abord nécessaire de calculer un champ d’identifiants uniques ID à la couche vecteur. En effet, lors du calcul des statistiques zonales, les résultats sont retournés dans un dataframe où chaque entité est simplement repérée via son ID. Dans un deuxième temps, nous calculons la statistique zonale désirée, ce qui nous retourne un simple dataframe. Enfin, nous joignons ce dataframe à la couche vecteur initiale afin de retrouver une couche vecteur.

library(terra)
library(sf)
library(tidyverse) # pour renommer les champs d'un dataframe

# import de la couche de polygones des communes
communes <- sf::st_read('communes_bv_Roya.gpkg')

# import du MNT
mnt_roya <- terra::rast('srtm_roya_L93.tif')

# on ajoute un champ d'ID à la couche des communes
communes <- cbind(ID = 1:nrow(communes), communes)

# conversion de la couche de polygones en SpatVector
communes <- terra::vect(communes)

# on calcule l'altitude médiane par commune (on obtient un simle dataframe)
communes_alti <- terra::extract(mnt_roya, communes, fun=median)

# on renomme le champ qu'on vient de calculer
communes_alti <- rename(communes_alti, 'Altitude médiane' = 'srtm_roya_L93')

# on joint ce df à la couche initiale des communes afin de retrouver une couche SpatVector
communes_alti <- merge(x=communes, y=communes_alti, by="ID")

# facultatif : on transforme ce SpatVector en sf
communes_alti <- sf::st_as_sf(communes_alti)

Comme nous le voyons dans le code ci-dessous, le calcul se fait grâce à la fonction extract() de la librairie Terra. Cette fonction prend en entrée un SpatRaster et un SpatVector, c’est pour cette raison que nous avons converti le vecteur initialement chargé en sf.

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.

alternate text

Fig. 290 Exemple d’application d’histogramme zonal.

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

alternate text

Fig. 291 Calcul d’un histogramme zonal dans QGIS.

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

alternate text

Fig. 292 Résultat d’un calcul d’histogramme zonal.

É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. 293).

alternate text

Fig. 293 Points d’échantillonnage.

É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. 294).

alternate text

Fig. 294 Échantillonnage d’un raster selon une couche de points dans QGIS.

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, …

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

alternate text

Fig. 295 Conversion d’un raster en vecteur dans QGIS.

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. 296 B) possédant une table attributaire associée.

alternate text

Fig. 296 Le raster initial (A) et le résultat de la conversion en vecteur (B).

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. 297). 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 icone_actualiser. 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.

alternate text

Fig. 297 Conversion d’un raster en vecteur avec SCP.

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

Cas pour des polygones

alternate text

Fig. 298 Conversion d’un raster en vecteur avec GRASS GIS dans QGIS.

À 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. 299).

alternate text

Fig. 299 Amincissement d’un raster de lignes avec GRASS GIS dans QGIS.

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

alternate text

Fig. 300 Raster initial (A) et raster après amincissement (B).

Nous constatons sur la figure précédente (Fig. 300) 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. 301).

alternate text

Fig. 301 Conversion d’un raster de lignes en vecteur avec GRASS GIS dans QGIS.

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 de la librairie terra
library(terra)

# 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)

# facultatif : transformation du résultat en objet sf
seuil_1500_sf <- sf::st_as_sf(seuil_1500_vect)

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. 302). À l’issue du traitement, nous disposerons ainsi de 25 dalles rasters correspondant exactement à notre image Landsat initiale par assemblage.

alternate text

Fig. 302 Image Landsat 8 à découper selon la grille rouge.

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

alternate text

Fig. 303 Grille de découpage individualisée carreau par carreau.

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

alternate text

Fig. 304 Sélection des carreaux de découpage.

Astuce

Il peut être plus rapide de choisir Sélectionner tout et d’enlever les couches à ne pas utiliser plutôt que de cliquer sur chacune des couches à utiliser.

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

alternate text

Fig. 305 Remplissage des colonnes Couche source et couche de masque.

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 icone_browse 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. 306).

alternate text

Fig. 306 Remplissage d’un paramètre automatique pour nommer les dalles produites.

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

alternate text

Fig. 307 Remplissage de la colonne d’export.

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

alternate text

Fig. 308 Dalles rasters correspondant aux carreaux de la grille de découpage.