biomod2 not plotting data

Question

I am trying to run an elementary species distribution model, and biomod2 is not plotting my presence data. I can attempt to recreate the problem here, but obviously I am unable to upload the tif file I am working with. To test it, I simply created a matrix of longitude/latitude values corresponding to a straight vertical line in Colombia, as follows:

tmp <- seq(-1,10,0.01) # latitude
tmp2 <- rep(-75,length(tmp)) # longitude
testVec <- cbind(tmp2,tmp)
myRespXY <- testVec

These are the coordinates of my presence data. I then only use presence data, so I have

myResp <- as.numeric(rep(1,nrow(myRespXY)))

Then, using the formatting function for biomod2,

myBiomodData <- BIOMOD_FormatingData(resp.var = myResp, # presence data (all 1's)
                                 expl.var = myExpl, # RasterStack
                                 resp.xy = myRespXY, # coordinates of presences, corresponding to a vertical line in Colombia
                                 resp.name = myRespName) # name of species
plot(myBiomodData)

What am I doing wrong? biomod2 is reading the presences, and sometimes it will plot a few points, but it doesn't plot most of them. How can I get the output to plot all of the presence data?

Answer 1

Well this is not a bug of the package but more a graphical issue that comes when you are working with high resolution raster of environmental variables.

By default the plotting function colorate pixels according to there class (presence/absence/undefined). Because pixels are so little here (and the presences/absences are so few compare to undefined area) we can't really see different colors in the plot.

A work around is to extract presences/absences from the BIOMOD_FormatingData() output use points function from raster package to display them on top of the map.

Here is a simple reproducible example:

##' This script is designed to show how you can simply extract presences, absences and pseudo absences
##' from a `data.formatted.Biomod.object` object (output of `BIOMOD_FormatingData()`).
##' 
##' This can be particularly useful to produce custom plots when high resolution rasters are used within `biomod2` framework.
##' 
##' 
##' ### Data formatting
##' 
##' Here we sill create a toy `data.formatted.Biomod.object` object using `BIOMOD_FormatingData()`
##' function.
##' 
##' Here are couple parameters used for the data formatting procedure:
##' 
##'   - *species* : Gulo Gulo
##'   - *explanatory variables*: [worlclim bioclimatic variables](http://www.worldclim.org/bioclim) 
##'   (bio_3, bio_4, bio_7, bio_11 & bio_12)
##'   - *numbeer of absences*: 0
##'   - *number of pseudo absences*: 500
##'   - *number of pseudo absences sampling repetiton*: 2
##'   - *pseudo absences selection algorithm*: `random`
##'   

library(biomod2)

# species occurrences
DataSpecies <- read.csv(system.file("external/species/mammals_table.csv",
                                    package="biomod2"), row.names = 1)
head(DataSpecies)

# the name of studied species
myRespName <- 'GuloGulo'

# the presence/absences data for our species
myRespDF <- DataSpecies[DataSpecies[, myRespName] == 1 ,c("X_WGS84","Y_WGS84", myRespName)]

myResp <- as.numeric(myRespDF[,myRespName])

# the XY coordinates of species data
myRespXY <- myRespDF[,c("X_WGS84","Y_WGS84")]


# Environmental variables extracted from Worldclim (bio_3, bio_4, bio_7, bio_11 & bio_12)
myExpl = raster::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"))

# Formatting data in biomod2 friendly sahpe
myBiomodData <- BIOMOD_FormatingData(resp.var = myResp,
                                     expl.var = myExpl,
                                     resp.xy = myRespXY,
                                     resp.name = myRespName,
                                     PA.nb.rep = 2,
                                     PA.nb.absences = 500)

##'
##' ## The extraction funcrtions
##'
##' Here are define 2 small function that will help you to extract the points locations of a 
##' `data.formatted.Biomod.object`
##' 
##' **note**: This function will be integreted in a next release of `biomod2`
##'  

library(dplyr)

## function to get PA dataset
get_PAtab <- function(bfd){
  dplyr::bind_cols(
    x = bfd@coord[, 1],
    y = bfd@coord[, 2],
    status = bfd@data.species,
    bfd@PA
  )
}

## function to get background mask
get_mask <- function(bfd){
  bfd@data.mask
}

##'
##' ### Extract presences, absences and pseudo-absences tables
##'

##'
##' When you have extracted the PA table from `data.formatted.Biomod.object` you can easily select 
##' presences (`filter(status == 1)`), pseudo-absences (`filter(is.na(status))`) or absences 
##' (`filter(status == 0)`) even if no absences are dfined in our case
##'

## get the coordiantes of presences
(pres.xy <- get_PAtab(myBiomodData) %>% 
    filter(status == 1) %>%
    select(x, y))

## get the coordiantes of absences ==> No absences here
(abs.xy <- get_PAtab(myBiomodData) %>% 
    filter(status == 0) %>%
    select(x, y))

## get the coordiantes of pseudo - absences
## all repetition of pseudo absences sampling merged 
(pa.all.xy <- get_PAtab(myBiomodData) %>% 
    filter(is.na(status)) %>%
    select(x, y)) %>%
    distinct()
## pseudo absences sampling for the first repetition only 
(pa.1.xy <- get_PAtab(myBiomodData) %>% 
    filter(is.na(status) & PA1 == TRUE) %>%
    select(x, y)) %>%
    distinct()

##'
##' ### Plot example
##'

##' 
##' Here is an example of plot that could be usefull when resolution is too high and presences 
##' could not be diplayed correctly using the `data.formatted.Biomod.object` `plot()` function
##' 


## plot the first PA selection and add the presences on top
plot(get_mask(myBiomodData)[['PA1']])
points(pres.xy, pch = 11)

This should give you a plot like this one: