Mean by month on R stacked raster

I have found on stack overflow a more generic way with the raster package using stackApply().

#get the date from the names of the layers and extract the month
indices <- format(as.Date(names(ndvi.stack), format = "X%Y.%m.%d"), format = "%m")
indices <- as.numeric(indices)

#sum layers
MonthNDVI<- stackApply(ndvi.stack, indices, fun = mean)
names(MonthNDVI) <- month.abb

## Set up color gradient with 100 values between 0.0 and 1.0
breaks <- seq(0, 1, by=0.01)
cols <- colorRampPalette(c("red", "yellow", "lightgreen"))(length(breaks)-1)
levelplot(MonthNDVI,at=breaks, col.regions=cols)

Et voilà

The most accurate solution to create monthly composites from these 16-day best value images would be to take into consideration the accompanying 'composite_day_of_the_year' scientific data set (see also MOD13A1 V006 product description). For a rather straightforward solution, please have a look at temporalComposite from MODIS and, in particular, the example included in the documentation. Using MOD13A1.006 from 2016, the code to create monthly mean value composites goes like this:

library(MODIS)

## download and extract required layers
runGdal("MOD13A1", collection = getCollection("MOD13A1", forceCheck = TRUE),
        begin = "2016001", end = "2016366", extent = "Luxembourg",
        job = "temporalComposite", SDSstring = "100000000010")

## import 16-day ndvi
ndvi <- list.files(paste0(getOption("MODIS_outDirPath"), "/temporalComposite"),
                   pattern = "NDVI.tif", full.names = TRUE)

## import corresponding composite day of the year
cdoy <- list.files(paste0(getOption("MODIS_outDirPath"), "/temporalComposite"),
                   pattern = "day_of_the_year.tif", full.names = TRUE)

## create monthly mean value composites
mmvc <- temporalComposite(ndvi, cdoy, fun = function(x) mean(x, na.rm = TRUE))

plot(mmvc[[1:4]] / 10000, zlim = c(-.1, .95))

Remember to set 'localArcPath' and 'outDirPath' before downloading images (see also ?MODISoptions).