add compute_ht_metrics()

Author: ctoney

NAMESPACE

@@ -3,6 +3,7 @@
 export(.get_fia_plot_geom)
 export(.get_tree_list_xy)
 export(DEFAULT_TREE_COLUMNS)
+export(compute_ht_metrics)
 export(create_fia_owin)
 export(create_fia_ppp)
 export(load_tree_data)

R/compute_ht_metrics.R

@@ -9,32 +9,38 @@
 #' * `meanTreeHt`: mean height of trees `>= 5.0` in. (`12.7` cm) diameter
 #' * `meanTreeHtBAW`: basal-area weighted mean height of trees `>= 5.0` in.
 #'   (`12.7` cm) diameter
-#' * `meanTreeHtDom`: mean height of canopy dominant/codominant trees
+#' * `meanTreeHtDom`: mean height of canopy dominant/co-dominant trees
 #'   `>= 5.0` in. (`12.7` cm) diameter
 #' * `meanTreeHtDomBAW`: basal-area weighted mean height of canopy
-#'   dominant/codominant trees `>= 5.0` in. (`12.7` cm) diameter
+#'   dominant/co-dominant trees `>= 5.0` in. (`12.7` cm) diameter
 #' * `maxTreeHt`: height of the tallest tree `>= 5.0` in. (`12.7` cm) diameter
 #' * `predomTreeHt`: predominant tree height, as the mean height of the tallest
-#'   trees `>= 5.0` in. (`12.7` cm) diameter comprising `16` trees per acre
-#'   (`39.5` trees per hectare)
+#'   trees `>= 5.0` in. (`12.7` cm) diameter comprising up to `16` trees per
+#'   acre (`39.5` trees per hectare)
 #' * `meanSapHt`: mean height of saplings (trees `>= 1.0` in. diameter but
 #'   `< 5.0` in. diameter, i.e., `>= 2.54` cm but `< 12.7` cm)
 #' * `maxSapHt`: height of the tallest sapling
 #'
-#' For the purpose of height calculations, canopy dominant/codominant include
-#' "open grown" trees, i.e., include trees with FIA crown class codes (`CCLCD`)
-#' of `1` (open grown), `2` (dominant) and `3` (codominant), but exclude trees
-#' with `CCLCD` of `4` (intermediate) and `5` (overtopped).
+#' For the purpose of height calculations, canopy dominant/co-dominant include
+#' "open grown" trees, i.e., include trees with FIA crown class code (`CCLCD`)
+#' of `1` (open grown), `2` (dominant) or `3` (co-dominant), but exclude trees
+#' with `CCLCD` of `4` (intermediate) or `5` (over-topped).
 #'
 #' @param tree_list A data frame with tree records for one FIA plot.  Must have
 #' columns `DIA` (tree diameter), `HT` (tree height), `ACTUALHT` (tree actual
-#' height, `ACTUALHT < HT` indicating a broken top), `CCLCD1` (FIA crown class
+#' height, `ACTUALHT < HT` indicating a broken top), `CCLCD` (FIA crown class
 #' code), `TPA_UNADJ` (trees per acre).
+#' @param digits Optional integer indicating the number of digits to keep in the
+#' return values (defaults to `1`).
 #' @return
-#' A named list of height metrics computed for the input tree list, as described
+#' A named list of computed height metrics for the input tree list, as described
 #' in Details.
 #'
-compute_ht_metrics <- function(tree_list) {
+#' @examples
+#' compute_ht_metrics(plantation)
+#'
+#' @export
+compute_ht_metrics <- function(tree_list, digits = 1) {
 
     if (missing(tree_list) || is.null(tree_list))
         stop("'tree_list' is required", call. = FALSE)
@@ -47,6 +53,9 @@ compute_ht_metrics <- function(tree_list) {
     if (!all(required_cols %in% colnames(tree_list)))
         stop("'tree_list' is missing required columns", call. = FALSE)
 
+    if (is.null(digits))
+        digits <- 1
+
     # separate "trees" vs "saplings" here
     if ("STATUSCD" %in% colnames(tree_list)) {
         trees_in <- tree_list[tree_list$STATUSCD == 1 & tree_list$DIA >= 5, ]
@@ -65,42 +74,52 @@ compute_ht_metrics <- function(tree_list) {
         warning("one or more sapling heights are missing", call. = FALSE)
 
     ht_metrics <- vector(mode = "list", length = 8)
+    ht_metrics[1:length(ht_metrics)] <- 0  # by definition
     names(ht_metrics) <- c("meanTreeHt", "meanTreeHtBAW", "meanTreeHtDom",
                            "meanTreeHtDomBAW", "maxTreeHt", "predomTreeHt",
                            "meanSapHt", "maxSapHt")
 
-    basal_area <- pi * (trees_in$DIA / 2)^2
-
-    ht_metrics$meanTreeHt <- mean(tree_ht, na.rm = TRUE)
-    ht_metrics$meanTreeHtBAW <- stats::weighted.mean(tree_ht, basal_area,
-                                                     na.rm = TRUE)
-
-    tree_ht_doms <- tree_ht[trees_in$CCLCD %in% c(1, 2, 3)]
-    ht_metrics$meanTreeHtDom <- mean(tree_ht_doms, na.rm = TRUE)
-
-    basal_area_doms <- basal_area[trees_in$CCLCD %in% c(1, 2, 3)]
-    ht_metrics$meanTreeHtDomBAW <-
-        stats::weighted.mean(tree_ht_doms, basal_area_doms, na.rm = TRUE)
-
-    ht_metrics$maxTreeHt <- max(tree_ht, na.rm = TRUE)
-
-    tree_ht_tpa <- data.frame(tree_ht, trees_in$TPA_UNADJ, check.names = FALSE)
-    tree_ht_tpa <- tree_ht_tpa[order(tree_ht, decreasing = TRUE), ]
-    tpa <- 0
-    sum_ht <- 0
-    n <- 0
-    for (i in seq_len(nrow(tree_ht_tpa))) {
-        sum_ht <- sum_ht + tree_ht_tpa[i, 1]
-        tpa <- tpa + tree_ht_tpa[i, 2]
-        n <- n + 1
-        if (tpa > 16)
-            break
+    if (nrow(trees_in) > 0) {
+        basal_area <- pi * (trees_in$DIA / 2)^2
+        ht_metrics$meanTreeHt <-
+            round(mean(tree_ht, na.rm = TRUE), digits = digits)
+
+        ht_metrics$meanTreeHtBAW <-
+            round(stats::weighted.mean(tree_ht, basal_area, na.rm = TRUE),
+                  digits = digits)
+
+        tree_ht_doms <- tree_ht[trees_in$CCLCD %in% c(1, 2, 3)]
+        basal_area_doms <- basal_area[trees_in$CCLCD %in% c(1, 2, 3)]
+        ht_metrics$meanTreeHtDom <-
+            round(mean(tree_ht_doms, na.rm = TRUE), digits = digits)
+
+        ht_metrics$meanTreeHtDomBAW <-
+            round(stats::weighted.mean(tree_ht_doms, basal_area_doms,
+                                       na.rm = TRUE),
+                  digits = digits)
+
+        ht_metrics$maxTreeHt <- max(tree_ht, na.rm = TRUE)
+
+        tree_ht_tpa <-
+            data.frame(tree_ht, trees_in$TPA_UNADJ, check.names = FALSE)
+        tree_ht_tpa <- tree_ht_tpa[order(tree_ht, decreasing = TRUE), ]
+        tpa <- sum_ht <- n <- 0
+        for (i in seq_len(nrow(tree_ht_tpa))) {
+            sum_ht <- sum_ht + tree_ht_tpa[i, 1]
+            tpa <- tpa + tree_ht_tpa[i, 2]
+            n <- n + 1
+            if (tpa > 16)
+                break
+        }
+        ht_metrics$predomTreeHt <- round(sum_ht / n, digits = digits)
     }
 
+    if (nrow(saplings_in) > 0) {
+        ht_metrics$meanSapHt <-
+            round(mean(sapling_ht, na.rm = TRUE), digits = digits)
 
-
-
-
+        ht_metrics$maxSapHt <- max(sapling_ht, na.rm = TRUE)
+    }
 
     return(ht_metrics)
 }

_pkgdown.yml

@@ -14,6 +14,7 @@ home:
 reference:
 - title: Functions
 - contents:
+  - compute_ht_metrics
   - load_tree_data
   - overlay_crowns
   - plot_crowns

man/compute_ht_metrics.Rd

@@ -0,0 +1,15 @@
+test_that("compute_ht_metrics works", {
+    expected <- vector(mode = "list", length = 8)
+    names(expected) <- c("meanTreeHt", "meanTreeHtBAW", "meanTreeHtDom",
+                         "meanTreeHtDomBAW", "maxTreeHt", "predomTreeHt",
+                         "meanSapHt", "maxSapHt")
+    expected$meanTreeHt <- 45
+    expected$meanTreeHtBAW <- 45.4
+    expected$meanTreeHtDom <- 45
+    expected$meanTreeHtDomBAW <- 45.4
+    expected$maxTreeHt <- 51
+    expected$predomTreeHt <- 50.3
+    expected$meanSapHt <- 33.5
+    expected$maxSapHt <- 42
+    expect_equal(compute_ht_metrics(plantation), expected, tolerance = 0.1)
+})