I work with MODIS NDVI rasters in 2016. I have 23 rasters stacked in one object. I have 2 raster by month. I would like the average by months and conserve a raster for each month.

ndvi.stack <- stack(result)

# attribute name for each raster stacked
idx <- seq(as.Date('2016-01-17'), as.Date('2017-01-03'), '16 day')
names(ndvi.stack) <- idx

#[1] 302 268  23
## 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(ndvi.stack,at=breaks, col.regions=cols, main="NDVI 2016")

NDVI for each 16 days

I would like to get something like that rasterViz package visualization


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à enter image description here


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:


## 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).

  • I have ever downloaded 5 years of modis images ... I would like to use it! But thank you! – delaye Apr 19 '17 at 17:26
  • 2
    I've added data download and layer extraction via runGdal for the sake of reproducibility only. As long as you have the initial .hdf files available, simply extract all 'composite_day_of_the_year' layers and put them in temporalComposite() together with the raw NDVI. – fdetsch Apr 19 '17 at 17:36
  • Really great to see this! Are weekly/monthly means possible with the MODIS package for other products like MCD19A2 that do not provide a daily composite? – philiporlando Sep 17 '18 at 22:18
  • 3
    Unfortunately, no. For daily products like MCD19A2, the raster::stackApply() based solution provided by @delaye is totally sufficient for creating weekly, monthly, etc. composites. – fdetsch Sep 18 '18 at 6:41
  • Thank you for the info! I'm making some progress with raster::stackApply(). – philiporlando Sep 18 '18 at 19:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.