The opening of large archives of satellite data such as Landsat, MODIS and the Sentinels has given researchers unprecedented access to data, allowing them to better quantify and understand local and global land change. The need to analyse such large data sets has led to the development of automated and semi-automated methods for satellite image time series analysis. However, few of the proposed methods for remote sensing time series analysis are available as open source software. The package dtwSat provides an implementation of the Time-Weighted Dynamic Time Warping (TWDTW) method for land cover mapping using multi-band satellite image time series (Maus et al. 2016). Methods based on dynamic time warping are flexible to handle irregular sampling and out-of-phase time series, and they have achieved significant results in time series analysis (Velichko and Zagoruyko 1970; Hiroaki Sakoe and Chiba 1971; H. Sakoe and Chiba 1978; Rabiner and Juang 1993; Berndt and Clifford 1994; Keogh and Ratanamahatana 2005; Müller 2007). dtwSat is also available from the Comprehensive R Archive Network (CRAN). dtwSat provides full cycle of land cover classification using image time series, ranging from selecting temporal patterns to visualising, and assessing the results. Bellow we show a quick demo of the package usage.

The GitHub version requires the package devtools

install.packages("devtools")
devtools::install_github("vwmaus/dtwSat")

In this quick demo we will perform a TWDTW analysis for a single time series. The data for the analysis are a set of temporal patterns in MOD13Q1.patterns.list and an example of time series in MOD13Q1.ts. These time series are in zoo format and come with the package installation. Suppose that we want to know the crop type of each subinterval in following time series:

library(dtwSat)
# Create and plot object time series 
ts <- twdtwTimeSeries(MOD13Q1.ts)
class(ts)
plot(ts, type = "timeseries")

Fig. 1. Example time series which we want to classify.

We know that in the region where the time series was observed we have soybean, cotton, and maize, whose typical temporal pattern are:

# Create and plot object time series 
patt <- twdtwTimeSeries(MOD13Q1.patterns.list)
class(patt)
plot(patt, type = "patterns") 

Fig. 2. Typical temporal patterns of *soybean*, *cotton*, and *maize*.

Using these temporal patterns we run the TWDTW analysis, such that

# Define logistic time-weight, see Maus et al. (2016)
log_fun <- logisticWeight(alpha = -0.1, beta = 100) 
# Run TWDTW analysis 
matches <- twdtwApply(x = ts, y = patt, weight.fun = log_fun, keep = TRUE) 

The result is a twdtwMatches object with all possible matches of the patterns to the time series

class(matches)

## [1] "twdtwMatches"
## attr(,"package")
## [1] "dtwSat"

show(matches)

## An object of class "twdtwMatches"
## Number of time series: 1 
## Number of Alignments: 27 
## Patterns labels: Soybean Cotton Maize

We can use several plot methods to visualize the results of the analysis in the twdtwMatches object, for example, to plot the alignments

plot(x = matches, type = "alignments")

Fig. 3. TWDTW alignments over time and cost (distance) in y-axis.

to plot matching point

plot(x = matches, type = "matches", attr = "evi", patterns.labels = "Soybean", k <- 4) 

Fig. 4. The four best matches of *soybean*.

to plot minimum cost paths

plot(x = matches, type = "paths", k <- 1:4) 

Fig. 1. The minimum cost path of the TWDTW alignment for each crop type.

and, finally to classify the subintervals of the time series. The plot will select the best match for each period of 6 months, i.e. the class for each period.

plot(x = matches, type = "classification",
     from = "2009-09-01", to = "2013-09-01", 
     by = "6 month", overlap = 0.5) 

Fig. 2. Classification using the best match for each subinterval.

The next example shows how to classify a raster time series, i.e. the same as we did in the quick demo but now for each pixel location. For that we use a set of MODIS (MOD13Q1 product) images from 2007 to 2013 for a region in the Brazilian Amazon. These data is included in the package installation. Load raster time series:

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"))

Load the dates of the MODIS images:

timeline <- scan(system.file("lucc_MT/data/timeline", package = "dtwSat"), what = "date")

Build raster time series:

rts <- twdtwRaster(evi, ndvi, red, blue, nir, mir, timeline = timeline, doy = doy)

Load the set of ground truth samples and projection information:

field_samples <- read.csv(system.file("lucc_MT/data/samples.csv", package = "dtwSat"))
proj_str <- scan(system.file("lucc_MT/data/samples_projection", package = "dtwSat"), what = "character")

We use the package caret to split the samples into training (10%) and validation (90%)

library(caret)
set.seed(1)
I <- unlist(createDataPartition(field_samples$label, p = 0.1))
training_samples <- field_samples[I, ]
validation_samples <- field_samples[-I, ]

Extract training time series from raster time series

training_ts <- getTimeSeries(rts, y = training_samples, proj4string = proj_str)
validation_ts <- getTimeSeries(rts, y = validation_samples, proj4string = proj_str)

Create temporal patterns using training samples

temporal_patterns <- createPatterns(training_ts, freq = 8, formula = y ~ s(x))

plot(temporal_patterns, type = "patterns") 

Fig. 3. Typical temporal patterns of *Cotton-fallow*, *Forest*, *Soybean-cotton*, *Soybean-maize*, and *Soybean-millet*.

Apply TWDTW analysis:

# Define logistic time-weight, see Maus et al. (2016)
log_fun <- logisticWeight(-0.1, 50)

# Run serial TWDTW analysis 
r_twdtw <- twdtwApply(x = rts, y = temporal_patterns, weight.fun = log_fun, progress = 'text')

# or Run parallel TWDTW analysis
beginCluster()
r_twdtw <- twdtwApplyParallel(x = rts, y = temporal_patterns, weight.fun = log_fun, progress = 'text')
endCluster()

Classify raster raster time series using the results from the TWDTW analysis

r_lucc <- twdtwClassify(r_twdtw, progress = 'text')

Visualising the results.

Land cover maps

plot(x = r_lucc, type = "maps")

Fig. 4. Land cover maps based on TWDTW analysis.

Land cover area for each class over time

plot(x = r_lucc, type = "area")

Fig. 5. Land cover area based on TWDTW analysis.

Land cover changes over time (gains and losses from/to classes)

plot(x = r_lucc, type = "changes")

Fig. 6. Land cover changes based on TWDTW analysis.

We use the validation samples to compute the metrics for accuracy assessment.

twdtw_assess <- twdtwAssess(object = r_lucc, y = validation_samples, 
  proj4string = proj_str, conf.int = .95) 
show(twdtw_assess)

## An object of class "twdtwAssessment"
## Number of classification intervals: 6 
## Accuracy metrics summary
## 
## Overall
## Accuracy      Var       sd      ci* 
##  9.8e-01  5.8e-05  7.6e-03  1.5e-02 
## 
## User's
##                Accuracy     Var    sd   ci*
## Cotton-fallow      0.95 0.00071 0.027 0.052
## Forest             1.00 0.00000 0.000 0.000
## Soybean-cotton     1.00 0.00000 0.000 0.000
## Soybean-maize      0.95 0.00036 0.019 0.037
## Soybean-millet     1.00 0.00000 0.000 0.000
## unclassified       1.00 0.00000 0.000 0.000
## 
## Producer's
##                Accuracy    Var    sd  ci*
## Cotton-fallow      1.00 0.0000 0.000 0.00
## Forest             1.00 0.0000 0.000 0.00
## Soybean-cotton     0.72 0.0044 0.067 0.13
## Soybean-maize      1.00 0.0000 0.000 0.00
## Soybean-millet     1.00 0.0000 0.000 0.00
## unclassified       1.00 0.0000 0.000 0.00
## 
## Area and uncertainty
##                 Mapped Adjusted     ci*
## Cotton-fallow  4.8e+07  4.5e+07 2484480
## Forest         7.5e+07  7.5e+07       0
## Soybean-cotton 1.9e+07  2.6e+07 4806920
## Soybean-maize  1.1e+08  1.0e+08 4115074
## Soybean-millet 7.0e+07  7.0e+07       0
## unclassified   0.0e+00  0.0e+00       0
## 
## * 95 % confidence interval

Visualizing User's and Producer's accuracy

plot(twdtw_assess, type = "accuracy")

Fig. 7. User's and Producer's accuracy.

Visualizing area uncertainty

plot(twdtw_assess, type = "area")

Fig. 8. Area uncertainty.

For further discussion on the package see the vignettes and if you want to learn more about the TWDTW method (see, Maus et al. 2016).

dtwSat: Time-Weighted Dynamic Time Warping for Satellite Image Time Series Analysis in R

Berndt, Donald J., and James Clifford. 1994. “Using Dynamic Time Warping to Find Patterns in Time Series.” In KDD Workshop, edited by Usama M. Fayyad and Ramasamy Uthurusamy, 359–70. AAAI Press.

Keogh, Eamonn, and Chotirat Ann Ratanamahatana. 2005. “Exact Indexing of Dynamic Time Warping.” Knowledge Information Systems 7 (3): 358–86.

Maus, Victor, Gilberto Camara, Ricardo Cartaxo, Alber Sanchez, Fernando M. Ramos, and Gilberto R. de Queiroz. 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 9 (8): 3729–39. doi:10.1109/JSTARS.2016.2517118.

Müller, Meinard. 2007. Information Retrieval for Music and Motion. London: Springer-Verlag.

Rabiner, Lawrence, and Biing-Hwang Juang. 1993. Fundamentals of Speech Recognition. New Jersey: Prentice-Hall International, Inc.

Sakoe, H., and S. Chiba. 1978. “Dynamic Programming Algorithm Optimization for Spoken Word Recognition.” IEEE Transactions on Acoustics, Speech, and Signal Processing 26 (1): 43–49. doi:10.1109/TASSP.1978.1163055.

Sakoe, Hiroaki, and Seibi Chiba. 1971. “A Dynamic Programming Approach to Continuous Speech Recognition.” In Proceedings of the Seventh International Congress on Acoustics, Budapest, 3:65–69. Budapest: Akadémiai Kiadó.

Velichko, V.M., and N.G. Zagoruyko. 1970. “Automatic Recognition of 200 Words.” International Journal of Man-Machine Studies 2 (3): 223–34. doi:10.1016/S0020-7373(70)80008-6.