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. I use ArcGIS 10.5.1. The two raster files and the vector file should have the same extent, because I used the outline surrounding all polygons in the vector file to clip the raster files.

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 initially tried zonal statistics in ArcGIS but had issues with loading the SOC Raster as a raster input (likely because it is a floating point rather than signed integer). So I decided to move everything into R environment.

In R, I 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 if moving from ArcGIS into R was required, or whether my question could be tackled in ArcGIS directly. Second, if I were to go ahead with R, what code could I use to now extract the mean SOC values for grassland for each polygon.

Answer 1

You can use a mask to limit the mean to a specific landcover type.

First, add our libraries and create example data. The polys are the state boundaries of North Carolina, lulc is a random raster with (1,2,3,4) values and soc is a random raster (0.5-3.0). I believe that this somewhat emulates your data.

library(sf)
library(terra)

polys <- st_read(system.file("shape/nc.shp", package="sf"))
  polys <- st_transform(polys, crs="ESRI:102008")

lulc <- rast(ext(polys), resolution = 1000)
  lulc[] <- sample(1:4, ncell(lulc), replace=TRUE) 
soc <- rast(ext(polys), resolution = 1000)
  soc[] <- sample(seq(0.5,3,0.01), ncell(soc), replace=TRUE) 

plot(c(lulc, soc))

Assuming that the grasslands class is 3 we turn everything else into nodata (NA) and then mask the soc raster so that the only soc cells that are not NA correspond with the grassland class.

grass <- ifel(lulc == 3, 1, NA)
  grass.soc <- mask(soc, grass)
    plot(grass.soc, col="black",legend=FALSE)

Now we can extract values, for each polygon, for the soc.grass and soc rasters and take the mean. We then assign the results back to the polygons and plot the results.

mean.soc <- extract(c(soc,grass.soc), vect(polys), 
           fun=function(x) {mean(x,na.rm=TRUE)})
  names(mean.soc) <- c("ID","SOC","GRASS_SOC")  
    head(mean.soc)
    
polys$soc <- mean.soc[,2]
polys$grass.soc <- mean.soc[,3]
  plot(polys[c("soc","grass.soc")])

Keep in mind that these examples are random rasters so, the means should be somewhat equally distributed. An example of code efficiency you can mask the data on the fly,

extract(mask(soc, grass), vect(polys), 
  fun=function(x) {mean(x,na.rm=TRUE)})