Computing-Spectral-Indexes.Rmd

---
title: "Computing Spectral Indices"
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(prismaread)
```

`prismaread` allows automatically computing spectral indices starting from 
either an original PRISMA hdf image, or from an hyperspectral cube already 
processed with `pr_convert()`.

## Computing spectral indices from a predefined list 

Spectral indices to be computed can be selected from a list of predefined ones. 
To do so, either specify a vector of desired indices as the `indexes` argument 
of `pr_convert()` to compute them directly from the hdf prisma data:

```{r index_1, echo=TRUE, message=FALSE, warning=FALSE}
l2d_he5_path <- file.path(
  system.file("testdata", package = "prismaread"),
  "PRS_L2D_STD_20200524103704_20200524103708_0001.he5"
)

# Download and unzip using piggyback if necessary
if (!file.exists(l2d_he5_path)) {
  message("Downloading test data - this may need a long time...")
  if (!requireNamespace("piggyback", quietly = TRUE)) {
    install.packages("piggyback")
  }
  l2d_zip_path <- file.path(
    system.file("testdata", package = "prismaread"),
    "PRS_L2D_STD_20200524103704_20200524103708_0001.zip"
  )
  piggyback::pb_download(
    "PRS_L2D_STD_20200524103704_20200524103708_0001.zip",
    repo = "irea-cnr-mi/prismaread",
    dest = dirname(l2d_zip_path)
  )
  piggyback::pb_track(glob = "inst/testdata/*.zip, inst/testdata/*.he5")
  unzip(l2d_zip_path, exdir = dirname(l2d_he5_path))
  unlink(l2d_zip_path)
}
idx_out_dir <- file.path(tempdir(), "prismaread/indices")
dir.create(dirname(idx_out_dir))
pr_convert(
  in_file = l2d_he5_path,
  out_format = "ENVI",
  out_folder = idx_out_dir, 
  indexes = c("GI", "MSAVI")
)
```

Alternatively, use function `pr_compute_indexes()` to compute them starting from 
an hyperspectral cube already processed (note that this is a bit slower and you
need to be sure that bands required to compute the index are available in the
file you are using as input).

```{r index_2, echo=TRUE, message=FALSE, warning=FALSE}
in_tif_path <- system.file(
  "testdata/prismaread_test_HCO_FULL.tif",
  package = "prismaread"
)
idx2_out_dir <- file.path(tempdir(), "prismaread/indices2")
dir.create(idx2_out_dir, recursive = TRUE)
idx2_out_path <- tempfile(fileext = ".tif", tmpdir = idx2_out_dir)
pr_compute_indexes(
  in_file = in_tif_path,
  out_file = idx2_out_path,
  indexes = c("GI", "MSAVI")
)
```

Output file names are created by adding a suffix corresponding to the indices
names to the base output filename. In addition, an ancillary file containing 
the formulas used to compute each index is saved (extension = `.formulas`) 
for reference. The file shows the formulas used, giving reference to the true
wavelengths of the PRISMA bands used in the computation.

```{r echo=TRUE, message=FALSE, warning=FALSE}
list.files(idx2_out_dir, full.names = TRUE)

idx2_out_MSAVI <- raster::raster(gsub(".tif", "_MSAVI.tif", idx2_out_path))
# Remove NA areas for better visualisation
idx2_out_MSAVI[idx2_out_MSAVI == -0.5 ] <- NA

mapview::mapview(idx2_out_MSAVI)
```


The list of available indices, shown below, was derived from the list of 
spectral indices available in package 
[`hsdar`](https://cran.r-project.org/web/packages/hsdar/hsdar.pdf),
as well as from the list reported
[HERE](https://cubert-gmbh.com/applications/vegetation-indices/).

```{r tbl, echo=FALSE, message=FALSE, warning=FALSE}
pr_listindexes()
```

**NOTE:** although a check was done on indices formulas, we do not guarantee
that all of them are correct. Please check the formulas in the table to be sure.

**IMPORTANT NOTE:** computation of spectral indices can be done alongside 
extraction of spectral cubes/ancillary info. For example, commands below
export the VNIR cube, LATLON and ANGLES datasets along with
images relative to "GI" and "MTCI" indices. 

```{r index2, eval=FALSE}
idx3_out_dir <- file.path(tempdir(), "prismaread/indices3")
pr_compute_indexes(
  in_file = l2d_he5_path,
  out_folder = idx3_out_dir,
  VNIR = TRUE, 
  LATLON = TRUE, 
  ANGLES = TRUE,
  out_format = "ENVI",
  indices = c("GI", "MTCI")
)
```


### Adding a new index to the list

Additional spectral indices can be added to the aforementioned list using 
function `pr_addindex()`, as in the following example:

```{r l1example4, eval=FALSE}
pr_addindex(
  Name = "myindex", 
  Formula = "R600 / R700", 
  Description = "My custom Index", 
  Reference = "Me (2020)"
)
```

## Computing custom spectral indices

Custom spectral indices can also be computed on the fly by specifying the
`cust_indexes` argument, as a __named list__, such as in:

```{r l1example2, eval=FALSE, message=FALSE, warning=FALSE}
idx4_out_dir <- file.path(tempdir(), "prismaread/indices4")
pr_convert(
  l2d_he5_path,
  out_format = "ENVI",
  out_folder = idx4_out_dir, 
  indices = c("GI", "MTCI"), 
  cust_indices = list(myindex1 = "R500 / R600",
                      myindex2 = "(R800 - R680) / (R800 + R680)")
)
```