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)