Grid vs. Polygon Overlay in R gives different results than in QGIS

Question

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

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

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.

RangeShape<-readOGR(".","overflowExample07162015")

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) {
    LatitudeRows<-length(seq(YMin,YMax,BinSize))
    LongitudeColumns<-length(seq(XMin,XMax,BinSize))
    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")
    Grid[]<-1:ncell(Grid)
    Grid<-as(Grid,"SpatialPixelsDataFrame")
    return(Grid)
    }

# Create a GridShape object
GridShape<-generateGrid(XMin=-180,XMax=180,YMin=-90,YMax=90,BinSize=5)

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) {
    Occupancies<-vector("numeric",length=nrow(RangeShape))
    for (i in 1:nrow(RangeShape)) {
        Intersect<-over(GridShape,RangeShape[i,])
        Present<-apply(is.na(Intersect),1,all) # identify grids that do not intersect range polygons
        Occupancies[i]<-nrow(Intersect[Present!=TRUE,])
        }
    return(Occupancies)
    }

 Occupancy<-shapeOccupancy(RangeShape,GridShape)

################## 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/

Answer 1

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.

library(raster)
library(sp)
library(rgdal)

setwd("D:/TMP/spp")
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) 
plot(spp.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)

par(mfrow=c(2,1))
  plot(spp.sum)
  plot(spp.pct.sum)