Add dominant diameter #1

Author: Cidree

NAMESPACE

@@ -1,11 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
 export(eq_biomass_dieguez_aranda_2009)
+export(eq_biomass_manrique_2017)
 export(eq_biomass_montero_2005)
 export(eq_biomass_ruiz_peinado_2011)
 export(eq_biomass_ruiz_peinado_2012)
 export(eq_hd_vazquez_veloso_2025)
-export(eq_manrique_2017)
 export(lid_fcov)
 export(lid_lhdi)
 export(plot)
@@ -26,6 +26,7 @@ export(silv_sample_size_stratified)
 export(silv_spacing_index)
 export(silv_sqrmean_diameter)
 export(silv_stand_basal_area)
+export(silv_stand_dominant_diameter)
 export(silv_stand_dominant_height)
 export(silv_stand_lorey_height)
 export(silv_stand_qmean_diameter)

NEWS.md

@@ -1,5 +1,5 @@
 
-# silviculture 0.2.0 (dev)
+# silviculture 0.1.9000 (dev)
 
 This new version brings new naming conventions that will be useful for sorting the package into "modules" of related functions.
 
@@ -17,6 +17,12 @@ The old functions are now deprecated and will be eliminated in a future release.
 
 * `silv_predict_height()`: estimates height from diameter, using the so-called h-d curves. The argument `equation` allows to choose which equations to use. Currently, only `eq_hd_aitor2025()` available.
 
+* `silv_stand_dominant_diameter()`: calculates dominant diameter using two methods:
+
+    - `Assman`: the mean diameter of the 100 thickest trees per hectare
+
+    - `Weise`: the quadratic mean diameter of the 20% thickest trees per hectare
+
 * `eq_biomass_*()`: equations to be used inside the `model` argument of `silv_predict_biomass()`.
 
 ## Bug Fixes

R/deprecated-funs.R

@@ -193,7 +193,7 @@ silv_dominant_height <- function(diameter,
         .cumtrees = cumsum(nt),
         .nmax     = which(.cumtrees >= 100)[1],
         .nmax     = if (is.na(.nmax[1])) which.max(.cumtrees) else .nmax,
-        .do       = calc_dominant_height(.nmax, nt, h)
+        .do       = calc_dominant_metric(.nmax, nt, h)
       ) |>
       dplyr::pull(.do)
   } else {
@@ -204,7 +204,7 @@ silv_dominant_height <- function(diameter,
         .cumtrees = cumsum(nt),
         .nmax     = which(.cumtrees >= 100)[1],
         .nmax     = if (is.na(.nmax[1])) which.max(.cumtrees) else .nmax,
-        .do       = calc_dominant_height(.nmax, nt, h)
+        .do       = calc_dominant_metric(.nmax, nt, h)
       ) |>
       dplyr::pull(.do)
   }

R/metrics-stand-level.R

@@ -142,32 +142,129 @@ silv_stand_dominant_height <- function(diameter,
     )
   }
 
-
   # 2. Calculate dominant height
   if (tolower(which) == "assman") {
-    d0 <- data |>
+    h0 <- data |>
       ## sort descending by diameter class
       dplyr::arrange(dplyr::desc(d)) |>
       dplyr::mutate(
         .cumtrees = cumsum(nt),
         .nmax     = which(.cumtrees >= 100)[1],
         .nmax     = if (is.na(.nmax[1])) which.max(.cumtrees) else .nmax,
-        .do       = calc_dominant_height(.nmax, nt, h)
+        .do       = calc_dominant_metric(.nmax, nt, h)
       ) |>
       dplyr::pull(.do)
   } else {
-    d0 <- data |>
+    h0 <- data |>
       ## sort descending by height
       dplyr::arrange(dplyr::desc(h)) |>
       dplyr::mutate(
         .cumtrees = cumsum(nt),
         .nmax     = which(.cumtrees >= 100)[1],
         .nmax     = if (is.na(.nmax[1])) which.max(.cumtrees) else .nmax,
-        .do       = calc_dominant_height(.nmax, nt, h)
+        .do       = calc_dominant_metric(.nmax, nt, h)
       ) |>
       dplyr::pull(.do)
   }
 
+  # 3. If it's not vectorized, retrieve just one value
+  if (is.null(ntrees)) h0[1] else h0
+
+}
+
+
+
+
+
+#' Calculates the dominant diameter
+#'
+#' Calculates the dominant diameter using Assman and Friedrich method, or
+#' Weise method
+#'
+#' @param diameter Numeric vector with diameter classes
+#' @param ntrees Optional. Numeric vector with number of trees per hectare.
+#' Use this argument when you have aggregated data by diametric classes (see details).
+#' @param which The method to calculate the dominant diameter (see details)
+#' @param quiet if \code{TRUE}, messages will be supressed
+#'
+#' @details
+#' The dominant diameter \eqn{D_0} is the mean diameter of the 100 thickest trees per
+#' hectare. Therefore, `diameter` and `ntrees` should be vectors of the same length.
+#' 
+#' - \bold{Assman}: calculates the \eqn{D_0} as the mean diameter of the 100 thickest
+#' trees per hectare
+#'
+#' - \bold{Weise}: calculates the \eqn{D_0} as the quadratic mean diameter of the
+#' 20% thickest trees per hectare
+#'
+#' @return A numeric vector
+#' @export
+#' @include utils-not-exported.R
+#'
+#' @examples
+#' ## calculate d0 for inventory data grouped by plot_id and species
+#' library(dplyr)
+#' inventory_samples |>
+#' mutate(dclass = silv_tree_dclass(diameter)) |>
+#'   summarise(
+#'     height = mean(height, na.rm = TRUE),
+#'     ntrees = n(),
+#'     .by    = c(plot_id, species, dclass)
+#'   ) |>
+#'   mutate(
+#'     ntrees_ha = silv_density_ntrees_ha(ntrees, plot_size = 10),
+#'     d0        = silv_stand_dominant_diameter(dclass, ntrees_ha),
+#'     .by       = c(plot_id, species)
+#'   )
+silv_stand_dominant_diameter <- function(diameter,
+                                        ntrees = NULL,
+                                        which = "assman",
+                                        quiet = FALSE) {
+
+  # 0. Handle errors and setup
+  ## 0.1. Errors
+  if (!tolower(which) %in% c("assman", "weise")) cli::cli_abort("`which` must be either <assman> or <weise>.")
+  if (!is.numeric(diameter)) cli::cli_abort("`diameter` must be a numeric vector")
+  ## 0.2. Invalid values
+  if (any(diameter <= 0, na.rm = TRUE)) cli::cli_warn("Any value in `diameter` is less than 0. Review your data.")
+
+  # 1. Create a data frame with input variables
+  if (is.null(ntrees)) {
+    data <- data.frame(
+      d  = diameter,
+      nt = 1
+    )
+  } else {
+    data <- data.frame(
+      d  = diameter,
+      nt = ntrees
+    )
+  }
+
+  if (tolower(which) == "assman") {
+    d0 <- data |> 
+        ## sort descending by diameter class
+        dplyr::arrange(dplyr::desc(d)) |>
+        dplyr::mutate(
+          .cumtrees = cumsum(nt),
+          .nmax     = which(.cumtrees >= 100)[1],
+          .nmax     = if (is.na(.nmax[1])) which.max(.cumtrees) else .nmax,
+          .do       = calc_dominant_metric(.nmax, nt, d)
+        ) |>
+        dplyr::pull(.do)
+  } else {
+    n_tickest_trees <- 0.2 * sum(data$nt)
+    d0 <- data |> 
+      ## sort descending by diameter class
+      dplyr::arrange(dplyr::desc(d)) |>
+      dplyr::mutate(
+        .cumtrees = cumsum(nt),
+        nt_sel = calc_accumulated_trees(nt, .cumtrees, n_tickest_trees),
+        .do = silv_stand_qmean_diameter(d, nt_sel)
+      ) |> 
+      dplyr::pull(.do)
+  }
+
   # 3. If it's not vectorized, retrieve just one value
   if (is.null(ntrees)) d0[1] else d0
 

R/utils-not-exported.R

@@ -1,15 +1,16 @@
 
 
-#' Calculates dominant height
+#' Calculates dominant height or dominant diameter
 #'
 #' @param nmax Index of first diametric class with > 100 trees; if there are no
 #'    100 trees, it is the index of the maximum
 #' @param ntress Number of trees per hectare
-#' @param height Height of the diametric class
+#' @param metric Height of the diametric class or diameter
+#' @param max_trees Number to trees to calculate dominant metric
 #'
 #' @return A numeric vector
 #' @keywords internal
-calc_dominant_height <- function(nmax, ntress, height) {
+calc_dominant_metric <- function(nmax, ntress, metric, max_trees = 100) {
 
   # initialize n and empty list
   n <- 0
@@ -20,14 +21,14 @@ calc_dominant_height <- function(nmax, ntress, height) {
     ## sum previous trees plus new trees
     n <- n + ntress[i]
     ## are we over 100 trees already?
-    if (n > 100) {
-      new_trees <- ntress[i] - n + 100
+    if (n > max_trees) {
+      new_trees <- ntress[i] - n + max_trees
       ## add to list and exit loop
-      l[[i]] <- c(new_trees, height[i])
+      l[[i]] <- c(new_trees, metric[i])
       break
     } else {
       ## add to list
-      l[[i]] <- c(ntress[i], height[i])
+      l[[i]] <- c(ntress[i], metric[i])
     }
   }
 
@@ -46,6 +47,39 @@ calc_dominant_height <- function(nmax, ntress, height) {
 
 
 
+#' Calculates number of trees until reaching a maximum number of trees
+#'
+#' @param ntrees Number of trees per hectare
+#' @param cumtrees Accumulated trees and sorted from thickest to thinner diameter
+#' @param max_trees Number to trees to calculate dominant metric
+#'
+#' @return A numeric vector
+#' @keywords internal
+calc_accumulated_trees <- function(ntrees, cumtrees, max_trees) {
+
+  ## row with with accumulated max_trees
+  row_with_target <- which(cumtrees >= max_trees)[1]
+
+  ## calculate trees needed from each diameter class
+  trees_needed <- rep(0, length(ntrees))
+
+  ## for rows before the target row, take all trees
+  if (row_with_target > 1) {
+    trees_needed[1:(row_with_target-1)] <- ntrees[1:(row_with_target-1)]
+  }
+
+  ## for the target row, calculate remaining trees needed
+  trees_from_previous <- ifelse(row_with_target == 1, 0, cumtrees[row_with_target-1])
+  trees_needed[row_with_target] <- max_trees - trees_from_previous
+
+  ## return
+  return(trees_needed)
+}
+
+
+
+
+
 #' Calculates weighted mean
 #'
 #' @param var An object containing the values whose weighted median is to be computed

R/zzz.R

@@ -15,6 +15,7 @@ utils::globalVariables(
     "h0",
     "n",
     "nt",
+    "nt_sel",
     "ntrees",
     "ntrees_ha",
     "remaining_to_extract",