leca-dev · MayaGueguen · Oct 29, 2021 · Nov 23, 2021 · Nov 30, 2021 · Nov 30, 2021
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: RFate
 Type: Package
 Title: FATE with R
-Version: 1.0.1
-Date: 2021-06-14
+Version: 1.0.2
+Date: 2021-09-27
 Authors@R: c(person("Wilfried", "Thuiller", role = "aut"
                     , email = "wilfried.thuiller@univ-grenoble-alpes.fr"),
              person("Damien", "Georges", role = "aut"),
@@ -35,7 +35,7 @@ Description: Wrapper of the C++ model FATE. FATE is a vegetation model based on
         1) gather, prepare and format data to be used with FATE ;
         2) run FATE simulation(s)
         3) process and analyze data produced by FATE.
-RoxygenNote: 7.1.1
+RoxygenNote: 7.1.2
 Encoding: UTF-8
 NeedsCompilation: yes
 License:
@@ -59,6 +59,6 @@ Suggests:
     testthat,
     covr,
     knitr, markdown, rmarkdown
-LinkingTo: Rcpp, RcppThread, BH (>= 1.75.0-0)
+LinkingTo: Rcpp, RcppThread, BH (>= 1.75.0-0), RangeShiftR
 VignetteBuilder: knitr, rmarkdown
 Language: en-US
diff --git a/NAMESPACE b/NAMESPACE
@@ -64,6 +64,7 @@ importFrom(R.utils,gunzip)
 importFrom(R.utils,gzip)
 importFrom(RColorBrewer,brewer.pal)
 importFrom(Rcpp,sourceCpp)
+importFrom(ade4,dudi.hillsmith)
 importFrom(ade4,dudi.pca)
 importFrom(ade4,dudi.pco)
 importFrom(ade4,inertia.dudi)

diff --git a/R/FATE_RS.R b/R/FATE_RS.R
@@ -0,0 +1,305 @@
+### HEADER #####################################################################
+##' @title S
+##' 
+##' @name FATE_RS
+##'
+##' @author Maya Guéguen
+##' 
+##' 
+##' @description ...
+##'              
+##'  
+##' @keywords RangeShifter
+##'  
+##' @export
+##' 
+##' @importFrom RangeShiftR RunRS readPop
+##' @importFrom raster stack rasterFromXYZ writeRaster rasterToPolygons disaggregate rasterize
+##' @importFrom sf st_as_sf st_buffer st_union st_area st_write
+##' @importFrom rmapshaper ms_explode
+##' @importFrom spatstat im predict
+##' @importFrom stringr str_count
+##' @importFrom scales rescale
+##' 
+## END OF HEADER ###############################################################
+
+
+FATE_RS = function(name.simulation, file.simulParam, opt.no_CPU = 1, verbose.level = 2
+                   , PPM.names_PFG, PPM.ras_env, PPM.xy, PPM.mod
+                   , PATCH.threshold = 80, PATCH.buffer = 500, PATCH.min_m2 = 300000
+                   , COST.ras_elev, COST.ras_barr, COST.wei = c(0.43, 0.53, 0.04), COST.range = c(1, 10)
+                   , RS.name_simul, RS.dens_dep, RS.sim, RS.demo, RS.disp, RS.init
+                   , POP.carrying_capacity = 0.115
+                   , year.start, year.end, year.step)
+{
+  ## + argument pour savoir si on commence par FATE ou RS
+
+  ## LOOP over simulation years
+  for (ye in seq(year.start, year.end, 1))
+  {
+    .setParam(params.lines = file.simulParam
+              , flag = "SAVING_DIR"
+              , flag.split = "^--.*--$"
+              , value = paste0(name.simulation, "/RESULTS/SIMUL_YEAR_", ye, ""))
+
+    ## Run FATE ###################################################################################
+    FATE(simulParam = file.simulParam, no_CPU = opt.no_CPU, verboseLevel = verbose.level)
+
+    ## Get resulting FATE vegetation maps #########################################################
+    ## Get results directories ------------------------------------------------
+    GLOB_DIR = .getGraphics_results(name.simulation  = name.simulation
+                                    , abs.simulParam = file.simulParam)
+
+    ## Get raster mask --------------------------------------------------------
+    GLOB_MASK = .getGraphics_mask(name.simulation  = name.simulation
+                                  , abs.simulParam = file.simulParam)
+
+    ## UNZIP the raster saved -------------------------------------------------
+    raster.perPFG.allStrata = .getRasterNames(year.step, "allStrata", "ABUND", GLOB_DIR)
+    .unzip(folder_name = GLOB_DIR$dir.output.perPFG.allStrata
+           , list_files = raster.perPFG.allStrata
+           , no_cores = opt.no_CPU)
+
+    ## GET PFG abundance maps (all strata) ------------------------------------
+    file_name = paste0(GLOB_DIR$dir.output.perPFG.allStrata
+                       , "Abund_YEAR_"
+                       , year.step
+                       , "_"
+                       , PPM.names_PFG
+                       , "_STRATA_all.tif")
+    gp = PPM.names_PFG[which(file.exists(file_name))]
+    file_name = file_name[which(file.exists(file_name))]
+
+    if (length(file_name) > 0)
+    {
+      ras.PFG = stack(file_name) * GLOB_MASK$ras.mask
+      names(ras.PFG) = gp
+    } else
+    {
+      stop(paste0("Missing data!\n The folder "
+                  , GLOB_DIR$dir.output.perPFG.allStrata
+                  , " does not contain adequate files"))
+    }
+
+    ## Combine with other variables ###############################################################
+    ras.all = stack(ras.PFG, PPM.ras_env)
+
+    env.data = as.data.frame(ras.all)
+    env.data = na.exclude(env.data)
+    for(pfg in PPM.names_PFG) {
+      env.data[which(env.data[, pfg] < 0), pfg] = 0
+    }
+
+    ### Create all explanatory variables (with interactions between vegetation and environment)
+    mess = paste0(c(paste0("env.data$", names(PPM.ras_env))
+                    , paste0("sqrt(env.data$", PPM.names_PFG, ")")), collapse = ", ")
+    eval(parse(text = paste0("X.des = poly(", mess, ", degree = 2, raw = TRUE)")))
+    col_power_1 = which(str_count(colnames(X.des), "1") == 1)
+    col_power_2 = which(str_count(colnames(X.des), "1") == 0)
+    col_inter_var1 = grep("^1.", colnames(X.des))
+    col_inter_var2 = grep("^0.1.", colnames(X.des))
+    col_toKeep = sort(unique(c(
+      col_power_1, col_power_2, col_inter_var1, col_inter_var2
+    )))
+    X.des = X.des[, col_toKeep]
+
+    ### Change explanatory variables into images ------------------------------
+    ux = sort(unique(PPM.xy[, 1])) #x ordinates
+    uy = sort(unique(PPM.xy[, 2])) #y ordinates
+    nx = length(ux)
+    ny = length(uy)
+    col.ref = match(PPM.xy[, 1], ux) # indice de l'abscisse de quad dans ux
+    row.ref = match(PPM.xy[, 2], uy)
+
+    int.list = list()
+    for (n in 1:dim(X.des)[2]) {
+      all.vec = rep(NA, max(row.ref) * max(col.ref))
+      vec.ref = (col.ref - 1) * max(row.ref) + row.ref
+      all.vec[vec.ref] = X.des[, n]
+      int.list[[n]] = im(matrix(all.vec, max(row.ref), max(col.ref), dimnames = list(uy, ux)),
+                         xcol = ux,
+                         yrow = uy)
+    }
+    names(int.list) = paste0("V", 1:dim(X.des)[2])
+
+    ## Transform them into RS species distribution map ############################################
+    mod.pred = predict(PPM.mod, covariates = int.list, ngrid = c(ny, nx))
+    ras.pred = raster(mod.pred)
+    ras.quality = ras.pred / max(na.exclude(values(ras.pred))) * 100
+    tmp = log(ras.quality) + abs(min(log(ras.quality)[], na.rm = TRUE))
+    ras.log = tmp / max(tmp[], na.rm = TRUE) * 100
+
+
+    ## Update RS patch maps #######################################################################
+    ras.patch = ras.log > PATCH.threshold # defining patches
+    ras.patch[ras.patch[] < 1] = NA # keep only patches
+
+    pol.patch = rasterToPolygons(ras.patch, dissolve = TRUE) # converting into polygons
+    pol.patch = disaggregate(pol.patch) # separating into multiple polygons
+
+    pol.patch = st_as_sf(pol.patch)
+    pol.patch_buffer = st_buffer(pol.patch, dist = PATCH.buffer)
+    pol.patch_buffer = st_union(pol.patch_buffer) # merging overlapping polygons
+    pol.patch_buffer = ms_explode(pol.patch_buffer) # separating into multiple polygons
+
+    pol.patch_buffer$area_m2 = as.numeric(st_area(pol.patch_buffer))
+    pol.patch_buffer = pol.patch_buffer[which(pol.patch_buffer$area_m2 > PATCH.min_m2)] # excluding too small patches
+
+    ## Change into raster
+    ras.patch_buffer = rasterize(pol.patch_buffer, ras.log)
+    ## Set Habitat Suitability out of patches to NA
+    ras.log[which(is.na(ras.patch_buffer[]))] = NA
+
+    ## SAVE PATCH into ASCII
+    val.buf = ras.patch_buffer[]
+    val.buf[is.na(val.buf)] = -9
+    val.buf = as.character(val.buf)
+    seq_linestarts <- seq(from = nrow(ras.patch_buffer)
+                          , to = ncell(ras.patch_buffer) - nrow(ras.patch_buffer) + 1
+                          , by = nrow(ras.patch_buffer))
+    val.buf[seq_linestarts] = paste0(val.buf[seq_linestarts], "\n") # delineating rows for future character string
+    cha.buf = paste(val.buf, collapse = ' ') # concatenating in 1 character
+    cha.buf = gsub("\n ", "\n", cha.buf) # removing blank from line starts
+
+    file.rename(from = paste0(RS.name_simul, "Inputs/PatchFile.txt")
+                , to = paste0(RS.name_simul, "Inputs/PatchFile_YEAR_", ye - 1, ".txt"))
+
+    fileConn <- file(paste0(RS.name_simul, "Inputs/PatchFile.txt"))
+    writeLines(c(paste("ncols", ncol(ras.patch_buffer), sep = " "),
+                 paste("nrows", nrow(ras.patch_buffer), sep = " "),
+                 paste("xllcorner", extent(ras.patch_buffer)[1], sep = " "),
+                 paste("yllcorner", extent(ras.patch_buffer)[3], sep = " "),
+                 paste("cellsize", resolution(ras.patch_buffer)[1], sep = " "),
+                 "NODATA_value -9",
+                 cha.buf), 
+               fileConn)
+    close(fileConn)
+
+
+    ## SAVE HS RASTER MAP AS ASCII ################################################################
+    val.log = ras.log[]
+    val.log[is.na(val.log)] = -9
+    val.log = as.character(val.log)
+    seq_linestarts <- seq(from = nrow(ras.log), to = ncell(ras.log) - nrow(ras.log) + 1, by = nrow(ras.log))
+    val.log[seq_linestarts] = paste0(val.log[seq_linestarts], "\n") # delineating rows for future character string
+    cha.log = paste(val.log, collapse = ' ') # concatenating in 1 character
+    cha.log = gsub("\n ", "\n", cha.log) # removing blank from line starts
+
+    file.rename(from = paste0(RS.name_simul, "Inputs/LandscapeFile_HabSuit.txt")
+                , to = paste0(RS.name_simul, "Inputs/LandscapeFile_HabSuit_YEAR_", ye - 1, ".txt"))
+
+    fileConn <- file(paste0(RS.name_simul, "Inputs/LandscapeFile_HabSuit.txt"))
+    writeLines(c(paste("ncols", ncol(ras.log), sep = " "),
+                 paste("nrows", nrow(ras.log), sep = " "),
+                 paste("xllcorner", extent(ras.log)[1], sep = " "),
+                 paste("yllcorner", extent(ras.log)[3], sep = " "),
+                 paste("cellsize", resolution(ras.log)[1], sep = " "),
+                 "NODATA_value -9",
+                 cha.log), 
+               fileConn)
+    close(fileConn)
+
+
+    ## Update RS disp-cost maps ###################################################################
+    ## Function to rescale habitat suitability maps
+    fun_rescale = function(a = 1, b = 100, ras, ext) {
+      mini = min(values(ras), na.rm = TRUE)
+      maxi = max(values(ras), na.rm = TRUE)
+      res = (((b - a) * (ras - mini)) / (maxi - mini)) + a
+      res = mask(res, ext)
+      return(res)
+    }
+
+    ras.log_ = fun_rescale(a = a, b = b, ras = ras.log, ext = ras.log)
+    ras.HS = fun_rescale(a = a, b = b, ras = 100 + (-1 * ras.log_) ^ 1, ext = ras.log_)
+
+    # Compute cost map --------------------------------------------------------
+    ras.cost <- COST.wei[1] * ras.HS + COST.wei[2] * COST.ras_elev + COST.wei[3] * COST.ras_barr
+    ras.cost[] <- rescale(ras.cost[], c(COST.range[1], COST.range[2]))
+    ras.cost[] <- round(ras.cost[]) # round for SMS module
+
+    ## SAVE COST RASTER MAP AS ASCII ##############################################################
+    val.cost = ras.cost[]
+    val.cost[is.na(val.cost)] = -9
+    val.cost = as.character(val.cost)
+    seq_linestarts <- seq(from = nrow(ras.cost), to = ncell(ras.cost) - nrow(ras.cost) + 1, by = nrow(ras.cost))
+    val.cost[seq_linestarts] = paste0(val.cost[seq_linestarts], "\n") # delineating rows for future character string
+    cha.cost = paste(val.cost, collapse = ' ') # concatenating in 1 character
+    cha.cost = gsub("\n ", "\n", cha.cost) # removing blank from line starts
+
+    file.rename(from = paste0(RS.name_simul, "Inputs/LandscapeFile_CostMap.txt")
+                , to = paste0(RS.name_simul, "Inputs/LandscapeFile_CostMap_YEAR_", ye - 1, ".txt"))
+
+    fileConn <- file(paste0(RS.name_simul, "Inputs/LandscapeFile_CostMap.txt"))
+    writeLines(c(paste("ncols", ncol(ras.cost), sep = " "),
+                 paste("nrows", nrow(ras.cost), sep = " "),
+                 paste("xllcorner", extent(ras.cost)[1], sep = " "),
+                 paste("yllcorner", extent(ras.cost)[3], sep = " "),
+                 paste("cellsize", resolution(ras.cost)[1], sep = " "),
+                 "NODATA_value -9",
+                 cha.cost), 
+               fileConn)
+    close(fileConn)
+
+
+    ## Run RangeShifter ###########################################################################
+    RS.land = ImportedLandscape(LandscapeFile = paste0(RS.name_simul, "Inputs/LandscapeFile_HabSuit.txt"),
+                                  Resolution = resolution(ras.log), 
+                                  HabPercent = TRUE, # continuous values for habitat quality (and not discrete)
+                                  K_or_DensDep = RS.dens_dep,
+                                  PatchFile = paste0(RS.name_simul, "Inputs/PatchFile.txt"),
+                                  CostsFile = paste0(RS.name_simul, "Inputs/LandscapeFile_CostMap.txt"))
+
+    RS.params = RSsim(simul = RS.sim, land = RS.land, demog = RS.demo, dispersal = RS.disp, init = RS.init)
+
+    RunRS(RSparams = RS.params, dirpath = RS.name_simul)
+
+
+    ## Get resulting RS species density ###########################################################
+    # read population output file into a dataframe
+    pop_df = readPop(s = RS.params, dirpath = RS.name_simul)
+    pop_df = pop_df[which(pop_df$Rep == 0 & pop_df$Year == 1), ]
+
+    file.rename(from = paste0(RS.name_simul, "Outputs/Batch1_Sim1_Land1_Pop.txt")
+                , to = paste0(RS.name_simul, "Outputs/Batch1_Sim1_Land1_Pop_YEAR_", ye, ".txt"))
+
+    # Make stack of different raster layers for each year and for only one repetition (Rep==0):
+    stack_years_rep0 = ras.patch_buffer
+    stack_years_rep0[which(!is.na(stack_years_rep0[]))] = 0
+    for (pp in 1:nrow(pop_df)) {
+      stack_years_rep0[which(ras.patch_buffer[] == pop_df$PatchID[pp])] = pop_df$NInd[pp] / pop_df$Ncells[pp]
+    }
+    stack_years_rep0 = stack_years_rep0 / POP.carrying_capacity
+    stack_years_rep0[which(stack_years_rep0[] > 1)] = 1
+    name.dist = paste0(name.simulation, "/DATA/MASK/MASK_DIST_YEAR_", ye + 1, ".tif")
+    writeRaster(stack_years_rep0, filename = name.dist)
+
+    ## Update FATE disturbance maps ###############################################################
+    dist.files = .getParam(params.lines = file.simulParam
+                           , flag = "DIST_MASK"
+                           , flag.split = "^--.*--$"
+                           , is.num = FALSE)
+    dist.files[1] = paste0(name.simulation, "/DATA/MASK/MASK_DIST_YEAR_", ye + 1, ".tif")
+
+    .setParam(params.lines = file.simulParam
+              , flag = "DIST_MASK"
+              , flag.split = "^--.*--$"
+              , value = dist.files)
+
+
+    ## Update FATE simulation to load #############################################################
+    year.save = readLines(.getParam(params.lines = file.simulParam
+                                    , flag = "SAVING_YEARS_OBJECTS"
+                                    , flag.split = "^--.*--$"
+                                    , is.num = FALSE))
+    file.rename(from = paste0(GLOB_DIR$dir.save, "SimulMap_", year.save[1], ".sav")
+                , to = paste0(GLOB_DIR$dir.save, "SimulMap_STEP_YEAR_", ye, ".sav"))
+
+    .setParam(params.lines = file.simulParam
+              , flag = "SAVED_STATE"
+              , flag.split = "^--.*--$"
+              , value = paste0(GLOB_DIR$dir.save, "SimulMap_STEP_YEAR_", ye, ".sav"))
+
+  }
+}
+
diff --git a/R/PRE_FATE.params_globalParameters.R b/R/PRE_FATE.params_globalParameters.R
@@ -124,8 +124,6 @@
 ##' @param doAliens default \code{FALSE}. \cr If \code{TRUE}, invasive plant 
 ##' introduction is activated in the \code{FATE} simulation, and associated 
 ##' parameters are required
-##' @param ALIEN.no (\emph{optional}) \cr an \code{integer} corresponding to the 
-##' number of introductions
 ##' @param ALIEN.freq (\emph{optional}) \cr a \code{vector} of \code{integer} 
 ##' corresponding to the frequency of each introduction (\emph{in years})
 ##' @param doFire default \code{FALSE}. \cr If \code{TRUE}, fire 
@@ -435,7 +433,6 @@
 ##'   introduction areas. \cr If the habitat suitability filter is on, 
 ##'   suitability maps will also be needed for these new groups.
 ##'   \describe{
-##'     \item{ALIEN.no}{the number of different introductions}
 ##'     \item{ALIEN.freq}{the frequency of each introduction (\emph{in years})}
 ##'   }
 ##'   }
@@ -642,7 +639,6 @@
 ##' 
 ##' \itemize{
 ##'   \item DO_ALIENS_INTRODUCTION
-##'   \item ALIENS_NO
 ##'   \item ALIENS_FREQ
 ##' }
 ##' 
@@ -768,7 +764,6 @@ PRE_FATE.params_globalParameters = function(
   , doDrought = FALSE
   , DROUGHT.no_sub = 4
   , doAliens = FALSE
-  , ALIEN.no
   , ALIEN.freq = 1
   , doFire = FALSE
   , FIRE.no
@@ -862,11 +857,10 @@ PRE_FATE.params_globalParameters = function(
   }
   if (doAliens)
   {
-    .testParam_notInteger.m("ALIEN.no", ALIEN.no)
     .testParam_notInteger.m("ALIEN.freq", ALIEN.freq)
-    if (length(ALIEN.freq) != ALIEN.no){
+    if (length(ALIEN.freq) != required.no_PFG){
       stop(paste0("Wrong type of data!\n `ALIEN.freq` must contain as many "
-                  , "values as the number of introductions (`ALIEN.no`)"))
+                  , "values as the number of PFG (`required.no_PFG`)"))
     }
   }
   if (doFire)
@@ -1012,10 +1006,8 @@ PRE_FATE.params_globalParameters = function(
   if (doAliens)
   {
     params.ALIEN = list(as.numeric(doAliens)
-                        , ALIEN.no
                         , ALIEN.freq)
     names.params.list.ALIEN = c("DO_ALIENS_INTRODUCTION"
-                                , "ALIENS_NO"
                                 , "ALIENS_FREQ")
   } else
   {