来自 "smwrQW" 包的带有 KendallATS 函数的栅格计算

Question

这是我在 Whosebug 上的第一个问题。我正在尝试对由 15 层（计数的时间序列）组成的栅格块使用 Mann-Kendall 趋势检验来制作三个栅格（sen 斜率、tau-b 和 p 值）。我需要使用 tau-b 而不是 tau-a，因为我的一些数据基本上是绑定的并且使用 tau-a 产生的输出（例如 Envstats 包中的函数 kendallTrendTest(x) ）有时是异常的（例如，显着的零斜率）。我在 R 包中找到的唯一函数，它完全满足我的需要，是 USGS 包 smwrQW 中的 kendallATS.test(x,y)。问题是，我无法使此函数与 raster 包中的 calc 函数一起使用，尽管我删除了研究区域之外的所有 NA 并且所有系列完全由相同的值组成（例如 c(0,0,0,0,0)).我对操作进行了不同的编码，但徒劳无功。以下是我使用的一些示例：

library(raster)
library(smwrQW)

setwd("C:R/WorkingDirectory/Annual_SCD")
list <- list.files(pattern = "*.tif")
RasterStack <- stack(list)
rb <- brick(RasterStack)

# Remove NA outside study area
j <- calc(rb, function(x) {x[is.na(x)] <- -999; return(x)}) 

# MK tau-b Slope (PROBLEM)
m <- calc(j, function(x) { y <- 1:15; X <- all(x == x[1]); if (X == FALSE) {
kendallATS.test(c(x[]),y)$estimate[1] } else return(-999)}, forcefun=TRUE)

获得错误信息：

Error in x@.Data[i, , drop = FALSE] : 
  (subscript) logical subscript too long

我研究区域内的所有单元格都包含砖块中 15 层的数值数据。我使用 c(x[]) 而不是 x 来获得向量，否则，我会收到另一条错误消息：

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘as.lcens’ for signature
 ‘"matrix", "missing", "missing"’.

as.lcens 函数属于 smwrQW。我已经在 kendallATS.test 和 calc 中使用了 na.rm=TRUE，但没有改变。然而，该函数与指定的行和列一起完美地工作：

kendallATS.test(c(rb[250,300]),1:15)

我的编码知识非常基础，如果能帮助我解决这个问题，我将不胜感激。

谢谢

Answer 1

这是一个可行的示例，也许可以作为起点：

library(raster)
library(smwrQW)

b <- brick(ncol=10, nrow=10, nl=15)
b <- setValues(b, matrix(runif(1200), 100, 15))
b[1:5] <- NA

y <- 1:nlayers(b)
f <- function(x) { 
    if (any(is.na(x))) {
        return(c(NA, NA, NA))
    } else if (any(x != x[1])) {
        kendallATS.test(x , y)$estimate 
    } else {
        return(c(NA, NA, NA))
    }
}

m <- calc(b, f)

如果您想 return 单个变量，请注意下面如何编辑函数：

fp <- function(x) { 
    if (any(is.na(x))) {
        return(NA)
    } else if (any(x != x[1])) {
        kendallATS.test(x , y)$p.value
    } else {
        return(NA)
    }
}

n <- calc(b, fp)

来自 "smwrQW" 包的带有 KendallATS 函数的栅格计算

Raster Calc with KendallATS function from "smwrQW" package

r

raster

calc

subscript

r-raster