18  Classification of raster data cubes

Configurations to run this chapter

# load packages "torch" and "luz"
library(torch)
library(luz)
torch::install_torch()
# load packages "sits" and "sitsdata"
library(sits)
# set tempdir if it does not exist 
tempdir_r <- "~/sitsbook/tempdir/R/cl_rasterclassification"
dir.create(tempdir_r, showWarnings = FALSE)
# Create a directory to store the data
tempdir_data <- "~/sitsbook/tempdir/R/data/Rondonia20LMR"
dir.create(tempdir_data, recursive = TRUE, showWarnings = FALSE)
# load "pysits" library
from pysits import *
from pathlib import Path
# set tempdir if it does not exist 
tempdir_py = Path.home() / "sitsbook/tempdir/Python/cl_rasterclassification"
tempdir_py.mkdir(parents=True, exist_ok=True)
# Create a directory to store the data
tempdir_data = Path.home() /"sitsbook/tempdir/data/Rondonia20LMR"
tempdir_data.mkdir(parents=True, exist_ok=True)

This chapter discusses how to classify data cubes by providing a step-by-step example. Our study area is the state of Rondonia, Brazil, which underwent substantial deforestation in the last decades. The objective of the case study is to detect deforested areas.

18.1 Data cube for case study

The examples of this chapter use a pre-built data cube of Sentinel-2 images, available in the package sitsdata. These images are from the SENTINEL-2-L2A collection in Microsoft Planetary Computer (MPC). The data consists of bands B02, B8A, and B11, and indexes NDVI, EVI, and NBR in a small area of \(1200 \times 1200\) pixels in the state of Rondonia. As explained in Chapter Data cubes from local files, we need to inform sits how to parse these file names to obtain tile, date, and band information. Image files are named according to the convention “satellite_sensor_tile_band_date” (e.g., SENTINEL-2_MSI_20LKP_B02_2020_06_04.tif) which is the default format in sits.

# Read data cube from Hugging Face
rondonia_20LMR <- sits_from_hf(
    repo = "e-sensing/book_Rondonia20LMR",
    type = "dataset",
    output_dir = tempdir_data
)

# Plot the cube
plot(rondonia_20LMR, date = "2022-07-16", band = "NDVI")

Figure 18.1: Color composite image of the cube for date 2022-07-16.
# Read data cube from Hugging Face
rondonia_20LMR <- sits_from_hf(
    repo = "e-sensing/book_Rondonia20LMR",
    type = "dataset",
    output_dir = tempdir_data
)

# Plot the cube
plot(rondonia_20LMR, date = "2022-07-16", band = "NDVI")

Color composite image of the cube for date 2022-07-16.

18.2 Training data for the case study

This case study uses the training dataset samples_deforestation_rondonia, available in package sitsdata. This dataset consists of 6,007 samples collected from Sentinel-2 images covering the state of Rondonia. There are nine classes: Clear_Cut_Bare_Soil, Clear_Cut_Burned_Area, Mountainside_Forest, Forest, Riparian_Forest, Clear_Cut_Vegetation, Water, Wetland, and Seasonally_Flooded. Each time series contains values from Sentinel-2/2A bands B02, B03, B04, B05, B06, B07, B8A, B08, B11 and B12, from 2022-01-05 to 2022-12-23 in 16-day intervals. The samples are intended to detect deforestation events and have been collected by remote sensing experts using visual interpretation.

# Get the samples from Hugging Face
samples_deforestation_rondonia <- sits_from_hf(
  repo = "e-sensing/samples_book",
  file = "samples_deforestation_rondonia.parquet",
  type = "dataset"
)

# Show the contents of the samples
summary(samples_deforestation_rondonia)
# A tibble: 9 × 3
  label                 count   prop
  <chr>                 <int>  <dbl>
1 Clear_Cut_Bare_Soil     944 0.157 
2 Clear_Cut_Burned_Area   983 0.164 
3 Clear_Cut_Vegetation    603 0.100 
4 Forest                  964 0.160 
5 Mountainside_Forest     211 0.0351
6 Riparian_Forest        1247 0.208 
7 Seasonally_Flooded      731 0.122 
8 Water                   109 0.0181
9 Wetland                 215 0.0358
# Obtain the samples 
samples_deforestation_rondonia = sits_from_hf(
  repo = "e-sensing/samples_book",
  file = "samples_deforestation_rondonia.parquet",
  type = "dataset"
)

# Show the contents of the samples
summary(samples_deforestation_rondonia)
# A tibble: 9 × 3
  label                 count   prop
  <chr>                 <int>  <dbl>
1 Clear_Cut_Bare_Soil     944 0.157 
2 Clear_Cut_Burned_Area   983 0.164 
3 Clear_Cut_Vegetation    603 0.100 
4 Forest                  964 0.160 
5 Mountainside_Forest     211 0.0351
6 Riparian_Forest        1247 0.208 
7 Seasonally_Flooded      731 0.122 
8 Water                   109 0.0181
9 Wetland                 215 0.0358

It is helpful to plot the basic patterns associated with the samples to understand the training set better. The function sits_patterns() uses a generalized additive model (GAM) to predict a smooth, idealized approximation to the time series associated with each class for all bands. Since the data cube used in the classification has 10 bands, we obtain the indexes NDVI, EVI, and NBR before showing the patterns.

# Generate indexes 
samples_deforestation_indices <- samples_deforestation_rondonia |> 
    sits_apply(NDVI = (B08 - B04)/(B08 + B04)) |> 
    sits_apply(NBR = (B08 - B12) / (B08 + B12)) |> 
    sits_apply(EVI = 2.5 * (B08 - B04) / ((B08 + 6.0 * B04 - 7.5 * B02) + 1.0)) 

# Generate and plot patterns
samples_indices <- sits_select(
  samples_deforestation_indices,
  bands = c("NDVI", "EVI", "NBR")
)
patterns_20LMR <- sits_patterns(samples_indices)
plot(patterns_20LMR)

Figure 18.2: Time series patterns for deforestation study in Rondonia.
# Generate indexes 
samples_deforestation_indices = sits_apply(
    data = samples_deforestation_rondonia,
    NDVI = "(B08 - B04)/(B08 + B04)",
    NBR = "(B08 - B12) / (B08 + B12)",
    EVI = "2.5 * (B08 - B04) / ((B08 + 6.0 * B04 - 7.5 * B02) + 1.0)"
)

# Generate patterns
patterns_20LMR = sits_patterns(
    sits_select(
        data = samples_deforestation_indices, 
        bands = ("NDVI", "EVI", "NBR")
    )
)

# Plot patterns
plot(patterns_20LMR)

Time series patterns for deforestation study in Rondonia.

The patterns show different temporal responses for the selected classes. They match the typical behavior of deforestation in the Amazon. In most cases, the forest is cut at the start of the dry season (May/June). At the end of the dry season, some clear-cut areas are burned to clean the remains; this action is reflected in the steep fall of the response of B11 values of burned area samples after August. The areas where native trees have been cut but some vegetation remains (“Clear_Cut_Vegetation”) have values in the B8A band that increase during the period.

18.3 Training machine learning models

The next step is to train a machine learning model to illustrate CPU-based classification. We build a Random Forest model using sits_train() and then create a plot to find out what are the most important variables for the model.

# set the seed to get the same result
set.seed(03022024)

# Train model using Random Forest algorithm
rfor_model <- sits_train(
  samples_deforestation_rondonia,
  ml_method = sits_rfor()
)

# plot the model results
plot(rfor_model)

Figure 18.3: Most relevant variables of the Random Forest model.
# set the seed to get the same result
r_set_seed(3022024)

# Train model using Random Forest algorithm
rfor_model = sits_train(
  samples_deforestation_rondonia,
  ml_method = sits_rfor()
)

# plot the model results
plot(rfor_model)

Most relevant variables of the Random Forest model.

The figure shows bands and dates that represent relevant inflection points in the image time series.

18.4 Classification of machine learning models in CPUs

By default, all classification algorithms in sits use CPU-based parallel processing, done internally by the package. The algorithms are adaptable; the only requirement for users is to inform the configuration of their machines. To achieve efficiency, sits implements a fault-tolerant multitasking procedure, using a cluster of independent workers linked to a virtual machine. To avoid communication overhead, all large payloads are read and stored independently; direct interaction between the main process and the workers is kept at a minimum. Details of CPU-based parallel processing in sits can be found in the Chapter How parallel processing works in SITS.

To classify both data cubes and sets of time series, use sits_classify(), which uses parallel processing to speed up the performance, as described at the end of this chapter. Its parameters are:

Parameter Meaning
data either a data cube or a set of time series (mandatory).
ml_model a trained model that uses one of the machine learning methods provided (mandatory).
impute_fn imputation function used to interpolate NA values (default = impute_linear).
roi Region of interest where pixels will be classified. See details in the function documentation.
exclusion_mask Areas to be excluded from the classification, defined by either an sf object or by a shapefile (optional).
start_date Start date for the classification process (optional).
end_date End date for the classification process (optional).
multicores number of CPU cores that will be used for processing (default = 2).
memsize RAM memory available for classification in GB (default = 4).
gpu_memory if a GPU is available, defines the graphics memory to be used for classification (default = 4).
batch_size size of samples to be sent to GPU for classification (default is 1000 * gpu_memory).
output_dir Directory where results will be stored (mandatory for raster data cubes).
version For version control (default = “v1”).
verbose Detailed information about processing (default = FALSE).
progress progress bar to track processing steps (default = TRUE).

The version parameter distinguishes different variants of the same classification, as in the case of outputs produced by different models.

sits also supports a recovery mode. If a classification is interrupted, running sits_classify() again with the same parameters, including the same output_dir, makes the software reuse the results already computed and resume from the point of interruption.

The classification result is a data cube with a set of probability layers, one for each output class. Each probability layer contains the model’s assessment of how likely each pixel belongs to the related class. The probability cube can be visualized with plot(). In this example, we show only the probabilities associated with label “Forest”.

# Classify data cube to obtain a probability cube
rondonia_20LMR_probs <- sits_classify(
    data     = rondonia_20LMR,
    ml_model = rfor_model,
    output_dir = tempdir_r,
    version = "rf-raster",
    multicores = 4,
    memsize = 16)

plot(rondonia_20LMR_probs, labels = "Forest", palette = "YlGn")

Figure 18.4: Probability map for class Forest for Random Forest model.
# Classify data cube to obtain a probability cube
rondonia_20LMR_probs = sits_classify(
    data     = rondonia_20LMR,
    ml_model = rfor_model,
    output_dir = tempdir_py,
    version = "rf-raster",
    multicores = 4,
    memsize = 16)

plot(rondonia_20LMR_probs, labels = "Forest", palette = "YlGn")

Probability map for class Forest for Random Forest model.

The probability cube provides information on the output values of the algorithm for each class. Most probability maps contain outliers or misclassified pixels. The labeled map generated from the pixel-based time series classification method exhibits several misclassified pixels, which are small patches surrounded by a different class. This occurrence of outliers is a common issue that arises due to the inherent nature of this classification approach. Regardless of their resolution, mixed pixels are prevalent in images, and each class exhibits considerable data variability. As a result, these factors can lead to outliers that are more likely to be misclassified. To overcome this limitation, sits employs post-processing smoothing techniques that leverage the spatial context of the probability cubes to refine the results. These techniques will be discussed in the Chapter Bayesian smoothing for post-processing. In what follows, we will generate the smoothed cube to illustrate the procedure.

The parameters for sits_smooth() are:

Parameter Meaning
cube A data cube with a set of probability matrices (mandatory).
window_size Size of the spatial neighborhood (default = 9 for a 9 x 9 filter).
neigh_fraction Fraction of neighbors with high probabilities to be used in Bayesian inference (default = 0.5).
smoothness Bayesian smoothing parameter (default = 20).
exclusion_mask Areas to be excluded from the classification, defined by either an sf object or by a shapefile (optional).
multicores number of CPU cores that will be used for processing (default = 2).
memsize RAM memory available for classification (default = 4).
output_dir Directory where results will be stored (mandatory).
version For version control (default = “v1”).
progress progress bar to track processing steps (default = TRUE). 
# Smoothen a probability cube
rondonia_20LMR_bayes <- sits_smooth(
    cube    = rondonia_20LMR_probs,
    output_dir = tempdir_r,
    version = "rf-raster",
    multicores = 4,
    memsize = 16)

plot(rondonia_20LMR_bayes, labels = c("Forest"),  palette = "YlGn")

Figure 18.5: Smoothed probability map for class Forest.
# Smoothen a probability cube
rondonia_20LMR_bayes = sits_smooth(
    cube    = rondonia_20LMR_probs,
    output_dir = tempdir_py,
    version = "rf-raster",
    multicores = 4,
    memsize = 16)

plot(rondonia_20LMR_bayes, labels = ("Forest"),  palette = "YlGn")

Smoothed probability map for class Forest.

In general, users should perform a post-processing smoothing after obtaining the probability maps in raster format. After the post-processing operation, we apply sits_label_classification() to obtain a map with the most likely class for each pixel. For each pixel, the sits_label_classification() function takes the label with highest probability and assigns it to the resulting map. The output is a labelled map with classes.

The parameters for sits_label_classification() for raster data cubes are:

Parameter Meaning
cube A data cube with a set of probability matrices (mandatory).
multicores number of CPU cores that will be used for processing (default = 2).
memsize RAM memory available for classification (default = 4).
output_dir Directory where results will be stored (mandatory).
version For version control (default = “v1”).
progress progress bar to track processing steps (default = TRUE). 
# Generate a thematic map
rondonia_20LMR_class <- sits_label_classification(
    cube = rondonia_20LMR_bayes,
    multicores = 4,
    memsize = 12,
    output_dir = tempdir_r,
    version = "rf-raster")

# Plot the thematic map
plot(rondonia_20LMR_class, legend_text_size = 0.7, legend_position = "inside")

Figure 18.6: Final map of deforestation obtained by Random Forest model.
# Generate a thematic map
rondonia_20LMR_class = sits_label_classification(
    cube = rondonia_20LMR_bayes,
    multicores = 4,
    memsize = 12,
    output_dir = tempdir_py,
    version = "rf-raster")

# Plot the thematic map
plot(rondonia_20LMR_class, legend_text_size = 0.7, legend_position = "inside")

Final map of deforestation obtained by Random Forest model.

18.5 Training and running deep learning models

The next examples show how to run deep learning models in sits. The case study uses the Temporal CNN model [1], which is described in Chapter Machine learning algorithms for image time series. Deep learning time series classification methods in sits, which include sits_lighttae(), sits_lstm_fcn(), sits_mlp(), sits_tempcnn(), and sits_resnet(), are written using the torch package, which is an adaptation of PyTorch to the R environment.

Deep learning algorithms use a CUDA-compatible NVIDIA GPU if one is available and has been properly configured. Please refer to the torch installation guide for details on how to configure torch to use CUDA GPUs. Macs with an integrated MPS architecture are also supported by sits; in this case, no special configuration is required.

If no GPU is available, these algorithms will run on regular CPUs, using the same parallelization methods described in the traditional machine learning methods. There is a relevant performance increase when running torch-based methods in GPUs relative to their processing time on CPUs.

We take the same data cube and training data used in the previous examples and use a Temporal CNN method. The first step is to obtain a deep learning model using the sits_tempcnn() algorithm. We use the hyperparameters obtained by sits_tuning() in the previous chapter. Based on our experience with sits_tempcnn(), we prefer smaller kernels due to the size of the time series and thus we chose the second-best parameter combination, which uses kernels of size 3.

tcnn_model <- sits_train(
  samples_deforestation_rondonia,
  sits_tempcnn(
    cnn_layers = c(128, 128, 128),
    cnn_kernels = c(3, 3, 3),
    cnn_dropout_rates = c(0.30, 0.30, 0.30),
    opt_hparams = list(
      lr = 7.18e-04,
      eps = 1.0e-08,
      weight_decay = 1.1e-07
    )
  )
)
tcnn_model = sits_train(
  samples_deforestation_rondonia,
  sits_tempcnn(
    cnn_layers = (128, 128, 128),
    cnn_kernels = (3, 3, 3),
    cnn_dropout_rates = (0.30, 0.30, 0.30),
    opt_hparams = dict(
      lr = 7.18e-04,
      eps = 1.0e-08,
      weight_decay = 1.1e-07
      )
  )
)

After training the model, we classify the data cube. If a GPU is available, users need to provide the additional parameters gpu_memory and batch_size. Batch size is the number of training data samples processed by the neural network at each pass before the model updates its internal weights. Instead of feeding the whole dataset in at once or processing just one sample at a time, data is split into small groups or batches. This information will be used by sits to optimize access to the GPU and speed up processing. In the classification phase, the model weights are fixed, so we recommend using a large value for batch_size.

rondonia_20LMR_probs_tcnn <- sits_classify(
    rondonia_20LMR,
    ml_model = tcnn_model,
    output_dir = tempdir_r,
    version = "tcnn-raster",
    gpu_memory = 16,
    batch_size = 2048,
    multicores = 6,
    memsize = 24
)

plot(rondonia_20LMR_probs_tcnn, labels = c("Forest"),  palette = "YlGn")

Figure 18.7: Probability map for class Forest for TempCNN model.
# ensure model and data cube use the same bands for classification
rondonia_20LMR = sits_select(rondonia_20LMR, bands = sits_bands(tcnn_model))
#| 
# Classify the data cube
rondonia_20LMR_probs_tcnn = sits_classify(
    rondonia_20LMR,
    ml_model = tcnn_model,
    output_dir = tempdir_py,
    version = "tcnn-raster",
    gpu_memory = 16,
    batch_size = 2048,
    multicores = 6,
    memsize = 24
)

plot(rondonia_20LMR_probs_tcnn, labels = ("Forest"),  palette = "YlGn")

Probability map for class Forest for TempCNN model.

After classification, we can smooth the probability cube. It is useful to compare the smoothed map for class Forest resulting from the TempCNN model with that produced by the Random Forest algorithm. The TempCNN model tends to show more confidence in its predictions than the Random Forest one. This is a feature of the model more than an intrinsic property of the training data or the data cube.

# Smoothen the probability map 
rondonia_20LMR_bayes_tcnn <- sits_smooth(
    rondonia_20LMR_probs_tcnn,
    output_dir = tempdir_r,
    version = "tcnn-raster",
    multicores = 6,
    memsize = 24
)

plot(rondonia_20LMR_bayes_tcnn, labels = c("Forest"),  palette = "YlGn")

Figure 18.8: Smoothed probability map for class Forest for TempCNN model.
# Smoothen the probability map 
rondonia_20LMR_bayes_tcnn = sits_smooth(
    rondonia_20LMR_probs_tcnn,
    output_dir = tempdir_py,
    version = "tcnn-raster",
    multicores = 6,
    memsize = 24
)

plot(rondonia_20LMR_bayes_tcnn, labels = ("Forest"),  palette = "YlGn")

Smoothed probability map for class Forest for TempCNN model.

We then label the resulting smoothed probabilities to obtain a classified map.

# Obtain the final labelled map
rondonia_20LMR_class_tcnn <- sits_label_classification(
    rondonia_20LMR_bayes_tcnn,
    output_dir = tempdir_r,
    version = "tcnn-raster",
    multicores = 6,
    memsize = 24
)

plot(rondonia_20LMR_class_tcnn, legend_text_size = 0.7, legend_position = "inside")

Figure 18.9: Classified map of data cube for TempCNN model.
# Obtain the final labelled map
rondonia_20LMR_class_tcnn = sits_label_classification(
    rondonia_20LMR_bayes_tcnn,
    output_dir = tempdir_py,
    version = "tcnn-raster",
    multicores = 6,
    memsize = 24
)

plot(rondonia_20LMR_class_tcnn)

Classified map of data cube for TempCNN model.

18.6 Comparison between classifications using TempCNN and RandomForest

A detailed accuracy assessment of the maps produced by TempCNN and RandomForest is outside the scope of this chapter. Please see the chapter on “Map accuracy assessment” for a discussion on the subject. A qualitative assessment of the two maps (Figure 18.6 and Figure 18.9) shows that TempCNN is better at classifying subtle changes. For example, TempCNN detects more areas labelled as Clear_Cut_Burned_Area. This comparison is compatible with the published literature [2].

18.7 Summary

This chapter presents a detailed example of how to train models and apply them to raster classification in sits. The procedure is simple and direct, using a workflow that combines sits_train(), sits_classify(), sits_smooth() and sits_label_classification(). For traditional machine learning models, such as Random Forest, the code is optimized for CPU processing and users will see good performance. In the case of deep learning models, sits is optimized for GPU processing. The typical time for classifying a 10-band Sentinel-2 tile in a CPU using Random Forest is 20 minutes in a moderately sized machine (16 cores, 64 GB RAM). The same performance is expected when running deep learning models in a standard GPU with 16 GB memory.

References

[1]
C. Pelletier, G. I. Webb, and F. Petitjean, “Temporal Convolutional Neural Network for the Classification of Satellite Image Time Series,” Remote Sensing, vol. 11, no. 5, 2019.
[2]
M. Rußwurm, C. Pelletier, M. Zollner, S. Lefèvre, and M. Körner, “BreizhCrops: A Time Series Dataset for Crop Type Mapping,” 2020, [Online]. Available: http://arxiv.org/abs/1905.11893.