add process_tree_data()

Author: ctoney

NAMESPACE

@@ -11,3 +11,4 @@ export(create_fia_owin)
 export(create_fia_ppp)
 export(load_tree_data)
 export(plot_crowns)
+export(process_tree_data)

R/calc_crwidth.R

@@ -19,13 +19,13 @@
 #' equation presented by Gill et al. (2000) to estimate crown width for nine
 #' tree species when their diameter is greater than 50 in. (127 cm).
 #'
-#' @param tree_list A data frame containing tree records. Must have columns
+#' @param tree_table A data frame containing tree records. Must have columns
 #' `SPCD` (FIA integer species code), `STATUSCD` (FIA integer tree status code,
 #' 1 = live) and `DIA` (FIA tree diameter in inches).
 #' @param digits Optional integer indicating the number of digits to keep in the
 #' return values (defaults to `1`).
 #' @return
-#' A numeric vector of length `nrow(tree_list)` with predicted crown width in
+#' A numeric vector of length `nrow(tree_table)` with predicted crown width in
 #' feet for live trees. `NA` is returned for trees with `STATUSCD != 1`.
 #'
 #' @references
@@ -47,87 +47,87 @@
 #' @examples
 #' calc_crwidth(plantation)
 #' @export
-calc_crwidth <- function(tree_list, digits = 1) {
-    if (missing(tree_list) || is.null(tree_list))
-        stop("'tree_list' is required", call. = FALSE)
+calc_crwidth <- function(tree_table, digits = 1) {
+    if (missing(tree_table) || is.null(tree_table))
+        stop("'tree_table' is required", call. = FALSE)
 
-    if (!is.data.frame(tree_list))
-        stop("'tree_list' must be a data frame", call. = FALSE)
+    if (!is.data.frame(tree_table))
+        stop("'tree_table' must be a data frame", call. = FALSE)
 
     required_cols <- c("SPCD", "STATUSCD", "DIA")
-    if (!all(required_cols %in% colnames(tree_list)))
-        stop("'tree_list' is missing required columns", call. = FALSE)
+    if (!all(required_cols %in% colnames(tree_table)))
+        stop("'tree_table' is missing required columns", call. = FALSE)
 
-    if (!is.numeric(tree_list$SPCD))
-        stop("'tree_list$SPCD' must be numeric or integer", call. = FALSE)
-    if (any(is.na(tree_list$SPCD)))
-        stop("'tree_list$SPCD' cannot have missing values", call. = FALSE)
+    if (!is.numeric(tree_table$SPCD))
+        stop("'SPCD' must be numeric or integer", call. = FALSE)
+    if (any(is.na(tree_table$SPCD)))
+        stop("'SPCD' cannot have missing values", call. = FALSE)
 
-    if (!is.numeric(tree_list$STATUSCD))
-        stop("'tree_list$STATUSCD' must be numeric or integer", call. = FALSE)
-    if (any(is.na(tree_list$STATUSCD)))
-        stop("'tree_list$STATUSCD' cannot have missing values", call. = FALSE)
+    if (!is.numeric(tree_table$STATUSCD))
+        stop("'STATUSCD' must be numeric or integer", call. = FALSE)
+    if (any(is.na(tree_table$STATUSCD)))
+        stop("'STATUSCD' cannot have missing values", call. = FALSE)
 
-    if (!is.numeric(tree_list$DIA))
-        stop("'tree_list$DIA' must be numeric", call. = FALSE)
-    if (any(is.na(tree_list$DIA)))
-        stop("'tree_list$DIA' cannot have missing values", call. = FALSE)
+    if (!is.numeric(tree_table$DIA))
+        stop("'DIA' must be numeric", call. = FALSE)
+    if (any(is.na(tree_table$DIA[tree_table$STATUSCD == 1])))
+        stop("'DIA' has missing values for live trees", call. = FALSE)
 
     if (is.null(digits))
         digits <- 1
 
-    cw <- rep_len(NA_real_, nrow(tree_list))
+    cw <- rep_len(NA_real_, nrow(tree_table))
 
     # define a default equation to use in case a species-specific one is missing
     # SPCD == 807, blue oak
     b_default <- cw_coef[cw_coef$SPCD == 807, c("b0", "b1", "b2")]
 
     # special case for large trees of certain species in the PNW region:
     # use the "old growth" equation from Gill et al. (2000)
-    old_growth_trees <- tree_list$DIA > 50 & tree_list$STATUSCD == 1 &
-        tree_list$SPCD %in% c(11, 98, 108, 119, 122, 202, 242, 263, 264)
+    old_growth_trees <- tree_table$DIA > 50 & tree_table$STATUSCD == 1 &
+        tree_table$SPCD %in% c(11, 98, 108, 119, 122, 202, 242, 263, 264)
 
-    cw[old_growth_trees] <- 16.449 + 0.4067 * tree_list$DIA[old_growth_trees]
+    cw[old_growth_trees] <- 16.449 + 0.4067 * tree_table$DIA[old_growth_trees]
 
     # apply species-specific equations
     # NB: crwidth of trees with DIA < 5 in. (i.e. "saplings") is predicted for
     # DIA = 5 and then sapling crwidth adjustment factors are applied afterward
-    for (spcd in unique(tree_list$SPCD)) {
+    for (spcd in unique(tree_table$SPCD)) {
         b <- cw_coef[cw_coef$SPCD == spcd, c("b0", "b1", "b2")]
         if (nrow(b) == 0)
             b <- b_default
 
         this_subset <-
-            tree_list$SPCD == spcd & tree_list$STATUSCD == 1 & is.na(cw)
+            tree_table$SPCD == spcd & tree_table$STATUSCD == 1 & is.na(cw)
 
         cw[this_subset] <-
-            b$b0 + b$b1 * pmax(5, tree_list$DIA[this_subset]) +
-            b$b2 * pmax(5, tree_list$DIA[this_subset])^2
+            b$b0 + b$b1 * pmax(5, tree_table$DIA[this_subset]) +
+            b$b2 * pmax(5, tree_table$DIA[this_subset])^2
     }
 
     # apply sapling crown width adjustment factors
-    saplings <- tree_list$DIA < 5 & tree_list$STATUSCD == 1
-    sapling_spp <- unique(tree_list$SPCD[saplings])
+    saplings <- tree_table$DIA < 5 & tree_table$STATUSCD == 1
+    sapling_spp <- unique(tree_table$SPCD[saplings])
     # species-specific adjustment factors if any (based on Bragg 2001)
     spcd_adj <- intersect(sapling_spp, cw_sapling_adj$SPCD)
     for (spcd in spcd_adj) {
         rowid <- which(cw_sapling_adj$SPCD == spcd)
         # adjustment factors at 1, 2, 3, 4, 5 inches DIA:
         adj_factors <- c(as.numeric(cw_sapling_adj[rowid, 2:5]), 1)
-        this_subset <- saplings & tree_list$SPCD == spcd
+        this_subset <- saplings & tree_table$SPCD == spcd
         # interpolated adjustment factors for the actual sapling diameters:
-        n <- trunc(tree_list$DIA[this_subset])
-        cw_adj <- (tree_list$DIA[this_subset] - n) *
+        n <- trunc(tree_table$DIA[this_subset])
+        cw_adj <- (tree_table$DIA[this_subset] - n) *
                   (adj_factors[n + 1] - adj_factors[n]) + adj_factors[n]
         cw[this_subset] <- cw[this_subset] * cw_adj
     }
     # otherwise use avarage adjustment factors based on Bragg (2001) data
     # average adjustment factors at 1, 2, 3, 4, 5 inches DIA:
     adj_factors <- c(0.509, 0.644, 0.767, 0.885, 1.0)
-    this_subset <- saplings & !(tree_list$SPCD %in% cw_sapling_adj$SPCD)
+    this_subset <- saplings & !(tree_table$SPCD %in% cw_sapling_adj$SPCD)
     # interpolated adjustment factors for the actual sapling diameters:
-    n <- trunc(tree_list$DIA[this_subset])
-    cw_adj <- (tree_list$DIA[this_subset] - n) *
+    n <- trunc(tree_table$DIA[this_subset])
+    cw_adj <- (tree_table$DIA[this_subset] - n) *
               (adj_factors[n + 1] - adj_factors[n]) + adj_factors[n]
     cw[this_subset] <- cw[this_subset] * cw_adj
 

R/calc_tcc_metrics.R

@@ -97,6 +97,8 @@
 #' @param digits Optional integer indicating the number of digits to keep in the
 #' return values (defaults to `1`). May be passed to `calc_crwidth()` and
 #' `calc_ht_metrics()`.
+#' @param ... Optional arguments passed to `create_fia_ppp()` if
+#' `stem_map = TRUE`.
 #' @return
 #' If `full_output = TRUE`, a named list with element `model_tcc` containing
 #' the plot-level predicted tree canopy cover as percent (`0:100`), and
@@ -142,7 +144,7 @@
 #' calc_tcc_metrics(plantation, stem_map = FALSE, full_output = FALSE)
 #' @export
 calc_tcc_metrics <- function(tree_list, stem_map = TRUE, full_output = TRUE,
-                             digits = 1) {
+                             digits = 1, ...) {
 
     if (!(is.logical(stem_map) && length(stem_map) == 1))
         stop("'stem_map' must be a single logical value", call. = FALSE)
@@ -155,7 +157,7 @@ calc_tcc_metrics <- function(tree_list, stem_map = TRUE, full_output = TRUE,
         # validate the input tree list for stem-mapping and get an estimate
         # of the L-function (square root transform of Ripley's K)
         # r = 0:12 feet
-        L <- create_fia_ppp(tree_list) |>
+        L <- create_fia_ppp(tree_list, ...) |>
             spatstat.explore::Lest(r = 0:12, correction = "isotropic")
 
         # mean of L at r = 6, 8, 10, 12 ft (Ripley's isotropic edge corrected)

R/process_tree_data.R

@@ -0,0 +1,115 @@
+#' Generate plot-level stand structure metrics for a tree dataset
+#'
+#' `process_tree_data()` takes a table of tree records for a set of forest
+#' inventory plots as input, and generates selected plot-level stand structure
+#' metrics.
+#'
+#' @param tree_table A data frame containing tree records for a set of forest
+#' inventory plots. Must have column `PLT_CN` containing the plot unique
+#' identifier for each tree. Other required columns are those of
+#' `calc_crwidth()` (if column `CRWIDTH` is not included), `calc_ht_metrics()`
+#' and `calc_tcc_metrics()`, depending on values given for `stem_map` and
+#' `full_output`.
+#' @param stem_map A logical value indicating whether to map individual tree
+#' stems explicitly, using coordinates specified in terms of distance and
+#' azimuth from subplot/microplot centers. The default is `TRUE`, in which case
+#' the input `tree_table` must contain columns `"DIST"` and `"AZIMUTH"`. This
+#' argument may be set to `FALSE` if individual tree locations are not
+#' available, in which case TCC will be predicted assuming a random arrangement
+#' of tree locations (see Details for [calc_tcc_metrics()]).
+#' @param full_output A logical value indicating whether to include the full set
+#' of components used to derive the plot-level TCC prediction. By default, the
+#' output data includes subplot-level TCC estimates, live tree and sapling
+#' counts, stand height metrics, and point pattern statistics, depending on the
+#' value given for `stem_map` (see Details for [calc_tcc_metrics()]).
+#' @param digits Optional integer indicating the number of digits to keep in the
+#' return values (defaults to `1`). May be passed to `calc_crwidth()` and
+#' `calc_ht_metrics()`.
+#' @return
+#' A data frame with one row for each unique `PLT_CN` in the input `tree_table`,
+#' and additional columns containing the output of `calc_tcc_metrics()`
+#' conditional on the values given for `stem_map` and `full_output`.
+#'
+#' @seealso
+#' [calc_ht_metrics()], [calc_tcc_metrics()], [load_tree_data()]
+#'
+#' @examples
+#' # Lolo NF, single-condition forest plots, INVYR 2022, from public FIADB
+#' f <- system.file("extdata/mt_lnf_2022_1cond_tree.csv", package="FIAstemmap")
+#' tree_table <- load_tree_data(f)
+#'
+#' process_tree_data(tree_table, stem_map = FALSE, full_output = TRUE)
+#' @export
+process_tree_data <- function(tree_table, stem_map = TRUE, full_output = TRUE,
+                              digits = 1) {
+
+    if (missing(tree_table) || is.null(tree_table))
+        stop("'tree_table' is required", call. = FALSE)
+
+    if (!is.data.frame(tree_table))
+        stop("'tree_table' must be a data frame", call. = FALSE)
+
+    if (!("PLT_CN" %in% colnames(tree_table))) {
+        stop("'tree_table' must have column 'PLT_CN' with unique plot IDs",
+             call. = FALSE)
+    }
+
+    if (!(is.logical(stem_map) && length(stem_map) == 1))
+        stop("'stem_map' must be a single logical value", call. = FALSE)
+
+    if (!(is.logical(full_output) && length(full_output) == 1))
+        stop("'full_output' must be a single logical value", call. = FALSE)
+
+    if (is.null(digits))
+        digits <- 1
+
+    plot_id_dt <- storage.mode(tree_table$PLT_CN)
+    if (!(plot_id_dt %in% c("character", "numeric", "integer", "integer64"))) {
+        stop("'PLT_CN' must be character, numeric, integer or integer64",
+             call. = FALSE)
+    }
+    plot_ids <- unique(tree_table$PLT_CN)
+    num_plots <- length(plot_ids)
+
+    # avoid creating a new owin on every call, pass through to create_fia_ppp()
+    w <- create_fia_owin()
+
+    # get the output for one plot, this validates input columns and defines
+    # the output data structure
+    tree_list <- tree_table[tree_table$PLT_CN == plot_ids[1], ]
+    x <- calc_tcc_metrics(tree_list, stem_map, full_output, digits, window = w)
+
+    out <- vector("list", 1 + length(x))
+    names(out) <- c("PLT_CN", names(x))
+
+    if (plot_id_dt == "character")
+        out$PLT_CN <- character(num_plots)
+    else if (plot_id_dt == "numeric")
+        out$PLT_CN <- rep(NA_real_, num_plots)
+    else if (plot_id_dt == "integer")
+        out$PLT_CN <- rep(NA_integer_, num_plots)
+    else if (plot_id_dt == "integer64")
+        out$PLT_CN <- rep(bit64::NA_integer64_, num_plots)
+
+    for (j in 2:length(out)) {
+        if (storage.mode(x[[j - 1]]) == "integer")
+            out[[j]] <- rep(NA_integer_, num_plots)
+        else
+            out[[j]] <- rep(NA_real_, num_plots)
+    }
+
+    cli::cli_progress_bar("Processing tree data", total = num_plots)
+    for (i in seq_along(plot_ids)) {
+        tree_list <- tree_table[tree_table$PLT_CN == plot_ids[i], ]
+        x <- calc_tcc_metrics(tree_list, stem_map, full_output, digits,
+                              window = w)
+        out$PLT_CN[i] <- plot_ids[i]
+        for (j in 2:length(out)) {
+            out[[j]][i] <- x[[j - 1]]
+        }
+        cli::cli_progress_update()
+    }
+    cli::cli_progress_done()
+
+    as.data.frame(out)
+}

README.Rmd

@@ -129,15 +129,14 @@ calc_tcc_metrics(plantation, stem_map = FALSE, full_output = FALSE)
 ### Data processing
 
 ```{r data-proc}
-## load tree data from a file or database connection
+# load tree data from a file or database connection
+# Lolo NF, single-condition forest plots, INVYR 2022, from public FIADB
 f <- system.file("extdata/mt_lnf_2022_1cond_tree.csv", package="FIAstemmap")
-tree <- load_tree_data(f)
+tree_table <- load_tree_data(f)
 
-head(tree)
+head(tree_table)
 
-## process tree data
-
-# TODO...
+process_tree_data(tree_table, stem_map = FALSE, full_output = TRUE)
 ```
 
 ## References

README.md

@@ -294,16 +294,17 @@ calc_tcc_metrics(plantation, stem_map = FALSE, full_output = FALSE)
 ### Data processing
 
 ``` r
-## load tree data from a file or database connection
+# load tree data from a file or database connection
+# Lolo NF, single-condition forest plots, INVYR 2022, from public FIADB
 f <- system.file("extdata/mt_lnf_2022_1cond_tree.csv", package="FIAstemmap")
-tree <- load_tree_data(f)
+tree_table <- load_tree_data(f)
 #> ! The data source does not have DIST and/or AZIMUTH
 #> ℹ Fetching tree data...
-#> ✔ Fetching tree data... [15ms]
+#> ✔ Fetching tree data... [14ms]
 #> 
 #> ℹ 910 tree records returned
 
-head(tree)
+head(tree_table)
 #>            PLT_CN SUBP TREE STATUSCD SPCD DIA HT ACTUALHT CCLCD TPA_UNADJ
 #> 1 670951075126144    1    1        2  108  NA NA       NA    NA        NA
 #> 2 670951075126144    1    2        1  108   1  9        9     3  74.96528
@@ -312,9 +313,53 @@ head(tree)
 #> 5 670951075126144    2    3        2  108  NA NA       NA    NA        NA
 #> 6 670951075126144    2    4        2  108  NA NA       NA    NA        NA
 
-## process tree data
-
-# TODO...
+process_tree_data(tree_table, stem_map = FALSE, full_output = TRUE)
+#>             PLT_CN model_tcc numTrees meanTreeHt meanTreeHtBAW meanTreeHtDom
+#> 1  670951075126144       1.2        0        0.0           0.0           0.0
+#> 2  670950940126144      38.4       24       61.4          66.4          64.5
+#> 3  670950992126144       3.4        1       43.0          43.0          43.0
+#> 4  670950609126144      17.2        4      102.2         102.6         102.2
+#> 5  670950600126144      34.9       16       62.1          79.0          69.9
+#> 6  670951118126144      20.6        9       24.6          28.0          30.0
+#> 7  670950964126144      37.9       16       58.8          67.1          64.1
+#> 8  670951031126144      51.2       29       70.2          72.8          72.4
+#> 9  670950608126144      70.8       32       73.7          94.8          86.6
+#> 10 670950599126144      66.4       44       61.8          66.4          64.1
+#> 11 670950967126144      57.4       23       86.0         100.7          96.1
+#> 12 670950732126144      34.4       12       64.3          91.8          72.6
+#> 13 670950725126144      66.5       69       66.3          87.0          73.9
+#> 14 670950598126144      55.8       20       65.7          89.6          83.9
+#> 15 670950965126144      81.3       74       53.1          55.1          54.5
+#> 16 670951032126144      32.5        5       15.0          14.2          15.0
+#> 17 670951034126144      16.4        7       40.7          45.0          40.7
+#> 18 670950625126144      44.5       23       42.0          61.9          42.9
+#> 19 670951029126144      55.1       33       64.7          68.5          64.7
+#> 20 670951035126144      97.6       54       44.9          50.6          45.7
+#> 21 670951089126144      21.1        7       79.9          83.0          79.9
+#> 22 670951152126144       5.3        3       21.3          21.7          21.3
+#>    meanTreeHtDomBAW maxTreeHt predomTreeHt numSaplings meanSapHt maxSapHt
+#> 1               0.0         0          0.0           1       9.0        9
+#> 2              67.9        85         81.7           1      16.0       16
+#> 3              43.0        43         43.0           0       0.0        0
+#> 4             102.6       114        106.3           0       0.0        0
+#> 5              83.4       104         99.3           0       0.0        0
+#> 6              34.8        47         39.7           0       0.0        0
+#> 7              69.0        80         78.0           0       0.0        0
+#> 8              73.6        85         83.0           0       0.0        0
+#> 9              98.3       120        112.7          19      15.8       38
+#> 10             67.2        84         81.7           1      14.0       14
+#> 11            103.7       123        117.0           2      13.0       16
+#> 12             94.2       109         93.7           0       0.0        0
+#> 13             89.6       118        116.3           2      12.5       16
+#> 14             97.3       128        112.0           5      11.4       18
+#> 15             56.1        72         67.0           3      40.3       45
+#> 16             14.2        22         18.3          15      12.3       20
+#> 17             45.0        53         48.3           0       0.0        0
+#> 18             63.0       104         70.0           2      20.5       25
+#> 19             68.5        87         83.3           3      19.0       22
+#> 20             51.2        74         66.0          27      23.5       39
+#> 21             83.0        92         85.3           0       0.0        0
+#> 22             21.7        24         21.3           1      14.0       14
 ```
 
 ## References

_pkgdown.yml

@@ -21,6 +21,7 @@ reference:
   - spatstat_helpers
   - load_tree_data
   - plot_crowns
+  - process_tree_data
 
 - title: Data
 - contents: