The aim of this document is to show different options to produce biomod2 projections plots. We will use the example with some single models projections but the scripts and functions should be easily adapted to ensemble models case.

Let’s load biomod2 package and execute BIOMOD_Projection example.

rm(list=ls())
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-64 loaded.
## 
## Type browseVignettes(package='biomod2') to access directly biomod2 vignettes.
example(BIOMOD_Projection)
## 
## BIOMOD> # species occurrences
## BIOMOD> DataSpecies <- read.csv(system.file("external/species/mammals_table.csv",
## BIOMOD+                                     package="biomod2"), row.names = 1)
## 
## BIOMOD> head(DataSpecies)
##   X_WGS84 Y_WGS84 ConnochaetesGnou GuloGulo PantheraOnca PteropusGiganteus
## 1   -94.5      82                0        0            0                 0
## 2   -91.5      82                0        1            0                 0
## 3   -88.5      82                0        1            0                 0
## 4   -85.5      82                0        1            0                 0
## 5   -82.5      82                0        1            0                 0
## 6   -79.5      82                0        1            0                 0
##   TenrecEcaudatus VulpesVulpes
## 1               0            0
## 2               0            0
## 3               0            0
## 4               0            0
## 5               0            0
## 6               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=1, 
## BIOMOD+                                        DataSplit=70,
## BIOMOD+                                        models.eval.meth = c('TSS'),
## BIOMOD+                                        do.full.models = FALSE)
## 
## 
## 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 : 1
## Models selected : SRE RF 
## 
## Total number of model runs : 2 
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
## 
## 
## -=-=-=- 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...
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
## 
## BIOMOD> # 4.1 Projection on current environemental conditions
## BIOMOD> 
## BIOMOD> myBiomodProjection <- BIOMOD_Projection(modeling.output = myBiomodModelOut,
## BIOMOD+                                           new.env = myExpl,
## BIOMOD+                                           proj.name = 'current',
## BIOMOD+                                           selected.models = 'all',
## BIOMOD+                                           binary.meth = 'TSS',
## BIOMOD+                                           compress = FALSE,
## BIOMOD+                                           build.clamping.mask = FALSE)
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-= Do Models Projections -=-=-=-=-=-=-=-=-=-=-=-=-=
## 
##  > Projecting GuloGulo_AllData_RUN1_SRE ...
##  > Projecting GuloGulo_AllData_RUN1_RF ...
## 
##  > Building TSS binaries
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
## 
## BIOMOD> ## Not run: 
## BIOMOD> ##D # 4.2 Projection on future environemental conditions
## BIOMOD> ##D myExplFuture = stack(system.file("external/bioclim/future/bio3.grd",package="biomod2"),
## BIOMOD> ##D                      system.file("external/bioclim/future/bio4.grd",package="biomod2"),
## BIOMOD> ##D                      system.file("external/bioclim/future/bio7.grd",package="biomod2"),
## BIOMOD> ##D                      system.file("external/bioclim/future/bio11.grd",package="biomod2"),
## BIOMOD> ##D                      system.file("external/bioclim/future/bio12.grd",package="biomod2"))
## BIOMOD> ##D 
## BIOMOD> ##D myBiomodProjectionFuture <- BIOMOD_Projection(modeling.output = myBiomodModelOut,
## BIOMOD> ##D                                               new.env = myExplFuture,
## BIOMOD> ##D                                               proj.name = 'future',
## BIOMOD> ##D                                               selected.models = 'all',
## BIOMOD> ##D                                               binary.meth = 'TSS',
## BIOMOD> ##D                                               compress = FALSE,
## BIOMOD> ##D                                               build.clamping.mask = TRUE)
## BIOMOD> ##D 
## BIOMOD> ##D # print summary and plot projections
## BIOMOD> ##D myBiomodProjectionFuture
## BIOMOD> ##D plot(myBiomodProjectionFuture)
## BIOMOD> ## End(Not run)
## BIOMOD>                                        
## BIOMOD> 
## BIOMOD> 
## BIOMOD>

Lets check what we have in our work space.

ls()
## [1] "DataSpecies"        "myBiomodData"       "myBiomodModelOut"  
## [4] "myBiomodOption"     "myBiomodProjection" "myExpl"            
## [7] "myResp"             "myRespName"         "myRespXY"

The BIOMOD_Projection output object is myBiomodProjection.

myBiomodProjection
## 
## -=-=-=-=-=-=-=-=-=-=-=-= 'BIOMOD.projection.out' -=-=-=-=-=-=-=-=-=-=-=-=
## 
## Projection directory : GuloGulo/current
## 
## 
## sp.name : GuloGulo
## 
## expl.var.names : bio3 bio4 bio7 bio11 bio12
## 
## 
## modeling id : 1418059580 ( GuloGulo/GuloGulo.1418059580.models.out )
## 
## models projected : GuloGulo_AllData_RUN1_SRE, GuloGulo_AllData_RUN1_RF
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

Here only 2 models have been produced and the outputs are stored into “GuloGulo/proj_current” directory.

list.files("GuloGulo/proj_current/", recursive = T)
## [1] "GuloGulo.current.projection.out"  "proj_current_GuloGulo.grd"       
## [3] "proj_current_GuloGulo.gri"        "proj_current_GuloGulo_TSSbin.grd"
## [5] "proj_current_GuloGulo_TSSbin.gri"

Here we see a copy of myBiomodProjection object called GuloGulo.current.projection.out and 2 .grd files (with associated .gri files) corresponding to the save of the projections and binary projections under R rasterStack format (each model projection saved as a layer).

Note : If you have use the option do.stack = F in BIOMOD_Projection then each models projections (resp. binary projection) will be stored in individual_projection directory. In this case object will be rasterLayers object.

Note : If you have use the option output.format = '.img' in BIOMOD_Projection then you will have ‘.img’ outputs instead of ‘.grd’ (and ‘.gri’).

Now the projections maps are built we should be interested in plotting them. There is several way to do that. Here some (not exhaustive) are exposed.

  1. using biomod2 utilities

Plot all models projection using an unified color scale.

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

plot of chunk unnamed-chunk-5

You should also make some selection of models projections you want to plot :

## all RF projections
plot(myBiomodProjection, str.grep = "RF")

plot of chunk unnamed-chunk-6

## all projection built with "AllData" dataset
plot(myBiomodProjection, str.grep = "AllData")

plot of chunk unnamed-chunk-6

## plot a define models here "GuloGulo_AllData_RUN1_SRE"
plot(myBiomodProjection, str.grep = "GuloGulo_AllData_RUN1_SRE")

plot of chunk unnamed-chunk-6

## you can also make whatever other models assembling using grep utilities
( all_proj <- get_projected_models(myBiomodProjection) )
## [1] "GuloGulo_AllData_RUN1_SRE" "GuloGulo_AllData_RUN1_RF"
( sel_proj <- grep("RF", grep("AllData", all_proj, value = T), value=T ) ) ## RF and AllData
## [1] "GuloGulo_AllData_RUN1_RF"
plot(myBiomodProjection, str.grep = sel_proj)

plot of chunk unnamed-chunk-6

If this representations not suits you then you should load directly the rasterStack (resp rasterLayer) object in your work space. You will be then free to produce whatever graph you want.

## you should use get_predictions wrapper 
mod_proj <- get_predictions(myBiomodProjection)
mod_proj
## class       : RasterStack 
## dimensions  : 47, 120, 5640, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 3, 3  (x, y)
## extent      : -180, 180, -57.5, 83.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : GuloGulo_AllData_RUN1_SRE, GuloGulo_AllData_RUN1_RF 
## min values  :                         0,                        0 
## max values  :                      1000,                     1000
## or search directly in your hard drive for the files that interest you
mod_proj_b <- raster::stack("GuloGulo/proj_current/proj_current_GuloGulo.grd")
mod_proj_b
## class       : RasterStack 
## dimensions  : 47, 120, 5640, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 3, 3  (x, y)
## extent      : -180, 180, -57.5, 83.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : GuloGulo_AllData_RUN1_SRE, GuloGulo_AllData_RUN1_RF 
## min values  :                         0,                        0 
## max values  :                      1000,                     1000

Then you are able to use all raster, rasterVis, … packages function to make your plots.

## classical plot
plot(mod_proj)

plot of chunk unnamed-chunk-8

## plot only the 2 layer and rename it
plot(subset(mod_proj,2), main = "RF projections")

plot of chunk unnamed-chunk-8

## use some rasterVis functions
rasterVis::levelplot(mod_proj)

plot of chunk unnamed-chunk-8

## or sp
sp::spplot(mod_proj)

plot of chunk unnamed-chunk-8

Now you have the full control on what and how you plot your biomod2 projections.

This document is experimental.. any comment or amelioration is more than welcome.