In ReseauBiodiversiteQuebec/tableaudelta: Tableau-delta

knitr::opts_chunk$set(echo = TRUE)

library(sf)
library(terra)
library(tidyverse)

On charge une shapefile du Québec, notre zone de référence.

mask <- terra::vect("./data/quebec_nad83.shp")
study_area_ascott <- terra::vect("./data/Contour_boise.shp")
study_area_ascott <- terra::project(study_area_ascott, mask)

On charge les données du réseau routier du Québec

routes <- terra::vect("./data/ESRI(SHP)/Reseau_routier.shp")
routes <- terra::project(routes, mask)

routes_sf <- read_sf("./data/ESRI(SHP)/Reseau_routier.shp")
routes_sf <- st_transform(routes_sf, st_crs(mask))

0. Rasteriser les routes

# Avec gadal_rasterize()
## load the packages
require("rgdal")
devtools:::install_github("gearslaboratory/gdalUtils")
require("gdalUtils")

## Prep le raster 1km/1km
mask_r <- terra::rast(mask)
terra::res(mask_r) <- c(1000, 1000)
terra::values(mask_r) <- 1
mask_r <- terra::mask(mask_r, mask)

## Sauver le raster mask
raster::writeRaster(raster::raster(mask_r), filename="./data/test.tif", overwrite=TRUE)

## Resterize
routes_r <- gdalUtils::gdal_rasterize(burn=c(0,0,128), src_datasource="./data/ESRI(SHP)/Reseau_routier.shp", dst_filename="./data/test.tif", output_Raster=TRUE)

## Cellules 0 = pas de routes
routes_r <- 1 - routes_r

# Sauver le raster
raster::writeRaster(routes_r, filename="./data/routes_pres_absence.tif", overwrite=TRUE)

1. Calcul de la densité de routes sur la zone de référence

table(routes$ClsRte)

# Selection des routes achalandées
# - Ne comprends pas les routes forestières et isolées
routesOfI <- c("Artère", "Autoroute", "Collectrice de transit", "Collectrice municipale", "Locale", "Nationale", "Régionale", "Rue piétonne")
i <- which(routes$ClsRte %in% routesOfI)
routesSel <- routes[i,]

2. Calcul du stresseur sur une zone d’étude

Le calcul de la densité de route par cellule est computationellement lourd. Pour le bien de l'exploration, les calaculs sont limités à la zone d'étude.

study_area <- terra::crop(routesSel, study_area_ascott)
terra::plot(study_area)

Calcul de la densité de route

# Création de la grille
# - cellules de 1000m x 1000m
library(dplyr)
grid <- sf::st_make_grid(sf::st_bbox(study_area_ascott),
                         cellsize = 1000, square = TRUE) %>%
  sf::st_as_sf() %>%
  dplyr::mutate(cell = 1:nrow(.)) %>%
  sf::st_as_sf()

grid |> terra::plot()
terra::plot(study_area, add = TRUE)

split_count <- sf::st_intersection(sf::st_as_sf(study_area), grid) %>%
                  group_by(cell) %>%
                  count() %>%
                  arrange(desc(n))

split_count$geometry <- NULL
RteDst <- dplyr::left_join(grid, split_count) %>% 
  sf::st_sf(sf_column_name = "x")

RteDst$n[is.na(RteDst$n)] <- 0

# Normaliser la densité
RteDst$n <- (RteDst$n - min(RteDst$n, na.rm = TRUE))/(max(RteDst$n, na.rm = TRUE) - min(RteDst$n, na.rm = TRUE))
RteDst <- terra::vect(RteDst)

# Plot
color <- colorRampPalette(c("grey90", "steelblue4", "steelblue2", 
                            "steelblue1", "gold", "red1", "red4"),
                          bias = 1)
fixBreak <- seq(0, 1, length = 99) |> round(2)

terra::plot(RteDst, "n", color(length(fixBreak)), breaks = fixBreak)
terra::plot(study_area, add = TRUE)