Search code examples
rrasterterra

How can I calculate how many pixels share the same value in two raster layers?


I have made a stacked raster of two layers with the same extent, which contain information about the presence of two different things across the area. Both are binary; the values are either 1 (presence) or 0 (absence) for each pixel.

This is the raster's description:

> stacked
class       : SpatRaster 
dimensions  : 166, 1622, 2  (nrow, ncol, nlyr)
resolution  : 0.1666667, 0.1666667  (x, y)
extent      : -131.5, 138.8333, 36.33333, 64  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : memory  
              memory  
names       : lyr1, lyr1 
min values  :    0,    0 
max values  :    1,    1 

How can I use conditions to extract the number of pixels where value = 1 in both layers?

I tried a few things when the two areas were separate rasters, but had no success with that, so I used extend() to enlarge both their extents so that they could share a common area. I think this means there will be NAs in there as well as 1s and 0s.

I have tried the following:

> freq(stacked, value = 1)
     layer value count
[1,]     1     1 12243
[2,]     2     1 14804

but this just counts how many pixels have a value of 1 in each of the layers separately, whereas what I need is the number of pixels in which the values for BOTH layers =1, so they match.

Many thanks for any tips!


Solution

  • If I fully understand your request, please find below one possible solution.

    Reprex

    library(terra)
    
    # Build a dummy Spatraster with two layers containing only 1 and 0
    set.seed(0) # for reproducibility
    r <- rast(nrows=10, ncols=10, nlyrs=2)
    values(r) <- sample(c(1,0), 200, replace = TRUE)
    
    r
    #> class       : SpatRaster 
    #> dimensions  : 10, 10, 2  (nrow, ncol, nlyr)
    #> resolution  : 36, 18  (x, y)
    #> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    #> coord. ref. : lon/lat WGS 84 
    #> source      : memory 
    #> names       : lyr.1, lyr.2 
    #> min values  :     0,     0 
    #> max values  :     1,     1
    
    
    # If you use terra::freq(), you get the number of cells with 0 and 1 for each layer
    terra::freq(r)
    #>      layer value count
    #> [1,]     1     0    52
    #> [2,]     1     1    48
    #> [3,]     2     0    47
    #> [4,]     2     1    53
    
    # So, one approach is to sum the two layers and retrieve the number of cells with 2
    # (here 24 cells have the value 1 in the two layers)
    
    terra::freq(sum(r))
    #>      layer value count
    #> [1,]     1     0    23
    #> [2,]     1     1    53
    #> [3,]     1     2    24
    
    # OR more directly:
    terra::freq(sum(r), value = 2)[, 'count']
    #> count 
    #>    24
    

    Created on 2022-03-15 by the reprex package (v2.0.1)