Merge branch 'develop'

DESCRIPTION

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

NAMESPACE

@@ -3,7 +3,9 @@
 export(annualSummary)
 export(areaSieve)
 export(bfmPixel)
+export(bfmSpOver)
 export(bfmSpatial)
+export(bfmZoo)
 export(changeMonth)
 export(cleanBrick)
 export(cleanMODIS)
@@ -22,10 +24,13 @@ export(subsetRasterTS)
 export(summaryBrick)
 export(timeStack)
 export(timeStackMODIS)
+export(zooExtract)
 import(bfast)
 import(bitops)
 import(gdalUtils)
 import(parallel)
 import(raster)
 import(rgdal)
+import(sp)
 import(stringr)
+import(zoo)

R/LandsatVIs.R

@@ -59,70 +59,72 @@
                 fun=fun))
 }
 
-# Tasseled Cap Components -------------------------------------------------
+# Tasseled Cap Components ------------------------------------------------- 
 
-.tcbright <- function(sensor) {
-    ind <- c('band1','band2','band3','band4','band5','band7') 
-    # make compatible with getSceneinfo() output
-    if(sensor %in% c("ETM+", "ETM+ SLC-on", "ETM+ SLC-off"))
-        sensor <- 7
-    if(sensor == "TM")
-        sensor <- 5
+.tcbright <- function(sensor) { 
+    ind <- c('band1','band2','band3','band4','band5','band7')  
+    # make compatible with getSceneinfo() output 
+    if(sensor %in% c("ETM+", "ETM+ SLC-on", "ETM+ SLC-off")) 
+        sensor <- 7 
+    if(sensor == "TM") 
+        sensor <- 5 
 
-    if(sensor == 5) {
-        tc_coef <- c(0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863)
-    } else if (sensor == 7) {
+    if(sensor == 5) { 
+        tc_coef <- c(0.2043, 0.4158, 0.5524, 0.5741, 0.3124, 0.2303) 
+    } else if (sensor == 7) { 
         tc_coef <- c(0.3561, 0.3972, 0.3904, 0.6966, 0.2286, 0.1596)
-    }
+    } 
 
-    fun <- function(x1, x2, x3, x4, x5, x7) {
-        tcbright <- sum(c(x1, x2, x3, x4, x5, x7) * tc_coef)
-    }
+    fun <- function(x1, x2, x3, x4, x5, x7) { 
+        tcbright <- sum(c(x1, x2, x3, x4, x5, x7) * tc_coef) 
+    } 
 
-    return(list(ind=ind,
-                fun=fun))
-}
+    return(list(ind=ind, 
+                fun=fun)) 
+} 
+
 
-.tcgreen <- function(sensor) {
-    ind <- c('band1','band2','band3','band4','band5','band7')
-    # make compatible with getSceneinfo() output
-    if(sensor %in% c("ETM+", "ETM+ SLC-on", "ETM+ SLC-off"))
-        sensor <- 7
-    if(sensor == "TM")
-        sensor <- 5
+.tcgreen <- function(sensor) { 
+    ind <- c('band1','band2','band3','band4','band5','band7') 
+    # make compatible with getSceneinfo() output 
+    if(sensor %in% c("ETM+", "ETM+ SLC-on", "ETM+ SLC-off")) 
+        sensor <- 7 
+    if(sensor == "TM") 
+        sensor <- 5 
 
-    if(sensor == 5) {
-        tc_coef <- c(-0.2848, -0.2435, -0.5436, 0.7243, 0.0840, -0.1800)
-    } else if (sensor == 7) {
-        tc_coef <- c(-0.3344, -0.3544, -0.4556, 0.6966, -0.0242, -0.263)
-    }
+    if(sensor == 5) { 
+        tc_coef <- c(-0.1603, -0.2819, -0.4934,  0.7940, 0.0002, -0.1446) 
+    } else if (sensor == 7) { 
+        tc_coef <- c(-0.3344, -0.3544, -0.4556,  0.6966, -0.0242,-0.2630)
+    } 
 
-    fun <- function(x1, x2, x3, x4, x5, x7) {
-        tcbright <- sum(c(x1, x2, x3, x4, x5, x7) * tc_coef)
-    }
+    fun <- function(x1, x2, x3, x4, x5, x7) { 
+        tcbright <- sum(c(x1, x2, x3, x4, x5, x7) * tc_coef) 
+    } 
 
-    return(list(ind=ind,
-                fun=fun))
-}
+    return(list(ind=ind, 
+                fun=fun)) 
+} 
 
-.tcwet <- function(sensor) {
-    ind <- c('band1','band2','band3','band4','band5','band7')
-    # make compatible with getSceneinfo() output
-    if(sensor %in% c("ETM+", "ETM+ SLC-on", "ETM+ SLC-off"))
-        sensor <- 7
-    if(sensor == "TM")
-        sensor <- 5
+
+.tcwet <- function(sensor) { 
+    ind <- c('band1','band2','band3','band4','band5','band7') 
+    # make compatible with getSceneinfo() output 
+    if(sensor %in% c("ETM+", "ETM+ SLC-on", "ETM+ SLC-off")) 
+        sensor <- 7 
+    if(sensor == "TM") 
+        sensor <- 5 
 
-    if(sensor == 5) {
-        tc_coef <- c(0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572)
-    } else if (sensor == 7) {
-        tc_coef <- c(0.2626, 0.2141, 0.0926, 0.0656, -0.7629, -0.5388)
-    }
+    if(sensor == 5) { 
+        tc_coef <- c(0.0315,  0.2021,  0.3102,  0.1594, 0.6806, -0.6109) 
+    } else if (sensor == 7) { 
+        tc_coef <- c(0.2626, 0.2141, 0.0926, 0.0656, -0.7629, -0.5388) 
+    } 
 
-    fun <- function(x1, x2, x3, x4, x5, x7) {
-        tcbright <- sum(c(x1, x2, x3, x4, x5, x7) * tc_coef)
-    }
+    fun <- function(x1, x2, x3, x4, x5, x7) { 
+        tcbright <- sum(c(x1, x2, x3, x4, x5, x7) * tc_coef) 
+    } 
 
-    return(list(ind=ind,
-                fun=fun))
-}
+    return(list(ind=ind, 
+                fun=fun)) 
+} 

R/bfmSpOver.R

@@ -0,0 +1,83 @@
+#' Runs bfastmonitor for a spatial subset with aggregation
+#' 
+#' @description Runs \link{bfastmonitor} on a rasterBrick object for a set of locations, determined by an object of class \link{Spatial-class}.
+#' 
+#' @param x A rasterBrick or rasterStack, ideally with time written to the z dimension. In case time is not written to the z dimension, the \code{dates=} argument has to be supplied (see \link{zooExtract})
+#' @param y A SpatialPoints, SpatialPointsDataFrame, SpatialPolygons, SpatialPolygonsDataFrame, SpatialLines, SpatialLinesDataFrame, or extent. \link{bfastmonitor} will be ran at these locations. In case each feature of the object covers several pixels (typically SpatialPolygons(DataFrames), SpatialLines(DataFrames) and extent), an aggregation function (\code{fun=}) has to be supplied (see \link{extract}).
+#' @param start See \code{\link{bfastmonitor}}
+#' @param formula See \code{\link{bfastmonitor}}
+#' @param order See \code{\link{bfastmonitor}}
+#' @param lag See \code{\link{bfastmonitor}}
+#' @param slag See \code{\link{bfastmonitor}}
+#' @param history See \code{\link{bfastmonitor}}
+#' @param type See \code{\link{bfastmonitor}}
+#' @param h See \code{\link{bfastmonitor}}
+#' @param level See \code{\link{bfastmonitor}}
+#' @param mc.cores Numeric NUmber of cores to use (for parallel processing)
+#' @param ... Arguments to be passed to \link{zooExtract}
+#' 
+#' @author Loic Dutrieux
+#' 
+#' @examples
+#' # Load data
+#' data(tura)
+#' 
+#' # 1- SpatialPoints case
+#' # Generate SpatialPoints
+#' sp <- sampleRegular(x = tura, size = 20, sp=TRUE)
+#' 
+#' # Run bfmSpOver with monitoring period starting year 2005 and all other default parameters of bfastmonitor
+#' out <- bfmSpOver(tura, y = sp, start=c(2005,1))
+#' 
+#' # Visualize the results
+#' plot(tura, 166)
+#' 
+#' # Build color palette
+#' colfunc <- colorRampPalette(c("yellow", "red"))
+#' colList <- colfunc(2013 - 2005)
+#' points(out, col= colList[out$breakpoint - 2005], pch=16, cex = abs(out$magnitude/max(out$magnitude)))
+#' # Color corresponds to timing of break and size to magnitude
+#' 
+#' # 2 - SpatialPolygons case
+#' data(turaSp)
+#' # Run bfmSpOver with monitoring period starting year 2002 and mean spatial aggregation function
+#' out2 <- bfmSpOver(tura, y = turaSp, fun = mean, start=c(2002,1))
+#' 
+#' # Visualize
+#' plot(tura, 166)
+#' # Build color palette
+#' colfunc <- colorRampPalette(c("yellow", "red"))
+#' colList <- colfunc(2013 - 2002)
+#' plot(out2, col = colList[out2$breakpoint - 2002], add = TRUE)
+#' # Interpretation: The redder the latter the break was detected. If transparent, no break detected in spatially aggregated polygon time-series.
+#'
+#' 
+#' @import raster
+#' @import sp
+#' 
+#' @export
+
+bfmSpOver <- function(x, y, start, 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, ...) {
+
+    ts <- zooExtract(x, y, ...)
+    bfm <- bfmZoo(ts, mc.cores = mc.cores, start = start,
+                  formula=formula,
+                  order=order, lag=lag, slag=slag,
+                  history=history,
+                  type=type, h=h,
+                  end=end, level=level)
+
+    if(inherits(y, 'SpatialPoints')) {
+        y <- SpatialPoints(y) # Not sure if that is necessary
+        out <- SpatialPointsDataFrame(y, data = bfm)
+    } else if(inherits(y, 'SpatialPolygons')) {
+        out <- SpatialPolygonsDataFrame(y, data = bfm, match.ID = FALSE)
+    } else if(inherits(y, 'SpatialLines')) {
+        out <- SpatialLinesDataFrame(y, data = bfm, match.ID = FALSE)
+    } else if(class(y) == 'extent') {
+        y <- polygonFromExtent(y)
+        out <- SpatialPolygonsDataFrame(y, data = bfm, match.ID = FALSE)
+        proj4string(out) <- CRS(projection(x))
+    }
+    return(out)
+}

R/bfmSpatial.R

@@ -14,14 +14,14 @@
 #' @param slag See \code{\link{bfastmonitor}}
 #' @param history See \code{\link{bfastmonitor}}
 #' @param type See \code{\link{bfastmonitor}}
-#' @param n See \code{\link{bfastmonitor}}
+#' @param h See \code{\link{bfastmonitor}}
 #' @param level See \code{\link{bfastmonitor}}
 #' @param mc.cores Numeric. Number of cores to be used for the job.
 #' @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 layers depending on what has been supplied to \code{returnLayers}. See details for more information.
+#' @return A rasterBrick with layers depending on what has been supplied to \code{returnLayers}. By default, 3 layers are returned: (1) breakpoint: timing of breakpoints detected for each pixel; (2) magnitude: the median of the residuals within the monitoring period; (3) error: a value of 1 for pixels where an error was encountered by the algorithm (see \code{\link{try}}), and NA where the method was successfully run. See \code{\link{bfastmonitor}} for more information on the other possible layers.
 #' 
 #' @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.

R/bfmZoo.R

@@ -0,0 +1,47 @@
+#' Runs bfastmonitor on a zoo object
+#' 
+#' @description This function is analog to the \link{bfastmonitor} function, but differs in terms of inputs. In \link{bfastmonitor} an object of class ts needs to be provided, usually pre-processed using \link{bfastts}; the present function accepts directly zoo time-series. The return also differs from \link{bfastmonitor}; instead of returning an object of class 'bfastmonitor'the present function returns a dataframe. The zoo object may contain several time-series, which results in a return dataframe containing several rows.
+#' 
+#' @param x A zoo object that may contain several time-series
+#' @param mc.cores Numeric For parallel processing. Number of workers.
+#' @param ... Arguments to be passed to \link{bfastmonitor}
+#' 
+#' @return A dataframe
+#' 
+#' @author Loic Dutrieux
+#' 
+#' @import bfast
+#' 
+#' @export
+#' 
+
+
+bfmZoo <- function(x, mc.cores = 1, ...) {
+
+    bfm2df <- function(x) {
+        if(class(x) == 'bfastmonitor') {
+            data.frame(breakpoint = x$breakpoint,
+                       magnitude = x$magnitude,
+                       history = (x$history[2] - x$history[1]),
+                       rsq = summary(x$model)$r.squared,
+                       adj_rsq = summary(x$model)$adj.r.squared)
+        } else if(class(x) == 'try-error') {
+            data.frame(breakpoint = NA,
+                       magnitude = NA,
+                       history = NA,
+                       rsq = NA,
+                       adj_rsq = NA)
+        }
+
+    }
+
+    bfastmonitorFun <- function(x, ...) {
+        bfm <- try(bfastmonitor(x, ...))
+        bfm2df(bfm)
+    }
+
+    ts <- bfastts(x, index(x), 'irregular')
+    out <- mclapply(X = ts, FUN = bfastmonitorFun, mc.cores = mc.cores, ...)
+    # Convert list to df
+    do.call(rbind, out)
+}

R/cleanMODIS.R

@@ -6,7 +6,7 @@
 #' @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 QC_val Numeric or vector of numerics. Quality control values to \bold{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}}.
 #' 

R/processLandsat.R

@@ -3,7 +3,7 @@
 #' 
 #' @description Processes a single Landsat scene, from tarball or zip archive (or hdf/tiff if untar is set to FALSE) to vegetation index. Easy to batch using sapply or mclapply for parallel implementation. Data obtained from espi may already contained pre-processed indices layers, in which case they are directly used.
 #' @param x Character. filename of the tarball or zip archive of the hdf/tiff file.
-#' @param vi Character. Vegetation index to be computed or extracted from the archive. Can be either 'ndvi', 'evi', 'savi', 'ndmi'*, 'nbr', 'nbr2'* or 'msavi'*. Indices with * need to be present in the archive.
+#' @param vi Character. Vegetation index to be computed or extracted from the archive. Can be either 'ndvi', 'evi', 'savi', 'ndmi'*, 'nbr', 'nbr2'* or 'msavi'*. Indices with * need to be present in the archive. Note that it is also possible to extract single bands using the \code{vi=} argument. \code{vi='sr_band1'} for instance will extract surface reflectance band 1 from the archive and perform the same pre-processing steps as if it was a vegetation index layer.
 #' @param srdir Character. Directory where the tarball should be uncompressed. Can be ommited if \code{untar} is set to \code{FALSE}
 #' @param outdir Character. Directory where the vegetation index rasterLayer should be written.
 #' @param untar Logical. Is there a need to untar data, or have they been previously unpacked.

R/test_data.R

@@ -8,3 +8,18 @@
 NULL
 
 
+#' SpatialPolygons object
+#' @aliases turaSp
+#' @name turaSp
+#' @description Three polygons located in the tura area (ethiopia)
+#' @docType data
+#' @usage data(turaSp)
+#' @author Ben DeVries
+#' @keywords data
+#' @examples
+#' data(turaSp)
+#' data(tura)
+#' plot(tura, 2)
+#' plot(turaSp, add = TRUE)
+
+NULL

R/timeStack.R

@@ -59,8 +59,9 @@ timeStack <- function(x, pattern=NULL, orderChrono = TRUE, ...) {
     }
 
     s <- stack(x)
+    names(s) <- row.names(getSceneinfo(x))
 
-    time <- getSceneinfo(str_extract(string=basename(x), '(LT4|LT5|LE7|LC8)\\d{13}'))$date
+    time <- getSceneinfo(x)$date
     s <- setZ(x=s, z=time)
 
     if(hasArg(filename)) {