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?
raster
intersection looks fine with some simple test data. I suspect your species data is wrongly coded. Try the old trick of making a 0-width buffer:mdau = gBuffer(mdau,width=0)
and then doing the raster intersection - this can sometimes clean up the geometry.