CODEWITHSUNDEEP
rrasterquantilecalcmulti-layer
Marcos
Marcos

Reputation: 51

Compute quantiles for each pixel across a number of raster layers probs 0.05

I want to compute quantiles for each pixel across many layers

library(raster)

r <- raster(ncols=36, nrows=18)
r[] <- 1:ncell(r)
s <- stack(r, r*2, sqrt(r))
s[10:20] <- NA

beginCluster()
clusterR(s, calc, args=list(fun = quantile, na.rm= TRUE, prob = c(0.05, 0.5, 0.95)))
endCluster()

I would expect three layers, but instead it produces the default probs (see stats::quantile())

I also tried generating the function quantile with the specific arguments, as

function(x){quantile(x, probs = c(0.05, 0.5, 0.95)}

but got the following error

Error in is.infinite(v) : default method not implemented for type 'list'

Upvotes: 0

Views: 560

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47696

The terra package has a quantile method for SpatRaster.

library(terra)
rr <- rast(ncols=36, nrows=18)
values(rr) <- 1:ncell(rr)
ss <- c(rr, rr*2, sqrt(rr))
ss[10:20] <- NA

q <- quantile(ss, c(0.05, 0.5, 0.95))
q
#class       : SpatRaster 
#dimensions  : 18, 36, 3  (nrow, ncol, nlyr)
#resolution  : 10, 10  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#names       :      q0.05,       q0.5,      q0.95 
#min values  :        1.0,        1.0,        1.9 
#max values  :   87.71026,  648.00000, 1231.20000

That should be faster than with raster; but you can get the same result with raster like this

library(raster)
r <- raster(ncols=36, nrows=18)
values(r) <- 1:ncell(r)
s <- stack(r, r*2, sqrt(r))
s[10:20] <- NA

# improved version of your function
qfun <- function(x){
  if(is.na(x)){
    c(NA, NA, NA)
  }else{
    quantile(x, c(0.05, 0.5, .95))
  }
}

z <- calc(s, qfun)

But I would do this instead:

qf <- function(x){
    quantile(x, c(0.05, 0.5, .95), na.rm=TRUE)
}
z <- calc(s, qf)

z
#class      : RasterBrick 
#dimensions : 18, 36, 648, 3  (nrow, ncol, ncell, nlayers)
#resolution : 10, 10  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      :        X5.,       X50.,       X95. 
#min values :        1.0,        1.0,        1.9 
#max values :   87.71026,  648.00000, 1231.20000

And you could use the same function with terra::app.

zz <- app(ss, qf)

Upvotes: 0

Marcos
Marcos

Reputation: 51

I was playing around and did not realize that it was as simple as this

library(raster)

r <- raster(ncols=36, nrows=18)
r[] <- 1:ncell(r)
s <- stack(r, r*2, sqrt(r))
s[10:20] <- NA

q95 <- function(x){
  if(is.na(x)){
    NA
  }else{
    stats::quantile(x, .95)
  }
}

beginCluster()
clusterR(r, calc, args=list(fun = q95), verbose=TRUE)
endCluster()

And then I repeated it for prob = 0.05

Upvotes: 0

Related Questions