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.

1 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.


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")  
polys$soc <- mean.soc[,2]
polys$grass.soc <- mean.soc[,3]

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)})  
  • I notice that lulc and soc are SpatRaster files. Currently, my land cover and soc files are RasterLayers. Is it necessary to convert to SpatRaster? I tried this using rast(), but I get names for min and max value (i.e. land cover type names), rather than numbers. Does the ifel function also work with names? If not, is there a way to get the number associated with the land cover type name?
    – tabtimm
    Commented Mar 6, 2023 at 23:49
  • The RasterLayer attributes have 'from:' and 'to:' as rows and 'ID' (numbers) and 'lulc' (names) as columns.
    – tabtimm
    Commented Mar 6, 2023 at 23:53
  • The raster library is being phased out so, you should start using terra, terra::rast is used to read rasters. Raster data is always numeric! It may have attributes associated but, the actual cell values are numeric so, you must know what number is associate with a given attribute. You can use levels to see the attribute table of the raster and the associated cell values represented in the "ID" column. Commented Mar 7, 2023 at 0:26
  • Ok, thanks, understand. I got this error (Error: [mask] extents do not match), despite having clipped both rasters in ArcGIS earlier based on the same shape file outline. Is there a way to align the extents in R, or do I have to go back into ArcGIS for that?
    – tabtimm
    Commented Mar 7, 2023 at 0:39
  • Due to rounding and origin offsets this tends to not be as straight forward as it seems in ESRI products. You can use mask(crop(soc, ext(grass)), grass) to fix the extent within the mask call. Commented Mar 7, 2023 at 2:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.