diff --git a/NAMESPACE b/NAMESPACE index fcc5e98..a6879cc 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_sdimax) export(silv_diametric_class) export(silv_dominant_height) export(silv_lorey_height) diff --git a/R/data.R b/R/data.R index b739df3..f34b562 100644 --- a/R/data.R +++ b/R/data.R @@ -107,3 +107,29 @@ #' 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" + + +#' Maximum stand density index (SDImax) models +#' +#' Coefficients for calculating maximum stand density index (SDImax) from +#' Rodríguez de Prado (2020). +#' +#' @format A `tibble` with 88 rows and 13 variables: +#' \describe{ +#' \item{article_id}{Character. Identifier of the article.} +#' \item{title}{Character. Title of the article.} +#' \item{doi_url}{Character. DOI URL of the article.} +#' \item{country}{Character. Country where the study was conducted.} +#' \item{species}{Character. Tree species scientific name.} +#' \item{model_name}{Character. Name of the model/equation variant (e.g. "basic", "P1", "MXT3").} +#' \item{a0}{Numeric. Coeffient a0.} +#' \item{a1}{Numeric. Coeffient a1 (0 if not used/applicable).} +#' \item{b0}{Numeric. Coeffient b0.} +#' \item{b1}{Numeric. Coeffient b1 (0 if not used/applicable).} +#' \item{aic}{Numeric. Akaike Information Criterion.} +#' \item{pseudo_r2}{Numeric. Pseudo R-squared value.} +#' \item{q_index}{Numeric. Q index value.} +#' } +#' @references +#' Rodríguez-de-Prado, M., et al. (2020). Potential climatic influence on maximum stand carrying capacity for 15 Mediterranean coniferous and broadleaf species. Forest Ecology and Management, 458, 117824. +"sdimax_models" diff --git a/R/metrics-stand-density.R b/R/metrics-stand-density.R index 6af2fb1..e010d10 100644 --- a/R/metrics-stand-density.R +++ b/R/metrics-stand-density.R @@ -180,5 +180,113 @@ silv_density_hart <- function( } +#' Calculates the Maximum Stand Density Index (SDImax) +#' +#' The Maximum Stand Density Index (SDImax) represents the maximum stand carrying capacity, +#' calculated using coefficients from Rodríguez de Prado (2020) by default. +#' +#' @param species Character vector. Scientific names of the tree species. +#' @param model Character. The source article or model database (default is \code{"rodriguez-prado-2020"}). +#' @param climatic_model Character. The specific climate-dependent model name (e.g. \code{"P1"}, \code{"MXT3"}). +#' Required if \code{clim_value} is provided, and must not be \code{"basic"}. +#' @param clim_value Numeric vector. Values of the climatic variable corresponding to the selected +#' climate model. If \code{NULL} (default), the reference model (\code{"basic"}) is calculated. +#' +#' @return A numeric vector representing the SDImax for each species. +#' @export +#' +#' @details +#' If \code{clim_value} is \code{NULL}, the function computes the reference SDImax (SDImaxREF) +#' based on the "basic" model parameters: +#' \deqn{SDImaxREF = exp(a0 + b0 * log(25.4))} +#' If \code{clim_value} is provided, a climate-dependent model must be specified in \code{climatic_model}, +#' and the climate-dependent SDImax is calculated as: +#' \deqn{SDImax(Clim) = exp((a0 + a1 * log(clim_value)) + (b0 + b1 * clim_value) * log(25.4))} +#' +#' @references +#' Rodríguez-de-Prado, M., et al. (2020). Potential climatic influence on maximum stand carrying capacity for 15 Mediterranean coniferous and broadleaf species. Forest Ecology and Management, 458, 117824. +#' +#' @examples +#' ## Calculate reference SDImax for Pinus sylvestris +#' silv_density_sdimax("Pinus sylvestris") +#' +#' ## Calculate climate-dependent SDImax for Pinus canariensis using model P1 +#' silv_density_sdimax("Pinus canariensis", climatic_model = "P1", clim_value = 400) +silv_density_sdimax <- function( + species, + model = "rodriguez-prado-2020", + climatic_model = NULL, + clim_value = NULL +) { + # 0. Validate inputs + if (!is.character(species)) { + cli::cli_abort("{.arg species} must be a character vector.") + } + if (!is.character(model) || length(model) != 1) { + cli::cli_abort("{.arg model} must be a single character string.") + } + if (!is.null(climatic_model) && (!is.character(climatic_model) || length(climatic_model) != 1)) { + cli::cli_abort("{.arg climatic_model} must be a single character string.") + } + if (!is.null(clim_value) && !is.numeric(clim_value)) { + cli::cli_abort("{.arg clim_value} must be a numeric vector.") + } + + # 1. Determine target climatic model name + if (is.null(clim_value)) { + if (!is.null(climatic_model) && climatic_model != "basic") { + cli::cli_abort("Argument {.arg clim_value} is required when using climate-dependent models.") + } + target_model <- "basic" + } else { + if (is.null(climatic_model) || climatic_model == "basic") { + cli::cli_abort("Argument {.arg climatic_model} must be specified and cannot be 'basic' when {.arg clim_value} is provided.") + } + target_model <- climatic_model + } + + # If clim_value is provided, align lengths with species vector (if necessary) + if (!is.null(clim_value)) { + if (length(clim_value) == 1 && length(species) > 1) { + clim_value <- rep(clim_value, length(species)) + } + if (length(species) != length(clim_value)) { + cli::cli_abort("{.arg species} and {.arg clim_value} must have the same length.") + } + } + + # 2. Get coefficients from internal dataset + coefs_tbl <- sdimax_models[sdimax_models$article_id == model & sdimax_models$model_name == target_model, ] + + if (nrow(coefs_tbl) == 0) { + cli::cli_abort("No coefficients found for model {.val {model}} and climatic_model {.val {target_model}}.") + } + + # Check that all requested species are supported + missing_species <- species[!species %in% coefs_tbl$species] + if (length(missing_species) > 0) { + cli::cli_abort("The following species are not supported by this model: {.val {unique(missing_species)}}.") + } + + # Match species to extract coefficients + matched_indices <- match(species, coefs_tbl$species) + a0 <- coefs_tbl$a0[matched_indices] + a1 <- coefs_tbl$a1[matched_indices] + b0 <- coefs_tbl$b0[matched_indices] + b1 <- coefs_tbl$b1[matched_indices] + + # 3. Calculate SDImax + if (target_model == "basic") { + # Reference SDImax + sdimax <- exp(a0 + (b0 * log(25.4))) + } else { + # Climate-dependent SDImax + sdimax <- exp((a0 + (a1 * log(clim_value))) + ((b0 + (b1 * clim_value)) * log(25.4))) + } + + return(sdimax) +} + + diff --git a/R/zzz.R b/R/zzz.R index 1df35d8..f3b312a 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -19,6 +19,7 @@ utils::globalVariables( "ntrees", "ntrees_ha", "remaining_to_extract", + "sdimax_models", "weighted.mean", ".cumtrees", ".data", diff --git a/_pkgdown.yml b/_pkgdown.yml index 2b7b908..b175d4d 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -77,6 +77,7 @@ reference: - silv_density_ntrees_ha - silv_density_hart - silv_density_sdi + - silv_density_sdimax - silv_spacing_index - silv_ntrees_ha @@ -103,6 +104,7 @@ reference: contents: - biomass_models - carbon_models + - sdimax_models - snfi3_volume_coefficients - snfi4_volume_coefficients - inventory_samples diff --git a/data/sdimax_models.rda b/data/sdimax_models.rda new file mode 100644 index 0000000..129f8a7 Binary files /dev/null and b/data/sdimax_models.rda differ diff --git a/inst/2020_RodriguezdePrado_sdi_models_spain.csv b/inst/2020_RodriguezdePrado_sdi_models_spain.csv new file mode 100644 index 0000000..d3bdbd5 --- /dev/null +++ b/inst/2020_RodriguezdePrado_sdi_models_spain.csv @@ -0,0 +1,89 @@ +Species,Model,a0,a1,b0,b1,AIC,pseudoR2,Q_index +Pinus canariensis,basic,12.672,,-1.8226,,2616.9,0.3378, +Pinus canariensis,P1,3.639,2.448,-2.0891,,2320.8,0.4178,0.305 +Pinus canariensis,PWM,4.176,2.059,-1.9567,,2347.4,0.4111,0.251 +Pinus canariensis,P,13.161,,-2.6082,0.0015,2364.4,0.4067,0.245 +Pinus canariensis,P2,12.989,,-2.3961,0.0075,2420.9,0.3921,0.189 +Pinus canariensis,M1,11.738,1.061,-1.85,,2431.5,0.3893,0.189 +Pinus halepensis,basic,11.982,,-1.776,,12622.5,0.3388, +Pinus halepensis,M,9.241,0.886,-1.5559,-0.0095,12325.5,0.3549,0.063 +Pinus halepensis,MXT3,96.948,-14.977,-1.7045,,12368.1,0.3526,0.079 +Pinus halepensis,MXT4,105.595,-16.445,-1.7171,,12383.7,0.3517,0.077 +Pinus halepensis,MXTWM,100.504,-15.542,-1.7134,,12394.4,0.3512,0.073 +Pinus halepensis,PWM,8.722,0.784,-1.6057,-0.0026,12401.9,0.3509,0.051 +Pinus nigra,basic,12.756,,-1.8346,,5117.9,0.2965, +Pinus nigra,MXT3,140.953,–22.536,-1.9324,,5010.9,0.3128,0.123 +Pinus nigra,MXT,154.667,-24.995,-1.9154,,5028.5,0.3102,0.119 +Pinus nigra,MXT4,104.61,-16.094,-1.9119,,5045.9,0.3076,0.091 +Pinus nigra,MXT2,13.019,,5.7005,-0.0268,5046.7,0.3075,0.119 +Pinus nigra,P2,11.821,0.29,-1.8973,,5047.8,0.3073,0.098 +Pinus pinaster,basic,13.096,,-1.9063,,10593,0.2716, +Pinus pinaster,MXT,13.446,,4.177,-0.0213,10229,0.3011,0.129 +Pinus pinaster,MXT3,13.365,,3.5759,-0.019,10241.6,0.3001,0.128 +Pinus pinaster,T3,13.324,,3.911,-0.0206,10296.4,0.2958,0.121 +Pinus pinaster,MXT4,13.462,,2.6955,-0.0159,10307.5,0.2949,0.114 +Pinus pinaster,MXT2,13.389,,3.3318,-0.0187,10317.1,0.2941,0.114 +Pinus pinea,basic,13.562,,-2.1855,,3270.9,0.3887, +Pinus pinea,P4,15.072,-0.46,-2.4379,0.0093,3139.5,0.4185,0.262 +Pinus pinea,M4,13.531,-0.467,-2.4556,0.2919,3144,0.4176,0.257 +Pinus pinea,P,13.213,,-2.2271,0.0003,3210.7,0.4026,0.131 +Pinus pinea,TAR,77.368,-11.127,-2.279,,3213.2,0.402,0.143 +Pinus pinea,M,13.304,,-2.2518,0.0077,3216.7,0.4013,0.155 +Pinus radiata,basic,12.947,,-1.8254,,1432.8,0.3723, +Pinus radiata,PET3,110.968,-21.507,-8.049,0.0652,1402.4,0.3845,0.058 +Pinus radiata,PET4,88.959,-16.269,-6.5496,0.0441,1409.2,0.3821,0.062 +Pinus radiata,PET1,6.92,1.675,-1.3894,-0.0119,1421.2,0.3778,0.02 +Pinus sylvestris,basic,12.685,,-1.7524,,7718.9,0.368, +Pinus sylvestris,TAR,66.47,-9.442,-1.7478,,7594.7,0.3777,0.078 +Pinus sylvestris,MNTCM,617.791,-108.147,-40.0934,0.1425,7630.1,0.3751,0.109 +Pinus sylvestris,MXTWM,74.54,-10.872,-1.7675,,7637.6,0.3744,0.075 +Pinus sylvestris,MXT4,71.686,-10.376,-1.7699,,7643.9,0.3739,0.073 +Pinus sylvestris,MXT3,58.945,-8.154,-1.7767,,7653,0.3732,0.064 +Pinus uncinata,basic,12.519,,-1.7336,,556.6,0.4414, +Pinus uncinata,PET3,12.918,,-1.6378,-0.0031,534.6,0.4586,0.068 +Pinus uncinata,PET4,16.777,-0.838,-1.8979,,535.5,0.458,0.063 +Pinus uncinata,PET,12.899,,-1.6288,-0.0004,535.6,0.4578,0.108 +Pinus uncinata,PET2,12.908,,-1.6784,-0.0077,536.7,0.4571,0.062 +Pinus uncinata,P2,11.386,0.364,-1.9112,,538.1,0.4561,0.052 +Fagus sylvatica,basic,13.17,,-1.9471,,1577.1,0.5137, +Fagus sylvatica,MXT3,12.87,,2.088,-0.0137,1507.5,0.529,0.085 +Fagus sylvatica,T3,12.813,,2.0872,-0.0138,1510.2,0.5285,0.085 +Fagus sylvatica,MXT2,75.624,-11.138,-1.836,,1512.2,0.5281,0.07 +Fagus sylvatica,PET1,12.911,,-1.5935,-0.0085,1514.5,0.5276,0.061 +Fagus sylvatica,M1,12.133,0.671,-2.0013,,1514.9,0.5275,0.135 +Quercus faginea,basic,12.097,,-1.7055,,2003.5,0.1811, +Quercus faginea,MXTWM,247.037,-41.233,-1.7874,,1883.7,0.2508,0.315 +Quercus faginea,TAR,12.606,,12.9044,-0.0495,1886.9,0.249,0.35 +Quercus faginea,MXT4,254.074,-42.519,-1.7485,,1899.6,0.242,0.315 +Quercus faginea,T4,271.627,-45.75,-1.6856,,1910.6,0.2359,0.303 +Quercus faginea,M,9.667,0.812,-1.8657,,1915.9,0.2329,0.188 +Quercus ilex,basic,12.508,,-2.0951,,8099.8,0.5025, +Quercus ilex,PET3,11.777,,-1.3094,-0.0044,7398.6,0.5487,0.211 +Quercus ilex,PET,11.773,,-1.405,-0.0004,7449.7,0.5455,0.207 +Quercus ilex,MXT3,11.899,,5.0064,-0.0234,7474.1,0.544,0.215 +Quercus ilex,MXTWM,11.969,,4.7651,-0.0223,7484.2,0.5433,0.172 +Quercus ilex,PET2,11.865,,-1.5025,-0.0087,7491.1,0.5429,0.159 +Quercus petraea,basic,12.277,,-1.6777,,431.3,0.3877, +Quercus petraea,MXT,-489.861,88.759,36.5003,-0.1334,357.6,0.4954,0.242 +Quercus petraea,MXT4,12.593,,9.0312,-0.037,358.5,0.4917,0.247 +Quercus petraea,MXT3,12.615,,7.5139,-0.0323,360,0.4899,0.23 +Quercus petraea,MXTWM,12.382,,8.8624,-0.036,360.8,0.4889,0.227 +Quercus petraea,T4,12.674,,11.0925,-0.0446,361.6,0.4878,0.24 +Quercus pyrenaica,basic,12.271,,-1.7203,,4718.4,0.2962, +Quercus pyrenaica,T4,-187.581,35.255,17.946,-0.0679,4537.2,0.33,0.213 +Quercus pyrenaica,MNT4,12.312,,7.1163,-0.0309,4566.5,0.3244,0.186 +Quercus pyrenaica,MXTWM,12.335,,5.632,-0.025,4570,0.3238,0.191 +Quercus pyrenaica,MXT3,-310.973,57.023,24.1039,-0.0892,4577.6,0.3228,0.204 +Quercus pyrenaica,MXT4,12.328,,5.5596,-0.0248,4578.1,0.3223,0.182 +Quercus robur,basic,12.043,,-1.6698,,1017.7,0.4394, +Quercus robur,MNT3,-795.789,143.317,49.1578,-0.1812,974.7,0.4624,0.12 +Quercus robur,MNT,-820.659,147.74,51.1787,-0.1885,981.1,0.4594,0.125 +Quercus robur,MNT2,-605.574,109.939,37.8316,-0.1435,985.5,0.4572,0.123 +Quercus robur,MNT4,-1112.201,198.611,70.2864,-0.2505,989.2,0.4554,0.131 +Quercus robur,MNT1,-624.82,113.08,39.0364,-0.1458,993.6,0.4533,0.115 +Quercus suber,basic,12.704,,-1.9674,,1340.2,0.4839, +Quercus suber,PET3,11.948,,-1.2349,-0.0043,1233.6,0.5231,0.176 +Quercus suber,MXTWM,12.097,,9.7879,-0.0385,1235.9,0.5223,0.208 +Quercus suber,PET4,11.846,,-1.3656,-0.0025,1239.2,0.5211,0.15 +Quercus suber,MXT4,-670.091,119.608,43.6583,-0.1515,1239.5,0.5217,0.147 +Quercus suber,MXT3,12.343,,9.4775,-0.0384,1243.8,0.5195,0.185 diff --git a/man/sdimax_models.Rd b/man/sdimax_models.Rd new file mode 100644 index 0000000..0f86c31 --- /dev/null +++ b/man/sdimax_models.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{sdimax_models} +\alias{sdimax_models} +\title{Maximum stand density index (SDImax) models} +\format{ +A \code{tibble} with 88 rows and 13 variables: +\describe{ +\item{article_id}{Character. Identifier of the article.} +\item{title}{Character. Title of the article.} +\item{doi_url}{Character. DOI URL of the article.} +\item{country}{Character. Country where the study was conducted.} +\item{species}{Character. Tree species scientific name.} +\item{model_name}{Character. Name of the model/equation variant (e.g. "basic", "P1", "MXT3").} +\item{a0}{Numeric. Coeffient a0.} +\item{a1}{Numeric. Coeffient a1 (0 if not used/applicable).} +\item{b0}{Numeric. Coeffient b0.} +\item{b1}{Numeric. Coeffient b1 (0 if not used/applicable).} +\item{aic}{Numeric. Akaike Information Criterion.} +\item{pseudo_r2}{Numeric. Pseudo R-squared value.} +\item{q_index}{Numeric. Q index value.} +} +} +\usage{ +sdimax_models +} +\description{ +Coefficients for calculating maximum stand density index (SDImax) from +Rodríguez de Prado (2020). +} +\references{ +Rodríguez-de-Prado, M., et al. (2020). Potential climatic influence on maximum stand carrying capacity for 15 Mediterranean coniferous and broadleaf species. Forest Ecology and Management, 458, 117824. +} +\keyword{datasets} diff --git a/man/silv_density_sdimax.Rd b/man/silv_density_sdimax.Rd new file mode 100644 index 0000000..50d4c72 --- /dev/null +++ b/man/silv_density_sdimax.Rd @@ -0,0 +1,49 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/metrics-stand-density.R +\name{silv_density_sdimax} +\alias{silv_density_sdimax} +\title{Calculates the Maximum Stand Density Index (SDImax)} +\usage{ +silv_density_sdimax( + species, + model = "rodriguez-prado-2020", + climatic_model = NULL, + clim_value = NULL +) +} +\arguments{ +\item{species}{Character vector. Scientific names of the tree species.} + +\item{model}{Character. The source article or model database (default is \code{"rodriguez-prado-2020"}).} + +\item{climatic_model}{Character. The specific climate-dependent model name (e.g. \code{"P1"}, \code{"MXT3"}). +Required if \code{clim_value} is provided, and must not be \code{"basic"}.} + +\item{clim_value}{Numeric vector. Values of the climatic variable corresponding to the selected +climate model. If \code{NULL} (default), the reference model (\code{"basic"}) is calculated.} +} +\value{ +A numeric vector representing the SDImax for each species. +} +\description{ +The Maximum Stand Density Index (SDImax) represents the maximum stand carrying capacity, +calculated using coefficients from Rodríguez de Prado (2020) by default. +} +\details{ +If \code{clim_value} is \code{NULL}, the function computes the reference SDImax (SDImaxREF) +based on the "basic" model parameters: +\deqn{SDImaxREF = exp(a0 + b0 * log(25.4))} +If \code{clim_value} is provided, a climate-dependent model must be specified in \code{climatic_model}, +and the climate-dependent SDImax is calculated as: +\deqn{SDImax(Clim) = exp((a0 + a1 * log(clim_value)) + (b0 + b1 * clim_value) * log(25.4))} +} +\examples{ +## Calculate reference SDImax for Pinus sylvestris +silv_density_sdimax("Pinus sylvestris") + +## Calculate climate-dependent SDImax for Pinus canariensis using model P1 +silv_density_sdimax("Pinus canariensis", climatic_model = "P1", clim_value = 400) +} +\references{ +Rodríguez-de-Prado, M., et al. (2020). Potential climatic influence on maximum stand carrying capacity for 15 Mediterranean coniferous and broadleaf species. Forest Ecology and Management, 458, 117824. +} diff --git a/scripts/regenerate_sdimax_models.R b/scripts/regenerate_sdimax_models.R new file mode 100644 index 0000000..f3899ff --- /dev/null +++ b/scripts/regenerate_sdimax_models.R @@ -0,0 +1,56 @@ +# Script to regenerate sdimax_models.rda from inst/2020_RodriguezdePrado_sdi_models_spain.csv + +library(dplyr) +library(readr) + +# Read the raw CSV +raw_data <- read_csv("inst/2020_RodriguezdePrado_sdi_models_spain.csv", show_col_types = FALSE) + +# Helper function to clean numeric fields (e.g. replace en-dashes with standard minus signs) +clean_num <- function(x) { + if (is.character(x)) { + x <- gsub("–", "-", x) + x <- gsub("—", "-", x) + } + suppressWarnings(as.numeric(x)) +} + +# Clean, format, and add metadata +sdimax_models <- raw_data %>% + transmute( + article_id = "rodriguez-prado-2020", + title = "Potential climatic influence on maximum stand carrying capacity for 15 Mediterranean coniferous and broadleaf species", + doi_url = "https://doi.org/10.1016/j.foreco.2019.117824", + country = "Spain", + species = as.character(Species), + model_name = as.character(Model), + a0 = clean_num(a0), + a1 = clean_num(a1), + b0 = clean_num(b0), + b1 = clean_num(b1), + aic = clean_num(AIC), + pseudo_r2 = clean_num(pseudoR2), + q_index = clean_num(Q_index) + ) %>% + mutate( + # Fill NAs in climate slope/interaction terms with 0 + a1 = ifelse(is.na(a1), 0, a1), + b1 = ifelse(is.na(b1), 0, b1) + ) %>% + as_tibble() + +# Validation checks +cat("=== VALIDATING sdimax_models ===\n") +cat("Dimensions:", nrow(sdimax_models), "x", ncol(sdimax_models), "\n") +cat("Columns:", paste(names(sdimax_models), collapse = ", "), "\n") +cat("Unique species count:", length(unique(sdimax_models$species)), "\n") +cat("Checking for NAs in critical columns:\n") +cat(" - a0 NAs:", sum(is.na(sdimax_models$a0)), "\n") +cat(" - b0 NAs:", sum(is.na(sdimax_models$b0)), "\n") +cat(" - a1 NAs:", sum(is.na(sdimax_models$a1)), "\n") +cat(" - b1 NAs:", sum(is.na(sdimax_models$b1)), "\n") + +# Save dataset to data/ +cat("\nSaving sdimax_models.rda...\n") +usethis::use_data(sdimax_models, overwrite = TRUE) +cat("Done!\n") diff --git a/tests/testthat/test-stand-density.R b/tests/testthat/test-stand-density.R index c54fd2c..d627e18 100644 --- a/tests/testthat/test-stand-density.R +++ b/tests/testthat/test-stand-density.R @@ -79,3 +79,61 @@ test_that("Errors work", { expect_error(silv_density_hart(17.8, "400")) expect_error(silv_density_hart(c(17.8, 20.5), 400)) }) + + +# 3. silv_density_sdimax -------------------------------------------------- + +test_that("SDImax reference model calculations are correct", { + # Pinus sylvestris reference: a0 = 12.685, b0 = -1.7524 + # exp(12.685 - 1.7524 * log(25.4)) = 1114.77 + expect_equal( + silv_density_sdimax("Pinus sylvestris"), + 1114.77, + tolerance = 0.01 + ) + + # Pinus canariensis reference: a0 = 12.672, b0 = -1.8226 + # exp(12.672 - 1.8226 * log(25.4)) = 876.86 + expect_equal( + silv_density_sdimax("Pinus canariensis"), + 876.86, + tolerance = 0.01 + ) + + # Vectorized calculation + expect_equal( + silv_density_sdimax(c("Pinus sylvestris", "Pinus canariensis")), + c(1114.77, 876.86), + tolerance = 0.01 + ) +}) + +test_that("SDImax climate-dependent calculations are correct", { + # Pinus canariensis model P1 with clim_value = 400: + # a0 = 3.639, a1 = 2.448, b0 = -2.0891, b1 = 0 + # exp((3.639 + 2.448 * log(400)) + (-2.0891) * log(25.4)) = 103600.4 + expect_equal( + silv_density_sdimax("Pinus canariensis", climatic_model = "P1", clim_value = 400), + 103600.4, + tolerance = 0.1 + ) +}) + +test_that("SDImax error handling works", { + # Non-character species + expect_error(silv_density_sdimax(123)) + + # Unsupported species + expect_error(silv_density_sdimax("Pinus nonexistus")) + + # Missing clim_value when climate model is requested + expect_error(silv_density_sdimax("Pinus canariensis", climatic_model = "P1")) + + # clim_value provided but climatic_model is basic/NULL + expect_error(silv_density_sdimax("Pinus canariensis", clim_value = 400)) + expect_error(silv_density_sdimax("Pinus canariensis", climatic_model = "basic", clim_value = 400)) + + # Length mismatch between species and clim_value + expect_error(silv_density_sdimax(c("Pinus canariensis", "Pinus sylvestris"), climatic_model = "P1", clim_value = c(400, 500, 600))) +}) +