---
title: "rOPTRAM: Deriving Soil Moisture from Sentinel-2 Imagery"
output: rmarkdown::html_document
vignette: >
  %\VignetteIndexEntry{rOPTRAM-1}
  %\VignetteEngine{knitr::rmarkdown_notangle}
  %\VignetteEncoding{UTF-8}
bibliography: "references.bib"
---



## Introduction  <img align="right" width="100" height="100" src="images/rOPTRAM_logo.jpg">

`rOPTRAM` implements The OPtical TRapezoid Model (OPTRAM) to derive soil moisture based on the linear relation between a vegetation index, i.e. NDVI, and Land Surface Temperature (LST). The Short Wave Infra-red (SWIR) band is used as a proxy for LST. The SWIR band is transformed to Swir Transformed Reflectance (STR).

A scatterplot of NDVI vs. STR is used to produce wet and dry linear regression lines, and the slope/intercept coefficients of these lines comprise the trapezoid. These coefficients are then used on a new satellite image to determine soil moisture.

  See:  @sadeghi_optical_2017, @burdun_satellite_2020, @ambrosone_retrieving_2020


### Prerequisites

Only a small number of commonly used R packages are required to use {rOPTRAM}. This includes:

 - base packages {tools} and {utils}
 - spatial packages {sf} and {terra}
 - data.frame and plotting {dplyr}, {ggplot2}

Also, to allow {rOPTRAM} to download Sentinel-2 images, clip to a study area, and prepare the necessary vegetation index and STR products, the R package {CDSE} (see @karaman_cdse_2023) as well as {jsonlite} are required.

## Workflows

Users can download Sentinel-2 tiles from the Copernicus manually, and run thru the steps one by one to produce the OPTRAM trapezoid, and predicted soil moisture maps. However, this approach is not optimal. The complete workflow can be initiated with a single function call to download, clip to area of interest, and produce the trapezoid coefficients. This all-inclusive approach is highly recommended since processing of the Sentinel-2 data is performed "in the cloud" and only the final products are downloaded, **greatly** reducing the download file sizes.

That recommended R package {CDSE} interfaces with the Copernicus DataSpace Ecosystem in one of two ways:

 - Thru the [Scihub API](https://shapps.dataspace.copernicus.eu/dashboard/#/).
 - Thru the [openEO platform](https://openeo.dataspace.copernicus.eu/)

Both methods require [registering](https://dataspace.copernicus.eu/) on the Copernicus DataSpace



``` r
remotes::install_github("ropensci/rOPTRAM")
library(rOPTRAM)
ls(getNamespace('rOPTRAM'))
# The {CDSE} and {jsonlite} packages are required for downloading Sentinel imagery from Copernicus DataSpace.
library("CDSE")
library("jsonlite")
```

## Package options

The {rOPTRAM} package is controlled by several options, set automatically when the package first loads. These defaults can be viewed by running the function `optram_options()` without arguments:


``` r
optram_options()
#> [1] "edge_points = TRUE"
#> [1] "feature_col = ID"
#> [1] "max_cloud = 12"
#> [1] "max_tbl_size = 1e+06"
#> [1] "period = seasonal"
#> [1] "plot_colors = colors"
#> [1] "remote = scihub"
#> [1] "rm.hi.str = FALSE"
#> [1] "rm.low.vi = FALSE"
#> [1] "SWIR_band = 11"
#> [1] "trapezoid_method = linear"
#> [1] "veg_index = NDVI"
#> [1] "vi_step = 0.005"
#> NULL
```

Users can override any of the defaults by calling this function with alternative values for the options. For example, Sentinel images are downloaded, by default for all available image dates between the specified `from_date` and `to_date`. Users can, instead, choose to download images only for a specific season as follows:


``` r
optram_options("period", "seasonal")
#> [1] "edge_points = TRUE"
#> [1] "feature_col = ID"
#> [1] "max_cloud = 12"
#> [1] "max_tbl_size = 1e+06"
#> [1] "period = seasonal"
#> [1] "plot_colors = colors"
#> [1] "remote = scihub"
#> [1] "rm.hi.str = FALSE"
#> [1] "rm.low.vi = FALSE"
#> [1] "SWIR_band = 11"
#> [1] "trapezoid_method = linear"
#> [1] "veg_index = NDVI"
#> [1] "vi_step = 0.005"
```
Now, only images between the day/month of `from_date` and day/month of `to_date` will be acquired, but for all years in the date range.


Three trapezoid fitting functions are implemented in {rOPTRAM}. The default is "linear", but users can change as follows:


``` r
optram_options("trapezoid_method", "polynomial")
#> [1] "edge_points = TRUE"
#> [1] "feature_col = ID"
#> [1] "max_cloud = 12"
#> [1] "max_tbl_size = 1e+06"
#> [1] "period = seasonal"
#> [1] "plot_colors = colors"
#> [1] "remote = scihub"
#> [1] "rm.hi.str = FALSE"
#> [1] "rm.low.vi = FALSE"
#> [1] "SWIR_band = 11"
#> [1] "trapezoid_method = polynomial"
#> [1] "veg_index = NDVI"
#> [1] "vi_step = 0.005"
```

The default trapezoid plot shows a scatter plot cloud of vegetation index values vs SWIR Transformed Reflectance, where all points have a uniform color. This plot can be improved by coloring the points by their density within the scatterplot. This is accomplished by choosing to set the option "plot_colors" to "colors":


``` r
optram_options("plot_colors", "colors")
#> [1] "edge_points = TRUE"
#> [1] "feature_col = ID"
#> [1] "max_cloud = 12"
#> [1] "max_tbl_size = 1e+06"
#> [1] "period = seasonal"
#> [1] "plot_colors = colors"
#> [1] "remote = scihub"
#> [1] "rm.hi.str = FALSE"
#> [1] "rm.low.vi = FALSE"
#> [1] "SWIR_band = 11"
#> [1] "trapezoid_method = polynomial"
#> [1] "veg_index = NDVI"
#> [1] "vi_step = 0.005"
```

## Main wrapper function

### Run the full OPTRAM model procedure with a single function call

Registration on Copernicus DataSpace was done in advance, and the OAuth credentials (*clientid*, and *secret*) have been saved to the user's home directory.

Downloaded Sentinel-2 images are saved to `S2_output_dir`. For this example, outputs are saved to `tempdir()`


``` r
from_date <- "2022-05-01"
to_date <- "2023-04-30"
output_dir <- tempdir()
aoi <- sf::st_read(system.file("extdata",
                               "lachish.gpkg", package = "rOPTRAM"))
optram_options("veg_index", "NDVI")
optram_options("trapezoid_method", "linear")
```


``` r
rmse <- optram(aoi,
              from_date, to_date,
              S2_output_dir = output_dir,
              data_output_dir = output_dir)
#> [1] "RMSE for fitted trapezoid:"
#>    RMSE.wet  RMSE.dry
#> 1 0.2747341 0.1407656
```

``` r
knitr::kable(rmse, caption = "Trapezoid coefficients")
```



Table: Trapezoid coefficients

|  RMSE.wet|  RMSE.dry|
|---------:|---------:|
| 0.2747341| 0.1407656|



### Show trapezoid plot


``` r
edges_df <- read.csv(file.path(output_dir, "trapezoid_edges_lin.csv"))
df_file <- file.path(output_dir, "VI_STR_data.rds")
full_df <- readRDS(df_file)
pl <- plot_vi_str_cloud(full_df, edges_df)
```

![plot of chunk plot-prepare](figure/plot-prepare-1.png)

``` r
pl <- pl + ggplot2::ggtitle("Lachish area trapezoid plot",
                            subtitle = "linear fitted")
ggplot2::ggsave(file.path(output_dir, "trapezoid_lachish_linear.png"),
                width = 18, height = 12, units = "cm")
```


``` r
knitr::include_graphics("images/trapezoid_lachish_linear.png")
```

![plot of chunk plot-trap](images/trapezoid_lachish_linear.png)

## Step by step

### The same procedure as the wrapper function, but in explicit steps

  - Acquire Sentinel 2 images within a date range, and crop to AOI;
  - Prepare the SWIR Transformed Reflectance;
  - Prepare a data.frame of Vegetation Index and STR values;
  - Activate the options to remove VI values below zero
  - and the option to remove outlier STR values
  - Get trapezoid coefficients from the scatterplot of VI-STR pixels


``` r
# Using the same dates, and aoi as previous example
s2_file_list <- optram_acquire_s2(aoi,
                            from_date, to_date,
                            output_dir = output_dir)
STR_list <- list.files(file.path(output_dir, "STR"),
                      pattern = ".tif$", full.names = TRUE)
VI_list <- list.files(file.path(output_dir, "NDVI"),
                      pattern = ".tif$", full.names = TRUE)
full_df <- optram_ndvi_str(STR_list, VI_list,
                           output_dir = output_dir)
rmse <- optram_wetdry_coefficients(full_df,
                                   output_dir = output_dir)
knitr::kable(rmse, caption = "RMSE of fitted trapezoid edges")
```



Table: RMSE of fitted trapezoid edges

|  RMSE.wet|  RMSE.dry|
|---------:|---------:|
| 0.2823248| 0.1338994|



## Soil Moisture Estimate

### Use trapezoid coefficients, VI, and STR rasters to derive soil moisture grid


``` r
img_date <- "2023-01-10"   # After a rain
VI_dir <- file.path(output_dir, "NDVI")
STR_dir <- file.path(output_dir, "STR")
data_dir <- output_dir
SM <- optram_calculate_soil_moisture(img_date, VI_dir, STR_dir,
                                     data_dir = data_dir)
```

### Soil moisture plot


``` r
library("viridis")
library("terra")
names(SM) <- "Lachish soil moisture"
SM[terra::values(SM) < 0] <- 0

terra::plot(SM, col = rev(viridis::viridis(25)), smooth = TRUE,
            main = "Soil moisture map")
```

![plot of chunk sm-plot](figure/sm-plot-1.png)


``` r
knitr::include_graphics("images/sm-plot.png")
```

![plot of chunk sm-plot-1](images/sm-plot.png)