Back in 2013 I published a post with a little function for creating presence-absence rasters from point data. I was in the midst of creating ~1000 species distribution models for Australian species and I needed something that I could easily automate in R.
Both R and my approach to working with spatial data have changed a lot in the last nine years, particularly with excellent new packages such as
fasterize. This post is a brief update to a new option for creating presence-absence rasters, prompted by a reader’s question.
The steps are essentially the same but I’ve adopted some new packages (
fasterize) that make life easier:
- Create an sf point data set from a dataframe of point locations using
- Convert the points to polygons using
st_buffer. This step is necessary as
fasterizecurrently doesn’t work with points. For this simple example, I’ve elected to leave my points in WGS84, which means I get a warning that st_buffer doesn’t work properly with lat/long data. I could have used
st_transformto convert it to a different projection before using
- Convert the polygons to a raster1 using
maskout the ocean.
Step 1: Load packages and get the point data
require(biomod2) # for the spatial data require(dplyr) # for tidy data wrangling require(sf) # for working with spatial data require(raster) # for masking rasters require(fasterize) # for creating the raster # set select to come from dplyr rather than raster select <- dplyr::select # read in a raster of the world world <- raster(system.file( "external/bioclim/current/bio3.grd",package="biomod2")) # read in species point data and extract data for foxes mammals <- read.csv(system.file("external/species/mammals_table.csv", package="biomod2"), row.names = 1) head(mammals) X_WGS84 Y_WGS84 ConnochaetesGnou GuloGulo PantheraOnca PteropusGiganteus 1 -94.5 82.00001 0 0 0 0 2 -91.5 82.00001 0 1 0 0 3 -88.5 82.00001 0 1 0 0 4 -85.5 82.00001 0 1 0 0 5 -82.5 82.00001 0 1 0 0 6 -79.5 82.00001 0 1 0 0 TenrecEcaudatus VulpesVulpes 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0
Step 2: Convert the point data to polygons
# extract fox data from larger dataset and keep only the x and y coordinates fox_data <- mammals %>% # keep only the spatial data and foxes select(X_WGS84,Y_WGS84,VulpesVulpes) %>% # keep only the presence points for foxes filter(VulpesVulpes %in% 1) %>% # convert to an sf points object projected in WGS84 st_as_sf(coords=c("X_WGS84","Y_WGS84"), crs = 4326) %>% # convert to an sf polygon object by buffering the points st_buffer(1) Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle = endCapStyle, : st_buffer does not correctly buffer longitude/latitude data head(fox_data) Simple feature collection with 6 features and 1 field Geometry type: POLYGON Dimension: XY Bounding box: xmin: -95.5 ymin: 72.00001 xmax: -75.5 ymax: 74.00001 Geodetic CRS: WGS 84 VulpesVulpes geometry 1 1 POLYGON ((-93.5 73.00001, -... 2 1 POLYGON ((-87.5 73.00001, -... 3 1 POLYGON ((-84.5 73.00001, -... 4 1 POLYGON ((-81.5 73.00001, -... 5 1 POLYGON ((-78.5 73.00001, -... 6 1 POLYGON ((-75.5 73.00001, -...
Step 3. Create the presence-absence raster
fox_raster <- fasterize(fox_data, world, # select the column to use in the raster field = "VulpesVulpes", # set the value of the background points background = 0) %>% # only retain values in non-NA cells in the original raster mask(world)
In this plot, the presences (1) are shown in green and the absences (0) in light grey. Now I could have kept the absences in the original dataset and created the raster without using the mask but I wanted to demonstrate what you would need to do if your dataset doesn’t include absence data (which seems to be the most typical scenario).
Pro tip: don’t forget to check that your species point data and raster are in the same projection and that they actually overlap before you get started.
Hopefully this is helpful. But, as always, there are many ways to do things in R and this may not be the most efficient method. I’d be interested to know how else folks are doing this sort of thing.
1 Note that you need to use a single layer raster for this approach to work as expected. If you have a raster brick or raster stack, then you can either calculate some summary statistics using the
fun argument or you can pull out a single layer of the raster brick/stack by using
my_raster[[layer]] where layer could be a name or index value.