Skip to content

Commit

Permalink
More robust (pseudo)RGB plotting
Browse files Browse the repository at this point in the history
  • Loading branch information
mzarowka committed Mar 22, 2024
1 parent e95138d commit 776c473
Show file tree
Hide file tree
Showing 54 changed files with 673 additions and 1,302 deletions.
1 change: 1 addition & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ export(plot_raster_rgb)
export(prepare_core)
export(raster_crop)
export(remove_continuum)
export(roi_to_vect)
export(run_core)
export(spectra_position)
export(spectra_sub)
Expand Down
138 changes: 10 additions & 128 deletions R/normalization.R
Original file line number Diff line number Diff line change
@@ -1,70 +1,7 @@
#' Find position of selected spectra
#'
#' @param raster a terra SpatRaster.
#' @param spectra vector with choice of desired spectra.
#'
#' @return positions (indices) of desired spectra in SpatRaster
#' @export
#'
#' @description find index position of the nearest spectra (band) in the dataset.
#' Match for the lowest difference between integer band and actual SpatRaster band.
#' This will produce duplicates with multiple bands. Drop.
spectra_position <- function(raster, spectra) {
# Check if correct class is supplied.
if (!inherits(raster, what = "SpatRaster")) {
rlang::abort(message = "Supplied data is not a terra SpatRaster.")
}

# Find index (position) of selected spectra by comparing choice and names
spectraIndex <- purrr::map(spectra, \(x) which.min(abs(x - as.numeric(names(raster))))) |>
# Get positions
purrr::as_vector()

# Create tibble with spectra of choice and respective position
spectraIndex <- dplyr::tibble(
spectra = spectra,
position = spectraIndex) |>
# Keep second observation if duplicates are present
# From experience closer to desired product
dplyr::slice_tail(by = .data$position)

# Return values
return(spectraIndex)
}

#' Subset SpatRaster by spectra
#'
#' @param raster a terra SpatRaster to be subset.
#' @param spectra_tbl a tibble with spectra positions from spectra_position.
#'
#' @return SpatRaster subset to contain only required spectral bands.
#' @export
#'
#' @description subset SpatRaster with spectra (bands) positions.
spectra_sub <- function(raster, spectra_tbl) {
# Check if correct class is supplied.
if (!inherits(raster, what = "SpatRaster")) {
rlang::abort(message = "Supplied data is not a terra SpatRaster.")
}

# Get spectra from tibble
spectra <- dplyr::pull(spectra_tbl, 1)

# Get positions from tibble
position <- dplyr::pull(spectra_tbl, 2)

# Subset raster by position
raster <- terra::subset(raster, position)

# Set raster names to match spectra
names(raster) <- as.character(spectra)

# Return raster
return(raster)
}

#' Crop SpatRaster
#'
#' @family Normalization
#'
#' @param raster terra SpatRaster to be cropped.
#' @param type either data raster or reference raster.
#' @param dir directory.
Expand Down Expand Up @@ -179,71 +116,10 @@ raster_crop <- function(raster, type, dir = NULL, roi) {
return(raster)
}

#' #' Create reference SpatRaster
#' #'
#' #' @param raster terra SpatRaster of the captured reference.
#' #' @param roi Region Of Interest: extent to match data raster.
#' #' @param ref_type type of reference, one of "whiteref" or "darkref".
#' #' @param ... additional arguments.
#' #'
#' #' @return a terra SpatRaster of reference matching the data raster extent.
#' #' @export
#' #'
#' #' @description Creating reference SpatRaster covering core extent
#' #' Create one mean reference row SpatRaster by averaging data every column by aggregation
#' #' Create reference SpatRaster matching capture SpatRaster extent by disaggregation.
#' create_reference_raster_disagg <- function(raster, roi, ref_type, ...) {
#' # Store additional parameters
#' params <- rlang::list2(...)
#'
#' # Check if correct class is supplied.
#' if (!inherits(raster, what = "SpatRaster")) {
#' rlang::abort(message = "Supplied data is not a terra SpatRaster.")
#' }
#'
#' # Check number of rois and prepare ids.
#' if (length(roi) == 1) {
#' roi_id <- NULL
#' } else {
#' roi_id <- paste0("ROI_", seq(1:length(roi)))
#' }
#'
#' if (ref_type == "whiteref") {
#' name <- "WHITEREF"
#'
#' } else {
#' name <- "DARKREF"
#' }
#'
#' # Aggregate data into one row SpatRaster, divide by number of rows
#' cli::cli_alert("Aggregate { name }")
#'
#' raster <- terra::aggregate(
#' raster,
#' fact = c(terra::nrow(raster), 1),
#' fun = "mean",
#' overwrite = TRUE,
#' steps = terra::ncell(raster) * terra::nlyr(raster))
#'
#' # Set new extent to match extent of capture SpatRaster
#' terra::ext(raster) <- roi
#'
#' # Disaggregate data over entire extent to match capture SpatRaster extent, multiply by ymax
#' cli::cli_alert("Disaggregate { name }")
#'
#' raster <- terra::disagg(
#' raster,
#' fact = c(terra::ymax(raster), 1),
#' filename = paste0(params$path, "/products/", name, "_", basename(params$path), "_disaggregated.tif"),
#' overwrite = TRUE,
#' steps = terra::ncell(raster) * terra::nlyr(raster))
#'
#' # Return raster
#' return(raster)
#' }

#' Create reference SpatRaster
#'
#' @family Normalization
#'
#' @param raster terra SpatRaster of the captured reference.
#' @param roi Region Of Interest: extent to match data raster.
#' @param ref_type type of reference, one of "whiteref" or "darkref".
Expand Down Expand Up @@ -312,6 +188,8 @@ create_reference_raster <- function(raster, roi, ref_type, ...) {

#' Raster normalization: calculation
#'
#' @family Normalization
#'
#' @param capture a terra SpatRaster of captured data.
#' @param whiteref a terra SpatRaster of the white reference matching capture extent.
#' @param darkref a terra SpatRaster of the dark reference matching capture extent.
Expand Down Expand Up @@ -342,6 +220,8 @@ normalization <- function(capture = capture, whiteref = whiteref, darkref = dark

#' Raster normalization
#'
#' @family Normalization
#'
#' @param capture terra SpatRaster of captured data.
#' @param whiteref terra SpatRaster of the white reference matching capture extent.
#' @param darkref terra SpatRaster of the dark reference matching capture extent.
Expand All @@ -365,6 +245,8 @@ create_normalized_raster <- function(capture = capture, whiteref = whiteref, dar
# Create terra dataset
terra::sds()

# Do the while loop for file naming

# Apply function over the dataset and write to file
raster <- terra::lapp(
x = dataset,
Expand Down
83 changes: 69 additions & 14 deletions R/plotting.R
Original file line number Diff line number Diff line change
@@ -1,11 +1,16 @@
#' Stretch and optionally save full RGB preview of SpatRaster
#'
#' @family Plotting
#' @param raster a SpatRaster, preferably reflectance file.
#' @param ext character, a graphic format extension.
#' @param write logical, should resulting SpatRaster be written to file.
#'
#' @export
stretch_raster_full <- function(raster, ext = NULL, write = TRUE) {
stretch_raster_full <- function(
raster,
ext = NULL,
write = TRUE) {

# Check if correct class is supplied.
if (!inherits(raster, what = "SpatRaster")) {
rlang::abort(message = "Supplied data is not a terra SpatRaster.")
Expand All @@ -26,29 +31,39 @@ stretch_raster_full <- function(raster, ext = NULL, write = TRUE) {
fs::path_file() |>
fs::path_ext_remove()

cli::cli_h1("{raster_name}")

# Prepare name
filename <- paste0(raster_src, "/RGB_", raster_name, ".", ext)

# Check if there are values close to RGB, within the 25 nm.
if (all(purrr::list_c(purrr::map(c(450, 550, 650), \(i) dplyr::near(i, as.numeric(names(raster)), tol = 25)))) == TRUE) {
spectra <- c(450, 550, 650)
} else {
rlang::warn("No layers matching the RGB. Using the first, middle and last available layers.")
spectra <- c(
min(1:terra::nlyr(raster)),
terra::median(1:terra::nlyr(raster)),
max(terra::nlyr(raster))) |>
(\(i) as.numeric(names(1:terra::subset(raster, i))))()
}

# Check if raster is written to file
if (write == FALSE) {
# Subset and write to RGB
raster <- HSItools::spectra_position(
raster,
spectra = c(650, 550, 450)
spectra = spectra
) |>
HSItools::spectra_sub(
raster = raster,
spectra_tbl = _
) |>
terra::stretch()
} else {
cli::cli_alert("Writing RGB SpatRaster to {filename}")

# Subset and write to RGB
raster <- HSItools::spectra_position(
raster,
spectra = c(650, 550, 450)
spectra = spectra
) |>
HSItools::spectra_sub(
raster = raster,
Expand All @@ -66,8 +81,9 @@ stretch_raster_full <- function(raster, ext = NULL, write = TRUE) {

#' Plot spatial map plots of calculated proxies, and optionally save to file
#'
#' @family Plotting
#' @param raster a SpatRaster with calculated hyperspectral indices and RGB layers.
#' @param calibration result of pixel_to_distance or actuall call to pixel_to_distance with appropriate input.
#' @param calibration result of pixel_to_distance or actual call to pixel_to_distance with appropriate input.
#' @param hsi_index a character indicating hyperspectral index layer to plot.
#' @param palette a character indicating one of \pkg{viridis} palettes of choice: "viridis", "magma", "plasma", "inferno", "civids", "mako", "rocket" and "turbo”.
#' @param extent an extent or SpatVector used to subset SpatRaster. Defaults to the entire SpatRaster.
Expand All @@ -78,7 +94,15 @@ stretch_raster_full <- function(raster, ext = NULL, write = TRUE) {
#'
#' @return a plot with color map of selected hyperspectral index.
#' @export
plot_raster_proxy <- function(raster, hsi_index, calibration = NULL, palette = c("viridis", "magma", "plasma", "inferno", "civids", "mako", "rocket", "turbo"), extent = NULL, ext = NULL, write = FALSE, ...) {
plot_raster_proxy <- function(
raster,
hsi_index,
calibration = NULL,
palette = c("viridis", "magma", "plasma", "inferno", "civids", "mako", "rocket", "turbo"),
extent = NULL,
ext = NULL,
write = FALSE,
...) {
# Check if correct class is supplied.
if (!inherits(raster, what = "SpatRaster")) {
rlang::abort(message = "Supplied data is not a terra SpatRaster.")
Expand Down Expand Up @@ -141,7 +165,7 @@ plot_raster_proxy <- function(raster, hsi_index, calibration = NULL, palette = c
ggplot2::labs(
x = hsi_index,
y = "Depth (px)",
fill = "RABD"
fill = "Value"
)

} else {
Expand Down Expand Up @@ -176,7 +200,7 @@ plot_raster_proxy <- function(raster, hsi_index, calibration = NULL, palette = c
ggplot2::labs(
x = hsi_index,
y = "Depth (mm)",
fill = "RABD"
fill = "Value"
)
}

Expand All @@ -199,15 +223,21 @@ plot_raster_proxy <- function(raster, hsi_index, calibration = NULL, palette = c

#' Spatial map plots of RGB image
#'
#' @family Plotting
#' @param raster a SpatRaster with calculated hyperspectral indices and RGB layers or just RGB layers.
#' @param calibration result of pixel_to_distance or actuall call to pixel_to_distance with appropriate input.
#' @param calibration result of pixel_to_distance or actual call to pixel_to_distance with appropriate input.
#' @param extent an extent or SpatVector used to subset SpatRaster. Defaults to the entire SpatRaster.
#' @param ext character, a graphic format extension.
#' @param write logical, should resulting SpatRaster be written to file.
#'
#' @return a plot with color map of selected hyperspectral index.
#' @export
plot_raster_rgb <- function(raster, calibration = NULL, extent = NULL, ext = NULL, write = FALSE) {
plot_raster_rgb <- function(
raster,
calibration = NULL,
extent = NULL,
ext = NULL,
write = FALSE) {
# Check if correct class is supplied.
if (!inherits(raster, what = "SpatRaster")) {
rlang::abort(message = "Supplied data is not a terra SpatRaster.")
Expand All @@ -224,11 +254,25 @@ plot_raster_rgb <- function(raster, calibration = NULL, extent = NULL, ext = NUL
fs::path_file() |>
fs::path_ext_remove()

# Prepare filename
filename <- paste0(raster_src, "/RGB_GG_", raster_name, ".", ext)

# Check if there are values close to RGB, within the 25 nm.
if (all(purrr::list_c(purrr::map(c(450, 550, 650), \(i) dplyr::near(i, as.numeric(names(raster)), tol = 25)))) == TRUE) {
spectra <- c(450, 550, 650)
} else {
rlang::warn("No layers matching the RGB. Using the first, middle and last available layers.")
spectra <- c(
min(1:terra::nlyr(raster)),
terra::median(1:terra::nlyr(raster)),
max(1:terra::nlyr(raster))) |>
(\(i) as.numeric(names(terra::subset(raster, i))))()
}

# Prepare SpatRaster
raster <- HSItools::spectra_position(
raster,
spectra = c(650, 550, 450)) |>
spectra = spectra) |>
HSItools::spectra_sub(
raster = raster,
spectra_tbl = _)
Expand Down Expand Up @@ -317,6 +361,7 @@ plot_raster_rgb <- function(raster, calibration = NULL, extent = NULL, ext = NUL

#' Overlay color plot of proxy on RGB
#'
#' @family Plotting
#' @param raster raster a SpatRaster with calculated hyperspectral indices and RGB layers.
#' @param hsi_index a character indicating hyperspectral index layer to plot.
#' @param palette a character indicating one of \pkg{viridis} palettes of choice: "viridis”, “magma”, “plasma”, “inferno”, “civids”, “mako”, “rocket” and “turbo”.
Expand All @@ -327,7 +372,14 @@ plot_raster_rgb <- function(raster, calibration = NULL, extent = NULL, ext = NUL
#'
#' @return a plot with color map of selected hyperspectral index overlain on RGB image.
#' @export
plot_raster_overlay <- function(raster, hsi_index, palette = c("viridis”, “magma”, “plasma”, “inferno”, “civids”, “mako”, “rocket”, “turbo"), alpha = 0.5, extent = NULL, ext = NULL, write = FALSE) {
plot_raster_overlay <- function(
raster,
hsi_index,
palette = c("viridis”, “magma”, “plasma”, “inferno”, “civids”, “mako”, “rocket”, “turbo"),
alpha = 0.5,
extent = NULL,
ext = NULL,
write = FALSE) {
# Check if correct class is supplied.
if (!inherits(raster, what = "SpatRaster")) {
rlang::abort(message = "Supplied data is not a terra SpatRaster.")
Expand Down Expand Up @@ -432,6 +484,7 @@ plot_raster_overlay <- function(raster, hsi_index, palette = c("viridis”, “m
#' Composite hyperspectral indices plots
#' Can composite line profiles and SpatRasters
#'
#' @family Plotting
#' @param raster a SpatRaster with REFLECTANCE file. Used for correct placement.
#' @param plots a list of plots.
#' @param ext character, a graphic format extension.
Expand Down Expand Up @@ -489,6 +542,7 @@ plot_composite <- function(raster, plots, ext = NULL, write = FALSE) {

#' Line plots of calculated proxies series
#'
#' @family Plotting
#' @param raster a SpatRaster with calculated hyperspectral indices and RGB layers.
#' @param hsi_index a character indicating hyperspectral index layer to plot.
#' @param calibration result of pixel_to_distance or actual call to pixel_to_distance with appropriate input.
Expand Down Expand Up @@ -638,6 +692,7 @@ plot_profile_spectral_series <- function(raster, hsi_index, calibration = NULL,

#' Line plot of spectral profile from the ROI
#'
#' @family Plotting
#' @param raster Reflectance SpatRaster.
#' @param extent an extent or SpatVector used to subset SpatRaster. Defaults to the entire SpatRaster.
#' @param ext character, a graphic format extension.
Expand Down
Loading

0 comments on commit 776c473

Please sign in to comment.