Localisation des zones d’exploitation optimale du Pin blanc avec R

Tutoriel pédagogique décrivant une méthode SIG sous R pour localiser les zones d’exploitation optimale du Pin blanc à partir d’un raster d’occupation du sol, en intégrant des contraintes de distance hydrologique et de superficie minimale.

Auteur·rice

NEBIE G.C.

Date de publication

26 novembre 2025

Résumé

Ce tutoriel présente une démarche complète, étape par étape, permettant d’identifier les zones optimales d’exploitation du Pin blanc (Pin de Weymouth) à partir d’un raster d’occupation du sol, en utilisant exclusivement le langage de programmation R.

L’analyse repose sur l’application successive de trois critères spatiaux : la nature de l’essence forestière (Pin blanc), la distance minimale de 700 mètres par rapport aux plans d’eau, et une superficie minimale de 1000 hectares par parcelle. La méthodologie combine des traitements raster et vectoriels, incluant la conversion des données, la création de zones tampons (buffer), les opérations de découpage spatial, le calcul des superficies et la cartographie thématique finale.

Le tutoriel met en évidence la capacité de R à fonctionner comme un véritable Système d’Information Géographique (SIG), en proposant une approche reproductible, transparente et accessible pour la planification forestière et l’aide à la décision spatiale.

Introduction

Dans ce tutoriel, nous allons montrer étape par étape comment identifier, à l’aide de R, les zones optimales d’exploitation du Pin blanc (Pin de Weymouth) à partir d’un raster d’occupation du sol.

Nous appliquerons successivement trois critères :

Essence forestière : seules les zones de Pin blanc sont retenues.
Distance aux plans d’eau : aucune parcelle ne doit se situer à moins de 700 m de l’eau.
Superficie minimale : seules les parcelles d’au moins 1000 hectares sont considérées exploitables.

À la fin du tutoriel, nous disposerons :

d’un jeu de polygones représentant les blocs exploitables,
d’une carte finale combinant l’occupation du sol et les parcelles éligibles.

Objectif du tutoriel

Ce document guide pas à pas l’utilisateur dans la mise en oeuvre d’une analyse spatiale sous R visant à identifier les zones optimales d’exploitation du Pin blanc à partir de données raster. Il s’adresse aussi bien aux étudiants débutants en SIG qu’aux praticiens souhaitant automatiser leurs traitements spatiaux.

Public cible

Étudiants en géographie, environnement, aménagement du territoire
Gestionnaires forestiers
Utilisateurs débutants en SIG sous R
Analystes souhaitant reproduire un workflow SIG sans logiciel propriétaire

Préparation de l’environnement R

Nous commençons par charger les packages nécessaires pour :

manipuler les rasters (terra),
travailler avec les objets vectoriels (sf),
faire des graphiques (ggplot2, ggspatial),
utiliser la syntaxe tidyverse.

Afficher le code

library(terra)
library(sf)
library(ggspatial)
library(tidyterra)
library(tidyverse)

Import du raster et définition du système de coordonnées

Nous chargeons le raster d’occupation du sol landuse1.tif. Ce raster contient les classes suivantes :

1 : Pin blanc
2 : Pin gris
3 : Plans d’eau

Nous procédons ensuite au chargement du fichier raster représentant l’occupation du sol, nommé landuse, qui constitue la base de toutes les analyses spatiales ultérieures.

Afficher le code

landuse <- rast("landuse/landuse1.tif")

landuse

class       : SpatRaster 
dimensions  : 77, 205, 1  (nrow, ncol, nlyr)
resolution  : 350, 350  (x, y)
extent      : 0, 71750, 350, 27300  (xmin, xmax, ymin, ymax)
coord. ref. :  
source      : landuse1.tif 
name        : landuse1

Vérifions le système de coordonnées (CRS) :

Afficher le code

crs(landuse)

[1] ""

Dans notre cas, aucun CRS n’était défini dans le fichier d’origine. Nous le définissons donc selon les consignes du travail pratique : WGS84 / UTM Zone 17N (EPSG:32617).

Afficher le code

crs(landuse) <- "EPSG:32617"

crs(landuse)

[1] "PROJCRS[\"WGS 84 / UTM zone 17N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 17N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-81,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n        AREA[\"Between 84°W and 78°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Ecuador - north of equator. Canada - Nunavut; Ontario; Quebec. Cayman Islands. Colombia. Costa Rica. Cuba. Jamaica. Nicaragua. Panama. United States (USA).\"],\n        BBOX[0,-84,84,-78]],\n    ID[\"EPSG\",32617]]"

Nous créons ensuite un objet de travail landuse_proj :

Afficher le code

landuse_proj <- landuse

plot(landuse_proj, main = "Occupation du sol - landuse")

Comprendre les classes d’occupation du sol

Pour rappel :

1 = Pin blanc
2 = Pin gris
3 = Eau

Nous allons maintenant extraire deux couches raster :

une pour le Pin blanc,
une pour les plans d’eau.

Afficher le code

# 1 = Pin blanc, 2 = Pin gris, 3 = Eau

pin_blanc <- landuse_proj == 1
pin_blanc[pin_blanc == 0] <- NA

eau <- landuse_proj == 3
eau[eau == 0] <- NA

par(mfrow = c(1, 2))
plot(pin_blanc, main = "Pixels de Pin blanc")
plot(eau, main = "Pixels d'eau")

Afficher le code

par(mfrow = c(1, 1))

Pourquoi une approche vectorielle ?

Bien que les données initiales soient raster, le passage en format vectoriel permet un découpage plus précis des parcelles intersectant partiellement la zone tampon. Contrairement au masquage raster, qui supprime des cellules entières, l’approche vectorielle autorise une géométrie plus fidèle des limites exploitables.

Création d’un buffer de 700 m autour des plans d’eau

L’objectif est maintenant de créer une zone de protection autour des plans d’eau : toute exploitation est interdite dans un rayon de 700 m.

Conversion de la couche “eau” en polygones

Afficher le code

# Raster eau -> polygones (fusion des zones contiguës)

eau_vect_terra <- as.polygons(eau, dissolve = TRUE, na.rm = TRUE)

# Passage au format sf

eau_vect <- st_as_sf(eau_vect_terra)

# S'assurer que les géométries sont valides

eau_vect <- st_make_valid(eau_vect)

# Création du buffer de 700 m autour de l'eau

eau_buffer <- eau_vect |>
st_union() |>
st_buffer(
dist        = 700,
endCapStyle = "SQUARE",
joinStyle   = "MITRE"
)

plot(st_geometry(eau_buffer),
main = "Buffer 700 m autour des zones d'eau",
col  = "lightblue")

À ce stade, eau_buffer représente toute la zone située à moins de 700 m de l’eau.

Extraction des parcelles de Pin blanc candidates

Nous allons maintenant appliquer le buffer aux parcelles de Pin blanc.

Afficher le code

# Raster Pin blanc -> polygones (parcelles continues)

pin_vect_terra <- as.polygons(pin_blanc, dissolve = TRUE, na.rm = TRUE)

# Passage en sf

pin_vect <- st_as_sf(pin_vect_terra)

# Correction éventuelle des géométries

pin_vect    <- st_make_valid(pin_vect)
eau_buffer  <- st_make_valid(eau_buffer)

# Harmonisation des systèmes de coordonnées

st_crs(pin_vect) <- st_crs(eau_buffer)

# Découpage : retrait de la partie à l'intérieur du buffer

pin_hors_buffer <- st_difference(pin_vect, eau_buffer)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Afficher le code

# Suppression des géométries vides

pin_hors_buffer <- pin_hors_buffer[!st_is_empty(pin_hors_buffer), ]

# Conversion en POLYGON si nécessaire

pin_hors_buffer <- st_cast(pin_hors_buffer, "POLYGON")

Warning in st_cast.sf(pin_hors_buffer, "POLYGON"): repeating attributes for all
sub-geometries for which they may not be constant

Afficher le code

plot(st_geometry(pin_hors_buffer),
main = "Pin blanc hors zone buffer 700 m")

Les polygones obtenus représentent les parcelles de Pin blanc qui se trouvent à plus de 700 m des plans d’eau. Deux des trois critères sont donc déjà respectés.

Calcul de la surface et sélection des blocs exploitables

Il reste à appliquer le troisième critère : superficie minimale de 1000 hectares.

Afficher le code

# Surface en m² puis conversion en hectares

pin_hors_buffer$area_ha <- as.numeric(st_area(pin_hors_buffer)) / 10000

# Filtrage des blocs exploitables (≥ 1000 ha)

blocs_exploitables <- pin_hors_buffer |>
dplyr::filter(area_ha >= 1000)

# Attribution d'un identifiant unique de bail

if (nrow(blocs_exploitables) > 0) {
blocs_exploitables$id_bail <- paste0(
"B",
sprintf("%03d", seq_len(nrow(blocs_exploitables)))
)
}

blocs_exploitables

Simple feature collection with 5 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 15750 ymin: 350 xmax: 51800 ymax: 27300
Projected CRS: WGS 84 / UTM zone 17N
     landuse1                       geometry   area_ha id_bail
1.2         1 POLYGON ((38910.05 27300, 3... 10393.003    B001
1.3         1 POLYGON ((51100 27300, 5110...  3122.795    B002
1.17        1 POLYGON ((35000 2450, 34650...  1774.813    B003
1.19        1 POLYGON ((20650 2450, 20650...  1578.931    B004
1.26        1 POLYGON ((16100 22460.05, 1...  1474.927    B005

Chaque ligne de blocs_exploitables représente une parcelle candidate qui :

est du Pin blanc,
se situe à plus de 700 m de l’eau,
possède une superficie ≥ 1000 ha,
dispose d’un identifiant unique id_bail.

Cartographie finale simple

Nous produisons d’abord une carte simple avec le raster en fond et les blocs exploitables en rouge.

Afficher le code

# Palette pour l'occupation du sol

cols_landuse <- c(
"1" = "darkgreen",   # Pin blanc
"2" = "palegreen3",  # Pin gris
"3" = "lightblue"    # Plan d'eau
)

plot(landuse_proj, col = cols_landuse, main = "Blocs exploitables (vecteur)")
if (nrow(blocs_exploitables) > 0) {
plot(st_geometry(blocs_exploitables),
add   = TRUE,
border= "red",
lwd   = 2)
}

Cartographie finale stylisée avec ggplot2

Pour une carte plus présentable (par exemple pour un rapport ou une présentation), nous pouvons utiliser ggplot2 et ggspatial.

Afficher le code

# Conversion raster -> data.frame pour ggplot

landuse_df <- as.data.frame(landuse_proj, xy = TRUE) %>%
rename(valeur = 3) %>%  # Adapter si nécessaire au nom de la colonne
mutate(
categorie = case_when(
valeur == 1 ~ "Pin blanc",
valeur == 2 ~ "Pin gris",
valeur == 3 ~ "Plan d'eau"
)
)

plot_parcelles <- ggplot() +

# Fond raster

geom_tile(
data = landuse_df,
aes(x = x, y = y, fill = categorie)
) +
scale_fill_manual(
name   = "Occupation du sol",
values = c(
"Pin blanc" = "#2E8B57",
"Pin gris"  = "#90EE90",
"Plan d'eau"= "#4682B4"
),
na.value = "grey90"
) +

# Parcelles exploitables

geom_sf(
data      = blocs_exploitables,
aes(color = id_bail),
fill      = "red",
linewidth = 0.9,
alpha     = 0.5
) +
scale_color_viridis_d(
name  = "Parcelles",
option = "plasma",
guide = guide_legend(
keywidth  = unit(0.5, "cm"),
keyheight = unit(0.5, "cm"),
override.aes = list(linewidth = 1.5)
)
) +

# Flèche nord et échelle graphique

annotation_north_arrow(
location = "tl",
style    = north_arrow_minimal(
text_size = 10,
line_width = 1
),
pad_x = unit(0.3, "cm"),
pad_y = unit(0.3, "cm")
) +
annotation_scale(
location   = "br",
width_hint = 0.25,
bar_cols   = c("black", "white"),
text_cex   = 0.9,
height     = unit(0.15, "cm")
) +
labs(
title = "Occupation du sol et parcelles exploitables",
x = NULL,
y = NULL
) +
theme_light() +
theme(
plot.title = element_text(
hjust  = 0.5,
face   = "bold",
size   = 14,
margin = margin(b = 10)
),
legend.position   = "right",
legend.box        = "vertical",
legend.spacing.y  = unit(0.3, "cm"),
legend.text       = element_text(size = 9),
legend.title      = element_text(size = 10, face = "bold"),
axis.text         = element_text(size = 8),
panel.grid        = element_line(color = "grey85", linewidth = 0.3),
panel.background  = element_rect(fill = "white"),
plot.margin       = margin(10, 10, 10, 10)
) +
guides(
fill  = guide_legend(order = 1),
color = guide_legend(order = 2)
)

print(plot_parcelles)

La carte obtenue montre clairement que toutes les parcelles de Pin blanc ne sont pas automatiquement exploitables. Plusieurs sont éliminées en raison de leur proximité avec les plans d’eau ou de leur superficie insuffisante. Seules les parcelles éloignées des zones hydrologiques sensibles et présentant une superficie suffisante sont retenues pour l’analyse finale.

Résumé de la démarche

En résumé, la démarche sous R a consisté à :

Importer le raster d’occupation du sol et définir le CRS.
Extraire les classes Pin blanc et Eau.
Convertir la couche d’eau en polygones et créer un buffer de 700 m.
Convertir le Pin blanc en polygones et retirer la partie à l’intérieur du buffer (distance > 700 m).
Calculer la surface des parcelles restantes et sélectionner celles de ≥ 1000 ha.
Attribuer un identifiant de bail à chaque bloc exploitable.
Produire des cartes finales (simple et stylisée) montrant les parcelles éligibles.

Ce flux de travail illustre comment R peut être utilisé comme un SIG complet pour réaliser des analyses spatiales complexes de manière transparente, reproductible et documentée.