R MODIS MOD13Q1 - How to handle data quality issues?

Question

I'm running into some issues with 250m 16-day MOD13Q1 NDVI data using the MODIS package in R.

Looking at the following raster plot, there are certain days where the NDVI does not match the true surface characteristics. Specifically ndvi.2017.11.01 and ndvi.2018.01.17. Even ndvi.2018.03.06 and ndvi.2018.02.18 are questionable. Is it possible to programmatically exclude these layers? Is there something to gain by using the following sds layer: MODIS_Grid_16DAY_250m_500m_VI:250m 16 days VI Quality?

Below is the code to produce the above plot:

library(raster)
library(MODIS)
library(parallel)
library(pbapply)
library(magrittr)
library(dplyr)
library(stringr)

# Final Study Area Boundary
study_area <- readRDS(gzcon(url('http://web.pdx.edu/~porlando/study_area.RDS')))  
# Study CRS
epsg_26910 <- "+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs "

# Download Data
runGdal(product = "MOD13Q1", collection = "006"
           ,tileH = 9, tileV = 4
           ,begin = str_replace_all("2017-09-01", "-", ".")
           ,end = str_replace_all("2018-09-30", "-", ".")
           ,overwrite = TRUE
           ,job="temporalComposite"
           #,SDSstring="" # not sure about which SDS string to use...
)


# Process data
ndvi_path <- "./MODIS/MOD13Q1.006" # 250 m

ndvi_files <- list.files(path = ndvi_path, pattern = "\\.hdf$"
                     ,all.files = FALSE, full.names = TRUE
                     ,recursive = TRUE
                     ,ignore.case = FALSE)


processNDVI <- function(file_path, study_area = study_area, proj = epsg_26910) {

  cat("\n")

  # extract date from file path
  date <- stringr::str_extract(file_path, '[0-9][0-9][0-9][0-9].[0-9][0-9].[0-9][0-9]')
  date_clean <- stringr::str_replace_all(date, "\\.", "-")

  # coerce 16 day average to daily composite?
  dates <- seq.Date(from = as.Date(date_clean)
                    #,to = as.Date(date_clean) + lubridate::days(15)
                    ,to = as.Date(date_clean)
                    ,by = "1 day")

  sds <- get_subdatasets(file_path)

  ndvi <- sds[grepl("16DAY_250m_500m_VI:250m 16 days NDVI", sds)] %>% readGDAL %>% raster
  evi <- sds[grepl("16DAY_250m_500m_VI:250m 16 days EVI", sds)] %>% readGDAL %>% raster
  # VI Quality may be useful?
  vi_quality <- sds[grepl("16DAY_250m_500m_VI:250m 16 days VI Quality", sds)] %>% readGDAL %>% raster

  # combine NDVI and EVI into a single stack
  r <- stack(ndvi, evi)
  names(r) <- c(paste0("ndvi.", date), paste0("evi.", date))

  r <- raster::projectRaster(from = r, res = 250, method = 'ngb', crs = proj)
  study_area <- spTransform(study_area, CRSobj = proj)

  if(identical(crs(study_area), crs(r))) {

    m <- mask(r, study_area)
    cr <- crop(m, study_area)

    pblapply(dates, function(x) {
      pblapply(1:nlayers(cr), function(y) {
        layer_name <- names(cr[[y]])
        var <- gsub("\\..*$", "", layer_name)
        date <- gsub("-", "\\.", x)
        writeRaster(cr[[y]]
                    ,filename = paste0("./output/", var, "/", var, ".", date, ".tif")
                    ,format = "GTiff"
                    ,overwrite = TRUE)
      })
    })
}

pblapply(ndvi_files, function(x) {
  processNDVI(file_path = x
              ,study_area = study_area
              )
  }
  ,cl = detectCores()-1
    )

ndvi <- stack(list.files(path = "./output/ndvi/", pattern = ".tif$", full.names = T))
plot(ndvi)