Here is example code. It is fairly straight forward to adapt this code to work in a loop for processing all of your rasters. If your rasters share a common extent and resolution you can create a raster stack and loop through the bands in the stack. To create a vector containing all rasters in a directory, in a specific format, you can use "list.files" and then pass this vector to stack.

Example:

rlist=list.files(getwd(), pattern="img$", full.names=TRUE)

r <- stack(rlist)

# Add required libraries
require(raster)
require(sp)
require(rgdal)

# Set working directory, raster, in and out shapefiles
setwd("C:/test")
inshp="MyPolys"
outshp="PolyMeans"
rdata <- "Year2012.img"

# Read polygon feature class shapefile
sdata <- readOGR(dsn=getwd(), layer=inshp)

# Read raster
r <- raster(rdata)

# Extract raster values to list object
r.vals <- extract(r, sdata)

# Use list apply to calculate mean for each polygon
r.mean <- lapply(r.vals, FUN=mean)

# Join mean values to polygon data
sdata@data <- data.frame(sdata@data, m2012=r.mean)

# Write results
writeOGR(sdata, getwd(), outshp, driver="ESRI Shapefile", check_exists=TRUE, 
     overwrite_layer=TRUE)

Here is example code. It is fairly straight forward to adapt this code to work in a loop for processing all of your rasters. If your rasters share a common extent and resolution you can create a raster stack and loop through the bands in the stack. To create a vector containing all rasters in a directory, in a specific format, you can use "list.files" and then pass this vector to stack.

Example:

rlist=list.files(getwd(), pattern="img$", full.names=TRUE) 

r <- stack(rlist)   


    # Add required libraries
    require(raster)
    require(sp)
    require(rgdal)
    
    # Set working directory, raster, in and out shapefiles
    setwd("C:/test")
    inshp="MyPolys"
    outshp="PolyMeans"
    rdata <- "Year2012.img"
    
    # Read polygon feature class shapefile
    sdata <- readOGR(dsn=getwd(), layer=inshp)
    
    # Read raster
    r <- raster(rdata)
    
    # Extract raster values to list object
    r.vals <- extract(r, sdata)
    
    # Use list apply to calculate mean for each polygon
    r.mean <- lapply(r.vals, FUN=mean)
    
    # Join mean values to polygon data
    sdata@data <- data.frame(sdata@data, m2012=r.mean)
    
    # Write results
    writeOGR(sdata, getwd(), outshp, driver="ESRI Shapefile", check_exists=TRUE, 
         overwrite_layer=TRUE)