library(MODISTools)
library(raster)

# Download modis data
LC <- mt_subset(product = "MCD12Q1",
                lat = 48.383662,
                lon = 2.610250,
                band = "LC_Type1",
                start = "2005-01-01",
                end = "2005-12-30",
                km_lr = 10,
                km_ab = 10,
                site_name = "testsite",
                internal = TRUE,
                progress = FALSE)

VI <- mt_subset(product = "MOD13Q1",
                lat = 48.383662,
                lon = 2.610250,
                band = "250m_16_days_NDVI",
                start = "2005-01-01",
                end = "2005-12-31",
                km_lr = 10,
                km_ab = 10,
                site_name = "testsite",
                internal = TRUE,
                progress = FALSE)

QA <- mt_subset(product = "MOD13Q1",
                lat = 48.383662,
                lon = 2.610250,
                band = "250m_16_days_pixel_reliability",
                start = "2005-01-01",
                end = "2005-12-31",
                km_lr = 10,
                km_ab = 10,
                site_name = "testsite",
                internal = TRUE,
                progress = FALSE)

mt_to_raster <- function(df = subset){

  # find unique dates for which data should exist
  dates <- unique(df$calendar_date)

  # convert scale to 1 if not available
  # should not change the numeric value of a band
  df$scale[df$scale == "Not Available"] <- 1

  # loop over all dates, format rasters and return
  r <- do.call("stack",
    lapply(dates, function(date){

     # stuff values into raster
     m <- matrix(df$value[df$calendar_date == date] *
                   as.numeric(df$scale[df$calendar_date == date]),
                df$nrows[1],
                df$ncols[1],
                byrow = TRUE)

     # convert to raster and return
     return(raster(m))
    })
   )

  # get bounding box
  bb <- mt_bbox(
    xllcorner = df$xllcorner[1],
    yllcorner = df$yllcorner[1],
    cellsize = df$cellsize[1],
    nrows = df$nrows[1],
    ncols = df$ncols[1])

  # convert to Spatial object (easier to get extent)
  bb <- as(bb, 'Spatial')

  # assign extent + projection bb to raster
  extent(r) <- extent(bb)
  projection(r) <- projection(bb)
  names(r) <- as.character(dates)

  # return the data
  return(r)
}

# convert df to raster
VI_r <- mt_to_raster(df = VI)
QA_r <- mt_to_raster(df = QA)
LC_r <- mt_to_raster(df = LC)

# upsample due to lower resolution
LC_r_2x <- resample(LC_r, QA_r, method = "ngb")

# create mask on pixel reliability flag
m <- (QA_r < 0 | QA_r > 1 )

# create land cover mask
lc_m <- !(LC_r_2x == 4 | LC_r_2x == 5)

# combine both masks
VI_m <- mask(VI_r, m, maskvalue = 1)
VI_m <- mask(VI_r, lc_m, maskvalue = 1)

# plot masked data combining pixel quality
# and showing deciduous and mixed forest pixels
# only (for the Fontainbleau forest area)
par(oma= rep(4,4))
plot(VI_r, zlim = c(0,1))
plot(VI_m, zlim = c(0,1))
original data

Screenshot from 2019-03-13 14-34-15
masked values

Screenshot from 2019-03-13 15-36-50
ropensci-archive / MODISTools

mt_to_raster function #13

original data

masked values
