0

I want to use rasterize to get how much of a raster cell is covered by polygons. This can be done like this:

library(raster)

p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
pols <- SpatialPolygons( list(  Polygons(list(Polygon(p1), Polygon(hole)), 1), Polygons(list(Polygon(p2)), 2), Polygons(list(Polygon(p3)), 3)))
pols@polygons[[1]]@Polygons[[2]]@hole <- TRUE

r <- raster(ncol=180, nrow=90)

r <- rasterize(pols, r, getCover=T)

plot(r)

But is it possible to make getCover account for overlaps, so that the cover can exceed 100 if there are more polygons overlapping?

1

Do a rasterize for each polygon and sum the stack:

> st = stack(lapply(1:length(pols), function(i){rasterize(pols[i], r, getCover=TRUE)}))
> plot(st)
> plot(sum(st))

enter image description here

  • I was thinking in the same line. Its just... I have 400.000 polygons. But might still be the way to go. – Jeppe Olsen Mar 19 at 12:38
  • With 400,000 you might want to do one at a time and add as you go using maybe Reduce or a for loop instead of an apply loop that makes a stack of 400,000 layers! – Spacedman Mar 19 at 12:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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