Computing Spectral Indices

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:

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).

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.

list.files(idx2_out_dir, full.names = TRUE)
#> [1] "/tmp/Rtmpo7jOQM/prismaread/indices2/file2e77421948fc_GI.tif"   
#> [2] "/tmp/Rtmpo7jOQM/prismaread/indices2/file2e77421948fc_MSAVI.tif"

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, as well as from the list reported HERE.

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.

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:

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:

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)")
)