diff --git a/R/Climate.R b/R/Climate.R
index dc02b9cdf0e11ddcd13670aff9b8219641bfef14..6de9529b4bf4938e0f5be43ad16e1f65fdd18428 100644
--- a/R/Climate.R
+++ b/R/Climate.R
@@ -11,6 +11,7 @@
#' @param ... additional param not use for the moment
#' @import archive raster
download_climatologies <- function(name = "Bio01",path = '', overwrite = FALSE, ...) {
+require(archive)
if(length(name)>1) stop("\nname must be of length one")
if(!name %in% paste0("Bio", c(paste0("0", 1:9), 10:19))) stop(paste0("\nArgument name must be in ",
paste(paste0("Bio", c(paste0("0", 1:9), 10:19)),
@@ -109,6 +110,27 @@ return(res)
}
+# Functions to extract soil data from SoilGrid
+## TODO change for ftp://ftp.soilgrids.org/data/aggregated/1km/ for sl1 and sl2 and take the mean
+
+download_ph <- function(dir_temp = "soil_temp"){
+
+require(raster)
+require(R.utils)
+# download raster
+url_ph <- "ftp://ftp.soilgrids.org/data/aggregated/1km/PHIHOX_M_sl2_1km_ll.tif"
+raster_name <- "PHIHOX_M_sl1_1km_ll.tif"
+if(!file.exists(file.path(dir_temp, raster_name))){
+ dir.create(dir_temp)
+ download.file(url_ph, file.path(dir_temp, raster_name),mode="wb")
+}
+# load raster
+raster_ph <- raster(file.path(dir_temp, raster_name))
+e1 <- extent(-30, 60, 30, 80)
+res <- crop(raster_ph, e1)
+return(res)
+}
+
#######################
## Function to extract annual AET from http://www.cgiar-csi.org/data/global-high-resolution-soil-water-balance
download_AETy <- function(dir_temp = "clim_temp"){
@@ -134,11 +156,11 @@ return(res)
## stack data
stack_list <- function(Bio01, Bio04, Bio10, Bio11, Bio15, Bio16,
- aridity, alpha, AETy){
+ aridity, alpha, AETy, Ph){
listt <- list(Bio01, Bio04, Bio10, Bio11, Bio15, Bio16)
require(raster)
e <- extent(-30, 60, 30, 80)
res <- crop(raster::stack(listt), e)
-res <- raster::stack(res, aridity, alpha, AETy)
+res <- raster::stack(res, aridity, alpha, AETy, Ph)
return(res)
}
diff --git a/R/FormatData.R b/R/FormatData.R
index 03a6ab74f3734471b4d04d2e843c4963c17d6bd0..99541bfae0d17b206b39c29612b34a312f34bbd2 100644
--- a/R/FormatData.R
+++ b/R/FormatData.R
@@ -26,7 +26,7 @@ return(res)
extract_clim <- function(points, l_stack){
val <- raster::extract(x=l_stack,y=points) # Extraction of the values
colnames(val) <- c("Bio01", "Bio04", "Bio10", "Bio11",
- "Bio15", "Bio16", "aridity", "alpha", "AETy")
+ "Bio15", "Bio16", "aridity", "alpha", "AETy", "Ph")
df <- cbind(data.frame(points),val)
return(df)
}
diff --git a/R/SDM.R b/R/SDM.R
index aba74d06eeb4e660760beefc529c184f23906cb7..714904b0e5507753b8bfafa151b18ec56d0f72fc 100644
--- a/R/SDM.R
+++ b/R/SDM.R
@@ -2,7 +2,12 @@
### FUNCTION TO FIT SDM
BIOMOD_DATA_FOR_SP <- function(sps="Fagus.sylvatica", df){
+ require(dplyr)
+ print(paste("start species ", sps))
+ wd <- getwd()
+ dir.create("biomod2", showWarnings = FALSE)
if (!file.exists(file.path("output", paste0(sps,".rds")))){
+ setwd(file.path(wd, "biomod2"))
## todo sample plot to limit spatial autocorrelation
## caleval <- ecospat.caleval(data = as.numeric(df[[sps]]),
## xy = df[2:3], row.num = 1:nrow(df),
@@ -12,53 +17,52 @@ BIOMOD_DATA_FOR_SP <- function(sps="Fagus.sylvatica", df){
expl.var = df[c("Bio01", "Bio04",
"Bio15",
"Bio16", "aridity",
- "AETy")],
+ "AETy", "Ph")],
resp.xy = coordinates(df),
resp.name = sps)
myBiomodOption <- BIOMOD_ModelingOptions()
myBiomodModelOut <- BIOMOD_Modeling(myBiomodData,
- models = c('GLM','GBM','RF'),
+ models = c('GLM','GBM','RF'),#'GLM','GBM','RF',
models.options = myBiomodOption,
- NbRunEval=1,
+ NbRunEval=10,
DataSplit=70,
- Prevalence=0.5,
+ Prevalence=NULL,
VarImport=0,
models.eval.meth = c('TSS','ROC'),
SaveObj = TRUE,
rescal.all.models = TRUE,
do.full.models = FALSE,
modeling.id = paste(sps,
- "FirstModeling",sep=""))
+ "FirstModeling",sep=""))
myBiomodModelEval <- get_evaluations(myBiomodModelOut)
- #get_variables_importance(myBiomodModelOut)
- myBiomodEM <- BIOMOD_EnsembleModeling(modeling.output = myBiomodModelOut,
- chosen.models = 'all',
- em.by='all',
- eval.metric = c('TSS'),
- eval.metric.quality.threshold = c(0.55),
- prob.mean = T,
- prob.cv = T,
- prob.ci = T,
- prob.ci.alpha = 0.05,
- prob.median = T,
- committee.averaging = T,
- prob.mean.weight = T,
- prob.mean.weight.decay = 'proportional' )
+ ## TODO AFTER PROCESSING ALL SPECIES
+ ## myBiomodEM <- BIOMOD_EnsembleModeling(modeling.output = myBiomodModelOut,
+ ## chosen.models = 'all',
+ ## em.by='all',
+ ## eval.metric = c('TSS'),
+ ## eval.metric.quality.threshold = c(0.55),
+ ## prob.mean = T,
+ ## prob.cv = T,
+ ## prob.ci = T,
+ ## prob.ci.alpha = 0.05,
+ ## prob.median = T,
+ ## committee.averaging = T,
+ ## prob.mean.weight = T,
+ ## prob.mean.weight.decay = 'proportional' )
print(paste(sps, " DONE"))
-
- saveRDS(list(myBiomodModelOut, myBiomodModelEval, myBiomodEM),
+ setwd(wd)
+ saveRDS(list(out = myBiomodModelOut, eval = myBiomodModelEval),
file = file.path("output", paste0(sps,".rds")))
- return(list(myBiomodModelOut, myBiomodModelEval, myBiomodEM))
- }else{
- list_res <- readRDS(file.path("output", paste0(sps,".rds")))
- return(list_res)
}
+ print(paste("done species ", sps))
+
}
BIOMOD2_ALL_SP <- function(df){
sp_vec <- c("Pinus.sylvestris", "Pinus.pinaster", "Picea.abies", "Pinus.nigra",
- "Pinus.halepensis", "Quercus.ilex", "Fagus.sylvatica", "Quercus.robur",
+ "Pinus.halepensis", "Quercus.ilex", "Fagus.sylvatica",
+ "Quercus.robur",
"Betula.pendula", "Betula.pubescens", "Quercus.petraea", "Quercus.pubescens", "Quercus.pyrenaica",
"Quercus.suber", "Pinus.pinea", "Abies.alba", "Carpinus.betulus",
"Pinus.mugo", "Fraxinus.excelsior", "Quercus.faginea", "Alnus.glutinosa",
diff --git a/README.md b/README.md
index 9afa22f6f348529f55f54bbdedf1a2e5458a04d1..d8d2fe87b3b39a1a17a284d2ea9c548f21865463 100644
--- a/README.md
+++ b/README.md
@@ -2,10 +2,11 @@
Based on harmonized NFI data from Mauri et al. (2017), with climatic data extracted from Chelsa and CGIAR
-http://www.cgiar-csi.org/data/global-aridity-and-pet-database
-http://www.cgiar-csi.org/data/global-high-resolution-soil-water-balance
+* http://chelsa-climate.org/
+* http://www.cgiar-csi.org/data/global-aridity-and-pet-database
+* http://www.cgiar-csi.org/data/global-high-resolution-soil-water-balance
-The analysis are done using the R package remake see https://github.com/richfitz/remake to facilitate replicability.
+The analysis are done using the R package 'remake' see https://github.com/richfitz/remake to facilitate replicability.
# References
diff --git a/remake.yml b/remake.yml
index 7a6cb7b79a37310d134c5ab63c1ce141f8607328..f7d5d231a25d396708bbf30ffb397723e251ad15 100644
--- a/remake.yml
+++ b/remake.yml
@@ -5,6 +5,7 @@ packages:
- ggplot2
- raster
- sp
+ - archive
sources:
- R
@@ -54,14 +55,17 @@ targets:
depends:
- download/EUForest.csv
+ sdm_all:
+ command: BIOMOD2_ALL_SP(EUForestClim)
+
EUForest:
command: fomat_species_data("download/EUForest.csv")
EUForestClim:
command: extract_clim(EUForest, BioStack)
- sdm_all:
- command: BIOMOD2_ALL_SP(EUForestClim)
+ Ph:
+ command: download_ph()
Bio01:
command: download_climatologies(I(name = "Bio01"))
@@ -82,7 +86,7 @@ targets:
command: download_climatologies(I(name = "Bio16"))
BioStack:
- command: stack_list(Bio01, Bio04, Bio10, Bio11, Bio15, Bio16, aridity, alpha, AETy)
+ command: stack_list(Bio01, Bio04, Bio10, Bio11, Bio15, Bio16, aridity, alpha, AETy, Ph)
alpha:
command: download_alpha()