Remove dependency from day of the year

DESCRIPTION

@@ -3,8 +3,9 @@ Type: Package
 Title: Time-Weighted Dynamic Time Warping for Satellite Image Time Series
     Analysis
 Version: 0.2.2.9000
-Date: 2016-07-19
-Authors@R: c(person('Victor', 'Maus', role = c('aut', 'cre'), email = '[email protected]'))
+Date: 2016-11-20
+Authors@R: c(person('Victor', 'Maus', role = c('aut', 'cre'), email = '[email protected]'), 
+             person('Marius', 'Appel', role = c('ctb')))
 Description: Provides an implementation of the Time-Weighted Dynamic Time
     Warping (TWDTW) method for land use and land cover mapping using satellite
     image time series. TWDTW is based on the Dynamic Time Warping technique and

R/class-twdtwRaster.R

@@ -119,20 +119,22 @@ setMethod("initialize",
   definition = 
     function(.Object, timeseries, timeline, doy, layers, labels, levels){
 
-      .Object@timeseries = list(doy=brick(), Layer0=brick())
+      .Object@timeseries = list(Layer0=brick())
       .Object@timeline = as.Date(0)
       .Object@labels = as.character()
       .Object@levels = numeric()
       if(!missing(timeseries)){
         if(is(timeseries, "RasterBrick") | is(timeseries, "RasterStack") | is(timeseries, "RasterLayer") ) 
-           timeseries = list(timeseries)
-        if(is.null(names(timeseries))) names(timeseries) = paste0("Layer",seq_along(timeseries))
-        .Object@timeseries = c(.Object@timeseries[1], timeseries)
+          timeseries = list(timeseries)
+        .Object@timeseries = timeseries
+        if(is.null(names(.Object@timeseries)))
+          names(.Object@timeseries) = paste0("Layer", seq_along(.Object@timeseries)-1)
+      } else {
+        if(!missing(layers))
+          names(.Object@timeseries) = layers
       }
-      if(missing(doy)) doy = creat.doy(.Object@timeseries[[2]], timeline)
-      .Object@timeseries$doy = doy 
-      if(!missing(layers))
-        names(.Object@timeseries) = c("doy", layers)
+      if(!missing(doy)) 
+        .Object@timeseries = c(doy = doy, .Object@timeseries)
       .Object@layers = names(.Object@timeseries)
       if(!missing(labels))
         .Object@labels = as.character(labels)
@@ -142,7 +144,7 @@ setMethod("initialize",
       if(!missing(timeline))
         .Object@timeline = as.Date(timeline)
       validObject(.Object)
-      names(.Object@timeline) = paste0("date.",format(.Object@timeline,"%Y.%m.%d"))
+      names(.Object@timeline) = paste0("date.", format(.Object@timeline,"%Y.%m.%d"))
       .Object@timeseries = lapply(.Object@timeseries, function(x) { names(x)=names(.Object@timeline); x})
       return(.Object)
   }
@@ -195,39 +197,48 @@ setMethod(f = "twdtwRaster",
                                 layers=layers, labels=labels, levels=levels, filepath=filepath, dotargs=dotargs)
           })
 
+
 creat.twdtwRaster = function(timeseries, timeline, doy, layers, labels, levels, filepath, dotargs){
 
-  if(is.null(doy)) doy = creat.doy(timeseries[[1]], timeline)
-
   # Check timeline 
-  nl = sapply(c(timeseries, doy), nlayers)
+  nl = sapply(c(timeseries), nlayers)
+  if(!is.null(doy))
+    nl = c(nlayers(doy), nl)
   if(any(nl!=length(timeline)))
     stop("raster objects do not have the same length as the timeline")
 
   res = timeseries
   # Save a single file (complete time series) for each raster attribute 
   if (!is.null(filepath)) {
-#     print("Saving raster objects. It might take some minutes depending on the objects size and format.")
     dir.create(filepath, showWarnings = FALSE)
     write(as.character(timeline), file = paste(filepath, "timeline", sep="/"))
-    aux = c(doy=doy, res)
+    aux = res
+    if(!is.null(doy))
+      aux = c(doy=doy, res)
     res_brick = lapply(names(aux), function(i){
       filename = paste(filepath, i, sep="/")
       dotargs = c(x = aux[[i]], filename = filename, dotargs)
       r = do.call(writeRaster, dotargs)
       r
     })
     names(res_brick) = names(aux)
-    res = res_brick[-1]
-    doy = res_brick[[1]]
+    doy = NULL
+    res = res_brick
+    if(any(names(res)=="doy")){
+      res = res_brick[-1]
+      doy = res_brick[[1]]
+    }
   }
   if(is.null(layers)) layers = names(res)
+  if(is.null(doy)) 
+    return(new("twdtwRaster", timeseries = res, timeline = timeline, layers = layers, labels = labels, levels=levels))
   new("twdtwRaster", timeseries = res, timeline = timeline, doy=doy, layers = layers, labels = labels, levels=levels)
 }
 
-creat.doy = function(x, timeline){
+.creat.doy = function(x, timeline){
     array_data = rep(as.numeric(format(as.Date(timeline), "%j")), each=ncell(x))
     e = extent(x)
     brick(array(array_data, dim = dim(x)), xmn=e[1], xmx=e[2], ymn=e[3], ymx=e[4], crs = projection(x))
 }
 
+

R/class-twdtwTimeSeries.R

@@ -40,7 +40,8 @@
 #'
 #' @examples 
 #' # Creating new object of class twdtwTimeSeries  
-#' ptt = new("twdtwTimeSeries", timeseries = MOD13Q1.patterns.list, labels = names(MOD13Q1.patterns.list))
+#' ptt = new("twdtwTimeSeries", timeseries = MOD13Q1.patterns.list, 
+#'            labels = names(MOD13Q1.patterns.list))
 #' class(ptt)
 #' labels(ptt)
 #' levels(ptt)

R/dtw.R

@@ -88,7 +88,7 @@
                 XM = matrix(as.integer(c(as.numeric(x[[1]]$from), as.numeric(x[[1]]$to))), ncol = 2),
                 AM = matrix(as.double(fill), nrow = n, ncol = m), 
                 DM = as.double(x[[1]]$distance),
-                DP  = as.integer(as.numeric(breaks)),
+                DP = as.integer(as.numeric(breaks)),
                 X  = as.integer(match(x[[1]]$label, levels)),
                 IM = matrix(as.integer(0), nrow = n, ncol = 3),
                 A  = as.integer(x[[1]]$Alig.N),

R/getTimeSeries.R

@@ -167,9 +167,12 @@ extractTimeSeries.twdtwRaster = function(x, y){
   # Check if the sample time interval overlaps the raster time series 
   from  = try(as.Date(s$from))
   to    = try(as.Date(s$to))
-  doy   = c(x$doy[p,])
-  year  = format(timeline, "%Y")
-  dates = getDatesFromDOY(year = year, doy = doy)
+  dates = timeline
+  if(any(names(x)=="doy")){
+    doy   = c(x$doy[p,])
+    year  = format(timeline, "%Y")
+    dates = getDatesFromDOY(year = year, doy = doy)
+  }
   if(is.null(from) | is(from, "try-error")) from = dates[1]
   if(is.null(to) | is(to, "try-error")) to = tail(dates,1)
   # layer =  which( from - dates <= 0 )[1]
@@ -181,8 +184,7 @@ extractTimeSeries.twdtwRaster = function(x, y){
     return(NULL)
   }
   # Extract raster values 
-  I = which(names(x) %in% c("doy"))
-  ts = data.frame(sapply(x[-I], function(x) x[p,layer:(layer+nl-1)]))
+  ts = data.frame(sapply(x[!names(x)%in%c("doy")], function(x) x[p,layer:(layer+nl-1)]))
   dates = dates[layer:(layer+nl-1)]
   k = !duplicated(dates)
   zoo(data.frame(ts[k,, drop=FALSE]), dates[k])

R/methods.R

@@ -309,7 +309,9 @@ setMethod("[", "twdtwRaster", function(x, i) {
   if(any(i > nlayers(x)))
     stop("subscript out of bounds")
   if(any(is.na(i))) stop("NA index not permitted")
-  new("twdtwRaster", timeseries=x@timeseries[i], timeline = x@timeline, doy = x@timeseries[[1]])
+  if(any("doy"==layers(x)))
+    return(new("twdtwRaster", timeseries=x@timeseries[i], timeline = x@timeline, doy = x@timeseries[[1]]))
+  new("twdtwRaster", timeseries=x@timeseries[i], timeline = x@timeline)
 })
 
 #' @inheritParams twdtwRaster-class
@@ -401,7 +403,7 @@ setMethod("crop",
 setMethod("extent", 
           signature = signature("twdtwRaster"),
           definition = function(x, y, ...){
-            extent(x@timeseries$doy)
+            extent(x@timeseries[[1]])
           }
 )
 

R/plot.R

@@ -121,7 +121,8 @@ setMethod("plot", signature(x = "twdtwRaster"), function(x, type="maps", ...) .P
           time.labels=NULL
         }
         else {
-          layers = coverages(x)[-1]
+          layers = coverages(x)
+          layers = layers[!layers%in%"doy"]
           time.levels = time.levels[1]
           time.labels = time.labels[1]
           labels = layers
@@ -162,7 +163,7 @@ setMethod("plot", signature(x = "twdtwRaster"), function(x, type="maps", ...) .P
         class.levels = levels(x)
 
       if(length(class.levels)<1)
-        class.levels = sort(unique(as.numeric(x[[2]][])))
+        class.levels = sort(unique(as.numeric(x[["Class"]][])))
 
       if( is.null(class.labels))  
         class.labels = labels(x)
@@ -185,7 +186,7 @@ setMethod("plot", signature(x = "twdtwRaster"), function(x, type="maps", ...) .P
       names(class.levels) = class.labels
       names(class.labels) = class.labels
 
-      x = subset(x=x[[2]], subset=time.levels)
+      x = subset(x=x[["Class"]], subset=time.levels)
 
       pt = pmatch(type, c("maps","area","changes"))
 

R/plotDistance.R

@@ -58,9 +58,7 @@
 #' r_twdtw = twdtwApply(x=rts, y=patt, weight.fun=log_fun, breaks=time_interval, 
 #'           filepath="~/test_twdtw", overwrite=TRUE, format="GTiff", mc.cores=3)
 #' 
-#' r_lucc = twdtwClassify(r_twdtw, format="GTiff")
-#' 
-#' plotMaps(r_lucc)
+#' plotDistance(r_twdtw)
 #' 
 #' }
 #' @export

R/plotMaps.R

@@ -64,7 +64,7 @@
 #' r_twdtw = twdtwApply(x=rts, y=patt, weight.fun=log_fun, breaks=time_interval, 
 #'           filepath="~/test_twdtw", overwrite=TRUE, format="GTiff", mc.cores=3)
 #' 
-#' r_lucc = twdtwClassify(r_twdtw, format="GTiff")
+#' r_lucc = twdtwClassify(r_twdtw, format="GTiff", overwrite=TRUE)
 #' 
 #' plotMaps(r_lucc)
 #' 

R/subset.R

@@ -134,9 +134,10 @@ setMethod("subset", "twdtwRaster", function(x, e=NULL, layers=NULL)
           subset.twdtwRaster(x=x, e=e, layers=layers) )
 
 subset.twdtwRaster = function(x, e, layers){
-    if(is.null(layers)) layers = names(x)
-    layers = c("doy", layers[layers!="doy"])
-    if(is.null(e)) e = extent(x@timeseries$doy)
+    if(is.null(layers)) 
+      layers = names(x)
+    if(is.null(e))
+      e = extent(x)
     res = x
     res@layers = layers
     res@timeseries = res@timeseries[layers]

R/twdtwApply.R

@@ -227,7 +227,7 @@ twdtwApply.twdtwRaster = function(x, y, weight.fun, dist.method, step.matrix, n,
 
     # Set blocks Multi-thread parameters
     minblocks = round(nrow(x)*ncol(x) / chunk.size)
-    blocks = blockSize(x@timeseries$doy, minblocks = minblocks)
+    blocks = blockSize(x@timeseries[[1]], minblocks = minblocks)
     threads = seq(1, blocks$n)
 
     # Match raster bands to pattern bands
@@ -245,7 +245,7 @@ twdtwApply.twdtwRaster = function(x, y, weight.fun, dist.method, step.matrix, n,
 
     # Open raster fiels for results 
     if(is.null(filepath)) filepath = paste0(getwd(), "/twdtw_results")
-    r_template = brick(x@timeseries$doy, nl=length(breaks)-1)
+    r_template = brick(x@timeseries[[1]], nl=length(breaks)-1)
     dir.create(filepath, showWarnings = FALSE, recursive = TRUE)
     filename = paste0(filepath,"/twdtw_distance_from_",levels)
     names(filename) = levels
@@ -269,7 +269,7 @@ twdtwApply.twdtwRaster = function(x, y, weight.fun, dist.method, step.matrix, n,
       ts_list = mclapply(as.list(x), FUN=getValuesBlock, row=blocks$row[i], nrows=blocks$nrows[i], mc.preschedule = mc.preschedule, mc.set.seed = mc.set.seed, mc.silent = mc.silent, mc.cores = mc.cores, mc.cleanup = mc.cleanup)
 
       # Create a dummy array 
-      nts = seq(1, nrow(ts_list$doy))
+      nts = seq(1, nrow(ts_list[[1]]))
       m = length(levels)
       n = length(breaks)-1
 
@@ -332,7 +332,12 @@ twdtwApply.twdtwRaster = function(x, y, weight.fun, dist.method, step.matrix, n,
 .bulidZoo = function(p, x, timeline){
   # Get time series for each band 
   datasets = lapply(x, function(x) x[p,])
-  datasets$doy = getDatesFromDOY(doy=datasets$doy, year=format(timeline, "%Y"))
+  if(any("doy"==names(datasets))){
+    datasets$doy = getDatesFromDOY(doy=datasets$doy, year=format(timeline, "%Y"))
+  } else {
+    datasets$doy = timeline
+  }
+
   idoy = which(names(datasets) %in% c("doy"))
 
   # Remove invalid values 

R/twdtwClassify.R

@@ -144,7 +144,8 @@ setMethod("twdtwClassify", "twdtwMatches",
 #' }
 setMethod("twdtwClassify", "twdtwRaster",
           function(x, patterns.labels=NULL, thresholds=Inf, fill=255, filepath, ...){
-                  if(is.null(patterns.labels)) patterns.labels = coverages(x)[-1]
+                  if(is.null(patterns.labels)) patterns.labels = coverages(x)
+                  patterns.labels = patterns.labels[!patterns.labels%in%"doy"]
                   if(missing(filepath)) filepath = if(fromDisk(x[[2]])){dirname(filename(x[[2]]))}else{NULL}
                   twdtwClassify.twdtwRaster(x, patterns.labels=patterns.labels, thresholds=thresholds, fill=fill, filepath=filepath, ...)
            })

src/computecost.f

@@ -10,6 +10,8 @@
 C                                                             C
 C   Efficient computation of DTW cost matrix  - 2015-10-16    C
 C                                                             C
+C   This function was adpted from the C function 'computeCM'  C
+C   implemented in the R package 'dtw' by Toni Giorgino.      C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C     CM - Input local cost and output cumulative cost matrix