library(raster)
## Loading required package: sp
library(rasterVis)
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
# Read data
r = brick("~/Downloads/US_NDVI_nyc_summer_mean-max-median-count_2000_2017.tif")

# Layer names
yrs = 2000:2017
funs = c("mean", "max", "median", "count")
n = rep(funs, each = length(yrs))
n = paste0(n, yrs)
n
##  [1] "mean2000"   "mean2001"   "mean2002"   "mean2003"   "mean2004"  
##  [6] "mean2005"   "mean2006"   "mean2007"   "mean2008"   "mean2009"  
## [11] "mean2010"   "mean2011"   "mean2012"   "mean2013"   "mean2014"  
## [16] "mean2015"   "mean2016"   "mean2017"   "max2000"    "max2001"   
## [21] "max2002"    "max2003"    "max2004"    "max2005"    "max2006"   
## [26] "max2007"    "max2008"    "max2009"    "max2010"    "max2011"   
## [31] "max2012"    "max2013"    "max2014"    "max2015"    "max2016"   
## [36] "max2017"    "median2000" "median2001" "median2002" "median2003"
## [41] "median2004" "median2005" "median2006" "median2007" "median2008"
## [46] "median2009" "median2010" "median2011" "median2012" "median2013"
## [51] "median2014" "median2015" "median2016" "median2017" "count2000" 
## [56] "count2001"  "count2002"  "count2003"  "count2004"  "count2005" 
## [61] "count2006"  "count2007"  "count2008"  "count2009"  "count2010" 
## [66] "count2011"  "count2012"  "count2013"  "count2014"  "count2015" 
## [71] "count2016"  "count2017"
names(r) = n

# Plot
for(i in funs) {
  
  # Subset
  j = paste0(i, 2000:2003)
  s = r[[j]]
  
  # Plot
  print(
    levelplot(
      s,
      par.settings = RdBuTheme,
      contour = FALSE,
      margin = FALSE,
      layout = c(2, 2)
    )
  )
  
}