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?