与多个栅格相关的嵌套条件语句

Question

假设我有四个光栅文件：x、y、w 和 z，并且我想根据所有这些文件的关系（即网格单元的层次顺序）使用多个嵌套条件创建一个新的光栅文件下面的函数：

x <- raster(nrows=100, ncols=100)
  x[] <- runif(ncell(x), min = 0, max = 100)
y <- raster(nrows=100, ncols=100)
  y[] <- runif(ncell(y), min = 0, max = 100)
w <- raster(nrows=100, ncols=100)
  w[] <- runif(ncell(w), min = 0, max = 100)
z <- raster(nrows=100, ncols=100)
  z[] <- runif(ncell(z), min = 0, max = 100)

My.fun <- function(x,y,z,w){
ifelse(x > y && x > w && x > z,
1,
ifelse(y > x && y > w && y > z,
2,
ifelse(w > x && w > y && w > z,
3,
ifelse(z > x && z > y && z > w,
4,NA),NA),NA),NA)
}

res <- overlay(x, y, z, w, fun = Vectorize(My.fun))

这是我到目前为止想出的。然而，它似乎不起作用。想知道是否有人可以帮我解答一下我该如何解决？

错误信息：

(函数 (x, fun, filename = "", recycle = TRUE, forcefun = FALSE, : 不能用这个公式，可能是因为没有向量化”

非常感谢。

Answer 1

这里有几个问题。首先，您的 ifelse 语句有太多参数（前 3 个在末尾有一个不必要的 NA。其次，您需要使用例如 x[] < y[] 而不是 x < y 比较每个栅格内的数据。第三，您应该有一个 & 而不是 &&。

这样，你的函数已经被向量化了：

My.fun <- function(x,y,z,w){
  ifelse(x[] > y[] & x[] > w[] & x[] > z[], 1,
         ifelse(y[] > x[] & y[] > w[] & y[] > z[], 2,
                ifelse(w[] > x[] & w[] > y[] & w[] > z[], 3,
                       ifelse(z[] > x[] & z[] > y[] & z[] > w[], 4, NA))))
}

所以你可以这样做：

My.fun(x, y, w, z)
#>    [1] 3 1 3 2 3 4 2 3 1 3 4 4 3 4 1 2 2 1 2 3 4 3 4 4 1 1 3 2 3 1 2 4 2 4 2 4
#>   [37] 1 4 4 2 2 4 1 2 3 2 2 3 1 1 1 2 3 4 4 4 3 4 4 2 3 2 2 4 2 2 3 1 1 4 1 4
#>   [73] 2 2 2 3 3 2 1 4 1 3 3 4 4 4 4 1 4 3 2 1 4 1 1 2 3 1 3 3 1 3 4 4 1 1 4 3
#>  [109] 4 3 1 4 4 1 2 3 1 3 2 4 1 4 1 2 2 3 2 3 2 1 3 1 3 2 2 3 3 1 3 1 3 3 3 2
#>  [145] 1 4 4 3 2 1 3 1 3 1 1 1 4 2 3 1 1 4 2 3 1 3 2 2 2 2 3 1 1 3 4 1 2 1 1 2
#>  [181] 2 3 4 1 4 1 3 3 1 4 3 2 3 1 1 2 4 3 1 3 2 1 4 2 4 3 2 1 1 1 1 1 1 1 2 3
#>  [217] 4 4 2 1 2 2 1 3 2 3 3 3 4 3 2 1 2 2 4 2 4 4 2 2 4 3 4 3 1 2 4 4 4 3 2 2
#>  [253] 2 4 3 4 4 3 1 3 2 4 3 2 2 2 4 3 3 4 4 3 3 4 3 4 3 1 1 1 1 3 2 3 3 3 1 2
#>  [289] 1 4 4 4 3 4 3 4 3 2 2 3 1 3 1 1 1 1 3 2 4 4 4 1 1 3 4 4 4 3 4 1 2 1 1 4
#>  [325] 4 4 2 4 2 3 4 4 2 3 1 1 1 4 3 2 3 4 4 1 3 3 4 1 3 1 2 4 1 1 2 1 2 4 2 4
#>  [361] 3 3 2 2 1 1 4 2 1 3 4 4 3 1 2 2 3 4 2 4 2 3 1 4 3 3 3 4 2 2 1 2 2 1 4 1
#>  [397] 1 1 3 4 3 1 2 1 1 2 3 3 4 2 1 1 1 3 3 1 2 2 1 3 1 4 1 4 2 3 2 2 1 4 1 4
#>  [433] 4 1 3 3 1 1 1 1 3 1 2 1 1 4 2 3 4 4 2 1 3 4 4 4 4 4 2 4 2 3 4 2 2 3 1 4
#>  [469] 3 1 3 3 3 2 1 4 1 4 2 2 3 1 3 2 1 4 3 3 2 4 4 4 3 3 1 2 1 4 4 1 3 1 2 3
#>  [505] 2 1 4 2 4 2 4 4 4 4 2 3 2 3 2 2 1 3 2 3 1 3 1 3 2 1 3 1 2 2 2 4 3 4 2 4
#>  [541] 1 4 1 2 3 3 3 2 1 1 3 4 2 1 1 4 4 2 4 4 2 2 3 4 1 4 1 3 2 4 3 3 1 2 2 1
#>  [577] 2 2 3 1 3 1 1 3 4 2 1 3 3 3 3 1 1 2 3 2 1 1 1 4 1 2 1 3 4 4 4 1 1 3 1 1
#>  [613] 3 3 1 1 3 3 1 2 1 1 4 4 1 2 1 3 2 3 1 4 2 3 3 4 3 4 1 2 3 4 3 3 2 4 3 2
#>  [649] 1 3 2 4 2 4 3 1 1 1 2 4 3 2 2 4 4 2 4 2 4 1 3 2 2 3 2 2 2 3 4 1 1 3 2 3
#>  [685] 1 2 4 3 2 2 4 4 3 2 1 2 4 3 2 4 4 2 2 1 2 3 1 1 3 3 3 1 4 2 2 2 2 3 2 3
#>  [721] 3 3 2 1 4 3 1 1 3 4 4 2 1 3 4 3 3 4 4 1 3 2 4 3 2 3 1 2 3 4 3 4 3 3 3 2
#>  [757] 4 3 2 3 1 2 4 4 4 3 4 3 4 2 3 2 4 4 4 4 2 4 2 1 1 3 4 3 3 1 3 4 3 4 1 1
#>  [793] 4 1 4 3 4 3 4 4 4 1 3 1 2 3 3 3 4 4 3 4 3 2 3 3 4 1 1 4 4 3 1 1 2 1 2 2
#>  [829] 2 1 1 3 2 4 4 1 3 2 3 4 2 3 2 3 2 3 4 4 4 4 1 3 2 3 2 2 3 2 4 1 1 4 4 2
#>  [865] 3 2 2 3 3 1 2 3 3 1 2 4 4 1 3 2 3 2 4 4 2 1 2 1 3 4 1 2 1 1 3 3 1 1 1 3
#>  [901] 3 4 2 1 1 4 4 4 2 1 3 3 3 1 4 3 1 2 3 1 1 1 3 2 1 4 1 1 4 1 4 1 1 3 1 2
#>  [937] 2 2 2 3 2 2 2 2 3 2 4 3 2 2 2 3 2 4 2 4 1 1 4 2 2 4 4 1 4 4 4 4 2 4 2 4
#>  [973] 4 1 4 3 1 2 3 1 3 2 3 3 2 1 3 1 2 1 2 3 1 3 4 3 3 2 3 2
#>  [ reached getOption("max.print") -- omitted 9000 entries ]

和

wxyz <- raster(nrows=100, ncols=100)
wxyz[] <- My.fun(x,y,z,w)
wxyz
#> class      : RasterLayer 
#> dimensions : 100, 100, 10000  (nrow, ncol, ncell)
#> resolution : 3.6, 1.8  (x, y)
#> extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 
#> source     : memory
#> names      : layer 
#> values     : 1, 4  (min, max)