I have a shapefile of range polygons. Each colored set of polygons represents a different group of species/genera/whatever.

Iniital Polygon Shapefile

I next overlay a 5x5 grid onto this map. And then use join attributes by location to find the number of grids that intersect each polygon

Join Grid and Polygon Layer

The final output gives me a column in the attribute table COUNT that indicates the number of grids intersected by the polygons.

EDIT: In other words, it gives me the site occupancy of each species.

################## NEXT PHASE

I want to duplicate this same process entirely within R. First I load the same shapefile of range polygons into R using the rgdal package.


I then generate an analogous 5x5 grid in R using the following function and the raster package.

# Generate lat/lng grid for use with mammals dataset
generateGrid<-function(XMin=-180,XMax=180,YMin=-90,YMax=90,BinSize=5) {
    Grid<-raster(extent(matrix(c(XMin,YMin,XMax,YMax),nrow=2)),nrow=LatitudeRows,ncol=LongitudeColumns,crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

# Create a GridShape object

Next I attempt to overlay them using the over() function in the sp package with the following function.

# Find Occupancy of taxa using shapefile
shapeOccupancy<-function(RangeShape,GridShape) {
    for (i in 1:nrow(RangeShape)) {
        Present<-apply(is.na(Intersect),1,all) # identify grids that do not intersect range polygons

################## The problem

The counts given by my R process are not even close to the counts given by the QGIS process. I could understand slight differences, but these are quite large.

# Sample of differences
                     RResults  QGISResult
Abrothrix                7         23     
Alouatta                31         69     
Aplodontia               1         6     
Artibeus                47         88     
Ateles                  24         48     
Avahi                    1         7     
Balantiopteryx           3         17

My best guess currently is that this is a result of the fact that some of the range polygons are discontinuous and R has a difficult time handling that, but I'm not sure. Based on a few manual spotchecks (i.e., manually counting grid and polygon intersections on the map), the QGIS results seem to be accurate and it is the R results that are wrong.

EDIT: A subset of the range polygons is included at the following github repository. https://github.com/aazaff/StackOverflowExamples/

  • I would need to work though your problem. Could you subset the data so it can be posted? I have no idea how QGIS is coded for this but one thing that comes to mind is if the pixels are being treated as cell centroids it could make a huge difference. The SpatialPixelsDataFrame class represents cells as polygons whereas SpatialGridDataFrame are cell centroid points. Because of this, one could test the differences in results based on different input vector topologies. You could also use rasterToPoints to directly test cell-centroid based representation of the raster. – Jeffrey Evans Jul 16 '15 at 16:22
  • @JeffreyEvans Hi, I added a link to a subset of the data. – Andy Jul 16 '15 at 16:44
  • @JeffreyEvans If you have problem accessing it let me know, I can put it up on one of our servers. – Andy Jul 16 '15 at 21:40

If I understand your desired outcome correctly, you would like a count (richness) of species for each grid cell in the defined raster. I cannot speak to the differences between R and QGIS but I came up with a much more optimized and faster way to conduct your analysis. I leverage the raster package and use a raster stack to accumulate species. The workflow is as follows.

Load required libraries and read the species distribution polygons.


spp <- readOGR(getwd(), "overflowSubset07162015")

Create raster with desired resolution or rows/columns. This is the reference raster that will be used to rasterize the species polygons.

xmin=-180; xmax=180; ymin=-90; ymax=90; bin=5
bins <- raster(extent(matrix(c(xmin,ymin,xmax,ymax),nrow=2)),
                 nrow=length(seq(ymin, ymax, bin)), 
                 ncol=length(seq(xmin, xmax, bin)),
                 crs= "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
bins[] <- 1:ncell(bins)

Now we can create binary rasters for each species corresponding to the reference raster. We first create an empty stack and then populate it by subsetting each species and rasterizing it to the reference raster. For clarity, we assign names of each species to each raster in the stack.

spp.pa <- stack()
  for(i in unique(spp@data$genus_name) ) {
    s <- as(spp[spp$genus_name == i,], "SpatialPolygons")
    spp.pa <- addLayer(spp.pa, rasterize(s, bins, fun = 'count', background = 0))
names(spp.pa) <- unique(spp@data$genus_name)

Now we can calculate the sum (species counts at each pixel) by using the calc function on the raster stack.

spp.sum <- calc(spp.pa, fun = sum) 

If you want to get at fractional cover of a species in a given pixel you can use the getCover argument in the rasterize function. In this way you could create a conditional argument in the for loop that if a species is less than a defined percent coverage the pixel is 0. Here is a quick example that sets pixels < 30% coverage to 0 else 1.

p = 30
spp.pct <- stack()
  for(i in unique(spp@data$genus_name) ) {
    s <- as(spp[spp$genus_name == i,], "SpatialPolygons")
    sfp <- rasterize(s, bins, getCover = TRUE, background = 0)
    sfp[sfp < p] <- 0
    sfp[sfp >= p] <- 1
    spp.pct <- addLayer(spp.pct, sfp)
names(spp.pct) <- unique(spp@data$genus_name)

spp.pct.sum <- calc(spp.pct, fun = sum)

  • I am looking for the site occupancy of each species, not the richness of each site. It's implicit in your answer with apply(spp.pa,2,sum). This gives an identical answer to the function I posted above. Unfortunately the answer still doesn't work out correctly. Some species have site occupancies of 0, which is definitely not correct. Is it possible that readOGR round small polygons down or something like that, making them disappear? – Andy Jul 17 '15 at 18:14
  • If you look at the resulting raster stack each species has an associated binary raster that indicates occupancy ie., plot(spp.pa). With a quick look at results I would recommend using the getCover argument and remove the [0,1] assignment in my example thus, directly using the resulting percents. One genus "Nyctiellus" exhibits a 1-13% pixel coverage at bin=10 and could potentially get dropped during rasterization with large pixels. It is much more stable when I drop bin to 1. You could use my last example and set p=0 and then make the assignments: sfp[sfp <= p] <- 0; sfp[sfp > p] <- 1 – Jeffrey Evans Jul 17 '15 at 18:32
  • I think I am confused about something fundamental here, sorry. I don't know how to explain the problem better. In QGIS, the Nyctiellus polygon intersects with six unique 5x5 lat/lng cells. In R, it crosses zero cells. This doesn't make sense as R gives a polygon extent of: -84.9574, -74.1318, 19.82804, 25.56344 (xmin, xmax, ymin, ymax). Which means that it should cross more than one 5x5 cell (if not exactly six). – Andy Jul 17 '15 at 18:53
  • Reproducing your subset extent, I am getting exactly 6 cells intersecting the Nyctiellus polygons with 1%, 4%, 8%, 13%, 30%, 38% polygon cell coverage. So, I am not sure what is going on. Try: plot(s); plot(as(bins, "SpatialPolygons"), add=T) to insure that there is overlap and please use my example that returns percent coverage as it seems to be more stable. – Jeffrey Evans Jul 17 '15 at 21:00
  • Okay, I think I am getting closer to understanding. Taking the instructions from your first comment, I set p = 0; sfp[sfp <= p] <- 0; sfp[sfp > p] <- 1. I then convert spp.pa<-as.data.frame(spp.pct) and sum the columns apply(spp.pa,2,sum). Nyctiellus now equals 5, which is much closer to the correct answer than before. What additional steps are you taking to get 6? – Andy Jul 20 '15 at 16:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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