0
votes

I am having three rasters.

library(raster)

a <- raster(ncol=100, nrow=100)
set.seed(2)     
values(a) = runif(10000, min=-1, max=1)    # define the range between -1 and 1

b <- raster(ncol=100, nrow=100)
set.seed(2)     
values(b) = runif(10000, min=0.97, max=0.99)

c <- raster(ncol=100, nrow=100)
set.seed(2)     
values(c) = runif(10000, min=0, max=1)

Now I want to create a new raster based on multiple conditions. The conditions are like

when a < 0.2, the new raster pixels should have the values from raster c;

when a > 0.5, the new raster pixels should have the values of 1.85;

when 0.2 <= a <= 0.5, the new raster pixels should have the values from raster b. I am using the following code to achieve that

# create empty copy
E <- raster(a)

E[(a < 0.2)] <- NA
E <- cover(E, c)
E[(a >= 0.2) & (a <= 0.5)] <- NA
E <- cover(E, b)
E[a > 0.5] <- 1.85

Whether it is the correct method to do that? Or there are any better alternative to this approach.

2
Does your code do what you want it to? Are you asking if this works or if this is the best method? It can be easier to tell if your code is working if you make your rasters smaller - much smaller - so you can tell by looking if the output raster is as expected. Maybe start with 4x4 rasters? - Spacedman
In your text you say "when a > 0.5, the new raster pixels should have the values of 1.85" but your code looks like E[a > 0.5] <- 0.99 - should that 0.99 be 1.85? - Spacedman
In your text you say "when 0.2 <= a <= 0.5, the new raster pixels should have the values from raster c." but your code looks like E <- cover(E, b) so did you mean "the values from raster b"? - Spacedman
Sorry for those mistakes, I have corrected the question. Is there a better approach to do that? - Bappa Das

2 Answers

1
votes

I'd do it with some ifelse replacement of the values in the raster:

Start with a blank raster:

> F = raster(a)

Then where a<0.2, take the value from c otherwise use NA:

> F[] = ifelse(a[]<0.2, c[], NA)

Then where a is over 0.5 set F to 1.85, otherwise use whatever was in F:

> F[] = ifelse(a[]>0.5, 1.85, F[])

And if a is between 0.2 and 0.5 take the value from b otherwise use whatever was in F:

> F[] = ifelse(a[]>=0.2 & a[]<=0.5, b[], F[])

This gives the same values as your E:

> all(F[]==E[])
[1] TRUE

Or you can do it by conditional replacement in the vector of raster values:

> F=raster(a)
> F[a[]<0.2] <- c[a[]<0.2]
> F[a[]>0.5] <- 1.85
> F[a[]>=0.2 & a[]<=0.5] <- b[a[]>=0.2 & a[]<=0.5]
> all(F[]==E[])
[1] TRUE

I'm not sure which is better in terms of speed or readability or flexibility. I'd wrap this all into a function anyway:

replacer = function(a, b, c,
  low_thresh=0.2, high_thresh=0.5,
  high_value=1.85){
...
}

and then benchmark it if speed is a concern. Otherwise use what you think looks best and is easiest for you to maintain.

Here's the three approaches (I've optimised the third one to avoid repeating tests):

replacer_1 = function(a, b, c,
  low_thresh=0.2, high_thresh=0.5,
  high_value=1.85){
    E <- raster(a)
    E[(a < low_thresh)] <- NA
    E <- cover(E, c)
    E[(a >= low_thresh) & (a <= high_thresh)] <- NA
    E <- cover(E, b)
    E[a > high_thresh] <- high_value
    return(E)
}

replacer_2 = function(a, b, c,
  low_thresh=0.2, high_thresh=0.5,
  high_value=1.85){
    F = raster(a)
    F[] = ifelse(a[]<0.2, c[], NA)
    F[] = ifelse(a[]>0.5, 1.85, F[])    
    F[] = ifelse(a[]>=0.2 & a[]<=0.5, b[], F[])
    return(F)
}

replacer_3 = function(a, b, c,
  low_thresh=0.2, high_thresh=0.5,
  high_value=1.85){
    low = a[]<low_thresh
    F[low] <- c[low]
    F[a[]>high_thresh] <- high_value
    mid =a[]>=low_thresh & a[]<=high_thresh 
    F[mid] <- b[mid]
    return(F)
}

Check they all return the same answer:

> E1 = replacer_1(a,b,c)
> E2 = replacer_2(a,b,c)
> E3 = replacer_3(a,b,c)
> all(E1[]==E2[])
[1] TRUE
> all(E1[]==E3[])
[1] TRUE
> all(E2[]==E3[])
[1] TRUE

Then use the microbenchmark package to compare:

> microbenchmark(replacer_1(a,b,c), replacer_2(a,b,c), replacer_3(a,b,c))
Unit: microseconds
                expr       min        lq       mean    median         uq
 replacer_1(a, b, c) 50687.567 52486.878 54089.3265 53721.770 55221.1125
 replacer_2(a, b, c)  2359.977  2507.661  2667.1080  2581.568  2642.9790
 replacer_3(a, b, c)   485.086   504.377   543.6963   542.970   565.7025

which shows me that the third method is about 100 times as fast as the first one. Enjoy all the time I've saved you!

1
votes

Please keep your example data as simple as possible making it easier to track what is going on

library(raster)
a <- b <- c <- raster(ncol=10, nrow=10)
set.seed(2)    
values(a) = sort(runif(100, min=-1, max=1))
values(b) = 3
values(c) = 5

One way to get there would be

m <- reclassify(a, c(-Inf, 0.2, 1, 0.2, Inf, NA))
x <- mask(c, m)
m <- reclassify(a, c(-Inf, 0.2, NA, 0.2, 0.5, 1, 0.5, Inf, NA))
y <- mask(b, m)
z <- init(a, 1.85)

And now

r <- merge(x, y, z)

or

r <- merge(x, y)
r <- cover(r, z)

To follow @Spacedman's lead on ifelse; you can also use the hidden .ifel method

s <- raster:::.ifel(a < 0.2, c, raster:::.ifel(a > 0.5, 1.85, b))

Or use terra::ifel --- should be faster too

library(terra)
aa <- rast(a)
bb <- rast(b)
cc <- rast(c)

ss <- ifel(aa < 0.2, cc, ifel(aa > 0.5, 1.85, bb))