1

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!

lovalery
  • 4,524
  • 3
  • 14
  • 28
Mairi
  • 35
  • 7
  • Hi Mairi, without your data and code, it is almost impossible to help you as we are limited to guessing. A question on StackOverflow should always include a [minimal reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). Please consider this in your future posts on SO. That said, I think you are looking for the `raster::freq()` or `terra::freq()` function. Hope this helps. Cheers. – lovalery Mar 15 '22 at 20:32
  • I think you're looking for something like `terra::selectRange` as seen [Hijmans answer](https://stackoverflow.com/questions/69712166/extract-values-from-a-certain-layer-in-a-stack-based-on-pixel-value-of-another-r), – Chris Mar 15 '22 at 21:08
  • 1
    @lovalery Sorry, yes, good point. I'll edit my question now to make it clearer! And thank you – Mairi Mar 15 '22 at 22:01
  • @Chris Thank you for that; I've had a look, but unfortunately I don't think that's quite the right thing. I think terra::selectRange would tell me all the values in layer 2 that correspond to the cells of a particular value in layer 1, for example. What I want to do though is effectively subset all the pixels, and identify which have a '1' down for both layer 1 and 2. Have I got that right? Apologies if terra::selectRange is in fact the right thing; but if not, do you have any other suggestions? Thanks for your help! – Mairi Mar 15 '22 at 22:14
  • 1
    @lovalery I have tried using freq(), and included it in my question just now, but I'd be really grateful for any further advice! Cheers. – Mairi Mar 15 '22 at 22:28
  • The important thing in the [Hijmans answer](https://stackoverflow.com/questions/69712166/extract-values-from-a-certain-layer-in-a-stack-based-on-pixel-value-of-another-r) above is the approach to building a small, generalized [MRE](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example), as constructing helps you think about the problem at hand and how to describe it, shows others what's working or not(errors), and greatly increases your success here with any class of problem. If raster or terra, seek out Hijmans answers. He wrote it. – Chris Mar 15 '22 at 23:41

1 Answers1

1

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)

lovalery
  • 4,524
  • 3
  • 14
  • 28