Assume I have four raster files: x, y, w and z and I want to create a new raster file based on the relationship of all of them (i.e., the hierarchical order of grid cells) using multiple nested conditionals as per function below:
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))
This is what I came up with so far. Yet, it does not seem to work. Wondering if anyone could give me a hand and shed some light on how I could sort it out?
Error message:
Error in (function (x, fun, filename = "", recycle = TRUE, forcefun = FALSE, : cannot use this formula, probably because it is not vectorized"
Many thanks.
There are a couple of problems here. Firstly, your ifelse statements have too many arguments (the first 3 have an unnecessary NA at the end. Secondly, you need to compare the data within each raster using e.g. x[] < y[] rather than x < y. Thirdly, you should have a single & rather than &&.
This way, your function is already vectorised:
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))))
}
So you can do:
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 ]
and
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)