2. Land-cover change analysis with TWDTW

Victor Maus

This vignette provides a short guide on how to perform a Time-Weighted Dynamic Time Warping (TWDTW) analysis on raster image time series using dtwSat. For more details about TWDTW read Maus et al. (2016) and Maus et al. (2019).

Create multi-band image time series

dtwSat provides a set of images extracted from MOD13Q1 dataset for a region within the Brazilian Amazon. You can load these images with

library(dtwSat)  
  
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"))

The tif files above have been pre-processed, so that a single file has the complete time series for each band. One can also built the time series for each band from independent files using the function stack from the R package raster.

Note that every band must have the same extend, resolution, and number of images, i.e. the same number of observations across space and time. To create the multi-band image time series we need the dates of each observation, which for the set of images above are

timeline <- scan(system.file("lucc_MT/data/timeline", package = "dtwSat"), what="date")
timeline
##   [1] "2007-09-14" "2007-09-30" "2007-10-16" "2007-11-01" "2007-11-17"
##   [6] "2007-12-03" "2007-12-19" "2008-01-01" "2008-01-17" "2008-02-02"
##  [11] "2008-02-18" "2008-03-05" "2008-03-21" "2008-04-06" "2008-04-22"
##  [16] "2008-05-08" "2008-05-24" "2008-06-09" "2008-06-25" "2008-07-11"
##  [21] "2008-07-27" "2008-08-12" "2008-08-28" "2008-09-13" "2008-09-29"
##  [26] "2008-10-15" "2008-10-31" "2008-11-16" "2008-12-02" "2008-12-18"
##  [31] "2009-01-01" "2009-01-17" "2009-02-02" "2009-02-18" "2009-03-06"
##  [36] "2009-03-22" "2009-04-07" "2009-04-23" "2009-05-09" "2009-05-25"
##  [41] "2009-06-10" "2009-06-26" "2009-07-12" "2009-07-28" "2009-08-13"
##  [46] "2009-08-29" "2009-09-14" "2009-09-30" "2009-10-16" "2009-11-01"
##  [51] "2009-11-17" "2009-12-03" "2009-12-19" "2010-01-01" "2010-01-17"
##  [56] "2010-02-02" "2010-02-18" "2010-03-06" "2010-03-22" "2010-04-07"
##  [61] "2010-04-23" "2010-05-09" "2010-05-25" "2010-06-10" "2010-06-26"
##  [66] "2010-07-12" "2010-07-28" "2010-08-13" "2010-08-29" "2010-09-14"
##  [71] "2010-09-30" "2010-10-16" "2010-11-01" "2010-11-17" "2010-12-03"
##  [76] "2010-12-19" "2011-01-01" "2011-01-17" "2011-02-02" "2011-02-18"
##  [81] "2011-03-06" "2011-03-22" "2011-04-07" "2011-04-23" "2011-05-09"
##  [86] "2011-05-25" "2011-06-10" "2011-06-26" "2011-07-12" "2011-07-28"
##  [91] "2011-08-13" "2011-08-29" "2011-09-14" "2011-09-30" "2011-10-16"
##  [96] "2011-11-01" "2011-11-17" "2011-12-03" "2011-12-19" "2012-01-01"
## [101] "2012-01-17" "2012-02-02" "2012-02-18" "2012-03-05" "2012-03-21"
## [106] "2012-04-06" "2012-04-22" "2012-05-08" "2012-05-24" "2012-06-09"
## [111] "2012-06-25" "2012-07-11" "2012-07-27" "2012-08-12" "2012-08-28"
## [116] "2012-09-13" "2012-09-29" "2012-10-15" "2012-10-31" "2012-11-16"
## [121] "2012-12-02" "2012-12-18" "2013-01-01" "2013-01-17" "2013-02-02"
## [126] "2013-02-18" "2013-03-06" "2013-03-22" "2013-04-07" "2013-04-23"
## [131] "2013-05-09" "2013-05-25" "2013-06-10" "2013-06-26" "2013-07-12"
## [136] "2013-08-13" "2013-08-29"

Finally, we can create the multi-band image time series with

rts <- twdtwRaster(evi, ndvi, red, blue, nir, mir, timeline = timeline, doy = doy)
rts
## An object of class "twdtwRaster"
## Time series layers: doy evi ndvi red blue nir mir 
## Time range: 2007-09-14 ... 2013-08-29 
## Dimensions: 7 27 37 137 (nlayers, nrow, ncol, length)
## Resolution: 231.6564 231.6564  (x, y)
## Extent    : -6089551 -6080979 -1339205 -1332951 (xmin, xmax, ymin, ymax)
## Coord.ref.: +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs

The multi-band image time series has 7 bands and 137 observations ranging from 2007-09-14 to 2013-08-29.

Create a library of vegetation profiles

To crate a library of profiles we need a set of samples. dtwSat provides a set of samples of land cover classes for the study area. We can load the samples with

field_samples <- read.csv(system.file("lucc_MT/data/samples.csv", package = "dtwSat"))
head(field_samples)
##   longitude  latitude       from         to         label
## 1 -55.98819 -12.03646 2011-09-01 2012-09-01 Cotton-fallow
## 2 -55.99118 -12.04062 2011-09-01 2012-09-01 Cotton-fallow
## 3 -55.98606 -12.03646 2011-09-01 2012-09-01 Cotton-fallow
## 4 -55.98562 -12.03437 2011-09-01 2012-09-01 Cotton-fallow
## 5 -55.98475 -12.03021 2011-09-01 2012-09-01 Cotton-fallow
## 6 -55.98393 -12.03646 2011-09-01 2012-09-01 Cotton-fallow

We split this samples into a training and a validation set, such that

library(caret)

set.seed(1) # set for reproducibility 

I <- unlist(createDataPartition(field_samples$label, p = 0.1))

training_samples <- field_samples[I,]

validation_samples <- field_samples[-I,]

To get the time series for each sample in the training set, we can run

training_ts <- getTimeSeries(rts, 
                             y = training_samples, 
                             proj4string = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

Finally, we can use the training samples to crate the profiles library with

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

plot(profiles_library, type = "patterns")

The function createPatterns using generalized additive models (GAMs) to create the profiles. For details see ?gam from mgcv package.

Classify image time series

system.time(
  twdtw_lucc <- twdtwApply(x = rts, 
                           y = profiles_library, 
                           alpha = -0.1,
                           beta = 50,
                           progress = 'text', 
                           minrows = 30,
                           legacy = FALSE,
                           time.window = TRUE)
)
## 
  |                                                                            
  |                                                                      |   0%
##    user  system elapsed 
##   7.061   0.020   7.084

Note the argument legacy = FALSE. Setting this argument to FALSE considerably improves the performance of the processing as it does not keep any intermediate dataset. Using time.window = TRUE implements a change in the TWDTW algorithm to reduce processing.

One can also run TWDTW in parallel. For a small area, the performance is not much better than the sequential processing. However, for larger areas, the parallel processing can reduce processing time. To run twdtwApply in parallel setup and register a cluster before calling the function, such that

library(doParallel)
library(parallel)
library(foreach)

cl <- makeCluster(detectCores(), type = "FORK")
registerDoParallel(cl)

system.time(
  twdtw_lucc <- twdtwApply(x = rts, 
                           y = profiles_library, 
                           alpha = -0.1,
                           beta = 50,
                           progress = 'text', 
                           minrows = 30,
                           legacy = FALSE,
                           time.window = TRUE)
)

registerDoSEQ()
stopCluster(cl)

Assess classifiaction results

dtwSat provides a set of functions to asses the classification results. For example with

# Plot TWDTW distances for the first year 
  plot(twdtw_lucc, type = "distance", time.levels = 1)

  
# Plot TWDTW classification results 
  plot(twdtw_lucc, type = "map")


# Plot mapped area time series 
  plot(twdtw_lucc, type = "area")

  
# Plot land-cover changes
  plot(twdtw_lucc, type = "changes")

The package also offers a set of methods to assess the classification accuracy and visualize the results.

# Assess classification 
  twdtw_assess <- 
    twdtwAssess(twdtw_lucc, 
                y = validation_samples, 
                proj4string = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",
                conf.int = .95)
  
# Plot map accuracy 
  plot(twdtw_assess, type = "accuracy")

  
# Plot area uncertainty 
  plot(twdtw_assess, type = "area")

  
# Plot misclassified samples  
  plot(twdtw_assess, type = "map", samples = "incorrect")

  
# Get latex table with error matrix 
  twdtwXtable(twdtw_assess, table.type = "matrix")
## % latex table generated in R 4.2.1 by xtable 1.8-4 package
## % Tue Oct 11 17:31:08 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrrr}
##   \hline
##   &\multicolumn{6}{c}{Reference class}&&\\
## \multicolumn{1}{c}{Map class} & Cotton-fallow & Forest & Soybean-cotton & Soybean-maize & Soybean-millet & unclassified & Total & Area & w\\
##  \hline
## Cotton-fallow & 61 & 0 & 3 & 0 & 0 & 0 & 64 & 43039064.00 & 0.134 \\ 
##   Forest & 0 & 124 & 0 & 0 & 0 & 0 & 124 & 73520595.60 & 0.229 \\ 
##   Soybean-cotton & 0 & 0 & 62 & 0 & 0 & 0 & 62 & 22968478.04 & 0.071 \\ 
##   Soybean-maize & 0 & 0 & 6 & 120 & 2 & 0 & 128 & 116505994.93 & 0.362 \\ 
##   Soybean-millet & 0 & 0 & 0 & 0 & 163 & 0 & 163 & 65631889.36 & 0.204 \\ 
##   unclassified & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.00 & 0.000 \\ 
##   Total & 61 & 124 & 71 & 120 & 165 & 0 & 541 & 321666021.93 & 1.000 \\ 
##    \hline
## \end{tabular}
## \end{table}
  
# Get latex table with error accuracy 
  twdtwXtable(twdtw_assess, table.type = "accuracy")
## % latex table generated in R 4.2.1 by xtable 1.8-4 package
## % Tue Oct 11 17:31:08 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{rllllllllll}
##   \hline
##   &\multicolumn{6}{c}{Reference class}&&&&\\
## \multicolumn{1}{c}{Map class} & Cotton-fallow & Forest & Soybean-cotton & Soybean-maize & Soybean-millet & unclassified & Total & User's* & Producers's* & Overall*\\
##  \hline
## Cotton-fallow & 0.13 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.13 & 0.95$\pm$0.05 & 1.00$\pm$0.00 & 0.97$\pm$0.02 \\ 
##   Forest & 0.00 & 0.23 & 0.00 & 0.00 & 0.00 & 0.00 & 0.23 & 1.00$\pm$0.00 & 1.00$\pm$0.00 &  \\ 
##   Soybean-cotton & 0.01 & 0.00 & 0.07 & 0.02 & 0.00 & 0.00 & 0.09 & 1.00$\pm$0.00 & 0.75$\pm$0.12 &  \\ 
##   Soybean-maize & 0.00 & 0.00 & 0.00 & 0.34 & 0.00 & 0.00 & 0.34 & 0.94$\pm$0.04 & 1.00$\pm$0.00 &  \\ 
##   Soybean-millet & 0.00 & 0.00 & 0.00 & 0.01 & 0.20 & 0.00 & 0.21 & 1.00$\pm$0.00 & 0.97$\pm$0.04 &  \\ 
##   unclassified & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 1.00$\pm$-0.00 & 1.00$\pm$0.00 &  \\ 
##   Total & 0.13 & 0.23 & 0.07 & 0.36 & 0.20 & 0.00 & 1.00 &  &  &  \\ 
##    \hline 
## \multicolumn{10}{l}{* 95\% confidence interval.}
## \end{tabular}
## \end{table}
  
# Get latex table with area uncertainty 
  twdtwXtable(twdtw_assess, table.type = "area")
## % latex table generated in R 4.2.1 by xtable 1.8-4 package
## % Tue Oct 11 17:31:08 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlll}
##   \hline
##   \multicolumn{1}{c}{Class} & Mapped area & Adjusted area & Margin of error*\\
##  \hline
## Cotton-fallow & 43039064.00 & 41021607.87 & $\pm$2246395.42 \\ 
##   Forest & 73520595.60 & 73520595.60 & $\pm$0.00 \\ 
##   Soybean-cotton & 22968478.04 & 30447152.68 & $\pm$4836290.48 \\ 
##   Soybean-maize & 116505994.93 & 109224370.25 & $\pm$4904786.98 \\ 
##   Soybean-millet & 65631889.36 & 67452295.53 & $\pm$2512955.54 \\ 
##   unclassified & 0.00 & 0.00 & $\pm$0.00 \\ 
##    \hline 
## \multicolumn{3}{l}{* 95\% confidence interval.}
## \end{tabular}
## \end{table}

This short introduction showed how to use dtwSat for land-cover change analyse. To learn TWDTW read Maus et al. (2016) and Maus et al. (2019).

References

Maus, Victor, Gilberto Camara, Marius Appel, and Edzer Pebesma. 2019. “dtwSat: Time-Weighted Dynamic Time Warping for Satellite Image Time Series Analysis in R.” Journal of Statistical Software 88 (5): 1–31. https://doi.org/10.18637/jss.v088.i05.
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. https://doi.org/10.1109/JSTARS.2016.2517118.