Extracting values from one raster based on condition in another raster and distinguished by polygons in a vector file

Question

I have one raster file on land cover (30x30m pixels), one raster file on soil organic carbon (SOC) (250x250m pixels), and one vector file with contiguous polygons of different sizes.

What I require is a mean SOC value for each polygon. But only based on pixels that have a select land cover, here grassland.

I first made sure that the two raster files have the same extent and resolution, so I scaled the SOC map to 30x30m. I then turned the vector map into a SpatRaster (perhaps not needed), and realised that the extent is slightly different from the two raster files (the vector file when loaded into R has a 'bounding box' rather than extent, perhaps it relates to that). I also stacked the land cover and SOC files.

I am stuck here now. First, I am not sure whether it is an issue that the vector file has a different extent, and if so, how to align it with the raster files. Second, I am not sure whether I had to to convert the vector file to a SpatRaster in the first place, or if keeping it as a vector file is just fine. Third, what code could I use to now extract the mean SOC values for grassland for each polygon (i.e. a conditional if clause).

Answer 1

Here is some example data

library(terra)
pols <- vect(system.file("ex/lux.shp", package="terra"))
carbon <- rast(system.file("ex/elev.tif", package="terra"))
names(carbon) <- "carbon"
landuse <- round(carbon/100)

Here is a solution. Assuming that grassland has value "3" in landuse

grass <- subst(landuse, 3, 1, others=NA)
grass_carbon <- mask(carbon, grass)
e <- extract(grass_carbon, pols, mean, na.rm=TRUE, ID=FALSE)
data.frame(pols$NAME_2, e)
#        pols.NAME_2   carbon
#1          Clervaux 344.6667
#2          Diekirch 305.2462
#3           Redange 300.0202
#4           Vianden 295.0000
#5             Wiltz 329.2917
#6        Echternach 312.9183
#7            Remich 281.1053
#8      Grevenmacher 294.5238
#9          Capellen 321.3462
#10 Esch-sur-Alzette 300.4255
#11       Luxembourg 299.0872
#12           Mersch 301.0226

This is good enough if differences in cell sizes within polygons are relatively small. Otherwise, you could use zonal with cell sizes as weights. You could also use that if the land cover classes were expressed as fractions of the grid cells.

To plot the results, you can do

pols$carbon <- e
plot(pols, "carbon")