enable parallel analysis on pixels from rasterbrick

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

enable parallel analysis on pixels from rasterbrick

Jackson Rodrigues
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)

--

Jackson M. Rodrigues

"In order to succeed, we must first believe that we can."

Nikos Kazantzakis

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: enable parallel analysis on pixels from rasterbrick

Loïc Dutrieux-3
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)
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo