基于另一图层中的值进行栅格图层的条件分析

Question

如果我有一个栅格堆栈：

require(raster)
r_test <- stack(
  raster(ncols = 10, nrows = 10, vals = sample(c(-1, 0, 1), 100, TRUE)),
  raster(ncols = 10, nrows = 10, vals = rnorm(100, 7, 0.4)

我想在layer.1中的同一单元格中基于其值将一个函数应用于layer.2，请问该如何操作？

举个简单的例子，将layer.1中等价单元格值为-1的layer.2中的值相乘。

- sinclairjesse

---

嗯，也许可以使用overlay来实现这个？ overlay(x = r_test[[1]], y = r_test[[2]], fun = function(x, y) ifelse(x == -1, x * y, y)) - Roman Luštrik

谢谢@RomanLuštrik。我也是这么做的。它有效，所以请将其写成答案！ - sinclairjesse

Answer 1

一般来说，overlay() 函数在许多情况下都能很好地发挥作用。与我首次回答时相矛盾的是，在我的经验中，对于更大的光栅数据，ifelse() 可能要更高效或更低效。对于一个包含 1000 行/列的光栅数据，简单的光栅代数会更快。如果我将其增加到 10000 行/列，则调用 ifelse() 更好。

library(raster)
r_test <- stack(
  raster(ncols = 1000, nrows = 1000, vals = sample(c(-1, 0, 1), 1000000, TRUE)),
  raster(ncols = 1000, nrows = 1000, vals = rnorm(1000000, 7, 0.4))
)

##  While overlay() works this is a more general solution and is
##      more efficient for large raster data sets:
system.time(
r_out <- r_test[[1]] * r_test[[2]] * (r_test[[1]] == -1) + r_test[[2]] * (r_test[[1]] != -1)
)
# N=1000x1000 cells:
#   user  system elapsed 
#   0.05    0.01    0.06
# N=10000x10000 cells:
#   user  system elapsed 
#   48.36   19.60   77.56 

system.time(
r_out <- overlay(x = r_test[[1]], y = r_test[[2]], fun = function(x, y) ifelse(x == -1, x * y, y))
)
# N=1000x1000 cells:
#   user  system elapsed 
#   0.53    0.08    0.64
# N=10000x10000 cells:
#   user  system elapsed 
#   19.77    5.82   26.53  

- Forrest R. Stevens

---

奇怪，不是吗？为什么会这样呢？ - sinclairjesse

我猜测这可能与硬件和数据相关的行为有关。但请注意，我在一台内存对于所提供的数据没有任何显著限制的机器上运行以上代码，因此我不认为这会是一个巨大的问题。话虽如此，如果不深入挖掘细节，我无法确定。 - Forrest R. Stevens