In ReseauBiodiversiteQuebec/tableaudelta: Tableau-delta

knitr::opts_chunk$set(echo = TRUE)
library("gdalcubes")
library("rstac")
library("tibble")
library("sp")
library("sf")
library("dplyr")
library("rgbif")
library("tidyr")
library("stars")
library("terra")
library("stacatalogue")
#devtools::install_github("ReseauBiodiversiteQuebec/stac-catalogue")

shp_to_bbox <- function(shp, proj_from = NULL, proj_to = NULL) {
  if(is.na(sf::st_crs(shp)) && is.null(proj_from)) {
    stop("proj.fom is null and shapefile has no crs.")
  }

  if(is.na(sf::st_crs(shp))) {
    crs(shp) <- proj_from
    shp <- shp %>% sf::st_set_crs(proj_from)
  }

  if (!is.null(proj_to) ) {
    shp <- shp %>% sf::st_as_sf() %>%
      sf::st_transform(crs = sp::CRS(proj_to))
  }


  bbox <- sf::st_bbox(shp, crs = proj)

  bbox
}
fast_crop <- function(predictors,
                       mask) {

  # convert into terra raster for increased speed
  if (!inherits(predictors, "SpatRaster")) {
    predictors <- terra::rast(predictors)
  }

  # convert into a SpatVector
 if (!inherits(mask, "SpatVector")) {
    mask <- terra::vect(mask)
  }
predictors <- terra::crop(predictors, mask)
  predictors <- terra::mask(predictors, mask, touches = FALSE)

  # convert to raster format for later processing
     return(predictors)
}

climate_velocity <- function(current_tmean, future_tmean, time_span,
                             type = "local",  opt="slope", units = "meters", 
                             neighbors=8) { 
  if (type == "local") {

    # 1. Spatial gradient
    # Neighborhood Slope Algorithm, average maximum technique
    spatial_gradient <- raster::terrain(current_tmean,
                                        opt = opt,
                                        units = units,
                                        neighbors=neighbors  # (queen case) 
    )
    # Truncating zero values
    spatial_gradient[spatial_gradient <= 0.00001] <- 0.00001
  }

  # 2. Temporal gradient
  temporal_gradient <- (future_tmean - current_tmean)/time_span

  # 3. Velocity
  local_velocity <- temporal_gradient/spatial_gradient

  return(local_velocity)
}

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

mask <- terra::vect("C:/GitHub/tableaudelta/data/shape_study_area_nolakes_nad83.shp")
study_area_ascott <- terra::vect("C:/GitHub/tableaudelta/data/Contour_boise.shp")
study_area_ascott <- terra::project(study_area_ascott, mask)
bbox <- shp_to_bbox(mask)
srs.cube <- "EPSG:6623"

1. Calcul du stresseur sur la zone de référence

Exemple de la vélocité climatique = difference between future and current temperature divided by number of years (°C/year) / slope of temperature across spatial neighborhood (meters/°C). Ici on utilise la définition de la vélocité climatique locale Sandel et al. 2011.

Température moyenne actuelle

Pour les fonctions permettant d'intéragir avec le stac catalogue, voir le package stacatalogue.

layers <- c("bio1")

cube <- 
  load_cube(stac_path = "https://io.biodiversite-quebec.ca/stac/",
            limit = 5000, 
            collections = c("chelsa-clim"), 
            use.obs = F,
            obs = obs.coords.proj,
            bbox = bbox,
            buffer.box = 0,
            layers = layers,
            srs.cube = srs.cube,
            t0 = "1981-01-01",
            t1 = "1981-01-01",
            spatial.res = 1000, # in meters
            temporal.res = "P1Y",
            aggregation = "mean",
            resampling = "near") 

current_tmean <- cube_to_raster(cube, format = "terra")
current_tmean <- fast_crop(current_tmean, mask)
current_tmean <- (current_tmean$bio1/10) - 273.15 #on convertit en degrés celsius
plot(current_tmean)

Température moyenne actuelle

Ici, on choisit un scénario d'émissions "pessimiste" (ssp585) et la période 2040-2070. Les autres options possibles sont indiquées en commentaires. On utilise la moyenne d'un ensemble de modèles climatiques.

cube.future <- load_cube_projection(stac_path = 
                                      "https://io.biodiversite-quebec.ca/stac/",
                                    limit = 5000, 
                                    collections = c('chelsa-clim-proj'), 
                                    use.obs = F,
                                    obs = NULL,
                                    buffer.box = 0,
                                    rcp = 'ssp585', #ssp126, ssp370, ssp585
                                    bbox = bbox,
                                    layers = NULL,
                                    variable = "bio1",
                                    srs.cube = srs.cube, 
                                    time.span = "2041-2070", #"2011-2040", 2041-2070 or 2071-2100
                                    spatial.res = 1000, 
                                    temporal.res  = "P1Y",
                                    aggregation = "mean",
                                    resampling = "near") 
future_tmean <- cube_to_raster(cube.future, format = "raster")
future_tmean <- (raster::calc(future_tmean, mean)/10) - 273.15
future_tmean <- fast_crop(future_tmean, mask)
plot(future_tmean)

Local climate velocity

c_velocity <- climate_velocity(raster::raster(current_tmean), raster::raster(future_tmean), time_span = 2060-1980, type = "local",  opt="slope", units = "meters", 
                             neighbors=8)

plot(c_velocity)

Normalisation de la distribution entre 0-1

c_velocity <- terra::rast(c_velocity)
c_velocity <- (c_velocity -terra::minmax(c_velocity)[1,])/(terra::minmax(c_velocity)[2,]-terra::minmax(c_velocity)[1,])

plot(c_velocity)

On récupère les valeurs.

c_velocity_values <- as.data.frame(c_velocity, row.names=NULL, optional=FALSE, xy=FALSE, na.rm=T)$layer
summary(c_velocity_values)

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

#study_area <- as(raster::extent(-388691.1, 141649.6,   565217.2,  1037335.6), 'SpatialPolygons')
#crs(study_area) <- srs.cube
study_area <- crop(c_velocity, study_area_ascott, snap = "out")
plot(study_area)

On calcule la moyenne du stresseur sur la zone.

study_area_val <- as.data.frame(study_area, row.names=NULL, optional=FALSE, xy=FALSE, na.rm=T)
study_area_mean <- mean(study_area_val$layer)
study_area_mean

On regarde son positionnement relatif par rapport au Québec.

ecdf_fun <- function(x,perc) ecdf(x)(perc)
ecdf_fun(c_velocity_values, study_area_mean)

ReseauBiodiversiteQuebec/tableaudelta documentation built on May 25, 2022, 5:26 p.m.