BIOMOD2 : How to do projection with models built in a previous stage.
Damien G.
02/12/2014
The aim of this script is to show how to do projections in biomod2 with models coming from a previous modelling campain.
This document is based on a simple code example. All comments, suggestions,… are more than welcome!
library(biomod2)## Loading required package: sp
## Loading required package: raster
## Loading required package: parallel
## Loading required package: reshape
## Loading required package: ggplot2
## biomod2 3.1-63 loaded.
##
## Type browseVignettes(package='biomod2') to access directly biomod2 vignettes.## model GuloGulo niche
example(BIOMOD_Modeling)##
## BIOMOD> # species occurrences
## BIOMOD> DataSpecies <- read.csv(system.file("external/species/mammals_table.csv",
## BIOMOD+ package="biomod2"))
##
## BIOMOD> head(DataSpecies)
## X X_WGS84 Y_WGS84 ConnochaetesGnou GuloGulo PantheraOnca
## 1 1 -94.5 82 0 0 0
## 2 2 -91.5 82 0 1 0
## 3 3 -88.5 82 0 1 0
## 4 4 -85.5 82 0 1 0
## 5 5 -82.5 82 0 1 0
## 6 6 -79.5 82 0 1 0
## PteropusGiganteus TenrecEcaudatus VulpesVulpes
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0
##
## BIOMOD> # the name of studied species
## BIOMOD> myRespName <- 'GuloGulo'
##
## BIOMOD> # the presence/absences data for our species
## BIOMOD> myResp <- as.numeric(DataSpecies[,myRespName])
##
## BIOMOD> # the XY coordinates of species data
## BIOMOD> myRespXY <- DataSpecies[,c("X_WGS84","Y_WGS84")]
##
## BIOMOD> # Environmental variables extracted from BIOCLIM (bio_3, bio_4, bio_7, bio_11 & bio_12)
## BIOMOD> myExpl = stack( system.file( "external/bioclim/current/bio3.grd",
## BIOMOD+ package="biomod2"),
## BIOMOD+ system.file( "external/bioclim/current/bio4.grd",
## BIOMOD+ package="biomod2"),
## BIOMOD+ system.file( "external/bioclim/current/bio7.grd",
## BIOMOD+ package="biomod2"),
## BIOMOD+ system.file( "external/bioclim/current/bio11.grd",
## BIOMOD+ package="biomod2"),
## BIOMOD+ system.file( "external/bioclim/current/bio12.grd",
## BIOMOD+ package="biomod2"))
##
## BIOMOD> # 1. Formatting Data
## BIOMOD> myBiomodData <- BIOMOD_FormatingData(resp.var = myResp,
## BIOMOD+ expl.var = myExpl,
## BIOMOD+ resp.xy = myRespXY,
## BIOMOD+ resp.name = myRespName)
##
## -=-=-=-=-=-=-=-=-=-=-=-= GuloGulo Data Formating -=-=-=-=-=-=-=-=-=-=-=-=
##
## > No pseudo absences selection !
## ! No data has been set aside for modeling evaluation
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
##
## BIOMOD> # 2. Defining Models Options using default options.
## BIOMOD> myBiomodOption <- BIOMOD_ModelingOptions()
##
## BIOMOD> # 3. Doing Modelisation
## BIOMOD>
## BIOMOD> myBiomodModelOut <- BIOMOD_Modeling( myBiomodData,
## BIOMOD+ models = c('SRE','RF'),
## BIOMOD+ models.options = myBiomodOption,
## BIOMOD+ NbRunEval=2,
## BIOMOD+ DataSplit=80,
## BIOMOD+ VarImport=0,
## BIOMOD+ models.eval.meth = c('TSS','ROC'),
## BIOMOD+ do.full.models=FALSE,
## BIOMOD+ modeling.id="test")
##
##
## Loading required library...
##
## Checking Models arguments...
##
## Creating suitable Workdir...
##
## > No weights : all observations will have the same weight
##
##
## -=-=-=-=-=-=-=-=-=-=-=-= GuloGulo Modeling Summary -=-=-=-=-=-=-=-=-=-=-=-=
##
## 5 environmental variables ( bio3 bio4 bio7 bio11 bio12 )
## Number of evaluation repetitions : 2
## Models selected : SRE RF
##
## Total number of model runs : 4
##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
##
##
## -=-=-=- Run : GuloGulo_AllData
##
##
## -=-=-=--=-=-=- GuloGulo_AllData_RUN1
##
## Model=Surface Range Envelop
## Evaluating Model stuff...
## Model=Breiman and Cutler's random forests for classification and regression
## Evaluating Model stuff...
##
## -=-=-=--=-=-=- GuloGulo_AllData_RUN2
##
## Model=Surface Range Envelop
## Evaluating Model stuff...
## Model=Breiman and Cutler's random forests for classification and regression
## Evaluating Model stuff...
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
##
## BIOMOD> ## print a summary of modeling stuff
## BIOMOD> myBiomodModelOut
##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-= BIOMOD.models.out -=-=-=-=-=-=-=-=-=-=-=-=-=-=
##
## Modeling id : test
##
## Species modeled : GuloGulo
##
## Considered variables : bio3 bio4 bio7 bio11 bio12
##
##
## Computed Models : GuloGulo_AllData_RUN1_SRE GuloGulo_AllData_RUN1_RF
## GuloGulo_AllData_RUN2_SRE GuloGulo_AllData_RUN2_RF
##
##
## Failed Models : none
##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=## list of all objects in the workspace
ls()## [1] "DataSpecies" "metadata" "myBiomodData"
## [4] "myBiomodModelOut" "myBiomodOption" "myExpl"
## [7] "myResp" "myRespName" "myRespXY"## print the modelling summary
myBiomodModelOut##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-= BIOMOD.models.out -=-=-=-=-=-=-=-=-=-=-=-=-=-=
##
## Modeling id : test
##
## Species modeled : GuloGulo
##
## Considered variables : bio3 bio4 bio7 bio11 bio12
##
##
## Computed Models : GuloGulo_AllData_RUN1_SRE GuloGulo_AllData_RUN1_RF
## GuloGulo_AllData_RUN2_SRE GuloGulo_AllData_RUN2_RF
##
##
## Failed Models : none
##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=## print the list of files saved on hard drive
list.files(myRespName, recursive = T)## [1] "GuloGulo.test.models.out"
## [2] "models/test/GuloGulo_AllData_RUN1_RF"
## [3] "models/test/GuloGulo_AllData_RUN1_SRE"
## [4] "models/test/GuloGulo_AllData_RUN2_RF"
## [5] "models/test/GuloGulo_AllData_RUN2_SRE"
## [6] "proj_current/GuloGulo.current.projection.out"
## [7] "proj_current/proj_current_GuloGulo.grd"
## [8] "proj_current/proj_current_GuloGulo.gri"We saw that a 2 types of files have been created :
xxx.models.out : which is the copy of your
BIOMOD_Modelling()output object. This is the file we will use after to do projectionsmodels/… : all single models produced. Should be useful if you want to see more accurately your models but it is not the point here.
So lets now clear the work space ( to simulate the case where you stop your modelling procedure here ).
rm(list=ls())
ls()## character(0)** IMPORTANT NOTE ** : To be able to do your projection, you will have to keep the relative path unchanged. That is to say to be in a directory where your species directory is. The easiest way to do that is to keep the same working directory for the whole analysis.
Here we have to reload the BIOMOD_Modelling() output. The file will be load with a specific name, here "GuloGulo.test.models.out". I will show you a trick to store it under the name you want.
## keep a track
(bm_out_file <- load("GuloGulo/GuloGulo.test.models.out"))## [1] "GuloGulo.test.models.out"## BIOMOD_MOdelling output is loaded with a predefined name
ls()## [1] "bm_out_file" "GuloGulo.test.models.out"## lets store this object with another name and free workspace
my_bm_model <- get(bm_out_file)
rm(list = c(bm_out_file, 'bm_out_file'))
ls()## [1] "my_bm_model"## check our model is well loaded
my_bm_model ##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-= BIOMOD.models.out -=-=-=-=-=-=-=-=-=-=-=-=-=-=
##
## Modeling id : test
##
## Species modeled : GuloGulo
##
## Considered variables : bio3 bio4 bio7 bio11 bio12
##
##
## Computed Models : GuloGulo_AllData_RUN1_SRE GuloGulo_AllData_RUN1_RF
## GuloGulo_AllData_RUN2_SRE GuloGulo_AllData_RUN2_RF
##
##
## Failed Models : none
##
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=Now we have our BIOMOD_Modelling() output correctly load, we just have to load the environmental predictors and apply BIOMOD_Projection() function.
** IMPORTANT NOTE ** : you have to ensure that your explanatory variables names match exactly with the ones you used at BIOMOD_FormatingData() step.
## here we use exacly the same data than in BIOMOD_FormatingData() step
new_env <- stack( system.file( "external/bioclim/current/bio3.grd",
package="biomod2"),
system.file( "external/bioclim/current/bio4.grd",
package="biomod2"),
system.file( "external/bioclim/current/bio7.grd",
package="biomod2"),
system.file( "external/bioclim/current/bio11.grd",
package="biomod2"),
system.file( "external/bioclim/current/bio12.grd",
package="biomod2"))
## then call Projection function
myBiomodProjection <- BIOMOD_Projection(modeling.output = my_bm_model,
new.env = new_env,
proj.name = 'current',
selected.models = 'all',
compress = FALSE,
build.clamping.mask = FALSE)##
## -=-=-=-=-=-=-=-=-=-=-=-=-= Do Models Projections -=-=-=-=-=-=-=-=-=-=-=-=-=
##
## > Projecting GuloGulo_AllData_RUN1_SRE ...
## > Projecting GuloGulo_AllData_RUN1_RF ...
## > Projecting GuloGulo_AllData_RUN2_SRE ...
## > Projecting GuloGulo_AllData_RUN2_RF ...
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=## check that it worked properly
plot(myBiomodProjection)## Loading required package: rasterVis
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
##
## Attaching package: 'latticeExtra'
##
## The following object is masked from 'package:ggplot2':
##
## layer
##
## Loading required package: hexbin
## Loading required package: grid