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

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!

CodePudding user response：

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

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