R - update raster values which fall inside a polygon

Question

I'm working with environmental layers in R using terra. I want to modify the values of a raster that fall inside a polygon. I know the classify function that can update all the values of a raster, but I don't know how to apply this only to the area of the polygon. I also found a nice and neat way to change all values inside a polygon, like this:

myraster[mypolygon] <- 7

but I just need to update one particular value, not all, for example, change all cells containg the value 21 for the value 7.

I unsuccessfully tried:

myraster[mypoligon][myraster[mypoligon] == 21] <- 7

getting the message:

Error: [`[`] value is too long

As an alternative I'm thinking of getting the raster that intersects the polygon, then classify the values, and merge back this small raster with the original, but I wonder if a quicker process exists.

PS. I posted a reproducible example below as an answer.

Answer 1

Here is one way to get there using cover.

Example data

library(terra)
p <- vect(system.file("ex/lux.shp", package="terra"))[3,]
r <- rast(system.file("ex/elev.tif", package="terra"))

Let's update the area of r that is covered by p to 400, but only if that area has values > 400.

x <- mask(r, p)
y <- clamp(x, -Inf, 400)
z <- cover(y, r)

If you wanted to change a particular value, let's say 482, to 1000 you could do

yy <- classify(x, cbind(482, 1000))
zz <- cover(yy, r)

And a range of values

yyy <- classify(x, cbind(480, 500, 1000))
zzz <- cover(yyy, r)

Answer 2

Not exactly a one-liner, and probably there are more elgant approaches, but this gets the job done (basically chaining terra's mask, classify and merge):

library(terra)
#> terra 1.7.71

# sample data
r <- system.file("ex/elev.tif", package = "terra") |> terra::rast()
v <- system.file("ex/lux.shp", package = "terra") |> terra::vect()

# subset raster
mask <- v[12, ]
r_sub <- terra::mask(r, mask, touches = FALSE)

# reclassify, here: 350 -> 1000
m <- c(350, 350, 1000) |> matrix(ncol = 3, byrow = TRUE)

r_rcl <- terra::classify(r_sub, rcl = m, include.lowest = TRUE)

# merge again, pay attention to the order of the arguments 
r2 <- terra::merge(r_rcl, r)
r2
#> class       : SpatRaster 
#> dimensions  : 90, 95, 1  (nrow, ncol, nlyr)
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> varname     : elev 
#> name        : elevation 
#> min value   :       141 
#> max value   :      1000

terra::plot(r2)
terra::plot(mask, add = TRUE)

^{Created on 2024-03-25 with reprex v2.1.0}

Edit: Oh, should have read your complete question. Just noticed I proposed exactly the alternative approach you described :(

Answer 3

I found myself a possible solution following the idea posted here: updating raster values with a vector or matrix of values based on polygon field in R terra.

The thing is accessing the values directly, eg.:

myraster[cells_id] <- 7

Here is a reproducible example:

library(terra)
#terra 1.7.55

# sample data
r <- rast(system.file("ex/elev.tif", package = "terra"))
v <- vect(system.file("ex/lux.shp", package = "terra"))[6,]
plot(r)
polys(v)

# Get a raster with the shape of the polygon and the same resolution of r
raster_polygon <- rasterize(v, r)

# Get the cell ids and values
cells_id <- as.data.frame(raster_polygon, cells = TRUE)
cells_id <- cells_id$cell
cell_values <- r[cells_id]

# Change the values keeping then in a vector
cell_values[cell_values == 300] <- 500

# Assign those values to the raster
r[cells_id] <- cell_values

plot(r)
polys(v)

Here are the maps, before and after the value update (notice the green dots):