Need help extracting mean raster data into quadrats during spatial analysis

This is what I'm trying to achieve. I'm not 100% sure this is the right path to do it, but I've used a very similar process to generate count data for each quadrat which was tallied from lat/long points throughout the study area.

  • break my study area in to a series of 20m and 50m square quadrats
  • then for each quadrat I want to extract the mean/average values from a raster for both 20m and 50m - for where data isn't available simply insert '0' THIS IS WHERE I'M HAVING TROUBLE AT THE MOMENT AND SEAKING GUIDANCE
  • Eventually I want to be able to export this as a spatial raster that I can use to extract data from for each quadrat and compare to other data for similar quadrats
turtles <- read.csv("Outputs/2015/2015-nov-tsdata-1m.csv")

seagrassbio <- raster("J:/Masters/GIS/Interpolations/IDWBio/15NovIDWBio.tif")

#######################
# set class for id and date #
#######################

turtles$id <- as.factor(turtles$id)
turtles$DateTime <- as.POSIXct(turtles$DateTime, format = "%Y-%m-%d %H:%M:%S")

##############################################
# Set Spatial Points and Coordinate Reference System #
##############################################

turtles.sp <- SpatialPoints(coords=cbind(turtles$lon, turtles$lat), 
                            proj4string = CRS("+init=epsg:4326"))

proj4string(seagrassbio) <- projection(turtles.sp)

turtles.sp$id <- turtles$id    #import Turtle Tag ID

###########################################################
# Import ShapeFile of Focus Study Area and set CRS same as sat data #
###########################################################

PBanks <- readOGR("J:/Masters/GIS/ShapeFiles/Landscape/PelicanBanks.shp")

proj4string(PBanks) <- projection(turtles.sp)

############################
# Turn the Shapefile into Raster #
############################

rPBanks20 <- raster(PBanks)
rPBanks50 <- raster(PBanks)

############################################################
# Set Resolution / Cell Size of Raster, this will impact Quadrat size later #
############################################################

res(rPBanks20) <- 0.0002
res(rPBanks50) <- 0.0005

#### 0.0005 = 50m2,  0.0002 = 20m2

##################################################################
# Add ShapeFile of PBanks to above Raster, defining Rasters spatial attributes #
##################################################################

rPBanks20 <- rasterize(PBanks, rPBanks20)
rPBanks50 <- rasterize(PBanks, rPBanks50)


###########################################
# Build 'quadrats' around our predefined Cell Size #
###########################################

quadrat20 <- as(rPBanks20, "SpatialPolygons")
quadrat50 <- as(rPBanks50, "SpatialPolygons")

# orginial seagrassbio raster is over a very large area, so to minimise calculation time I'm creating a subset of a # smaller area while I'm trying to work out the code

seagrassbio.mask <- mask(seagrassbio, PBanks)

################################################
# Summarise the mean seagrassbio within each quadrat #
################################################

bioPBanks20 <- rasterize(seagrassbio.mask, rPBanks20, fun='mean', background=0)
bioPBanks50 <- rasterize(seagrassbio.mask, rPBanks50, fun='mean', background=0)

### this is where my train of thought falls apart and I receive the following error 

** Unable to find an inherited method for function ‘rasterize’ for signature ‘"RasterLayer", "RasterLayer"’ **