In vwmaus/dtwSat: Time-Weighted Dynamic Time Warping for Satellite Image Time Series Analysis

knitr::opts_chunk$set(echo = TRUE, collapse = TRUE, dev = "png")

This vignette offers a concise guide for using version 1.0.0 or higher of the dtwSat package to generate a land-use map. The package utilizes Time-Weighted Dynamic Time Warping (TWDTW) along with a 1-Nearest Neighbor (1-NN) classifier. The subsequent sections will walk you through the process of creating a land-use map based on a set of training samples and a multi-band satellite image time series.

Reading Training Samples

First, let's read a set of training samples that come with the dtwSat package installation.

library(dtwSat)

samples <- st_read(system.file("mato_grosso_brazil/samples.gpkg", package = "dtwSat"), quiet = TRUE)

Preparing the Satellite Image Time Series

The dtwSat package supports satellite images read into R using the stars package. The installation comes with a set of MOD13Q1 images for a region within the Brazilian Amazon. Note that timing is crucial for the TWDTW distance metric. To create a consistent image time series, we start by extracting the date of acquisition from the MODIS file names [@Didan:2015].

tif_files <- dir(system.file("mato_grosso_brazil", package = "dtwSat"), pattern = "\\.tif$", full.names = TRUE)

acquisition_date <- as.Date(regmatches(tif_files, regexpr("[0-9]{8}", tif_files)), format = "%Y%m%d")

print(acquisition_date)

Side note: The date in the file name is not the true acquisition date for each pixel. MOD13Q1 is a 16-day composite product, and the date in the file name is the first day of this 16-day period.

With the files and dates in hand, we can construct a stars satellite image time series for dtwSat.

# read data-cube
dc <- read_stars(tif_files,
                 along = list(time = acquisition_date),
                 RasterIO = list(bands = 1:6))

# set band names
dc <- st_set_dimensions(dc, 3, c("EVI", "NDVI", "RED", "BLUE", "NIR", "MIR"))

# convert band dimension to attribute
dc <- split(dc, c("band"))

print(dc)

Note that it's important to set the date for each observation using the parameter along. This will produce a 4D array (data-cube) that will be collapsed into a 3D array by converting the 'band' dimension into attributes. This prepares the data for training the TWDTW-1NN model.

Create TWDTW-KNN1 model

There several wayis to buld the TWDTW-1NN model. The default options is to keep all samples as part of the model. However, this implies in higher computational time for prediction. To speed the precessing, here I show how to reduce the samples to a single pattern for land use class using Generalized Additive Models (GAM) with cubic regression splines. That can be achieved by defining a smoothing function, such as

library(mgcv)

twdtw_model <- twdtw_knn1(x = dc,
                          y = samples,
                          resampling_fun = function(data) mgcv::gam(y ~ s(x), data = data),
                          cycle_length = 'year',
                          time_scale = 'day',
                          time_weight = c(steepness = 0.1, midpoint = 50))

print(twdtw_model)

In addition to the mandatory arguments x (satellite data-cube) and y (training samples), the TWDTW distance calculation also requires setting cycle_length, time_scale, and time_weight. For more details, refer to the documentation using ?twdtw. The model will resample the training set using the defined function, collapsing all samples with the same land-use label into a single sample. This reduces computational demands. The sample in the model can be visualized as follows:

plot(twdtw_model)

Land Use Prediction

Finally, we predict the land-use classes for each pixel location in the data-cube:

lu_map <- predict(dc, model = twdtw_model)
print(lu_map)

The 'time' dimension was reduced to a single map. We can now visualize it using ggplot:

ggplot() +
  geom_stars(data = lu_map) +
  theme_minimal()

Note that some pixels (in a 3x3 box) in the northeast part of the map have NA values due to a missing value in the blue band recorded on 2011-11-17. This limitation will be addressed in future versions of the dtwSat package.

Ultimately, we can write the map to a TIFF file we can use:

write_stars(lu_map, "lu_map.tif")

Further Reading

This introduction outlined the use of dtwSat for land-use mapping. For more in-depth information, refer to the papers by @Maus:2016 and @Maus:2019 and the twdtw R package documentation.

For additional details on how to manage input and output satellite images, check the stars documentation.

References

vwmaus/dtwSat documentation built on Jan. 28, 2024, 9 a.m.