0

I am using the {terra} package of R-Stats. I have some quite complicated conditional maths to run on a set of rasters; a different formula depending on the value in the cell. I have tried to make a simplified version below to use as an example. But it does not work.

Currently I get the message " [lapp] I do not like 'fun' :( "

# Make an INT1U raster, where the value 255 signifies NA
r <- rast(ncol = 2, nrow = 2, vals = c(1,2,3,255))
r <- writeRaster(r, 
                 filename = "/databricks/driver/r.tif",
                 overwrite = TRUE,
                 wopt = list(datatype = "INT1U"))

# Make another raster with no NA values
s <- rast(ncol = 2, nrow = 2, vals = c(1,2,3,4))
s <- writeRaster(s, 
                 filename = "/databricks/driver/s.tif",
                 overwrite = TRUE,
                 wopt = list(datatype = "INT1U"))

# Now I want to do conditional functions on the combination of cells of each layer
t <- lapp(c(r,s),
          fun = function(a,b){
            
            if (is.na(a)) { 55 } 
            if (a > 10 & a < 20) {a + b}
            if (a > 20 & a < 30) {a / b}
            if (a == 50) { 99 }
            
          },
          filename = "/databricks/driver/t.tif",
          overwrite = TRUE,
          wopt = list(datatype = "INT1U"))

 [lapp] I do not like 'fun' :(

2 Answers 2

4

As @Spacedman says, a vectorized solution should work this function should give you the expected result.

t <- lapp(c(r,s),
          fun = function(a,b){
            a[is.na(a)] <- 55
            a[a>10 & a <20] <- a[a>10 & a <20] + b[a>10 & a <20]
            a[a>20 & a <30] <- a[a>20 & a <30]/b[a>20 & a <30]
            a[a==50] <- 99
            return(a)
            
          },
          filename = "t.tif",
          overwrite = TRUE,
          wopt = list(datatype = "INT1U"))
5
  • Wow. Nice one Hugh. I wouldn't have considered writing the function in that style. Yes, it seems to give the correct response in my minimal example. Thanks! I'll now try and do my 'proper' work and see if it works there also. May 12 at 14:07
  • I wonder if this is going to struggle with memory issues on larger rasters though? Hmm. Will have to wait and see I guess. May 12 at 14:19
  • No worries - be interesting to hear if this scales at all... I guess tiling might be needed... I don't often use this lapp function so thanks for the chance to dig into it a little! May 12 at 14:34
  • I've already split my massive raster (5m of Italy) into tiles, so touch wood it will run ok. Getting a new error on my 'real' data unfortunately. See EDIT1 above. Any thoughts please? May 12 at 16:50
  • Ignore me. It was a typo in my code. May 12 at 16:57
2

This from the docs:

 Before you use the function, test it to make sure that it is
 vectorized. That is, it should work for vectors longer than one,
 not only for single numbers. 

Let's test:

> fun(c(1,2,3),c(4,5,6))
Warning messages:
1: In if (is.na(a)) { :
  the condition has length > 1 and only the first element will be used
2: In if (a > 10 & a < 20) { :
  the condition has length > 1 and only the first element will be used
3: In if (a > 20 & a < 30) { :
  the condition has length > 1 and only the first element will be used
4: In if (a == 50) { :
  the condition has length > 1 and only the first element will be used

You have to rewrite the code so that it deals with a and b longer than 1.

I think you'll have to write this as a series of conditional subsets, or use Vectorize to make this function work (but you'll also have to explicitly use return(55) inside the if statements.

Here's the Vectorize solution (someone's already shown the chain-of-conditions solution, which will be faster) just for "fun"...

fun_scalar = function(a,b){
    if (is.na(a)) { return(55) } 
    if (a > 10 & a < 20) {return(a + b)}
    if (a > 20 & a < 30) {return(a / b)}
    if (a == 50) { return(99) }
    return(a) #?
}

fun = Vectorize(fun_scalar)

The conditional version is:

fun_cond = function(a,b){
    a[is.na(a)] <- 55
    a[a>10 & a <20] <- a[a>10 & a <20] + b[a>10 & a <20]
    a[a>20 & a <30] <- a[a>20 & a <30]/b[a>20 & a <30]
    a[a==50] <- 99
    return(a)    
}

And now:

> a = c(NA,10,15,20,25,30); b=c(10,15,20,25,25,25)
> all(fun(a,b) == fun_cond(a,b))
[1] TRUE
1
  • Thanks very much @Spacedman ! Every day is a school day :-) I've implemented Hugh's code now, and as you reckon it'll be faster I'll stick with it. But good to understand how to use Vectorise. May 12 at 17:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.