Below you’d find a simple R code that calculates the average value of a raster within polygons:
# load the required packages
library(raster)
library(sp)
# read your raster and polygon datasets
raster_data <- raster("directory/raster.tif")
polygon_data <- readOGR("directory/polygons.shp")
# crop your raster to the extent of your polygon dataset
raster_data_crop <- crop(raster_data, extent(polygon_data))
# convert the polygon dataset to the same coordinate reference system as the raster data
proj4string(polygon_data) <- proj4string(raster_data_crop)
# extract values from raster within each polygon
raster_values <- extract(raster_data_crop, polygon_data)
# calculate the average value of the raster within each polygon
polygon_data$average_value <- sapply(raster_values, mean)
# view the results
head(polygon_data)