1 Délimitation des bassins versants

Ce document présente les étapes permettant d’effectuer la délimitation d’un bassin versant dans R, à partir d’un Modèle Numérique d’Élévation (MNE) (ou en anglais Digital Elevation Model, DEM) et d’un point d’écoulement précis.

Nous utiliserons des librairies spécialisées pour l’analyse géospatiale et hydrologique, notamment Terra et WhiteboxTools.

L’analyse inclut le prétraitement du MNE, le calcul des directions et de l’accumulation d’écoulement, le traçage du réseau hydrographique et la délimitation finale du bassin versant à partir d’un point d’exutoire prédéterminé.

1.1 Données nécessaires

Pour exécuter ce tutoriel, vous aurez besoin de :

Un Modèle Numérique d’Élévation de la région d’intérêt en format raster (p. ex. .tif)
Les coordonnées GPS du point d’exutoire à partir duquel vous souhaitez délimiter votre bassin versant

Préparation de vos fichiers

Nous vous conseillons de bien structurer vos fichiers dès le départ en suivant les étapes suivantes :

Créer un projet R. Celui-ci vous servira de répertoire de travail (working directory) tout au long de votre analyse.
Y ajouter un dossier Inputs/ dans lequel vous insérerez vos fichiers sources : MNE, point GPS d’exutoire, données hydrographiques, etc.
Y ajouter un dossier Outputs/ dans lequel vous sauvegarderez tous les fichiers produits lors de votre analyse.

1.2 Chargement des librairies

library(tidyverse)
library(tidyterra)
library(terra)      # Manipulation de données raster
library(sf)         # Manipulation d'objets vectoriels
library(mapview)    # Cartographie interactive

#whitebox::install_whitebox() # Installer une fois
library(whitebox)   # Interface avec WhiteboxTools
whitebox::wbt_init()

WhiteboxTools — qu’est-ce que c’est ?

WhiteboxTools est une librairie d’analyse géospatiale développée par le Dr. John Lindsay (Université de Guelph). Elle offre plus de 500 outils spécialisés en analyse hydrologique, géomorphométrie et traitement d’images raster. Dans ce tutoriel, nous l’utilisons principalement pour le calcul des directions d’écoulement, de l’accumulation de flux et de la délimitation du bassin versant.

wbt_version()

2 Importation et préparation du MNE

Cette section importe le MNE et effectue les premiers traitements nécessaires pour éliminer les valeurs aberrantes, puis corrige les dépressions topographiques.

Au Québec, les données LiDAR peuvent être téléchargées depuis le portail Forêt Ouverte, sous l’onglet Catalogue → LiDAR → Modèle de terrain.

Pourquoi remplir les dépressions du MNE ?

Les dépressions (ou « puits ») dans un MNE sont des cellules dont l’élévation est inférieure à toutes leurs voisines. Elles peuvent être des artefacts de l’interpolation ou représenter de véritables cuvettes. Dans les deux cas, elles interrompent le cheminement de l’eau simulée, car celle-ci s’y accumule sans pouvoir en sortir. La fonction wbt_fill_depressions() comble ces puits afin de garantir un écoulement continu à travers le MNE.

# Lecture du MNE (format .tif)
MNE = rast("Inputs/MNE_SBL_clip.tif")

# Nettoyage des valeurs négatives (souvent des artefacts d'interpolation)
MNE[MNE < 0] = NA

# Sauvegarde du MNE nettoyé (étape intermédiaire requise par WhiteboxTools)
terra::writeRaster(MNE, "Inputs/MNE_cleaned.tif", overwrite = TRUE)

# Remplissage des dépressions
wbt_fill_depressions(
  dem    = "Inputs/MNE_cleaned.tif",
  output = "Outputs/MNE_filled.tif"
)

# Rechargement du MNE rempli dans l'environnement R
MNE_filled = rast("Outputs/MNE_filled.tif")

# Visualisation du MNE rempli
ggplot(as.data.frame(MNE_filled, xy = TRUE)) +
  geom_raster(aes(x = x, y = y, fill = MNE_filled)) +
  coord_equal() +
  theme_minimal()

# Résolution spatiale du MNE
cat("Résolution spatiale du MNE (en mètres) :"); print(res(MNE_filled))

Résolution spatiale du MNE (en mètres) :

[1] 1 1

3 Calcul de la direction d’écoulement

L’algorithme D8

La méthode D8 (Deterministic 8) détermine la direction d’écoulement de chaque cellule vers l’une de ses 8 cellules voisines (horizontale, verticale ou diagonale) selon la pente maximale. Chaque direction est encodée par une valeur numérique (1, 2, 4, 8, 16, 32, 64 ou 128). Pour en savoir plus, consultez la documentation ArcGIS sur la direction d’écoulement.

wbt_d8_pointer(
  dem    = "Outputs/MNE_filled.tif",
  output = "Outputs/flowdir.tif"
)

flowdir = rast("Outputs/flowdir.tif")

# Le raster de direction ne devrait contenir que les valeurs : 1, 2, 4, 8, 16, 32, 64, 128
knitr::kable(freq(flowdir),
             caption = "Nombre de cellules par code de direction (algorithme D8)")

Nombre de cellules par code de direction (algorithme D8)
layer	value	count
1	0	718
1	1	886845
1	2	1354633
1	4	1201186
1	8	1507142
1	16	786024
1	32	1287452
1	64	1277246
1	128	1687134

4 Calcul de l’accumulation d’écoulement

L’accumulation de flux représente le nombre de cellules qui drainent vers chaque cellule du raster. Les cellules à forte accumulation correspondent aux cours d’eau principaux.

wbt_d8_flow_accumulation(
  input  = "Outputs/flowdir.tif",
  output = "Outputs/flow_accD8.tif",
  pntr   = TRUE   # Indique que l'entrée est un raster de direction (et non un MNE)
)

flow_acc = rast("Outputs/flow_accD8.tif")

knitr::kable(summary(values(flow_acc), na.rm = TRUE),
             caption = "Distribution des valeurs d'accumulation de flux")

Distribution des valeurs d’accumulation de flux
	flow_accD8
	Min. : 1
	1st Qu.: 1
	Median : 3
	Mean : 1497
	3rd Qu.: 11
	Max. :2382036

5 Localisation de l’exutoire

L’exutoire est le point par lequel l’eau quitte le bassin versant. Nous le définissons ici manuellement à partir de ses coordonnées GPS, puis nous ajustons sa position pour qu’il tombe exactement sur une cellule à forte accumulation de flux.

# Coordonnées de l'exutoire (système WGS84, lon/lat)
ppoints = tribble(
  ~Lon,      ~Lat,
  -74.01528, 45.98928
)

# Conversion en objet sf
ppointsSF = st_as_sf(ppoints, coords = c("Lon", "Lat"), crs = 4326)

# Reprojection dans le même SCR que le MNE
ppointsSF_proj = st_transform(ppointsSF, crs(MNE_filled))

# Export en shapefile (requis par WhiteboxTools)
st_write(ppointsSF_proj, "Inputs/pourpoints.shp", delete_layer = TRUE, quiet = TRUE)

Pourquoi ajuster la position du point d’exutoire ?

En pratique, les coordonnées GPS peuvent ne pas tomber exactement sur la cellule de fort débit dans le raster d’accumulation. La fonction wbt_snap_pour_points() déplace le point vers la cellule à la plus forte accumulation dans un rayon donné, garantissant ainsi que la délimitation se fait à partir d’un point hydrologique cohérent.

wbt_snap_pour_points(
  pour_pts   = "Inputs/pourpoints.shp",
  flow_accum = "Outputs/flow_accD8.tif",
  output     = "Outputs/pourpoints_snapped.shp",
  snap_dist  = 5   # Distance maximale de déplacement, en mètres
)

snapped_pt = st_read("Outputs/pourpoints_snapped.shp", quiet = TRUE)

6 Délimitation du bassin versant

À partir du point d’exutoire ajusté et du raster de direction d’écoulement, wbt_watershed() remonte le réseau hydrologique pour identifier toutes les cellules qui drainent vers ce point.

wbt_watershed(
  d8_pntr  = "Outputs/flowdir.tif",
  pour_pts = "Outputs/pourpoints_snapped.shp",
  output   = "Outputs/watershed.tif"
)

watershed_raster = rast("Outputs/watershed.tif")

# Conversion du raster en polygone vectoriel
watershed_vector = as.polygons(watershed_raster) |> st_as_sf()
st_crs(watershed_vector) = st_crs(MNE_filled)

# Sauvegarde
st_write(watershed_vector, "Outputs/watershed.shp", delete_layer = TRUE, quiet = TRUE)

# Visualisation du bassin versant et point d'écoulement
ggplot() +
  geom_spatraster(data = MNE_filled) +
  geom_sf(data = watershed_vector, fill = "blue", alpha = 0.4, color = NA) +
  geom_sf(data = snapped_pt, color = "red", size = 3) +
  scale_fill_gradientn(colors = terrain.colors(256))+
  coord_sf() +
  theme_minimal()

#mapview(MNE_filled) +
#  mapview(watershed_vector, alpha.regions = 0.2) +
#  mapview(snapped_pt)

# Calcul de la superficie
watershed_vector$area_km2 = as.numeric(st_area(watershed_vector)) / 1e6
cat("Superficie du bassin versant (km²) :"); print(round(watershed_vector$area_km2, 2))

Superficie du bassin versant (km²) :

[1] 0.22

7 Extraction du réseau hydrographique

Pour extraire les cours d’eau, nous appliquons un seuil sur le raster d’accumulation de flux : toute cellule dont la valeur d’accumulation dépasse ce seuil est considérée comme appartenant au réseau hydrographique.

Choisir le bon seuil

Le seuil (threshold) détermine la densité du réseau extrait : une valeur élevée produit un réseau peu dense (seulement les cours d’eau principaux), tandis qu’une valeur faible produit un réseau très dense (incluant les petits ruisseaux, souvent intermittents). Il n’existe pas de valeur universelle — elle dépend de la résolution du MNE, des caractéristiques de la zone d’étude et des conditions hydrologiques que l’on tente de représenter.

Comparez le résultat avec des données hydrographiques de référence pour valider votre choix de seuil. Ces données peuvent provenir d’observations terrain ou de données cartographiques du Québec tels que:

Les lits d’écoulement potentiel issue du LIDAR (cours d’eau de tête)
Cours d’eau intermittents (Rivières, ruisseaux, etc.)
Cours d’eau permanents (Fleuve, lacs, rivières, ruisseaux, etc.)
Plans d’eau (Fleuve, lacs, rivières, réservoirs, etc.)

tous disponible depuis le portail Forêt Ouverte, sous l’onglet Catalogue → Hydrographie.

À l’échelle globale, les données HydroSHEDS (WWF / Université McGill) constituent également une référence utile.

wbt_raster_streams_to_vector vs wbt_raster_to_vector_lines

Il existe deux façons de convertir un raster de cours d’eau en vecteur :

wbt_raster_to_vector_lines() effectue une conversion directe cellule par cellule, sans tenir compte de la topologie du réseau.
wbt_raster_streams_to_vector() utilise le raster de direction d’écoulement (D8) pour reconstituer les segments dans le bon ordre amont–aval. Chaque segment relie une source à une confluence, ce qui produit un réseau topologiquement cohérent, mieux adapté aux analyses hydrologiques.

Nous utilisons ici la deuxième approche.

wbt_extract_streams(
  flow_accum = "Outputs/flow_accD8.tif",
  output     = "Outputs/raster_streams.tif",
  threshold  = 7000  # Seuil ajusté pour représenter le réseau au pic printanier
)

# Conversion vectorielle en suivant la direction d'écoulement (D8)
wbt_raster_streams_to_vector(
  streams  = "Outputs/raster_streams.tif",
  d8_pntr  = "Outputs/flowdir.tif",
  output   = "Outputs/stream_vector.shp"
)

stream_vector = st_read("Outputs/stream_vector.shp", quiet = TRUE)
st_crs(stream_vector) = st_crs(MNE_filled)  # Assigner le SCR du MNE

# Visualisation finale en mode interactif
# Visualisation du bassin versant et point d'écoulement

ggplot() +
  geom_spatraster(data = MNE_filled) +
  geom_sf(data = watershed_vector, fill = "blue", alpha = 0.4, color = NA) +
  geom_sf(data=stream_vector, color="turquoise", size=0.5)+
  geom_sf(data = snapped_pt, color = "red", size = 3) +
  scale_fill_gradientn(colors = terrain.colors(256))+
  coord_sf() +
  theme_minimal()

Visualisation interactive

Ce code permet de visualiser la carte en mode interactif 🤩

#mapview(MNE_filled) + 
# mapview(stream_vector, color = "turquoise") + 
# mapview(watershed_vector, alpha.regions = 0.2) + 
# mapview(snapped_pt)