将栅格值（从 Stack 中）提取到 for 循环中的点
Extract raster values (from Stack) to points in for loop
我有一个光栅堆栈和 100 个点。对于每个栅格，我想提取值并使用三个不同的 scales/buffers.
首先，这里是三个栅格组合成一个堆栈
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。下面的代码定义了一个矩阵（输出），它在循环中填充，每一列都是 r1 在两个不同的 Scales 处的提取值。然后在循环外标记这些列。
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 number"。例如，layer.1_60, layer.2_60, ..., layer.2_500, layer.3_500.
我也可以post新问题首选
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
为了结束，我发布了最适合我的解决方案。鉴于光栅包错误，我没有使用 0 缓冲区将值提取到点。
Scales <- c(60, 500)
然后，使用前 10 个点，
Points <- Points[1:10]
我使用以下代码为每个缓冲区级别创建了一个列表。
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 = "_")
}
然后，在 之后，我使用以下代码将数据帧列表合并为一个数据帧。
do.call(cbind,lapply(names(output),function(x){
res <- output[[x]]
colnames(res) <- paste(colnames(res),x,sep="_")
res
}))
返回的df的head如下
layer.1_Buffer_60 layer.2_Buffer_60 layer.3_Buffer_60 layer.1_Buffer_500
[1,] 50.67339 408657.5 3.013623 50.42280
[2,] 50.42401 449786.5 2.990888 50.42335
[3,] 49.96709 968829.9 2.995279 50.44288
[4,] 50.65492 119448.9 3.009086 50.52634
[5,] 50.60678 141819.5 2.998585 50.43535
[6,] 50.52477 394303.5 2.984253 50.48277
layer.2_Buffer_500 layer.3_Buffer_500
[1,] 435485.7 2.999983
[2,] 460519.9 2.999632
[3,] 775273.5 3.002715
[4,] 273116.8 3.000364
[5,] 289803.0 2.999054
[6,] 426887.0 3.000055
我有一个光栅堆栈和 100 个点。对于每个栅格，我想提取值并使用三个不同的 scales/buffers.
首先，这里是三个栅格组合成一个堆栈
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。下面的代码定义了一个矩阵（输出），它在循环中填充，每一列都是 r1 在两个不同的 Scales 处的提取值。然后在循环外标记这些列。
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 number"。例如，layer.1_60, layer.2_60, ..., layer.2_500, layer.3_500.
我也可以post新问题首选
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
为了结束，我发布了最适合我的解决方案。鉴于光栅包错误，我没有使用 0 缓冲区将值提取到点。
Scales <- c(60, 500)
然后，使用前 10 个点，
Points <- Points[1:10]
我使用以下代码为每个缓冲区级别创建了一个列表。
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 = "_")
}
然后，在
do.call(cbind,lapply(names(output),function(x){
res <- output[[x]]
colnames(res) <- paste(colnames(res),x,sep="_")
res
}))
返回的df的head如下
layer.1_Buffer_60 layer.2_Buffer_60 layer.3_Buffer_60 layer.1_Buffer_500
[1,] 50.67339 408657.5 3.013623 50.42280
[2,] 50.42401 449786.5 2.990888 50.42335
[3,] 49.96709 968829.9 2.995279 50.44288
[4,] 50.65492 119448.9 3.009086 50.52634
[5,] 50.60678 141819.5 2.998585 50.43535
[6,] 50.52477 394303.5 2.984253 50.48277
layer.2_Buffer_500 layer.3_Buffer_500
[1,] 435485.7 2.999983
[2,] 460519.9 2.999632
[3,] 775273.5 3.002715
[4,] 273116.8 3.000364
[5,] 289803.0 2.999054
[6,] 426887.0 3.000055