Merge branch 'develop' -- New release (MODIS functions, bfmSpatial up…

…dated)

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: bfastSpatial
 Type: Package
 Title: Utilities to monitor for change on satellite image time-series
-Version: 0.4.0
+Version: 0.5.0
 Date: 2014-04-06
 Author: Loic Dutrieux, Ben DeVries, Jan Verbesselt
 Maintainer: <[email protected]>
@@ -16,5 +16,6 @@ Depends:
     bfast,
     gdalUtils,
     stringr,
-    rgdal
+    rgdal,
+    bitops
 License: tbd

NAMESPACE

@@ -2,25 +2,28 @@
 
 export(annualSummary)
 export(areaSieve)
-export(bfmChange)
-export(bfmMagn)
 export(bfmPixel)
 export(bfmSpatial)
 export(changeMonth)
 export(cleanBrick)
+export(cleanMODIS)
 export(clumpSize)
 export(countObs)
 export(flattenBrick)
+export(getMODISinfo)
 export(getSceneinfo)
 export(mc.calc)
 export(processLandsat)
 export(processLandsatBatch)
+export(processMODISbatch)
 export(saveExtent)
 export(sr2vi)
 export(subsetRasterTS)
 export(summaryBrick)
 export(timeStack)
+export(timeStackMODIS)
 import(bfast)
+import(bitops)
 import(gdalUtils)
 import(parallel)
 import(raster)

R/bfmSpatial.R

@@ -17,19 +17,22 @@
 #' @param n See \code{\link{bfastmonitor}}
 #' @param level See \code{\link{bfastmonitor}}
 #' @param mc.cores Numeric. Number of cores to be used for the job.
-#' @param sensor Character. Optional: Limit analysis to a particular sensor. Can be one or more of \code{c("ETM+", "ETM+ SLC-on", "ETM+ SLC-off", "TM", or "OLI")}
+#' @param returnLayers Character. Result layers to be returned. Can be any combination of \code{c("breakpoint", "magnitude", "error", "history", "r.squared", "adj.r.squared", "coefficients")}. By default, \code{breakpoint}, \code{magnitude} and \code{error} are returned by the function. See \code{details} for more information.
+#' @param sensor Character. Optional: Limit analysis to one or more particular sensors. Can be any combintation of \code{c("ETM+", "ETM+ SLC-on", "ETM+ SLC-off", "TM", or "OLI")}
 #' @param ... Arguments to be passed to \code{\link{mc.calc}}
 #' 
-#' @return A rasterBrick, with 3 layers. (1) Breakpoints (time of change); (2) change magnitude; and (3) error flag (1, NA). See \code{\link{bfastmonitor}}
+#' @return A rasterBrick with layers depending on what has been supplied to \code{returnLayers}. See details for more information.
 #' 
 #' @details
 #' \code{bfmSpatial} applies \code{\link{bfastmonitor}} over a raster time series. For large raster datasets, processing times can be long. Given the number of parameters that can be set, it is recommended to first run \code{\link{bfmPixel}} over some test pixels or \code{bfmSpatial} over a small test area to gain familiarity with the time series being analyzed and to test several parameters.
 #' 
 #' Note that there is a difference between the \code{monend} argument included here and the \code{end} argument passed to \code{\link{bfastmonitor}}. Supplying a date in the format \code{c(year, Julian day)} to \code{monend} will result in the time series being trimmed \emph{before} running \code{\link{bfastmonitor}}. While this may seem identical to trimming the resulting \code{bfastmonitor} object per pixel, trimming the time series before running \code{bfastmonitor} will have an impact on the change magnitude layer, which is calculated as the median residual withint the entire monitoring period, whether or not a breakpoint is detected.
 #' 
-#' While \code{bfmSpatial} can be applied over any raster time series with a time dimension (implicit or externally supplied), an additional feature relating to the type of Landsat sensor is also included here. This feature allows the user to specify data from a particular sensor, excluding all others. This can be useful if bias in a particular sensor is of concern, and can be tested without re-making the input RasterBrick. The \code{sensor} argument accepts any combination of the following characters (also see \code{\link{getSceneinfo}}): "all" - all layers; "TM" - Landsat 5 Thematic Mapper; "ETM+" - Landsat 7 Enhanced Thematic Mapper Plus (all); "ETM+ SLC-on" - ETM+ data before failure of Scan Line Corrector; ETM+ data after failure of the Scan Line Corrector.
+#' While \code{bfmSpatial} can be applied over any raster time series with a time dimension (implicit or externally supplied), an additional feature relating to the type of Landsat sensor is also included here. This feature allows the user to specify data from a particular sensor, excluding all others. The \code{sensor} argument accepts any combination of the following characters (also see \code{\link{getSceneinfo}}): "all" - all layers; "TM" - Landsat 5 Thematic Mapper; "ETM+" - Landsat 7 Enhanced Thematic Mapper Plus (all); "ETM+ SLC-on" - ETM+ data before failure of Scan Line Corrector; ETM+ data after failure of the Scan Line Corrector.
 #' 
-#' Note that \code{names(x)} must correspond to Landsat sceneID's (see \code{\link{getSceneinfo}}), otherwise any value passed to \code{sensor} will be ignored with a warning.
+#' Note that for the \code{sensor} argument to work, \code{names(x)} must correspond to Landsat sceneID's (see \code{\link{getSceneinfo}}), otherwise any values passed to \code{sensor} will be ignored with a warning.
+#' 
+#' \code{returnLayers} can be used to specify which \code{bfasmonitor} results to return. Regardless of which parameters are assigned, the output layers will always follow the order: \code{c("breakpoint", "magnitude", "error", "history", "r.squared", "adj.r.squared", "coefficients")}. This is important if \code{mc.cores} is set to be greater than 1, since this causes the layer names in the output brick to be lost, so it is important to know which layers have been requested and in which order they will be exported. Note that if "coefficients" is included, the output will include the following: "(Intercept)" and any trend and/or harmonic coefficients depending on the values of \code{formula} and \code{order}.
 #' 
 #' @author Loic Dutrieux and Ben DeVries
 #' @import bfast
@@ -65,7 +68,7 @@
 bfmSpatial <- function(x, dates=NULL, pptype='irregular', start, monend=NULL,
                        formula = response ~ trend + harmon, order = 3, lag = NULL, slag = NULL,
                        history = c("ROC", "BP", "all"),
-                       type = "OLS-MOSUM", h = 0.25, end = 10, level = 0.05, mc.cores=1, sensor=NULL, ...) {
+                       type = "OLS-MOSUM", h = 0.25, end = 10, level = 0.05, mc.cores=1, returnLayers = c("breakpoint", "magnitude", "error"), sensor=NULL, ...) {
 
     if(is.character(x)) {
         x <- brick(x)
@@ -98,6 +101,15 @@ bfmSpatial <- function(x, dates=NULL, pptype='irregular', start, monend=NULL,
             s <- s[which(s$sensor %in% sensor), ]
         }
     }
+
+    # determine length of coefficient vector
+    # = intercept [+ trend] [+ harmoncos*order] [+ harmonsin*order]
+    coef_len <- 1 # intercept
+    modterms <- attr(terms(formula), "term.labels")
+    if("trend" %in% modterms)
+        coef_len <- coef_len + 1
+    if("harmon" %in% modterms)
+        coef_len <- coef_len + (order * 2) # sin and cos terms
 
     fun <- function(x) {
         # subset x by sensor
@@ -120,19 +132,42 @@ bfmSpatial <- function(x, dates=NULL, pptype='irregular', start, monend=NULL,
                                     type=type, h=h,
                                     end=end, level=level), silent=TRUE)
 
-            # assign NA to bkpt and magn if an error is encountered
+            # assign 1 to error and NA to all other fields if an error is encountered
             if(class(bfm) == 'try-error') {
-                res <- cbind(NA, NA, 1)
+                bkpt <- NA
+                magn <- NA
+                err <- 1
+                history <- NA
+                rsq <- NA
+                adj_rsq <- NA
+                coefficients <- rep(NA, coef_len)
             } else {
-                res <- cbind(bfm$breakpoint, bfm$magnitude, NA)
+                bkpt <- bfm$breakpoint
+                magn <- bfm$magnitude
+                err <- NA
+                history <- bfm$history[2] - bfm$history[1]
+                rsq <- summary(bfm$model)$r.squared
+                adj_rsq <- summary(bfm$model)$adj.r.squared
+                coefficients <- coef(bfm$model)
             }
         } else {
-            res <- cbind(NA, NA, NA)
+            bkpt <- NA
+            magn <- NA
+            err <- NA
+            history <- NA
+            rsq <- NA
+            adj_rsq <- NA
+            coefficients <- rep(NA, coef_len)
         }
-        names(res) <- c("breakpoint", "magnitude", "error") # these are lost if mc.cores > 1
+        res <- c(bkpt, magn, err, history, rsq, adj_rsq)
+        names(res) <- c("breakpoint", "magnitude", "error", "history", "r.squared", "adj.r.squared")
+        res <- res[which(names(res) %in% returnLayers)]
+        if("coefficients" %in% returnLayers)
+            res <- c(res, coefficients)
         return(res)
     }
 
     out <- mc.calc(x=x, fun=fun, mc.cores=mc.cores, ...)
+
     return(out)
 }

R/changeMonth.R

@@ -2,15 +2,15 @@
 #' 
 #' @description Returns one RasterLayer for each year represented by the change results. Months are represted by integers from 1 to 12.
 #' 
-#' @param change. RasterLayer output from \code{\link{bfmChange}} representing pixel breakpoints.
+#' @param change. RasterLayer representing pixel breakpoints extracted from \code{\link{bfmSpatial}} output.
 #' @param ... Additional arguments to pass to \code{\link{writeRaster}}
 #' 
 #' @return either a RasterLayer with values between 1 to 12 (representing month of change), or if multiple years are represented in the input change RasterLayer, a RasterBrick with one layer for each year, and values of 1 to 12 representing change months for each year.
 #' 
-#' @seealso \code{\link{bfmChange}}
-#' 
 #' @author Ben DeVries
 #' 
+#' @seealso \code{\link{bfmSpatial}}
+#' 
 #' @import raster
 #' @export
 

R/cleanMODIS.R

@@ -0,0 +1,73 @@
+#' Clean MODIS
+#' 
+#' @description Uses MODIS quality layer to perform cleaning over a single data layer. Quality band can be encoded bitewise or in numerics. This function needs rgdal to be configured with HDF4 driver in order to work.
+#' 
+#' @param x Character. Filename (full path) of a hdf dataset containing data and QC layers
+#' @param data_SDS Numeric. Index of the SDS containing the data in the hdf infrastructure
+#' @param QC_SDS Numeric. Index of the SDS containing the Quality control in the hdf infrastructure
+#' @param bit Logical. Is QC information provided bitwise?
+#' @param QC_val Numeric or vector of numerics. Quality control values to \textbf{keep} in the data. If \code{bit} is set to \code{TRUE} \code{QC_val} is a BYTE (hexadecimal or decimal) where each bit refers to a bit in the QC layer element (i.e.: bitpos = 0xA1 targets bits 7, 5 and 0 --- 1010 0001). When bits targetted by \code{QC_val} are activated, the corresponding observation in the data layer is filtered out.
+#' @param fill (vector of) Numeric or NULL. Fill values in the data layer to be filtered out. (important to do before reprojection with resampling method other than Nearest Neighbours)
+#' @param ... Arguments to be passed to \link{\code{writeRaster}}.
+#' 
+#' @return A rasterLayer
+#' 
+#' @author Loic Dutrieux
+#' 
+#' @seealso \code{\link{processMODISbatch}} for batcher
+#' 
+#' @examples
+#' \dontrun{
+#' library(raster)
+#' dir <- file.path(tmpDir(), 'MODIStest'); dir.create(path = dir)
+#' 
+#' fileDL <- 'http://e4ftl01.cr.usgs.gov/MOLT/MOD17A2.055/2003.07.12/MOD17A2.A2003193.h19v07.055.2011269202445.hdf'
+#' modis <- file.path(dir, basename(fileDL)); download.file(url = fileDL, destfile = modis)
+#'
+#' # Now We've just downloaded a MODIS file that is stored in a subdirectory of the raster tmp directory
+#' sprintf('These data have been acquired around %s', getMODISinfo(modis)$date)
+#' 
+#' # Clean dataset and replace fill values by NAs
+#' # Product details at https://lpdaac.usgs.gov/products/modis_products_table/mod17a2
+#' MODISclean <- cleanMODIS(x=modis, data_SDS=1, QC_SDS=3, bit=TRUE, QC_val=0x19, fill = (32761:32767))
+#' plot(MODISclean)
+#' 
+#' # In that case we did not write the result back to disk, but note that this is possible
+#' # For that use the filename= argument, and the datatype= argument (recommended)
+#' }
+#' 
+#' @import raster
+#' @import gdalUtils
+#' @import bitops
+#' @import rgdal
+#' @export
+#' 
+#' 
+
+cleanMODIS <- function(x, data_SDS, QC_SDS, bit=FALSE, QC_val, fill=NULL, ...){
+
+    if(!'HDF4' %in% gdalDrivers()$name){
+        stop('This command requires rgdal to be configured with HDF4 driver')
+    }
+
+
+    sds <- get_subdatasets(x)
+    data <- raster(readGDAL(sds[data_SDS], as.is=TRUE))
+    QC <- raster(sds[QC_SDS])
+
+    clean <- function(x, y) {
+        if(bit) {
+            x[bitAnd(y, QC_val) != 0] <- NA
+        } else {
+            x[!y %in% QC_val] <- NA
+        }
+        if(!is.null(fill)) {
+            x[x %in% fill] <- NA
+        }
+        return(x)
+    }
+
+
+    overlay(x = data, y = QC, fun = clean,...)
+
+}

R/getMODISinfo.R

@@ -0,0 +1,33 @@
+#' Get MODIS infos
+#' 
+#' @description Parses through MODIS file names and retrieve information on acquisition date. Future versions will include retrieval of sensor information, tile ID, etc... The function is vectorized.
+#' 
+#' @param x Character. MODIS file name
+#' 
+#' @return
+#' A dataframe with one column ($date) of class \code{Date}
+#' 
+#' @author Loic Dutrieux
+#' 
+#' @examples
+#' # Simple case
+#' fileName <- '/path/to/file/MOD17A2.A2008097.h11v09.005.2008142064019.hdf'
+#' getMODISinfo(fileName)
+#' 
+#' # Multile files
+#' fileNames <- c('/path/to/file/MOD17A2.A2008097.h11v09.005.2008142064019.hdf', 'MOD17A2.A2008097.h11v09.005.2008164064019.tif')
+#' getMODISinfo(fileNames)
+#' 
+#' # Get only a vector of dates
+#' getMODISinfo(fileNames)$date
+#' 
+#' 
+#' @import stringr
+#' @export
+#' 
+
+
+getMODISinfo <- function(x) {
+    date <- as.Date(str_match(pattern='(\\.A)(\\d{7})(\\.)', basename(x))[,3], format="%Y%j")
+    data.frame(date = date) # For consistency with getSceneinfo (Landsat)    
+}