From 2ccbdfafe69310f9fd31d4f611cdc92c1fd48d57 Mon Sep 17 00:00:00 2001 From: aitorvv96 Date: Tue, 30 Jun 2026 12:41:19 +0200 Subject: [PATCH 1/6] feat: enhance SDI calculations with beta parameter and add SDI density classification --- NAMESPACE | 1 + R/metrics-stand-density.R | 121 ++++++++++++++++++---------- man/silv_density_sdi.Rd | 30 +++---- man/silv_density_sdi_class.Rd | 43 ++++++++++ tests/testthat/test-stand-density.R | 55 ++++++++++++- 5 files changed, 192 insertions(+), 58 deletions(-) create mode 100644 man/silv_density_sdi_class.Rd diff --git a/NAMESPACE b/NAMESPACE index fcc5e98..1a1b235 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,6 +16,7 @@ export(silv_biomass) export(silv_density_hart) export(silv_density_ntrees_ha) export(silv_density_sdi) +export(silv_density_sdi_class) export(silv_diametric_class) export(silv_dominant_height) export(silv_lorey_height) diff --git a/R/metrics-stand-density.R b/R/metrics-stand-density.R index 6af2fb1..aa5c358 100644 --- a/R/metrics-stand-density.R +++ b/R/metrics-stand-density.R @@ -54,68 +54,105 @@ silv_density_ntrees_ha <- function(ntrees, #' Calculates the Stand Density Index #' -#' The Stand Density Index (SDI) is relationship between the average tree size and +#' The Stand Density Index (SDI) is the relationship between the average tree size and #' density of trees per hectare. #' #' @template ntrees #' @template dg -#' @param classify whether to classify the values using USDA thresholds -#' @param max_sdi used when \code{classify = TRUE}. The maximum SDi, which depends -#' on the species, stand type, and site -#' -#' @return A numeric vector +#' @param beta The Stand Density Index exponent (default is \code{1.605}). +#' +#' @return A numeric vector representing the absolute SDI. #' @export -#' +#' #' @details -#' The SDI has different interpretation depending on the species, location, and also +#' The SDI has different interpretations depending on the species, location, and also #' the management type (even-aged, uneven-aged...). The value of maximum SDI must -#' be determined from the literature and used carefully. The option \code{classify = TRUE} -#' will use this value to classify the SDI in low density (<24%), moderate density (24-35%), -#' high density (34-55%), and extremely high density (>55%). +#' be determined from the literature and used carefully. The \code{beta} exponent allows +#' adjustments for different species or mixed stands. +#' +#' @references Reineke, L. H. (1933). Perfecting a stand-density index for even-aged forests. +#' Journal of Agricultural Research, 46(7), 627-638. #' #' @examples -#' ## calculate SDI for a Pinus sulvestris stand (max 990) -#' silv_density_sdi(ntrees = 800, dg = 23.4, max_sdi = 990) -#' -#' ## check base classification (other can be used) -#' silv_density_sdi(ntrees = 800, dg = 23.4, classify = TRUE, max_sdi = 990) +#' ## calculate SDI for a Pinus sylvestris stand (beta = 1.605) +#' silv_density_sdi(ntrees = 800, dg = 23.4) +#' +#' ## calculate SDI with custom beta +#' silv_density_sdi(ntrees = 800, dg = 23.4, beta = 1.7) silv_density_sdi <- function( - ntrees, - dg, - classify = FALSE, - max_sdi = NULL + ntrees, + dg, + beta = 1.605 ) { - - # 0. Validate inputs + # 0. validate inputs assert_positive_numeric(ntrees, "ntrees") assert_positive_numeric(dg, "dg") - assert_logical(classify, "classify") + if (!is.numeric(beta)) cli::cli_abort("{.arg beta} has to be a numeric vector.") assert_same_length(ntrees, dg, names = c("ntrees", "dg")) + # 1. calculate sdi + sdi <- ntrees * ((25.4 / dg) ** -abs(beta)) # note: abs() avoids errors with signs + return(sdi) +} - # 1. Calculate SDI - sdi <- ntrees * ((dg / 25.4)) ** 1.605 - # 2. Classify? - if (classify) { +#' Classifies the Stand Density Index +#' +#' Classifies the Stand Density Index (SDI) into density classes or calculates the relative SDI +#' percentage based on USDA thresholds. +#' +#' @param sdi A numeric vector representing the Stand Density Index. +#' @param max_sdi A numeric vector representing the maximum SDI for the species/site. +#' @param classify A logical value indicating whether to classify the values into density classes +#' (default is \code{TRUE}). If \code{FALSE}, it returns the relative SDI as a percentage. +#' +#' @return A character vector with the density classes if \code{classify = TRUE}, or a numeric vector +#' with the relative SDI percentage if \code{classify = FALSE}. +#' @export +#' +#' @details +#' The option \code{classify = TRUE} will use the \code{max_sdi} value to classify the SDI into +#' four competitive and growth conditions: low density (<24%), moderate density (24-35%), +#' high density (34-55%), and extremely high density (>55%). +#' +#' @references USDA Forest Service. (n.d.). Stand Density Index. +#' https://www.fs.usda.gov/Internet/FSE_DOCUMENTS/stelprdb5270993.pdf +#' +#' @examples +#' ## calculate SDI for a Pinus sylvestris stand (max 990) +#' sdi_val <- silv_density_sdi(ntrees = 800, dg = 23.4) +#' +#' ## check base classification +#' silv_density_sdi_class(sdi = sdi_val, max_sdi = 990) +#' +#' ## get relative SDI percentage +#' silv_density_sdi_class(sdi = sdi_val, max_sdi = 990, classify = FALSE) +silv_density_sdi_class <- function( + sdi, + max_sdi, + classify = TRUE +) { + # 0. validate inputs + assert_positive_numeric(sdi, "sdi") + assert_positive_numeric(max_sdi, "max_sdi") + assert_logical(classify, "classify") + assert_same_length(sdi, max_sdi, names = c("sdi", "max_sdi")) - ## assert inputs - if (is.null(max_sdi)) cli::cli_abort("You must specify when ") - assert_positive_numeric(max_sdi, "max_sdi") - - ## calculate - sdi <- (sdi / max_sdi) * 100 - sdi <- dplyr::case_when( - sdi <= 24 ~ "Low density", - sdi > 24 & sdi <= 34 ~ "Moderate density", - sdi > 34 & sdi <= 55 ~ "High density", - sdi > 55 ~ "Extremely high density" + # 1. calculate relative sdi + rel_sdi <- (sdi / max_sdi) * 100 + + # 2. classify or return percentage + if (classify) { + res <- dplyr::case_when( + rel_sdi <= 24 ~ "Low density", + rel_sdi > 24 & rel_sdi <= 34 ~ "Moderate density", + rel_sdi > 34 & rel_sdi <= 55 ~ "High density", + rel_sdi > 55 ~ "Extremely high density" ) - } else if (!is.null(max_sdi)) { - sdi <- (sdi / max_sdi) * 100 + } else { + res <- rel_sdi } - return(sdi) - + return(res) } diff --git a/man/silv_density_sdi.Rd b/man/silv_density_sdi.Rd index 422bc24..2c77973 100644 --- a/man/silv_density_sdi.Rd +++ b/man/silv_density_sdi.Rd @@ -4,7 +4,7 @@ \alias{silv_density_sdi} \title{Calculates the Stand Density Index} \usage{ -silv_density_sdi(ntrees, dg, classify = FALSE, max_sdi = NULL) +silv_density_sdi(ntrees, dg, beta = 1.605) } \arguments{ \item{ntrees}{Numeric vector with number of trees of the diameter class per @@ -13,29 +13,29 @@ corresponds to only one tree} \item{dg}{Numeric vector of quadratic mean diameters} -\item{classify}{whether to classify the values using USDA thresholds} - -\item{max_sdi}{used when \code{classify = TRUE}. The maximum SDi, which depends -on the species, stand type, and site} +\item{beta}{The Stand Density Index exponent (default is \code{1.605}).} } \value{ -A numeric vector +A numeric vector representing the absolute SDI. } \description{ -The Stand Density Index (SDI) is relationship between the average tree size and +The Stand Density Index (SDI) is the relationship between the average tree size and density of trees per hectare. } \details{ -The SDI has different interpretation depending on the species, location, and also +The SDI has different interpretations depending on the species, location, and also the management type (even-aged, uneven-aged...). The value of maximum SDI must -be determined from the literature and used carefully. The option \code{classify = TRUE} -will use this value to classify the SDI in low density (<24\%), moderate density (24-35\%), -high density (34-55\%), and extremely high density (>55\%). +be determined from the literature and used carefully. The \code{beta} exponent allows +adjustments for different species or mixed stands. } \examples{ -## calculate SDI for a Pinus sulvestris stand (max 990) -silv_density_sdi(ntrees = 800, dg = 23.4, max_sdi = 990) +## calculate SDI for a Pinus sylvestris stand (beta = 1.605) +silv_density_sdi(ntrees = 800, dg = 23.4) -## check base classification (other can be used) -silv_density_sdi(ntrees = 800, dg = 23.4, classify = TRUE, max_sdi = 990) +## calculate SDI with custom beta +silv_density_sdi(ntrees = 800, dg = 23.4, beta = 1.7) +} +\references{ +Reineke, L. H. (1933). Perfecting a stand-density index for even-aged forests. +Journal of Agricultural Research, 46(7), 627-638. } diff --git a/man/silv_density_sdi_class.Rd b/man/silv_density_sdi_class.Rd new file mode 100644 index 0000000..f92a509 --- /dev/null +++ b/man/silv_density_sdi_class.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/metrics-stand-density.R +\name{silv_density_sdi_class} +\alias{silv_density_sdi_class} +\title{Classifies the Stand Density Index} +\usage{ +silv_density_sdi_class(sdi, max_sdi, classify = TRUE) +} +\arguments{ +\item{sdi}{A numeric vector representing the Stand Density Index.} + +\item{max_sdi}{A numeric vector representing the maximum SDI for the species/site.} + +\item{classify}{A logical value indicating whether to classify the values into density classes +(default is \code{TRUE}). If \code{FALSE}, it returns the relative SDI as a percentage.} +} +\value{ +A character vector with the density classes if \code{classify = TRUE}, or a numeric vector +with the relative SDI percentage if \code{classify = FALSE}. +} +\description{ +Classifies the Stand Density Index (SDI) into density classes or calculates the relative SDI +percentage based on USDA thresholds. +} +\details{ +The option \code{classify = TRUE} will use the \code{max_sdi} value to classify the SDI into +four competitive and growth conditions: low density (<24\%), moderate density (24-35\%), +high density (34-55\%), and extremely high density (>55\%). +} +\examples{ +## calculate SDI for a Pinus sylvestris stand (max 990) +sdi_val <- silv_density_sdi(ntrees = 800, dg = 23.4) + +## check base classification +silv_density_sdi_class(sdi = sdi_val, max_sdi = 990) + +## get relative SDI percentage +silv_density_sdi_class(sdi = sdi_val, max_sdi = 990, classify = FALSE) +} +\references{ +USDA Forest Service. (n.d.). Stand Density Index. +https://www.fs.usda.gov/Internet/FSE_DOCUMENTS/stelprdb5270993.pdf +} diff --git a/tests/testthat/test-stand-density.R b/tests/testthat/test-stand-density.R index c54fd2c..75e4952 100644 --- a/tests/testthat/test-stand-density.R +++ b/tests/testthat/test-stand-density.R @@ -77,5 +77,58 @@ test_that("Errors work", { ) expect_error(silv_density_hart("17.8", 400)) expect_error(silv_density_hart(17.8, "400")) - expect_error(silv_density_hart(c(17.8, 20.5), 400)) + expect_error( + silv_density_hart(c(17.8, 20.5), 400) + ) +}) + + +# 3. silv_density_sdi ---------------------------------------------------- + +## Tests +test_that("Stand Density Index is well calculated", { + # default beta (1.605) + expect_equal( + silv_density_sdi(ntrees = 800, dg = 23.4), + 702.40, + tolerance = 0.01 + ) + + # custom beta + expect_equal( + silv_density_sdi(ntrees = 800, dg = 23.4, beta = 1.7), + 692.68, + tolerance = 0.01 + ) + + # negative beta + expect_equal( + silv_density_sdi(ntrees = 800, dg = 23.4, beta = -1.605), + 702.40, + tolerance = 0.01 + ) + + # with max_sdi (returns percentage) using silv_density_sdi_class + sdi_val <- silv_density_sdi(ntrees = 800, dg = 23.4) + expect_equal( + silv_density_sdi_class(sdi = sdi_val, max_sdi = 990, classify = FALSE), + 70.95, + tolerance = 0.01 + ) + + # with classification using silv_density_sdi_class + expect_equal( + silv_density_sdi_class(sdi = sdi_val, max_sdi = 990), + "Extremely high density" + ) +}) + +## Test errors +test_that("Errors work in silv_density_sdi", { + expect_error(silv_density_sdi(800, 23.4, beta = "1.605")) +}) + +test_that("Errors work in silv_density_sdi_class", { + expect_error(silv_density_sdi_class(700, "990")) + expect_error(silv_density_sdi_class(700, 990, classify = "TRUE")) }) From 2f8e14c7ef89da2e4468909d9cea8e8153f706f3 Mon Sep 17 00:00:00 2001 From: aitorvv96 Date: Tue, 30 Jun 2026 12:43:45 +0200 Subject: [PATCH 2/6] docs: add direct PDF download URL for Reineke (1933) paper --- R/metrics-stand-density.R | 28 +++++++--------------------- man/silv_density_ntrees_ha.Rd | 2 +- man/silv_density_sdi.Rd | 2 +- 3 files changed, 9 insertions(+), 23 deletions(-) diff --git a/R/metrics-stand-density.R b/R/metrics-stand-density.R index aa5c358..da792bf 100644 --- a/R/metrics-stand-density.R +++ b/R/metrics-stand-density.R @@ -1,5 +1,3 @@ - - #' Calculates number of trees per hectare #' #' Calculates number of trees per hectare for a given plot size and shape @@ -28,12 +26,11 @@ #' n, #' plot_size = c(10, 15), #' plot_shape = "rectangular" -#' ) +#' ) #' ) silv_density_ntrees_ha <- function(ntrees, - plot_size, - plot_shape = "circular") { - + plot_size, + plot_shape = "circular") { # 0. Handle errors stopifnot(plot_shape %in% c("circular", "rectangular")) if (length(plot_size) == 1 && plot_size <= 0) cli::cli_abort("`plot_size` has to be greater than 0") @@ -44,8 +41,6 @@ silv_density_ntrees_ha <- function(ntrees, } else { ntrees * 10000 / prod(plot_size) } - - } @@ -71,7 +66,7 @@ silv_density_ntrees_ha <- function(ntrees, #' adjustments for different species or mixed stands. #' #' @references Reineke, L. H. (1933). Perfecting a stand-density index for even-aged forests. -#' Journal of Agricultural Research, 46(7), 627-638. +#' Journal of Agricultural Research, 46(7), 627-638. URL: https://research.fs.usda.gov/download/treesearch/60134.pdf #' #' @examples #' ## calculate SDI for a Pinus sylvestris stand (beta = 1.605) @@ -91,7 +86,7 @@ silv_density_sdi <- function( assert_same_length(ntrees, dg, names = c("ntrees", "dg")) # 1. calculate sdi - sdi <- ntrees * ((25.4 / dg) ** -abs(beta)) # note: abs() avoids errors with signs + sdi <- ntrees * ((25.4 / dg)**-abs(beta)) # note: abs() avoids errors with signs return(sdi) } @@ -111,8 +106,8 @@ silv_density_sdi <- function( #' @export #' #' @details -#' The option \code{classify = TRUE} will use the \code{max_sdi} value to classify the SDI into -#' four competitive and growth conditions: low density (<24%), moderate density (24-35%), +#' The option \code{classify = TRUE} will use the \code{max_sdi} value to classify the SDI into +#' four competitive and growth conditions: low density (<24%), moderate density (24-35%), #' high density (34-55%), and extremely high density (>55%). #' #' @references USDA Forest Service. (n.d.). Stand Density Index. @@ -156,9 +151,6 @@ silv_density_sdi_class <- function( } - - - #' Hart or Hart-Becking spacing index #' #' Calculates the Hart Index or the Hart-Becking Index for even-aged stands @@ -200,7 +192,6 @@ silv_density_hart <- function( ntrees, which = c("hart", "hart-becking") ) { - # 0. Validate inputs assert_positive_numeric(h0, "h0") assert_positive_numeric(ntrees, "ntrees") @@ -213,9 +204,4 @@ silv_density_hart <- function( "hart-becking" = sqrt(20000 / (ntrees * sqrt(3))) / h0 * 100, cli::cli_abort("`which` must be either or ") ) - } - - - - diff --git a/man/silv_density_ntrees_ha.Rd b/man/silv_density_ntrees_ha.Rd index a9361c9..35dfd4b 100644 --- a/man/silv_density_ntrees_ha.Rd +++ b/man/silv_density_ntrees_ha.Rd @@ -37,6 +37,6 @@ inventory_samples |> n, plot_size = c(10, 15), plot_shape = "rectangular" - ) + ) ) } diff --git a/man/silv_density_sdi.Rd b/man/silv_density_sdi.Rd index 2c77973..6800fb2 100644 --- a/man/silv_density_sdi.Rd +++ b/man/silv_density_sdi.Rd @@ -37,5 +37,5 @@ silv_density_sdi(ntrees = 800, dg = 23.4, beta = 1.7) } \references{ Reineke, L. H. (1933). Perfecting a stand-density index for even-aged forests. -Journal of Agricultural Research, 46(7), 627-638. +Journal of Agricultural Research, 46(7), 627-638. URL: https://research.fs.usda.gov/download/treesearch/60134.pdf } From b30f8fc34a0c3fa3fac8f6409a7a59ea414e7681 Mon Sep 17 00:00:00 2001 From: aitorvv96 Date: Tue, 30 Jun 2026 12:51:29 +0200 Subject: [PATCH 3/6] feat: add SDI beta autoselector and sdi_models dataset --- NAMESPACE | 1 + R/data.R | 18 +++ R/metrics-stand-density.R | 203 ++++++++++++++++++++++++++++ data/sdi_models.rda | Bin 0 -> 1065 bytes man/sdi_models.Rd | 26 ++++ man/silv_density_sdi_auto.Rd | 56 ++++++++ tests/testthat/test-stand-density.R | 25 ++++ 7 files changed, 329 insertions(+) create mode 100644 data/sdi_models.rda create mode 100644 man/sdi_models.Rd create mode 100644 man/silv_density_sdi_auto.Rd diff --git a/NAMESPACE b/NAMESPACE index 1a1b235..6e61043 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,6 +16,7 @@ export(silv_biomass) export(silv_density_hart) export(silv_density_ntrees_ha) export(silv_density_sdi) +export(silv_density_sdi_auto) export(silv_density_sdi_class) export(silv_diametric_class) export(silv_dominant_height) diff --git a/R/data.R b/R/data.R index b739df3..e24ee5c 100644 --- a/R/data.R +++ b/R/data.R @@ -107,3 +107,21 @@ #' MITECO. 4th Spanish National Forest Inventory - SIG database codes. #' \url{https://www.miteco.gob.es/content/dam/miteco/es/biodiversidad/temas/inventarios-nacionales/documentador_sig_tcm30-536622.pdf} "snfi4_volume_coefficients" + + +#' SDI beta models +#' +#' Specific beta coefficients for Reineke's Stand Density Index (SDI) per +#' species and region. +#' +#' @format A `tibble` +#' \describe{ +#' \item{article_id}{Character. Short identifier of the source article.} +#' \item{title}{Character. Full title of the source article.} +#' \item{doi_url}{Character. DOI URL of the source article.} +#' \item{country}{Character. Country where the study was conducted.} +#' \item{region}{Character. Region within the country.} +#' \item{species}{Character. Scientific name of the tree species.} +#' \item{beta}{Numeric. Beta coefficient for SDI calculation.} +#' } +"sdi_models" diff --git a/R/metrics-stand-density.R b/R/metrics-stand-density.R index da792bf..b9832da 100644 --- a/R/metrics-stand-density.R +++ b/R/metrics-stand-density.R @@ -90,6 +90,209 @@ silv_density_sdi <- function( return(sdi) } +# --- Internal Auto Selector Helper --- + +#' @noRd +.auto_select_sdi_beta <- function(species, region = NULL) { + + sdi_models <- silviculture::sdi_models + + # 1. Try exact species + region match + if (!is.null(region)) { + sel <- sdi_models[sdi_models$species == species & sdi_models$region == region, ] + if (nrow(sel) > 0) { + return(list( + model = sel$article_id[1], + beta = sel$beta[1], + matched_species = species, + matched_region = region, + is_fallback = FALSE, + fallback_type = NA_character_, + model_desc = paste0(sel$article_id[1], " (", region, ")") + )) + } + } + + # 2. Try species + "all" regions fallback + sel <- sdi_models[sdi_models$species == species & sdi_models$region == "all", ] + if (nrow(sel) > 0) { + return(list( + model = sel$article_id[1], + beta = sel$beta[1], + matched_species = species, + matched_region = "all", + is_fallback = !is.null(region), + fallback_type = "region", + model_desc = paste0(sel$article_id[1], " (all regions)") + )) + } + + # 3. Try genus fallback (genus spp.) + genus <- strsplit(species, " ")[[1]][1] + genus_spp <- paste0(genus, " spp.") + + if (!is.null(region)) { + sel <- sdi_models[sdi_models$species == genus_spp & sdi_models$region == region, ] + if (nrow(sel) > 0) { + return(list( + model = sel$article_id[1], + beta = sel$beta[1], + matched_species = genus_spp, + matched_region = region, + is_fallback = TRUE, + fallback_type = "genus", + model_desc = paste0(sel$article_id[1], " (genus fallback)") + )) + } + } + + sel <- sdi_models[sdi_models$species == genus_spp & sdi_models$region == "all", ] + if (nrow(sel) > 0) { + return(list( + model = sel$article_id[1], + beta = sel$beta[1], + matched_species = genus_spp, + matched_region = "all", + is_fallback = TRUE, + fallback_type = "genus", + model_desc = paste0(sel$article_id[1], " (genus fallback)") + )) + } + + # 4. Total fallback (default) + sel <- sdi_models[sdi_models$species == "default", ] + if (nrow(sel) == 0) { + # Absolute safety fallback in case dataset is malformed + default_beta <- 1.605 + } else { + default_beta <- sel$beta[1] + } + + return(list( + model = "default", + beta = default_beta, + matched_species = "default", + matched_region = "default", + is_fallback = TRUE, + fallback_type = "default", + model_desc = paste0("default (", default_beta, ")") + )) +} + +#' Predict Stand Density Index automatically +#' +#' @description +#' `silv_density_sdi_auto()` is a vectorized function that automatically selects +#' the best available Stand Density Index exponent (\code{beta}) for each row +#' based on a provided species and region from the internal \code{sdi_models} database. +#' +#' If an exact species and region match is not found, the function falls back to a +#' country-wide species model (\code{region = "all"}), then to a genus-level fallback +#' (e.g., "Pinus spp."), and finally to the default SDI exponent (\code{beta = 1.605}). +#' +#' @template ntrees +#' @template dg +#' @param species A character string or vector of tree species (e.g., `"Pinus sylvestris"`). +#' @param region A character string or vector of the region (e.g., `"Castilla y Leon"`). +#' Defaults to `NULL` (no region specified). +#' @param quiet Logical. If `FALSE`, informs the user about fallbacks to genus or default models. +#' +#' @return A `data.frame` with two columns: +#' - `sdi`: The computed absolute Stand Density Index. +#' - `sdi_model`: The model used (e.g., `"del-rio-2006 (Castilla y Leon)"`, +#' `"aguirre-2017 (genus fallback)"`, or `"default (1.605)"`). +#' +#' @name silv_density_sdi_auto +#' +#' @examples +#' # Calculate SDI with automatic selection +#' silv_density_sdi_auto( +#' ntrees = 800, +#' dg = 23.4, +#' species = "Pinus sylvestris", +#' region = "Castilla y Leon" +#' ) +#' +#' # Fallback to default +#' silv_density_sdi_auto( +#' ntrees = 800, +#' dg = 23.4, +#' species = "Unknown species" +#' ) +#' +#' @export +silv_density_sdi_auto <- function( + ntrees, + dg, + species, + region = NULL, + quiet = FALSE +) { + # Validations + n_trees <- length(ntrees) + assert_positive_numeric(ntrees, "ntrees") + assert_positive_numeric(dg, "dg") + assert_same_length(ntrees, dg, names = c("ntrees", "dg")) + + if (length(species) == 1) { + species <- rep(species, n_trees) + } else if (length(species) != n_trees) { + cli::cli_abort("{.arg species} must be of length 1 or the same length as {.arg ntrees}.") + } + + if (is.null(region)) { + region <- rep(NA_character_, n_trees) + } else if (length(region) == 1) { + region <- rep(region, n_trees) + } else if (length(region) != n_trees) { + cli::cli_abort("{.arg region} must be of length 1 or the same length as {.arg ntrees}.") + } + + sdi_values <- rep(NA_real_, n_trees) + sdi_models_used <- rep(NA_character_, n_trees) + + unique_combos <- unique(data.frame( + species = species, + region = region, + stringsAsFactors = FALSE + )) + + for (i in seq_len(nrow(unique_combos))) { + sp <- unique_combos$species[i] + reg <- unique_combos$region[i] + + # Handle NA in region + if (is.na(reg)) { + idx <- which(species == sp & is.na(region)) + reg_arg <- NULL + } else { + idx <- which(species == sp & region == reg) + reg_arg <- reg + } + + best_model_info <- .auto_select_sdi_beta(sp, reg_arg) + + if (best_model_info$is_fallback && !quiet) { + if (best_model_info$fallback_type == "default") { + cli::cli_alert_info("Exact model for {.val {sp}} not found. Using default beta {.val {best_model_info$beta}}.") + } else if (best_model_info$fallback_type == "genus") { + cli::cli_alert_info("Exact species {.val {sp}} not found. Using genus fallback {.val {best_model_info$matched_species}} from {.val {best_model_info$model}}.") + } else if (best_model_info$fallback_type == "region") { + cli::cli_alert_info("Exact region {.val {reg_arg}} not found for {.val {sp}}. Using fallback region {.val {best_model_info$matched_region}} from {.val {best_model_info$model}}.") + } + } + + sdi_values[idx] <- silv_density_sdi(ntrees[idx], dg[idx], beta = best_model_info$beta) + sdi_models_used[idx] <- best_model_info$model_desc + } + + return(data.frame( + sdi = sdi_values, + sdi_model = sdi_models_used, + stringsAsFactors = FALSE + )) +} + #' Classifies the Stand Density Index #' diff --git a/data/sdi_models.rda b/data/sdi_models.rda new file mode 100644 index 0000000000000000000000000000000000000000..63497e55069409d867c9ac91bf94e448de148f92 GIT binary patch literal 1065 zcmV+^1lIdPT4*^jL0KkKS$IHYh5!O&fB*l#O+-K^`_L`^e?Y(I|L{;0NCG1>8bHW^ z5CJd<;kVEPo&dZVMF@)_B`K)T&;jZI0g<3IGypw700006jWiENru9RB0B8Vc13&-< zj6eX;00000&;x(~XaHyfKmZ1eKmgDH0002c0ThTTc#Njg6w}HX27mwn0gwYA4H$-i z002c11e+w&Q!uIRNwkkpVF92Fk5Jl;0B8n)00E(*&+=*+F{bVn0Hh{*$S^SHg*y@0 zf=*sOPG?`8a;K%-?CznK2a~tM)2JwWS<)~Is;R>ZW?-&-wF3SVm$SKMO@pGICHxl7 zS4envnCs5W=MIhqbOv|@=ZdhZAB=>Qf?yOT1w&&1LJ*yJ9AG~9Q?&_}1^;jRy1mcn z(DNMu2|8V)M2CgyG4G|#yFirZBvBfr>BJOEfVn18#DOeIB+CgSYXd0&WR_N=b8r1G zTZ{6x*7W4Ty2Qcd4kqNaEbUS75?(+rbEj1xOn^d5=!c-d5Dz27g<5gSOIM)FYIWp8 zMI2;A15I%zCMS%c#K#uYz5W<-U*MlWi2`6*7#(I$k2-|6&m{E7Wim}jb!DN{>mlNg zVp-17!KB5;)csNfS?o50@)4MFWO;B(xo=3i!`;9XdSTv z-U6)w%6>_D#6=*g0DtOKt$YgZj)>@t ji{xX8&{MY{iRE(JNE`Chn#)oqOa3n8ig2MI@PNw Date: Tue, 30 Jun 2026 13:04:10 +0200 Subject: [PATCH 4/6] feat: rename sdi_models to sdi_coefficients and add country param and beta column --- R/data.R | 6 +- R/metrics-stand-density.R | 184 ++++++++++++++++----- data/sdi_coefficients.rda | Bin 0 -> 1160 bytes data/sdi_models.rda | Bin 1065 -> 0 bytes man/{sdi_models.Rd => sdi_coefficients.Rd} | 10 +- man/silv_density_sdi_auto.Rd | 32 ++-- tests/testthat/test-stand-density.R | 13 +- 7 files changed, 184 insertions(+), 61 deletions(-) create mode 100644 data/sdi_coefficients.rda delete mode 100644 data/sdi_models.rda rename man/{sdi_models.Rd => sdi_coefficients.Rd} (84%) diff --git a/R/data.R b/R/data.R index e24ee5c..f8098c8 100644 --- a/R/data.R +++ b/R/data.R @@ -109,10 +109,10 @@ "snfi4_volume_coefficients" -#' SDI beta models +#' SDI beta coefficients #' #' Specific beta coefficients for Reineke's Stand Density Index (SDI) per -#' species and region. +#' species, country, and region. #' #' @format A `tibble` #' \describe{ @@ -124,4 +124,4 @@ #' \item{species}{Character. Scientific name of the tree species.} #' \item{beta}{Numeric. Beta coefficient for SDI calculation.} #' } -"sdi_models" +"sdi_coefficients" diff --git a/R/metrics-stand-density.R b/R/metrics-stand-density.R index b9832da..28e3856 100644 --- a/R/metrics-stand-density.R +++ b/R/metrics-stand-density.R @@ -93,89 +93,166 @@ silv_density_sdi <- function( # --- Internal Auto Selector Helper --- #' @noRd -.auto_select_sdi_beta <- function(species, region = NULL) { +.auto_select_sdi_beta <- function(species, country = NULL, region = NULL) { - sdi_models <- silviculture::sdi_models + sdi_coefficients <- silviculture::sdi_coefficients - # 1. Try exact species + region match - if (!is.null(region)) { - sel <- sdi_models[sdi_models$species == species & sdi_models$region == region, ] + # 1. Try exact species + country + region match + if (!is.null(country) && !is.null(region)) { + sel <- sdi_coefficients[sdi_coefficients$species == species & + sdi_coefficients$country == country & + sdi_coefficients$region == region, ] if (nrow(sel) > 0) { return(list( model = sel$article_id[1], beta = sel$beta[1], matched_species = species, + matched_country = country, matched_region = region, is_fallback = FALSE, fallback_type = NA_character_, - model_desc = paste0(sel$article_id[1], " (", region, ")") + model_desc = paste0(sel$article_id[1], " (", country, ", ", region, ")") )) } } - # 2. Try species + "all" regions fallback - sel <- sdi_models[sdi_models$species == species & sdi_models$region == "all", ] + # 2. Try species + country + "all" regions match + if (!is.null(country)) { + sel <- sdi_coefficients[sdi_coefficients$species == species & + sdi_coefficients$country == country & + sdi_coefficients$region == "all", ] + if (nrow(sel) > 0) { + return(list( + model = sel$article_id[1], + beta = sel$beta[1], + matched_species = species, + matched_country = country, + matched_region = "all", + is_fallback = !is.null(region), + fallback_type = "region", + model_desc = paste0(sel$article_id[1], " (", country, ", all regions)") + )) + } + } + + # 3. Fallback when country is not found/specified, but we find the species in some other country + sel <- sdi_coefficients[sdi_coefficients$species == species, ] if (nrow(sel) > 0) { + if (!is.null(region)) { + sel_reg <- sel[sel$region == region, ] + if (nrow(sel_reg) > 0) { + return(list( + model = sel_reg$article_id[1], + beta = sel_reg$beta[1], + matched_species = species, + matched_country = sel_reg$country[1], + matched_region = region, + is_fallback = TRUE, + fallback_type = "country", + model_desc = paste0(sel_reg$article_id[1], " (", sel_reg$country[1], ", ", region, ")") + )) + } + } + + sel_all <- sel[sel$region == "all", ] + if (nrow(sel_all) > 0) { + return(list( + model = sel_all$article_id[1], + beta = sel_all$beta[1], + matched_species = species, + matched_country = sel_all$country[1], + matched_region = "all", + is_fallback = TRUE, + fallback_type = "region", + model_desc = paste0(sel_all$article_id[1], " (", sel_all$country[1], ", all regions)") + )) + } + return(list( model = sel$article_id[1], beta = sel$beta[1], matched_species = species, - matched_region = "all", - is_fallback = !is.null(region), + matched_country = sel$country[1], + matched_region = sel$region[1], + is_fallback = TRUE, fallback_type = "region", - model_desc = paste0(sel$article_id[1], " (all regions)") + model_desc = paste0(sel$article_id[1], " (", sel$country[1], ", ", sel$region[1], ")") )) } - # 3. Try genus fallback (genus spp.) + # 4. Try genus fallback (genus spp.) genus <- strsplit(species, " ")[[1]][1] genus_spp <- paste0(genus, " spp.") - if (!is.null(region)) { - sel <- sdi_models[sdi_models$species == genus_spp & sdi_models$region == region, ] + if (!is.null(country)) { + if (!is.null(region)) { + sel <- sdi_coefficients[sdi_coefficients$species == genus_spp & + sdi_coefficients$country == country & + sdi_coefficients$region == region, ] + if (nrow(sel) > 0) { + return(list( + model = sel$article_id[1], + beta = sel$beta[1], + matched_species = genus_spp, + matched_country = country, + matched_region = region, + is_fallback = TRUE, + fallback_type = "genus", + model_desc = paste0(sel$article_id[1], " (genus fallback: ", country, ", ", region, ")") + )) + } + } + + sel <- sdi_coefficients[sdi_coefficients$species == genus_spp & + sdi_coefficients$country == country & + sdi_coefficients$region == "all", ] if (nrow(sel) > 0) { return(list( model = sel$article_id[1], beta = sel$beta[1], matched_species = genus_spp, - matched_region = region, + matched_country = country, + matched_region = "all", is_fallback = TRUE, fallback_type = "genus", - model_desc = paste0(sel$article_id[1], " (genus fallback)") + model_desc = paste0(sel$article_id[1], " (genus fallback: ", country, ", all regions)") )) } } - sel <- sdi_models[sdi_models$species == genus_spp & sdi_models$region == "all", ] + sel <- sdi_coefficients[sdi_coefficients$species == genus_spp, ] if (nrow(sel) > 0) { return(list( model = sel$article_id[1], beta = sel$beta[1], matched_species = genus_spp, - matched_region = "all", + matched_country = sel$country[1], + matched_region = sel$region[1], is_fallback = TRUE, fallback_type = "genus", - model_desc = paste0(sel$article_id[1], " (genus fallback)") + model_desc = paste0(sel$article_id[1], " (genus fallback: ", sel$country[1], ", ", sel$region[1], ")") )) } - # 4. Total fallback (default) - sel <- sdi_models[sdi_models$species == "default", ] + # 5. Total fallback (default) + sel <- sdi_coefficients[sdi_coefficients$species == "default", ] if (nrow(sel) == 0) { - # Absolute safety fallback in case dataset is malformed - default_beta <- 1.605 + default_beta <- -1.605 + default_desc <- "reineke-1933 (-1.605)" } else { default_beta <- sel$beta[1] + default_desc <- paste0(sel$article_id[1], " (", default_beta, ")") } return(list( - model = "default", + model = sel$article_id[1], beta = default_beta, matched_species = "default", + matched_country = "default", matched_region = "default", is_fallback = TRUE, fallback_type = "default", - model_desc = paste0("default (", default_beta, ")") + model_desc = default_desc )) } @@ -184,23 +261,27 @@ silv_density_sdi <- function( #' @description #' `silv_density_sdi_auto()` is a vectorized function that automatically selects #' the best available Stand Density Index exponent (\code{beta}) for each row -#' based on a provided species and region from the internal \code{sdi_models} database. +#' based on a provided species, country, and region from the internal \code{sdi_coefficients} database. #' -#' If an exact species and region match is not found, the function falls back to a -#' country-wide species model (\code{region = "all"}), then to a genus-level fallback -#' (e.g., "Pinus spp."), and finally to the default SDI exponent (\code{beta = 1.605}). +#' If an exact species, country, and region match is not found, the function falls back to a +#' country-wide species model (\code{region = "all"}), then searches other countries, then falls +#' back to a genus-level fallback (e.g., "Pinus spp."), and finally to the default SDI exponent +#' (\code{beta = -1.605} from Reineke 1933). #' #' @template ntrees #' @template dg #' @param species A character string or vector of tree species (e.g., `"Pinus sylvestris"`). -#' @param region A character string or vector of the region (e.g., `"Castilla y Leon"`). +#' @param country A character string or vector of the country (e.g., `"Spain"`). +#' Defaults to `NULL` (no country specified). +#' @param region A character string or vector of the region (e.g., `"Castilla y León"`). #' Defaults to `NULL` (no region specified). #' @param quiet Logical. If `FALSE`, informs the user about fallbacks to genus or default models. #' -#' @return A `data.frame` with two columns: +#' @return A `data.frame` with three columns: #' - `sdi`: The computed absolute Stand Density Index. -#' - `sdi_model`: The model used (e.g., `"del-rio-2006 (Castilla y Leon)"`, -#' `"aguirre-2017 (genus fallback)"`, or `"default (1.605)"`). +#' - `beta`: The beta exponent used for the calculation. +#' - `sdi_model`: The model used (e.g., `"del-rio-2006 (Spain, Castilla y León)"`, +#' `"reineke-1933 (-1.605)"`, etc.). #' #' @name silv_density_sdi_auto #' @@ -210,7 +291,7 @@ silv_density_sdi <- function( #' ntrees = 800, #' dg = 23.4, #' species = "Pinus sylvestris", -#' region = "Castilla y Leon" +#' region = "Castilla y León" #' ) #' #' # Fallback to default @@ -225,6 +306,7 @@ silv_density_sdi_auto <- function( ntrees, dg, species, + country = NULL, region = NULL, quiet = FALSE ) { @@ -240,6 +322,14 @@ silv_density_sdi_auto <- function( cli::cli_abort("{.arg species} must be of length 1 or the same length as {.arg ntrees}.") } + if (is.null(country)) { + country <- rep(NA_character_, n_trees) + } else if (length(country) == 1) { + country <- rep(country, n_trees) + } else if (length(country) != n_trees) { + cli::cli_abort("{.arg country} must be of length 1 or the same length as {.arg ntrees}.") + } + if (is.null(region)) { region <- rep(NA_character_, n_trees) } else if (length(region) == 1) { @@ -249,28 +339,42 @@ silv_density_sdi_auto <- function( } sdi_values <- rep(NA_real_, n_trees) + beta_values <- rep(NA_real_, n_trees) sdi_models_used <- rep(NA_character_, n_trees) unique_combos <- unique(data.frame( species = species, + country = country, region = region, stringsAsFactors = FALSE )) for (i in seq_len(nrow(unique_combos))) { sp <- unique_combos$species[i] + cnt <- unique_combos$country[i] reg <- unique_combos$region[i] - # Handle NA in region + idx_sp <- species == sp + + if (is.na(cnt)) { + idx_cnt <- is.na(country) + cnt_arg <- NULL + } else { + idx_cnt <- country == cnt + cnt_arg <- cnt + } + if (is.na(reg)) { - idx <- which(species == sp & is.na(region)) + idx_reg <- is.na(region) reg_arg <- NULL } else { - idx <- which(species == sp & region == reg) + idx_reg <- region == reg reg_arg <- reg } - best_model_info <- .auto_select_sdi_beta(sp, reg_arg) + idx <- which(idx_sp & idx_cnt & idx_reg) + + best_model_info <- .auto_select_sdi_beta(sp, cnt_arg, reg_arg) if (best_model_info$is_fallback && !quiet) { if (best_model_info$fallback_type == "default") { @@ -279,15 +383,19 @@ silv_density_sdi_auto <- function( cli::cli_alert_info("Exact species {.val {sp}} not found. Using genus fallback {.val {best_model_info$matched_species}} from {.val {best_model_info$model}}.") } else if (best_model_info$fallback_type == "region") { cli::cli_alert_info("Exact region {.val {reg_arg}} not found for {.val {sp}}. Using fallback region {.val {best_model_info$matched_region}} from {.val {best_model_info$model}}.") + } else if (best_model_info$fallback_type == "country") { + cli::cli_alert_info("Exact country {.val {cnt_arg}} not found for {.val {sp}}. Using fallback country {.val {best_model_info$matched_country}} from {.val {best_model_info$model}}.") } } sdi_values[idx] <- silv_density_sdi(ntrees[idx], dg[idx], beta = best_model_info$beta) + beta_values[idx] <- best_model_info$beta sdi_models_used[idx] <- best_model_info$model_desc } return(data.frame( sdi = sdi_values, + beta = beta_values, sdi_model = sdi_models_used, stringsAsFactors = FALSE )) diff --git a/data/sdi_coefficients.rda b/data/sdi_coefficients.rda new file mode 100644 index 0000000000000000000000000000000000000000..a8b8fb23bf13ca167ac6355a127ece45e58102b7 GIT binary patch literal 1160 zcmV;31b6#FT4*^jL0KkKSvhhfQ~&~U|NsBLO+-K^|L`sTe?Y(Q|L{;0NCG1>8bHW^ z5CJd<;kVEP-aJv003=mb21t|-Q`7(e00x1e2AT~IQ`Be!Y5)MypdO7)>W2UT007Vc z000JqKmY&$$OAwC-~a#s8UO$Q0MKXv000>PXaESL5@|Iyk&`1vnq&h_4FCWGAOHX~ zXaE2JMI;&oH8L?N_^Ij|Hls~6003ehk?H}p01W^JgFsdN%8<7u8L*@Rh)ijc5HMo^ z^2aO^I^9l99DGv|>i@jikiG}0+S%Yz3Lf4mJOav0MIG2Twqt;~=lnMo5ysdX~KX=+~bSwm`7g$ZgvmZd2vFp?XX z86*KznG}S9z-<3X=4E^>xp9(n?_l7!ZY3U&Y@1TVuX;3S?bDk#4?PRalHOWKO&u^3 znqf)jhdGCJ2_$xT>>*Zpn|%!mrCQu>vuSnA*IJZI*HU;WS)N}gP`WgN$}I0pNwq4RTN5|9nn-$GLSB-5nbe;t z>m?H>*ZIa{p)pRvY@sMkaF9AyR{_&Qf05MF*8`8<-Z)272R$|;MFY&5l9xL1BZ5E~ zJIcH7F*?cP+i_i~Mbg(waPxr+(6zlNPU(c_UN`pkLSKSic75BnhO~eK>d>1~mXWF) zxFjY)ocC)QCI%8GKpF$9ouMv`NdpIvoNjA}5H}TQCtM{=lOjS%HY7G{fHMgRoN)qH zJ}p_Z3@w0YNy~pyPk|#Te|q-XnwxG(=mM0nCYTgM^rxxWam!J7mZp^6(~iWoHQKfT zG@2gDQ;Hiv!8fJ9IVQmxos_mP!n;#7Kz0U6RCcREhV2PG*+8HW1h5kNScJw(+iY04 z=~hc5DwuYGoCMN|0gTwsRFp}5QgQ-Pa6&-E_H3ABz#Dn>fjCPb>7-VbvXWqxrSC0C zdf1l88N`N>875THniEA#wT>un=}odrrvxZTHY{B_{?sc>?st4|8|g||84Lj?FQKq~ z^)VrgX*bhS!we-g6G|F%&cM>mzTa}<S}0(MhMbXmS&#$@2U97q?6z6~G>z9@#UweUrKkpaYT z8ju7ZQsuk>u8bHW^ z5CJd<;kVEPo&dZVMF@)_B`K)T&;jZI0g<3IGypw700006jWiENru9RB0B8Vc13&-< zj6eX;00000&;x(~XaHyfKmZ1eKmgDH0002c0ThTTc#Njg6w}HX27mwn0gwYA4H$-i z002c11e+w&Q!uIRNwkkpVF92Fk5Jl;0B8n)00E(*&+=*+F{bVn0Hh{*$S^SHg*y@0 zf=*sOPG?`8a;K%-?CznK2a~tM)2JwWS<)~Is;R>ZW?-&-wF3SVm$SKMO@pGICHxl7 zS4envnCs5W=MIhqbOv|@=ZdhZAB=>Qf?yOT1w&&1LJ*yJ9AG~9Q?&_}1^;jRy1mcn z(DNMu2|8V)M2CgyG4G|#yFirZBvBfr>BJOEfVn18#DOeIB+CgSYXd0&WR_N=b8r1G zTZ{6x*7W4Ty2Qcd4kqNaEbUS75?(+rbEj1xOn^d5=!c-d5Dz27g<5gSOIM)FYIWp8 zMI2;A15I%zCMS%c#K#uYz5W<-U*MlWi2`6*7#(I$k2-|6&m{E7Wim}jb!DN{>mlNg zVp-17!KB5;)csNfS?o50@)4MFWO;B(xo=3i!`;9XdSTv z-U6)w%6>_D#6=*g0DtOKt$YgZj)>@t ji{xX8&{MY{iRE(JNE`Chn#)oqOa3n8ig2MI@PNw Date: Tue, 30 Jun 2026 16:06:03 +0200 Subject: [PATCH 5/6] fix: configure Quarto and install R package in pkgdown workflow, and remove duplicates in _pkgdown.yml --- .github/workflows/pkgdown.yaml | 6 +++++- _pkgdown.yml | 3 --- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index bfc9f4d..0715b94 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -31,10 +31,14 @@ jobs: with: use-public-rspm: true + - name: Set up Quarto + uses: quarto-dev/quarto-actions/setup@v2 + - uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: any::pkgdown, local::. + extra-packages: any::pkgdown, any::quarto, local::. needs: website + install-quarto: true - name: Build site run: pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE) diff --git a/_pkgdown.yml b/_pkgdown.yml index 2b7b908..4ea65d1 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -27,8 +27,6 @@ reference: - silv_predict_biomass_components - silv_predict_carbon - silv_predict_carbon_auto - - silv_predict_biomass_auto - - silv_predict_biomass_components - silv_predict_height - silv_predict_snfi_volume - silv_snfi_provinces @@ -78,7 +76,6 @@ reference: - silv_density_hart - silv_density_sdi - silv_spacing_index - - silv_ntrees_ha - title: "Inventory Sample Size" desc: Calculate and explore forest inventory sample size From fcf83736c0e8d8683282b38086b7bd2ecc691c68 Mon Sep 17 00:00:00 2001 From: aitorvv Date: Tue, 30 Jun 2026 16:55:50 +0200 Subject: [PATCH 6/6] fix: add sdi_coefficients, silv_density_sdi_auto, and silv_density_sdi_class to _pkgdown.yml --- _pkgdown.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index 4ea65d1..e93c1f8 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -75,6 +75,8 @@ reference: - silv_density_ntrees_ha - silv_density_hart - silv_density_sdi + - silv_density_sdi_auto + - silv_density_sdi_class - silv_spacing_index - title: "Inventory Sample Size" @@ -100,6 +102,7 @@ reference: contents: - biomass_models - carbon_models + - sdi_coefficients - snfi3_volume_coefficients - snfi4_volume_coefficients - inventory_samples