高效地從R中的柵格堆棧訪問數據

我需要從grib2文件（我將其放入柵格堆棧）中有效地提取溫度數據。堆棧中的每個柵格圖層都表示一個時間點。高效地從R中的柵格堆棧訪問數據

現在我需要爲每個觀察值（x，y，t）提取一個單值。下面的代碼可以完成這項工作，但這需要花費太多時間。任何提高效率的建議非常感謝。

files <- list.files(path="Weather/NCEP/temperature_3hour_forecast", full.names = TRUE) 

s <- stack(files) 

userdata$x <- sample(1:ncol(s), nrow(userdata), replace=T) 
userdata$y <- sample(1:nrow(s), nrow(userdata), replace=T) 

smalldata <- userdata[1:100,] 
tic() 
smalldata$temp1morning <- getValues(s[[smalldata$t]], smalldata$y)[smalldata$x] 
toc()

編輯：loki的答案是非常有用的。但是，當我將這種方法應用於我的溫度數據時，它確實很慢。我懷疑這是由我的溫度數據結構引起的，並且有許多時間段。請參閱下面的建議方法與getValues的可比試驗之間的比較。任何想法，爲什麼這是這種情況，或者我可以如何改進代碼？

> files <- list.files(path="Weather/NCEP/temperature_3hour_forecast", full.names = TRUE, pattern = glob2rx("*06.f003.grib*")) 
> 
> s <- stack(files) 
> s 
class  : RasterStack 
dimensions : 197, 821, 161737, 971 (nrow, ncol, ncell, nlayers) 
resolution : 0.25, 0.25 (x, y) 
extent  : 190.875, 396.125, 22.875, 72.125 (xmin, xmax, ymin, ymax) 
coord. ref. : +proj=longlat +a=6371229 +b=6371229 +no_defs 
names  : gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, gfs.0p25.//mayr258302, ... 

> 
> userdata$x <- sample(1:ncol(s), nrow(userdata), replace=T) 
> userdata$y <- sample(1:nrow(s), nrow(userdata), replace=T) 
> 
> smalldata <- data.frame(x = userdata$x[1:2], 
+ y = userdata$y[1:2], 
+ t = userdata$t[1:2]) 
> 
> smalldata 
    x y t 
1 142 67 547 
2 779 14 829 
> 
> tic("apply") 
> smalldata$temp1morning <- apply(smalldata, 1, function(x){s[x[2], x[1]][x[3]]}) 
> toc() 
apply: 305.41 sec elapsed 
> 
> tic("getValues") 
> smalldata$temp2morning <- apply(smalldata, 1, function(x){getValues(s[[x[3]]], x[2])[x[1]]}) 
> toc() 
getValues: 0.33 sec elapsed 
> 
> smalldata 
    x y t temp1morning temp2morning 
1 142 67 547 13.650018 13.650018 
2 779 14 829 -1.750006 -1.750006 
>

來源

2017-09-14 Ha Ma

請拿出一個可重複的例子。 [這篇文章]（https://stackoverflow.com/q/5963269/3250126），特別是[這個答案]（https://stackoverflow.com/a/5963379/3250126）將幫助你。 – loki

讓我們先從一個可重複的例子：

library(raster) 
r <- raster(ncol = 100, nrow = 100) 
r[] <- runif(ncell(r)) 

s <- stack(r, r, r) 
s

現在，我們假設你的userdata持有以下結構：

x是像素行的索引
y是像素列的索引
t是該層的指標，我們在

尋找像素讓我們創建一個可重複的用戶數據：

userdata <- data.frame(x = sample(1:100, 10), 
         y = sample(1:100, 10), 
         t = sample(1:3, 10, replace = T))

然後我們可以使用apply通過所有的線userdata工作，使用行，列，和層的索引來提取值：

userdata$pixelvalue <- apply(userdata, 1, function(x){s[x[1], x[2]][x[3]]})

在apply的像素是本身的每次迭代通過其所有層的柵格中的x和y位置來引導。 x[3]然後只返回相應圖層的值。

在此之前的邏輯：

stack[*row*, *column*][*layer*]

你的方法的好處是，你不必整個光柵轉換爲矢量（這基本上是什麼getValues一樣），而是直接訪問數據作爲RasterStack的矩陣結構。

來源

2017-09-14 07:55:07 loki

感謝您對loki的快速反應。不幸的是，你提出的代碼對我的數據非常慢（見上面的編輯）。你有一個想法是如何發生的？ –

我發現了一個適用於我的簡單解決方案。首先，我使用as.array將我的溫度數據輸入到一個數組中。然後，我在陣列上使用apply爲洛基建議：

files <- list.files(path="Weather/NCEP/temperature_3hour_forecast", full.names = TRUE, pattern = glob2rx("*06.f003.grib*")) 

s <- stack(files) 
a <- as.array(s) 

userdata$x <- sample(1:ncol(s), nrow(userdata), replace=T) 
userdata$y <- sample(1:nrow(s), nrow(userdata), replace=T) 

smalldata <- data.frame(x = userdata$x[1:nrow(userdata)], 
      y = userdata$y[1:nrow(userdata)], 
      t = userdata$t[1:nrow(userdata)]) 

tic("array") 
userdata$temp1morning <- apply(smalldata, 1, function(x){a[x[2], x[1], x[3]]}) 
toc()

這是很容易足夠快，我的目的。洛基，謝謝你的幫助！

來源

2017-09-15 06:59:23

高效地從R中的柵格堆棧訪問數據

回答

相關問題