From 2e1ec09ea237e38357c551a7493d43aa04bad729 Mon Sep 17 00:00:00 2001 From: aitorvv Date: Tue, 30 Jun 2026 16:39:19 +0200 Subject: [PATCH 1/3] feat: implement silv_density_sdimax and add sdimax_models dataset --- NAMESPACE | 1 + R/data.R | 26 +++++ R/metrics-stand-density.R | 108 ++++++++++++++++++ R/zzz.R | 1 + data/sdimax_models.rda | Bin 0 -> 4261 bytes ...2020_RodriguezdePrado_sdi_models_spain.csv | 89 +++++++++++++++ man/sdimax_models.Rd | 35 ++++++ man/silv_density_sdimax.Rd | 49 ++++++++ scripts/regenerate_sdimax_models.R | 56 +++++++++ tests/testthat/test-stand-density.R | 58 ++++++++++ 10 files changed, 423 insertions(+) create mode 100644 data/sdimax_models.rda create mode 100644 inst/2020_RodriguezdePrado_sdi_models_spain.csv create mode 100644 man/sdimax_models.Rd create mode 100644 man/silv_density_sdimax.Rd create mode 100644 scripts/regenerate_sdimax_models.R 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/data/sdimax_models.rda b/data/sdimax_models.rda new file mode 100644 index 0000000000000000000000000000000000000000..129f8a73e8a87c2707bf64fcceef41add6d056bb GIT binary patch literal 4261 zcmV;W5L)j-T4*^jL0KkKS)6QA*8ml1fB*mg|NsC0|NsC0|NsB@|NsC0|NsC0|NsC0 z|NsC0|Nr0)?|ghokTIUcZA_k;n?p&bp{7&P zV@8n3%}++AN0NG*L6Dm@P3sX#Aw~m`05>G{Ax438qXn9-wKbq{zci z$j|`Mv`mbKhpB=98V^u1G}A_nq6~vT0BFz}10Z^jQMCq|WW)f{F)*7znr#{Z7)%N2 zF)$+lO)_8vzzrh=^a5Z2hS1Q+G$sH5000vsMggD=27pWe00fAk1gY&)^*sq2)74LC zQ`Ge`X^78M#;59sr>Y*Bk5gqmLFuT|^-US3r|O54+ImyTq|?zhq3Ul*9;c!PKzg2| zO$Vu_^+%}D>ODb-dW?*kG-;4}42bmuYCI|XRXwUKY%aJxvV&^nfw|(0Yb| z4^z|t8UPvq02%dQT!|niDlZ z{VC#~sp%(@Xnv^uL}EQdK=jHyk)Ua(q|-xBMA7O7nh!!62c-2H10zO&{X}|9j0EzB z>VWkRP#&NdnHqW_Ug;4`0a+GydKD`A+}RzL^{>Z0)V+9<#8_NJ>3ygEl)Uc&O86}H zce482$crs%Vkajc&JK-NW>ae)7k8aT>EojrC7=iisgQ_36=+uL5LYJzVweStQH&U? z5X8d7^t?S8;r1R!ljb^p*GX3C*>;|n(rxI{7S;tYFe!-Tww5MgPq+W?AEZ-H67IPE z!|AiCk)ou`!V)uV&#^6*N8WW2sno?j?t(28MOhol?y6gLQB`+L?CIE1W=iau9b#qN zMKcshRTNQBQx2z14v8{}As1F%ikUGL6*>#kXT@hYs{M9Cwe4J4m_}WS(i$NmrkQ8N zGvZ}TVs2!&^(I@olNIM@MoI0pGP2EIG3x3om=emvUKUiaea5P21%q0D@rrZ!u8omI@#M zhlogomP$+@0|bCeN%T~3E5co(%H9}cQ6N^NOqH{vACz0N4Wh;ZOH{M^uf!uBCA`0t&GLh{tSG9*;9=_x4(hCZ+HJhSC$(hYxZ| z?JGZ02re$sATsJOkuQ8yTObv`%z>-xvir9GwG`?%`ie?nadJ#+^$rSQ}nG09JeLfyX73(nRJXflJMqMT7LY@FZXNZOYG$tZ)c)I|V z+$Q_LUias4yFSe+9l?H6J30%lWe2wKDxsmU%D(Z335)dN3{Zxq3r!0wm}!$2t-Cx- zpF@^^gm~M!;tF#v@^+Ya!k&B03VR(VuH9pYx8=nBQ})9m__h2Cp^PQQ1Zj%~HCjHd8Y zt_po1r$Q;b6u_p7htKv5+kLmcw>IW6j>n0tDY7Zt3VDG}{p|Yx(Z;ZPLhdjCKG&!q zp`7%&*GJE-wDoB6_Rs1htVslO!{ex2(u$*he5VvZ6!E^$zBBJz_JK~X>Ob2{l>W^j zOmtYU#kkymf!>OG-#xq6_nP9E(W5B(+#UCu6!;DAf%ACTz3+FGb8+sw+=+Byt99q) z$)Eu+u%_1I#|gR!2rZa*;;e{QJ8Bftu3(t2vGS`PtV-l_0%5uvOm@di!l2#xv8i4E zdX_HCC3+H*?$VpMtsG3(Ih#yU%!f+7_MP+PTUS@Fy4;e;2p!f{)ijVhC^( z%J9NG9J*(t44VETizfrvfUd&`t8LUg3uOzz(o~?8qA|^cH;#{S8<*Vvc}%HPt3yDe zM! zqaH*6J+Y`bj3-MgUG-l4;1lNj(sN3k_EVT$)S z$ZOHSt(V-=poxi}q0*92ZoPwpUNW~Djj8|&D;vC}OO;if9%;)V0G`J5I&WS{&vRY!B z!B5p%ne9PMg~k-v6woP{Q$oO}E3_oPt_$fvQjtRd5+=Z*0Fc^mBM-|`?M2R&)>c6X z4r=#neBX+5>sSAHxZa20;8V-!$h6d^{8NP#!tjMQG-wpiDS=Lx#gW`-)Czchz2msq zddydh-n9M3+2C!8WE9|Xg*P+4h!~ALMk_p=cD?S%R6Rflq6Yrl_X6(jwB%soVSSvF(3+?)rT(P6($2 zQ<0YqRPqc8Wfa;J&?(o2JV2*PDVs_u((*QTIp^|!KZz9c;+m08d2J*6g*c*}9%b2j z{_7q;F|)Ok$^17to2}*L%k_V0yL6Ki7@B-xO^*%9Tf@f+ZFf-n#-Iix)kAHQ&K6XtSPhM{%36(eFhzV1v-Rn1vjFd);22g1slMdr;0kQ5HK*bQGzxW7yZLO7EsZs)?1erF9Qdcl z3UAi$pj=#qHrH6Nu&}VOSXeMxYa;WP{q=TT$O#nR1GIqvtOSq-J*81+kJD0OzLW9+ zAq|70j($-#$QHn-&S%Ui#4bAgev9^byMmn8=wmiI*E=2;u;pg27l1vJ^SXf)`*^wt z&|o0Il_6q^>>$b~^`nk*k`^IWe0lu@*5XSPw~;kZi?QGS@G=`Ore|jg03sEyCK&R> z0S2>fsf-05D!S;%561s!03sR3`#%i>U|^ZMAfL!K<+?~w(O2zJaAaFT#KK6dn^D8= zi^Bww+bp=h>;?ChQN}{aG=kF`)29+8#V`@e9>kxPT*vtF5!fC6${ z>QjV3!ITO5=%i9m|AHi>s6SntcP^=%YJLn?50; z&h8mSKDuD*VaMT9PkQew?d##p6A}tBpki_6VN@=?##E32wMOtN?5H3GrOg~JgJOCu zCFTF#`@BwPKXv#6RdypK<7y>?V%Inr`WgRQyFKFC*@PJE%w>_Nvu^xPAT813Bn4wS zYsE$9&xL*kc)8{TsHWQMlq<8;uZnU8f)${dZENsZ3Pj&_Z)KsvbrwxwISw3ng6U`lPgP5n|KE z8VIqtBn5$^4XunIwhf02L>UlTGyp^nCW^8=JU9cteC9P!3m3E4HK6tUJB*#jEhemH zsUhI$qu}kZIO2{v^SXeXk}Wh#%9a_BIyFgXNF#3;Ufk@Ns`1sHoSDa$xn{%pT)J1>c4NkOx)X1C z@5%1j=jeHKeeXs&G?w6&)ZvLAYr}f?ydw@sBl9R@T9M$*sse?-@|#6$m=zny&B+6r6$8 z$e;sB2=@!~xuNIKk5G#z(sC3O@tD}-MAOpMhUqph!*Hv6{~JSFNR3{wodW#kau8LX z%3=30PWf0ZnqNN#9lnyRog@vHGRPT=5xExY+(3D(@`E`MS{OM1!d$m92{RdiNgx0y ztz(do04pBE2ND_%FQ1gjy>Q?wr2FnUnW^nIMoVXaX!HeNgBTbQ0HJ_(N<+{BM4brV zCg>iXK@eu3;YqsRYt-i1*>>@?Pd%_phoHHV@$%KLE71mwt#7idk4eDt<5@N;3Y-r4 zn%1qa1_QX-kQ|Z+gtHJ&f2-_rI;B$Q8OkSoXWBYeOJ|(DAY8!C7}g)r2xgYv(uMnW zWz?8in7L=ZO4}E=wU(srAcK?EZ5kzvT&vLAC-N-kHS~(d#$B`U+ETe+F^e^+yX(pu z!0O0GI%v`LfnoGN{lW@Sn9^6R)lum3v%W2@+Y!{8C#-ki1A0$NP*jtt}B!Vf^WdUH} zM5G82Lwei2(?=3iNN)E(6Klq0v)QPSBh)-zlKMU$rye* zsPhQSW@00p-x%@+>2z#UxZ_y{GRt6Z%_ujCj4ImwJ|znDDeHOy2i=<8Hn8$nRHJcR zjBFm8C`{MX`G}r8W)Y5Pd}% + 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))) +}) + From 0a03668e29ff067bab001639f481de2ead04c7ab Mon Sep 17 00:00:00 2001 From: aitorvv Date: Tue, 30 Jun 2026 16:50:20 +0200 Subject: [PATCH 2/3] fix: add silv_density_sdimax and sdimax_models to _pkgdown.yml --- _pkgdown.yml | 2 ++ 1 file changed, 2 insertions(+) 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 From d1ca5a2bae1287bf15b8c65c4baf496dae22a4c1 Mon Sep 17 00:00:00 2001 From: aitorvv96 Date: Wed, 1 Jul 2026 12:09:03 +0200 Subject: [PATCH 3/3] feat: implement Stand Density Index (SDI) autoselector and classification --- NAMESPACE | 2 + R/data.R | 18 + R/metrics-stand-density.R | 493 ++++++++++++++++++++++++---- _pkgdown.yml | 3 + data/sdi_coefficients.rda | Bin 0 -> 1160 bytes man/sdi_coefficients.Rd | 26 ++ man/silv_density_ntrees_ha.Rd | 2 +- man/silv_density_sdi.Rd | 30 +- man/silv_density_sdi_auto.Rd | 82 +++++ man/silv_density_sdi_class.Rd | 43 +++ tests/testthat/test-stand-density.R | 113 +++++-- 11 files changed, 706 insertions(+), 106 deletions(-) create mode 100644 data/sdi_coefficients.rda create mode 100644 man/sdi_coefficients.Rd create mode 100644 man/silv_density_sdi_auto.Rd create mode 100644 man/silv_density_sdi_class.Rd diff --git a/NAMESPACE b/NAMESPACE index a6879cc..e992583 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,6 +16,8 @@ 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_density_sdimax) export(silv_diametric_class) export(silv_dominant_height) diff --git a/R/data.R b/R/data.R index f34b562..e10ec85 100644 --- a/R/data.R +++ b/R/data.R @@ -109,6 +109,24 @@ "snfi4_volume_coefficients" +#' SDI beta coefficients +#' +#' Specific beta coefficients for Reineke's Stand Density Index (SDI) per +#' species, country, 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_coefficients" + + #' Maximum stand density index (SDImax) models #' #' Coefficients for calculating maximum stand density index (SDImax) from diff --git a/R/metrics-stand-density.R b/R/metrics-stand-density.R index e010d10..cc5474b 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) } - - } @@ -54,72 +49,465 @@ 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. URL: https://research.fs.usda.gov/download/treesearch/60134.pdf #' #' @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 +# --- Internal Auto Selector Helper --- + +#' @noRd +.auto_select_sdi_beta <- function(species, country = NULL, region = NULL) { + + sdi_coefficients <- silviculture::sdi_coefficients + + # 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], " (", country, ", ", region, ")") + )) + } + } + + # 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)") + )) + } + } - # 2. Classify? - if (classify) { + # 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_country = sel$country[1], + matched_region = sel$region[1], + is_fallback = TRUE, + fallback_type = "region", + model_desc = paste0(sel$article_id[1], " (", sel$country[1], ", ", sel$region[1], ")") + )) + } + + # 4. Try genus fallback (genus spp.) + genus <- strsplit(species, " ")[[1]][1] + genus_spp <- paste0(genus, " spp.") + + 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_country = country, + matched_region = "all", + is_fallback = TRUE, + fallback_type = "genus", + model_desc = paste0(sel$article_id[1], " (genus fallback: ", country, ", all regions)") + )) + } + } + + 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_country = sel$country[1], + matched_region = sel$region[1], + is_fallback = TRUE, + fallback_type = "genus", + model_desc = paste0(sel$article_id[1], " (genus fallback: ", sel$country[1], ", ", sel$region[1], ")") + )) + } + + # 5. Total fallback (default) + sel <- sdi_coefficients[sdi_coefficients$species == "default", ] + if (nrow(sel) == 0) { + 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 = sel$article_id[1], + beta = default_beta, + matched_species = "default", + matched_country = "default", + matched_region = "default", + is_fallback = TRUE, + fallback_type = "default", + model_desc = default_desc + )) +} - ## assert inputs - if (is.null(max_sdi)) cli::cli_abort("You must specify when ") - assert_positive_numeric(max_sdi, "max_sdi") +#' 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, country, and region from the internal \code{sdi_coefficients} database. +#' +#' 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 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 classify A logical value indicating whether to automatically calculate `SDImax` +#' and classify the values (default is `FALSE`). +#' @param climatic_model Character. The specific climate-dependent model name (e.g. `"P1"`, `"MXT3"`). +#' Passed to \code{\link{silv_density_sdimax}} when \code{classify = TRUE}. +#' @param clim_value Numeric vector. Values of the climatic variable corresponding to the selected +#' climate model. Passed to \code{\link{silv_density_sdimax}} when \code{classify = TRUE}. +#' @param quiet Logical. If `FALSE`, informs the user about fallbacks to genus or default models. +#' +#' @return A `data.frame` with the columns: +#' - `sdi`: The computed absolute Stand Density Index. +#' - `beta`: The beta exponent used for the calculation. +#' - `sdi_model`: The model used for beta exponent. +#' - `sdimax`: (If `classify = TRUE`) The maximum SDI for the species. +#' - `sdi_class`: (If `classify = TRUE`) The density classification. +#' +#' @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 León" +#' ) +#' +#' # With automatic classification +#' silv_density_sdi_auto( +#' ntrees = 800, +#' dg = 23.4, +#' species = "Pinus sylvestris", +#' classify = TRUE +#' ) +#' +#' @export +silv_density_sdi_auto <- function( + ntrees, + dg, + species, + country = NULL, + region = NULL, + classify = FALSE, + climatic_model = NULL, + clim_value = 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")) - ## 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" - ) - } else if (!is.null(max_sdi)) { - sdi <- (sdi / max_sdi) * 100 + 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}.") } - return(sdi) + 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) { + 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) + beta_values <- rep(NA_real_, n_trees) + sdi_models_used <- rep(NA_character_, n_trees) + + if (classify) { + sdimax_values <- rep(NA_real_, n_trees) + sdi_class_values <- 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] + + 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_reg <- is.na(region) + reg_arg <- NULL + } else { + idx_reg <- region == reg + reg_arg <- reg + } + + 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") { + 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}}.") + } 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}}.") + } + } + + current_sdi <- silv_density_sdi(ntrees[idx], dg[idx], beta = best_model_info$beta) + sdi_values[idx] <- current_sdi + beta_values[idx] <- best_model_info$beta + sdi_models_used[idx] <- best_model_info$model_desc + + if (classify) { + # Handle SDImax per species + tryCatch({ + current_clim_value <- if(!is.null(clim_value)) clim_value[idx] else NULL + current_sdimax <- silv_density_sdimax( + species = rep(sp, length(idx)), + climatic_model = climatic_model, + clim_value = current_clim_value + ) + sdimax_values[idx] <- current_sdimax + sdi_class_values[idx] <- silv_density_sdi_class(current_sdi, current_sdimax, classify = TRUE) + }, error = function(e) { + if (!quiet) { + cli::cli_warn("SDImax classification failed for {.val {sp}}: {e$message}") + } + sdimax_values[idx] <- NA_real_ + sdi_class_values[idx] <- NA_character_ + }) + } + } + + if (classify) { + return(data.frame( + sdi = sdi_values, + beta = beta_values, + sdi_model = sdi_models_used, + sdimax = sdimax_values, + sdi_class = sdi_class_values, + stringsAsFactors = FALSE + )) + } else { + return(data.frame( + sdi = sdi_values, + beta = beta_values, + sdi_model = sdi_models_used, + stringsAsFactors = FALSE + )) + } } +#' 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")) + + # 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 { + res <- rel_sdi + } + return(res) +} #' Hart or Hart-Becking spacing index @@ -163,7 +551,6 @@ silv_density_hart <- function( ntrees, which = c("hart", "hart-becking") ) { - # 0. Validate inputs assert_positive_numeric(h0, "h0") assert_positive_numeric(ntrees, "ntrees") @@ -176,10 +563,8 @@ silv_density_hart <- function( "hart-becking" = sqrt(20000 / (ntrees * sqrt(3))) / h0 * 100, cli::cli_abort("`which` must be either or ") ) - } - #' Calculates the Maximum Stand Density Index (SDImax) #' #' The Maximum Stand Density Index (SDImax) represents the maximum stand carrying capacity, @@ -286,7 +671,3 @@ silv_density_sdimax <- function( return(sdimax) } - - - - diff --git a/_pkgdown.yml b/_pkgdown.yml index b175d4d..aa0cd2b 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -77,7 +77,9 @@ reference: - silv_density_ntrees_ha - silv_density_hart - silv_density_sdi + - silv_density_sdi_auto - silv_density_sdimax + - silv_density_sdi_class - silv_spacing_index - silv_ntrees_ha @@ -104,6 +106,7 @@ reference: contents: - biomass_models - carbon_models + - sdi_coefficients - sdimax_models - snfi3_volume_coefficients - snfi4_volume_coefficients 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>u n, plot_size = c(10, 15), plot_shape = "rectangular" - ) + ) ) } diff --git a/man/silv_density_sdi.Rd b/man/silv_density_sdi.Rd index 422bc24..6800fb2 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. URL: https://research.fs.usda.gov/download/treesearch/60134.pdf } diff --git a/man/silv_density_sdi_auto.Rd b/man/silv_density_sdi_auto.Rd new file mode 100644 index 0000000..d0bfdd9 --- /dev/null +++ b/man/silv_density_sdi_auto.Rd @@ -0,0 +1,82 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/metrics-stand-density.R +\name{silv_density_sdi_auto} +\alias{silv_density_sdi_auto} +\title{Predict Stand Density Index automatically} +\usage{ +silv_density_sdi_auto( + ntrees, + dg, + species, + country = NULL, + region = NULL, + classify = FALSE, + climatic_model = NULL, + clim_value = NULL, + quiet = FALSE +) +} +\arguments{ +\item{ntrees}{Numeric vector with number of trees of the diameter class per +hectare. If \code{ntrees = NULL}, the function will assume that each diameter +corresponds to only one tree} + +\item{dg}{Numeric vector of quadratic mean diameters} + +\item{species}{A character string or vector of tree species (e.g., \code{"Pinus sylvestris"}).} + +\item{country}{A character string or vector of the country (e.g., \code{"Spain"}). +Defaults to \code{NULL} (no country specified).} + +\item{region}{A character string or vector of the region (e.g., \code{"Castilla y León"}). +Defaults to \code{NULL} (no region specified).} + +\item{classify}{A logical value indicating whether to automatically calculate \code{SDImax} +and classify the values (default is \code{FALSE}).} + +\item{climatic_model}{Character. The specific climate-dependent model name (e.g. \code{"P1"}, \code{"MXT3"}). +Passed to \code{\link{silv_density_sdimax}} when \code{classify = TRUE}.} + +\item{clim_value}{Numeric vector. Values of the climatic variable corresponding to the selected +climate model. Passed to \code{\link{silv_density_sdimax}} when \code{classify = TRUE}.} + +\item{quiet}{Logical. If \code{FALSE}, informs the user about fallbacks to genus or default models.} +} +\value{ +A \code{data.frame} with the columns: +\itemize{ +\item \code{sdi}: The computed absolute Stand Density Index. +\item \code{beta}: The beta exponent used for the calculation. +\item \code{sdi_model}: The model used for beta exponent. +\item \code{sdimax}: (If \code{classify = TRUE}) The maximum SDI for the species. +\item \code{sdi_class}: (If \code{classify = TRUE}) The density classification. +} +} +\description{ +\code{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, country, and region from the internal \code{sdi_coefficients} database. + +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). +} +\examples{ +# Calculate SDI with automatic selection +silv_density_sdi_auto( + ntrees = 800, + dg = 23.4, + species = "Pinus sylvestris", + region = "Castilla y León" +) + +# With automatic classification +silv_density_sdi_auto( + ntrees = 800, + dg = 23.4, + species = "Pinus sylvestris", + classify = TRUE +) + +} 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 d627e18..366c2df 100644 --- a/tests/testthat/test-stand-density.R +++ b/tests/testthat/test-stand-density.R @@ -77,63 +77,108 @@ 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_sdimax -------------------------------------------------- +# 3. silv_density_sdi ---------------------------------------------------- -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 +## Tests +test_that("Stand Density Index is well calculated", { + # default beta (1.605) expect_equal( - silv_density_sdimax("Pinus sylvestris"), - 1114.77, + silv_density_sdi(ntrees = 800, dg = 23.4), + 702.40, tolerance = 0.01 ) - # Pinus canariensis reference: a0 = 12.672, b0 = -1.8226 - # exp(12.672 - 1.8226 * log(25.4)) = 876.86 + # custom beta expect_equal( - silv_density_sdimax("Pinus canariensis"), - 876.86, + silv_density_sdi(ntrees = 800, dg = 23.4, beta = 1.7), + 692.68, tolerance = 0.01 ) - # Vectorized calculation + # negative beta expect_equal( - silv_density_sdimax(c("Pinus sylvestris", "Pinus canariensis")), - c(1114.77, 876.86), + silv_density_sdi(ntrees = 800, dg = 23.4, beta = -1.605), + 702.40, 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 + # with max_sdi (returns percentage) using silv_density_sdi_class + sdi_val <- silv_density_sdi(ntrees = 800, dg = 23.4) expect_equal( - silv_density_sdimax("Pinus canariensis", climatic_model = "P1", clim_value = 400), - 103600.4, - tolerance = 0.1 + 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("SDImax error handling works", { - # Non-character species - expect_error(silv_density_sdimax(123)) +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")) +}) - # 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")) +# 4. silv_density_sdi_auto ----------------------------------------------- + +test_that("silv_density_sdi_auto calculates and falls back correctly", { + + # Exact match + res1 <- silv_density_sdi_auto(800, 23.4, "Pinus sylvestris", country = "Spain", region = "Castilla y León", quiet = TRUE) + expect_equal(res1$sdi_model, "del-rio-2006 (Spain, Castilla y León)") + expect_equal(res1$beta, -1.75) + expect_true(is.numeric(res1$sdi)) + + # Region fallback ("all") + res2 <- silv_density_sdi_auto(800, 23.4, "Pinus pinaster", country = "Spain", region = "Unknown region", quiet = TRUE) + expect_equal(res2$sdi_model, "aguirre-2017 (Spain, all regions)") + expect_equal(res2$beta, -1.9477) + + # Default fallback + res3 <- silv_density_sdi_auto(800, 23.4, "Unknown species", quiet = TRUE) + expect_equal(res3$sdi_model, "reineke-1933 (-1.605)") + expect_equal(res3$beta, -1.605) + +}) - # 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)) +test_that("Errors work in silv_density_sdi_auto", { + expect_error(silv_density_sdi_auto(800, 23.4, species = 123)) + expect_error(silv_density_sdi_auto(c(800, 700), c(23.4, 25.1), species = c("Pinus sylvestris", "Pinus pinaster", "Quercus robur"))) +}) - # 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))) +test_that("silv_density_sdi_auto handles classification correctly", { + # With classification and SDImax + res_class <- silv_density_sdi_auto(800, 23.4, "Pinus sylvestris", classify = TRUE, quiet = TRUE) + expect_true("sdimax" %in% names(res_class)) + expect_true("sdi_class" %in% names(res_class)) + expect_true(is.numeric(res_class$sdimax)) + expect_true(is.character(res_class$sdi_class)) }) +# 5. silv_density_sdimax ----------------------------------------------- + +test_that("silv_density_sdimax calculates reference and climate models", { + sdimax_ref <- silv_density_sdimax("Pinus sylvestris") + expect_true(is.numeric(sdimax_ref)) + + sdimax_clim <- silv_density_sdimax("Pinus canariensis", climatic_model = "P1", clim_value = 400) + expect_true(is.numeric(sdimax_clim)) + + expect_error(silv_density_sdimax("Pinus canariensis", climatic_model = "P1")) + expect_error(silv_density_sdimax("Unknown Species")) +})