I have data for species distribution modelling of two different species. I need to plot the predicted distributions of the two species in the same raster plot to be able to find the areas where they occur in allopatry(alone) and sympatry(together).

How do I plot such a map in R?

My data is a point data used with other environmental for the species distribution model. I have done the distribution modelling and have the predicted raster layer for the two species in two different raster plots. But what I want is to plot the two rasters in the same plot, such that I could show areas where the two species occur together(sympatry), where they occur differently (allopatry) and where they do not occur at all.

  • 1
    what is the type your data? movement data? or only spatial data? Oct 8, 2015 at 14:02
  • 2
    You could classify your estimated distribution rasters as [0,1] and then add the rasters together. Areas of sympatry will have a value of two, allopatry will be one and unoccupied, zero. Then you do not have to muck aground with things like raster transparencies to see where the overlap is occurring. Oct 8, 2015 at 14:51
  • Thank you immensely Jeffrey Evans, as you perfectly understood what I was asking for. I am very grateful. Please, could you let me know which functions or codes or package(s) I could use to do exactly what you have mentioned.
    – Emmanuel
    Oct 8, 2015 at 18:54
  • Welcome to GIS.se! Click 'edit' and update your original question with some more detail, including the type of file you're using (raster or point data?) and some example code of what you've tried already; the question can be re-opened.
    – Simbamangu
    Oct 9, 2015 at 4:44
  • If you have two rasters based on the same grid then you can just do arithmetic on them: both = (r1 + r2)==2 will create a raster where the sum of the two rasters r1 and r2 is 2. Is that the sort of things you want to do?
    – Spacedman
    Oct 11, 2015 at 11:27

1 Answer 1


Here is a dummy example of what you are after.

# Create some dummy data

p1 <- SpatialPolygons(list(Polygons(list(Polygon(rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20)))), 1)))
p2 <- SpatialPolygons(list(Polygons(list(Polygon(rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)))), 2)))
p1 <- SpatialPolygonsDataFrame(p1, data=data.frame(val=1), match.ID=F)
p2 <- SpatialPolygonsDataFrame(p2, data=data.frame(val=2), match.ID=F)

r <- raster(xmn=min(c(extent(p1)[1],extent(p2)[1])), 

rp1 <- rasterize(p1, r)
rp2 <- rasterize(p2, r)

# Add cell values and reclassify rasters to [0,1] to emulate species distribution
rp1[!is.na(rp1)] <- 1:length(rp1[!is.na(rp1)])
rp1[rp1 < cellStats(rp1, "mean")] <- 0 
rp1[rp1 > 0] <- 1

rp2[!is.na(rp2)] <- 1:length(rp2[!is.na(rp2)])
rp2[rp2 < cellStats(rp2, "mean")] <- 0 
rp2[rp2 > 0] <- 1

  plot(rp2, add=TRUE)  

# Intersect by summing binary rasters 
i <- rp1 + rp2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.