ts_basics.qmd

---
title: "Basic operations on image time series"
format: html
---

<a href="https://www.kaggle.com/code/esensing/basic-operations-on-image-time-series" 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_basics"
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_basics"
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)
```
:::

## Data structure for image time series

In `sits`, time series are stored in a data frame data structure, both in R and Python. The following code shows the first three lines of a time series containing 1,882 labeled samples of land classes in Mato Grosso state of Brazil. The samples have time series extracted from the MODIS MOD13Q1 product from 2000 to 2016, provided every 16 days at 250 m resolution in the Sinusoidal projection. Based on ground surveys and high-resolution imagery, it includes samples of seven classes: `Forest`, `Cerrado`, `Pasture`, `Soy_Fallow`, `Soy_Cotton`, `Soy_Corn`, and `Soy_Millet`. 

:::{.panel-tabset}
## R
```{r}
# Samples
data("samples_matogrosso_mod13q1")
samples_matogrosso_mod13q1[1:4,]
```
## Python
```{python}
# Samples
samples_matogrosso_mod13q1 = load_samples(
    name = "samples_matogrosso_mod13q1", 
    package = "sitsdata"
)

samples_matogrosso_mod13q1[0:4]
```
:::

The time series data frame contains data and metadata. The first six columns contain spatial and temporal information, the label assigned to the sample, and the data cube from where the data has been extracted. The first sample has been labeled `Pasture` at location (-58.5631, -13.8844), being valid for the period (2006-09-14, 2007-08-29). Informing the dates where the label is valid is crucial for correct classification. In this case, the researchers labeling the samples used the agricultural calendar in Brazil. The relevant dates for other applications and other countries will likely differ from those used in the example. The `time_series` column contains the time series data for each spatiotemporal location. This data is also organized as a data frame, with a column with the dates and the other columns with the values for each spectral band. 

## Utilities for handling time series

The package provides functions for data manipulation and displaying information for time series data frames. For example, `summary()` shows the labels of the sample set and their frequencies.

```{python}
#| echo: false
samples_matogrosso_mod13q1 = load_samples(
    name = "samples_matogrosso_mod13q1", 
    package = "sitsdata"
)
```

:::{.panel-tabset}
## R
```{r}
summary(samples_matogrosso_mod13q1)
```
## Python
```{python}
summary(samples_matogrosso_mod13q1)
```
:::

In many cases, it is helpful to relabel the dataset. For example, there may be situations where using a smaller set of labels is desirable because samples in one label on the original set may not be distinguishable from samples with other labels. We then could use `sits_labels()<-` to assign new labels. The example below shows how to do relabeling on a time series set shown above; all samples associated with crops are grouped in a single `Cropland` label.

:::{.panel-tabset}
## R
```{r}
# Copy the sample set for Mato Grosso 
samples_new_labels <- samples_matogrosso_mod13q1

# Show the current labels
sits_labels(samples_new_labels)

# Update the labels
sits_labels(samples_new_labels) <- c("Cerrado",   "Forest",    
                                     "Pasture",   "Cropland", 
                                     "Cropland", "Cropland",
                                     "Cropland")

summary(samples_new_labels)
```
## Python
```{python}
# Copy the sample set for Mato Grosso 
samples_new_labels = samples_matogrosso_mod13q1.copy(deep = True)

# Get the current labels
sits_labels(samples_new_labels)

# Update the labels
samples_new_labels.sits.labels = (
    "Cerrado", "Forest", "Pasture", "Cropland", 
    "Cropland", "Cropland","Cropland"
)

summary(samples_new_labels)
```
:::

In R, the functions from `dplyr`, `tidyr`, and `purrr` packages of the `tidyverse` [@Wickham2017] can be used to process the data. In Python, time series data frames can be processed using `pandas`. For example, the following code uses `sits_select()` to get a subset of the sample dataset with two bands (NDVI and EVI). In R, it then uses the `dplyr::filter()` to select the samples labeled as `Cerrado`, while in Python the code uses the `query` function from the `pandas` package.

:::{.panel-tabset}
## R
```{r}
#| message: false
# Select NDVI band
samples_ndvi <- sits_select(samples_matogrosso_mod13q1, 
                            bands = "NDVI")

# Select only samples with Cerrado label
samples_cerrado <- dplyr::filter(samples_ndvi, 
                                 label == "Cerrado")
```
## Python
```{python}
#| message: false
# Select NDVI band
samples_ndvi = sits_select(samples_matogrosso_mod13q1, 
                            bands = "NDVI")

# Select only samples with Cerrado label
samples_cerrado = samples_ndvi.query("label == 'Cerrado'")
```
:::

## Time series visualisation

The default time series plot in `sits` aligns all samples onto a single reference time interval (based on the first sample) for display purposes, even when the samples' actual observations span different calendar years. This plot shows the spread of values for the time series of each band. The strong red line in the plot indicates the median of the values, while the two orange lines are the first and third interquartile ranges. See `?sits::plot` for more details on data visualization in  `sits`.

:::{.panel-tabset}
## R
```{r}
#| label: fig-ts-cerrado
#| fig.align: "center"
#| out.width: "70%"
#| fig.cap: "Plot of Cerrado samples for NDVI band."
# Plot cerrado samples together
plot(samples_cerrado)
```
## Python
```{python}
#| label: fig-ts-cerrado-py
#| fig.align: "center"
#| out.width: "70%"
#| fig.cap: "Plot of Cerrado samples for NDVI band."
# Plot cerrado samples together
plot(samples_cerrado)
```
:::

To see the spatial distribution of the samples, use `sits_view()` to create an interactive plot. The spatial visulisation is useful to show where the data has been collected. 

:::{.panel-tabset}
## R
```{r}
#| eval: false
sits_view(samples_matogrosso_mod13q1)
```
## Python
```{python}
#| eval: false
sits_view(samples_matogrosso_mod13q1)
```
:::

```{r}
#| label: fig-ts-view
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: "View location of training samples"
knitr::include_graphics("./images/view_samples.png")
```

## Patterns of training samples

When dealing with large time series, it is useful to obtain a single plot that captures the essential temporal variability of each class. Following the work on the `dtwSat` R package [@Maus2019], we use a generalized additive model (GAM) to obtain a single time series based on statistical approximation. In a GAM, the predictor depends linearly on a smooth function of the predictor variables. 

$$
y = \beta_{i} + f(x) + \epsilon, \epsilon \sim N(0, \sigma^2).
$$ 

The function `sits_patterns()` uses a GAM to predict an idealized approximation to the time series associated with each class for all bands. The resulting patterns can be viewed using `plot()`.

:::{.panel-tabset}
## R
```{r}
#| label: fig-ts-pat
#| out.width: "100%"
#| fig.align: "center"
#| fig.cap: "Patterns for Mato Grosso MOD13Q1 samples."
# Estimate the patterns for each class and plot them
samples_matogrosso_mod13q1 |>  
    sits_patterns() |> 
    plot()
```
## Python
```{python}
#| label: fig-ts-pat-py
#| out.width: "100%"
#| fig.align: "center"
#| fig.cap: "Patterns for Mato Grosso MOD13Q1 samples."
# Estimate the patterns for each class and plot them
plot(
    sits_patterns(
        samples_matogrosso_mod13q1
    )
)
```
:::

The resulting patterns provide some insights over the time series behaviour of each class. The response of the Forest class is quite distinctive.  They also show that it should be possible to separate between the single and double cropping classes. There are similarities between the double-cropping classes (`Soy_Corn` and `Soy_Millet `) and between the `Cerrado` and `Pasture` classes. The subtle differences between class signatures provide hints at possible ways by which machine learning algorithms might distinguish between classes. One example is the difference between the middle-infrared response during the dry season (May to September) to differentiate between `Cerrado` and `Pasture`. 

## Geographical variability of training data

When working with machine learning classification of Earth observation data, it is important to evaluate if the training samples are well distributed in the study area. Training data often comes from ground surveys made at chosen locations. In large areas, ideally representative samples need to capture spatial variability. In practice, however, ground surveys or other means of data collection are limited to selected areas. In many cases, the geographical distribution of the training data does not cover the study area equally. Such mismatch can be a problem for achieving a good quality classification. As stated by Meyer and Pebesma [@Meyer2022]: *"large gaps in geographic space do not always imply large gaps in feature space"*. 

Meyer and Pebesma propose using a spatial distance distribution plot, which displays two distributions of nearest-neighbor distances: sample-to-sample and prediction-location-to-sample [@Meyer2022]. The difference between the two distributions reflects the degree of spatial clustering in the reference data. Ideally, the two distributions should be similar. Cases where the sample-to-sample distance distribution does not match prediction-location-to-sample distribution indicate possible problems in training data collection. 

`sits` implements spatial distance distribution plots with the `sits_geo_dist()` function. This function gets a training data in the `samples` parameter, and the study area in the `roi` parameter expressed as an `sf` object. Additional parameters are `n` (maximum number of samples for each distribution) and `crs` (coordinate reference system for the samples). By default, `n` is 1000, and `crs` is "EPSG:4326". The example below shows how to use `sits_geo_dist()`. 

:::{.panel-tabset}
## R
```{r}
#| echo: false
set.seed(777)
```
```{r}
#| label: fig-ts-geodist
#| out.width: "100%"
#| message: false
#| warning: false
#| fig.cap: | 
#|    Distribution of sample-to-sample and sample-to-prediction distances.
# Read a shapefile for the state of Mato Grosso, Brazil
mt_shp <- system.file("extdata/shapefiles/mato_grosso/mt.shp",
                      package = "sits")

# Convert to an sf object
mt_sf <- sf::read_sf(mt_shp)

# Calculate sample-to-sample and sample-to-prediction distances
distances <- sits_geo_dist(
    samples = samples_modis_ndvi,
    roi = mt_sf)

# Plot sample-to-sample and sample-to-prediction distances
plot(distances)
```
## Python
```{python}
#| echo: false
r_set_seed(777)
```
```{python}
#| label: fig-ts-geodist-py
#| out.width: "100%"
#| message: false
#| warning: false
#| fig.cap: | 
#|    Distribution of sample-to-sample and sample-to-prediction distances.
# Import GeoPandas
import geopandas as gpd

# Read a shapefile for the state of Mato Grosso, Brazil
mt_shp = r_package_dir("extdata/shapefiles/mato_grosso/mt.shp",
                      package = "sits")

# Convert to an GeoDataFrame object
mt_sf = gpd.read_file(mt_shp)

# Calculate sample-to-sample and sample-to-prediction distances
distances = sits_geo_dist(
    samples = samples_modis_ndvi,
    roi = mt_sf)

# Plot sample-to-sample and sample-to-prediction distances
plot(distances)
```
:::

The plot shows a mismatch between the sample-to-sample and the sample-to-prediction distributions. Most samples are closer to each other than they are close to the location where values need to be predicted. In this case, there are many areas where few or no samples have been collected and where the prediction uncertainty will be higher. In this and similar cases, improving the distribution of training samples is always welcome. If that is not possible, areas with insufficient samples could have lower accuracy. This information must be reported to potential users of classification results. 

## Time series from data cubes

To get a set of time series in `sits`, first create a regular data cube and then request one or more time series from the cube using `sits_get_data()`.  This function uses two mandatory parameters: `cube` and `samples`. The `cube` indicates the data cube from which the time series will be extracted.  The `samples` parameter accepts the following data types: 

- A data.frame with information on `latitude` and `longitude` (mandatory), `start_date`, `end_date`, and `label` for each sample point.
- A csv file with columns `latitude`, `longitude`, `start_date`, `end_date`, and `label`.
- A shapefile containing either `POINT`or `POLYGON` geometries. See details below.
- An `sf` object (from the `sf` package in R) with `POINT` or `POLYGON` geometry.
- A `geopandas` object from the `geopandas` package in Python with `POINT` or `POLYGON` geometry.

In the example below, given a data cube, the user provides the latitude and longitude of the desired location. Since the bands, start date, and end date of the time series are missing, `sits` obtains them from the data cube. The result is a tibble with one time series that can be visualized using `plot()`.

:::{.panel-tabset}
## R
```{r} 
#| label: fig-ts-get
#| fig.align: "center"
#| out.width: "90%" 
#| fig.cap: "NDVI and EVI time series fetched from local raster cube."
# Obtain a raster cube based on local files
data_dir <- system.file("extdata/sinop", package = "sitsdata")
raster_cube <- sits_cube(
    source     = "BDC",
    collection = "MOD13Q1-6.1",
    data_dir   = data_dir,
    parse_info = c("satellite", "sensor", "tile", "band", "date"))

# Obtain a time series from the raster cube from a point
sample_latlong <- tibble::tibble(
    longitude = -55.57320, 
    latitude  = -11.50566)

series <- sits_get_data(cube    = raster_cube,
                        samples = sample_latlong)

plot(series)
```
## Python
```{python}
#| label: fig-ts-get-py
#| fig.align: "center"
#| out.width: "90%" 
#| fig.cap: "NDVI and EVI time series fetched from local raster cube."
# Import pandas
import pandas as pd

# Obtain a raster cube based on local files
data_dir = r_package_dir("extdata/sinop", package = "sitsdata")
raster_cube = sits_cube(
    source     = "BDC",
    collection = "MOD13Q1-6.1",
    data_dir   = data_dir,
    parse_info = ("satellite", "sensor", "tile", "band", "date"))

# Obtain a time series from the raster cube from a point
sample_latlong = pd.DataFrame([dict(
    longitude = -55.57320,
    latitude  = -11.50566
)])

series = sits_get_data(cube = raster_cube,
                        samples = sample_latlong)

plot(series)
```
:::

A useful case is when a set of labeled samples can be used as a training dataset. In this case, trusted observations are usually labeled and commonly stored in plain text files in comma-separated values (csv) or using shapefiles (shp). 

:::{.panel-tabset}
## R
```{r}
#| message: false
# Retrieve a list of samples described by a csv file
samples_csv_file <- system.file("extdata/samples/samples_sinop_crop.csv",
                                package = "sits")

# Read the csv file into an R object
samples_csv <- read.csv(samples_csv_file)

# Print the first three samples
samples_csv[1:3,]
```
## Python
```{python}
#| message: false
# Retrieve a list of samples described by a csv file
samples_csv_file = r_package_dir("extdata/samples/samples_sinop_crop.csv",
                                package = "sits")

# Read the csv file into an R object
samples_csv = pd.read_csv(samples_csv_file)

# Print the first three samples
samples_csv[0:3]
```
:::

To retrieve training samples for time series analysis, users must provide the temporal information (`start_date` and `end_date`). In the simplest case, all samples share the same dates. That is not a strict requirement. It is possible to specify different dates as long as they have a compatible duration. For example, the dataset `samples_matogrosso_mod13q1` provided with the `sitsdata` package contains samples from different years covering the same duration. These samples are from the MOD13Q1 product, which contains the same number of images per year. Thus, all time series in the dataset `samples_matogrosso_mod13q1` have the same number of dates. 

Given a suitably built csv sample file, `sits_get_data()` requires two parameters: (a) `cube`, the name of the R object that describes the data cube; (b) `samples`, the name of the CSV file.

:::{.panel-tabset}
## R
```{r}
# Get the points from a data cube in raster brick format
points <- sits_get_data(cube = raster_cube, 
                        samples = samples_csv_file)

# Show the tibble with the first three points
points[1:3,]
```
## Python
```{python}
# Get the points from a data cube in raster brick format
points = sits_get_data(cube = raster_cube, 
                        samples = samples_csv_file)

# Show the tibble with the first three points
points[1:3]
```
:::

Users can also specify samples by providing shapefiles, `sf` objects (in R), or `geopandas` objects (in Python) containing `POINT` or `POLYGON` geometries. The geographical location is inferred from the geometries associated with the shapefile or `sf` (or `geopandas`) object. For files containing points, the geographical location is obtained directly. For polygon geometries, the parameter `n_sam_pol` (defaults to 20) determines the number of samples to be extracted from each polygon. The temporal information can be provided explicitly by the user; if absent, it is inferred from the data cube. If label information is available in the shapefile or `sf` object, the parameter `label_attr` is indicates which column contains the label associated with each time series.

In what follows, we provide a shapefile with location of forested areas in the Cerrado biome in Brazil. The shapefile is first used to define the region of interest to retrieve a data cube. Then the point locations are used to retrieve the time series.

:::{.panel-tabset}
## R
```{r}
#| eval: false
# Obtain a set of points inside the state of Mato Grosso, Brazil
shp_file <- system.file("extdata/shapefiles/cerrado/cerrado_forested.shp", 
                        package = "sits")
# Read the shapefile into an "sf" object
sf_shape <- sf::st_read(shp_file)

# Create a data cube based on MOD13Q1 collection from BDC
modis_cube <- sits_cube(
    source      = "BDC",
    collection  = "MOD13Q1-6.1",
    bands       = c("NDVI", "EVI"),
    roi         = sf_shape,
    start_date  = "2020-06-01", 
    end_date    = "2021-08-29")

# Read the points from the cube and produce a tibble with time series
samples_cerrado_forested <- sits_get_data(
    cube         = modis_cube, 
    samples      = shp_file, 
    start_date   = "2020-06-01",
    end_date     = "2021-08-29", 
    label        = "Woody Savanna",
    multicores   = 4)

# Display the time series for the locations of Woody Savanna
samples_cerrado_forested
```
## Python
```{python}
#| eval: false
# Obtain a set of points inside the state of Mato Grosso, Brazil
shp_file = r_package_dir("extdata/shapefiles/cerrado/cerrado_forested.shp", 
                        package = "sits")
# Read the shapefile into an "geopandas" object
sf_shape = gpd.read_file(shp_file)

# Create a data cube based on MOD13Q1 collection from BDC
modis_cube = sits_cube(
    source      = "BDC",
    collection  = "MOD13Q1-6.1",
    bands       = ("NDVI", "EVI"),
    roi         = sf_shape,
    start_date  = "2020-06-01", 
    end_date    = "2021-08-29")

# Read the points from the cube and produce a tibble with time series
samples_cerrado_forested = sits_get_data(
    cube         = modis_cube, 
    samples      = shp_file, 
    start_date   = "2020-06-01",
    end_date     = "2021-08-29", 
    label        = "Woody Savanna",
    multicores   = 4)

# Display the time series for the locations of Woody Savanna
samples_cerrado_forested
```
:::

```{r}
#| echo: false
samples_cerrado_forested <- readRDS("./etc/samples_cerrado_forested.rds")

# Display the time series for the locations of Woody Savanna
samples_cerrado_forested
```

```{python}
#| echo: false
# read it to have it in memory
samples_cerrado_forested = read_rds("./etc/samples_cerrado_forested.rds")
```

## Summary

In this chapter, we presented the data structure used  by `sits` to store pixel-based image time series. The text also shows how to retrieve time series from data cubes, as well as utilities available for visualisation and extracting patterns and geographical distribution. In the next chapters, we discuss how to improve the quality of training samples. 

## References{-}