1

I'm trying to rasterize a shapefile in parallel using the solution from this post: Processing vector to raster faster with R

However, this fails to give the correct answer when the GetCover = T in the rasterize`function.

Has anyone successfully rasterized in parallel with the GetCover option enabled?

closed as off-topic by PolyGeo Oct 24 '18 at 22:58

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "When seeking help to debug/write/improve code always provide the desired behavior, a specific problem/error and the shortest code (as formatted text, not pictures) needed to reproduce it in the question body. Providing a clear problem statement and a code attempt helps others to help you." – PolyGeo
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 1
    How is the answer not correct when you use getCover=T? Is the raster nonsense? Wrong in places? Upside down? Please provide a full reproducible example, even if its mostly a copy of the one in the linked Q - that would really help. – Spacedman Oct 24 '18 at 21:31
1

It is not clear what is not working, but this works for me:

library(raster)
p1 <- rbind(c(-180,-20), c(-140,55), c(-10, 0), c(-140,-60), c(-180,-20))
p2 <- rbind(c(0,0), c(140,60), c(160,0), c(140,-55), c(0,0))
pols <- spPolygons(p1, p2)
r <- raster(ncol=90, nrow=45)

library('parallel')
cl <- makeCluster(2)
z <- clusterEvalQ(cl, library("raster"))
clusterExport(cl, c("r", "pols"))
clusterExport(cl, "r")
p <- parLapply(cl = cl, X = 1:length(pols), fun = function(i) rasterize(pols[i, ], r, getCover=TRUE))

x <- sum(stack(p))
x <- reclassify(x, cbind(1,Inf,1))
plot(x)

I would just do

y <- rasterize(pols, r, getCover=TRUE)

In the current version of raster that is pretty fast. And it also works if you have overlapping polygons, while the parallelized version would be close but not exactely the same.

  • thank you @RobertHijmams. You are correct. The problem was not the rasterize part but the combination of tiles afterward. I was erroneously doing do.call(merge, p) which takes only the first value when rasters overlap. Thanks again for your help. – JuanMayorga Oct 25 '18 at 22:52

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