#> R-Script "analysis.R" for #> a) making predictions with the SIM and SDM on new climate data, and #> b) converting the Site index into top height, as well as #> c) converting occurrence probabilities into climate suitability classes, #> according to the table of SDM thresholds #> #> INPUT data from Woehlbrandt et al. (in prep.): #> 1) Site index model predictions #> 2) Table nlrq-parameters #> 3) Species distribution model predictions #> 4) Table SDM thresholds # 00. packages, paths, functions ############################################### require(mgcv) require(terra) ## paths --------------------------------------------------------------------- # path <- "E:/00_Benutzer/01_Projekte/03_EVA/EVA Ergebnisse/DIB_paper" # file.choose() path_SIM_predictions <- paste0(path, "/SIM_predictions") path_SIM_predictions <- "//10.46.52.95/Datenpool3/Waldbau/03_Projekte/03_Modellierung/W1 Eva Evidenzbas Anbauempf/2_Daten/03_EVA_AW/results/02_SIM/GAMs/Run_20240523/data/EUpred/tifs" path_SDM_predictions <- paste0(path, "/SDM_predictions") path_SIMs <- paste0(path, "/SIMs") path_SDMs <- paste0(path, "/SDMs") path_nlrq_file <- paste0(path, "/nlrq-parameters.xlsx") path_SDM_thresholds <- paste0(path, "/SDM thresholds.xlsx") dir(path_SIM_predictions) # check ## functions ----------------------------------------------------------------- # quantile_boundaries <- function(nlrq_file, species = "Fagus sylvtica", age = 50) { #browser() nlrq_file <- nlrq_file[nlrq_file$species == species, ] H0pred05_age = nlrq_file$A_low * (1 - exp(- nlrq_file$k_low * age)) ^ nlrq_file$p_low # age = i H0pred95_age = nlrq_file$A_up * (1 - exp(- nlrq_file$k_up * age)) ^ nlrq_file$p_up # age = i return(data.frame(H0pred05 = H0pred05_age, H0pred95 = H0pred95_age)) } convert_SI_into_top_height <- function(tif_file, nlrq_file, species = "Fagus sylvatica", age = 100) { names(tif_file) <- "predict" df <- quantile_boundaries(nlrq_file, species = species, age = age) tif_file <- (df$H0pred95 - df$H0pred05) * tif_file + df$H0pred05 tif_file[tif_file$predict < 0] <- 0 return(tif_file) } create_thresholds_SDMs <- function(SDM_prediction, species = "Fagus sylvatica", SDM_quantiles = c(.05, .3), path_SDMs) { # load SDM load(file = paste0(path_SDMs, "/", species, ".RData")) # load model # transform SDM_quantiles in SDM_thresholds Trans = function(x){exp(x)/(1 + exp(x))} myPred = Trans(predict(themod)[which(themod$model[1] == 1)]) # perform a prediction using only the presences of the calibration data rcl = c(0, quantile(myPred, SDM_quantiles), 1) # rcl = reclassify predictions using quantiles return(rcl) } # 01. Analysis ################################################################# nlrq_file <- readxl::read_excel(path_nlrq_file) species <- "Fagus sylvatica" time_period <- "1981_2010" # or 2011_2040, 2041_2070, 2071_2100 scenario <- "ref" # or rcp26, rcp45, rcp85 ## a) SIM and SDM predictions on new climate data ------------------------------ # load SIM or SDM themod <- paste(path_SIMs, # SIMs dir(path_SIMs)[grepl(species, dir(path_SIMs))], sep = "/") load(paste(path_SDMs, dir(path_SDMs)[1], sep = "/")) # SDM #> load your climate raster or polygons #> important here: you need the aggregated temperature and precipitation values #> of spring (Mar-May), summer (Jun-Aug) and winter (Dec-Feb), as well as the #> mean annual air temperature (bio1) and annual precipitation sum (bio12) #> => sadly, the variables are differently named within the SIMs and SDMs #> SIM SDM #> "reference_19812010_Bio." | "Bio." # 1 or 12 #> "reference_19812010_su_p" | "Bio.18" #> "reference_19812010_wi_t" | "Bio.11" #> "reference_19812010_su_t" | "Bio.10" #> "reference_19812010_sp_t" | "sp_t" #> "reference_19812010_sp_p" | "sp_p" #> "reference_19812010_wi_p" | "wi_p" # However, you have to name your climate variables exactly the same, such as themod <- readRDS(paste(path_SIMs, dir(path_SIMs)[grepl(species, dir(path_SIMs))], sep = "/")) # SIM head(themod$model[-1]) # then you use the predict.gam-function from "mgcv"-package pre <- mgcv::predict.gam(themod, newdata = head(themod$model[-1]), type = "response") #, se.fit = TRUE pre # and do further analysis or visualization with it .. cbind(head(themod$model[-1]), predict = pre) ## b) Site index into top height ----------------------------------------------- # load .tif file SIM_prediction_SI <- terra::rast(paste(path_SIM_predictions, dir(path_SIM_predictions)[grepl(paste(gsub(" ", "_", species), time_period, scenario, sep = "_"), dir(path_SIM_predictions))], sep = "/")) # caution: this function is written for raster files! SIM_prediction_top_height <- convert_SI_into_top_height(tif_file = SIM_prediction_SI, nlrq_file, species, age = 100) plot(SIM_prediction_top_height) ## c) occurence probability into climate suitability classes ------------------- ### thresholds --------------------------------------------------------------- # calculate your own thresholds SDM_prediction <- terra::rast(paste(path_SDM_predictions, dir(path_SDM_predictions)[grepl(paste(gsub(" ", "_", species), time_period, scenario, sep = "_"), dir(path_SDM_predictions))], sep = "/")) rcl <- create_thresholds_SDMs(SDM_prediction, species, SDM_quantiles = c(.05, .3), # here specify quantiles path_SDMs) # or get those from the table of SDM thresholds SDM_thresholds <- readxl::read_excel(path_SDM_thresholds) SDM_thresholds <- SDM_thresholds[SDM_thresholds$species == species, ] # filter for species rcl <- as.numeric(unlist(strsplit(SDM_thresholds$thresholds, " "))) # filter for thresholds and format accordingly ### apply thresholds --------------------------------------------------------- # and apply thresholds SDM_prediction[SDM_prediction < rcl[3]] <- NA # for example, to build a mask for the SIM_predictions .. SDM_prediction[!is.na(SDM_prediction)] <- 1 plot(SDM_prediction) # then mask the SIM_prediction_top_height: plot(SIM_prediction_top_height * SDM_prediction)