In lbusett/prismaread: Import PRISMA L1/L2 hyperspectral data and convert them to a more user friendly format

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:

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)

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.

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.

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

lbusett/prismaread documentation built on Feb. 22, 2022, 7:33 p.m.