Stack 2 decades EVI MODIS time series data into a single HDF file

Hi,

I am working with 16-day interval time-series MODIS data for 2 decades, available at https://lpdaac.usgs.gov/products/mod13a2v006/. The data is in hdf format. I would like to find out whether it is possible to strictly work with hdf files rather than convert the datasets to tif files. This would save time and disc space.

My main question is: Is it possible to stack EVI or NDVI subdatasets from the MODIS data into a single file (i.e stacked EVI for the whole time series as an hdf)?

Below is a script for a single input file but I would like to do it for the entire time series maybe as a loop or something, as long as it is automated...

Thanks in advance
Edward

library(sp)
library(rgdal)
library(raster)
library(gdalUtils)

setwd("C:/Users/Downloads/MOD/MOD")

Single input file

this_file <- "2000.05.24/MOD13A2.A2000145.h20v11.006.2015137094705.hdf"

List subdatasets

subdatasets <- get_subdatasets(this_file)

Print subdataset names

print(subdatasets)

Subdataset 1 is NDVI, subdataset 2 is EVI

Load each dataset

NDVI <- raster(subdatasets[1])
EVI <- raster(subdatasets[2])

Datasets are scaled to integer values. Multiply by 0.0001 to get back to index values.

NDVI_cal <- NDVI * 0.0001 * 0.0001
EVI_cal <- EVI * 0.0001 * 0.001

Display an EVI image

plot(EVI_cal)