Chapter 5 Consolidated Provincial Yield Table - NSYT and MSYT Combined Table

5.1 Overview

This workflow creates a single provincial yield table by combining Natural Stand Yield Tables (NSYT) generated from VDYP7 projections with Managed Stand Yield Tables (MSYT). The objective is to provide a standardized yield dataset for every forested feature in the provincial inventory using a common structure and age range. The final output contains conifer and deciduous merchantable volume (MVcon, MVdec) and dominant height (HTcon, HTdec) projections at 10-year intervals from age 0 to 250. The combined yield table provides a consistent provincial yield source for timber supply analysis, forest estate modelling, carbon accounting, and other downstream applications.

## Load Libraries
library(data.table)
library(dplyr)
library(tidyr)

5.2 Build MSYT component

5.2.1 Load MSYT tables

The provincial MSYT output and reference datasets are loaded to identify features assigned to managed stand yields and to retrieve the corresponding merchantable volume and dominant height projections.

msyt <- read.csv(
  "D://data//projects//MSYTDelivery//msyt_2025//prov//MSYT_prov_current_input_output.csv"
)

reference <- read.csv(
  "D://data//projects//MSYTDelivery//msyt_2025//prov//MSYT_prov_reference.csv"
)

5.2.2 Extend MSYT curves to age 250

Provincial MSYT yield curves are available to age 120. To align with the NSYT projection horizon, all volume and height attributes are extended to age 250 by carrying forward the age-120 values for each subsequent 10-year age interval. This produces a complete set of merchantable volume and dominant height fields from age 0 to 250 using the same structure as the NSYT outputs.

new_ages <- seq(130, 250, 10)

msyt_expanded <- msyt

for (age in new_ages) {

  msyt_expanded <- msyt_expanded %>%
    mutate(
      !!paste0("MVcon_", age) := MVcon_120,
      !!paste0("MVdec_", age) := MVdec_120,
      !!paste0("HTcon_", age) := HTcon_120,
      !!paste0("HTdec_", age) := HTdec_120
    )
}

### Reorder MSYT columns

ages <- seq(0, 250, 10)

ordered_cols <- c(
  "feature_id",
  paste0("MVcon_", ages),
  paste0("MVdec_", ages),
  paste0("HTcon_", ages),
  paste0("HTdec_", ages)
)

msyt_expanded <- msyt_expanded %>%
  select(all_of(ordered_cols))

5.2.3 Finalize MSYT output

The expanded MSYT yield curves are filtered to include only features that use managed stand yields. The resulting table is standardized to contain conifer and deciduous merchantable volume and dominant height projections for all age classes and is tagged with a source identifier of “MSYT”.

msyt_part <- reference %>%
  filter(!current_yield %in% c("VDYP", "Excluded")) %>%
  select("feature_id") %>%
  left_join(msyt_expanded, by = "feature_id") %>%
  mutate(source = "MSYT")

5.3 Build NSYT component

5.3.1 Load VDYP yield table

VDYP7 yield projections generated through the provincial NSYT workflow are loaded and converted to a data.table structure for processing. These projections contain age-specific merchantable volume, dominant height, and species-level yield information for natural stands.

yld <- readRDS(
  "D:\\data\\projects\\VDYP\\VDYP\\2024\\prov\\VDYP7RunParallel\\prov_lesschanges_nsyt_yieldtable.rds"
)

setDT(yld)

yld_trim<-yld%>%select(FEATURE_ID, PRJ_TOTAL_AGE, PRJ_VOL_WS, PRJ_VOL_DWB)
saveRDS(yld_trim, "yld_trim.rds")

5.3.2 Separate conifer volume and deciduous volume

5.3.2.1 Define deciduous species

Species in this list are treated as deciduous

decid <- c(
  "AC","ACB","ACT","AT","AD","AX",
  "D","DG","DM","DR","E","EA","EB","EE","EP","ES","EW","EXP",
  "G","GP","KC","M","MB","MN","MV","Q","QG",
  "WT","WS","WP","WD","WB","WA","W", 
  "VP","VB","V","U","RA","ZH","XH"
)

NoDecidVol_calculator() calculates conifer-only projected volume by removing deciduous species from VDYP projected species volumes.

NoDecidVol_calculator <- function (sp_list,
                                   vol_list,
                                   decidsp,
                                   multilayer = FALSE) {

  # Standardize species slot names
  sp_list[spcode %in% "SPECIES_1_CODE", `:=`(spcode, "sp1")]
  sp_list[spcode %in% "SPECIES_2_CODE", `:=`(spcode, "sp2")]
  sp_list[spcode %in% "SPECIES_3_CODE", `:=`(spcode, "sp3")]
  sp_list[spcode %in% "SPECIES_4_CODE", `:=`(spcode, "sp4")]
  sp_list[spcode %in% "SPECIES_5_CODE", `:=`(spcode, "sp5")]
  sp_list[spcode %in% "SPECIES_6_CODE", `:=`(spcode, "sp6")]

  # Remove deciduous species
  sp_list2 <- sp_list[
    !is.na(sp) &
      sp != "" &
      !sp %in% decidsp
  ]

  # Replace NA projected volumes with 0
  vol_list[is.na(PRJ_SP1_VOL_DWB), `:=`(PRJ_SP1_VOL_DWB, 0)]
  vol_list[is.na(PRJ_SP2_VOL_DWB), `:=`(PRJ_SP2_VOL_DWB, 0)]
  vol_list[is.na(PRJ_SP3_VOL_DWB), `:=`(PRJ_SP3_VOL_DWB, 0)]
  vol_list[is.na(PRJ_SP4_VOL_DWB), `:=`(PRJ_SP4_VOL_DWB, 0)]
  vol_list[is.na(PRJ_SP5_VOL_DWB), `:=`(PRJ_SP5_VOL_DWB, 0)]
  vol_list[is.na(PRJ_SP6_VOL_DWB), `:=`(PRJ_SP6_VOL_DWB, 0)]

  # Convert projected species volumes to long format
    vol_list2 <- melt(
      vol_list,
      id.vars = c("FEATURE_ID", "PRJ_TOTAL_AGE"),
      measure.vars = c(
        "PRJ_SP1_VOL_DWB",
        "PRJ_SP2_VOL_DWB",
        "PRJ_SP3_VOL_DWB",
        "PRJ_SP4_VOL_DWB",
        "PRJ_SP5_VOL_DWB",
        "PRJ_SP6_VOL_DWB"
      ),
      variable.name = "spcode",
      value.name = "vol"
    )

  # Standardize projected volume species names
  vol_list2[spcode %in% "PRJ_SP1_VOL_DWB", `:=`(spcode, "sp1")]
  vol_list2[spcode %in% "PRJ_SP2_VOL_DWB", `:=`(spcode, "sp2")]
  vol_list2[spcode %in% "PRJ_SP3_VOL_DWB", `:=`(spcode, "sp3")]
  vol_list2[spcode %in% "PRJ_SP4_VOL_DWB", `:=`(spcode, "sp4")]
  vol_list2[spcode %in% "PRJ_SP5_VOL_DWB", `:=`(spcode, "sp5")]
  vol_list2[spcode %in% "PRJ_SP6_VOL_DWB", `:=`(spcode, "sp6")]

  # Merge species and volume tables
    all <- merge(
      sp_list2,
      vol_list2,
      by = c("FEATURE_ID", "spcode"),
      all.x = TRUE
    )

    all2 <- all[
      ,
      .(PRJ_VOL_DWB_conifer = sum(vol, na.rm = TRUE)),
      by = c("FEATURE_ID", "PRJ_TOTAL_AGE")
    ]

  return(all2)
}

5.3.2.2 Extract species codes

Convert species fields into long format.

d <- unique(
  yld[, .(
    FEATURE_ID,
    SPECIES_1_CODE,
    SPECIES_2_CODE,
    SPECIES_3_CODE,
    SPECIES_4_CODE,
    SPECIES_5_CODE,
    SPECIES_6_CODE
  )]
)

d2 <- melt(
  d,
  id.vars = "FEATURE_ID",
  measure.vars = c(
    "SPECIES_1_CODE",
    "SPECIES_2_CODE",
    "SPECIES_3_CODE",
    "SPECIES_4_CODE",
    "SPECIES_5_CODE",
    "SPECIES_6_CODE"
  ),
  variable.name = "spcode",
  value.name = "sp"
)

5.3.2.3 Extract projected volumes

e <- unique(
  yld[, .(
    FEATURE_ID,
    PRJ_TOTAL_AGE,
    PRJ_SP1_VOL_DWB,
    PRJ_SP2_VOL_DWB,
    PRJ_SP3_VOL_DWB,
    PRJ_SP4_VOL_DWB,
    PRJ_SP5_VOL_DWB,
    PRJ_SP6_VOL_DWB
  )]
)

t <- unique(
  yld[, .(
    FEATURE_ID,
    PRJ_TOTAL_AGE,
    PRJ_VOL_DWB
  )]
)

5.3.2.4 Calculate conifer and deciduous Volume

f <- NoDecidVol_calculator(
  sp_list = d2,
  vol_list = e,
  decidsp = decid
)

final <- t %>%
  left_join(f, by = join_by(FEATURE_ID, PRJ_TOTAL_AGE)) %>%
  mutate(
    PRJ_VOL_DWB = round(PRJ_VOL_DWB, 2),
    PRJ_VOL_DWB_conifer = round(PRJ_VOL_DWB_conifer, 2),

    PRJ_VOL_DWB_decid =
      PRJ_VOL_DWB - coalesce(PRJ_VOL_DWB_conifer, 0)
  )

#convert to wide format
final_wide <- final %>%
  mutate(
    age_label = paste0("_", PRJ_TOTAL_AGE)
  ) %>%
  pivot_wider(
    id_cols = FEATURE_ID,
    names_from = age_label,
    values_from = c(
      PRJ_VOL_DWB_conifer,
      PRJ_VOL_DWB_decid
    ),
    names_glue =
      "{ifelse(.value == 'PRJ_VOL_DWB_conifer', 'MVcon', 'MVdec')}{age_label}"
  )

5.3.3 Separate conifer height and deciduous height

5.3.3.1 Develop Height Ratio Tables

The following develops a set of BEC-specific conifer-to-deciduous dominant height ratios that are subsequently used to separate VDYP projected dominant heights into conifer and deciduous components for the provincial yield table. From the MSYT outputs, only features containing valid conifer (HTcon_120) and deciduous (HTdec_120) dominant heights at age 120 are retained. These features are then joined to their corresponding Biogeoclimatic Ecosystem Classification (BEC) zone and subzone information, creating a unique BEC label for each feature. For every stand, the script calculates both the deciduous-to-conifer height ratio (HTdec_120 / HTcon_120) and the conifer-to-deciduous height ratio (HTcon_120 / HTdec_120). These ratios are then summarized by BEC zone and subzone, producing average height relationships for each ecological unit based on all available managed stand observations. The resulting table contains the number of contributing samples and mean height ratios for each BEC variant and is exported as bec_height_ratios.csv. These ratios are subsequently used within the provincial NSYT workflow to estimate separate conifer and deciduous dominant height projections from a single VDYP dominant height prediction, ensuring consistency between the provincial managed stand and natural stand yield datasets.

library(RPostgreSQL)
library(glue)

schema <- 'msyt_2025'
opt <- paste0("-c search_path=",schema)
user_name <- 'results'
database <- 'msyt'
db<- dbConnect('PostgreSQL',dbname=database,user=user_name,options = opt)

hts<-msyt

hts<-hts%>%select(feature_id,HTcon_120,HTdec_120)%>%
  filter(HTdec_120 > 0 & HTcon_120 > 0)
glimpse(hts)


sql<-glue_sql("select * from bec12",.con = db)

bec<-dbGetQuery(db,sql)
bec<-bec%>%mutate(
  bec_label = paste0(zone,"_",subzone))%>%
  select(feature_id,bec_label)

hts<-inner_join(hts,bec,join_by("feature_id"))

ht_ratio<-hts%>%
  mutate(
    dec_con_ratio = HTdec_120/HTcon_120,,
    con_dec_ratio = HTcon_120/HTdec_120)%>%
    group_by(bec_label)%>%
  summarise(
    N = n(),
    dec_con_ratio = round(mean(dec_con_ratio),2),
    con_dec_ratio = round(mean(con_dec_ratio),2)
  )

write_csv(ht_ratio,"./bec_height_ratios.csv")

5.3.3.2 Calculate conifer and deciduous height

VDYP dominant height projections are joined to the BEC-specific height ratio table and evaluated using the leading species of each stand. The resulting conifer and deciduous height projections are reshaped into a wide table containing HTcon_* and HTdec_* fields for all projected ages.

vdyp_ht <- unique(
  yld[, .(
    FEATURE_ID,
    PRJ_TOTAL_AGE,
    PRJ_DOM_HT,
    SPECIES_1_CODE
  )]
)

vdyp_ht <- vdyp_ht %>%
  left_join(
    reference %>%
      select(feature_id, bec_zone, bec_subzone) %>%
      mutate(
        bec_label = paste0(bec_zone, "_", bec_subzone)
      ),
    by = c("FEATURE_ID" = "feature_id")
  )

vdyp_ht <- vdyp_ht %>%
  left_join(ht_ratio, by = "bec_label")

vdyp_ht <- vdyp_ht %>%
  mutate(
    is_decid = SPECIES_1_CODE %in% decid,

    HTcon = case_when(
      is_decid ~ PRJ_DOM_HT * con_dec_ratio,
      TRUE ~ PRJ_DOM_HT
    ),

    HTdec = case_when(
      is_decid ~ PRJ_DOM_HT,
      TRUE ~ PRJ_DOM_HT * dec_con_ratio
    )
  )

vdyp_ht_wide <- vdyp_ht %>%
  mutate(
    age_label = paste0("_", PRJ_TOTAL_AGE)
  ) %>%
  pivot_wider(
    id_cols = FEATURE_ID,
    names_from = age_label,
    values_from = c(HTcon, HTdec),
    names_glue = "{.value}{age_label}"
  )

5.3.4 Finalize NSYT tables

The conifer and deciduous merchantable volume projections are combined with the conifer and deciduous height projections to create the final NSYT yield table. The output is restricted to features assigned to VDYP-based yields and is tagged with a source identifier of “NSYT”.

vdyp_part <- reference %>%
  filter(current_yield %in% c("VDYP", "Excluded")) %>%
  select("feature_id") %>%
  left_join(
    final_wide,
    by = c("feature_id" = "FEATURE_ID")
  ) %>%
  left_join(
    vdyp_ht_wide,
    by = c("feature_id" = "FEATURE_ID")
  ) %>%
  mutate(source = "NSYT")

5.4 Combine NSYT tables and MSYT tables and export

The finalized NSYT and MSYT yield tables are appended into a single provincial dataset containing standardized merchantable volume and dominant height projections for all forested features. The combined table is exported as yld_combined.csv and serves as the provincial yield library for downstream modelling applications.

final_yld <- bind_rows(
  msyt_part,
  vdyp_part
)

write.csv(
  final_yld,
  "D:\\data\\projects\\VDYP\\VDYP\\2024\\prov\\yld_combined.csv",
  row.names = FALSE
)

5.5 QA/QC checks

A final quality assurance review identifies records that contain zero merchantable volume across all projected age classes. These features are flagged for review because they may indicate missing yield information, invalid stand attributes, or yield assignment issues that require further investigation. Features with zero projected volume were identified, and the likely causes of missing yield estimates were summarized using preprocessing logs, validation checks, VDYP warnings, and VDYP error messages.

user_name <- 'postgres'
data_base <- 'vdyp'
 
con <- dbConnect('PostgreSQL',dbname=data_base,user=user_name)

err <- readRDS(
  "D:\\data\\projects\\VDYP\\VDYP\\2024\\prov\\VDYP7RunParallel\\prov_lesschanges_nsyt_errorTable.rds"
)
setDT(err)

warning <- readRDS(
  "D:\\data\\projects\\VDYP\\VDYP\\2024\\prov\\VDYP7RunParallel\\prov_lesschanges_nsyt_warningTable.rds"
)
setDT(warning)

query <- "select * from  prov_ba_age_ht_checks_2024"
prov_ba_age_ht_checks_2024<- dbGetQuery(con, query)

query <- "select * from  prov_updated_other_veg_2024"
prov_updated_other_veg_2024<- dbGetQuery(con, query)

query <- "select * from  prov_updated_ba_age_ht_2024"
prov_updated_ba_age_ht_2024<- dbGetQuery(con, query)

query <- "select * from  prov_updated_spp_2024"
prov_updated_spp_2024<- dbGetQuery(con, query)

other_veg<- prov_updated_other_veg_2024%>% select(feature_id, change_log)%>%rename(other_veg_change_log = change_log)
ah<- prov_updated_ba_age_ht_2024%>% select(feature_id, change_log)%>%rename(layer_attri_change_log = change_log)
spp <- prov_updated_spp_2024%>% select(feature_id, change_log)%>%rename(spp_change_log = change_log)
ah_check<- prov_ba_age_ht_checks_2024%>%rename(layer_attri_checks = change_log)
yld_si<- yld%>% select(FEATURE_ID, PRJ_SITE_INDEX)%>%rename(feature_id = FEATURE_ID)%>%distinct()
yld_err <- err %>%select(FEATURE_ID, MESSAGE) %>%
  rename(
    feature_id = FEATURE_ID,
    err_message = MESSAGE
  )%>%
  group_by(feature_id) %>%
  summarise(
    err_message = paste(unique(err_message), collapse = "; "),
    .groups = "drop"
  )%>%distinct()

yld_war<- warning%>% select(FEATURE_ID, MESSAGE) %>% 
  rename(
    feature_id = FEATURE_ID,
    war_message = MESSAGE
  )%>%
  group_by(feature_id) %>%
  summarise(
    war_message = paste(unique(war_message), collapse = "; "),
    .groups = "drop"
  )%>%distinct()

mv_cols <- grep(
  "^MV(con|dec)_",
  names(final_yld),
  value = TRUE
)

zero_features <- final_yld %>%
  mutate(
    total_mv = rowSums(
      across(all_of(mv_cols)),
      na.rm = TRUE
    )
  ) %>%
  filter(total_mv == 0) %>%
  select(feature_id, source)

zero_features_errorlog <- zero_features%>%
  left_join(other_veg, by = "feature_id")%>%
  left_join(ah, by = "feature_id")%>%
  left_join(spp, by = "feature_id")%>%
  left_join(ah_check, by = "feature_id")%>%
  left_join(yld_si, by = "feature_id")%>%
  left_join(yld_err, by = "feature_id")%>%
  left_join(yld_war, by = "feature_id")  

write_csv(zero_features_errorlog,"./zero_features_errorlog.csv")

zero_features_errorlog <- zero_features_errorlog %>%
  mutate(
    type = case_when(
      source =='MSYT' ~ "MSYT missing yield",
      grepl("NO_LAYER", other_veg_change_log, ignore.case = TRUE) ~ "No Treed Layer", 
      grepl("MISSING_BEC", other_veg_change_log, ignore.case = TRUE) ~ "Missing BEC",  
      grepl("NO_VALID_SPECIES", spp_change_log, ignore.case = TRUE) ~  "No Valid Species",
      grepl("missing or low est_age_spp1", layer_attri_checks, ignore.case = TRUE) ~  "Missing or Low Species1 Age Or Height",       
      grepl("missing or low est_height_spp1", layer_attri_checks, ignore.case = TRUE) ~   "Missing or Low Species1 Age Or Height",
      grepl("missing or low est_age_spp2 ", layer_attri_checks, ignore.case = TRUE) ~  "Missing or Low Species2 Age Or Height",
      grepl("missing or low est_height_spp2", layer_attri_checks, ignore.case = TRUE) ~  "Missing or Low Species2 Age Or Height",
      grepl("Very low or null site index", err_message, ignore.case = TRUE) ~  "Low Site Index",
      grepl("Very low or null site index", war_message, ignore.case = TRUE) ~  "Low Site Index",      
      grepl("VRI Inventory Standard but BA missing", war_message, ignore.case = TRUE) ~  "Missing or Low Basal Area",    
      grepl("Predicted Basal Area too small to be used", war_message, ignore.case = TRUE) ~  "Missing or Low Basal Area",
      grepl("BA null or lower than 1", layer_attri_checks, ignore.case = TRUE) ~  "Missing or Low Basal Area",  
      grepl("VRI Inventory Standard but TPH missing", war_message, ignore.case = TRUE) ~  "Missing or Low Density",   
      grepl("stems_per_ha_75 null or lower than 1", layer_attri_checks, ignore.case = TRUE) ~  "Missing or Low Density",     
      grepl("Unable to Project Stand backward beyond 400 years from stand reference age", war_message, ignore.case = TRUE) ~  "Reference Age > 400",  
      PRJ_SITE_INDEX <=5 ~ "Low Site Index",  
      grepl("Primary Layer height was too short to generate a stand description", war_message, ignore.case = TRUE) ~  "Ht L1 Too Low",  
      grepl("Breast Height Age too young to generate a stand description", war_message, ignore.case = TRUE) ~  "Age Too Low",  
      is.na(other_veg_change_log) ~ "Missing From VDYP input",
      TRUE ~ "Unknow"
    )
  )

err_summary_prov<-zero_features_errorlog %>%
  group_by(type) %>%
  summarise(
    polygon_count = n_distinct(feature_id),
    .groups = "drop"
  ) %>%
  arrange(desc(polygon_count))

saveRDS(err_summary_prov, "data/tables/err_summary_prov.rds")

Table 1. Summary of Possible Causes for Features Missing Yields

type	polygon_count
No Treed Layer	494183
Low Site Index	70123
Missing From VDYP input	47100
Missing or Low Basal Area	37911
No Valid Species	23393
Reference Age > 400	403
MSYT missing yield	88
Unknow	64
Missing BEC	18
Missing or Low Density	13