smoothing raster layer stack using R

Question

I'm trying to do smoothing for .tif stack layer of 331 bands. I used this code in R :

setwd("D:\\")    

raster_dir <-   "D:\\Modis VI quality\\"

VI <- brick(paste(raster_dir,"laikipia_VI_quality.tif", sep=""))

timeInfo <- orgTime(VI, nDays=16, begin="2001001", end="2015145")

system.time(whittaker.raster(vi,timeInfo=timeInfo,lambda=5000,overwrite="true"))

But, it is not working and giving,

Warning messages: 1: In min.default(numeric(0), na.rm = FALSE) : no non-missing arguments to min; returning Inf 2: In max.default(numeric(0), na.rm = FALSE) : no non-missing arguments to max; returning -Inf Errors: Error: dir.create(path, recursive = TRUE, showWarnings = showWarnings) is not TRUE

Any advice or help writing the correct code please?

Answer 1

A possible reason might be that vi (which you pass on to whittaker.raster) does not exist. Please provide the error message in order to clarify this issue.

Note that MODIS has been moved to GitHub a while ago. You can easily install the latest versions via

• (stable version)
  devtools::install_github("MatMatt/MODIS")

• (development version)
  devtools::install_github("MatMatt/MODIS", ref = "develop").

However, I noticed that there is something wrong with the most recent version of whittaker.raster. In my case, it finishes without throwing an error, but the smoothed rasters hold integer values only (i.e., 0, 1, 2, ...). A possible, yet quite inconvenient workaround until this issue is solved would be to temporarily install an older version of MODIS (you can download MODIS_0.10-18.tar.gz from my Dropbox) and employ the whittaker.raster algorithm included therein.

Here is some sample code that works on my machine

> sessionInfo()
R version 3.3.0 (2016-05-03)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS
(...)

using the 'develop' version from above and, temporarily, MODIS_0.10-18. The Kenyan shapefile I am using for demonstration purposes (since you did not provide any sample data) originates from GADM and can easily be downloaded via getData from raster. I omitted the usage of orgTime since (i) the required information is automatically extracted from the supplied 'vi' argument during whittaker.raster and (ii) I did not want to modify the default 'pillow' size (see ?orgTime).

## install and load 'MODIS' package
devtools::install_github("MatMatt/MODIS", ref = "develop")
library(MODIS)

## set 'MODIS' options
if (!dir.exists("data")) 
  dir.create("data")

MODISoptions(localArcPath = "data/MODIS_ARC", outProj = "+init=epsg:21037",
             outDirPath = "data/MODIS_ARC/PROCESSED", 
             MODISserverOrder = c("LPDAAC", "LAADS"))

## reference extent (kenya)
shp <- getData(country = "KEN", level = 0, path = "data")
shp <- spTransform(shp, CRS = CRS("+init=epsg:21037"))

## download tiles
collection <- getCollection("MYD13A2", forceCheck = TRUE)

runGdal("MYD13A2", collection = collection, job = "MYD13A2.006", 
        tileH = 21, tileV = 9, SDSstring = "101000000000", 
        begin = "2004001", end = "2005365")

## pre-processing (import, crop, mask)
lst <- lapply(c("NDVI", "VI_Quality"), function(i) {
  # list and import available files
  fls <- list.files(paste0(getOption("MODIS_outDirPath"), "/MYD13A2.006"), 
                    pattern = paste0(i, ".tif$"), full.names = TRUE)

  rst <- stack(fls)

  # crop reference extent and, if required, apply scaling factor
  rst_crp <- crop(rst, shp)     

  if (i == "NDVI") 
    rst_crp <- rst_crp * 0.0001

  # return processed images
  return(rst_crp)
})

## install and load 'MODIS' version 0.10-18 (assuming that the older version of 
## 'MODIS' is located inside your current working directory)
detach("package:MODIS", unload = TRUE)
install.packages("MODIS_0.10-18.tar.gz", repos = NULL, type = "source")
library(MODIS)

## apply whittaker smoother
rst_wht <- whittaker.raster(vi = lst[[1]], w = lst[[2]], overwrite = TRUE)

---