[R-sig-Geo] enable parallel analysis on pixels from rasterbrick

[Loïc Dutrieux]

Hi Jackson,

The strategy is to make your function work for raster::calc, and 
parallelize it with raster's built in parallelization utils (see ?clusterR).
Looping over pixels of a raster/matrix is nearly never the correct answer.

stars::st_apply would probably also be a good option.


library(raster)
library(trend)


r <- raster(ncol=96, nrow=63, vals=runif(63*96))
n <- 672  # number of copies
s <- stack(lapply(1:n, function(i) rz<-stack(r[i]<-(r^2/3))))

# Define the function to apply to each pixel; it should return a numeric 
or a vector of numerics with always the same lenght
fun <- function(x){
     out <- pettitt.test(x)[3]$estimate
     return(unname(out))
}

# Check that it works
fun(seq(200))

# Check that it works with calc
s_out <- calc(s, fun = fun)
s_out

# Parallelize with clusterR
beginCluster()
s_out <- clusterR(s, fun=calc, args=list(fun=fun))
endCluster()


Cheers,
Loïc


On 05/06/2020 05:16 PM, Jackson Rodrigues wrote:
> Dear all,
> 
> Could someone help me on performing parallel analysis on pixels from a
> rasterbrick?
> 
> I would like to extract and perform  Change-Point Detection from package
> trend on pixels from a rasterbrick. and write result as a netcdf file.
> However, when applied to my original data it is really time demanding
> (days).
> 
> So I would like some support on enabling parallel analysis, activating
> other cores from my computer aiming to save some time.
> I have already tried several techniques and got no satisfactory results, so
> far.
> 
> I am going through packages doParallel and foreach, but I think that my
> error is on selecting pixels from all layers as time series rather than
> from a single layer.
> 
> Below there is a copy of my function and some syntetic data created.
> 
> Thank you all in advance,
> 
> Jackson
> 
> ############
> library(raster)
> r <- raster(ncol=96, nrow=63, vals=runif(63*96))
> n <- 672  # number of copies
> s <- stack(lapply(1:n, function(i) rz<-stack(r[i]<-(r^2/3))))
> 
> Pettitt <- function(x, date, adjust = "none") {
>    require("raster")
>    require("trend")
>    require("zoo")
> 
> #Define how many cores you want to use
> UseCores <- detectCores() -1
> #Register CoreCluster
> cl<- makeCluster(UseCores)
> registerDoParallel(cl)
> 
>      # pre allocate memory
> ptt <- rep(NA, ncell(x))
>      # compute Pettitt  Change-Point Detection
>    #foreach(i=1:ncell(x)) %dopar% {}
> for(i in 1:ncell(x)) {
>      temporal <- x[i]
>      if(!any(is.na(temporal))) {
>        ptt[i] <- pettitt.test(x[i])[3] # Pettitt Test
>      }
>      cat("\r Processing :", round(i/ncell(x) * 100), "%")
>    }
> 
>    PTT <- raster(nrow=nrow(x), ncol =ncol(x), crs=crs(x))
>    extent(PTT) <- extent(x)
>    PTT[] <- unlist(ptt)
>      # write output
>    raster::writeRaster(PTT, "PTT.nc", overwrite = T)
> }
> 
> start <- Sys.time()
> Pettitt(s, dates)
> end <- Sys.time()
> difftime(end,start)
> 
> Ptt<-brick("PTT.nc")
> plot(Ptt)
>