Search code examples
rrasterr-raster

How to unify the number of non-NA cells in raster stack?


I have raster files that have the same resolution and extent but differ in the number of NA. I want to unify the number of NA between all of them. Is it possible to do it by considering a cell as non-NA if it's not NA in all the raster files?

Here an example :

   library(raster)
library(terra)

f <- system.file("external/test.grd", package="raster")
r1 <- raster(f)
r2 <- calc(r1, fun=function(x){ x[x < 500] <- NA; return(x)} )
r1 <- calc(r1, fun=function(x){ x[x > 1200] <- NA; return(x)} )
raste <-  rast(r1)
rNA <- terra:: global(!(is.na(raste)), sum, na.rm=TRUE)
print(paste0("Non-NA of r1", rNA))
raste <-  rast(r2)
rNA <- terra:: global(!(is.na(raste)), sum, na.rm=TRUE)
print(paste0("Non-NA of r2", rNA))

I want both r1 and r2 to have the same number of non-NA cells. I have more than two rasters, so I wonder if I can do it for a large number of files.


Solution

  • It can be a bit confusing to use raster and terra together, so I will just use terra (but you can do the same with raster, using stack in stead of c and cellStats in stead of global.

    Your example data

    library(terra)
    f <- system.file("external/test.grd", package="raster")
    r <- rast(f)
    r1 <- clamp(r, upper=1200, values=FALSE)
    r2 <- clamp(r, lower=500, values=FALSE)
    
    global(!(is.na(r1)), sum)
    #       sum
    #lyr.1 3145
    global(!(is.na(r2)), sum)
    #      sum
    #lyr.1 802
    

    Solution:

    r <- c(r1, r2)
    names(r) <- c("r1", "r2")
    
    m <- any(is.na(r))
    x <- mask(r, m, maskvalue=1)
    
    global(!(is.na(x)), sum, na.rm=TRUE)
    #   sum
    #r1 769
    #r2 769
    

    I like the use of any(is.na()) because it makes clear what the intent is.

    But you could combine the layers in one of many other ways. As long as you do not use na.rm=TRUE the cells with an NA in one of the layers will be NA in the output. For example with sum, diff, prod, mean or app.

    m <- sum(r)
    x <- mask(r, m)
    global(!(is.na(x)), sum, na.rm=TRUE)