Add tests comparing modeled TCC values with output from the v. 1.12 legacy code

DESCRIPTION

@@ -39,4 +39,4 @@ Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
 Language: en-US
-Config/roxygen2/version: 8.0.0
+Config/roxygen2/version: 8.1.0

R/calc_tcc_metrics.R

@@ -167,11 +167,11 @@ calc_tcc_metrics <- function(tree_list, stem_map = TRUE, full_output = TRUE,
     ht_metrics <- NULL
     if (stem_map || full_output) {
         # validate the input tree list for stand height calc and get metrics
-        ht_metrics <- calc_ht_metrics(tree_list)
+        ht_metrics <- calc_ht_metrics(tree_list, digits)
     }
 
     if (!("CRWIDTH" %in% colnames(tree_list)))
-        tree_list$CRWIDTH <- calc_crwidth(tree_list)
+        tree_list$CRWIDTH <- calc_crwidth(tree_list, digits)
 
     model_tcc <- NA_real_
     if (stem_map) {
@@ -189,7 +189,7 @@ calc_tcc_metrics <- function(tree_list, stem_map = TRUE, full_output = TRUE,
             if (nrow(tree_subp) == 0) {
                 subp_overlay[i] <- 0
             } else {
-                subp_overlay[i] <- calc_crown_overlay(tree_subp, 24)
+                subp_overlay[i] <- calc_crown_overlay(tree_subp, 24, digits)
             }
 
             # saplings in the microplot
@@ -200,7 +200,7 @@ calc_tcc_metrics <- function(tree_list, stem_map = TRUE, full_output = TRUE,
             if (nrow(sap_micr) == 0) {
                 micr_overlay[i] <- 0
             } else {
-                micr_overlay[i] <- calc_crown_overlay(sap_micr, 6.8)
+                micr_overlay[i] <- calc_crown_overlay(sap_micr, 6.8, digits)
             }
         }
 
@@ -230,12 +230,12 @@ calc_tcc_metrics <- function(tree_list, stem_map = TRUE, full_output = TRUE,
                 # and 926 omitted)
                 # coefficients given in Table 2, Toney et al. (2009)
                 sapling_component <-
-                    -8.036 +
-                    0.211 * micr_overlay_tcc +
-                    0.552 * ht_metrics$numSaplings +
-                    -0.131 * ht_metrics$numTrees +
-                    4.367 * log(ht_metrics$meanTreeHtBAW) +
-                    0.222 * L_mean
+                    -8.036498 +
+                    0.211119 * micr_overlay_tcc +
+                    0.552399 * ht_metrics$numSaplings +
+                    4.366521 * log(ht_metrics$meanTreeHtBAW) +
+                    -0.130718 * ht_metrics$numTrees +
+                    0.222308 * L_mean
 
                 # constrain sapling adjustment >= 0
                 sapling_component <- max(c(sapling_component, 0))

inst/extdata/test-tree_data.csv

@@ -0,0 +1,145 @@
+PLT_CN,SUBP,TREE,AZIMUTH,DIST,STATUSCD,SPCD,SPGRPCD,DIA,HT,ACTUALHT,CCLCD,TPA_UNADJ
+1,4,117,298,49.7,1,211,16,191,251,251,2,0.999188
+1,2,110,45,37.5,1,211,16,157.6,234,234,3,0.999188
+1,1,106,220,48,1,211,16,156.9,281,281,2,0.999188
+1,1,102,103,48.3,1,211,16,147.2,290,258,2,0.999188
+1,3,124,102,27.6,1,211,16,146.4,309,309,3,0.999188
+1,2,113,351,53.6,1,211,16,144.5,300,280,3,0.999188
+1,3,125,118,56.3,1,211,16,129.7,280,239,3,0.999188
+1,3,121,9,58.9,1,211,16,129.7,305,305,2,0.999188
+1,3,129,317,14.7,1,211,16,111.8,239,239,3,0.999188
+1,3,126,150,40.6,1,211,16,101.9,271,271,3,0.999188
+1,3,128,283,11.4,1,211,16,79.6,210,178,4,0.999188
+1,1,107,260,51.3,2,211,16,75.8,209,13,,0.999188
+1,3,120,360,54.6,1,211,16,74,201,201,4,0.999188
+1,3,122,102,22,1,211,16,58.8,160,134,4,0.999188
+1,4,114,104,41.8,1,263,13,45.8,170,170,3,0.999188
+1,1,105,201,37.4,1,263,13,37.5,190,190,3,0.999188
+1,2,112,263,51.2,2,211,16,34.1,135,11,,0.999188
+1,4,119,269,50.3,1,263,13,33.3,185,94,3,0.999188
+1,4,118,254,20.9,1,263,13,32.2,177,177,3,0.999188
+1,1,103,184,22.9,1,263,13,25.6,130,96,4,0.999188
+1,4,115,111,20.5,2,263,13,21.9,115,25,,6.018046
+1,3,130,354,18.7,1,263,13,20.5,140,140,5,6.018046
+1,4,116,190,12.6,1,263,13,17.7,102,102,4,6.018046
+1,2,111,157,21.1,1,263,13,15.5,120,83,4,6.018046
+1,3,127,169,9.9,1,263,13,10.8,78,78,5,6.018046
+1,1,101,94,11.5,1,263,13,9.2,61,61,5,6.018046
+1,1,109,343,12.1,1,263,13,7.9,57,57,5,6.018046
+1,1,100,84,21.1,1,263,13,7.5,56,56,5,6.018046
+1,3,123,91,14.8,1,211,16,6.8,35,35,5,6.018046
+1,1,108,284,13.9,1,263,13,6.5,61,61,5,6.018046
+1,1,104,185,10.4,1,263,13,6.3,28,28,5,6.018046
+2,3,137,173,15.5,1,202,10,15.6,84,84,2,6.018046
+2,1,108,209,19.7,1,202,10,15.4,90,90,3,6.018046
+2,1,109,216,23,1,202,10,15.3,84,84,3,6.018046
+2,1,103,138,21.3,1,202,10,15.3,90,90,3,6.018046
+2,2,114,120,21.6,1,202,10,13.9,81,81,3,6.018046
+2,2,117,162,21.6,1,202,10,12.3,85,85,3,6.018046
+2,2,126,323,20.3,1,73,19,11.8,99,99,2,6.018046
+2,3,132,28,16.1,1,108,21,11.2,79,79,3,6.018046
+2,4,156,180,20.8,1,202,10,11.1,79,79,3,6.018046
+2,2,118,177,16.5,1,202,10,10.9,75,75,3,6.018046
+2,2,113,108,11.4,1,202,10,10.6,79,79,3,6.018046
+2,2,123,315,2.9,2,108,21,10.3,85,85,,6.018046
+2,3,131,14,20.2,1,202,10,9.8,74,74,3,6.018046
+2,4,149,360,16.4,1,202,10,9.7,75,75,3,6.018046
+2,3,134,18,5.3,1,108,21,9.6,79,79,3,6.018046
+2,2,116,147,6.9,1,202,10,9.4,73,73,3,6.018046
+2,1,101,50,22.5,1,202,10,9.4,67,67,3,6.018046
+2,4,155,166,20.1,1,73,19,9.1,86,86,3,6.018046
+2,3,139,198,22.5,1,73,19,8.9,80,80,3,6.018046
+2,2,111,46,20.3,2,202,10,8.9,59,17,,6.018046
+2,2,129,7,11.1,1,73,19,8.7,84,84,3,6.018046
+2,3,147,334,19.6,1,73,19,8.5,85,85,3,6.018046
+2,4,153,76,14.6,1,73,19,8.5,88,88,3,6.018046
+2,2,119,184,16.3,1,202,10,8.4,73,73,3,6.018046
+2,1,100,19,13.1,1,202,10,8.3,51,51,3,6.018046
+2,4,161,260,21,1,73,19,8.1,71,71,3,6.018046
+2,4,157,187,1.9,1,73,19,8.1,79,79,3,6.018046
+2,1,106,169,15.4,1,202,10,7.8,57,57,3,6.018046
+2,3,148,336,22,1,108,21,7.7,72,72,3,6.018046
+2,2,128,360,14.4,1,202,10,7.5,65,65,3,6.018046
+2,2,120,186,8.3,2,202,10,7.5,27,24,,6.018046
+2,2,124,319,8.4,1,202,10,7.4,60,60,3,6.018046
+2,4,158,231,20.5,1,73,19,7.2,77,77,3,6.018046
+2,2,127,340,10.3,2,202,10,7.1,28,26,,6.018046
+2,2,122,235,22,1,202,10,7,55,55,4,6.018046
+2,4,154,166,15.3,1,202,10,6.9,50,46,4,6.018046
+2,3,146,316,15,1,108,21,6.9,77,77,3,6.018046
+2,4,150,13,8,1,73,19,6.8,68,68,3,6.018046
+2,2,121,228,19,1,202,10,6.8,46,46,4,6.018046
+2,1,107,197,15.6,1,202,10,6.6,58,58,4,6.018046
+2,3,142,257,17.2,2,108,21,6.5,46,46,,6.018046
+2,2,115,128,16,2,202,10,6.5,33,24,,6.018046
+2,2,139,352,7.2,1,202,10,6.4,57,57,4,6.018046
+2,4,159,233,16.2,1,73,19,6.4,76,76,3,6.018046
+2,4,160,258,18.5,1,73,19,6.2,72,72,3,6.018046
+2,4,177,98,23,1,73,19,5.9,64,64,4,6.018046
+2,1,102,84,22,1,202,10,5.9,36,36,4,6.018046
+2,3,138,179,10.2,1,202,10,5.8,45,45,4,6.018046
+2,4,162,303,13,1,73,19,5.7,64,64,3,6.018046
+2,3,145,274,9.4,2,108,21,5.7,59,59,,6.018046
+2,1,105,169,5.8,1,202,10,5.6,30,30,4,6.018046
+2,4,180,310,13.5,1,242,22,5.6,28,28,4,6.018046
+2,4,179,38,23,1,73,19,5.6,58,58,3,6.018046
+2,4,178,185,18,1,242,22,5.5,25,25,4,6.018046
+2,3,141,243,15,2,108,21,5.5,67,67,,6.018046
+2,1,104,140,7.7,2,202,10,5.3,37,37,,6.018046
+2,4,151,58,5.9,1,73,19,5.3,66,66,3,6.018046
+2,4,152,58,15,2,375,47,5.1,48,25,,6.018046
+2,4,181,270,19.4,1,202,10,5.1,32,32,4,6.018046
+2,4,174,138,1,1,73,19,4.5,51,51,4,74.965282
+2,4,173,125,2,1,202,10,3.4,33,33,4,74.965282
+2,4,172,7,3.5,1,242,22,2.4,12,12,5,74.965282
+2,4,176,278,3,1,242,22,1.9,8,8,5,74.965282
+2,1,110,294,2.9,1,202,10,1.6,11,11,5,74.965282
+2,1,120,25,5.9,1,202,10,1.1,9,9,5,74.965282
+2,4,175,243,4,1,242,22,1,8,8,5,74.965282
+2,2,125,322,14.7,2,108,21,,,,,
+2,3,140,230,13.2,2,108,21,,,,,
+2,3,143,269,16.1,2,108,21,,,,,
+2,2,112,100,22.5,2,202,10,,,,,
+2,3,144,272,14.9,2,108,21,,,,,
+2,3,136,116,9.5,2,108,21,,,,,
+2,3,135,86,19.6,2,108,21,,,,,
+2,3,133,30,10.4,2,108,21,,,,,
+3,1,1,299,19.8,1,922,41,6.9,36,36,3,6.018046
+3,1,2,315,12.7,1,742,37,7.2,56,56,3,6.018046
+3,1,3,313,21.2,1,922,41,7.3,45,45,3,6.018046
+3,1,4,279,5.9,1,313,41,4.6,28,28,4,74.965282
+3,1,5,280,6.2,1,6748,55,1.7,26,26,4,74.965282
+3,1,6,280,6.2,1,6748,55,1.1,21,21,5,74.965282
+3,2,1,18,7,1,731,41,7.1,48,48,3,6.018046
+3,2,2,132,6,1,313,41,6.9,42,42,3,6.018046
+3,2,3,151,5.6,1,731,41,6.9,53,53,3,6.018046
+3,2,4,165,16.6,1,313,41,6.1,35,35,3,6.018046
+3,2,5,190,16.3,1,731,41,8,54,54,3,6.018046
+3,2,6,223,22.6,1,313,41,7,35,35,3,6.018046
+3,2,7,238,18,1,313,41,8.5,38,38,3,6.018046
+3,2,8,292,15.2,1,922,41,13.2,58,58,3,6.018046
+3,2,9,359,3.6,1,731,41,1.7,38,38,5,74.965282
+3,2,10,56,6.3,1,6748,55,1.3,20,20,5,74.965282
+3,2,11,215,6.8,2,731,41,2.4,21,21,,74.965282
+3,2,12,325,5.1,2,731,41,1.1,19,19,,74.965282
+3,3,1,28,22.8,1,731,41,5.5,40,40,3,6.018046
+3,3,2,58,3.7,1,731,41,2.3,33,33,3,74.965282
+3,3,3,59,5.5,1,731,41,2.1,28,28,3,74.965282
+3,3,4,73,4.7,1,972,41,4.6,34,34,3,74.965282
+3,3,5,291,4.4,1,731,41,2,28,28,3,74.965282
+3,4,1,32,19.8,1,922,41,16.4,72,72,3,6.018046
+3,4,2,139,22.5,1,731,41,6.5,58,58,3,6.018046
+3,4,3,156,23.9,1,922,41,8.6,60,60,3,6.018046
+3,4,4,188,21.8,1,731,41,6.7,48,48,3,6.018046
+3,4,5,211,22.7,1,922,41,13.2,61,61,3,6.018046
+3,4,6,196,17.2,1,922,41,11,60,60,3,6.018046
+3,4,10,94,6.3,2,6748,55,1.1,15,7,,74.965282
+3,4,7,324,15.8,1,313,41,5.9,24,24,4,6.018046
+3,4,8,94,6.3,1,6748,55,1.4,16,16,5,74.965282
+3,4,9,94,6.3,1,6748,55,1.9,25,25,5,74.965282
+3,4,11,97,4.5,1,731,41,2.3,20,20,5,74.965282
+3,4,12,247,5,1,731,41,4.1,31,31,4,74.965282
+3,4,13,284,6.6,1,6748,55,2.6,19,19,5,74.965282
+3,4,14,291,6.2,1,6748,55,2.5,20,20,5,74.965282
+3,4,15,291,6.2,1,6748,55,1.8,18,18,5,74.965282
+3,4,16,353,6.3,1,313,41,3.1,25,25,5,74.965282

tests/testthat/test-calc_crown_overlay.R

@@ -15,4 +15,15 @@ test_that("calc_crown_overlay is calculated correctly", {
                 round((pi * 2^2) / (pi * 10^2) * 100 / 2, 0)
     res <- calc_crown_overlay(tree_list, sample_radius = 10, digits = 0)
     expect_equal(res, expected)
+
+    # test tree list used in the original subp_crcov.c
+    trees <- c(117,22.1,7.71,72,19.4,10.0,66,17.3,8.94,258,13.2,17.28) |>
+        matrix(nrow = 4, ncol = 3, byrow = TRUE) |>
+        as.data.frame()
+    colnames(trees) <- c("AZIMUTH", "DIST", "CRWIDTH")
+    expect_equal(calc_crown_overlay(trees, 24, digits = 0), 20)
+
+    # one tree with half the area of a subplot
+    tree_list <- data.frame(DIST = 1, AZIMUTH = 0, CRWIDTH = 33.94113)
+    expect_equal(calc_crown_overlay(tree_list, 24), 50, tolerance = 1e-3)
 })

tests/testthat/test-calc_tcc_metrics.R

@@ -1,6 +1,6 @@
-test_that("calc_tcc_metrics works", {
-    # one sapling in each microplot
-    # each sapling crown has half the area of the microplot
+test_that("calc_tcc_metrics returns correct values", {
+    ### one sapling in each microplot, each sapling crown has half the area of
+    ## the microplot
     subp <- c(1, 2, 3, 4)
     tree <- rep(1, 4)
     azimuth <- c(0, 90, 180, 270)
@@ -27,4 +27,55 @@ test_that("calc_tcc_metrics works", {
     tcc_pred <- calc_tcc_metrics(tree_list, stem_map = FALSE,
                                  full_output = FALSE)
     expect_equal(tcc_pred, 39.3, tolerance = 1e-3)
+
+    ### test against results from the legacy v. 1.12 Python code
+    f <- system.file("extdata/test-tree_data.csv", package="FIAstemmap")
+    tree_tbl <- load_tree_data(f)
+
+    tree_list_1 <- tree_tbl[tree_tbl$PLT_CN == "1", ]
+    expect_equal(nrow(tree_list_1), 31)
+    expect_warning(
+        tcc_pred <- calc_tcc_metrics(tree_list_1),
+        "13 points were rejected as lying outside the specified window")
+    expect_true(is.list(tcc_pred))
+    expect_equal(round(tcc_pred$model_tcc), 95)
+    expect_equal(round(tcc_pred$subp1_crown_overlay), 78)
+    expect_equal(round(tcc_pred$subp2_crown_overlay), 75)
+    expect_equal(round(tcc_pred$subp3_crown_overlay), 100)
+    expect_equal(round(tcc_pred$subp4_crown_overlay), 71)
+    expect_equal(round(tcc_pred$micr1_crown_overlay), 0)
+    expect_equal(round(tcc_pred$micr2_crown_overlay), 0)
+    expect_equal(round(tcc_pred$micr3_crown_overlay), 0)
+    expect_equal(round(tcc_pred$micr4_crown_overlay), 0)
+
+    tree_list_2 <- tree_tbl[tree_tbl$PLT_CN == "2", ]
+    expect_equal(nrow(tree_list_2), 74)
+    expect_no_error(
+        tcc_pred <- calc_tcc_metrics(tree_list_2, digits = 3))
+    expect_true(is.list(tcc_pred))
+    # model_tcc skipped here due to rounding error difference
+    # expect_equal(round(tcc_pred$model_tcc), 72)
+    expect_equal(round(tcc_pred$subp1_crown_overlay), 50)
+    expect_equal(round(tcc_pred$subp2_crown_overlay), 63)
+    expect_equal(round(tcc_pred$subp3_crown_overlay), 42)
+    expect_equal(round(tcc_pred$subp4_crown_overlay), 64)
+    expect_equal(round(tcc_pred$micr1_crown_overlay), 33)
+    expect_equal(round(tcc_pred$micr2_crown_overlay), 0)
+    expect_equal(round(tcc_pred$micr3_crown_overlay), 0)
+    expect_equal(round(tcc_pred$micr4_crown_overlay), 88)
+
+    tree_list_3 <- tree_tbl[tree_tbl$PLT_CN == "3", ]
+    expect_equal(nrow(tree_list_3), 39)
+    expect_no_error(
+        tcc_pred <- calc_tcc_metrics(tree_list_3, digits = 3))
+    expect_true(is.list(tcc_pred))
+    expect_equal(round(tcc_pred$model_tcc), 70)
+    expect_equal(round(tcc_pred$subp1_crown_overlay), 14)
+    expect_equal(round(tcc_pred$subp2_crown_overlay), 68)
+    expect_equal(round(tcc_pred$subp3_crown_overlay), 6)
+    expect_equal(round(tcc_pred$subp4_crown_overlay), 57)
+    expect_equal(round(tcc_pred$micr1_crown_overlay), 53)
+    expect_equal(round(tcc_pred$micr2_crown_overlay), 46)
+    expect_equal(round(tcc_pred$micr3_crown_overlay), 84)
+    expect_equal(round(tcc_pred$micr4_crown_overlay), 99)
 })