2

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

Species distribution and grid

# 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

Plot from both 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?

3
  • The 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.
    – Spacedman
    Mar 23, 2018 at 16:34
  • @Spacedman any idea why the functions are retreiving different outputs? I believe intersect was actually based on gIntersection, am I correct?
    – FAmorim
    Mar 23, 2018 at 16:39
  • By the way, your suggestion did the trick! Thank you @Spacedman
    – FAmorim
    Mar 23, 2018 at 16:45

1 Answer 1

3

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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