计算多个栅格层中每个像素的分位数 probs 0.05

如何解决计算多个栅格层中每个像素的分位数 probs 0.05

我想计算多个层中每个像素的分位数

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()

我希望有三层，但它会产生默认的概率（参见 stats::quantile()）

我还尝试使用特定参数生成函数分位数，如

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

但出现以下错误

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

解决方法

我在玩，没想到这么简单

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()

然后我重复了 prob = 0.05

,

terra 包具有 quantile 的 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 

这应该比 raster 更快；但是你可以像这样使用raster

得到相同的结果

library(raster)
r <- raster(ncols=36,nrows=18)
values(r) <- 1:ncell(r)
s <- stack(r,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,.95))
  }
}

z <- calc(s,qfun)

但我会这样做：

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

z
#class      : RasterBrick 
#dimensions : 18,648,ncell,nlayers)
#resolution : 10,y)
#extent     : -180,ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      :        X5.,X50.,X95. 
#min values :        1.0,1.9 
#max values :   87.71026,1231.20000 

您可以对 terra::app 使用相同的函数。

zz <- app(ss,qf)