提取柵格值於點中用於環路

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 

我想要這個相同的輸出，但而不是調用單個光柵（例如R 1以上），我想在RasterStack每個柵格執行此操作。最終的結果將是一個矩陣（或data.frame），其具有兩列用於每個光柵（R1：R3）。如例子中，標籤將對應於相應的規模，這樣的列被標記r1_60, r1_500, r2_60, ... , r3_500.

我覺得嵌套for循環將在哪裏工作我遍歷RasterStack並通過Scales但懷疑有可能是一個更好的辦法。

實際的數據我有20個光柵是由1541和1293 3萬點左右的位置。我也有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 = "_") 
} 

從這個輸出，我怎樣才能使其中每列將被標記爲「layer_buffer號碼」的單個6由100數據幀。例如，layer.1_60，layer.2_60，...，layer.2_500，layer.3_500。

我也可以張貼的首選的新問題。

Answer 1

對於關閉的緣故，我張貼的解決方案，爲我工作最好。根據柵格包錯誤，我沒有使用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 = "_") 
} 

然後，跟隨的post linked here，我用下面的代碼來的數據幀的列表合併成單個數據幀。

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 

Answer 2

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