Raster intersect and gIntersection give different results in R

Question

I'm trying to clip a polygon shapefile (species distribution) with another polygon grid. The goal is to calculate the area occupied by the species within each grid cell. I was able to do it with gIntersection from 'rgeos' package, however this function does not retain the arributes from the original shapefiles. I've also tried it using intersect from 'raster' package, and although it does keep all the attributes there is an obvious problem wiht the resulting shapefile.

Here is my code.

# libraries used
library(rgdal)
library(raster)
library(rgeos)

# Read shapefiles    
mdau <- readOGR(shp.dir, 'myo_dau') # species disribution data
grid <- readOGR(grid.dir, 'grid300km') # grid

plot(mdau, axes=T); plot(grid, add=T) # check shapefiles are correctly ploted

# Clip using intersect from 'raster' package
inter.raster <- intersect(mdau, grid)
# here I get warnings, which I suspect may be the cause of the unexpected result
# Example Warning messages:
# 1: In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
  Ring Self-intersection at or near point -354560.80502695002 8127597.0226868102

inter.raster$area <- area(inter.raster) # area calculation

# Clip using gIntersection from 'rgeos' package
inter.rgeos <- gIntersection(mdau, grid, byid = TRUE)
inter.rgeos$area <- area(inter.rgeos) # area calculation

# Plot results
par(mfrow=c(1,2))
plot(inter.raster, axes=T, main='Raster::Intersect'); plot(inter.rgeos, axes=T, main='Rgeos::gIntersection')

This is the compared output

I would be happy with gIntersection output, however it does not keep attributes from the original shapefiles.

head(inter.rgeos@data)
        area
1   82349136
2   82349136
3   41174568
4   82349136
5 1297047409
6  150000000

Any idea or suggestion of what might be going on?

Answer 1

Well, I am not sure why raster::intersect is failing, probably an error in the feature geometry of your data. In the sp instance, intersect is basically a wrapper function for rgeos::gIntersection so, from a topological perspective, it should be replicating the results from gIntersection. The trick to using the gIntersection function is the byid=TRUE argument, which defaults to FALSE.

Let's approximate the raster::intersect function. First we need some data and the required libraries.

library(sp)
library(spdep)
library(rgdal)
library(raster)
library(rgeos)

polys <- readOGR(system.file("shapes/columbus.shp", package="spData")[1])
r <- raster(extent(polys))
  r[] <- 1:ncell(r)
  r <- as(r, "SpatialPolygonsDataFrame")  

plot(r)
  plot(polys, add=TRUE) 

Now we can apply the gIntersection function with byid=TRUE to return the row.names from each data. The byid argument controls the return of the associated feature row.names.

m <- rgeos::gIntersection(polys, r, byid=TRUE)
  row.names(m)

We coerce the row.names into a usable vector, create new row.names and an associated data.frame that will match the merged data.

ids <- do.call(rbind, strsplit(row.names(m), ' '))
row.ids <- 1:length(ids[,1])
d <- data.frame(row.names=row.ids)  

Using the row.names, we can now create an index to match the data.frames from each corresponding feature.

idx.x <- match(ids[,1], rownames(polys@data))
idx.y <- match(ids[,2], rownames(r@data))

Then we coerce the intersected object into a SpatialPolygonsDataFrame object containing the attributes from each feature.

m <- SpatialPolygonsDataFrame(m, as.data.frame(cbind(polys@data[idx.x, ], r@data[idx.y, ])),
                              match.ID = FALSE )
head(m@data)                              
plot(m)   

At this point, if things are still going sideways on your data, it is time to look at the topology of the data.