diff --git a/.Rbuildignore b/.Rbuildignore index 3c1cbd1..c581ec9 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -7,3 +7,5 @@ ^README\.Rmd$ ^data-raw$ ^LICENSE\.md$ +^doc$ +^Meta$ diff --git a/.gitignore b/.gitignore index 758507c..61c0935 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,6 @@ .Ruserdata docs rando +inst/doc +/doc/ +/Meta/ diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..2591a4a --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,74 @@ +# ngr - New Graph Reporting + +R package with utility functions for dynamic reporting, spatial analysis, hydrology, and data wrangling. + +## Build & Test Commands + +```r +devtools::document() +devtools::test() +devtools::check() +devtools::install() +``` +Build documentation and run checks before committing. + +## Commit Style (fledge) + +This repo uses [fledge](https://github.com/cynkra/fledge) for changelog management. + +**NEWS-worthy commits** start with `- ` (bullet point): +``` +- initial commit of [ngr_function_name()] to close [#1](https://github.com/NewGraphEnvironment/ngr/issues/1) +``` + +**Function references** use square brackets for pkgdown auto-links: +- `[function_name()]` - creates link to function docs in pkgdown site + +**Issue references** use full markdown links: +- `[#19](https://github.com/NewGraphEnvironment/ngr/issues/19)` + +**Non-NEWS commits** (docs, chores) don't start with `- `: +``` +rebuild documentation and update build config +``` + +## Function Naming Conventions + +All exported functions use prefix `ngr_` followed by category: + +| Prefix | Category | Example | +|--------|----------|---------| +| `ngr_str_` | String manipulation | `ngr_str_extract_between()` | +| `ngr_spk_` | Spatial/raster (spacehakr) | `ngr_spk_stac_calc()` | +| `ngr_hyd_` | Hydrology | `ngr_hyd_q_daily()` | +| `ngr_dbqs_` | Database/SQL queries | `ngr_dbqs_filter_predicate()` | +| `ngr_tidy_` | Data frame tidying | `ngr_tidy_cols_rm_na()` | +| `ngr_xl_` | Excel operations | `ngr_xl_read_formulas()` | +| `ngr_fs_` | File system | `ngr_fs_copy_if_missing()` | +| `ngr_s3_` | S3/cloud storage | `ngr_s3_dl()` | +| `ngr_git_` | Git/GitHub | `ngr_git_issue()` | +| `ngr_chk_` | Validation/checking | `ngr_chk_dt_complete()` | +| `ngr_pkg_` | Package utilities | `ngr_pkg_detach()` | + +## Documentation Style + +- Use roxygen2 with markdown enabled +- Add `@family` tags to group related functions (e.g., `@family spacehakr`, `@family string`) +- Use `@importFrom` for all external functions +- Include `@examples` (use `\dontrun{}` for slow/network examples) +- Reference related functions with `@seealso` + +## Key Dependencies + +- **chk** - Argument validation +- **cli** - User messaging +- **sf/terra** - Spatial operations +- **tidyhydat** - Hydrometric data +- **dplyr/purrr/stringr** - Data wrangling +- **glue** - String interpolation + +## DESCRIPTION Management + +- Keep Imports alphabetized +- Don't duplicate packages in both Imports and Suggests +- Add vignette-only packages to Suggests (e.g., mapview, rstac) diff --git a/DESCRIPTION b/DESCRIPTION index c98b8c0..4912b9a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -9,37 +9,41 @@ License: MIT + file LICENSE Encoding: UTF-8 Roxygen: list(markdown = TRUE) RoxygenNote: 7.3.2.9000 -Suggests: +Suggests: testthat (>= 3.0.0), mockery, - sf, - purrr + knitr, + rmarkdown, + mapview, + rstac, + stars Config/testthat/edition: 3 URL: https://github.com/NewGraphEnvironment/ngr, https://newgraphenvironment.github.io/ngr/ BugReports: https://github.com/NewGraphEnvironment/ngr/issues -Imports: +Imports: chk, - processx, cli, - stringr, - fs, curl, - rvest, - lifecycle (>= 1.0.0), - janitor, - tidyxl, - sf, - terra, - purrr, - httr, dplyr, + fs, gh, + glue, + httr, + janitor, + lifecycle (>= 1.0.0), + processx, + purrr, + rvest, + sf, + stringr, + terra, tibble, tidyhydat, - xml2, - gh + tidyxl, + xml2 Remotes: nacnudus/tidyxl Depends: R (>= 2.10) LazyData: true +VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 5b20c9e..ceac108 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,7 @@ # Generated by roxygen2: do not edit by hand +export(ngr_chk_coerce_date) +export(ngr_chk_dt_complete) export(ngr_dbqs_filter_predicate) export(ngr_dbqs_ltree) export(ngr_dbqs_tbl_quote) @@ -7,6 +9,7 @@ export(ngr_fs_copy_if_missing) export(ngr_fs_id_missing) export(ngr_git_issue) export(ngr_git_issue_details) +export(ngr_hyd_q_daily) export(ngr_hyd_realtime) export(ngr_pkg_detach) export(ngr_s3_dl) @@ -25,10 +28,13 @@ export(ngr_spk_rast_ext) export(ngr_spk_rast_not_empty) export(ngr_spk_rast_rm_empty) export(ngr_spk_res) +export(ngr_spk_stac_calc) export(ngr_str_df_col_agg) export(ngr_str_df_detect_filter) +export(ngr_str_df_extract) export(ngr_str_dir_from_file) export(ngr_str_dir_from_path) +export(ngr_str_extract_between) export(ngr_str_link_url) export(ngr_str_replace_filenames) export(ngr_str_replace_in_files) @@ -38,6 +44,7 @@ export(ngr_tidy_type) export(ngr_xl_map_colnames) export(ngr_xl_map_formulas) export(ngr_xl_read_formulas) +importFrom(chk,abort_chk) importFrom(chk,chk_character) importFrom(chk,chk_data) importFrom(chk,chk_dir) @@ -73,6 +80,7 @@ importFrom(dplyr,across) importFrom(dplyr,all_of) importFrom(dplyr,arrange) importFrom(dplyr,bind_rows) +importFrom(dplyr,distinct) importFrom(dplyr,filter) importFrom(dplyr,group_by) importFrom(dplyr,if_else) @@ -121,21 +129,28 @@ importFrom(sf,st_read) importFrom(sf,st_sample) importFrom(sf,st_sf) importFrom(sf,st_transform) +importFrom(stringr,regex) importFrom(stringr,str_detect) importFrom(stringr,str_extract) +importFrom(stringr,str_match) importFrom(stringr,str_replace) importFrom(stringr,str_replace_all) +importFrom(stringr,str_squish) importFrom(stringr,str_which) importFrom(terra,crop) importFrom(terra,crs) importFrom(terra,ext) +importFrom(terra,mask) importFrom(terra,project) importFrom(terra,rast) importFrom(terra,res) importFrom(terra,union) importFrom(terra,values) +importFrom(terra,vect) importFrom(tibble,as_tibble) importFrom(tibble,tibble) +importFrom(tidyhydat,hy_daily_flows) +importFrom(tidyhydat,realtime_stations) importFrom(tidyhydat,realtime_ws) importFrom(tidyxl,xlsx_cells) importFrom(utils,glob2rx) diff --git a/R/ngr_chk_coerce_date.R b/R/ngr_chk_coerce_date.R new file mode 100644 index 0000000..782abf6 --- /dev/null +++ b/R/ngr_chk_coerce_date.R @@ -0,0 +1,43 @@ +#' Check and Validate Date-Coercible Input +#' +#' Checks whether an object is a [Date] or can be safely coerced to a [Date] +#' using [base::as.Date()]. If the input is coercible, it is returned +#' invisibly; otherwise, an informative error is thrown. +#' +#' @param x [any] An object expected to be a [Date] or a string coercible to a +#' [Date]. +#' @param x_name [character] Optional. A single string giving the name of `x` +#' used in error messages. If `NULL`, the name is inferred from the calling +#' expression. +#' +#' @details +#' This helper is intended for lightweight validation in functions that accept +#' date inputs but do not need to modify them. Coercion is attempted via +#' [base::as.Date()], and failure results in a call to +#' [chk::abort_chk()]. +#' +#' @returns +#' `x` invisibly if it is a [Date] or coercible to one. An error is thrown if +#' coercion fails. +#' +#' @seealso +#' [base::as.Date()], [chk::abort_chk()] +#' +#' @family chk +#' +#' @importFrom chk abort_chk +#' +#' @export +ngr_chk_coerce_date <- function(x, x_name = NULL) { + if (!inherits(try(as.Date(x), silent = TRUE), "try-error")) { + return(invisible(x)) + } + if (is.null(x_name)) { + x_name <- deparse(substitute(x)) + } + chk::abort_chk( + x_name, " must be a Date or a string coercible to Date", + x = x + ) +} + diff --git a/R/ngr_chk_dt_complete.R b/R/ngr_chk_dt_complete.R new file mode 100644 index 0000000..0f6e034 --- /dev/null +++ b/R/ngr_chk_dt_complete.R @@ -0,0 +1,66 @@ +#' Check for Complete Date Sequences +#' +#' Checks whether a vector of dates forms a complete, uninterrupted sequence +#' between its minimum and maximum values. If missing dates are detected, +#' they are printed using `dput()` for easy reuse and the function returns +#' `FALSE` invisibly. +#' +#' @param x [Date] A vector of dates to check for completeness. +#' @param units [character] Optional. A single string giving the time unit +#' for the sequence passed to [base::seq.Date()]. Default is `"days"`. +#' @param dates_print [logical] Optional. Whether to print missing dates to the +#' console using `dput()` when gaps are detected. Default is `TRUE`. +#' @param dates_capture [logical] Optional. Whether to capture and return the +#' missing dates as an object when gaps are detected. Default is `FALSE`. +#' +#' @details +#' The function computes a full date sequence using [base::seq.Date()] from +#' `min(x)` to `max(x)` and compares it to the unique values of `x` using +#' [base::setdiff()]. Missing dates, if any, are printed with `dput()` so they +#' can be copied directly into code or tests. +#' +#' @return +#' If `dates_capture = FALSE` (default), returns `TRUE` invisibly if no dates +#' are missing and `FALSE` invisibly if one or more dates are missing. +#' +#' If `dates_capture = TRUE`, returns the missing dates as a [Date] vector +#' (invisibly) when gaps are detected, and `TRUE` invisibly when no dates +#' are missing. +#' +#' @examples +#' dates_ok <- as.Date("2024-01-01") + 0:4 +#' ngr_chk_dt_complete(dates_ok) +#' +#' dates_bad <- as.Date(c("2024-01-01", "2024-01-02", "2024-01-04")) +#' ngr_chk_dt_complete(dates_bad) +#' +#' @importFrom cli cli_alert_warning cli_alert_success +#' +#' @export +#' @family chk +ngr_chk_dt_complete <- function(x, units = "days", dates_print = TRUE, dates_capture = FALSE) { + ngr_chk_coerce_date(x) + + x <- as.Date(x) + + full_seq <- seq.Date( + from = min(x, na.rm = TRUE), + to = max(x, na.rm = TRUE), + by = units + ) + + dates_missing <- as.Date(setdiff(full_seq, unique(x)), origin = "1970-01-01") + + if (length(dates_missing) > 0) { + cli::cli_alert_warning("There are missing dates:") + if (dates_print) { + dput(as.character(dates_missing)) + } + if (dates_capture) { + return(invisible(dates_missing)) + } + return(invisible(FALSE)) + } + + invisible(TRUE) +} diff --git a/R/ngr_hyd_q_daily.R b/R/ngr_hyd_q_daily.R new file mode 100644 index 0000000..a376bb7 --- /dev/null +++ b/R/ngr_hyd_q_daily.R @@ -0,0 +1,121 @@ +#' Combine daily Streamflow Data from HYDAT and ECCC Realtime Sources +#' +#' Retrieve and amalgamate daily streamflow data for one or more HYDAT stations, +#' combining historical daily flows from HYDAT with available realtime daily flow +#' data where supported. The function warns when realtime data are unavailable and +#' checks for missing dates in the combined time series. +#' +#' @param id_station [character] A character vector of one or more HYDAT station IDs +#' @param date_rt_start [Date] Optional. A single start date (or string coercible to date) for realtime data retrieval. +#' Defaults to `Sys.Date() - 581` (18 months). +#' @param date_rt_end [Date] Optional. A single end date (or string coercible to date) for realtime data retrieval. +#' Defaults to `Sys.Date()`. +#' @param date_hydat_start [character] Optional. A single date (or string coercible to date) giving the start +#' date for historical HYDAT daily flows. Default is "1980-01-01". +#' @param ... Optional. Additional arguments passed to +#' [ngr::ngr_chk_dt_complete()] when checking for missing dates within each +#' `STATION_NUMBER` time series (ex print_missing = FALSE). +#' +#' @details +#' Historical daily flows are retrieved using +#' [tidyhydat::hy_daily_flows()]. Realtime daily flows are retrieved using +#' [tidyhydat::realtime_ws()] for stations that support realtime reporting, as +#' identified by [tidyhydat::realtime_stations()]. Realtime data are coerced to +#' daily resolution and typed to match the HYDAT data using +#' [ngr::ngr_tidy_type()]. +#' +#' Input dates are validated using [ngr::ngr_chk_coerce_date()]. Missing daily +#' timestamps are checked per station using [ngr::ngr_chk_dt_complete()]. +#' +#' Duplicate station-date records are removed, keeping historical HYDAT values +#' where overlaps occur. The function reports stations with missing dates using +#' cli messages. +#' +#' @returns +#' A data frame containing daily streamflow data for the requested stations. +#' Warnings are emitted when realtime data are unavailable for some stations or +#' when missing dates are detected. +#' +#' @examples +#' ngr_hyd_q_daily(id_station = c("08NH118", "08NH126")) +#' +#' @importFrom tidyhydat hy_daily_flows realtime_ws realtime_stations +#' @importFrom dplyr bind_rows distinct +#' @importFrom cli cli_alert_warning cli_alert_success +#' +#' @export +ngr_hyd_q_daily <- function( + id_station, + date_rt_start = Sys.Date() - 581, + date_rt_end = Sys.Date(), + date_hydat_start = "1980-01-01", + ... +){ + ngr_chk_coerce_date(date_rt_start) + ngr_chk_coerce_date(date_rt_end) + ngr_chk_coerce_date(date_hydat_start) + chk::chk_character(id_station) + + hy_hy <- tidyhydat::hy_daily_flows( + station_number = id_station, + start_date = date_hydat_start + ) + # below used for test + # dplyr::filter( + # Date < "2018-01-01" | Date > "2018-01-03" + # ) + + stations_all <- tidyhydat::realtime_stations()[["STATION_NUMBER"]] + stations_no_rt <- setdiff(id_station, stations_all) + if (!all(id_station %in% stations_all)) { + cli::cli_alert_warning("There is no realtime daily flow data for {stations_no_rt}") + } + + stations_rt <- setdiff(id_station, stations_no_rt) + + if (length(stations_rt) > 0) { + hy_rt <- tidyhydat::realtime_ws( + station_number = stations_rt, + parameters = 6, + start_date = date_rt_start, + end_date = date_rt_end + ) + + hy_rt$Date <- as.Date(hy_rt$Date) + + hy_rt <- ngr::ngr_tidy_type(dat_w_types = hy_hy, dat_to_type = hy_rt) + + hy_all <- dplyr::bind_rows(hy_hy, hy_rt) + } + + if (length(stations_rt) == 0) { + hy_all <- hy_hy + } + + hy_all <- dplyr::distinct( + hy_all, + STATION_NUMBER, + Date, + .keep_all = TRUE + ) + + by_station <- split(hy_all$Date, hy_all$STATION_NUMBER) + by_station <- by_station[order(names(by_station))] + + bad <- character() + + for (stn in names(by_station)) { + ok <- ngr_chk_dt_complete( + by_station[[stn]], + units = "days", + ... + ) + + if (!ok) { + cli::cli_alert_warning("Station {stn} missing dates printed above") + bad <- c(bad, stn) + } + } + + return(hy_all) +} diff --git a/R/ngr_spk_stac_calc.R b/R/ngr_spk_stac_calc.R new file mode 100644 index 0000000..b8f1fcc --- /dev/null +++ b/R/ngr_spk_stac_calc.R @@ -0,0 +1,241 @@ +#' Spectral index calculation from a STAC item +#' +#' Compute a spectral index for a single STAC item. By default, the function reads +#' red and near-infrared (NIR) assets to compute NDVI, but this is only the +#' default configuration. The function is designed to read any two or three +#' assets and apply different spectral index calculations as they are added in +#' the future. +#' +#' The default asset names (`"red"` and `"nir08"`) align with common conventions +#' in Landsat Collection 2 Level-2 STAC items (e.g. the `landsat-c2-l2` collection +#' on Planetary Computer), but are fully parameterized to support other +#' collections with different band naming schemes. +#' +#' @param feature [list] A single STAC Item (one element from an `items$features` +#' list). +#' @param aoi [sf] or [NULL] Optional. An AOI object used to restrict reads (by +#' bbox) and to crop and mask the output. If `NULL`, the full assets are read and +#' no crop/mask is applied. Default is `NULL`. +#' @param asset_a [character] Optional. A single string giving the asset name for +#' the first input band. For NDVI this corresponds to the red band. Default is +#' `"red"`. +#' @param asset_b [character] Optional. A single string giving the asset name for +#' the second input band. For NDVI this corresponds to the NIR band. Default is +#' `"nir08"`. +#' @param asset_c [character] or [NULL] Optional. A single string giving the asset +#' name for an optional third input band, used by some calculations. Default is +#' `NULL`. +#' @param vsi_prefix [character] Optional. A single string giving the GDAL VSI +#' prefix used by [terra::rast()] to enable streaming, windowed reads over HTTP. +#' Exposed to allow alternative VSI backends (e.g. `/vsis3/`). Default is +#' `"/vsicurl/"`. +#' @param quiet [logical] Optional. If `TRUE`, suppress CLI messages. Default is +#' `FALSE`. +#' @param timing [logical] Optional. If `TRUE`, emit simple elapsed-time messages +#' for reads. Default is `FALSE`. +#' +#' @details +#' This function expects `feature` to provide the required `assets` referenced by +#' `asset_a`, `asset_b`, and (if used) `asset_c`. When `aoi` is not `NULL`, +#' `feature` must also have `properties$"proj:epsg"` so the AOI can be transformed +#' for windowed reads. +#' +#' The calculation performed is controlled by `calc`. By default, `calc = "ndvi"` +#' computes the Normalized Difference Vegetation Index using `(b - a) / (b + a)`. +#' This structure allows additional spectral indices (e.g. NDWI, EVI) to be added +#' in the future without changing the data access or cropping logic. +#' +#' When `aoi` is provided, reading is limited to the AOI bounding box in the +#' feature's projected CRS, and then cropped/masked to the AOI. When `aoi = NULL`, +#' the full assets are read and the calculated raster is returned for the full +#' raster extent. +#' +#' @return A [terra::SpatRaster] with calculated index values. +#' +#' @section Production-ready alternative: +#' This function is a proof of concept demonstrating how to perform raster +#' calculations on STAC items using `terra`. For production workflows involving +#' time series reductions, multi-image composites, or data cubes with varying +#' pixel sizes and coordinate reference systems, consider the +#' [gdalcubes](https://github.com/appelmar/gdalcubes) package. Key functions +#' include: +#' \itemize{ +#' \item `stac_image_collection()`: Create an image collection directly from +#' STAC query results with automatic band detection and VSI prefix handling. +#' \item `cube_view()`: Define spatiotemporal extent, resolution, aggregation, +#' and resampling method in a single specification. +#' \item `raster_cube()`: Build a data cube from an image collection, +#' automatically reprojecting and resampling images to a common grid. +#' \item `apply_pixel()`: Apply arithmetic expressions (e.g., +#' `"(nir - red) / (nir + red)"`) across all pixels. +#' \item `reduce_time()`: Temporal reductions (mean, median, max, etc.) to +#' composite multi-date imagery. +#' } +#' +#' @seealso [terra::rast()], [terra::crop()], [terra::mask()], [sf::st_transform()] +#' +#' @importFrom chk chk_flag chk_string +#' @importFrom cli cli_abort cli_alert_info +#' @importFrom sf st_crs st_transform +#' @importFrom terra crop ext mask rast vect +#' +#' @family spacehakr +#' @export +#' @examples +#' \dontrun{ +#' # STAC query against Planetary Computer (Landsat C2 L2) +#' stac_url <- "https://planetarycomputer.microsoft.com/api/stac/v1" +#' y <- 2000 +#' date_time <- paste0(y, "-05-01/", y, "-07-30") +#' +#' # Define an AOI from a bounding box (WGS84) +#' bbox <- c( +#' xmin = -126.55350240037997, +#' ymin = 54.4430453753869, +#' xmax = -126.52422763064457, +#' ymax = 54.46001902038006 +#' ) +#' +#' aoi <- sf::st_as_sfc(sf::st_bbox(bbox, crs = 4326)) |> +#' sf::st_as_sf() +#' +#' stac_query <- rstac::stac(stac_url) |> +#' rstac::stac_search( +#' collections = "landsat-c2-l2", +#' datetime = date_time, +#' intersects = sf::st_geometry(aoi)[[1]], +#' limit = 200 +#' ) |> +#' rstac::ext_filter(`eo:cloud_cover` <= 10) +#' +#' items <- stac_query |> +#' rstac::post_request() |> +#' rstac::items_fetch() |> +#' rstac::items_sign_planetary_computer() +#' +#' ndvi_list <- items$features |> +#' purrr::map(ngr_spk_stac_calc, aoi = aoi) +#' +#' ndvi_list <- ndvi_list |> +#' purrr::set_names(purrr::map_chr(items$features, "id")) +#' } + +ngr_spk_stac_calc <- function( + feature, + aoi = NULL, + asset_a = "red", + asset_b = "nir08", + asset_c = NULL, + calc = "ndvi", + vsi_prefix = "/vsicurl/", + quiet = FALSE, + timing = FALSE +) { + chk::chk_string(asset_a) + chk::chk_string(asset_b) + if (!is.null(asset_c)) chk::chk_string(asset_c) + chk::chk_string(calc) + chk::chk_string(vsi_prefix) + chk::chk_flag(quiet) + chk::chk_flag(timing) + + if (!is.list(feature) || is.null(feature$assets) || !is.list(feature$assets)) { + cli::cli_abort("`feature` must be a single STAC item (a list with an `assets` element).") + } + if (!is.null(aoi) && !inherits(aoi, "sf")) { + cli::cli_abort("`aoi` must be an sf object or `NULL`.") + } + + id <- feature$id %||% "" + + if (!is.null(aoi)) { + f_epsg <- feature[["properties"]][["proj:epsg"]] + if (is.null(f_epsg) || is.na(f_epsg)) { + cli::cli_abort("Missing `properties$proj:epsg` for item: {id}.") + } + } + + a_href <- feature$assets[[asset_a]]$href + b_href <- feature$assets[[asset_b]]$href + if (is.null(a_href) || !nzchar(a_href)) { + cli::cli_abort("Missing asset href for `asset_a` `{asset_a}` in item: {id}.") + } + if (is.null(b_href) || !nzchar(b_href)) { + cli::cli_abort("Missing asset href for `asset_b` `{asset_b}` in item: {id}.") + } + + a_url <- paste0(vsi_prefix, a_href) + b_url <- paste0(vsi_prefix, b_href) + + if (is.null(aoi)) { + if (!quiet) cli::cli_alert_info("aoi is NULL: reading full assets: {id}") + + if (!quiet) cli::cli_alert_info("read asset_a: {id}") + t0 <- if (timing) base::proc.time()[[3]] else NA_real_ + a <- terra::rast(a_url) + if (timing && !quiet) cli::cli_alert_info("read asset_a elapsed (s): {format(base::proc.time()[[3]] - t0, digits = 3)}") + + if (!quiet) cli::cli_alert_info("read asset_b: {id}") + t1 <- if (timing) base::proc.time()[[3]] else NA_real_ + b <- terra::rast(b_url) + if (timing && !quiet) cli::cli_alert_info("read asset_b elapsed (s): {format(base::proc.time()[[3]] - t1, digits = 3)}") + + out <- .ngr_spk_calc(calc, a = a, b = b, c = NULL) + return(out) + } + + aoi_epsg <- sf::st_transform(aoi, f_epsg) + win <- terra::ext(aoi_epsg) + + if (!quiet) cli::cli_alert_info("read asset_a: {id}") + t0 <- if (timing) base::proc.time()[[3]] else NA_real_ + a <- terra::rast(a_url, win = win) + if (timing && !quiet) cli::cli_alert_info("read asset_a elapsed (s): {format(base::proc.time()[[3]] - t0, digits = 3)}") + + if (!quiet) cli::cli_alert_info("read asset_b: {id}") + t1 <- if (timing) base::proc.time()[[3]] else NA_real_ + b <- terra::rast(b_url, win = win) + if (timing && !quiet) cli::cli_alert_info("read asset_b elapsed (s): {format(base::proc.time()[[3]] - t1, digits = 3)}") + + out <- .ngr_spk_calc(calc, a = a, b = b, c = NULL) + + aoi_v <- aoi |> + sf::st_transform(sf::st_crs(out)) |> + terra::vect() + + terra::crop(out, aoi_v) |> + terra::mask(aoi_v) +} + +# ---- helpers + +#' Calculation dispatcher +#' +#' @param calc [character] Name of the calculation. +#' @param a [terra::SpatRaster] First input raster. +#' @param b [terra::SpatRaster] Second input raster. +#' @param c [terra::SpatRaster] or [NULL] Optional third input raster. +#' +#' @return A [terra::SpatRaster]. +#' +#' @noRd +.ngr_spk_calc <- function(calc, a, b, c = NULL) { + if (calc == "ndvi") { + return((b - a) / (b + a)) + } + + cli::cli_abort("Unknown calc: {calc}.") +} + + +#' Null-coalescing helper +#' +#' @param x An object. +#' @param y An object to return if `x` is `NULL`. +#' +#' @return `x` if not `NULL`, otherwise `y`. +#' +#' @noRd +`%||%` <- function(x, y) { + if (is.null(x)) y else x +} diff --git a/R/ngr_str_df_extract.R b/R/ngr_str_df_extract.R new file mode 100644 index 0000000..194e073 --- /dev/null +++ b/R/ngr_str_df_extract.R @@ -0,0 +1,68 @@ +#' Extract segments of text from a data-frame column +#' +#' Pull substrings **between** regex start/end delimiters from a string column. +#' Each element of `segment_pairs` is a length-2 character vector: +#' `c(start_regex, end_regex)`. The **list name** becomes the output column name. +#' +#' @param data [data.frame] or [tibble::tbl_df] A table containing the source text column. +#' @param col [character] A single string: the **name of the text column** (e.g., "details"). +#' @param segment_pairs [list] A **named** list where each element is a length-2 +#' [character] vector `c(start_regex, end_regex)`. List names become output +#' column names. Names must be non-empty and unique. +#' +#' @returns [data.frame] The input `data` with one new column per named segment. +#' Values are `NA` where no match is found. An error is raised if `col` is +#' missing from `data`, if `segment_pairs` is not a properly named list, or if +#' any element is not a length-2 [character] vector. +#' +#' @seealso [ngr_str_extract_between()] to extract a single pair from a character vector. +#' +#' @family string dataframe +#' @export +#' +#' @importFrom cli cli_abort +#' +#' @examples +#' df <- tibble::tibble(details = c( +#' "Grant Amount: $400,000 Intake Year: 2025 Region: Fraser Basin Project Theme: Restoration", +#' "Grant Amount: $150,500 Intake Year: 2024 Region: Columbia Basin Project Theme: Food Systems" +#' )) +#' +#' segs <- list( +#' amount = c("Grant\\\\s*Amount", "Intake\\\\s*Year|Region|Project\\\\s*Theme|$"), +#' year = c("Intake\\\\s*Year", "Region|Project\\\\s*Theme|$"), +#' region = c("Region", "Project\\\\s*Theme|$"), +#' theme = c("Project\\\\s*Theme", "$") +#' ) +#' +#' out <- ngr_str_df_extract(df, "details", segs) +#' out +ngr_str_df_extract <- function(data, col, segment_pairs) { + # sanity checks + chk::chk_data(data) + chk::chk_string(col) + if (!col %in% names(data)) { + cli::cli_abort("ngr_str_df_extract: column '{col}' not found in `data`.") + } + if (!is.list(segment_pairs)) { + cli::cli_abort("ngr_str_df_extract: `segment_pairs` must be a named list of length-2 character vectors.") + } + if (is.null(names(segment_pairs))) { + cli::cli_abort("ngr_str_df_extract: `segment_pairs` must have names; names become output columns.") + } + nms <- names(segment_pairs) + if (any(nms == "")) cli::cli_abort("ngr_str_df_extract: all `segment_pairs` names must be non-empty.") + if (any(duplicated(nms))) cli::cli_abort("ngr_str_df_extract: names in `segment_pairs` must be unique.") + bad <- !vapply(segment_pairs, function(p) is.character(p) && length(p) == 2L, logical(1)) + if (any(bad)) cli::cli_abort("ngr_str_df_extract: each element must be `c(start_regex, end_regex)`.") + + + # extraction + for (nm in names(segment_pairs)) { + p <- segment_pairs[[nm]] + data[[nm]] <- ngr_str_extract_between(data[[col]], p[1], p[2]) + } + + + data +} diff --git a/R/ngr_str_extract_between.R b/R/ngr_str_extract_between.R new file mode 100644 index 0000000..be05119 --- /dev/null +++ b/R/ngr_str_extract_between.R @@ -0,0 +1,68 @@ +#' Extract text between two regex delimiters +#' +#' Pull a substring found **between** a start and end regular-expression pattern +#' from each element of a character vector. Matching is case-insensitive and +#' dot-all by default, and an optional colon after the start pattern is ignored +#' (e.g., `"Label:"`). You may optionally normalize internal whitespace. +#' +#' @param x [character] A character vector to search. +#' @param reg_start [character] A single string: the **start** regex pattern. +#' Optional trailing colon and whitespace in the source text are ignored. +#' @param reg_end [character] A single string: the **end** regex pattern used +#' in a lookahead; the matched text will end **before** this pattern. +#' @param squish [logical] Optional. If `TRUE`, collapse and trim whitespace in +#' the extracted text via [stringr::str_squish()]. Default is `FALSE`. +#' +#' @returns [character] A character vector the same length as `x`, with the +#' extracted substrings. Elements are `NA` when no match is found. Errors may +#' be thrown by the underlying regex engine if `reg_start` or `reg_end` +#' contain invalid regular expressions. +#' +#' @section Matching details: +#' * Flags used: `(?i)` case-insensitive, `(?s)` dot matches newline. +#' * Pattern built: non-capturing `(?:reg_start)` then optional `:\s*`, then +#' the first **non-greedy** capture `(.*?)`, ending **just before** +#' `(?:reg_end)` via a lookahead. If `squish = TRUE`, surrounding and internal +#' whitespace is normalized. +#' +#' @seealso [ngr_str_df_extract()] for applying multiple start/end pairs to a +#' data-frame column. +#' +#' @family string +#' +#' @export +#' +#' @importFrom stringr regex str_match str_squish +#' @importFrom chk chk_character chk_string chk_flag +#' +#' @examples +#' x <- c( +#' "Grant Amount: $400,000 Intake Year: 2025", +#' "Grant Amount: $150,500 Intake Year: 2024" +#' ) +#' ngr_str_extract_between(x, +#' reg_start = "Grant\\\\s*Amount", +#' reg_end = "Intake\\\\s*Year|$" +#' ) +#' +#' # With whitespace normalization +#' ngr_str_extract_between( +#' x = "Region : Fraser Basin Project Theme: Something", +#' reg_start = "Region", +#' reg_end = "Project\\\\s*Theme|$", +#' squish = TRUE +#' ) +ngr_str_extract_between <- function(x, reg_start, reg_end, squish = FALSE) { + chk::chk_character(x) + chk::chk_string(reg_start) + chk::chk_string(reg_end) + chk::chk_flag(squish) + + + pat <- stringr::regex( + paste0("(?is)(?:", reg_start, ")\\s*(.*?)\\s*(?=(?:", reg_end, "))") + ) + out <- stringr::str_match(x, pat)[, 2] + if (squish) out <- stringr::str_squish(out) + ifelse(is.na(out), NA_character_, out) +} diff --git a/man/ngr_chk_coerce_date.Rd b/man/ngr_chk_coerce_date.Rd new file mode 100644 index 0000000..cfe9fb2 --- /dev/null +++ b/man/ngr_chk_coerce_date.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ngr_chk_coerce_date.R +\name{ngr_chk_coerce_date} +\alias{ngr_chk_coerce_date} +\title{Check and Validate Date-Coercible Input} +\usage{ +ngr_chk_coerce_date(x, x_name = NULL) +} +\arguments{ +\item{x}{\link{any} An object expected to be a \link{Date} or a string coercible to a +\link{Date}.} + +\item{x_name}{\link{character} Optional. A single string giving the name of \code{x} +used in error messages. If \code{NULL}, the name is inferred from the calling +expression.} +} +\value{ +\code{x} invisibly if it is a \link{Date} or coercible to one. An error is thrown if +coercion fails. +} +\description{ +Checks whether an object is a \link{Date} or can be safely coerced to a \link{Date} +using \code{\link[base:as.Date]{base::as.Date()}}. If the input is coercible, it is returned +invisibly; otherwise, an informative error is thrown. +} +\details{ +This helper is intended for lightweight validation in functions that accept +date inputs but do not need to modify them. Coercion is attempted via +\code{\link[base:as.Date]{base::as.Date()}}, and failure results in a call to +\code{\link[chk:abort_chk]{chk::abort_chk()}}. +} +\seealso{ +\code{\link[base:as.Date]{base::as.Date()}}, \code{\link[chk:abort_chk]{chk::abort_chk()}} + +Other chk: +\code{\link{ngr_chk_dt_complete}()} +} +\concept{chk} diff --git a/man/ngr_chk_dt_complete.Rd b/man/ngr_chk_dt_complete.Rd new file mode 100644 index 0000000..8a4536e --- /dev/null +++ b/man/ngr_chk_dt_complete.Rd @@ -0,0 +1,58 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ngr_chk_dt_complete.R +\name{ngr_chk_dt_complete} +\alias{ngr_chk_dt_complete} +\title{Check for Complete Date Sequences} +\usage{ +ngr_chk_dt_complete( + x, + units = "days", + dates_print = TRUE, + dates_capture = FALSE +) +} +\arguments{ +\item{x}{\link{Date} A vector of dates to check for completeness.} + +\item{units}{\link{character} Optional. A single string giving the time unit +for the sequence passed to \code{\link[base:seq.Date]{base::seq.Date()}}. Default is \code{"days"}.} + +\item{dates_print}{\link{logical} Optional. Whether to print missing dates to the +console using \code{dput()} when gaps are detected. Default is \code{TRUE}.} + +\item{dates_capture}{\link{logical} Optional. Whether to capture and return the +missing dates as an object when gaps are detected. Default is \code{FALSE}.} +} +\value{ +If \code{dates_capture = FALSE} (default), returns \code{TRUE} invisibly if no dates +are missing and \code{FALSE} invisibly if one or more dates are missing. + +If \code{dates_capture = TRUE}, returns the missing dates as a \link{Date} vector +(invisibly) when gaps are detected, and \code{TRUE} invisibly when no dates +are missing. +} +\description{ +Checks whether a vector of dates forms a complete, uninterrupted sequence +between its minimum and maximum values. If missing dates are detected, +they are printed using \code{dput()} for easy reuse and the function returns +\code{FALSE} invisibly. +} +\details{ +The function computes a full date sequence using \code{\link[base:seq.Date]{base::seq.Date()}} from +\code{min(x)} to \code{max(x)} and compares it to the unique values of \code{x} using +\code{\link[base:sets]{base::setdiff()}}. Missing dates, if any, are printed with \code{dput()} so they +can be copied directly into code or tests. +} +\examples{ +dates_ok <- as.Date("2024-01-01") + 0:4 +ngr_chk_dt_complete(dates_ok) + +dates_bad <- as.Date(c("2024-01-01", "2024-01-02", "2024-01-04")) +ngr_chk_dt_complete(dates_bad) + +} +\seealso{ +Other chk: +\code{\link{ngr_chk_coerce_date}()} +} +\concept{chk} diff --git a/man/ngr_hyd_q_daily.Rd b/man/ngr_hyd_q_daily.Rd new file mode 100644 index 0000000..e140558 --- /dev/null +++ b/man/ngr_hyd_q_daily.Rd @@ -0,0 +1,60 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ngr_hyd_q_daily.R +\name{ngr_hyd_q_daily} +\alias{ngr_hyd_q_daily} +\title{Combine daily Streamflow Data from HYDAT and ECCC Realtime Sources} +\usage{ +ngr_hyd_q_daily( + id_station, + date_rt_start = Sys.Date() - 581, + date_rt_end = Sys.Date(), + date_hydat_start = "1980-01-01", + ... +) +} +\arguments{ +\item{id_station}{\link{character} A character vector of one or more HYDAT station IDs} + +\item{date_rt_start}{\link{Date} Optional. A single start date (or string coercible to date) for realtime data retrieval. +Defaults to \code{Sys.Date() - 581} (18 months).} + +\item{date_rt_end}{\link{Date} Optional. A single end date (or string coercible to date) for realtime data retrieval. +Defaults to \code{Sys.Date()}.} + +\item{date_hydat_start}{\link{character} Optional. A single date (or string coercible to date) giving the start +date for historical HYDAT daily flows. Default is "1980-01-01".} + +\item{...}{Optional. Additional arguments passed to +\code{\link[=ngr_chk_dt_complete]{ngr_chk_dt_complete()}} when checking for missing dates within each +\code{STATION_NUMBER} time series (ex print_missing = FALSE).} +} +\value{ +A data frame containing daily streamflow data for the requested stations. +Warnings are emitted when realtime data are unavailable for some stations or +when missing dates are detected. +} +\description{ +Retrieve and amalgamate daily streamflow data for one or more HYDAT stations, +combining historical daily flows from HYDAT with available realtime daily flow +data where supported. The function warns when realtime data are unavailable and +checks for missing dates in the combined time series. +} +\details{ +Historical daily flows are retrieved using +\code{\link[tidyhydat:hy_daily_flows]{tidyhydat::hy_daily_flows()}}. Realtime daily flows are retrieved using +\code{\link[tidyhydat:realtime_ws]{tidyhydat::realtime_ws()}} for stations that support realtime reporting, as +identified by \code{\link[tidyhydat:realtime_stations]{tidyhydat::realtime_stations()}}. Realtime data are coerced to +daily resolution and typed to match the HYDAT data using +\code{\link[=ngr_tidy_type]{ngr_tidy_type()}}. + +Input dates are validated using \code{\link[=ngr_chk_coerce_date]{ngr_chk_coerce_date()}}. Missing daily +timestamps are checked per station using \code{\link[=ngr_chk_dt_complete]{ngr_chk_dt_complete()}}. + +Duplicate station-date records are removed, keeping historical HYDAT values +where overlaps occur. The function reports stations with missing dates using +cli messages. +} +\examples{ +ngr_hyd_q_daily(id_station = c("08NH118", "08NH126")) + +} diff --git a/man/ngr_spk_gdalwarp.Rd b/man/ngr_spk_gdalwarp.Rd index e2427b0..5c19766 100644 --- a/man/ngr_spk_gdalwarp.Rd +++ b/man/ngr_spk_gdalwarp.Rd @@ -92,7 +92,8 @@ Other spacehakr: \code{\link{ngr_spk_rast_ext}()}, \code{\link{ngr_spk_rast_not_empty}()}, \code{\link{ngr_spk_rast_rm_empty}()}, -\code{\link{ngr_spk_res}()} +\code{\link{ngr_spk_res}()}, +\code{\link{ngr_spk_stac_calc}()} Other processx: \code{\link{ngr_spk_odm}()} diff --git a/man/ngr_spk_join.Rd b/man/ngr_spk_join.Rd index 2b7d2c4..d480dd0 100644 --- a/man/ngr_spk_join.Rd +++ b/man/ngr_spk_join.Rd @@ -62,6 +62,7 @@ Other spacehakr: \code{\link{ngr_spk_rast_ext}()}, \code{\link{ngr_spk_rast_not_empty}()}, \code{\link{ngr_spk_rast_rm_empty}()}, -\code{\link{ngr_spk_res}()} +\code{\link{ngr_spk_res}()}, +\code{\link{ngr_spk_stac_calc}()} } \concept{spacehakr} diff --git a/man/ngr_spk_layer_info.Rd b/man/ngr_spk_layer_info.Rd index d8ce56d..1d874be 100644 --- a/man/ngr_spk_layer_info.Rd +++ b/man/ngr_spk_layer_info.Rd @@ -35,6 +35,7 @@ Other spacehakr: \code{\link{ngr_spk_rast_ext}()}, \code{\link{ngr_spk_rast_not_empty}()}, \code{\link{ngr_spk_rast_rm_empty}()}, -\code{\link{ngr_spk_res}()} +\code{\link{ngr_spk_res}()}, +\code{\link{ngr_spk_stac_calc}()} } \concept{spacehakr} diff --git a/man/ngr_spk_odm.Rd b/man/ngr_spk_odm.Rd index 0c7bcbb..fec02c7 100644 --- a/man/ngr_spk_odm.Rd +++ b/man/ngr_spk_odm.Rd @@ -90,7 +90,8 @@ Other spacehakr: \code{\link{ngr_spk_rast_ext}()}, \code{\link{ngr_spk_rast_not_empty}()}, \code{\link{ngr_spk_rast_rm_empty}()}, -\code{\link{ngr_spk_res}()} +\code{\link{ngr_spk_res}()}, +\code{\link{ngr_spk_stac_calc}()} Other processx: \code{\link{ngr_spk_gdalwarp}()} @@ -104,7 +105,8 @@ Other spacehakr: \code{\link{ngr_spk_rast_ext}()}, \code{\link{ngr_spk_rast_not_empty}()}, \code{\link{ngr_spk_rast_rm_empty}()}, -\code{\link{ngr_spk_res}()} +\code{\link{ngr_spk_res}()}, +\code{\link{ngr_spk_stac_calc}()} } \concept{odm} \concept{processx} diff --git a/man/ngr_spk_poly_to_points.Rd b/man/ngr_spk_poly_to_points.Rd index b190da6..541bc99 100644 --- a/man/ngr_spk_poly_to_points.Rd +++ b/man/ngr_spk_poly_to_points.Rd @@ -46,6 +46,7 @@ Other spacehakr: \code{\link{ngr_spk_rast_ext}()}, \code{\link{ngr_spk_rast_not_empty}()}, \code{\link{ngr_spk_rast_rm_empty}()}, -\code{\link{ngr_spk_res}()} +\code{\link{ngr_spk_res}()}, +\code{\link{ngr_spk_stac_calc}()} } \concept{spacehakr} diff --git a/man/ngr_spk_q_layer_info.Rd b/man/ngr_spk_q_layer_info.Rd index b2ea7f6..e71d7e6 100644 --- a/man/ngr_spk_q_layer_info.Rd +++ b/man/ngr_spk_q_layer_info.Rd @@ -49,6 +49,7 @@ Other spacehakr: \code{\link{ngr_spk_rast_ext}()}, \code{\link{ngr_spk_rast_not_empty}()}, \code{\link{ngr_spk_rast_rm_empty}()}, -\code{\link{ngr_spk_res}()} +\code{\link{ngr_spk_res}()}, +\code{\link{ngr_spk_stac_calc}()} } \concept{spacehakr} diff --git a/man/ngr_spk_rast_ext.Rd b/man/ngr_spk_rast_ext.Rd index c348f1a..95bae1b 100644 --- a/man/ngr_spk_rast_ext.Rd +++ b/man/ngr_spk_rast_ext.Rd @@ -61,6 +61,7 @@ Other spacehakr: \code{\link{ngr_spk_q_layer_info}()}, \code{\link{ngr_spk_rast_not_empty}()}, \code{\link{ngr_spk_rast_rm_empty}()}, -\code{\link{ngr_spk_res}()} +\code{\link{ngr_spk_res}()}, +\code{\link{ngr_spk_stac_calc}()} } \concept{spacehakr} diff --git a/man/ngr_spk_rast_not_empty.Rd b/man/ngr_spk_rast_not_empty.Rd index 13fff04..d0d30d4 100644 --- a/man/ngr_spk_rast_not_empty.Rd +++ b/man/ngr_spk_rast_not_empty.Rd @@ -25,6 +25,7 @@ Other spacehakr: \code{\link{ngr_spk_q_layer_info}()}, \code{\link{ngr_spk_rast_ext}()}, \code{\link{ngr_spk_rast_rm_empty}()}, -\code{\link{ngr_spk_res}()} +\code{\link{ngr_spk_res}()}, +\code{\link{ngr_spk_stac_calc}()} } \concept{spacehakr} diff --git a/man/ngr_spk_rast_rm_empty.Rd b/man/ngr_spk_rast_rm_empty.Rd index 4f1ebcb..4098034 100644 --- a/man/ngr_spk_rast_rm_empty.Rd +++ b/man/ngr_spk_rast_rm_empty.Rd @@ -39,6 +39,7 @@ Other spacehakr: \code{\link{ngr_spk_q_layer_info}()}, \code{\link{ngr_spk_rast_ext}()}, \code{\link{ngr_spk_rast_not_empty}()}, -\code{\link{ngr_spk_res}()} +\code{\link{ngr_spk_res}()}, +\code{\link{ngr_spk_stac_calc}()} } \concept{spacehakr} diff --git a/man/ngr_spk_res.Rd b/man/ngr_spk_res.Rd index be38e1f..698f449 100644 --- a/man/ngr_spk_res.Rd +++ b/man/ngr_spk_res.Rd @@ -41,6 +41,7 @@ Other spacehakr: \code{\link{ngr_spk_q_layer_info}()}, \code{\link{ngr_spk_rast_ext}()}, \code{\link{ngr_spk_rast_not_empty}()}, -\code{\link{ngr_spk_rast_rm_empty}()} +\code{\link{ngr_spk_rast_rm_empty}()}, +\code{\link{ngr_spk_stac_calc}()} } \concept{spacehakr} diff --git a/man/ngr_spk_stac_calc.Rd b/man/ngr_spk_stac_calc.Rd new file mode 100644 index 0000000..039b8cc --- /dev/null +++ b/man/ngr_spk_stac_calc.Rd @@ -0,0 +1,159 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ngr_spk_stac_calc.R +\name{ngr_spk_stac_calc} +\alias{ngr_spk_stac_calc} +\title{Spectral index calculation from a STAC item} +\usage{ +ngr_spk_stac_calc( + feature, + aoi = NULL, + asset_a = "red", + asset_b = "nir08", + asset_c = NULL, + calc = "ndvi", + vsi_prefix = "/vsicurl/", + quiet = FALSE, + timing = FALSE +) +} +\arguments{ +\item{feature}{\link{list} A single STAC Item (one element from an \code{items$features} +list).} + +\item{aoi}{\link[sf:sf]{sf::sf} or \link{NULL} Optional. An AOI object used to restrict reads (by +bbox) and to crop and mask the output. If \code{NULL}, the full assets are read and +no crop/mask is applied. Default is \code{NULL}.} + +\item{asset_a}{\link{character} Optional. A single string giving the asset name for +the first input band. For NDVI this corresponds to the red band. Default is +\code{"red"}.} + +\item{asset_b}{\link{character} Optional. A single string giving the asset name for +the second input band. For NDVI this corresponds to the NIR band. Default is +\code{"nir08"}.} + +\item{asset_c}{\link{character} or \link{NULL} Optional. A single string giving the asset +name for an optional third input band, used by some calculations. Default is +\code{NULL}.} + +\item{vsi_prefix}{\link{character} Optional. A single string giving the GDAL VSI +prefix used by \code{\link[terra:rast]{terra::rast()}} to enable streaming, windowed reads over HTTP. +Exposed to allow alternative VSI backends (e.g. \verb{/vsis3/}). Default is +\code{"/vsicurl/"}.} + +\item{quiet}{\link{logical} Optional. If \code{TRUE}, suppress CLI messages. Default is +\code{FALSE}.} + +\item{timing}{\link{logical} Optional. If \code{TRUE}, emit simple elapsed-time messages +for reads. Default is \code{FALSE}.} +} +\value{ +A \link[terra:SpatRaster-class]{terra::SpatRaster} with calculated index values. +} +\description{ +Compute a spectral index for a single STAC item. By default, the function reads +red and near-infrared (NIR) assets to compute NDVI, but this is only the +default configuration. The function is designed to read any two or three +assets and apply different spectral index calculations as they are added in +the future. +} +\details{ +The default asset names (\code{"red"} and \code{"nir08"}) align with common conventions +in Landsat Collection 2 Level-2 STAC items (e.g. the \code{landsat-c2-l2} collection +on Planetary Computer), but are fully parameterized to support other +collections with different band naming schemes. + +This function expects \code{feature} to provide the required \code{assets} referenced by +\code{asset_a}, \code{asset_b}, and (if used) \code{asset_c}. When \code{aoi} is not \code{NULL}, +\code{feature} must also have \code{properties$"proj:epsg"} so the AOI can be transformed +for windowed reads. + +The calculation performed is controlled by \code{calc}. By default, \code{calc = "ndvi"} +computes the Normalized Difference Vegetation Index using \code{(b - a) / (b + a)}. +This structure allows additional spectral indices (e.g. NDWI, EVI) to be added +in the future without changing the data access or cropping logic. + +When \code{aoi} is provided, reading is limited to the AOI bounding box in the +feature's projected CRS, and then cropped/masked to the AOI. When \code{aoi = NULL}, +the full assets are read and the calculated raster is returned for the full +raster extent. +} +\section{Production-ready alternative}{ + +This function is a proof of concept demonstrating how to perform raster +calculations on STAC items using \code{terra}. For production workflows involving +time series reductions, multi-image composites, or data cubes with varying +pixel sizes and coordinate reference systems, consider the +\href{https://github.com/appelmar/gdalcubes}{gdalcubes} package. Key functions +include: +\itemize{ +\item \code{\link[gdalcubes:stac_image_collection]{gdalcubes::stac_image_collection()}}: Create an image collection +directly from STAC query results with automatic band detection and VSI +prefix handling. +\item \code{\link[gdalcubes:cube_view]{gdalcubes::cube_view()}}: Define spatiotemporal extent, resolution, +aggregation, and resampling method in a single specification. +\item \code{\link[gdalcubes:raster_cube]{gdalcubes::raster_cube()}}: Build a data cube from an image +collection, automatically reprojecting and resampling images to a common +grid. +\item \code{\link[gdalcubes:apply_pixel]{gdalcubes::apply_pixel()}}: Apply arithmetic expressions (e.g., +\code{"(nir - red) / (nir + red)"}) across all pixels. +\item \code{\link[gdalcubes:reduce_time]{gdalcubes::reduce_time()}}: Temporal reductions (mean, median, max, +etc.) to composite multi-date imagery. +} +} + +\examples{ +\dontrun{ +# STAC query against Planetary Computer (Landsat C2 L2) +stac_url <- "https://planetarycomputer.microsoft.com/api/stac/v1" +y <- 2000 +date_time <- paste0(y, "-05-01/", y, "-07-30") + +# Define an AOI from a bounding box (WGS84) +bbox <- c( + xmin = -126.55350240037997, + ymin = 54.4430453753869, + xmax = -126.52422763064457, + ymax = 54.46001902038006 +) + +aoi <- sf::st_as_sfc(sf::st_bbox(bbox, crs = 4326)) |> + sf::st_as_sf() + +stac_query <- rstac::stac(stac_url) |> + rstac::stac_search( + collections = "landsat-c2-l2", + datetime = date_time, + intersects = sf::st_geometry(aoi)[[1]], + limit = 200 + ) |> + rstac::ext_filter(`eo:cloud_cover` <= 10) + +items <- stac_query |> + rstac::post_request() |> + rstac::items_fetch() |> + rstac::items_sign_planetary_computer() + +ndvi_list <- items$features |> + purrr::map(ngr_spk_stac_calc, aoi = aoi) + +ndvi_list <- ndvi_list |> + purrr::set_names(purrr::map_chr(items$features, "id")) +} +} +\seealso{ +\code{\link[terra:rast]{terra::rast()}}, \code{\link[terra:crop]{terra::crop()}}, \code{\link[terra:mask]{terra::mask()}}, \code{\link[sf:st_transform]{sf::st_transform()}} + +Other spacehakr: +\code{\link{ngr_spk_gdalwarp}()}, +\code{\link{ngr_spk_join}()}, +\code{\link{ngr_spk_layer_info}()}, +\code{\link{ngr_spk_odm}()}, +\code{\link{ngr_spk_poly_to_points}()}, +\code{\link{ngr_spk_q_layer_info}()}, +\code{\link{ngr_spk_rast_ext}()}, +\code{\link{ngr_spk_rast_not_empty}()}, +\code{\link{ngr_spk_rast_rm_empty}()}, +\code{\link{ngr_spk_res}()} +} +\concept{spacehakr} diff --git a/man/ngr_str_df_col_agg.Rd b/man/ngr_str_df_col_agg.Rd index 749814f..fc537d3 100644 --- a/man/ngr_str_df_col_agg.Rd +++ b/man/ngr_str_df_col_agg.Rd @@ -121,6 +121,7 @@ head(dat_subset) \code{\link[=mean]{mean()}}, \code{\link[=sum]{sum()}}, \code{\link[=median]{median()}}, \code{\link[=min]{min()}}, \code{\link[=max]{max()}} Other string dataframe: -\code{\link{ngr_str_df_detect_filter}()} +\code{\link{ngr_str_df_detect_filter}()}, +\code{\link{ngr_str_df_extract}()} } \concept{string dataframe} diff --git a/man/ngr_str_df_detect_filter.Rd b/man/ngr_str_df_detect_filter.Rd index 47fbe6f..0e5372f 100644 --- a/man/ngr_str_df_detect_filter.Rd +++ b/man/ngr_str_df_detect_filter.Rd @@ -21,6 +21,7 @@ This function filters rows in a data frame where a specified column's name match } \seealso{ Other string dataframe: -\code{\link{ngr_str_df_col_agg}()} +\code{\link{ngr_str_df_col_agg}()}, +\code{\link{ngr_str_df_extract}()} } \concept{string dataframe} diff --git a/man/ngr_str_df_extract.Rd b/man/ngr_str_df_extract.Rd new file mode 100644 index 0000000..424e3b6 --- /dev/null +++ b/man/ngr_str_df_extract.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ngr_str_df_extract.R +\name{ngr_str_df_extract} +\alias{ngr_str_df_extract} +\title{Extract segments of text from a data-frame column} +\usage{ +ngr_str_df_extract(data, col, segment_pairs) +} +\arguments{ +\item{data}{\link{data.frame} or \link[tibble:tbl_df-class]{tibble::tbl_df} A table containing the source text column.} + +\item{col}{\link{character} A single string: the \strong{name of the text column} (e.g., "details").} + +\item{segment_pairs}{\link{list} A \strong{named} list where each element is a length-2 +\link{character} vector \code{c(start_regex, end_regex)}. List names become output +column names. Names must be non-empty and unique.} +} +\value{ +\link{data.frame} The input \code{data} with one new column per named segment. +Values are \code{NA} where no match is found. An error is raised if \code{col} is +missing from \code{data}, if \code{segment_pairs} is not a properly named list, or if +any element is not a length-2 \link{character} vector. +} +\description{ +Pull substrings \strong{between} regex start/end delimiters from a string column. +Each element of \code{segment_pairs} is a length-2 character vector: +\code{c(start_regex, end_regex)}. The \strong{list name} becomes the output column name. +} +\examples{ +df <- tibble::tibble(details = c( +"Grant Amount: $400,000 Intake Year: 2025 Region: Fraser Basin Project Theme: Restoration", +"Grant Amount: $150,500 Intake Year: 2024 Region: Columbia Basin Project Theme: Food Systems" +)) + +segs <- list( +amount = c("Grant\\\\\\\\s*Amount", "Intake\\\\\\\\s*Year|Region|Project\\\\\\\\s*Theme|$"), +year = c("Intake\\\\\\\\s*Year", "Region|Project\\\\\\\\s*Theme|$"), +region = c("Region", "Project\\\\\\\\s*Theme|$"), +theme = c("Project\\\\\\\\s*Theme", "$") +) + +out <- ngr_str_df_extract(df, "details", segs) +out +} +\seealso{ +\code{\link[=ngr_str_extract_between]{ngr_str_extract_between()}} to extract a single pair from a character vector. + +Other string dataframe: +\code{\link{ngr_str_df_col_agg}()}, +\code{\link{ngr_str_df_detect_filter}()} +} +\concept{string dataframe} diff --git a/man/ngr_str_dir_from_path.Rd b/man/ngr_str_dir_from_path.Rd index 48fc7d6..9b9346f 100644 --- a/man/ngr_str_dir_from_path.Rd +++ b/man/ngr_str_dir_from_path.Rd @@ -25,6 +25,7 @@ ngr_str_dir_from_path("~/Projects/repo/ngr/data-raw/extdata.R", levels = 2) } \seealso{ Other string: +\code{\link{ngr_str_extract_between}()}, \code{\link{ngr_str_link_url}()} } \concept{string} diff --git a/man/ngr_str_extract_between.Rd b/man/ngr_str_extract_between.Rd new file mode 100644 index 0000000..ea9b48a --- /dev/null +++ b/man/ngr_str_extract_between.Rd @@ -0,0 +1,70 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ngr_str_extract_between.R +\name{ngr_str_extract_between} +\alias{ngr_str_extract_between} +\title{Extract text between two regex delimiters} +\usage{ +ngr_str_extract_between(x, reg_start, reg_end, squish = FALSE) +} +\arguments{ +\item{x}{\link{character} A character vector to search.} + +\item{reg_start}{\link{character} A single string: the \strong{start} regex pattern. +Optional trailing colon and whitespace in the source text are ignored.} + +\item{reg_end}{\link{character} A single string: the \strong{end} regex pattern used +in a lookahead; the matched text will end \strong{before} this pattern.} + +\item{squish}{\link{logical} Optional. If \code{TRUE}, collapse and trim whitespace in +the extracted text via \code{\link[stringr:str_trim]{stringr::str_squish()}}. Default is \code{FALSE}.} +} +\value{ +\link{character} A character vector the same length as \code{x}, with the +extracted substrings. Elements are \code{NA} when no match is found. Errors may +be thrown by the underlying regex engine if \code{reg_start} or \code{reg_end} +contain invalid regular expressions. +} +\description{ +Pull a substring found \strong{between} a start and end regular-expression pattern +from each element of a character vector. Matching is case-insensitive and +dot-all by default, and an optional colon after the start pattern is ignored +(e.g., \code{"Label:"}). You may optionally normalize internal whitespace. +} +\section{Matching details}{ + +\itemize{ +\item Flags used: \code{(?i)} case-insensitive, \code{(?s)} dot matches newline. +\item Pattern built: non-capturing \verb{(?:reg_start)} then optional \verb{:\\s*}, then +the first \strong{non-greedy} capture \verb{(.*?)}, ending \strong{just before} +\verb{(?:reg_end)} via a lookahead. If \code{squish = TRUE}, surrounding and internal +whitespace is normalized. +} +} + +\examples{ +x <- c( +"Grant Amount: $400,000 Intake Year: 2025", +"Grant Amount: $150,500 Intake Year: 2024" +) +ngr_str_extract_between(x, +reg_start = "Grant\\\\\\\\s*Amount", +reg_end = "Intake\\\\\\\\s*Year|$" +) + +# With whitespace normalization +ngr_str_extract_between( +x = "Region : Fraser Basin Project Theme: Something", +reg_start = "Region", +reg_end = "Project\\\\\\\\s*Theme|$", +squish = TRUE +) +} +\seealso{ +\code{\link[=ngr_str_df_extract]{ngr_str_df_extract()}} for applying multiple start/end pairs to a +data-frame column. + +Other string: +\code{\link{ngr_str_dir_from_path}()}, +\code{\link{ngr_str_link_url}()} +} +\concept{string} diff --git a/man/ngr_str_link_url.Rd b/man/ngr_str_link_url.Rd index a60caa5..36c3195 100644 --- a/man/ngr_str_link_url.Rd +++ b/man/ngr_str_link_url.Rd @@ -66,6 +66,7 @@ print(df$link) } \seealso{ Other string: -\code{\link{ngr_str_dir_from_path}()} +\code{\link{ngr_str_dir_from_path}()}, +\code{\link{ngr_str_extract_between}()} } \concept{string} diff --git a/tests/testthat/test-ngr_chk_dt_complete.R b/tests/testthat/test-ngr_chk_dt_complete.R new file mode 100644 index 0000000..a115492 --- /dev/null +++ b/tests/testthat/test-ngr_chk_dt_complete.R @@ -0,0 +1,20 @@ +# test data +dates_ok <- as.Date("2024-01-01") + 0:4 +dates_bad <- as.Date(c("2024-01-01", "2024-01-02", "2024-01-04")) + + +testthat::test_that("ngr_chk_dt_complete returns nothing for complete date sequences", { + res <- ngr_chk_dt_complete(dates_ok) + testthat::expect_true(res) +}) + +testthat::test_that("ngr_chk_dt_complete returns FALSE for incomplete date sequences", { + res <- ngr_chk_dt_complete(dates_bad) + testthat::expect_false(res) +}) + +testthat::test_that("ngr_chk_dt_complete returns dates when dates_capture = TRUE", { + res <- ngr_chk_dt_complete(dates_bad, dates_capture = TRUE) + testthat::expect_identical(res, as.Date("2024-01-03")) + +}) diff --git a/tests/testthat/test-ngr_spk_gdalwarp.R b/tests/testthat/test-ngr_spk_gdalwarp.R index ded6204..3b4b0dc 100644 --- a/tests/testthat/test-ngr_spk_gdalwarp.R +++ b/tests/testthat/test-ngr_spk_gdalwarp.R @@ -25,8 +25,8 @@ expected_args <- c( file_out ) -test_that("ngr_spk_gdalwarp constructs correct arguments", { - expect_equal(args, expected_args) +testthat::test_that("ngr_spk_gdalwarp constructs correct arguments", { + testthat::expect_equal(args, expected_args) }) diff --git a/tests/testthat/test-ngr_spk_stac_calc.R b/tests/testthat/test-ngr_spk_stac_calc.R new file mode 100644 index 0000000..864562a --- /dev/null +++ b/tests/testthat/test-ngr_spk_stac_calc.R @@ -0,0 +1,188 @@ +# Test input validation -------------------------------------------------------- + +testthat::test_that("ngr_spk_stac_calc errors when feature is not a list with assets", { + + testthat::expect_error( + ngr::ngr_spk_stac_calc(feature = "not_a_list"), + "must be a single STAC item" + ) + + testthat::expect_error( + ngr::ngr_spk_stac_calc(feature = list(no_assets = TRUE)), + "must be a single STAC item" + ) + + testthat::expect_error( + ngr::ngr_spk_stac_calc(feature = list(assets = "not_a_list")), + "must be a single STAC item" + ) +}) + +testthat::test_that("ngr_spk_stac_calc errors when aoi is not sf or NULL", { + mock_feature <- list( + id = "test_item", + assets = list( + red = list(href = "https://example.com/red.tif"), + nir08 = list(href = "https://example.com/nir.tif") + ) + ) + + testthat::expect_error( + ngr::ngr_spk_stac_calc(feature = mock_feature, aoi = "not_sf"), + "must be an sf object" + ) + + testthat::expect_error( + ngr::ngr_spk_stac_calc(feature = mock_feature, aoi = data.frame(x = 1)), + "must be an sf object" + ) +}) + +testthat::test_that("ngr_spk_stac_calc errors when required assets are missing", { + mock_feature <- list( + id = "test_item", + assets = list( + blue = list(href = "https://example.com/blue.tif") + ) + ) + + testthat::expect_error( + ngr::ngr_spk_stac_calc(feature = mock_feature, quiet = TRUE), + "Missing asset href.*asset_a.*red" + ) + + mock_feature2 <- list( + id = "test_item", + assets = list( + red = list(href = "https://example.com/red.tif") + ) + ) + + testthat::expect_error( + ngr::ngr_spk_stac_calc(feature = mock_feature2, quiet = TRUE), + "Missing asset href.*asset_b.*nir08" + ) +}) + +testthat::test_that("ngr_spk_stac_calc errors when asset href is empty", { + mock_feature <- list( + id = "test_item", + assets = list( + red = list(href = ""), + nir08 = list(href = "https://example.com/nir.tif") + ) + ) + + testthat::expect_error( + ngr::ngr_spk_stac_calc(feature = mock_feature, quiet = TRUE), + "Missing asset href.*asset_a.*red" + ) +}) + +testthat::test_that("ngr_spk_stac_calc errors when aoi provided but proj:epsg missing", { + mock_feature <- list( + id = "test_item", + properties = list(), + assets = list( + red = list(href = "https://example.com/red.tif"), + nir08 = list(href = "https://example.com/nir.tif") + ) + ) + + aoi <- sf::st_as_sfc(sf::st_bbox(c(xmin = 0, ymin = 0, xmax = 1, ymax = 1), crs = 4326)) |> + sf::st_as_sf() + + testthat::expect_error( + ngr::ngr_spk_stac_calc(feature = mock_feature, aoi = aoi, quiet = TRUE), + "Missing.*proj:epsg" + ) +}) + +# Test parameter validation ---------------------------------------------------- + +testthat::test_that("ngr_spk_stac_calc validates string parameters", { + mock_feature <- list( + id = "test_item", + assets = list( + red = list(href = "https://example.com/red.tif"), + nir08 = list(href = "https://example.com/nir.tif") + ) + ) + + testthat::expect_error(ngr::ngr_spk_stac_calc(feature = mock_feature, asset_a = 123)) + testthat::expect_error(ngr::ngr_spk_stac_calc(feature = mock_feature, asset_b = 123)) + testthat::expect_error(ngr::ngr_spk_stac_calc(feature = mock_feature, calc = 123)) + testthat::expect_error(ngr::ngr_spk_stac_calc(feature = mock_feature, vsi_prefix = 123)) +}) + +testthat::test_that("ngr_spk_stac_calc validates flag parameters", { + mock_feature <- list( + id = "test_item", + assets = list( + red = list(href = "https://example.com/red.tif"), + nir08 = list(href = "https://example.com/nir.tif") + ) + ) + + testthat::expect_error(ngr::ngr_spk_stac_calc(feature = mock_feature, quiet = "yes")) + testthat::expect_error(ngr::ngr_spk_stac_calc(feature = mock_feature, timing = "yes")) +}) + +# Test internal .ngr_spk_calc helper ------------------------------------------- + +testthat::test_that(".ngr_spk_calc computes NDVI correctly", { + testthat::skip_if_not_installed("terra") + + # Create synthetic rasters with known values + a <- terra::rast(nrows = 2, ncols = 2, vals = c(100, 200, 150, 50)) + b <- terra::rast(nrows = 2, ncols = 2, vals = c(300, 400, 150, 250)) + + # Calculate expected NDVI: (nir - red) / (nir + red) + red_vals <- c(100, 200, 150, 50) + nir_vals <- c(300, 400, 150, 250) + expected <- (nir_vals - red_vals) / (nir_vals + red_vals) + + result <- ngr:::.ngr_spk_calc("ndvi", a = a, b = b) + result_vals <- as.vector(terra::values(result)) + + testthat::expect_equal(result_vals, expected, tolerance = 1e-10) +}) + +testthat::test_that(".ngr_spk_calc returns values in expected NDVI range", { + testthat::skip_if_not_installed("terra") + + # NDVI should be between -1 and 1 for valid reflectance values + set.seed(42) + a <- terra::rast(matrix(stats::runif(100, 0, 1000), nrow = 10)) + b <- terra::rast(matrix(stats::runif(100, 0, 1000), nrow = 10)) + + result <- ngr:::.ngr_spk_calc("ndvi", a = a, b = b) + result_vals <- terra::values(result) + + testthat::expect_true(all(result_vals >= -1 & result_vals <= 1, na.rm = TRUE)) +}) + +testthat::test_that(".ngr_spk_calc errors on unknown calc type", { + testthat::skip_if_not_installed("terra") + + a <- terra::rast(matrix(1:4, nrow = 2)) + b <- terra::rast(matrix(5:8, nrow = 2)) + + testthat::expect_error( + ngr:::.ngr_spk_calc("unknown_index", a = a, b = b), + "Unknown calc" + ) +}) + +# Test null coalescing helper -------------------------------------------------- + +testthat::test_that("%||% returns first value when not NULL", { + testthat::expect_equal(ngr:::`%||%`("value", "default"), "value") + testthat::expect_equal(ngr:::`%||%`(0, "default"), 0) + testthat::expect_equal(ngr:::`%||%`(FALSE, TRUE), FALSE) +}) + +testthat::test_that("%||% returns second value when first is NULL", { + testthat::expect_equal(ngr:::`%||%`(NULL, "default"), "default") + testthat::expect_equal(ngr:::`%||%`(NULL, 42), 42) +}) diff --git a/tests/testthat/test-ngr_str_extract_between.R b/tests/testthat/test-ngr_str_extract_between.R new file mode 100644 index 0000000..1c3d974 --- /dev/null +++ b/tests/testthat/test-ngr_str_extract_between.R @@ -0,0 +1,23 @@ +# tests/testthat/test-ngr_str_extract_between.R + +# Test data (reusable across tests) +x <- c( + "Grant Amount: $400,000 Intake Year: 2025", + "Grant Amount: $150,500 Intake Year: 2024" +) + +# Colon required when explicitly included in start pattern +# Expect dollars without the leading colon/space + +testthat::test_that("extracts after colon", { + out <- ngr::ngr_str_extract_between(x, reg_start = "Amount:", reg_end = "Intake") + testthat::expect_equal(out, c("$400,000", "$150,500")) +}) + +# When start pattern is 'Amount' (no colon), current behavior keeps the literal ': ' +# before the captured value in these examples + +testthat::test_that("extracts including leading colon", { + out <- ngr::ngr_str_extract_between(x, reg_start = "Amount", reg_end = "Intake") + testthat::expect_equal(out, c(": $400,000", ": $150,500")) +}) diff --git a/vignettes/.gitignore b/vignettes/.gitignore new file mode 100644 index 0000000..097b241 --- /dev/null +++ b/vignettes/.gitignore @@ -0,0 +1,2 @@ +*.html +*.R diff --git a/vignettes/stac-spectral-indices.Rmd b/vignettes/stac-spectral-indices.Rmd new file mode 100644 index 0000000..2cd970e --- /dev/null +++ b/vignettes/stac-spectral-indices.Rmd @@ -0,0 +1,157 @@ +--- +title: "STAC spectral indices with ngr_spk_stac_calc" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{STAC spectral indices with ngr_spk_stac_calc} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + echo = FALSE, + collapse = TRUE, + comment = "#>" +) +``` + +```{r setup} +library(ngr) +``` + +```{r rstac-query-2000} +# Define an AOI from a bounding box (WGS84) +bbox <- c( + xmin = -126.55350240037997, + ymin = 54.4430453753869, + xmax = -126.52422763064457, + ymax = 54.46001902038006 +) + +aoi <- sf::st_as_sfc(sf::st_bbox(bbox, crs = 4326)) |> + sf::st_as_sf() + +stac_url <- "https://planetarycomputer.microsoft.com/api/stac/v1" +y <- 2000 +date_time <- paste0(y, "-05-01/", y, "-09-15") + +stac_query <- rstac::stac(stac_url) |> + rstac::stac_search( + collections = "landsat-c2-l2", + datetime = date_time, + intersects = sf::st_geometry(aoi)[[1]], + limit = 200 + ) |> + rstac::ext_filter(`eo:cloud_cover` <= 20) + +items <- stac_query |> + rstac::post_request() |> + rstac::items_fetch() |> + rstac::items_sign_planetary_computer() + +ndvi_list <- items$features |> + purrr::map(ngr_spk_stac_calc, aoi = aoi, timing = TRUE) |> + purrr::set_names(purrr::map_chr(items$features, "id")) + +mv <- purrr::imap( + ndvi_list, + function(x, nm) { + rng <- terra::minmax(x) + at <- seq( + floor(rng[1] * 10) / 10, + ceiling(rng[2] * 10) / 10, + by = 0.1 + ) + + mapview::mapview( + x, + layer.name = nm, + at = at, + col.regions = grDevices::hcl.colors(length(at) - 1L, "RdYlGn") + ) + } +) + + +``` + + +```{r map1, out.width="100%"} +# combine into one map with toggleable layers +purrr::reduce(mv, `+`) +``` + + +Create NDVI for multiple years - selecting the raster with the highest likely shrub/tree values. + +```{r multi-year, eval=TRUE} +#----------------------------------------------------------------------------------------------------- +years <- seq(2000, 2025, by = 5) + +ndvi_by_year <- purrr::set_names(years) |> + purrr::map(function(y) { + date_time <- paste0(y, "-06-01/", y, "-07-15") + + items <- rstac::stac(stac_url) |> + rstac::stac_search( + collections = "landsat-c2-l2", + datetime = date_time, + intersects = sf::st_geometry(aoi)[[1]], + limit = 200 + ) |> + rstac::ext_filter(`eo:cloud_cover` <= 50) |> + rstac::post_request() |> + rstac::items_fetch() |> + rstac::items_sign_planetary_computer() + + items$features |> + purrr::map(ngr_spk_stac_calc, aoi = aoi, timing = TRUE) |> + purrr::set_names(purrr::map_chr(items$features, "id")) + }) + +# ndvi_by_year[["2000"]] is your old ndvi_list for 2000 +ndvi_best_by_year <- ndvi_by_year |> + purrr::map(function(ndvi_list) { + + # Compute the proportion of pixels with NDVI >= 0.4 for each raster + prop_high <- sapply( + ndvi_list, + \(r) terra::global(r >= 0.4, mean, na.rm = TRUE)[[1]] + ) + + ndvi_list[which.max(prop_high)] + }) |> + purrr::keep(~ length(.x) > 0) + +ndvi_flat <- purrr::imap(ndvi_best_by_year, \(lst, yr) { + purrr::set_names(lst, paste0(yr, "__", names(lst))) +}) |> + purrr::flatten() + +mv <- purrr::imap( + ndvi_flat, + function(x, nm) { + rng <- terra::minmax(x) + at <- seq( + floor(rng[1] * 10) / 10, + ceiling(rng[2] * 10) / 10, + by = 0.1 + ) + + mapview::mapview( + x, + layer.name = nm, + at = at, + col.regions = grDevices::hcl.colors(length(at) - 1L, "RdYlGn") + ) + } +) + + +``` + +```{r map2, out.width="100%", eval = TRUE} +# combine into one map with toggleable layers +purrr::reduce(mv, `+`) +``` +