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

prismaread provides a very efficient function based on package exactextractr(https://github.com/isciences/exactextractr) for extracting and summarising data from the converted hyperspectral cubes over features of vector spatial files. The function allows computing several statistics, as well as extracting all pixel values, and saving them to RData, CSV or EXCEL files (see prisma_extract_spectra() documentation for info on available arguments).

For example, starting from a VNIR Cube obtained with pr_convert() and a vector polygon file:

library(prismaread)
library(ggplot2)
library(tmap)
in_full_path  <- file.path(
  system.file("testdata", package = "prismaread"),
  "prismaread_test_HCO_FULL.tif"
)
in_poly_path <- system.file(
  "extdata/testdata/polys_prisma.gpkg",
  package = "prismaread"
)

# Download and unzip using piggyback if necessary
if (!file.exists(in_full_path)) {
  message("Downloading test data - this may need a long time...")
  if (!requireNamespace("piggyback", quietly = TRUE)) {
    install.packages("piggyback")
  }
  in_full_zip_path <- gsub("tif$", "zip", in_full_path)
  piggyback::pb_download(
    basename(in_full_zip_path),
    repo = "irea-cnr-mi/prismaread",
    dest = dirname(in_full_zip_path)
  )
  piggyback::pb_track(glob = "inst/testdata/*.tif")
  unzip(in_full_zip_path, exdir = dirname(in_full_path))
  unlink(in_full_zip_path)
}

# extract base statistics, in "long" format
extr1 <- pr_extract_spectra(in_full_path, in_poly_path, id_field = "field_id")
extr1
# plot results using ggplot
library(ggplot2)
ggplot(extr1, aes(x = wvl, y = mean)) +
  geom_line(aes(color = ID, group = ID)) +
  facet_wrap(~ID) +
  theme_light()

Argument stats_format defines the format used for statistics' output. If "long" (default, as above), the output has one column for the ID of the feature, and one column for each statistic. If "wide", the output has one column for each ID/statistic couple (e.g., mean_id_1, stdev_id_1, mean_id_2, etcetera).

in_vnir_path  <- file.path(
  system.file("testdata", package = "prismaread"),
  "prisma_test_HCO_VNIR.tif"
)

# Download and unzip using piggyback if necessary
if (!file.exists(in_full_path)) {
  message("Downloading test data - this may need a long time...")
  in_vnir_zip_path <- gsub("tif$", "zip", in_vnir_path)
  piggyback::pb_download(
    basename(in_vnir_zip_path),
    repo = "irea-cnr-mi/prismaread",
    dest = dirname(in_vnir_zip_path)
  )
  piggyback::pb_track(glob = "inst/testdata/*.tif")
  unzip(in_vnir_zip_path, exdir = dirname(in_vnir_path))
  unlink(in_vnir_zip_path)
}

# Extract base statistics for VNIR and save results as excel file, in "wide" format
extr2 <- pr_extract_spectra(
  in_vnir_path, in_poly_path,
  out_file = tempfile(fileext = ".xlsx"),
  stats_format = "wide", id_field = "field_id"
)
extr2

Additional arguments allow specifying which statistics are needed, as well as selecting if info regarding single pixels contained in each polygon need to be returned:

# extract custom statistics
extr3 <- pr_extract_spectra(
  in_vnir_path, in_poly_path,
  selstats = c("mean", "coeffvar", "stdev", "min", "max"),
  id_field = "field_id"
)
extr3

# plot results using ggplot
ggplot(extr3, aes(x = wvl)) +
  geom_ribbon(aes(ymin = mean - stdev, ymax = mean + stdev, fill = ID), alpha = 0.2) +
  geom_line(aes(y = mean, color = ID, colour = ID)) +
  facet_wrap(~ID) +
  theme_light()

# extract custom statistics and quantiles
extr4 <- pr_extract_spectra(
  in_vnir_path, in_poly_path,
  quantiles = TRUE,
  selstats = c("mean", "stdev"),
  id_field = "field_id"
)
extr4


# extract also all pixels
extr5 <- pr_extract_spectra(
  in_vnir_path, in_poly_path,
  allpix = TRUE,
  selstats = c("mean", "stdev"),
  id_field = "field_id"
)

# stats are saved in the "stats" slot of the output
extr5$stats

# pixel values are saved in the "allpix" slot of the output
extr5$allpix

ggplot(extr5$allpix, aes(x = wvl)) +
  geom_line(aes(y = value, group = pixel, color = ID), lwd = 0.02) +
  facet_wrap(~ID) +
  theme_light()

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