7

I´m trying to parallelize a simple "extract" operation like this:

UseCores<-detectCores() -1
cl<- makeCluster(UseCores)
registerDoParallel(cl)
all.test.poly_period.test<-extract(stack.ts, test.poly)
stopCluster(cl)

stack.ts= stack of some hundred of images

test.poly= a shapefile with some polygons

But when I go to the task manager I can see that only one core is working.

I was wondering if anyone could give me a short advice about it?

The show() of the stack:

class       : RasterStack 
dimensions  : 1596, 436, 695856, 186  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent      : 284085, 297165, -2458245, -2410365  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=23 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names       : stack_MGM.1, stack_MGM.2, stack_MGM.3, stack_MGM.4, stack_MGM.5, stack_MGM.6, stack_MGM.7, stack_MGM.8, stack_MGM.9, stack_MGM.10, stack_MGM.11, stack_MGM.12, stack_MGM.13, stack_MGM.14, stack_MGM.15, ... 
min values  : -0.26170975, -0.26808712, -0.27280286, -0.27959326, -0.28714523, -0.25714213, -0.22633411, -0.21442896, -0.20916982,  -0.20526922,  -0.17409700,  -0.16982903,  -0.17428745,  -0.18175986,  -0.18332276, ... 
max values  :   0.9389138,   1.0646172,   0.8939760,   0.8914277,   0.8828803,   0.8622602,   0.8733163,   0.8952460,   0.9045494,    0.9184299,    1.0648371,    0.8974358,    0.8855312,    0.8626966,    0.8333160, ... 
  • It's already parallelized in one sense, it only builds the cell-index-to-polygon map once, then applies it. I'd suspect that the individual load for each layer in stack is a bit inefficient. Please describe how the stack is created, what the file format is, and the show the raster summary of one of them. Extract with cell numbers will still be slow if the file/s are natively tiled, for example, so you can readAll on individual layers and use cellFromPolygons (but iterating the list is a hassle). I have a package to make this easier but it's still dependent on the per layer load efficiency. – mdsumner Jul 4 '17 at 14:09
  • That´s the problem. It´s incredible slow.. taking 84 hours for each stack. The stack actually is an image time series of vegetation index. It comes from a fitting process, separately. It´s in ".tif" format. I´ll edit the question with the the show of one of them. But basically it´s: dimensions : 1596, 436, 695856, 186 (nrow, ncol, ncell, nlayers). The problem is that I need the number of the cell for each pixel time series, in order to assess them for each pixel label and relate it to other table. Thanks a lot! – Bindini Jul 4 '17 at 16:01
  • Found the solution here: www.r-bloggers.com/extract-values-from-numerous-rasters-in-less-time/ – Bindini Jul 4 '17 at 18:04
14

You need to use the beginCluster and endCluster functions of the raster package. See the example below.

library(raster)
library(snow)

# Make test data
# RasterStack
r <- raster(ncol=36, nrow=18)
r[] <- 1:ncell(r)
s <- stack(r, sqrt(r), r/r)

# SpatialPolygons
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- spPolygons(cds1, cds2)

# Visualize
plot(s, 1); plot(polys, add = TRUE)

# Extract
beginCluster(n=2)
extract(s, polys)
endCluster()

However, most of the processing time is probably spent on rasterizing the polygons, and that part is not parallelized and known to be quite inefficient. There are alternative packages to speed up that step. See velox and fasterize.

  • Velox package is amazing! Thanks a lot! Unfortunately it seems it doesn´t work for extract the time series of each pixel, right? Only the centroid pixel and the result of a function. But with the beginCluster() I was able to use more processors, so then although is still a little slow, it´s already much faster. Tks! – Bindini Jul 4 '17 at 17:34
2

Given that it's .tif we now need to know if it's tiled. It probably is, and raster extract is very slow in this situation (and is effectively in an un-maintained state with no known prospect of improvement).

I would lapply(filenames, function(x) extract(readAll(raster(x)), ts.poly)) - but that's still going to do the geometry look up every layer, so it's best to flatten the cell-index to one column with a grouping for each polygon. That's what this does:

https://github.com/hypertidy/tabularaster

I'm still guessing because we can't reproduce your situation, and it's very open-ended, but untested code I'd try is

library(tabularaster)  ## devtools::install_github("hypertidy/tabularaster")

 ## mapping between polygon `object_` and `cell_` number as per
 ## extract(..., cellnumbers = TRUE) / cellFrom* / extract(x, cells)
cells <- cellnumbers(stack.ts[[1]], ts.poly) 
exvalues <- lapply(seq_len(nlayers(stack.ts)), function(i) extract(readAll(stack.ts[[i]]), cells$cell_))

Then exvalues can be do.call(cbind, exvalues) into the matrix you would have otherwise gotten the high-level way.

I wouldn't normally write the code that way, I'd loop over the file names probably but it's too open-ended to cover all possibilities. I'm sorry not to explain everything in detail, this is sadly a topic not oft discussed and so the tools are capable, just not well understood and have a bunch of problems.

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.