# ============================================================================= # DATA PREPARATION # ============================================================================= # This script provides an example of processing the raw data to calculate # biomass fluxes at the crown-position × quadrat level. The BCI (Barro Colorado # Island) censuses 6 and 7 are used as sample data. # # WORKFLOW: # 1. Load example data (BCI, census 6 and 7) # 2. Apply taper correction to DBH measurements # 3. Calculate AGB using allometric equations # 4. Apply modulus transformation to correct unrealistic growth # 5. Classify trees by demographic status (alive, recruited, died) # 6. Calculate AWP and AWM at tree and quadrat levels # 7. Assign crown positions using spatial analysis # 8. Apply Kohyama correction # # OUTPUT: Quadrat-level productivity data with crown position assignments # ============================================================================= # ============================================================================= # R SESSION SETUP # ============================================================================= # Clear the environment to start fresh rm(list = ls()) # Configure R session display options for better numeric precision and output formatting options( digits = 15, # Display 15 significant digits for numeric values tibble.print_max = 10, # Maximum number of tibble rows to display tibble.print_min = 6, # Minimum number of tibble rows to display tibble.width = Inf, # No limit on tibble width (show all columns) pillar.sigfig = 10 # Number of significant figures for pillar package ) # ============================================================================= # LIBRARY LOADING # ============================================================================= # Load required libraries for data manipulation and spatial analysis library(tidyverse) # Core data manipulation and visualization (dplyr, ggplot2, tidyr, etc.) library(conflicted) # Handle function conflicts between packages # Resolve function name conflicts by specifying preferred packages conflict_prefer("lag", "dplyr") # Use dplyr::lag instead of stats::lag conflict_prefer("filter", "dplyr") # Use dplyr::filter instead of stats::filter # ============================================================================= # DATA ACQUISITION: BCI CENSUS DATA # ============================================================================= # Load BCI forest plot data from two censuses # Data source: https://datadryad.org/dataset/doi:10.15146/5xcp-0d46 # # Required files (download and place in working directory): # - bci.spptable.rdata: Species table with taxonomic information and wood density # - bci.stem/bci.stem6.rdata: Stem measurements from census 6 # - bci.stem/bci.stem7.rdata: Stem measurements from census 7 # # The analysis uses two consecutive censuses to calculate productivity changes # ============================================================================= # Store current working directory as the main path for file operations main_path <- getwd() # Load BCI species table containing taxonomic information and wood specific gravity (WSG) load(file.path(main_path, "bci.spptable.rdata")) # Load BCI stem data from census 6 (first census in analysis) load(file.path(main_path, "bci.stem", "bci.stem6.rdata")) # Load BCI stem data from census 7 (second census in analysis) load(file.path(main_path, "bci.stem", "bci.stem7.rdata")) # ============================================================================= # CENSUS DATA PREPARATION # ============================================================================= # Process census 6 data: convert to tibble and add census identifier bci.stem6 <- bci.stem6 %>% as_tibble() %>% mutate( census = 6 # Add census number for tracking ) # Process census 7 data: convert to tibble and add census identifier bci.stem7 <- bci.stem7 %>% as_tibble() %>% mutate( census = 7 # Add census number for tracking ) # ============================================================================= # SPECIES TABLE PREPARATION # ============================================================================= # Process species table to extract key taxonomic and wood density information sp_tab <- bci.spptable %>% as_tibble() %>% rename_all(tolower) %>% # Convert all column names to lowercase select(sp, family, genus, species, latin, wsg) %>% # Select relevant columns distinct() %>% # Remove duplicate entries rename(name = latin) # Rename latin name column for consistency # Note: BCI species table provides species names and wood density (wsg) values. # However, for our study, we standardized species names for all sites using # WorldFlora (v.1.14-3, Kindt 2020). We then assigned wood density values to # each tree based on standardized species identity using the Global Wood Density # Database (Chave et al. 2009; Zanne et al. 2009) with the getWoodDensity # function from the BIOMASS package (v.2.1.11, Réjou-Méchain et al. 2017). See # the main manuscript for further details. # References: # Kindt, R. (2020). WorldFlora: An R package for exact and fuzzy matching of # plant names against the World Flora Online taxonomic backbone data. Appl. # Plant Sci., 8, e11388. # Chave, J., Coomes, D., Jansen, S., Lewis, S.L., Swenson, N.G. & Zanne, A.E. # (2009). Towards a worldwide wood economics spectrum. Ecol Lett, 12, 351–66. # Zanne, A.E., Lopez-Gonzalez, G., Coomes, D.A., Ilic, J., Jansen, S., Lewis, # S.L., et al. (2009). Global wood density database. Data from: Towards a # worldwide wood economics spectrum. Dryad Digital Repository. # Réjou-Méchain, M., Tanguy, A., Piponiot, C., Chave, J. & Hérault, B. (2017). # biomass: an r package for estimating above-ground biomass and its uncertainty # in tropical forests. Methods Ecol. Evol., 8, 1163–1167. # ============================================================================= # TAXONOMIC CLASSIFICATION: LIFE FORMS TO EXCLUDE # ============================================================================= # Create lookup table for life forms that should be excluded from analysis. # Palms and strangling figs have different growth patterns and allometries # compared to typical trees, so they are filtered out tax <- data.frame( genus = c( "Ficus", "Attalea", "Astrocaryum", "Socratea", "Oenocarpus", "Bactris", "Elaeis", "Chamaedorea" ), tax = factor(c( "strangler", # Ficus classification "palm", # All palm classifications "palm", "palm", "palm", "palm", "palm", "palm" )) ) # ============================================================================= # TAPER CORRECTION FUNCTION # ============================================================================= # Correct DBH measurements for taper when measured at non-standard heights. # Equation derived from Cushman et al. 2021: https://doi.org/10.1111/btp.12994 # # Parameters: # dbh_cm: Diameter at breast height in centimeters # hom: Height of measurement in meters (standard is 1.3 m) # wsg: Wood specific gravity # # Returns: Corrected DBH as if measured at standard 1.3 m height taper_correction <- function(dbh_cm, hom, wsg) { # Replace missing hom values with standard breast height (1.3 m) hom[is.na(hom)] <- 1.3 # Calculate taper parameter b based on tree characteristics # Larger trees, higher measurement points, and denser wood have different taper patterns b <- 0.151 - 0.025 * log(dbh_cm) - 0.02 * log(hom) - 0.021 * log(wsg) # Apply taper correction to estimate DBH at standard height dbhc <- dbh_cm / (exp(-b * (hom - 1.3))) return(dbhc) } # ============================================================================= # BASE FGEO DATA CONSTRUCTION # ============================================================================= # Combine both census datasets and prepare for analysis # This creates the foundation dataset with corrected measurements base_fgeo_data <- bind_rows(bci.stem6, bci.stem7) %>% rename_all(tolower) %>% # Standardize column names to lowercase mutate( dbh = dbh / 10, # Convert DBH from mm to cm site = "BCI", # Add site identifier exactdate = lubridate::ymd(exactdate), # Convert date to standard format E = 0.05176398, # Environmental stress parameter for BCI year = lubridate::year(exactdate) # Extract year from date ) %>% rename(dbh_cm = dbh) %>% # Rename for clarity rename(quadrat_20m = quadrat) %>% # Rename quadrat column filter(!is.na(dbh_cm) & status == "A") %>% # Keep only alive stems with valid DBH arrange(census, year) %>% # Sort by census and year left_join(., sp_tab, by = "sp") %>% # Add species information mutate( dbhc_cm = taper_correction(dbh_cm, hom, wsg) # Apply taper correction to DBH ) # ============================================================================= # QUADRAT CENTROID CALCULATION # ============================================================================= # Calculate central coordinates for each 20x20m quadrat # These coordinates are used for spatial analyses # Helper function to find the closest value in a sequence # Used to snap coordinates to quadrat centers # Note: This function may not work for smaller quadrats or for quadrats with a # small number of trees displaced in a non-random pattern—always check results. find_closest_quadrat_center <- function(num, sequence) { differences <- abs(sequence - num) # Calculate distances to all sequence values closest <- sequence[which.min(differences)] # Find minimum distance return(closest) } # Create sequence of quadrat centers (every 20m from 10 to 1030) sequence <- seq(from = 10, to = 1030, by = 20) # Calculate quadrat centroids based on average tree coordinates within each quadrat coord <- base_fgeo_data %>% group_by(quadrat_20m) %>% reframe( gx = mean(gx, na.rm = TRUE), # Mean x-coordinate gy = mean(gy, na.rm = TRUE) # Mean y-coordinate ) %>% ungroup() %>% mutate( qx = sapply(gx, find_closest_quadrat_center, sequence), # Snap to nearest quadrat center x qy = sapply(gy, find_closest_quadrat_center, sequence) # Snap to nearest quadrat center y ) %>% select(quadrat_20m, qx, qy) # ============================================================================= # ABOVEGROUND BIOMASS (AGB) CALCULATION FUNCTION # ============================================================================= # Calculate aboveground biomass in kilograms using modified Chave et al. (2014) equation # Reference: Global Change Biology, https://doi.org/10.1111/gcb.12629 # Code sourced from: https://github.com/umr-amap/BIOMASS/blob/master/R/computeAGB.R # # Parameters: # dbh: Diameter at breast height (cm) # wd: Wood density (g/cm^3) # E: Environmental stress parameter (unitless) # # Returns: Aboveground biomass in kilograms calculate_agb_chave2014 <- function(dbh, wd, E) { # Modified Eq 7 from Chave et al. 2014 # This is a pantropical allometric equation that accounts for wood density and climate agb <- exp(-2.023977 - 0.89563505 * E + 0.92023559 * log(wd) + 2.79495823 * log(dbh) - 0.04606298 * (log(dbh)^2)) return(agb) # Returns biomass in kg } # ============================================================================= # FINALIZE BASE FGEO DATASET # ============================================================================= # Add spatial coordinates, filter out non-tree life forms, and calculate AGB base_fgeo_data <- base_fgeo_data %>% left_join(coord, by = "quadrat_20m") %>% # Add quadrat centroid coordinates left_join(., tax %>% select(genus, tax), by = "genus") %>% # Add taxonomic classification filter(!tax %in% c("palm", "strangler")) %>% # Exclude palms and strangling figs select(c( "quadrat_20m", "site", "census", "year", "treeid", "stemid", "tag", "stemtag", "sp", "gx", "gy", "dbhc_cm", "hom", "exactdate", "wsg", "name", "E", "qx", "qy" )) %>% arrange(census, year) %>% # Sort by census and year mutate( agb = calculate_agb_chave2014(dbhc_cm, wsg, E) / 1000 # Calculate AGB and convert to Mg (tons) ) # ============================================================================= # DBH GROWTH CORRECTION PARAMETERS # ============================================================================= # Define parameters for detecting and correcting unrealistic DBH growth rates # Define allowed range of DBH change (cm/year) # Values outside this range are considered unrealistic and will be corrected # Negative growth up to -0.1 cm/yr allows for small measurement errors # Maximum growth of 5 cm/yr is a reasonable upper limit Ddbh_range <- c(-0.1, 5) # Lambda parameter for modulus transformation # This controls the shape of the transformation applied to DBH changes lambda <- 0.5 # ============================================================================= # HELPER FUNCTION: GET MODE OF DATES # ============================================================================= # Calculate the most frequent date in a vector # Used to assign a consistent date to grouped observations get_mode_value <- function(x) { # Remove NA values from input vector x <- na.omit(x) # Find unique values unique_x <- unique(x) # Return the value that appears most frequently unique_x[which.max(tabulate(match(x, unique_x)))] } # ============================================================================= # CREATE BASE DATASET WITH CENSUS STANDARDIZATION # ============================================================================= # Standardize census labels and prepare data for wide-format conversion base_fgeo <- base_fgeo_data %>% # Standardize census labels (convert to census_1 and census_2) mutate( census = ifelse(census == min(unique(census)), "census_1", "census_2"), exactdate = ymd(exactdate) # Ensure date format is consistent ) %>% rename(wd = wsg) %>% # Rename wood specific gravity to wood density select( # Select relevant columns for subsequent analyses # Keep quadrat_20m for estimating mean recruitment dates site, quadrat_20m, name, treeid, stemid, census, year, exactdate, gx, gy, hom, dbhc_cm, agb, wd, E, ) # ============================================================================= # CONVERT TO LONG FORMAT: MEASUREMENTS # ============================================================================= # Reshape measurements (hom, dbhc_cm, agb, year) to long format # This creates one row per tree-census-metric combination base_fgeo_long <- base_fgeo %>% pivot_longer( cols = c("hom", "dbhc_cm", "agb", "year"), # Metrics to reshape names_to = "metric", # Column name for metric type values_to = "value" # Column name for metric value ) # ============================================================================= # CONVERT TO WIDE FORMAT: MEASUREMENTS # ============================================================================= # Reshape measurements to wide format for cross-census comparisons # Creates columns like dbhc_cm_census_1, dbhc_cm_census_2, etc. # This allows calculation of changes between census periods base_fgeo_wide_measurements <- base_fgeo_long %>% pivot_wider( id_cols = c( site, quadrat_20m, name, treeid, stemid, gx, gy, wd, E, # Tree identifiers and attributes ), names_from = c(metric, census), # Create columns from metric-census combinations values_from = value # Fill with measurement values ) # ============================================================================= # CONVERT TO LONG FORMAT: DATES # ============================================================================= # Reshape exactdate to long format # Note: Sometimes exact dates differ within trees because when date was NA, # the mean date was used. This is resolved by using tree-level dates. # Date variations are typically only 1-2 days and don't affect analyses. base_fgeo_long_date <- base_fgeo %>% pivot_longer( cols = "exactdate", # Date column to reshape names_to = "metric", # Column name for metric type values_to = "value" # Column name for date value ) # ============================================================================= # CONVERT TO WIDE FORMAT: DATES # ============================================================================= # Reshape dates to wide format for cross-census time interval calculations base_fgeo_wide_date <- base_fgeo_long_date %>% pivot_wider( id_cols = c( site, quadrat_20m, name, treeid, stemid, gx, gy, wd, E, # Tree identifiers and attributes ), names_from = c(metric, census), # Create columns from metric-census combinations values_from = value # Fill with date values ) # ============================================================================= # MERGE MEASUREMENTS AND DATES INTO SINGLE WIDE DATASET # ============================================================================= # Combine measurement and date wide-format datasets base_fgeo_wide <- base_fgeo_wide_measurements %>% left_join(., base_fgeo_wide_date, by = c( "site", "quadrat_20m", "name", "treeid", "stemid", "gx", "gy", "wd", "E" # Match by tree identifiers ) ) # ============================================================================= # CALCULATE TIME INTERVALS AND GROWTH METRICS # ============================================================================= # Compute time between censuses and DBH changes # Handle missing dates by using quadrat-level mode dates base_fgeo_wide <- base_fgeo_wide %>% group_by(quadrat_20m) %>% # For recruits and lost stems with missing dates, use the most common date # in the quadrat (mode) mutate( mode_date_census_1 = get_mode_value(exactdate_census_1), # Most common census 1 date in quadrat exactdate_census_1 = if_else(is.na(exactdate_census_1), mode_date_census_1, exactdate_census_1), mode_date_census_2 = get_mode_value(exactdate_census_2), # Most common census 2 date in quadrat exactdate_census_2 = if_else(is.na(exactdate_census_2), mode_date_census_2, exactdate_census_2) ) %>% ungroup() %>% mutate( # Calculate exact time interval between censuses in years dT = time_length(interval(exactdate_census_1, exactdate_census_2), unit = "years"), # Flag stems with change in height of measurement (HOM) # TRUE if HOM changed between censuses, indicating potential measurement issues dHOM = !is.na(hom_census_2 - hom_census_1) & (hom_census_2 - hom_census_1) != 0, # Calculate annualized DBH change (cm/year) dDbh = (dbhc_cm_census_2 - dbhc_cm_census_1) / dT ) # ============================================================================= # MODULUS TRANSFORMATION FOR DBH CORRECTIONS # ============================================================================= # Apply modulus transformation to correct unrealistic DBH growth rates # This method identifies and adjusts outlier growth rates while preserving # the sign and relative magnitude of growth # Assign size classes for modulus calculation # Size classes are based on initial DBH (census 1, or census 2 for recruits) base_fgeo_wide_size_std <- base_fgeo_wide %>% group_by(site) %>% mutate( # Create size classes size_std_1 = as.character(cut(dbhc_cm_census_1, breaks = c(0, 5, 10, 20, 30, 50, 500), include.lowest = TRUE, right = FALSE )), size_std_2 = as.character(cut(dbhc_cm_census_2, breaks = c(0, 5, 10, 20, 30, 50, 500), include.lowest = TRUE, right = FALSE )), # Use census 1 size class, or census 2 if recruited between censuses size_std_initial = ifelse(is.na(dbhc_cm_census_1), size_std_2, size_std_1) ) %>% select(-size_std_1, -size_std_2) %>% ungroup() # ============================================================================= # MODULUS TRANSFORMATION FUNCTION # ============================================================================= # Apply power transformation that preserves sign # Compresses the range of extreme values while maintaining direction of change # # Parameters: # d : Numeric value to transform # lambda : Exponent for transformation (typically 0.5 for square root-like behavior) # # Formula: sign(d) * |d|^lambda apply_modulus_power_transformation <- function(d, lambda) { sign(d) * abs(d)^lambda } # ============================================================================= # CALCULATE SIZE CLASS-SPECIFIC MODULUS STATISTICS # ============================================================================= # For each size class, calculate mean and SD of modulus-transformed DBH growth # Only use "valid" growth rates (within expected range, no HOM change) base_fgeo_grouped <- base_fgeo_wide_size_std %>% mutate( # Classify growth rates as valid or needing correction wtd = case_when( # Keep values: no HOM change AND within expected range !dHOM & dDbh > Ddbh_range[1] & dDbh < Ddbh_range[2] & !is.na(dDbh) ~ "keep_values", # Change values: HOM changed OR outside expected range (dHOM | dDbh <= Ddbh_range[1] | dDbh >= Ddbh_range[2]) & !is.na(dDbh) ~ "change_values" ) ) %>% # Process each size class separately group_by(size_std_initial) %>% nest() %>% mutate( # Calculate modulus transformation statistics for valid values only mod = map(data, ~ .x %>% filter(wtd == "keep_values") %>% # Only use valid growth rates mutate(mod = apply_modulus_power_transformation(dDbh, lambda = lambda)) %>% # Get mean and SD of modulus-transformed DBH change reframe( mu = mean(mod), # Mean of transformed values sd = sd(mod) # SD of transformed values )) ) %>% mutate( mu = map_dbl(mod, ~ .$mu), # Extract mean values sigma = map_dbl(mod, ~ .$sd) # Extract SD values ) %>% select(-mod) # ============================================================================= # MODULUS TRANSFORMATION HELPER FUNCTIONS # ============================================================================= # Functions for calculating corrected DBH growth using modulus transformation # Calculates the density function for modulus-transformed values # Used in integration to compute variance of corrected DBH values # # Parameters: # x : Value at which to evaluate density # lambda : Modulus transformation parameter # mu : Mean of modulus-transformed distribution # sigma : SD of modulus-transformed distribution modulus_density_function <- function(x, lambda, mu, sigma) { # Apply inverse modulus transformation and multiply by normal density apply_modulus_power_transformation(x, 1 / lambda) * dnorm(x, mu, sigma) } # Integrates modulus_density_function over infinite interval to calculate expected # corrected DBH value for outliers # # Parameters: # lambda : Modulus transformation parameter # mu : Mean of modulus-transformed distribution # sigma : SD of modulus-transformed distribution # # Returns: Expected DBH change for corrected values compute_expected_dbh_change <- function(lambda, mu, sigma) { result <- integrate(modulus_density_function, -Inf, Inf, lambda = lambda, mu = mu, sigma = sigma)$value return(result) } # ============================================================================= # APPLY DBH CORRECTIONS TO OUTLIERS # ============================================================================= # Replace outlier DBH growth rates with size class-specific corrected values base_fgeo_grouped <- base_fgeo_grouped %>% mutate( lambda = lambda, # Add lambda parameter to each size class # Calculate corrected DBH value (varD) for each size class varD = pmap_dbl(list(lambda, mu, sigma), compute_expected_dbh_change), # Apply corrections to the data data_dbh_corr = map2(data, varD, ~ .x %>% mutate(dDbh_corr = case_when( wtd == "change_values" ~ .y, # Replace outliers with corrected value TRUE ~ dDbh # Keep original value for valid growth ))) ) %>% select(-data) # Remove original nested data # Unnest the corrected data # Now each row has dDbh_corr with outliers replaced by size class-specific corrections base_fgeo_unnested <- base_fgeo_grouped %>% unnest(cols = c(data_dbh_corr, lambda, mu, sigma)) # ============================================================================= # MODULUS TRANSFORMATION FOR AGB CORRECTIONS # ============================================================================= # Apply similar correction approach to AGB changes # More complex than DBH because AGB is a nonlinear function of DBH # For valid growth, calculate AGB from corrected DBH # For outliers, use expected AGB based on integration over DBH distribution # # Parameters: # data : Nested data frame for one size class # lambda : Modulus transformation parameter # mu : Mean of modulus-transformed DBH distribution # sigma : SD of modulus-transformed DBH distribution correct_agb_using_integration <- function(data, lambda, mu, sigma) { data_cor_agb <- data %>% rowwise() %>% mutate( dAGB_corr = if (isTRUE(wtd == "change_values")) { # For outliers, calculate expected AGB using integration (calculate_expected_agb_change( Di = dbhc_cm_census_1, # Initial DBH lambda = lambda, mu_s_b = mu, sigma_s_b = sigma, wd = wd, # Wood density E = E # Environmental stress )) / 1000 # Convert kg to Mg } else { dAGB_corr_base # For valid values, use DBH-based AGB } ) %>% ungroup() return(data_cor_agb) } # Expected AGB change for corrected stems # Integrates AGB difference over truncated normal distribution of DBH changes # # Parameters: # Di : Initial DBH (cm) # lambda : Modulus transformation parameter # mu_s_b : Mean of modulus-transformed DBH distribution # sigma_s_b : SD of modulus-transformed DBH distribution # wd : Wood density (g/cm³) # E : Environmental stress parameter calculate_expected_agb_change <- function(Di, lambda, mu_s_b, sigma_s_b, wd, E) { # Define the integrand function integrand <- function(x) { # Probability density from truncated normal (truncated at -Di^lambda to prevent negative DBH) dens <- truncnorm::dtruncnorm(x, mean = mu_s_b, sd = sigma_s_b, a = -Di^lambda) # AGB change for this DBH change value AGB_diff <- calculate_agb_difference_from_dbh_change(x, Di, lambda, wd, E) # Weight AGB change by probability density AGB_diff * dens } # Integrate over all possible DBH changes (truncated at minimum) result <- integrate(integrand, lower = -Di^lambda, upper = Inf) return(result$value) } # compute_AGB_diff: Calculate AGB difference for a given DBH change # Applies inverse modulus transformation to get final DBH, then calculates # AGB difference using Chave equation # # Parameters: # x : Modulus-transformed DBH change # Di : Initial DBH (cm) # lambda : Modulus transformation parameter # wd : Wood density (g/cm³) # E : Environmental stress parameter calculate_agb_difference_from_dbh_change <- function(x, Di, lambda, wd, E) { # Back-transform to get final DBH Di_mod <- Di + apply_modulus_power_transformation(x, 1 / lambda) # Calculate AGB difference calculate_agb_chave2014(Di_mod, wd, E) - calculate_agb_chave2014(Di, wd, E) } # ============================================================================= # APPLY AGB CORRECTIONS # ============================================================================= # Calculate baseline corrected AGB (from corrected DBH) and handle recruitment base_fgeo_corrected_agb <- base_fgeo_unnested %>% mutate( # Calculate baseline AGB change from corrected DBH # Uses final DBH = initial DBH + corrected DBH growth dAGB_corr_base = (calculate_agb_chave2014(dbhc_cm_census_1 + dDbh, wd, E) - calculate_agb_chave2014(dbhc_cm_census_1, wd, E)) / 1000 ) %>% # Handle recruitment (stems that appear only in census 2) mutate(dAGB_corr_base = case_when( # For recruits: AGB change is the final AGB divided by time interval is.na(agb_census_1) & !is.na(agb_census_2) ~ agb_census_2 / dT, TRUE ~ dAGB_corr_base # Otherwise use calculated value )) %>% # Group by size class and transformation parameters group_by(size_std_initial, lambda, mu, sigma) %>% nest() %>% # Nest data for each size class mutate( # Apply AGB correction function to each size class data_corrected_agb = pmap(list(data, lambda, mu, sigma), correct_agb_using_integration) ) # ============================================================================= # FINALIZE CORRECTED DATASET # ============================================================================= # Unnest and clean the final dataset with corrected DBH and AGB values base_fgeo_final <- base_fgeo_corrected_agb %>% select(data_corrected_agb) %>% unnest(data_corrected_agb) %>% ungroup() %>% mutate(id = row_number()) %>% # Add unique row identifier # Round corrected AGB to 6 decimal places for consistency mutate(dAGB_corr = round(dAGB_corr, 6)) # ============================================================================= # DEMOGRAPHIC CLASSIFICATION # ============================================================================= # Classify stems and trees based on their presence/absence across two censuses # Categories: alive, recruited (ingrowth), or dead # Note: This step assumes all DBH data are correct, leading to an accurate # estimation of presence/absence data. It is important to use all available # census data, which allows the user to interpolate missing DBH values (e.g., # via linear interpolation) when individual stems were missed in the field. In # this example, we skip this step, as it can be performed by the user. base_fgeo_demo <- base_fgeo_final %>% # Categorize individual stems based on presence in censuses mutate( demo = case_when( # Stems present in both censuses are 'alive' !is.na(dbhc_cm_census_1) & !is.na(dbhc_cm_census_2) ~ "alive", # Stems absent in census 1 but present in census 2 are recruits ('ingrowth_stems') is.na(dbhc_cm_census_1) & !is.na(dbhc_cm_census_2) ~ "ingrowth_stems", # Stems present in census 1 but absent in census 2 have died !is.na(dbhc_cm_census_1) & is.na(dbhc_cm_census_2) ~ "died", # Other cases should not exist in properly filtered data TRUE ~ "dont_exist" ), # Create binary presence/absence indicators for stem census data # Used to determine tree-level demography c1 = ifelse(!is.na(dbhc_cm_census_1) & dbhc_cm_census_1 > 0, 1, 0), # Stem present in census 1 c2 = ifelse(!is.na(dbhc_cm_census_2) & dbhc_cm_census_2 > 0, 1, 0) # Stem present in census 2 ) %>% # Group by tree to determine tree-level demographic status # (trees can have multiple stems) group_by(treeid) %>% mutate( demo_tree = case_when( # Trees with at least one stem in census 1 but no stems in census 2 are dead any(c1 == 1) & all(c2 == 0) ~ "death_tree", # Trees with stems in both censuses are alive any(c1 == 1) & any(c2 == 1) ~ "alive_tree", # Trees with no stems in census 1 but stems in census 2 are recruits all(c1 == 0) & any(c2 == 1) ~ "recruited_tree", # Other cases should not exist TRUE ~ "dont_exist" ) ) %>% ungroup() # ============================================================================= # ABOVEGROUND WOOD PRODUCTIVITY: STANDING (ALIVE) STEMS # ============================================================================= # Calculate productivity for stems that were alive in both censuses # These stems contribute through diameter/biomass growth AWP_alive <- base_fgeo_demo %>% ungroup() %>% # Filter to stems alive in both censuses (no HOM changes or extreme growth already corrected) filter(demo == "alive" & demo_tree == "alive_tree") %>% mutate( status = "AWP_standing" # Label as standing productivity ) # ============================================================================= # ABOVEGROUND WOOD PRODUCTIVITY: INGROWTH STEMS # ============================================================================= # Calculate productivity for stems that recruited on existing trees # (new stems on trees that were already present in census 1) AWP_ingrowth_stems <- base_fgeo_demo %>% # Note: Use original quadrat census 1 date to estimate recruitment timing # For recruited stems, census 1 date is NA, so we use the quadrat's # mean census date as the assumed recruitment date filter(demo == "ingrowth_stems" & demo_tree == "alive_tree") %>% mutate( status = "AWP_ingrowth_stems", # Recalculate time interval using assigned census 1 date dT = time_length(interval(exactdate_census_1, exactdate_census_2), unit = "years"), # For recruited stems, productivity is final AGB divided by time interval dAGB_corr = agb_census_2 / dT ) %>% # Quality control: Remove ingrowth stems with unrealistic DBH # DBH should not exceed maximum possible growth over the interval filter((dbhc_cm_census_2) < (dT * Ddbh_range[2])) # ============================================================================= # ABOVEGROUND WOOD PRODUCTIVITY: RECRUITED TREES # ============================================================================= # Calculate productivity for entirely new trees (not present in census 1) AWP_ingrowth_tree <- base_fgeo_demo %>% # Filter to stems on newly recruited trees filter(demo == "ingrowth_stems" & demo_tree == "recruited_tree") %>% mutate( status = "AWP_ingrowth_trees", # Calculate time interval dT = time_length(interval(exactdate_census_1, exactdate_census_2), unit = "years"), # Productivity is final AGB divided by time interval dAGB_corr = agb_census_2 / dT ) %>% # Quality control: Remove recruited trees with unrealistic DBH filter((dbhc_cm_census_2) < (dT * Ddbh_range[2])) # ============================================================================= # COMBINE ALL PRODUCTIVITY COMPONENTS AT TREE LEVEL # ============================================================================= # Aggregate all productivity sources (standing, ingrowth stems, recruited trees) # Sum stem-level productivity to tree-level (trees can have multiple stems) AWP <- bind_rows(AWP_alive, AWP_ingrowth_stems, AWP_ingrowth_tree) %>% # Sum productivity across all stems of each tree group_by(site, quadrat_20m, name, treeid, gx, gy, wd, demo_tree) %>% # Note: sum with na.rm=TRUE converts NA to 0 reframe( agb_census_1 = sum(agb_census_1, na.rm = TRUE), # Initial AGB (Mg) dAGB_corr = sum(dAGB_corr, na.rm = TRUE), # Annual productivity (Mg/yr) dT = mean(dT, na.rm = TRUE) # Time interval (years) ) # ============================================================================= # ABOVEGROUND WOOD MORTALITY (AWM) # ============================================================================= # Calculate biomass loss from mortality # Includes both entire tree deaths and individual stem deaths on living trees AWM <- base_fgeo_demo %>% # Include dead stems from dead trees AND dead stems from trees that survived filter((demo == "died" & demo_tree == "death_tree") | (demo == "died" & demo_tree == "alive_tree")) %>% mutate( status = "AWM", # Calculate time interval dT = time_length(interval(exactdate_census_1, exactdate_census_2), unit = "years"), # Mortality rate is initial AGB divided by time interval dAGB_corr = agb_census_1 / dT ) %>% # Aggregate to tree level group_by(site, quadrat_20m, name, treeid, gx, gy, wd, demo_tree) %>% reframe( agb_census_1 = sum(agb_census_1, na.rm = TRUE), # Initial AGB (Mg) dAGB_corr = sum(dAGB_corr, na.rm = TRUE), # Annual mortality (Mg/yr) dT = mean(dT, na.rm = TRUE) # Time interval (years) ) # ============================================================================= # ALLOMETRIC FUNCTIONS: GENERAL FORMS # ============================================================================= # Three common allometric equation forms used for height and crown radius predictions # All include log-bias correction factor: exp(sigma^2 / 2) # Power Function (PF): Simple power law relationship # Parameters: # dbhc : Diameter at breast height (cm) # a : Scaling coefficient # b : Allometric exponent # sigma : Residual standard error (for bias correction) # log : Return log-transformed value if TRUE pf_function <- function(dbhc, a, b, sigma, log = FALSE) { if (missing(a) || missing(b) || missing(sigma)) { stop("Please provide values for 'a', 'b', and 'sigma'.") } if (log) { log(a) + (b * log(dbhc)) } else { (a * (dbhc^b)) * exp((sigma^2) / 2) # Bias-corrected power function } } # Generalized Michaelis-Menten (GMM): Asymptotic relationship with saturation # Parameters: # dbhc : Diameter at breast height (cm) # a : Maximum asymptotic value # b : Growth rate parameter # k : Half-saturation constant # sigma : Residual standard error (for bias correction) # log : Return log-transformed value if TRUE gmm_function <- function(dbhc, a, b, k, sigma, log = FALSE) { if (missing(a) || missing(b) || missing(k) || missing(sigma)) { stop("Please provide values for 'a', 'b', 'k', and 'sigma'.") } if (log) { log(a) + (b * log(dbhc)) - log(k + exp(b * log(dbhc))) } else { ((a * (dbhc^b)) / (k + (dbhc^b))) * exp((sigma^2) / 2) # Bias-corrected GMM } } # ============================================================================= # HEIGHT ALLOMETRY FOR BCI # ============================================================================= # Site-specific height-DBH allometry for BCI # Parameters from Martinez Cano et al. (2019) # Reference: https://doi.org/10.5194/bg-16-847-2019 height_bci <- function(dbhc) { gmm_function( dbhc = dbhc, a = 57.17, # Maximum height (m) b = 0.73, # Growth rate parameter k = 21.57, # Half-saturation constant sigma = 0.181, # Residual error log = FALSE ) } # ============================================================================= # CROWN AREA ALLOMETRY FOR BCI # ============================================================================= # Crown area-DBH allometry for BCI. Parameters estimated from ForestGEO Americas # sites and the Tallo database. Refer to the main study for details. See # `crown_allometry.R` for derivation details. crown_bci <- function(dbhc) { pf_function( dbhc = dbhc, a = 0.599, # Scaling coefficient b = 1.25, # Allometric exponent sigma = 0.765, # Residual error log = FALSE ) } # ============================================================================= # SPATIAL DATA PREPARATION: CREATE SF OBJECTS # ============================================================================= # Convert tree crowns to spatial features (SF) objects for spatial analyses # Uses predicted crown radii to create circular crown polygons # Calculate allometric dimensions for the initial census # Note: Only initial census to avoid double-counting; recruited trees will be # assigned to lowest crown position layer in each quadrat base_fgeo_allo <- base_fgeo_wide %>% mutate(dbhc = dbhc_cm_census_1) %>% # Filter to standing trees in census 1 only # Recruited trees (NA in census 1) will be added to understory layer later filter(!is.na(dbhc)) %>% # For multi-stemmed trees, use the largest stem for allometry group_by(site, name, treeid, gx, gy) %>% reframe(dbhc = max(dbhc, na.rm = TRUE)) %>% # Maximum DBH per tree mutate(heightm = height_bci(dbhc)) %>% # Predict height (m) mutate(caream2 = crown_bci(dbhc)) %>% # Predict crown area (m²) mutate(cradsm = sqrt(caream2 / pi)) # Calculate crown radius (m) # Split data into list for parallel processing # Each element is one tree fgeo_data_quad_allo <- base_fgeo_allo %>% mutate(group_var = paste(site, treeid, sep = "-")) %>% split(.$group_var) # ============================================================================= # PARALLEL PROCESSING SETUP FOR CROWN POLYGONS # ============================================================================= # Create circular crown polygons for each tree using parallel processing # This is computationally intensive for large datasets library(foreach) library(doSNOW) # Set up parallel backend num_cores <- parallel::detectCores(logical = TRUE) - 1 # Use all cores except one cl <- makeCluster(num_cores, type = "SOCK") registerDoSNOW(cl) # Define number of iterations (one per tree) iterations <- length(fgeo_data_quad_allo) # Set up progress bar for monitoring library(progress) pb <- progress_bar$new( format = "% = :perc [:bar] :elapsed | eta: :eta", total = iterations, width = 60 ) progress_perc <- round(seq(from = 0, to = 100, length.out = iterations), 2) progress <- function(n) { pb$tick(tokens = list(perc = progress_perc[n])) } opts <- list(progress = progress) # Initialize results list results <- vector("list", length = iterations) # ============================================================================= # HELPER FUNCTION: CREATE CIRCULAR CROWN POLYGONS # ============================================================================= # Generates an SF circular polygon for a tree crown # # Parameters: # radius : Crown radius (m) # centerPoint : Named vector with tree coordinates (x, y) # nptsPerimeter : Number of points to define circle (default 100) # spID : Unique identifier for the polygon create_spatial_circle <- function(radius, centerPoint = c(x = 0, y = 0), # Default center at origin nptsPerimeter = 100, spID = paste("circle", .StemEnv$randomID(), sep = ":"), ...) { # Validate centerPoint if (any(is.na(match(c("x", "y"), names(centerPoint))))) { stop("Please use names x and y for circle centerPoint vector") } if (length(centerPoint) != 2) { stop("Please supply one set of (x,y) coordinates for the plot center location.") } # Validate parameters if (radius <= 0) { stop("radius must be positive!") } if (nptsPerimeter < 20) { warning("Using less than 20 points for the circle perimeter is not recommended--set to 20") nptsPerimeter <- 20 } area <- pi * radius * radius # Generate circle points using parametric equations # x = center_x + r*cos(theta), y = center_y + r*sin(theta) circ <- seq(0, 2 * pi, len = nptsPerimeter) # Create circle outline coordinates circle <- matrix(c( centerPoint["x"] + radius * cos(circ), # X coordinates centerPoint["y"] + radius * sin(circ), # Y coordinates rep(1, nptsPerimeter) # Height component (unused) ), nrow = nptsPerimeter) # Close the polygon by adding first point at end # This ensures proper polygon creation in SF circle <- rbind(circle, circle[1, ]) dimnames(circle) <- list(NULL, c("x", "y", "hc")) # Create SF polygon sfCircle <- sf::st_sfc(sf::st_polygon(list(circle[, -3]))) return(sfCircle) } # ============================================================================= # CREATE CROWN POLYGONS IN PARALLEL # ============================================================================= # Process each tree to create crown polygon with progress reporting results <- foreach( i = 1:iterations, .packages = c("sf"), .options.snow = opts ) %dopar% { # Create circular crown for this tree sfcrown <- create_spatial_circle( radius = fgeo_data_quad_allo[[i]]$cradsm, # Crown radius (m) centerPoint = c(x = fgeo_data_quad_allo[[i]]$gx, y = fgeo_data_quad_allo[[i]]$gy), spID = fgeo_data_quad_allo[[i]]$treeid ) # Combine tree data with crown geometry sfcrownaream2 <- sf::st_sf(fgeo_data_quad_allo[[i]], geometry = sfcrown) sfcrownaream2 # Return SF object } # Stop parallel cluster stopCluster(cl) library(sf) # Combine all crown polygons into single SF dataframe fgeo_data_quad_allo_sf <- dplyr::bind_rows(results) # ============================================================================= # DEFINE PLOT BOUNDARIES AND BUFFERS # ============================================================================= # Create plot outline and buffer zones for edge effect handling # Determine plot dimensions from data plot_length <- ceiling(max(fgeo_data_quad_allo_sf$gx)) # West-East extent (m) plot_width <- ceiling(max(fgeo_data_quad_allo_sf$gy)) # South-North extent (m) # Buffer distance for edge effects # Matches quadrat size (20m) to handle crowns extending beyond plot edges buffer_distance <- 20 # Define plot corner coordinates in clockwise order # Polygon must close (first point = last point) corners <- matrix( c( 0, 0, # Bottom left (SW) plot_length, 0, # Bottom right (SE) plot_length, plot_width, # Top right (NE) 0, plot_width, # Top left (NW) 0, 0 # Close polygon ), byrow = TRUE, ncol = 2 ) # Create the plot polygon and define its coordinate reference system (CRS) plot_polygon <- st_polygon(list(corners)) plot_sf <- st_sfc(plot_polygon) # Temporary CRS definition st_crs(plot_sf) <- st_crs(fgeo_data_quad_allo_sf) # Set CRS from existing sf object # Create the inward buffer as a separate polygon (negative distance creates an # inward buffer) usable_area <- st_buffer(plot_sf, dist = -buffer_distance) # Calculate the actual usable area by subtracting the buffer from the original # plot area buffer <- st_difference(plot_sf, usable_area) # At this point, 'usable_area' represents the area within the plot excluding the # buffer 'buffer' represents the subtracted area around the plot edge defined by # 'buffer_distance' ## 30 BY 30 PPA # Define grid size grid_size <- 30 # Function to create a grid centered on the usable area create_centered_grid <- function(usable_area, grid_size) { # Calculate the centroid of the usable area centroid <- st_centroid(usable_area) # Extract the coordinates of the centroid centroid_coords <- st_coordinates(centroid) # Calculate half the width and height of the grid based on the grid size half_width <- grid_size * (floor((st_bbox(usable_area)$xmax - st_bbox(usable_area)$xmin) / grid_size) / 2) half_height <- grid_size * (floor((st_bbox(usable_area)$ymax - st_bbox(usable_area)$ymin) / grid_size) / 2) # Define the bounding box (bbox) for the grid centered around the centroid grid_bbox <- st_bbox(c( xmin = as.numeric(centroid_coords[1] - half_width), xmax = as.numeric(centroid_coords[1] + half_width), ymin = as.numeric(centroid_coords[2] - half_height), ymax = as.numeric(centroid_coords[2] + half_height) ), crs = st_crs(usable_area)) # Create a grid using the defined bbox grid <- st_make_grid(st_as_sfc(grid_bbox), cellsize = c(grid_size, grid_size), square = TRUE) # Convert the grid to a simple feature (sf) object grid_sf <- st_sf(geometry = st_sfc(grid, crs = st_crs(usable_area))) # Intersect the grid with the usable area to get only the cells within the usable area grid_intersected <- st_intersection(grid_sf, usable_area) # Calculate the centroids of each grid cell centroids <- st_centroid(grid_intersected) # Extract x and y coordinates of each centroid coords <- st_coordinates(centroids) # Add the centroid coordinates as columns to the grid data frame grid_intersected$X_center <- coords[, "X"] grid_intersected$Y_center <- coords[, "Y"] # Add two new columns for the minimum x and y coordinates of each grid cell grid_intersected$X_min <- purrr::map_dbl(grid_intersected$geometry, ~ min(st_bbox(.)$xmin)) grid_intersected$Y_min <- purrr::map_dbl(grid_intersected$geometry, ~ min(st_bbox(.)$ymin)) grid_intersected$X_max <- purrr::map_dbl(grid_intersected$geometry, ~ min(st_bbox(.)$xmax)) grid_intersected$Y_max <- purrr::map_dbl(grid_intersected$geometry, ~ min(st_bbox(.)$ymax)) return(grid_intersected) } # Create centered grid and plot it quadrats_sf_fgeo_site_30 <- create_centered_grid(usable_area, grid_size) %>% mutate(quadrats_ids_30 = paste(X_center, Y_center, sep = "_")) %>% rename( X_center_30 = X_center, Y_center_30 = Y_center, X_min_30 = X_min, Y_min_30 = Y_min ) # Pre-calculate unique minimum X and Y values unique_x_min <- unique(quadrats_sf_fgeo_site_30$X_min) unique_y_min <- unique(quadrats_sf_fgeo_site_30$Y_min) # This function sums crown areas for trees and resets the running total when it # exceeds a predefined threshold (typically the plot area). assign_trees_to_area_groups <- function(x, threshold) { cumsum <- 0 group <- 0 result <- numeric(length(x)) group_vector <- integer(length(x)) for (i in 1:length(x)) { if (cumsum + x[i] > threshold) { group <- group + 1 cumsum <- x[i] } else { if (group == 0) { group <- 1 } cumsum <- cumsum + x[i] } result[i] <- cumsum group_vector[i] <- group } return(list(cumulative_sum = result, group = group_vector)) } # Intersect trees with quadrats intersection_result_fgeo_site_quadrat_30 <- st_intersection(fgeo_data_quad_allo_sf, quadrats_sf_fgeo_site_30) %>% mutate(intersect_area = st_area(.)) %>% group_by(quadrats_ids_30) %>% arrange(desc(caream2), .by_group = TRUE) %>% mutate( group_with_neighbors_quadrat_30 = assign_trees_to_area_groups(intersect_area, grid_size^2)$group ) %>% ungroup() %>% select(site, treeid, group_with_neighbors_quadrat_30, quadrats_ids_30) %>% st_drop_geometry() # Drop geometry once and keep it dropped cdf_30m <- base_fgeo_demo %>% select(site, name, treeid, gx, gy) %>% distinct() %>% mutate( gx_check = ifelse(gx %in% unique_x_min, "yes", "no"), gy_check = ifelse(gy %in% unique_y_min, "yes", "no"), ## Shift trees placed in the border gx_2 = ifelse(gx %in% unique_x_min, gx + 0.025, gx), gy_2 = ifelse(gy %in% unique_y_min, gy + 0.025, gy) ) # Convert to sf object and join cdf_30m_sf_object <- st_as_sf(cdf_30m, coords = c("gx_2", "gy_2"), crs = st_crs(fgeo_data_quad_allo_sf)) cdf_30m_joined <- st_join(cdf_30m_sf_object, quadrats_sf_fgeo_site_30, left = TRUE) # Combine with intersection results cdf_full_quadrat_30 <- cdf_30m_joined %>% left_join(intersection_result_fgeo_site_quadrat_30, by = c("site", "treeid", "quadrats_ids_30")) %>% # The lowest layers (>= 3) have similar height, suggesting similar light environment—we group them together. mutate( group_with_neighbors_quadrat_30_3_levels = ifelse(group_with_neighbors_quadrat_30 >= 3, 3, group_with_neighbors_quadrat_30), ) # ============================================================================= # EXTRACT CROWN POSITION DATA FOR JOINING # ============================================================================= # Extract crown layer assignments and quadrat information for each tree # Helper function to extract relevant columns from SF object process_crown_layers_data <- function(df, col_name) { df %>% ungroup() %>% sf::st_drop_geometry() %>% # Remove spatial geometry select(site, name, treeid, gx, gy, all_of(col_name)) # Keep tree ID and crown info } # Extract crown layer assignments for all trees crown_layers <- process_crown_layers_data(cdf_full_quadrat_30, c( "group_with_neighbors_quadrat_30", # Full crown position (1, 2, 3, 4, ...) "group_with_neighbors_quadrat_30_3_levels", # Simplified (1=canopy, 2=mid, 3=understory) "quadrats_ids_30", # Quadrat identifier "X_center_30", "Y_center_30", "X_min_30", "Y_min_30" # Quadrat coordinates )) # ============================================================================= # JOIN CROWN POSITION TO PRODUCTIVITY DATA (AWP) # ============================================================================= # Add crown layer assignments to tree-level productivity estimates # Helper function to assign recruited trees to bottom crown layer # Recruited trees don't have crown position data (not in census 1) # Assumption: Recruits are small and belong to understory (bottom layer) assign_recruits_to_bottom_layer <- function(df, group_var, fill_var) { # Separate rows with missing group_var (these are outside quadrats) na_rows <- df %>% filter(is.na({{ group_var }})) # For trees within quadrats, fill missing crown positions with maximum value # Maximum layer number = bottom/understory layer df %>% filter(!is.na({{ group_var }})) %>% group_by({{ group_var }}) %>% mutate({{ fill_var }} := ifelse(is.na({{ fill_var }}), max({{ fill_var }}, na.rm = TRUE), {{ fill_var }})) %>% ungroup() %>% bind_rows(na_rows) # Re-add trees outside quadrats } # Join crown position data to productivity estimates # Recruited trees get assigned to bottom layer (understory) AWP_QUAD_CA <- AWP %>% mutate(id = row_number()) %>% # Preserve original order left_join(., crown_layers, by = c("site", "name", "treeid", "gx", "gy")) %>% assign_recruits_to_bottom_layer(quadrats_ids_30, group_with_neighbors_quadrat_30_3_levels) %>% arrange(id) %>% # Restore original order select(-id) # ============================================================================= # ADD INITIAL BIOMASS TO PRODUCTIVITY DATA # ============================================================================= # Calculate initial (census 1) biomass for each tree (sum across all stems) # GET INITIAL BIOMASS FOR ALL TREES PRESENT IN CENSUS 1 init_agb_ba_vol_whole_tree <- base_fgeo_demo %>% filter(!is.na(dbhc_cm_census_1)) %>% # Only trees present in census 1 group_by(treeid) %>% reframe( initial_agb_whole_tree = sum(agb_census_1, na.rm = TRUE) # Total tree biomass (Mg) ) AWP_QUAD_CA <- AWP_QUAD_CA %>% left_join(., init_agb_ba_vol_whole_tree, by = "treeid") # ============================================================================= # JOIN CROWN POSITION TO MORTALITY DATA (AWM) # ============================================================================= # Add crown layer assignments to tree-level mortality estimates AWM_QUAD_CA <- AWM %>% mutate(id = row_number()) %>% # Preserve original order left_join(., crown_layers, by = c("site", "name", "treeid", "gx", "gy")) %>% arrange(id) %>% # Restore original order select(-id) # ============================================================================= # ADD INITIAL BIOMASS TO MORTALITY DATA # ============================================================================= AWM_QUAD_CA <- AWM_QUAD_CA %>% left_join(., init_agb_ba_vol_whole_tree, by = "treeid") # ============================================================================= # PREPARE DATASET FOR KOHYAMA BIOMASS TURNOVER CORRECTION # ============================================================================= # AWP and AWM were constructed at the stem level, and added up for each tree. full_tree_list <- bind_rows( AWP_QUAD_CA %>% select( site, treeid, group_with_neighbors_quadrat_30:Y_min_30 ), AWM_QUAD_CA %>% select(site, treeid, group_with_neighbors_quadrat_30:Y_min_30) ) %>% distinct() # Remove duplicates (to be safe) # Prepare productivity data for Kohyama correction # Calculate census 2 AGB using corrected growth rates to_link_1 <- AWP_QUAD_CA %>% select( site, treeid, agb_census_1, # Initial biomass (Mg) dAGB_corr, # Annual productivity (Mg/yr) dT # Time interval (years) ) %>% mutate( # Project forward: census 2 AGB = census 1 AGB + growth over interval agb_census_2 = agb_census_1 + (dT * dAGB_corr) ) %>% rename_with(~ paste0("A_", .), c(agb_census_1:agb_census_2)) # Prefix with "A_" for productivity # Prepare mortality data for Kohyama correction # For dead trees/stems, census 2 AGB is zero to_link_2 <- AWM_QUAD_CA %>% select( site, treeid, agb_census_1, # Initial biomass (Mg) dAGB_corr, # Annual mortality (Mg/yr) dT # Time interval (years) ) %>% mutate( agb_census_2 = 0 # Dead stems have zero biomass in census 2 ) %>% rename_with(~ paste0("B_", .), c(agb_census_1:agb_census_2)) # Prefix with "B_" for mortality # ============================================================================= # COMBINE PRODUCTIVITY AND MORTALITY AT TREE LEVEL # ============================================================================= # Join AWP and AWM for each tree # Some trees appear in both datasets (e.g., lost branches but trunk grew) # This is essential for the Kohyama correction join_all_trees <- full_tree_list %>% full_join(., to_link_1, by = c("site", "treeid")) %>% # Add productivity data full_join(., to_link_2, by = c("site", "treeid")) %>% # Add mortality data ungroup() %>% select(site:Y_min_30, sort(tidyselect::peek_vars())) # Sort columns for clarity # Calculate combined census biomasses # For trees with both growth and mortality, sum the components join_all_trees <- join_all_trees %>% mutate( # Census 1 AGB: Sum productivity and mortality components # Handles trees with both growing stems (A_) and dead stems (B_) agb_census_1 = coalesce(A_agb_census_1, 0) + coalesce(B_agb_census_1, 0), # Census 2 AGB: Sum surviving stems only (dead stems = 0) agb_census_2 = coalesce(A_agb_census_2, 0) + coalesce(B_agb_census_2, 0), # Time interval: Average if tree has both components dT = rowMeans(cbind(A_dT, B_dT), na.rm = TRUE) ) %>% select( # Remove intermediate columns, keep only combined values -A_agb_census_1, -A_agb_census_2, -B_agb_census_1, -B_agb_census_2, -A_dT, -B_dT, -A_dAGB_corr, -B_dAGB_corr ) # Join with initial whole-tree biomass df_kohyama_raw <- join_all_trees %>% left_join(., init_agb_ba_vol_whole_tree, by = "treeid") # ============================================================================= # KOHYAMA BIOMASS TURNOVER CORRECTION FUNCTION # ============================================================================= # Implements the Kohyama et al. (2019) correction for biomass turnover # Reference: https://doi.org/10.1016/j.foreco.2018.11.010 # # Parameters: # x1 : Initial biomass (census 1) # x2 : Final biomass (census 2) # xs1 : Biomass of surviving stems at census 1 # intvl : Time interval between censuses (years) # Area : Plot area for density calculation # subpop : Grouping variable (e.g., crown layer, quadrat) # simple : Return simple arithmetic estimates if TRUE # annual : Return annualized estimates if TRUE turnover <- function(x1, x2, xs1, intvl, Area, subpop = NULL, simple = FALSE, annual = FALSE) { # Default subpopulation: treat all trees as one group if (length(subpop) == 0) subpop <- rep(1, length(x1)) # Calculate biomass components by subpopulation B0k <- tapply(x1, subpop, sum) # Initial total biomass BTk <- tapply(x2, subpop, sum) # Final total biomass Bs0k <- tapply(xs1, subpop, sum) # Initial biomass of survivors # Logarithmic mean biomass (Eq 10 from Kohyama et al. 2019) # More appropriate than arithmetic mean for exponential growth Bwk <- ifelse(BTk != B0k, (BTk - B0k) / log(BTk / B0k), B0k) # Annualized logarithmic mean biomass (Eq 14 from Kohyama et al. 2019) Bw_annk <- ifelse(BTk != B0k, (BTk - B0k) / ((BTk / B0k)^(1 / intvl) - 1) / intvl, B0k) # Convert to per-area densities B0k <- B0k / Area BTk <- BTk / Area Bs0k <- Bs0k / Area Bwk <- Bwk / Area Bw_annk <- Bw_annk / Area # Simple estimates (arithmetic averages, no logarithmic correction) L_simple <- sum(B0k - Bs0k) / intvl # Simple mortality rate P_simple <- sum(BTk - Bs0k) / intvl # Simple productivity rate Bw_simple <- sum(B0k + BTk) / 2 # Arithmetic mean biomass # Annualized estimates (geometric progression) L_ann <- sum(Bw_annk * (1 - (Bs0k / B0k)^(1 / intvl))) # Annualized mortality P_ann <- sum(Bw_annk * (BTk / B0k)^(1 / intvl) * (1 - (Bs0k / BTk)^(1 / intvl))) # Annualized productivity Bw_ann <- sum(Bw_annk) # Annualized mean biomass # Kohyama estimates (logarithmic turnover rates - RECOMMENDED) L <- sum(Bwk * log(B0k / Bs0k)) / intvl # Kohyama mortality rate P <- sum(Bwk * log(BTk / Bs0k)) / intvl # Kohyama productivity rate Bw <- sum(Bwk) # Logarithmic mean biomass # Return requested estimate type if (simple == TRUE) { return(list("P_simple" = P_simple, "L_simple" = L_simple, "Bw_simple" = Bw_simple)) } else if (annual == TRUE) { return(list("P_annual" = P_ann, "L_annual" = L_ann, "Bw_annual" = Bw_ann)) } else { return(list("P" = P, "L" = L, "Bw" = Bw)) } } # ============================================================================= # DATA PROCESSING FUNCTION FOR KOHYAMA CORRECTION # ============================================================================= # Applies Kohyama correction to dataset grouped by crown position and quadrat # # Parameters: # dataset : Input data (df_kohyama_raw) # quadrat_size : Size of quadrat (e.g., 30 for 30x30m) # grouping_crown_column : Column name for crown position grouping # case_name : Label for this analysis case apply_kohyama_correction_by_crown_layer <- function( dataset = df_raw, quadrat_size = quadrat_size, grouping_crown_column = grouping_crown_column, case_name = case_name ) { print(paste0("Preparing ", case_name, ", sit tight....")) # Convert column names to symbols for tidy evaluation group_column <- sym(as.character(grouping_crown_column)) quadrat_column <- sym(paste0("quadrats_ids_", quadrat_size)) X_center <- sym(paste0("X_center_", quadrat_size)) X_min <- sym(paste0("X_min_", quadrat_size)) Y_center <- sym(paste0("Y_center_", quadrat_size)) Y_min <- sym(paste0("Y_min_", quadrat_size)) print("Define grouping variables...") # Aggregate to tree level (sum stem-level biomasses) df_processed <- dataset %>% group_by( treeid, {{ group_column }}, {{ quadrat_column }}, {{ X_center }}, {{ Y_center }}, {{ X_min }}, {{ Y_min }} ) %>% summarise( across(c(agb_census_1, agb_census_2), sum), # Sum biomass across stems dT = mean(dT), # Average time interval .groups = "drop" ) print("Prepare columns for correction...") # Nest data by crown position and quadrat processed_data <- df_processed %>% select( treeid, agb_census_1:dT, {{ group_column }}, {{ quadrat_column }}, {{ X_center }}, {{ Y_center }}, {{ X_min }}, {{ Y_min }} ) %>% filter(!is.na({{ quadrat_column }})) %>% # Remove trees outside quadrats group_by( {{ group_column }}, {{ quadrat_column }}, {{ X_center }}, {{ Y_center }}, {{ X_min }}, {{ Y_min }} ) %>% nest() %>% mutate( # Create unique identifier for each crown layer-quadrat combination link = paste0("crown_", as.character({{ group_column }}), "-", as.character({{ quadrat_column }})), # Prepare data for turnover function dss = map2(data, link, ~ .x %>% select(agb_census_1, agb_census_2, dT) %>% rename( w1_awp = agb_census_1, # Census 1 biomass w2_awp = agb_census_2 # Census 2 biomass ) %>% mutate(plot_id = .y) %>% select(plot_id, dT, w1_awp, w2_awp)) ) print("Biomass fluxes corrections...") # Apply Kohyama correction to each crown layer-quadrat combination processed_data <- processed_data %>% mutate(area = quadrat_size * quadrat_size) %>% # Quadrat area (m²) mutate(dT = map_dbl(dss, ~ mean(.x$dT))) %>% rownames_to_column("id") %>% mutate( # Calculate survivor biomass (xs1) - trees alive in both censuses dss_bio = map(dss, ~ .x %>% mutate(xs1 = ifelse(w1_awp > 0 & w2_awp > 0, w1_awp, 0))), # Apply Kohyama turnover correction (annualized version) cor_annual_bio = pmap( list(dss_bio, dT, area), function(a, b, c) { as.data.frame(turnover( x1 = a$w1_awp, # Census 1 biomass x2 = a$w2_awp, # Census 2 biomass xs1 = a$xs1, # Survivor biomass intvl = b, # Time interval Area = c, # Quadrat area subpop = NULL, # No sub-grouping simple = FALSE, # Use Kohyama method annual = TRUE # Return annualized rates )) } ), # Calculate initial biomass metrics initial_biomass = map_dbl(dss_bio, ~ sum(.x$w1_awp, na.rm = TRUE)), initial_biomass_Mg_ha = initial_biomass * (10000 / area) # Convert to Mg/ha ) print("Select corrected data...") # Extract and format results processed_data <- processed_data %>% select( {{ group_column }}, {{ quadrat_column }}, {{ X_center }}, {{ Y_center }}, {{ X_min }}, {{ Y_min }}, link, area, dT, cor_annual_bio, initial_biomass_Mg_ha ) %>% unnest( c( cor_annual_bio # Expand nested Kohyama results ), names_sep = "_" ) %>% # Convert productivity and mortality to Mg/ha/yr (from Mg/m²/yr) mutate(across(c( cor_annual_bio_P_annual, cor_annual_bio_L_annual ), ~ . * 10000, .names = "{.col}_Mg_ha_Yr")) %>% # Multiply by 10000 m²/ha mutate(case = {{ case_name }}) %>% # Add case label rename(crown = {{ group_column }}, quadrat_id = {{ quadrat_column }}) return(processed_data %>% ungroup()) } # ============================================================================= # APPLY KOHYAMA CORRECTION TO DATA # ============================================================================= # Process data for 30x30m quadrats with 3 crown position levels. This produces # Kohyama-corrected productivity estimates stratified by: # - Crown position (1=top, 2=middle, 3=bottom/understory) # - 30x30m quadrat df_kohyama_30_3_levels <- apply_kohyama_correction_by_crown_layer( dataset = df_kohyama_raw, quadrat_size = 30, grouping_crown_column = "group_with_neighbors_quadrat_30_3_levels", case_name = "case_30_by_30_3_levels" ) # ============================================================================= # END OF DATA PREPARATION SCRIPT # =============================================================================