使用另一个栅格的像元值按索引对栅格堆栈进行子集化

Question

我有一个由温度数据组成的 12 层栅格堆栈 - 每个月一层。我还有一个栅格，每个栅格用于某种作物的种植和收获日期。我想对各种气候参数（例如温度标准偏差）进行一些计算，但仅限于 planting/harvest 日期之间的种植期，不包括该期以外的月份。温度栅格堆栈中不属于裁剪期的所有月份，我想为给定单元格指定为 NA。更具体地说，我想使用 planting/harvest 栅格中的单元格值按索引（即第 1 层到第 12 层）而不是单元格值对裁剪月份进行子集化。我对这个束手无策。下面是我尝试在示例数据集上使用的代码，但没有给我想要的输出：

set.seed(150)
t1<- raster(ncol=10, nrow=5) #Create rasters
values(t1) = round(runif(ncell(t1),1,20))
t2<- raster(ncol=10, nrow=5)
values(t2) = round(runif(ncell(t2),1,20))
t3<- raster(ncol=10, nrow=5)
values(t3) = round(runif(ncell(t3),1,20))
t4<- raster(ncol=10, nrow=5)
values(t4) = round(runif(ncell(t4),1,20))
t5<- raster(ncol=10, nrow=5)
values(t5) = round(runif(ncell(t5),1,20))
t6<- raster(ncol=10, nrow=5)
values(t6) = round(runif(ncell(t6),1,20))

Planting<- raster(ncol=10, nrow=5) #Create Planting date raster
values(Planting) = round(runif(ncell(Planting),3,5)) # All planting dates are between 3 and 5

t.stack<-stack(t1,t2,t3,t4,t5,t6) #Create raster stack
layer.names<-c(1:6) # Rename raster stack layers from 1 to 6
names(t.stack)<-as.numeric(as.character(layer.names)) # Attempt to coerce raster stack layer names to numeric
names(t.stack) #View new raster stack layer names (they don't seem to be numberic)

t.stack[t.stack[["X1"]]< Planting] <- NA #Attempt to use Planting raster cell values to coerce raster layer index X1 (1st layer) to NA because it lies outside of the range of 3 to 5

head(t.stack[["X1"]]) #View new values of raster layer X1

第一个栅格堆栈层 (X1) 的输出是：

#   1  2  3  4  5  6  7  8  9 10
# 1 NA  9 10 17  6  8  8  8  9 17
# 2 12 19 16  9 14 19 12  6  5  7
# 3 10 NA NA NA 10 15  6 15  6 12
# 4 20 16 12  5 18 13 13  5 NA 10
# 5 13 14 18 10 16  8 10 NA 20  7

如您所见，我尝试将图层名称更改为数值，并使用条件语句根据 t.stack 中的单元格值查询它们。相反，如果 X1 图层中的特定像元值小于 Planting 栅格中的相应像元值，则会替换它们。但是，我希望将此数据集中的整个 X1 层指定为 NA，因为种植日期仅从 3 到 5，不包括 1（即第一层）。如有任何想法，我们将不胜感激。

Answer 1

不确定，如果我理解正确的话，也许这有助于你入门：

for (i in 1:nlayers(t.stack)) {
  tt <- t.stack[[i]][]
  tt[i < Planting[]] <-  NA
  t.stack[[i]][] <- tt
}

Answer 2

这里有一个更正式（更安全，记忆更明智）的方法：

set.seed(150)
r <- raster(ncol=10, nrow=5)
tt <- sapply(1:6, function(x) setValues(r,  round(runif(ncell(r),1,20))))
t.stack <- stack(tt)
Planting <- setValues(r, round(runif(ncell(r), 3, 5)))

f <- function(p) {
    x <- matrix(0, nrow=length(p), ncol=12)
    x[cbind(1:nrow(x), p)] <- 1
    x <- t(apply(x, 1, cumsum))
    x[x==0] <- NA
    x
}

p <- calc(Planting, f, filename='planting12.grd')

tp <- t.stack * p

使用另一个栅格的像元值按索引对栅格堆栈进行子集化

subsetting raster stack by index using cell values of another raster

indexing

r

raster

subset