8.2 Exercices

Dans cette section, vous mettrez en pratique certains concepts vus dans la section leçon de ce module. Bien que la réponse à chaque question soit disponible, il est très important de tenter d’y répondre par vous-même!

Question 1

Identifier le point d’élévation maximale sur une carte du parc du Mont-Orford.

Réponse

Isoler le polygone du parc du Mont-Orford du shapefile parcs en utilisant la fonction subset:

parc_orford <- subset(parcs, TRQ_NM_== "Parc national du Mont-Orford")

Filtrer le raster dem_lcc pour ne retenir que les cellules situées à l’intérieure du polygone délimitant les limites du parc du Mont-Orford

dem_orford <- mask(dem_lcc, parc_orford)
mapview(dem_orford)

Déterminer l’indice de la cellule ou des cellules d’élévation maximale en utilisant les fonctions getValues() et which.max().

imax_orford<-which.max(getValues(dem_orford))

Trouver les coordonnées correspondant à cet indice en utilisant la fonction xyFromCell.

coordmax_orford <- xyFromCell(dem_orford,imax_orford)
coordmax_orford

           x      y
[1,] -293815 154024

Transformer ces coordonnées en une donnée vectorielle de type POINT en utilisant la fonction st_as_sf().

#Créer d'abord un data.frame à partir de coordmax_orford
coordmax_orford_df <- as.data.frame(coordmax_orford) #ou simplement data.frame(coordmax_orford)

#Créer un point
pt_max_orford <- st_as_sf(coordmax_orford_df,
                   coords = c("x", "y"),
                   crs = st_crs(dem_orford)) #attribuer le même SCR

Visualiser la carte du parc du Mont-Orford ainsi que le point d’élévalion maximale.

mapview(dem_orford) + mapview(pt_max_orford, col.regions="red")

Question 2

Reclassifier le raster dem_orford selon quatre classes correspondant aux valeurs comprises entre le zéro et le \(1^{er}\) quantile, le \(1^{er}\) et le \(2^{ième}\), le \(2^{ième}\) et \(3^{ième}\), et le \(3^{ième}\) et \(4^{ième}\).

Réponse

Trouver d’abord les quantiles des valeurs d’élévation dans le parc du Mont-Orford.

quantile_orford <- quantile(getValues(dem_orford), na.rm = TRUE)
quantile_orford

   0%   25%   50%   75%  100% 
270.1 338.6 394.2 536.6 816.6 

q_25 <- as.numeric(quantile_orford["25%"]) # as.numeric pour conserver seulement les chiffres
q_50 <- as.numeric(quantile_orford["50%"])
q_75 <- as.numeric(quantile_orford["75%"])
q_100 <- as.numeric(quantile_orford["100%"])

En se servant de ces valeurs, construire une matrice qui détermine la nouvelle classification.

classes_orford <- matrix(c(0, q_25, 1,       #1ere classe
                           q_25, q_50, 2,    #2e classe
                           q_50, q_75, 3,    #3e classe
                           q_75, q_100, 4),  #4e classe
                           nrow = 4, ncol = 3, byrow = TRUE)
# donner des titres aux colonnes de cette matrice
colnames(classes_orford) <- c("Limite_min", "Limite_max", "Classement_quantile")

Utiliser la fonction reclassify() pour assigner les nouvelles classes au raster dem_orford.

dem_orford_quantile <- reclassify(dem_orford, classes_orford, rigth = FALSE)

Confirmer visuellement la reclassification en utilisant la fonction mapview()

mapview(dem_orford_quantile)

Noter que nous aurions pu vouloir attribuer une nouvelle classification de type catégorique au raster dem_orford. Par exemple,

Categorie_elevation <- c("Faible", "Intermediare", "Forte", "Très forte")

# Ou encore
Categorie_quantile <- c("[0%, 25%[", "[25%, 50%[", "[50%, 75%[", "[75%, 100%]")

Dans ce cas, il faut ajouter ces catégories au raster dem_orford en utilisant la fonction levels().

levels(dem_orford_quantile) <- data.frame(ID=1:4, quantile = Categorie_quantile)

dem_orford_quantile

class      : RasterLayer 
dimensions : 178, 304, 54112  (nrow, ncol, ncell)
resolution : 521, 741  (x, y)
extent     : -336797, -178413, 110676, 242574  (xmin, xmax, ymin, ymax)
crs        : +proj=lcc +lat_0=44 +lon_0=-68.5 +lat_1=60 +lat_2=46 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 
source     : memory
names      : DEM 
values     : 1, 4  (min, max)
attributes :
 ID    quantile
  1   [0%, 25%[
  2  [25%, 50%[
  3  [50%, 75%[
  4 [75%, 100%]

La classification est maintenant traitée comme une variable catégorique.

mapview(dem_orford_quantile)

Question 3

À partir du raster dem_lcc, créer un seul raster composé de zones tampons circulaires, de 10 km de rayon, autour des points d’élévation maximale de chaque parc national de la région de Sherbrooke.

Réponse

Répéter les étapes de la Question 1 pour créer un shapefile correspondant au point d’élévation maximale de chaque parc.

Pour le Parc du Mont-Orford:

# Isoler les limites du parc
parc_orford <- subset(parcs, TRQ_NM_== "Parc national du Mont-Orford")
# Trouver dem à l'intérieur du parc
dem_orford <- mask(dem_lcc, parc_orford)
# Trouver l'indice de la cellule d'élévation maximale
imax_orford<-which.max(getValues(dem_orford))
# Trouver les coordonnées de cette cellule
coordmax_orford <- xyFromCell(dem_orford,imax_orford)
# Transformer en data.frame
coordmax_orford_df <- as.data.frame(coordmax_orford)
# Transformer en point
pt_max_orford <- st_as_sf(coordmax_orford_df,
                   coords = c("x", "y"),
                   crs = st_crs(dem_orford)) #attribuer le même SCR

Pour le Parc national du Mont-Mégantic:

parc_megantic <- subset(parcs, TRQ_NM_== "Parc national du Mont-Mégantic")
dem_megantic <- mask(dem_lcc, parc_megantic)
imax_megantic <- which.max(getValues(dem_megantic))
coordmax_megantic <- xyFromCell(dem_megantic,imax_megantic)
coordmax_megantic_df <- as.data.frame(coordmax_megantic)
pt_max_megantic <- st_as_sf(coordmax_megantic_df,
                   coords = c("x", "y"),
                   crs = st_crs(dem_orford)) 

Pour le Parc national de la Yamaska:

parc_yamaska <- subset(parcs, TRQ_NM_== "Parc national de la Yamaska")
dem_yamaska <- mask(dem_lcc, parc_yamaska)
imax_yamaska <- which.max(getValues(dem_yamaska))
coordmax_yamaska <- xyFromCell(dem_yamaska,imax_yamaska)
coordmax_yamaska_df <- as.data.frame(coordmax_yamaska)
pt_max_yamaska <- st_as_sf(coordmax_yamaska_df,
                   coords = c("x", "y"),
                   crs = st_crs(dem_orford)) 

Pour le Parc national de Frontenac:

parc_frontenac <- subset(parcs, TRQ_NM_== "Parc national de Frontenac")
dem_frontenac <- mask(dem_lcc, parc_frontenac)
imax_frontenac <- which.max(getValues(dem_frontenac))
coordmax_frontenac <- xyFromCell(dem_frontenac,imax_frontenac)
coordmax_frontenac_df <- as.data.frame(coordmax_frontenac)
pt_max_frontenac <- st_as_sf(coordmax_frontenac_df,
                   coords = c("x", "y"),
                   crs = st_crs(dem_orford)) 

Confirmer vos calculs en visualisant les points trouvés.

mapview(pt_max_orford, col.regions = "red") + 
  mapview(pt_max_megantic, col.regions = "blue") +
  mapview(pt_max_yamaska, col.regions = "green") +
  mapview(pt_max_frontenac, col.regions = "yellow")

m <- mapview(pt_max_orford, col.regions = "red") + 
  mapview(pt_max_megantic, col.regions = "blue") +
  mapview(pt_max_yamaska, col.regions = "green") +
  mapview(pt_max_frontenac, col.regions = "yellow")
m@map

Calculer une zone tampon de 10 km de rayon autour de chaque point en utilisant la fonction st_buffer().

tampon_orford <- st_buffer(pt_max_orford, dist = 10e3) #nous savons que le CRS est d'unité métrique
tampon_megantic <- st_buffer(pt_max_megantic, dist = 10e3) 
tampon_yamaska <- st_buffer(pt_max_yamaska, dist = 10e3) 
tampon_frontenac <- st_buffer(pt_max_frontenac, dist = 10e3) 

Filtrer le raster dem_lcc pour ne retenir que les cellules comprises à l’intérieur des zones tampons grâce à la fonction mask().

dem_tampon_orford <- mask(dem_lcc, tampon_orford)
dem_tampon_megantic <- mask(dem_lcc, tampon_megantic)
dem_tampon_yasmaka <- mask(dem_lcc, tampon_yamaska)
dem_tampon_frontenac <- mask(dem_lcc, tampon_frontenac)

Combiner les quatre raster en un seul en utilisant la fonction merge()

dem_tampon_max <- merge(dem_tampon_orford, 
                        merge(dem_tampon_megantic,
                              merge(dem_tampon_yasmaka, dem_tampon_frontenac)))

Noter que la fonction merge() s’utilise sur deux rasters à la fois. Il faut donc l’appliquer successivement pour combiner quatre rasters.

Visualiser le raster final contenant les quatre zones tampons:

mapview(dem_tampon_max)