提取栅格值(从堆栈)到for循环中的点 [英] Extract raster values (from Stack) to points in for loop
问题描述
我有一个栅格堆栈和100分.对于每个栅格，我要提取值，并使用三个不同的比例尺/缓冲区进行提取.
I have a raster stack and 100 points. For each raster I want to extract the value and do so using three different scales/buffers.
首先，这是将三个栅格合并为一个堆栈
First, here are three rasters combined into a stack
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 对象
I then generate 100 points as a SpatialPoints object
#make points
Points <- SpatialPoints(data.frame(xPoints = sample(1:1000, 100),
yPoints = sample(1:1000, 100)))
接下来，我定义要循环通过的三个缓冲区
Next, I define the three buffers that I want to loop through
Scales <- c(60, 500)
为了更好地描述期望的结果，我将首先仅使用单个栅格，而不使用RasterStack.下面的代码定义了一个矩阵(输出)，该矩阵在循环中填充，每一列都是在两个不同的 Scales 上提取的 r1 的值.然后将这些列标记在循环外部.
To better describe the desired outcome, I will first use only a single raster, not the RasterStack. The code below defines a matrix (output) which is populated in the loop with each column being the extracted values of r1 at the two different Scales. The columns are then labeled outside of the loop.
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.
I want this same output, but rather than calling a single raster (e.g. r1 above), I want to do this for each raster in the RasterStack. The final result would be a matrix (or data.frame) that has two columns for each raster (r1:r3). As in the example, labeling would correspond to the respective scale so that the columns were labeled r1_60, r1_500, r2_60, ... , r3_500.
我认为嵌套的 for 循环可以在通过 RasterStack 和通过 Scales 循环的地方工作，但我怀疑可能会有更好的方法
I think a nested for loop would work where I loop through the RasterStack and through the Scales but suspect there might be a better way.
对于真实数据，我有20个栅格，这些栅格分别是1293年和1541年的30,000个位置.我也有5种不同的音阶，因此嵌套的 for 循环将花费很长时间.
For the real data I have 20 rasters that are 1541 by 1293 and around 30,000 locations. I also have 5 different scales so a nested for loop will take a very long time to run.
添加采用不同的方法，我可以使用以下代码创建数据帧列表，每个数据帧对应于使用给定缓冲区提取的每一层的值.
Addition Taking a different approach, I can use the following code to create a list of data frames, each of which corresponds to the extracted values of each layer using a given buffer.
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.
From this output, how can I make a single 6 by 100 data frame where each column would be labeled as the "layer_buffer number". For example, layer.1_60, layer.2_60, ... , layer.2_500, layer.3_500.
我还可以发布一个首选问题.
I can also post a new question of preferred.
推荐答案
为方便起见，我发布了最适合我的解决方案.鉴于栅格数据包错误，我没有使用0缓冲区将值提取到点.
For the sake of closure, I am posting the solution that worked best for me. In light of the raster package bug, I did not extract values to points using the 0 buffer.
Scales <- c(60, 500)
然后，使用前10点，
Points <- Points[1:10]
我使用以下代码为每个缓冲区级别创建了一个列表.
I created a list for each buffer level using the following code.
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 = "_")
}
Then, following the post linked here, I used the following code to combine the list of data frames into a single data frame.
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
这篇关于提取栅格值(从堆栈)到for循环中的点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
