From 2e1ec09ea237e38357c551a7493d43aa04bad729 Mon Sep 17 00:00:00 2001 From: aitorvv Date: Tue, 30 Jun 2026 16:39:19 +0200 Subject: [PATCH 1/2] 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/2] 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