R: Annual composite raster based on median: 如何获取每个像素的原始图层的索引？

Question

我有一个 rasterstacks 列表，它们是不同波段的时间序列（300+ 层）。时间序列是不规则的，因此我想根据中值 ndvi 创建年度复合材料。因此，从一年的可用图像中，选择具有（接近）中值 ndvi 的像素。对于其他波段，我想基于此 ndvi 合成材料创建年度合成材料。我尝试制作一个掩码，其中每个像素都具有用于中值 ndvi 合成的图像索引值。我将把它应用到其他波段，所以我每年每个波段都有 'the same' 年度综合。

不幸的是，我一直在做面具。我创建了一些虚拟数据，并以某种方式 returns 两个索引栅格（一个值为 1-3，另一个为 2-4），而我期望有一个（值为 1-4）。此外，我的函数无法处理 NA 值，将 na.rm 添加到 calc 函数也无法解决此问题。我想知道我需要调整什么以获得一个值为 1-4 的输出层，以及如何让 'which'-函数处理 NA。

#dummy data:
r <- raster(nrow=5, ncol=5)
set.seed(20181114)
s <- stack(lapply(1:5, function(i) setValues(r, runif(25, max=50))))
names(s) <- c("X2001", "X2002a", "X2002b", "X2002c", "X2002d")
#s$X2002a[2] <- NA 

AnnualMask <- function(ts, year){
  year <- as.character(year)
  ts_year <- subset(ts, (grep(year, names(ts))))
  indexraster <- calc(ts_year, function(x){
    medianval <- median(x, na.rm = TRUE)
    which(abs(x - medianval) == min(abs(x - medianval)))
           })
  return(indexraster)
}

mask2002 <- AnnualMask(s, 2002)
plot(mask2002)

Answer 1

我还不确定值在 1-4 之间的部分，但是这里有一些关于 which 部分和处理 NA:

# here is some dummy data
x <- c(3,4,5,NA,1,-2-45,2,54)

# now here is the which part as found in your function
medianval <- median(x, na.rm = TRUE)
which(abs(x - medianval) == min(abs(x - medianval)))
# returning in this case
integer(0)

# now let's add na.rm = T in the min function
medianval <- median(x, na.rm = TRUE)
which(abs(x - medianval) == min(abs(x - medianval), na.rm = T))
# returning
[1] 1

Answer 2

感谢您提供很好的示例数据：

library(raster)
r <- raster(nrow=5, ncol=5)
set.seed(20181114)
s <- stack(lapply(1:5, function(i) setValues(r, runif(25, max=50))))
names(s) <- c("X2001", "X2002a", "X2002b", "X2002c", "X2002d")
s[[2]][1:3] <- NA

这里有两个功能相同。 #1 可能更简单，但效率可能较低。

annualIndex1 <- function(ts, year) {
    year <- as.character(year)
    ts_year <- subset(ts, (grep(year, names(ts))))
    medianval <- calc(ts_year, median, na.rm = TRUE)
    dif <- abs(ts_year - medianval)
    which.min(dif)
}
x1 <- annualIndex1(s, 2002)

这个使用 calc 组合多个步骤

annualIndex2 <- function(ts, year){
    year <- as.character(year)
    ts_year <- subset(ts, (grep(year, names(ts))))

    fun <- function(x){
        y <- median(x, na.rm=TRUE)
        which.min(abs(x - y))
    }
    calc(ts_year, fun)
}  
x2 <- annualIndex2(s, 2002)

结果是一样的（而且没有NA）

all(values(x1) == values(x2))
#[1] TRUE

但是，如果您在所有层中都设置一个单元格 NA

s[5] <- NA

功能 2 失败。可以这样调整

annualIndex2b <- function(ts, year){
    year <- as.character(year)
    ts_year <- subset(ts, (grep(year, names(ts))))

    fun <- function(x){
        if (all(is.na(x))) return( NA )
        y <- median(x, na.rm=TRUE)
        which.min(abs(x - y))
    }
    calc(ts_year, fun)
}

x2b <- annualIndex2b(s, 2002)

您使用了术语 "mask"，但实际上您是在创建索引，您可以像这样在 strackSelect 中使用适当的单元格吗：

ts_year <- subset(ts, (grep(year, names(ts))))
v <- stackSelect(ts_year, x2)

R: Annual composite raster based on median: 如何获取每个像素的原始图层的索引？

R: Annual composite raster based on the median: How to get the index of the original layer for each pixel?

r

mask

raster

median

na