ts_som.qmd

---
title: "Self-organized maps for sample quality control"
format: html
---

<a href="https://www.kaggle.com/code/esensing/self-organized-maps-for-sample-quality-control" target="_blank">
<img src="https://kaggle.com/static/images/open-in-kaggle.svg"/>
</a>

### Configurations to run this chapter{-}

:::{.panel-tabset}
## R
```{r}
#| echo: true
#| eval: true
#| output: false
# load package "tibble"
library(tibble)
# load packages "sits" and "sitsdata"
library(sits)
library(sitsdata)
# set tempdir if it does not exist 
tempdir_r <- "~/sitsbook/tempdir/R/ts_som"
dir.create(tempdir_r, showWarnings = FALSE)
```
## Python
```{python}
#| echo: true
#| eval: true
#| output: 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/ts_som"
tempdir_py.mkdir(parents=True, exist_ok=True)
```
```{python}
#| echo: false
import pandas as pd
pd.set_option('display.max_columns', 7)
pd.set_option('display.width', 600)
pd.set_option('display.max_rows', 10)
```
:::

## Introduction

The `sits` package provides a clustering technique based on self-organizing maps (SOM) as an alternative to hierarchical clustering for quality control of training samples. SOM is a dimensionality reduction technique [@Kohonen1990], where high-dimensional data is mapped into a two-dimensional map, keeping the topological relations between data patterns. As shown in @fig-ts-som2d, the SOM 2D map is composed of units called neurons. Each neuron has a weight vector, with the same dimension as the training samples. At the start, neurons are assigned a small random value and then trained by competitive learning. The algorithm computes the distances of each member of the training set to all neurons and finds the neuron closest to the input, called the best matching unit.

```{r}
#| label: fig-ts-som2d
#| out.width: "90%"
#| echo: false
#| fig.align: "center"
#| fig.cap: "SOM 2D map creation (source: [@Santos2021a])."
knitr::include_graphics("./images/som_structure.png")
```

The input data for quality assessment is a set of training samples, which are high-dimensional data; for example, a time series with 25 instances of 4 spectral bands has 100 dimensions. When projecting a high-dimensional dataset into a 2D SOM map, the units of the map (called neurons) compete for each sample. Each time series will be mapped to one of the neurons. Since the number of neurons is smaller than the number of classes, each neuron will be associated with many time series. The resulting 2D map will be a set of clusters. Given that SOM preserves the topological structure of neighborhoods in multiple dimensions, clusters that contain training samples with a given label will usually be neighbors in 2D space. The neighbors of each neuron of a SOM map provide information on intraclass and interclass variability, which is used to detect noisy samples. The methodology of using SOM for sample quality assessment is discussed in detail in the reference paper [@Santos2021a].

```{r}
#| label: fig-ts-sommet
#| out.width: "90%"
#| echo: false
#| fig.align: "center"
#| fig.cap: "Using SOM for class noise reduction (source: [@Santos2021a])"

knitr::include_graphics("./images/methodology_bayes_som.png")
```

## Dataset used in this chapter

The examples of this chapter use  `samples_cerrado_mod13q1`, a set of time series from the Cerrado region of Brazil. The data ranges from 2000 to 2017 and includes 50,160 samples divided into 12 classes (`Dense_Woodland`, `Dunes`, `Fallow_Cotton`,  `Millet_Cotton`, `Pasture`, `Rocky_Savanna`, `Savanna`, `Savanna_Parkland`, `Silviculture`, `Soy_Corn`, `Soy_Cotton`, and  `Soy_Fallow`). Each time series covers 12 months (23 data points) from MOD13Q1 product, and has 4 bands (EVI, NDVI, MIR, and NIR). We use bands NDVI and EVI for faster processing.

:::{.panel-tabset}
## R
```{r}
# Take only the NDVI and EVI bands
samples_cerrado_mod13q1_2bands <- sits_select(
    data = samples_cerrado_mod13q1, 
    bands = c("NDVI", "EVI"))

# Show the summary of the samples
summary(samples_cerrado_mod13q1_2bands)
```
## Python
```{python}
# Load samples
samples_cerrado_mod13q1 = load_samples(
    name = "samples_cerrado_mod13q1", 
    package = "sitsdata"
)

# Take only the NDVI and EVI bands
samples_cerrado_mod13q1_2bands = sits_select(
    data = samples_cerrado_mod13q1, 
    bands = ("NDVI", "EVI"))

# Show the summary of the samples
summary(samples_cerrado_mod13q1_2bands)
```
:::


## Creating the SOM map

To perform the SOM-based quality assessment, the first step is to run `sits_som_map()`, which uses the `kohonen` R package to compute a SOM grid [@Wehrens2018], controlled by five parameters. The grid size is given by `grid_xdim` and `grid_ydim`. The starting learning rate is `alpha`, which decreases during the interactions. To measure the separation between samples, use `distance` (either "dtw" or "euclidean"). The number of iterations is set by `rlen`. When using `sits_som_map()` in machines which have multiprocessing support for the OpenMP protocol, setting the learning mode parameter `mode` to "patch" improves processing time. In Windows, please use "online". 

We suggest using the Dynamic Time Warping ("dtw") metric as the distance measure. It is a technique used to measure the similarity between two temporal sequences that may vary in speed or timing [@Berndt1994]. The core idea of DTW is to find the optimal alignment between two sequences by allowing non-linear mapping of one sequence onto another. In time series analysis, DTW matches two series slightly out of sync. This property is useful in land use studies for matching time series of agricultural areas [@Maus2015].

```{r}
#| echo: false
# Recover SOM cluster from RDS file
library(kohonen)
som_cluster <- readRDS(file = "./etc/som_cluster.rds")
```

```{python}
#| echo: false
# Recover SOM cluster from RDS file
som_cluster = read_rds(file = "./etc/som_cluster.rds")
```

:::{.panel-tabset}
## R
```{r}
#| eval: false
# Clustering time series using SOM
som_cluster <- sits_som_map(samples_cerrado_mod13q1_2bands,
    grid_xdim = 15,
    grid_ydim = 15,
    alpha = 1.0,
    distance = "dtw",
    rlen = 20
)

# Plot the SOM map
plot(som_cluster)
```
```{r}
#| echo: false
#| label: fig-ts-sommap
#| out.width: "100%"
#| fig.cap: "SOM map for the Cerrado samples."
#| fig.align: "center"
# Plot the SOM map
plot(som_cluster)
```
## Python
```{python}
#| eval: false
# Clustering time series using SOM
som_cluster = sits_som_map(samples_cerrado_mod13q1_2bands,
    grid_xdim = 15,
    grid_ydim = 15,
    alpha = 1.0,
    distance = "dtw",
    rlen = 20
)

# Plot the SOM map
plot(som_cluster)
```
```{python}
#| echo: false
#| label: fig-ts-sommap-py
#| out.width: "100%"
#| fig.cap: "SOM map for the Cerrado samples."
#| fig.align: "center"
# Clustering time series using SOM
# Plot the SOM map
plot(som_cluster)
```
:::

The output of the `sits_som_map()` is a list with three elements: (a) `data`, the original set of time series with two additional columns for each time series: `id_sample` (the original id of each sample) and `id_neuron` (the id of the neuron to which it belongs); (b) `labelled_neurons`, a tibble with information on the neurons. For each neuron, it gives the prior and posterior probabilities of all labels which occur in the samples assigned to it; and (c) the SOM grid. To plot the SOM grid, use `plot()`. The neurons are labelled using majority voting.

The SOM grid shows that most classes are associated with neurons close to each other, although there are exceptions. Some Pasture neurons are far from the main cluster because the transition between open savanna and pasture areas is not always well defined and depends on climate and latitude. Also, the neurons associated with Soy_Fallow are dispersed in the map, indicating possible problems in distinguishing this class from the other agricultural classes. The SOM map can be used to remove outliers, as shown below.

## Measuring confusion between labels using SOM

The second step in SOM-based quality assessment is understanding the confusion between labels. The function `sits_som_evaluate_cluster()` groups neurons by their majority label and produces a tibble. Neurons are grouped into clusters, and there will be as many clusters as there are labels. The results shows the percentage of samples of each label in each cluster. Ideally, all samples of each cluster would have the same label. In practice, cluster contain samples with different label. This information helps on measuring the confusion between samples. 

:::{.panel-tabset}
## R
```{r}
# Produce a tibble with a summary of the mixed labels
som_eval <- sits_som_evaluate_cluster(som_cluster)

# Show the result
som_eval 
```
## Python
```{python}
# Produce a tibble with a summary of the mixed labels
som_eval = sits_som_evaluate_cluster(som_cluster)

# Show the result
som_eval 
```
:::

Many labels are associated with clusters where there are some samples with a different label. Such confusion between labels arises because sample labeling is subjective and can be biased. In many cases, interpreters use high-resolution data to identify samples. However, the actual images to be classified are captured by satellites with lower resolution. In our case study, a MOD13Q1 image has pixels with 250 m resolution. As such, the correspondence between labeled locations in high-resolution images and mid to low-resolution images is not direct. The confusion by sample label can be visualized in a bar plot using `plot()`, as shown below. The bar plot shows some confusion between the labels associated with the natural vegetation typical of the Brazilian Cerrado (`Savanna`, `Savanna_Parkland`, `Rocky_Savanna`). This mixture is due to the large variability of the natural vegetation of the Cerrado biome, which makes it difficult to draw sharp boundaries between classes. Some confusion is also visible between the agricultural classes. The `Fallow_Cotton` class is a particularly difficult one since many of the samples assigned to this class are confused with `Soy_Cotton` and `Millet_Cotton`. 

:::{.panel-tabset}
## R
```{r}
#| label: fig-ts-someval
#| out.width: "90%"
#| message: false
#| warning: false
#| fig.align: "center"
#| fig.cap: "Confusion between classes as measured by SOM."
# Plot the confusion between clusters
plot(som_eval)
```
## Python
```{python}
#| label: fig-ts-someval-py
#| out.width: "90%"
#| message: false
#| warning: false
#| fig.align: "center"
#| fig.cap: "Confusion between classes as measured by SOM."
# Plot the confusion between clusters
plot(som_eval)
```
:::

## Detecting noisy samples using SOM

The third step in the quality assessment uses the discrete probability distribution associated with each neuron, which is included in the `labeled_neurons` tibble produced by `sits_som_map()`. This approach associates probabilities with frequency of occurrence. More homogeneous neurons (those with one label has high frequency) are assumed to be composed of good quality samples. Heterogeneous neurons (those with two or more classes with significant frequencies) are likely to contain noisy samples. The algorithm computes two values for each sample:

- *prior probability*: the probability that the label assigned to the sample is correct, considering the frequency of samples in the same neuron. For example, if a neuron has 20 samples, of which 15 are labeled as `Pasture` and 5 as `Forest`, all samples labeled Forest are assigned a prior probability of 25%. This indicates that Forest samples in this neuron may not be of good quality.

- *posterior probability*: the probability that the label assigned to the sample is correct, considering the neighboring neurons. Take the case of the above-mentioned neuron whose samples labeled `Pasture` have a prior probability of 75%. What happens if all the neighboring neurons have `Forest` as a majority label? To answer this question, we use Bayesian inference to estimate if these samples are noisy based on the surrounding neurons [@Santos2021a]. 

To identify noisy samples, we take the result of the `sits_som_map()` function as the first argument to the function `sits_som_clean_samples()`. This function finds out which samples are noisy, which are clean, and which need to be further examined by the user. It requires the `prior_threshold` and `posterior_threshold` parameters according to the following rules:

-   If the prior probability of a sample is less than `prior_threshold`, the sample is assumed to be noisy and tagged as "remove";
-   If the prior probability is greater or equal to `prior_threshold` and the posterior probability calculated by Bayesian inference is greater or equal to `posterior_threshold`, the sample is assumed not to be noisy and thus is tagged as "clean";
-   If the prior probability is greater or equal to `prior_threshold` and the posterior probability is less than `posterior_threshold`, we have a situation when the sample is part of the majority level of those assigned to its neuron, but its label is not consistent with most of its neighbors. This is an anomalous condition and is tagged as "analyze". Users are encouraged to inspect such samples to find out whether they are in fact noisy or not.

The default value for both `prior_threshold` and `posterior_threshold` is 60%. The `sits_som_clean_samples()` has an additional parameter (`keep`), which indicates which samples should be kept in the set based on their prior and posterior probabilities. The default for `keep` is `c("clean", "analyze")`. As a result of the cleaning, about 900 samples have been considered to be noisy and thus to be possibly removed. We first show the complete distribution of the samples and later remove the noisy ones.

:::{.panel-tabset}
## R
```{r}
#| label: fig-ts-somremove
#| out-width: 100%
#| fig-align: "center"
#| fig-cap: "Distribution of samples using som evaluation."
#| message: false
#| warning: false
all_samples <- sits_som_clean_samples(
    som_map = som_cluster, 
    prior_threshold = 0.6,
    posterior_threshold = 0.6,
    keep = c("clean", "analyze", "remove"))

# Print the sample distribution based on evaluation
plot(all_samples)
```
## Python
```{python}
#| label: fig-ts-somremove-py
#| out-width: 100%
#| fig-align: "center"
#| fig-cap: "Distribution of samples using som evaluation."
#| message: false
#| warning: false
all_samples = sits_som_clean_samples(
    som_map = som_cluster, 
    prior_threshold = 0.6,
    posterior_threshold = 0.6,
    keep = ("clean", "analyze", "remove"))

# Print the sample distribution based on evaluation
plot(all_samples)
```
:::

We now remove the noisy samples to improve the quality of the training set.

:::{.panel-tabset}
## R
```{r}
#| message: false
#| warning: false
new_samples <- sits_som_clean_samples(
    som_map = som_cluster, 
    prior_threshold = 0.6,
    posterior_threshold = 0.6,
    keep = c("clean", "analyze"))

# Print the new sample distribution
summary(new_samples)
```
## Python
```{python}
#| message: false
#| warning: false
new_samples = sits_som_clean_samples(
    som_map = som_cluster, 
    prior_threshold = 0.6,
    posterior_threshold = 0.6,
    keep = ("clean", "analyze"))

# Print the new sample distribution
summary(new_samples)
```
:::

All samples of the classes which had the highest confusion with others(`Fallow_Cotton`, `Silviculture`, and `Millet_Cotton`) are marken as noisy been removed. Classes `Fallow_Cotton` and `Millet_Cotton` are not distinguishable from other crops. Samples of class `Silviculture` (planted forests) have  removed since they have been confused with natural forests and woodlands in the SOM map. Further analysis includes calculating the SOM map and confusion matrix for the new set, as shown in the following example. 

:::{.panel-tabset}
## R
```{r}
#| eval: false
#| label: fig-ts-newcluster
#| out-width: 90%
#| fig-align: "center"
#| fig-cap: "Cluster confusion plot for samples cleaned by SOM."
#| message: false
#| warning: false
# Produce a new SOM map with the cleaned samples
new_cluster <- sits_som_map(
   data = new_samples,
   grid_xdim = 15,
   grid_ydim = 15,
   alpha = 1.0,
   rlen = 20,
   distance = "dtw")

# Evaluate the mixture in the new SOM clusters
new_cluster_mixture <- sits_som_evaluate_cluster(new_cluster)

# Plot the mixture information.
plot(new_cluster_mixture)
```
```{r}
#| message: false
#| warning: false
#| echo: false
# Evaluate the mixture in the SOM clusters of new samples
library(sits)
library(kohonen)
new_cluster <- readRDS("./etc/new_cluster.rds")
plot(
    sits_som_evaluate_cluster(new_cluster)
)
```

## Python
```{python}
#| eval: false
#| label: fig-ts-newcluster-py
#| out-width: 90%
#| fig-align: "center"
#| fig-cap: "Cluster confusion plot for samples cleaned by SOM."
#| message: false
#| warning: false
# Produce a new SOM map with the cleaned samples
new_cluster = sits_som_map(
   data = new_samples,
   grid_xdim = 15,
   grid_ydim = 15,
   alpha = 1.0,
   rlen = 20,
   distance = "dtw")
   
# Evaluate the mixture in the new SOM clusters
new_cluster_mixture = sits_som_evaluate_cluster(new_cluster)

# Plot the mixture information.
plot(new_cluster_mixture)
```
```{python}
#| message: false
#| warning: false
#| echo: false
new_cluster = read_rds("./etc/new_cluster.rds")
plot(
    sits_som_evaluate_cluster(new_cluster)
)
```
:::

As expected, the new confusion map shows a significant improvement over the previous one. This result should be interpreted carefully since it may be due to different effects. The most direct interpretation is that `Millet_Cotton` and `Silviculture` cannot be easily separated from the other classes, given the current attributes (a time series of NDVI and EVI indices from MODIS images). In such situations, users should consider improving the number of samples from the less represented classes, including more MODIS bands, or working with higher resolution satellites. The results of the SOM method should be interpreted based on the users' understanding of the ecosystems and agricultural practices of the study region. 

The SOM-based analysis discards samples that can be confused with samples of other classes. After removing noisy samples or uncertain classes, the dataset obtains a better validation score since there is less confusion between classes. Users should analyse the results with care. Not all discarded samples are low-quality ones. Confusion between samples of different classes can result from inconsistent labeling or from the lack of capacity of satellite data to distinguish between chosen classes. When many samples are discarded, as in the current example, revising the whole classification schema is advisable. The aim of selecting training data should always be to match the reality on the ground to the power of remote sensing data to identify differences. No analysis procedure can replace actual user experience and knowledge of the study region. 

## Summary 

In this chapter, we discuss the use of SOM as a proven clustering method for removing noisy samples and those that cannot be easily distinguishable from other samples of other classes. Experience with `sits` indicates that using SOM is a good way to assess data quality. In the next section, we focus on a complementary method of removing sample imbalance.

## References{-}