使用另一个栅格的单元值按索引对栅格堆进行子集提取

Question

使用另一个栅格的单元值按索引对栅格堆进行子集提取

3

我有一个由12层温度数据组成的光栅堆叠，每个月有一层。我还有一个特定作物的种植和收获日期的光栅。我想计算各种气候参数（例如温度标准差），但仅在种植/收获日期之间的种植期内进行，排除此期间之外的月份。对于温度光栅堆叠中不属于种植期的所有月份，我想将其分配为给定单元格的NA。更具体地说，我想使用我的种植/收获光栅中的单元格值按它们的索引（即第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（即第一层）。欢迎提出任何想法。

- jlab

2个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Robert Hijmans · Answer 1

这里有一种更正式（更安全，更节省内存）的方法：

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

- johannes · Answer 2

不太确定我是否理解了您的意思，也许这可以帮助您入门：

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