From 19745a93fdfe41851b861ef63b1f9080e4412c76 Mon Sep 17 00:00:00 2001
From: vwmaus <vwmaus1@gmail.com>
Date: Tue, 30 May 2017 17:03:18 +0200
Subject: [PATCH] First version of parallel processing. Unstable.

---
 DESCRIPTION               |   1 +
 NAMESPACE                 |   2 +
 R/twdtwApplyParallel.R    | 345 ++++++++++++++++++++++++++++++++++++++
 dtwSat.Rproj              |   6 +-
 man/twdtwApplyParallel.Rd | 151 +++++++++++++++++
 5 files changed, 502 insertions(+), 3 deletions(-)
 create mode 100644 R/twdtwApplyParallel.R
 create mode 100644 man/twdtwApplyParallel.Rd

diff --git a/DESCRIPTION b/DESCRIPTION
index 361ac1f..dc8c867 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -92,6 +92,7 @@ Collate:
     'subset.R'
     'twdtw.R'
     'twdtwApply.R'
+    'twdtwApplyParallel.R'
     'twdtwAssess.R'
     'twdtwClassify.R'
     'twdtwCrossValidate.R'
diff --git a/NAMESPACE b/NAMESPACE
index 09bcdd4..2fcbd71 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -26,6 +26,7 @@ export(shiftDates)
 export(symmetric1)
 export(symmetric2)
 export(twdtwApply)
+export(twdtwApplyParallel)
 export(twdtwClassify)
 exportMethods("[")
 exportMethods("[[")
@@ -64,6 +65,7 @@ exportMethods(show)
 exportMethods(subset)
 exportMethods(summary)
 exportMethods(twdtwApply)
+exportMethods(twdtwApplyParallel)
 exportMethods(twdtwAssess)
 exportMethods(twdtwClassify)
 exportMethods(twdtwCrossValidate)
diff --git a/R/twdtwApplyParallel.R b/R/twdtwApplyParallel.R
new file mode 100644
index 0000000..4c94749
--- /dev/null
+++ b/R/twdtwApplyParallel.R
@@ -0,0 +1,345 @@
+#' @include methods.R
+#' @title Apply TWDTW analysis to twdtwRaster using parallel processing 
+#' @name twdtwApplyParallel 
+#' @author Victor Maus, \email{vwmaus1@@gmail.com}
+#' 
+#' @description This function performs a multidimensional Time-Weighted DTW 
+#' analysis and retrieves the matches between the temporal patterns and 
+#' a set of time series [1].
+#' 
+#' @inheritParams twdtwApply
+#' 
+#' @references 
+#' [1] Maus  V,  Camara  G,  Cartaxo  R,  Sanchez  A,  Ramos  FM,  de Queiroz, GR.
+#' (2016). A Time-Weighted Dynamic Time Warping method for land use and land cover 
+#' mapping. IEEE Journal of Selected Topics in Applied Earth Observations and Remote 
+#' Sensing, vol.9, no.8, pp.3729-3739.
+#' @references 
+#' [2] Giorgino, T. (2009). Computing and Visualizing Dynamic Time Warping Alignments in R: 
+#' The dtw Package. Journal of Statistical Software, 31, 1-24.
+#' @references 
+#' [3] Muller, M. (2007). Dynamic Time Warping. In Information Retrieval for Music 
+#' and Motion (pp. 79-84). London: Springer London, Limited.
+#' 
+#' @details The linear \code{linearWeight} and \code{logisticWeight} weight functions 
+#' can be passed to \code{twdtwApply} through the argument \code{weight.fun}. This will 
+#' add a time-weight to the dynamic time warping analysis. The time weight 
+#' creates a global constraint useful to analyse time series with phenological cycles
+#' of vegetation that are usually bound to seasons. In previous studies by [1] the 
+#' logistic weight had better results than the linear for land cover classification. 
+#' See [1] for details about the method. 
+#' 
+#' @return An object of class twdtwRaster.
+#' 
+#' @seealso 
+#' \code{\link[dtwSat]{twdtwRaster-class}}, and 
+#' \code{\link[dtwSat]{createPatterns}}
+#' 
+#' @export  
+setGeneric(name = "twdtwApplyParallel", 
+           def = function(x, y, resample=TRUE, length=NULL, weight.fun=NULL, 
+                          dist.method="Euclidean", step.matrix = symmetric1, n=NULL, 
+                          span=NULL, min.length=0, theta = 0.5, ...) standardGeneric("twdtwApplyParallel"))
+
+
+
+#' @rdname twdtwApplyParallel 
+#' @aliases twdtwApplyParallel-twdtwRaster
+#' @examples
+#' \dontrun{
+#' # Run TWDTW analysis for raster time series 
+#' patt = MOD13Q1.MT.yearly.patterns
+#' evi = brick(system.file("lucc_MT/data/evi.tif", package="dtwSat"))
+#' ndvi = brick(system.file("lucc_MT/data/ndvi.tif", package="dtwSat"))
+#' red = brick(system.file("lucc_MT/data/red.tif", package="dtwSat"))
+#' blue = brick(system.file("lucc_MT/data/blue.tif", package="dtwSat"))
+#' nir = brick(system.file("lucc_MT/data/nir.tif", package="dtwSat"))
+#' mir = brick(system.file("lucc_MT/data/mir.tif", package="dtwSat"))
+#' doy = brick(system.file("lucc_MT/data/doy.tif", package="dtwSat"))
+#' timeline = scan(system.file("lucc_MT/data/timeline", package="dtwSat"), what="date")
+#' rts = twdtwRaster(evi, ndvi, red, blue, nir, mir, timeline = timeline, doy = doy)
+#' 
+#' time_interval = seq(from=as.Date("2007-09-01"), to=as.Date("2013-09-01"), 
+#'                     by="12 month")
+#' log_fun = weight.fun=logisticWeight(-0.1,50)
+#' 
+#' r_twdtw = twdtwApply(x=rts, y=patt, weight.fun=log_fun, breaks=time_interval, 
+#'           filepath="~/test_twdtw", overwrite=TRUE, format="GTiff")
+#'
+#' plot(r_twdtw, type="distance")
+#' 
+#' r_lucc = twdtwClassify(r_twdtw, format="GTiff", overwrite=TRUE)
+#' 
+#' plot(r_lucc)
+#' 
+#' plot(r_lucc, type="distance")
+#' 
+#' }
+#' @export
+setMethod(f = "twdtwApplyParallel", "twdtwRaster",
+          def = function(x, y, resample, length, weight.fun, dist.method, step.matrix, n, span, min.length, theta, 
+                         breaks=NULL, from=NULL, to=NULL, by=NULL, overlap=0.5, chunk.size=1000, filepath=NULL, silent=FALSE, ...){
+            if(!is(step.matrix, "stepPattern"))
+              stop("step.matrix is not of class stepPattern")
+            if(is.null(weight.fun))
+              weight.fun = function(psi) 0 
+            if(!is(weight.fun, "function"))
+              stop("weight.fun is not a function")
+            if( overlap < 0 & 1 < overlap )
+              stop("overlap out of range, it must be a number between 0 and 1")
+            if(is.null(breaks))
+              if( !is.null(from) &  !is.null(to) ){
+                breaks = seq(as.Date(from), as.Date(to), by=by)    
+              } else {
+                patt_range = lapply(index(y), range)
+                patt_diff = trunc(sapply(patt_range, diff)/30)+1
+                min_range = which.min(patt_diff)
+                by = patt_diff[[min_range]]
+                from = patt_range[[min_range]][1]
+                to = from 
+                month(to) = month(to) + by
+                year(from) = year(range(index(x))[1])
+                year(to) = year(range(index(x))[2])
+                if(to<from) year(to) = year(to) + 1
+                breaks = seq(from, to, paste(by,"month"))
+              }
+            breaks = as.Date(breaks)
+            if(resample)
+              y = resampleTimeSeries(object=y, length=length)
+            twdtwApplyParallel.twdtwRaster(x, y, weight.fun, dist.method, step.matrix, n, span, min.length, theta, 
+                                   breaks, overlap, chunk.size, filepath, silent, ...)
+          })
+
+
+twdtwApplyParallel.twdtwRaster = function(x, y, weight.fun, dist.method, step.matrix, n, span, min.length, theta, 
+                                  breaks, overlap, chunk.size, filepath, silent, ...){
+  
+  # Set blocks Multi-thread parameters
+  minblocks = round(nrow(x)*ncol(x) / chunk.size)
+  blocks = blockSize(x@timeseries[[1]], minblocks = minblocks)
+  threads = seq(1, blocks$n)
+  
+  # Match raster bands to pattern bands
+  raster_bands = coverages(x)
+  pattern_names = names(y@timeseries[[1]])
+  matching_bands = pattern_names[pattern_names %in% raster_bands]
+  if(length(matching_bands)<1)
+    stop(paste0("Attributes (bands) of the raster and patterns do not match"))
+  if("doy"%in%coverages(x) & !"doy"%in%matching_bands)
+    matching_bands = c("doy", matching_bands)
+  x = subset(x, layers=matching_bands)
+  raster_bands = coverages(x)
+  
+  # Set raster levels and labels 
+  levels = levels(y)
+  names(levels) = levels
+  
+  # Open raster fiels for results 
+  if(is.null(filepath)) filepath = paste0(getwd(), "/twdtw_results")
+  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
+  # b_files = lapply(filename, function(i) writeStart(r_template, filename=i, format="GTiff", overwrite=TRUE))
+  b_files = lapply(filename, function(i) writeStart(r_template, filename=i, ...))
+  
+  # Get time line 
+  timeline = as.Date(index(x))
+  
+  # Raster info 
+  proj_str = projection(x)
+  coord    = data.frame(coordinates(x))
+  names(coord) = c("longitude","latitude")
+  
+  fun = function(i){
+    
+    if(!silent) print(paste0("Procesing chunk ",i,"/",threads[length(threads)]))
+    
+    # Get twdtwTimeSeries from raster 
+    cells = cellFromRow(x@timeseries[[1]], blocks$row[i]:blocks$nrows[i])
+    ts = getTimeSeries(x, y = coord[cells,], proj4string = proj_str)
+    
+    # Apply TWDTW analysis  
+    twdtw_results = twdtwApply(ts, y=y, weight.fun=weight.fun, dist.method=dist.method, 
+                               step.matrix=step.matrix, n=n, span=span, 
+                               min.length=min.length, theta=theta, keep=FALSE)
+    
+    # Get best mathces for each point, period, and pattern 
+    m = length(levels)
+    h = length(breaks)-1
+    A = lapply(as.list(twdtw_results), FUN=.lowestDistances, m=m, n=h, levels=levels, breaks=breaks, overlap=overlap, fill=9999)
+    
+    # Reshape list to array 
+    A = sapply(A, matrix, nrow=h, ncol=m, simplify = 'array')
+    
+    # Write raster files 
+    lapply(seq_along(levels), function(l) writeValues(b_files[[levels[l]]], matrix(t(A[,l,]),ncol=h), blocks$row[i]))
+  }
+  
+  # Apply TWDTW analysis 
+  out.list = lapply(threads, FUN = fun)
+  
+  # Close raster files 
+  b_files = lapply(b_files, writeStop)
+  
+  # Create brick list for output 
+  out = lapply(sapply(b_files, filename), brick)
+  
+  # Save output timeline 
+  timeline = breaks[-1]
+  write(as.character(timeline), file = paste(filepath, "timeline", sep="/"))
+  
+  new("twdtwRaster", timeseries = out, timeline=timeline, layers = names(out))
+}
+
+
+
+
+# Prallel raster processing 
+apply_raster_parallel <- function(x, fun, A, filepath="", ...) {
+  
+  # Set blocks Multi-thread parameters
+  # minblocks <- round(nrow(x) * ncol(x) / chunk.size)
+  # blocks <- blockSize(x@timeseries[[1]], minblocks = minblocks)
+  # threads <-  seq(1, blocks$n)
+  
+  # Match raster bands to pattern bands
+  raster_bands <- coverages(x)
+  pattern_names <- names(y@timeseries[[1]])
+  matching_bands <- pattern_names[pattern_names %in% raster_bands]
+  
+  if(length(matching_bands) < 1)
+    stop(paste0("Bands from twdtwRaster do not match the bands from patterns"))
+  
+  if("doy" %in% coverages(x) & !"doy" %in% matching_bands)
+    matching_bands <- c("doy", matching_bands)
+  
+  x <- subset(x, layers = matching_bands)
+  raster_bands <- coverages(x)
+  
+  # Set raster levels and labels 
+  levels <- levels(y)
+  names(levels) <- levels
+  
+  # Create output raster objects 
+  # if(is.null(filepath)) 
+  #   filepath = paste0(getwd(), "/twdtw_results")
+  r_template <- brick(x@timeseries[[1]], nl = length(breaks) - 1, values = FALSE)
+  out <- rep(list(r_template), length(levels))
+  names(out) <- names(levels)
+  # dir.create(filepath, showWarnings = FALSE, recursive = TRUE)
+  # filename = paste0(filepath,"/twdtw_distance_from_",levels)
+  # names(filename) = levels
+  # b_files = lapply(filename, function(i) writeStart(r_template, filename=i, format="GTiff", overwrite=TRUE))
+  # b_files = lapply(filename, function(i) writeStart(r_template, filename=i, ...))
+  
+  # out <- brick(x@timeseries[[1]], nl=length(breaks)-1, values=FALSE)
+  # names(out) <- names(x)
+  
+  cl <- getCluster()
+  on.exit( returnCluster() )
+  
+  nodes <- length(cl)
+  
+  bs <- blockSize(x, minblocks = nodes*4)
+  bs$array_rows <- cumsum(c(1, bs$nrows*out@ncols))
+  pb <- pbCreate(bs$n)
+  
+  clFun <- function(k){
+    
+    # Get raster data 
+    # v <- lapply(x@timeseries, getValues, row = 1, nrows = 2)
+    v <- lapply(x, getValues, row = bs$row[k], nrows = bs$nrows[k])
+    
+    # Create twdtwTimeSeries 
+    ts <- twdtwTimeSeries(lapply(1:nrow(v[[1]]), function(i){
+      # Get time series dates 
+      if(any(names(v) %in% "doy")){
+        timeline <- getDatesFromDOY(year = format(index(x), "%Y"), doy = v[["doy"]][i,])
+      } else {
+        timeline <- index(x)
+      }
+      # Tag duplicate dates for removal 
+      k = !duplicated(timeline)
+      # Create multi band time series 
+      zoo(sapply(v[!(names(v) %in% "doy")], function(v) v[i, k, drop = FALSE]), timeline[k])
+    }))
+    
+    # Apply TWDTW analysis  
+    twdtw_results = twdtwApply(ts, y=y, weight.fun=weight.fun, dist.method=dist.method, 
+                               step.matrix=step.matrix, n=n, span=span, 
+                               min.length=min.length, theta=theta, keep=FALSE)
+    
+    # Get best mathces for each point, period, and pattern 
+    m = length(levels)
+    h = length(breaks)-1
+    A = lapply(as.list(twdtw_results), FUN=.lowestDistances, m=m, n=h, levels=levels, breaks=breaks, overlap=overlap, fill=9999)
+    
+    t(data.frame(A))
+    
+  }
+  
+  # get all nodes going
+  for (k in 1:nodes) {
+    sendCall(cl[[k]], clFun, list(k, A), tag = k)
+  }
+  
+  filepath <- trim(filepath)
+  filename <- NULL
+  if (filepath != "") {
+    filename <- paste0(filepath, "/", names(out), ".grd")
+  } else if (!canProcessInMemory(r_template, n = length(breaks))) {
+    filename <- sapply(names(out), rasterTmpFile)
+  } 
+  
+  if (!is.null(filename)) {
+    # out <- writeStart(out, filename = filename, ... )
+    out <- lapply(names(out), function(i) writeStart(out[[i]], filename = filename[i], ...))
+  } else {
+    # vv <- matrix(ncol=nrow(out), nrow=ncol(out))
+    # vv <- matrix(out, ncol=nlayers(out))
+    vv <- lapply(names(out), function(i) matrix(out[[i]], ncol = nlayers(out[[i]])))
+  }
+  
+  for (k in 1:bs$n) {
+    # receive results from a node
+    d <- recvOneData(cl)
+    
+    # error?
+    if (! d$value$success) { 
+      stop('cluster error')
+    }
+    
+    # which block is this?
+    b <- d$value$tag
+    cat('received block: ',b," / ",bs$n,'\n'); flush.console();
+    
+    if (!is.null(filename)) {
+      # out <- writeValues(out, d$value$value, bs$row[b])
+      out <- lapply(seq_along(levels), function(l) writeValues(out[[l]], d$value$value[[l]], bs$row[b]))
+    } else {
+      # cols <- bs$row[b]:(bs$row[b] + bs$nrows[b]-1)	
+      # vv[,cols] <- matrix(d$value$value, nrow=out@ncols)
+      rows <- seq(from = bs$array_rows[b], by = 1, length.out = bs$nrows[b]*out@ncols)
+      vv <- lapply(seq_along(levels), function(l) vv[[l]][rows,] <- d$value$value[[l]])
+    }
+    
+    # need to send more data?
+    ni <- nodes + k
+    if (ni <= bs$n) {
+      sendCall(cl[[d$node]], clFun, list(ni, A), tag = ni)
+    }	
+    pbStep(pb)
+  }
+  if (!is.null(filename)) {
+    # out <- writeStop(out)
+    out <- lapply(out, writeStop)
+  } else {
+    # out <- setValues(out, as.vector(vv))
+    out <- lapply(names(out), function(i) setValues(out[[i]], values = as.vector(vv[[i]])))
+  }
+  pbClose(pb)
+  
+  return(out)
+}
+
diff --git a/dtwSat.Rproj b/dtwSat.Rproj
index c850748..14cd891 100644
--- a/dtwSat.Rproj
+++ b/dtwSat.Rproj
@@ -1,8 +1,8 @@
 Version: 1.0
 
-RestoreWorkspace: Default
-SaveWorkspace: Default
-AlwaysSaveHistory: Default
+RestoreWorkspace: No
+SaveWorkspace: No
+AlwaysSaveHistory: No
 
 EnableCodeIndexing: Yes
 UseSpacesForTab: Yes
diff --git a/man/twdtwApplyParallel.Rd b/man/twdtwApplyParallel.Rd
new file mode 100644
index 0000000..740d720
--- /dev/null
+++ b/man/twdtwApplyParallel.Rd
@@ -0,0 +1,151 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/twdtwApplyParallel.R
+\docType{methods}
+\name{twdtwApplyParallel}
+\alias{twdtwApplyParallel}
+\alias{twdtwApplyParallel,twdtwRaster-method}
+\alias{twdtwApplyParallel-twdtwRaster}
+\title{Apply TWDTW analysis to twdtwRaster using parallel processing}
+\usage{
+twdtwApplyParallel(x, y, resample = TRUE, length = NULL,
+  weight.fun = NULL, dist.method = "Euclidean", step.matrix = symmetric1,
+  n = NULL, span = NULL, min.length = 0, theta = 0.5, ...)
+
+\S4method{twdtwApplyParallel}{twdtwRaster}(x, y, resample, length, weight.fun,
+  dist.method, step.matrix, n, span, min.length, theta, breaks = NULL,
+  from = NULL, to = NULL, by = NULL, overlap = 0.5, chunk.size = 1000,
+  filepath = NULL, silent = FALSE, ...)
+}
+\arguments{
+\item{x}{an object of class twdtw*. This is the target time series. 
+Usually, it is a set of unclassified time series.}
+
+\item{y}{an object of class \link[dtwSat]{twdtwTimeSeries}. 
+The temporal patterns.}
+
+\item{resample}{resample the patterns to have the same length. Default is TRUE.
+See \link[dtwSat]{resampleTimeSeries} for details.}
+
+\item{length}{An integer. Patterns length used with \code{patterns.length}. 
+If not declared the length of the output patterns will be the length of 
+the longest pattern.}
+
+\item{weight.fun}{A function. Any function that receive and performs a 
+computation on a matrix. The function receives a matrix of time differences 
+in days and returns a matrix of time-weights. If not declared the time-weight 
+is zero. In this case the function runs the standard version of the dynamic 
+time warping. See details.}
+
+\item{dist.method}{A character. Method to derive the local cost matrix.
+Default is ''Euclidean'' see \code{\link[proxy]{dist}} in package 
+\pkg{proxy}.}
+
+\item{step.matrix}{see \code{\link[dtw]{stepPattern}} in package \pkg{dtw} [2].}
+
+\item{n}{An integer. The maximun number of matches to perform. 
+NULL will return all matches.}
+
+\item{span}{A number. Span between two matches, \emph{i.e.} the minimum
+interval between two matches, for details see [3]. If not declared it removes
+all overlapping matches of the same pattern. To include overlapping matches 
+of the same pattern use \code{span=0}.}
+
+\item{min.length}{A number between 0 an 1. This argument removes the over fittings.
+Minimum length after warping. Percentage of the original pattern length. Default is 0.5, 
+meaning that the matching cannot be shorter than half of the pattern length.}
+
+\item{theta}{numeric between 0 and 1. The weight of the time 
+for the TWDTW computation. Use \code{theta=0} to cancel the time-weight, 
+\emph{i.e.} to run the original DTW algorithm. Default is 0.5, meaning that 
+the time has the same weight as the curve shape in the TWDTW analysis.}
+
+\item{...}{arguments to pass to \code{\link[raster]{writeRaster}}}
+
+\item{breaks}{A vector of class \code{\link[base]{Dates}}. This replaces the arguments \code{from},
+\code{to}, and \code{by}.}
+
+\item{from}{A character or \code{\link[base]{Dates}} object in the format "yyyy-mm-dd".}
+
+\item{to}{A \code{\link[base]{character}} or \code{\link[base]{Dates}} object in the format "yyyy-mm-dd".}
+
+\item{by}{A \code{\link[base]{character}} with the intevals size, \emph{e.g.} "6 month".}
+
+\item{overlap}{A number between 0 and 1. The minimum overlapping 
+between one match and the interval of classification. Default is 0.5, 
+\emph{i.e.} an overlap minimum of 50\%.}
+
+\item{chunk.size}{An integer. Set the number of cells for each block, 
+see \code{\link[raster]{blockSize}} for details.}
+
+\item{filepath}{A character. The path to save the raster with results. If not informed the 
+function saves in the current work directory.}
+
+\item{silent}{if set to TRUE then hide chunk processing message. Default is FALSE.}
+}
+\value{
+An object of class twdtwRaster.
+}
+\description{
+This function performs a multidimensional Time-Weighted DTW 
+analysis and retrieves the matches between the temporal patterns and 
+a set of time series [1].
+}
+\details{
+The linear \code{linearWeight} and \code{logisticWeight} weight functions 
+can be passed to \code{twdtwApply} through the argument \code{weight.fun}. This will 
+add a time-weight to the dynamic time warping analysis. The time weight 
+creates a global constraint useful to analyse time series with phenological cycles
+of vegetation that are usually bound to seasons. In previous studies by [1] the 
+logistic weight had better results than the linear for land cover classification. 
+See [1] for details about the method.
+}
+\examples{
+\dontrun{
+# Run TWDTW analysis for raster time series 
+patt = MOD13Q1.MT.yearly.patterns
+evi = brick(system.file("lucc_MT/data/evi.tif", package="dtwSat"))
+ndvi = brick(system.file("lucc_MT/data/ndvi.tif", package="dtwSat"))
+red = brick(system.file("lucc_MT/data/red.tif", package="dtwSat"))
+blue = brick(system.file("lucc_MT/data/blue.tif", package="dtwSat"))
+nir = brick(system.file("lucc_MT/data/nir.tif", package="dtwSat"))
+mir = brick(system.file("lucc_MT/data/mir.tif", package="dtwSat"))
+doy = brick(system.file("lucc_MT/data/doy.tif", package="dtwSat"))
+timeline = scan(system.file("lucc_MT/data/timeline", package="dtwSat"), what="date")
+rts = twdtwRaster(evi, ndvi, red, blue, nir, mir, timeline = timeline, doy = doy)
+
+time_interval = seq(from=as.Date("2007-09-01"), to=as.Date("2013-09-01"), 
+                    by="12 month")
+log_fun = weight.fun=logisticWeight(-0.1,50)
+
+r_twdtw = twdtwApply(x=rts, y=patt, weight.fun=log_fun, breaks=time_interval, 
+          filepath="~/test_twdtw", overwrite=TRUE, format="GTiff")
+
+plot(r_twdtw, type="distance")
+
+r_lucc = twdtwClassify(r_twdtw, format="GTiff", overwrite=TRUE)
+
+plot(r_lucc)
+
+plot(r_lucc, type="distance")
+
+}
+}
+\references{
+[1] Maus  V,  Camara  G,  Cartaxo  R,  Sanchez  A,  Ramos  FM,  de Queiroz, GR.
+(2016). A Time-Weighted Dynamic Time Warping method for land use and land cover 
+mapping. IEEE Journal of Selected Topics in Applied Earth Observations and Remote 
+Sensing, vol.9, no.8, pp.3729-3739.
+
+[2] Giorgino, T. (2009). Computing and Visualizing Dynamic Time Warping Alignments in R: 
+The dtw Package. Journal of Statistical Software, 31, 1-24.
+
+[3] Muller, M. (2007). Dynamic Time Warping. In Information Retrieval for Music 
+and Motion (pp. 79-84). London: Springer London, Limited.
+}
+\seealso{
+\code{\link[dtwSat]{twdtwRaster-class}}, and 
+\code{\link[dtwSat]{createPatterns}}
+}
+\author{
+Victor Maus, \email{vwmaus1@gmail.com}
+}