Zonal statistics in GIS refers to the calculation of statistics on values of a raster within the zones of another dataset. In R the ‘extract’ function is used for this. The example below shows a zonal statistics calculation on a set of multiple rasters using a ‘for’ loop and a polygon shapefile (zones) . ‘fun=sum’ indicates that the sum of the raster cells within each polygon (zone) are being calculated, but a range of statistics can be performed. Calculations are then outputed as both a data frame and CSV file.

#import required libraries
library(maptools)
library(raster)

#list files (in this case raster TIFFs)
grids <- list.files("./path/to/data", pattern = "*.tif$")

#check the number of files in the raster list (grids)
length <- length(grids)

#read-in the polygon shapefile
poly <- readShapePoly("./path/to/data/shapefile.shp")

#create a raster stack
s <- stack(paste0("./path/to/data/", grids))

#extract raster cell count (sum) within each polygon area (poly)
for (i in 1:length(grids)){
  ex <- extract(s, poly, fun=sum, na.rm=TRUE, df=TRUE)
}

#write to a data frame
df <- data.frame(ex)

#write to a CSV file
write.csv(df, file = "./path/to/data/CSV.csv")