Amy Whitehead's Research

the ecological musings of a conservation biologist

Creating a presence-absence raster from point data

13 Comments

I’m working on generating species distribution models at the moment for a few hundred species. Which means that I’m trying to automate as many steps as possible in R to avoid having to click buttons hundreds of times in ArcView.

One of the tasks that I need to do is to convert presence-only latitude and longitude data into a presence-absence raster for each species. It seems like this would be something that relatively simple but it took me longer than it should have to figure it out. So I’m posting my code here so 1) I don’t forget how I did it; and 2) because I had someone ask me how to exactly this thing this afternoon and it took me ages to hunt through my poorly organised files to find this piece of code! So here it is:

Because I’m a function kinda girl, I wrote this as a function. It basically goes through three steps:

1. Take an existing raster of the area you are interested in mask.raster and set the background cells to zero (absences).

2. rasterize the presence points for your species species.data and set those cells to one (presences).

3. Label the new raster by your species names raster.label and save it as a new raster.

presence.absence.raster <- function (mask.raster,species.data,raster.label="") {
require(raster)

# set the background cells in the raster to 0
mask.raster[!is.na(mask.raster)] <- 0

#set the cells that contain points to 1
speciesRaster <- rasterize(species.data,mask.raster,field=1)
speciesRaster <- merge(speciesRaster,mask.raster)

#label the raster
names(speciesRaster) <- raster.label
return(speciesRaster)
}

Below is an example of how the function works using data on the global distribution of foxes data from the biomod2 package.

library(biomod2)

# read in species point data and extract data for foxes
mySpecies <- read.csv(system.file("external/species/mammals_table.csv", package="biomod2"), row.names = 1)
species <- "VulpesVulpes"

# extract fox data from larger dataset and keep only the x and y coordinates
fox.data <- mySpecies[,c("X_WGS84", "Y_WGS84",species)]
fox.data <- fox.data[fox.data$VulpesVulpes==1,c("X_WGS84", "Y_WGS84")]

# read in a raster of the world
myRaster <- raster(system.file( "external/bioclim/current/bio3.grd",package="biomod2"))

# create presence absence raster for foxes
pa.raster <- presence.absence.raster(mask.raster=myRaster, species.data=fox.data, raster.label=species)
plot(pa.raster, main=names(pa.raster))

run_function

In this plot, the presences (1) are shown in green and the absences (0) in light grey.

Helpful things to remember (or things I learnt the hard way)

  1. Make sure your species point data and raster are in the same projection and that they actually overlap!
  2. Set your desired raster extent and resolution in the mask.raster before you get started.
  3. The species point data that you feed into the function should just be a list of x and y co-ordinates – no species names or abundances or you’ll confuse the poor beast and it won’t work!

And yes, foxes are also present in Australia where they are a pest. I guess this map shows their natural range before people started doing silly things.

Bought to you by the powers of knitr & RWordpress 

(Well it was and then things went a little awry, so I had to do some tinkering!)

Advertisements

13 thoughts on “Creating a presence-absence raster from point data

  1. Pingback: Recent Qaecologist blog posts (May 2013, Weeks 3 & 4) | Quantitative & Applied Ecology Group

  2. Hi Amy,

    How do you then read presence-absence raster in BIOMOD_FormatingData? I thought resp.var had to be vector.

    Thanks

    Andre

    • Hi Andre,
      Sorry for the super-slow reply – your comment got lost in my spam filter. I’m actually using my presence-absence rasters in Zonation but was using the Biomod data as example data to demonstrate how the function worked. Sorry if that caused some confusion.

  3. Hi Amy,
    I am trying to do similar things to this. What I would like to do is plot the presence/absence of two separate species that overlap just slightly, and eliminate records in cells where both species are present. I basically want to remove presence locations from my data set where both species occur in the same raster grid cell – In other words, eliminating presence records in the “contact zone” between my two species for further analysis and distribution modeling. Any idea how I would do this based on the code you have here?

    Really great stuff!

    Mike

    • Hi Michael, I think that the easiest way to do this would be to first make separate presence-absence rasters for each of your species using the technique in the post and then sum the rasters. This would give you an output raster that would have ones where the species occur separately and twos where they occur together. You can then remove those cells that have a value of two. If you are interested in identifying which species occurs where, you may wish to change the values of the presence points for one species before you sum them together. For example:

      # create PA rasters for each species
      sp1.raster <- presence.absence.raster(mask.raster=myRaster, species.data=sp1, raster.label="species1")
      sp2.raster <- presence.absence.raster(mask.raster=myRaster, species.data=sp2, raster.label="species2")
      sp2.raster[sp2.raster==1] <- 2 #change all the presence points for species 2 to have a value of 2

      # sum the individual species rasters to identify areas where they co-occur
      sp.overlap <- sum(sp1.raster,sp2.raster)
      sp.overlap[sp.overlap==3] <- 0 #set all cells where both species occur to 0

      Hopefully that makes sense. Good luck with your analysis!

  4. Pingback: Extracting raster data using a shapefile | Amy Whitehead's Research

  5. Hi,

    this code was written a while ago, but is really useful for me, so thanks Amy!
    I ended up with an error at the end, and I cannot figure out why

    when I do the pa.raster, I have this error message:

    Error in .doCellFromXY(object@ncols, object@nrows, object@extent@xmin, :
    not compatible with requested type
    Called from: .doCellFromXY(object@ncols, object@nrows, object@extent@xmin,
    object@extent@xmax, object@extent@ymin, object@extent@ymax,
    x, y)

    my code is as follow:

    library(biomod2)
    # read in species point data and extract data for cats
    mySpecies <- read.csv("all_species.csv",header = TRUE, sep = ";", row.names = 1)
    head(mySpecies)

    looks like this:
    X Y Felis_catus Sarcophilus_harrisii Dasyurus_viverrinus
    1 496500,0019 5192499,998 1 1 0
    4 520500,0018 5192499,998 0 0 0
    7 526499,9996 5195499,997 1 0 0
    10 496500,0019 5201499,996 0 1 1
    13 526499,9996 5201499,996 0 0 0
    16 505499,9988 5204500 0 1 1
    Dasyurus_maculatus
    1 0
    4 0
    7 0
    10 0
    13 0
    16 0

    species <- "Felis_catus"
    # extract cat data from larger dataset and keep only the x and y coordinates
    cat.data <- mySpecies[,c("X", "Y",species)]
    cat.data <- cat.data[cat.data$Felis_catus==1,c("X", "Y")]
    cat.data
    # read in a raster with the extent and resolution I want
    myRaster <- raster("slope2.tif")
    # create presence absence raster for cats
    pa.raster error message there

    plot(pa.raster, main=names(pa.raster))

    Can you help me?
    I don’t know if it’s just because my excel is doing a “,” instead of a “.” for the numbers with decimals???

    • Hi Elodie,
      Glad to hear that you are finding the code helpful (even if it isn’t working!). I suspect that you are correct about the incorrect coding of the decimal places in the location data. You might like to look at using read.csv2 to read in your data – see this StackOverflow thread for more details.
      You could also try doing a test by manually making a dataframe with a few rows of your data with the correct decimal coding and see if that fixes the problem.
      Hopefully that helps but let me know if you continue to have problems and I’ll see if I can nut something out.

      • Hi thanks for the quick reply
        it indeed didn’t put the error anymore!
        BUT

        now, I have a plot of my background image with 0 everywhere (nice) but none of my presence points…

        with the 1st bit of code, it should rasterize my presence points and set them at 1 right? I don’t know what’s wrong here,
        I even created a matrix just to see, or adding coordinates and a projection to my dataframe (cats presence/absence)
        but nothing worked…

        If you have another idea to help, that would be greately appreciated 🙂

      • Hi Elodie,
        Glad we managed to fix the first problem. As for the second problem, I would double check that your points overlap with the raster. Using the example data in the post:
        coordinates(fox.data) <- ~ X_WGS84 + Y_WGS84
        plot(myRaster)
        plot(fox.data,add=T)

        If they don’t overlap, then that will explain the lack of ones in your raster. If not, then there is something else going on. Let me know how you get on.

  6. they overlapped, so finally I managed to create the raster with background points at 0 and presences at 1, but using a different method, hope it’s good!

    ### create a mask raster = raster of the background with 0 values everywhere
    mask.raster= raster (“new_slope.tif”) # a raster with the right ext and res
    mask.raster[!is.na(mask.raster)] <- 0
    mask=mask.raster
    plot(mask)
    writeRaster(mask, "mask.tif", overwrite=TRUE)
    mask = raster ("mask.tif")

    ### presences
    cats = read.table("cats_presences_absences.csv",header = TRUE, sep=";", dec=",")
    head(cats)
    class(cats)
    #data.frame
    rcats = rasterFromXYZ(cats)
    class(rcats)
    #RasterBrick
    proj4string(rcats) <- CRS("+proj=utm +zone=55 +south +ellps=GRS80 +units=m +no_defs")
    rcats
    rastercats = raster(rcats, layer=1)
    #only keep the presence/absence column
    rastercats
    new_cat = resample(rastercats, new_slope)
    new_cat
    plot(new_cat)
    writeRaster(new_cat, "new_cat.tif", overwrite=TRUE)

    # now, combine the cat distribution with the mask
    cuts=c(0,0.25,0.5,0.75,1)
    pal <- colorRampPalette(c("grey88","violetred3"))
    together=merge(new_cat, mask)
    plot(together, breaks=cuts, col=pal(6))
    writeRaster(together,"cat_map.tif", overwrite=TRUE)

    • Hi Elodie. That seems like a sensible approach and should give you the results you were after. Still puzzled by why my function didn’t work but glad you were able to find a solution. Good luck with the rest of your work.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s