我想计算多层中每个像素的分位数
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.5, 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
包有一个用于SpatRaster
的quantile
方法。
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, 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)
但我会这样做:
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
你可以用terra::app
使用同样的函数。
zz <- app(ss, qf)