我有一个由温度数据组成的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(即第一层)。任何想法将不胜感激。

最佳答案

这是一种更正式(更安全,更明智的存储方式)的方法:

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

关于r - 使用另一个栅格的像元值按索引对栅格堆栈进行设置,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/34025809/

10-12 22:31