9

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, ...
Improve this question
3
  • 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
    Commented Jul 4, 2017 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
    Commented Jul 4, 2017 at 16:01
  • Found the solution here: www.r-bloggers.com/extract-values-from-numerous-rasters-in-less-time/
    – Bindini
    Commented Jul 4, 2017 at 18:04

3 Answers 3

Reset to default
16

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.

Improve this answer
1
  • 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
    Commented Jul 4, 2017 at 17:34
3

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.

Improve this answer
0

Sometimes parallel is not the way to go. If the function is not written for it you get all the overheads and none of the benefits.

Often, just reducing the demand for memory allocation solves performance issues.

I have found that if the number of points is high it is much quicker to hit the raster layers one at a time (or in smaller stacks) rather than all at once.

It is relatively trivial to loop (yes, performance benefits from for loops in R!) through a list of raster filenames, extract, and cbind the results as you go. You can go the other way and subset the points and hit the complete stack and rbind the results similarly.

Further, if you test for multiple points hitting the same pixel there is potentially significant gain there as well (especially as your raster objects are quite small (I often deal with rasters in the 30k*20k and throw more points at them than you have in your rasters)

Improve this answer

Your Answer

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

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