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).
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_defsThe multi-band image time series has 7 bands and 137 observations ranging from 2007-09-14 to 2013-08-29.
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-fallowWe 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.
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.084Note 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)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).