Délinéation de bassins versants

Author

Joël Gaxotte

Ce document présente un workflow protoype et incomplet pour la délimitation de bassins versants à partir d’un Modèle Numérique d’Élévation (MNE) en utilisant les packages R whitebox, terra, raster et sf. Le processus inclut le prétraitement du MNE, trouver un point d’écoulement, la génération de réseaux hydrographiques et la délimitation finale des bassins versants.

Chargement des packages

#Packages de base pour la manipulation de données spatiales
library(tidyverse)
library(raster)
library(terra)      # Alternative moderne au package raster
library(sf)         # Pour les objets vectoriels
library(tmap)       # Cartographie interactive
library(stars)      # Autre format de données spatiales
library(whitebox)   # Interface avec WhiteboxTools

#Initialisation de WhiteboxTools
whitebox::wbt_init()
wbt_version()

#Configuration des thèmes et modes d'affichage
theme_set(theme_classic())
tmap_mode("view")

Importation et préparation du MNE

Cette section importe le Modèle Numérique d’Élévation (MNE) et effectue les premiers traitements pour éliminer les valeurs aberrantes.

# Lecture du fichier MNE (format .adf)
MNE <- rast("/Users/thejo/Desktop/TD/Donnees_Audrey/DEM/demwatershed")

# Nettoyage des valeurs négatives (souvent des artefacts)
MNE[MNE < 0] <- NA

# Visualisation interactive du MNE
tm_shape(MNE) +
  tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE) +
  tm_scale_bar()

2. Conversion et ombrage du MNE

Conversion au format .tif pour une meilleure compatibilité et génération d’un ombrage pour la visualisation.

# Conversion en format .tif
writeRaster(MNE, "/Users/thejo/Desktop/TD/trois/donnees/MNE.tif", overwrite = TRUE)

# Génération d'un ombrage (hillshade)
wbt_hillshade(
  dem = "/Users/thejo/Desktop/TD/trois/donnees/MNE.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/hillshademne.tif",
  azimuth = 115  # Angle d'éclairage en degrés
)

# Visualisation de l'ombrage
hillshade2 <- raster("/Users/thejo/Desktop/TD/trois/donnees/hillshademne.tif")
tm_shape(hillshade2) +
  tm_raster(style = "cont", palette = "-Greys", legend.show = FALSE) +
  tm_scale_bar()

Correction des dépressions

Les dépressions dans le MNE peuvent affecter l’analyse hydrologique. Cette section les corrige.

# Conversion en format float32 avec valeur NA spécifique
terra::writeRaster(
  MNE, 
  "/Users/thejo/Desktop/TD/trois/donnees/MNE_float32.tif", 
  datatype = "FLT4S", 
  overwrite = TRUE, 
  NAflag = -9999
)

# Correction des dépressions par la méthode du moindre coût
wbt_breach_depressions_least_cost(
  dem = "/Users/thejo/Desktop/TD/trois/donnees/MNE_float32.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/MNE_breach_dep.tif",
  dist = 3,  # Distance maximale de correction
  fill = TRUE
)

# Remplissage complémentaire des dépressions
wbt_fill_depressions_wang_and_liu(
  dem = "/Users/thejo/Desktop/TD/trois/donnees/MNE_breach_dep.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif"
)
MNE_clean <- rast("/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif")

# Vérification des résultats
summary(values(MNE_clean))

Analyse du réseau hydrographique

Calcul de l’accumulation de flux et de la direction d’écoulement.

# Accumulation de flux (méthode D8)
wbt_d8_flow_accumulation(
  input = "/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/flow_accD8.tif"
)

  # Transformation logarithmique pour meilleure visualisation sur ArcGIS pro
flow_acc <- rast("/Users/thejo/Desktop/TD/trois/donnees/flow_accD8.tif")
scaled_acc <- log1p(flow_acc)
terra::writeRaster(
  scaled_acc, 
  "/Users/thejo/Desktop/TD/trois/donnees/flow_accD8_scaled.tif", 
  overwrite = TRUE
)

# Direction d'écoulement (méthode D8)
wbt_d8_pointer(
  dem = "/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/pointer3.tif"
)

# Check flow direction values (should only contain powers of 2: 1,2,4,8,16,32,64,128)
freq(pointer3)

Notes

Ainsi, si vous voyez 0 dans freq(pointer3), cela signifie que WhiteboxTools n’a pas pu assigner une direction d’écoulement à ces pixels, et non pas que leur valeur d’élévation était de 0. Cela se produit généralement dans les cas suivants :

Le pixel se trouve dans un puits ou une zone plate qui n’a pas été complètement résolue.
Le pixel est entouré de NA (bord de la trame)
Le pixel est entouré de NA (bord de la trame).
Il reste un petit plat ou une dépression après le prétraitement.

Il faut faire des corrections plus robustes et aggresives, mais dépendemment du MNA et des outils disponible, ce n’est pas toujours possible.

Ceci vient donc bloquer le reste de la délimitation. Les étapes concernant les points d’écoulement peuvent tout de même être accomplies sans le Flow Direction (pointer).

Définition des points d’exutoire

Avec point(s) déjà défini(s)

Création et préparation des points d’exutoire pour la délimitation des bassins versants.

# Création du ou des points d'exutoires
ppoints <- tribble(
  ~Lon, ~Lat,
  -63.369243, 50.828091
)

# Conversion en objet spatial
ppointsSP <- SpatialPoints(ppoints, proj4string = CRS("+proj=longlat +datum=WGS84"))

# Option 1: Export en shapefile
shapefile(ppointsSP, filename = "pourpoints.shp", overwrite = TRUE)

Sans point défini (non complété)

Cette étape n’est également pas fini.

#Étapes à faire si pas de points d'écoulement définis

#Trouvé et appliquer le percentile pour isoler les zones de fortes accumulations (Choix à faire selon le cas)
#95th percentile → petit cours d'eau
#98th percentile → medium
#99th percentile → principaux cours d'eau

quantile(values(flow_acc), probs = c(0.95, 0.98, 0.99), na.rm = TRUE)

#95%  98%  99% 
#126 1052 5253

#99%
high_flow <- flow_acc >= 5253
terra::writeRaster(high_flow, 
                   "/Users/thejo/Desktop/TD/trois/donnees/high_flow.tif", 
                   overwrite = TRUE)

#98%
medium_flow <- flow_acc >= 1052
terra::writeRaster(medium_flow, 
                   "/Users/thejo/Desktop/TD/trois/donnees/medium_flow.tif", 
                   overwrite = TRUE)
#95%
low_flow <- flow_acc >= 126
terra::writeRaster(low_flow, 
                   "/Users/thejo/Desktop/TD/trois/donnees/low_flow.tif", 
                   overwrite = TRUE)

#Besoin : flow_acc, pointer(flow_direct), lacs et rivières shp, medium ou high flow

lacs_riv <- st_read("/Users/thejo/Desktop/TD/Donnees_Audrey/Hydrologie/water_bodies_all.shp")

#Identifier les cellules en bordure du raster

# Créer un masque de bordure
border_mask <- focal(flow_acc, w = matrix(1, 3, 3), fun = function(x) any(is.na(x)), na.policy = "omit")

# Convertir les valeurs TRUE/FALSE en 0/1
border_mask[!is.na(border_mask)] <- 1
border_mask[is.na(border_mask)] <- 0

#---Automatique

# Trouver la cellule avec l'accumulation maximale
lmax_acc <- which.max(values(low_flow))
coords_outlet <- xyFromCell(low_flow, lmax_acc)

# Convertir en objet sf
outlet_point <- st_as_sf(data.frame(coords_outlet), coords = c("x", "y"), crs = crs(flow_acc))

#Visualiser le point
library(tmap)
tmap_mode("view")  # Interactive mode (like Leaflet)

# Plot the flow accumulation raster + outlet point
tm_shape(flow_acc) + 
  tm_raster(style = "cont", palette = "Blues", alpha = 0.7, title = "Flow Accumulation") +
  tm_shape(outlet_point) + 
  tm_dots(col = "red", size = 0.5, title = "Outlet Point") +
  tm_layout(legend.position = c("left", "bottom")) +
  tm_scale_bar(position = c("right", "bottom"))

#Méthode manuelle

wbt_gaussian_filter(
  input = "/Users/thejo/Desktop/TD/trois/donnees/medium_flow.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/medium_flow_smoothed.tif",
  sigma = 2  # Ajuster pour plus/moins de lissage
)

medium_flow_smoothed <- rast("/Users/thejo/Desktop/TD/trois/donnees/medium_flow_smoothed.tif")

# Convert medium_flow to factor (forces binary treatment)
medium_flow_smoothed_fact <- as.factor(medium_flow_smoothed)

tm_shape(low_flow) +
  tm_raster(col.scale = tm_scale_categorical(values = c("white", "black"),
                                             labels = c("No Flow", "Flow Lines"))) +
  tm_raster(col.scale = tm_scale_continuous(values = "viridis"),
            col_alpha = 0.5)

Extraction du réseau hydrographique

Le(s) point(s) d’écoulement sont également « snapped » ensemble.

wbt_extract_streams(
  flow_accum = "/Users/thejo/Desktop/TD/trois/donnees/flow_accD8.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/raster_streams.tif",
  threshold = 1052  # Ajuster selon la résolution et la zone d'étude
)

wbt_jenson_snap_pour_points(
  pour_pts = "/Users/thejo/Desktop/TD/trois/donnees/pour_points.shp",
  streams = "/Users/thejo/Desktop/TD/trois/donnees/raster_streams.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/snappedpp.shp",
  snap_dist = 100  # Distance d'ajustement (unités du SCR)
)
# Vérifier concordance des systèmes de coordonnées
print("Pointer3 CRS:")
print(crs(pointer3))

print("Snappp CRS:")
print(st_crs(snappp))

# Visualisation
snappp <- shapefile("/Users/thejo/Desktop/TD/trois/donnees/snappedpp.shp")
streams <- raster("/Users/thejo/Desktop/TD/trois/donnees/raster_streams.tif")

tm_shape(streams) +
  tm_raster(legend.show = TRUE, palette = "Blues") +
  tm_shape(snappp) +
  tm_dots(col = "red")

Délimitation des bassins versants

Calcul final des bassins versants pour chaque point d’exutoire.

wbt_watershed(d8_pntr = "/Users/thejo/Desktop/TD/trois/donnees/pointer3.tif",
              pour_pts = "/Users/thejo/Desktop/TD/trois/donnees/snappedpp.shp",
              output = "/Users/thejo/Desktop/TD/trois/donnees/brush_watersheds.tif")

ws <- raster("/Users/thejo/Desktop/TD/deuxieme/brush_watersheds.tif")

tm_shape(hillshade2)+
  tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
  tm_shape(ws)+
  tm_raster(legend.show = TRUE, alpha = 0.5, style = "cat")+
  tm_shape(pp)+
  tm_dots(col = "red")