r - 将栅格值(从堆栈)提取到 for 循环中的点

Question

我有一个光栅堆栈和 100 个点。对于每个栅格，我想提取值并使用三种不同的比例/缓冲区来提取值。

首先，这是组合成一个堆栈的三个栅格

library(raster)
# Make rasters and combine into stack
set.seed(123)
r1 = raster(ncol=1000, nrow=1000, xmn=0, xmx=1000, ymn=0, ymx=1000)
values(r1) = round(runif(ncell(r1),1,100))

r2 = raster(ncol=1000, nrow=1000, xmn=0, xmx=1000, ymn=0, ymx=1000)
values(r2) = round(seq(1:ncell(r1)))

r3 = raster(ncol=1000, nrow=1000, xmn=0, xmx=1000, ymn=0, ymx=1000)
values(r3) = round(runif(ncell(r1),1,5))

RasterStack <- stack(r1, r2, r3)

然后我生成 100 个点作为 SpatialPoints 对象

#make points
Points <- SpatialPoints(data.frame(xPoints = sample(1:1000, 100),
                                   yPoints = sample(1:1000, 100)))

接下来，我定义要循环的三个缓冲区

Scales <- c(60, 500)

为了更好地描述所需的结果，我将首先仅使用单个栅格，而不是 RasterStack。下面的代码定义了一个矩阵(输出)，该矩阵填充在循环中，每一列都是在两个不同的Scales 处提取的r1 值。然后将这些列标记在循环外部。

output <- matrix(ncol = length(Scales), nrow = length(Points))
for( i in 1:length(Scales)) {
  output[, i] <- extract(r1, Points, method='simple', buffer=Scales[i], fun=mean)
}
colnames(output) <- paste("r1", Scales, sep = "_" )

   > head(output)
        r1_60   r1_500
[1,] 50.67339 50.42280
[2,] 50.42401 50.42335
[3,] 49.96709 50.44288
[4,] 50.65492 50.52634
[5,] 50.60678 50.43535
[6,] 50.52477 50.48277

我想要相同的输出，但我不想调用单个栅格(例如上面的 r1)，而是想对 RasterStack 中的每个栅格执行此操作。最终结果将是一个矩阵(或 data.frame)，每个栅格 (r1:r3) 有两列。如示例中所示，标签将对应于相应的比例，以便将列标记为 r1_60、r1_500、r2_60、...、r3_500。

我认为嵌套的for循环可以在我循环遍历RasterStack和Scales的地方工作，但怀疑可能有更好的方法.

对于真实数据，我有 20 个 1541 x 1293 的栅格和大约 30,000 个位置。我还有 5 个不同的比例，因此嵌套的 for 循环将需要很长时间才能运行。

添加 采用不同的方法，我可以使用以下代码创建数据帧列表，每个数据帧对应于使用给定缓冲区提取的每层的值。

output <- list()
for(i in 1:length(Scales)){
  output[[i]] <- extract(RasterStack, Points, method='simple', buffer = Scales[i], fun = mean)
  names(output)[[i]] <- paste("Buffer", Scales[i], sep = "_")
}

根据此输出，我如何制作单个 6 x 100 数据帧，其中每一列都被标记为“layer_buffer 编号”。例如，layer.1_60、layer.2_60、...、layer.2_500、layer.3_500。

我还可以发布一个新的首选问题。

Answer 1

raster 包中似乎存在一个错误，如果 buffer< 表示的距离从 RasterStack 中提取值时，会导致抛出错误 小于网格分辨率。这也可引用here .

例如，

extract(RasterStack, Points, buffer=0, fun=mean)
## Error in apply(x, 2, fun2) : dim(X) must have a positive length

解决方法有点困惑:

# Just the first 10 points, for the example
Points <- Points[1:10, ]

dat <- do.call(cbind, lapply(Scales, function(b) {
  out <- do.call(rbind, lapply(extract(RasterStack, Points, buffer=b), 
                               function(x) if(is.matrix(x)) colMeans(x) else x))
  colnames(out) <- paste(colnames(out), b, sep='_')
  out
}))

这会产生:

dat
##       layer.1_0 layer.2_0 layer.3_0 layer.1_60 layer.2_60 layer.3_60 layer.1_500 layer.2_500 layer.3_500
##  [1,]        48    409158         4   50.67339   408657.5   3.013623    50.42280    435485.7    2.999983
##  [2,]        80    450287         1   50.42401   449786.5   2.990888    50.42335    460519.9    2.999632
##  [3,]        89    987912         3   49.96709   968829.9   2.995279    50.44288    775273.5    3.002715
##  [4,]        65    119952         5   50.65492   119448.9   3.009086    50.52634    273116.8    3.000364
##  [5,]        99    142320         4   50.60678   141819.5   2.998585    50.43535    289803.0    2.999054
##  [6,]        64    394804         3   50.52477   394303.5   2.984253    50.48277    426887.0    3.000055
##  [7,]        61    580925         2   50.96037   580424.5   3.001769    50.50032    559294.6    2.999218
##  [8,]        47     84918         3   50.83050    84417.5   2.998585    50.51135    258470.6    2.999923
##  [9,]         8    750667         4   50.16003   750166.5   2.987969    50.41984    655768.4    3.000635
## [10,]        88    273369         5   50.30219   272868.5   2.981157    50.44709    354833.6    2.999274