CODEWITHSUNDEEP
rrasterterra
apple
apple

Reputation: 500

How can I subset a raster by conditional statement in R using `terra`?

I am trying to plot only certain values from a categorical land cover raster I am working with. I have loaded it in to R using the terra package and it plots fine. However, since the original data did not come with a legend, I am trying to find out which raster value corresponds to what on the map.

Similar to the answer provided here: How to subset a raster based on grid cell values

I have tried using the following line:

> landcover
class       : SpatRaster 
dimensions  : 20057, 63988, 1  (nrow, ncol, nlyr)
resolution  : 0.0005253954, 0.0005253954  (x, y)
extent      : -135.619, -102, 59.99989, 70.53775  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : spat_n5WpgzBuVAV3Ijm.tif 
name        : CAN_LC_2015_CAL_wgs 
min value   :                   1 
max value   :                  18 

> plot(landcover[landcover == 18])
                                      
Error: cannot allocate vector of size 9.6 Gb

However, this line takes a very long time to run and produces a vector memory error. The object is 1.3 kb in the global environment and the original tif is about 300 mb.

Upvotes: 0

Views: 1770

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47591

You can use cats to find out which values correspond to which categories.

library(terra)
set.seed(0)
r <- rast(nrows=10, ncols=10)
values(r) <- sample(3, ncell(r), replace=TRUE) - 1
cls <- c("forest", "water", "urban")
levels(r) <- cls
names(r) <- "land cover"

cats(r)[[1]]
#  ID category
#1  0   forest
#2  1    water
#3  2    urban

To plot a logical (Boolean) layer for one category, you can do

plot(r == "water")

And from from the above you can see that in this case that is equivalent to

plot(r == 1)

Upvotes: 2

apple
apple

Reputation: 500

I think I found the solution to write the conditional within the plot function as below:

plot(landcover == 18)

For those looking for a reproduceable example, just load the rlogo:

s <- rast(system.file("ex/logo.tif", package="terra"))
s <- s$red
plot(s == 255)

Upvotes: 1

Related Questions