Creating MaxEnt Bias Rasters in R

Photo by Arturo Rivera on Unsplash

Intro
I was recently asked to help a colleague with an issue she was having when trying to add a binary bias file to a MaxEnt model of mango distribution (see photo). Though I thought this would be an easy process, I quickly found out that MaxEnt is very picky about the layers being used and the error messages offer almost no insight into the nature of the problem. After a lot of headaches, we finally managed to come up with a working procedure. This document will walk you through how we did it. Hopefully, this will prove helpful to someone in the future who is struggling to get their bias files working.


Required Packages

require(tidyverse)
require(raster)
require(sp)
require(rgdal)
require(spatstat)
require(here)
require(tools)

Environmental Layers
First we have to begin by reading in our environmental layers as raster objects. Here I am using two environmental layers for demonstrations purposes. We will also being using WGS84 for the sake of global applicability.

b1 <- raster(here(‘Data’, ‘bio_01_30s.asc’))
crs(b1) <-CRS(‘+init=EPSG:4326’)
b2 <- raster(here(‘Data’, ‘bio_02_30s.asc’))
crs(b2) <- CRS(‘+init=EPSG:4326’)

Check your work
We also want to check that our two rasters have the same CRS and extent (they align).

compareRaster(b1, b2)

So long as that stays ‘TRUE’ we can keep moving forward with the creation of our bias raster.


Creating a Blank Raster
We need to use one of our environmental layers as a template to create a blank raster with the same extent as the environmental layers. Later, we will use this blank raster to keep track of which cells do and do not contain occurrences.

blank <- raster(b1)
crs(blank) <- CRS(‘+init=EPSG:4326’)

Species Occurrence Points
For MaxEnt to work we also need a set of occurrence points. Here we read them in as a csv and turn that table into a SpatialPointsDataFrame. Again, you will want to set or transform your projections as necessary to get everything into WGS84.

oc <- read.csv(here(‘Data’, ‘occurrences.csv’))
oc <- SpatialPoints(data.frame(lon = oc$Longitude, lat = oc$Latitude))
crs(oc) <- CRS(‘+init=EPSG:4326’)

Counting Occurrences
Now we have all of the necessary information to begin making a binary bias raster. All we need to do is figure out which cells of our raster contain at least one point. NA cells are then assigned a very small value. This is really just a placeholder, since we then use our first environmental layer to mask away any water and ensure that all layers use the same NA value. If the NA values are different it will cause a MaxEnt error.

# Get count of occurrence points in each cell
oc.ras.bin <- rasterize(oc, blank, fun = ‘count’, na.rm = F)
# All cells with count >= 1 get assigned a 1 value
oc.ras.bin[oc.ras.bin >= 1] <- 1
# Redundancy to enforce consistent NA values
oc.ras.bin[is.na(oc.ras.bin)] <- 0.00001
# Mask away any water/missing data
oc.ras.bin <- mask(oc.ras.bin, b1)
# Match bias NA values with env. layers
NAvalue(oc.ras.bin) <- NAvalue(b1)

Check your work (again)
Let’s make sure we haven’t messed up any of our rasters before we move forward

compareRaster(oc.ras.bin, b1)

Creating a Spatial Mask
Okay, this is where things begin to get complicated. We need to create a new Polygon object which has corner coordinates slightly larger than the bounding box of our occurrence points. Since we are using WGS84, the units of the expansion factor are degrees. I decided to expand everything by one degree around the bounding box of the occurrence points (this is essentially a buffer). Next, we create a Polygons class object which is in turn converted into a SpatialPolygons class object by assigning it some spatial information. This final SpatialPolygons object is the ‘cookie cutter’ we will use to force our output ASCII files to have exactly identical spatial information. What is actually happening is that we are making a single template for all of our environmental layers and the bias raster to follow. Expanding this template beyond the minimum bounding box of the occurrence points ensures that we don’t sacrifice any information about the spatial neighborhood of points on the edges of the minimum bounding box (of which there will always be at least 4).

poly_expand <- function(ocr, exp_factor){
 return(Polygon(rbind(c(matrix(bbox(ocr), 2, 2)[1,1] — exp_factor,
 matrix(bbox(ocr), 2, 2)[2,2] + exp_factor),
 c(matrix(bbox(ocr), 2, 2)[1,1] — exp_factor,
 matrix(bbox(ocr), 2, 2)[2,1] — exp_factor),
 c(matrix(bbox(ocr), 2, 2)[1,2] + exp_factor,
 matrix(bbox(ocr), 2, 2)[2,1] — exp_factor),
 c(matrix(bbox(ocr), 2, 2)[1,2] + exp_factor,
 matrix(bbox(ocr), 2, 2)[2,2] + exp_factor),
 c(matrix(bbox(ocr), 2, 2)[1,1] — exp_factor,
 matrix(bbox(ocr), 2, 2)[2,2] + exp_factor)))
 )
}
p <- poly_expand(ocr = oc, exp_factor = 1)
px <- Polygons(list(p), 1)
spat.mask <- SpatialPolygons(list(px),
 proj4string = CRS(“+init=EPSG:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0”))

Apply the Spatial Mask
We are finally ready to put together all of the things we have built up so far. We extract new versions of each of our environmental layers and the binary raster using the newly created spatial mask. Then we check our alignment, and write the new files out to disk, and that’s it! You can go ahead and use these new files in MaxEnt without any of those pesky errors that crop up when you try to add bias files.

# Extract new rasters by spatial mask
b1_mod <- crop(b1, spat.mask)
b2_mod <- crop(b2, spat.mask)
bias.bin <- crop(oc.ras.bin, spat.mask)
# Double-check alignment between environmental layers and bias raster
compareRaster(b1_mod, bias.bin)
# Write new environemntal raster information with modified names
writeRaster(b1_mod, “b1_modified.asc”, overwrite = T)
writeRaster(b2_mod, “b2_modified.asc”, overwrite = T)
# Export bias layer
writeRaster(bias.bin, “bias_mask_bin_modified.asc”, overwrite = T)

Batch Processing Multiple Environmental Layers
Now that we have all of the elements working, let’s write a function we can use to convert a bunch of environmental layer files at once. In the function below I append the tag ‘_modified’ to each of the file names, but that is just so that I can better keep track of them and not risk overwriting the original files.

reduce.rasters <- function(bio.clim.file){
 bc <- raster(here(‘bioclim’, bio.clim.file))
 out <- raster::crop(bc, spat.mask)
 writeRaster(out, paste(file_path_sans_ext(bio.clim.file), ‘modified.asc’, sep = ‘_’), overwrite = T)
}
in.files <- list.files(here(‘bioclim’), pattern = ‘bio_.*\\.asc$’)
lapply(in.files, reduce.rasters)

Acknowledgements
Special thanks to Dr. Jessica McCarty and Tony Francis for their help in developing this solution.