1  Workflow Overview

The {targets} R package provides a make-like data analysis pipeline toolkit. Here, we use it to track the progress and reproducibility of this R workflow. By using this approach, all aspects of the data analysis from raw inputs to delivered calculations are tracked; any changes that are made will be detected necessitating the recalculation of all dependent analytical steps.

The following interactive diagram displays all of the “targets” for this project. Targets are discrete sections of the analysis (functions) which fulfil a specific purpose. This demonstrates the flow of information between the various aspects of the pipeline.

1.1 Code

To view the code sections, simply click on a subheading which is preceded by an arrow.

The _targets.R script defines the data analysis pipeline. Each target, defined by the function tar_target can be considered a discrete step in the workflow. The functions called in this script are also referenced throughout this book.

Main script to control the pipeline

library(targets)
library(tarchetypes)

tar_option_set(trust_timestamps = TRUE)

# ------ Set output options -----------------------
results_version <- "4-3" # > 4-0 for mr2
make_gh_release <- FALSE
release_type <- "minor" # must be "major" or "minor"
gh_release_body <- "Kuamut pipeline release: enable save supression in ERRs to faciliate Auditor testing - for PRR."
# ----------- Load functions ----------------------
tar_source()

# ----------- Define Inputs ------------------------
DATA_VERSION <- "v0_3_4"
# Define input data
source("targets_components/input_data_paths.R")

# ------------------- Define ML inputs --------------------------
ML_VERSION_mr1 <- "v1-2"
ML_VERSION_mr2 <- "v2-0"

source("targets_components/ml_data_paths.R")
# -------------------- Set package deps ----------------------------

source("targets_components/targets_set_options.R")

# ---- Data download targets (object named `TARGETS_DOWNLOAD`) ---------------
source("targets_components/targets_download_pipe.R")

# ---- main pipeline (object named `TARGETS_MAIN`) ---------------------------
source("targets_components/targets_main_pipe.R")

# ---- leaflet map targets (object named `LEAF_MAPS`) ------------------------
source("targets_components/targets_leaflet_pipe.R")

# ---- ggplot targets (object named `GGPLOT_TARGETS`) ------------------------
source("targets_components/targets_ggplots_pipe.R")

# ---- saving outputs including `CREATE_GH_RELEASE` and `EXPORT_OUTPUTS` -----
source("targets_components/targets_outputs_pipe.R")

# ---- build the pipeline book/webpage including `BUILD_REFS`, `BUILD_BOOK`
source("targets_components/targets_build_book_pipe.R")

# and `BUILD_BOOK_MINIMAL` ---------------------------------------------------

list(
  TARGETS_DOWNLOAD,
  TARGETS_MAIN,
  LEAF_MAPS,
  GGPLOT_TARGETS,
  CREATE_GH_RELEASE,
  BUILD_REFS,
  BUILD_BOOK,
  BUILD_BOOK_MINIMAL
)

Sourced components of the pipeline

# Here, the input paths for the pipeline are defined.
data_root <- paste0("data_in_NOEDIT_", DATA_VERSION)

kuamut_plot_data16 <- ("PlotData/FullPlotDataKuamut_FINAL_CDP_rechecked.csv")
kuamut_plot_data22 <- ("PlotData/FullPlotDataKuamut_24.6.2022.csv")
kuamut_plot_locs <- "PlotData/KuamutLocations.fgb"

kuamut_bounds <- "boundaries/Kuamut_boundaries_dissolved.fgb"
kuamut_logging_coups <- "coups2/coups_modified2.fgb"

SabahC_acd <- "CAO_Sabah_ACD/CAO_Sabah_ACD_30m_Masked.tif"
kuamut_dtm <- "Finished_Surface_Maps/Kuamut_Cleaned_ground.tif"
kuamut_dsm <- "Finished_Surface_Maps/Kuamut_Cleaned_tch.tif"
gao_chm <- "GAO_TCH_Map/GAO_TCH_30m_unmasked.tif"

timb_prod <- "SFD_TimberProduction/SabahTimberProduction2010to2019.csv"
timb_prod_MR <- "SFD_TimberProduction/SabahTimberProduction2020to2021.csv"

roads_detail_src <- "riversRoads_all/kuamut_all_roads_detail.fgb"
roads_src <- "riversRoads_all/roads_all.fgb"
rivers_src <- "riversRoads_all/rivers_all.fgb"
protectedArea_src <- "forest_reserve/forest_reserve.fgb"
KuamutRiver_src <- "river_polygon/river_poly.fgb"

forest_reserves_src <- "ForestReserves/ForestReserves.fgb"
forest_res_names <- "ForestReserves/ForestReserveNames.csv"

reference_areas_src <- "reference_areas/Kuamut_Reference_areas_Final.fgb"
pinangah_5k_pdf_src <- paste0(
  "reference_areas/",
  "PD_ReferenceArea_FMU 16 Sg Pinangah 5000ha_Compressed.pdf"
)
pinangah_20k_pdf_src <- paste0(
  "reference_areas/",
  "PD_ReferenceArea_FMU 16 Sg Pinangah 20000ha_Compressed.pdf"
)
gunung_rara_pdf_src <- paste0(
  "reference_areas/",
  "PD_ReferenceArea_FMU 23 Gunung Rara 28000ha_Compressed.pdf"
)

protected_taxa <- "protected_tree_sp/sabah-forestry-prot-taxa.csv"

man_defor23 <- "manual_digitise_deforest/man-deforest-2023/man-defor-23.fgb"

# Here we define the file paths for the machine learning data.
ML_ROOT_mr1 <- paste0("ML_data_in_NOEDIT", ML_VERSION_mr1)
ML_ROOT_mr2 <- paste0("ML_data_in_NOEDIT", ML_VERSION_mr2)
# Define input data
MerchVol_10cm <- "MerchVol_maps/PLOT_VAR_CLEAN_PREDICT_MerchVol10_ha2016_30.tif"
MerchVol_31cm <- "MerchVol_maps/PLOT_VAR_CLEAN_PREDICT_MerchVol31_ha2016_30.tif"
MerchVol_40cm <- "MerchVol_maps/PLOT_VAR_CLEAN_PREDICT_MerchVol40_ha2016_30.tif"

ACD_map_2016 <- "ACD_maps/PLOT_VAR_CLEAN_PREDICT_ACD_mgC_ha2016_30.tif"
# below files not explicitly used in this pipeline.
ACD_map_2017 <- "ACD_maps/PLOT_VAR_CLEAN_PREDICT_ACD_mgC_ha2017_30.tif"
ACD_map_2018 <- "ACD_maps/PLOT_VAR_CLEAN_PREDICT_ACD_mgC_ha2018_30.tif"
ACD_map_2019 <- "ACD_maps/PLOT_VAR_CLEAN_PREDICT_ACD_mgC_ha2019_30.tif"
ACD_map_2020 <- "ACD_maps/PLOT_VAR_CLEAN_PREDICT_ACD_mgC_ha2020_30.tif"
ACD_map_2021 <- "ACD_maps/PLOT_VAR_CLEAN_PREDICT_ACD_mgC_ha2021_30.tif"
ACD_map_2022 <- "ACD_maps/PLOT_VAR_CLEAN_PREDICT_ACD_mgC_ha2022_30.tif"

# for mr2 we use a new ACD map.
ACD_map_2023 <- "ACD_maps/PLOT_VAR_CLEAN_PREDICT_ACD_mgC_ha2023_30.tif"

tar_option_set(
  packages = c(
    # general
    "here", "tidyr", "measurements", "purrr", "stringr", "units", "stringr",
    "units", "readr", "data.table", "furrr", "future",
    # ecological
    "BIOMASS",
    # spatial
    "rmapshaper", "sf", "exactextractr", "terra", "tmap", "ezwarp", "demsrc",
    # statistics
    "gt", "Rmisc",
    # visualisation
    "viridis", "ggplot2", "scico", "ggpmisc", "ggtext", "ggdist", "patchwork",
    # dplyr to avoid plyr conflcts
    "dplyr"
  )
)

TARGETS_DOWNLOAD <- list(

  # target required for quarto report.
  tar_target(
    Data_versions,
    list(
      in_data_v = DATA_VERSION,
      in_ml_data_v = ML_VERSION_mr1,
      in_ml_data_v_mr2 = ML_VERSION_mr2
    )
  ),
  tar_skip(Download_Input_Data,
    download_gh_release_asset(
      src_target = "kuamut_data_in",
      dir_target = data_root,
      data_version = DATA_VERSION,
      repo = "kuamut-in-data"
    ),
    skip = dir.exists(data_root)
  ),
  tar_skip(Download_Input_ML_Data,
    download_gh_release_asset(
      src_target = "Kuamut_ACD_MV_outputs",
      dir_target = ML_ROOT_mr1,
      data_version = ML_VERSION_mr1,
      repo = "Kuamut-Biomass-ML"
    ),
    skip = dir.exists(ML_ROOT_mr1)
  ),
  tar_skip(Download_Input_ML_Data_mr2,
    download_gh_release_asset(
      src_target = "Kuamut_ACD_MV_outputs",
      dir_target = ML_ROOT_mr2,
      data_version = ML_VERSION_mr2,
      repo = "Kuamut-Biomass-ML"
    ),
    skip = dir.exists(ML_ROOT_mr2)
  )
)

TARGETS_MAIN <- list(
  tar_target(
    CF,
    check_existence(
      kuamut_plot_data16,
      kuamut_plot_data22,
      kuamut_bounds,
      kuamut_plot_locs,
      kuamut_logging_coups,
      SabahC_acd,
      kuamut_dtm,
      kuamut_dsm,
      gao_chm,
      roads_detail_src,
      roads_src,
      rivers_src,
      protectedArea_src,
      KuamutRiver_src,
      timb_prod,
      timb_prod_MR,
      forest_reserves_src,
      forest_res_names,
      reference_areas_src,
      pinangah_5k_pdf_src,
      pinangah_20k_pdf_src,
      gunung_rara_pdf_src,
      protected_taxa,
      man_defor23,
      dr = data_root,
      skip_target = Download_Input_Data
    )
  ),
  tar_target(
    MLF,
    check_existence(
      MerchVol_10cm,
      MerchVol_31cm,
      MerchVol_40cm,
      ACD_map_2016,
      ACD_map_2017,
      ACD_map_2018,
      ACD_map_2019,
      ACD_map_2020,
      ACD_map_2021,
      ACD_map_2022,
      dr = ML_ROOT_mr1,
      skip_target = Download_Input_ML_Data
    )
  ),
  tar_target(
    MLF_mr2,
    check_existence(
      ACD_map_2023,
      dr = ML_ROOT_mr2,
      skip_target = Download_Input_ML_Data_mr2
    )
  ),
  tar_target(
    acd_path_list,
    c(MLF[grep("ACD", names(MLF))], MLF_mr2)
  ),
  tar_target(
    TreeAllometrics,
    Kuamut_TreeAllometrics(
      CF$kuamut_plot_data16,
      CF$kuamut_plot_data22
    )
  ),

  # from "R/geoid_tools.R"
  tar_target(
    DTM_geoid_correct,
    geoid_correct_dtm_gdal(CF$kuamut_dtm),
    format = "file"
  ),
  # from "R/mergeDEM_Sabah.R"
  tar_target(
    MergeDEM,
    merge_dem_srcs(CF$kuamut_bounds, DTM_geoid_correct),
    format = "file"
  ),

  # from "R/SlopeCorrection.R"
  tar_target(
    SlopeCorrection,
    slope_correction(
      dem_path = MergeDEM,
      points_path = CF$kuamut_plot_locs
    )
  ),
  # from "R/KuamutTree2Plot.R"
  tar_target(
    TreeToPlotMetrics,
    tree_to_plot(
      TreeAllometrics,
      SlopeCorrection$slp_correct_spat
    )
  ),
  # from 'R/slopeAnalysis.R'
  tar_target(SlopeAnalysis, slope_analysis(MergeDEM)),

  # from 'R/DefineStrata.R'
  tar_target(
    DefineStrata,
    define_strata(
      slope_src = SlopeAnalysis$slope_raster,
      boundary_src = CF$kuamut_bounds,
      flat_zones_src = SlopeAnalysis$flat_areas,
      logg_coups_src = CF$kuamut_logging_coups,
      CF$roads_src,
      CF$rivers_src,
      CF$protectedArea_src,
      CF$KuamutRiver_src
    )
  ),
  tar_target(
    Deforestation_Extract,
    command = deforest_extract(CF$roads_detail_src, outfile_version = "v0_5"),
    format = "file"
  ),
  tar_target(
    bind_auto_manual_deforest,
    command = rbind_auto_manual_deforest(
      Deforestation_Extract,
      CF$man_defor23,
      outfile_version = "v0_5"
    ),
    format = "file"
  ),
  tar_target(
    name = deforest_emission_calcs,
    command = {
      deforest_emissions(acd_path_list, bind_auto_manual_deforest)
    }
  ),
  tar_target(
    sp_conserv_factor,
    calc_sp_conserv_factor(
      tree_mets = TreeAllometrics,
      plot_locs = CF$kuamut_plot_locs,
      strata = DefineStrata$mgmt_strata_union,
      prot_taxa = CF$protected_taxa
    )
  ),

  # from 'R/TimberHarvestPlan_Mosaic_Coupes.R'
  tar_target(
    TimberHarvestPlanMosaic,
    timber_harvest_plan_mosaic(
      boundary_src = CF$kuamut_bounds,
      MosaicCoupesGo_src = DefineStrata$MosaicCoupesGo,
      MosaicPLANT_src = DefineStrata$MosaicPLANT,
      MosaicUNplanted_src = DefineStrata$MosaicUNplanted,
      dbh31_src = MLF$MerchVol_31cm,
      dbh10_src = MLF$MerchVol_10cm,
      acd_src = MLF$ACD_map_2016,
      sp_conserve_factor = sp_conserv_factor$spc_mosaic
    )
  ),

  # from "R/TimberHarvestPlan_NFM_Coupes.R"
  tar_target(
    TimberHarvestPlanNFM,
    timber_harvest_plan_NFM(
      boundary_src = CF$kuamut_bounds,
      NFMCoupes_src = DefineStrata$NFMCoupesGo,
      dbh40_src = MLF$MerchVol_40cm,
      dbh10_src = MLF$MerchVol_10cm,
      acd_src = MLF$ACD_map_2016,
      sp_conserve_factor = sp_conserv_factor$spc_nfm
    )
  ),

  #  from "R/ParcelCalculations_VM10v1.3.R"
  tar_target(
    ParcelCalculations,
    parcel_calculations(
      V_EX_31cmDBH_csv = TimberHarvestPlanMosaic$`31cmDBH`,
      V_EX_10cmDBH_csv = TimberHarvestPlanMosaic$`10cmDBH`,
      V_EX_ACD_10cm_csv = TimberHarvestPlanMosaic$ACD_10cm,
      NFM_40cmDBH_csv = TimberHarvestPlanNFM$`40cmDBH`,
      NFM_10cmDBH_csv = TimberHarvestPlanNFM$`10cmDBH`,
      NFM_ACD_10cm = TimberHarvestPlanNFM$ACD_10cm,
      NoGo_src = DefineStrata$NoGo,
      WD = TreeToPlotMetrics$mean_plotWD
    )
  ),

  # from 'R/MosaicAreaGrowthRate.R'
  tar_target(
    MosaicAreaGrowthRate,
    mosaic_area_growth_rate(
      MosaicPLANT_src = DefineStrata$MosaicPLANT,
      Mosaic_Parcels_csv = TimberHarvestPlanMosaic$`31cmDBH`,
      Mosaic10cm_csv = TimberHarvestPlanMosaic$`10cmDBH`,
      MosaicACD_csv = TimberHarvestPlanMosaic$ACD_10cm
    )
  ),
  # from R/timber-production-combine.R
  tar_target(
    combine_timber_production,
    combine_timber_production_files(
      tp10_19 = CF$timb_prod,
      tp20_21 = CF$timb_prod_MR
    ),
    format = "file"
  ),

  #  from 'R/EstimateMarketLeakage.R'
  tar_target(
    EstimateLeakage,
    estimate_market_leakage(
      tp = combine_timber_production,
      FRs_src = CF$forest_reserves_src,
      FRnames_src = CF$forest_res_names,
      CAO_Sabah_src = CF$SabahC_acd,
      Boundary_src = CF$kuamut_bounds,
      Mosaic_src = DefineStrata$MosaicCoupesGo,
      NFMCoupes_src = DefineStrata$NFMCoupesGo,
      Mosaic_Parcels_csv = TimberHarvestPlanMosaic$`31cmDBH`,
      NFM_Parcels_csv = TimberHarvestPlanNFM$`40cmDBH`
    )
  ),

  #  from 'R/ActivityShiftingLeakage.R'
  tar_target(
    ActivityShiftingLeakage,
    estimate_act_shift_leakage(
      tp = combine_timber_production
    )
  ),

  # from `R/proj-growth-rate-multi.R`
  tar_target(
    calculate_bio_growth_rate,
    calculate_proj_growth(
      agb_rasters = acd_path_list[1:7], # 2016 - 2022
      aoi_path = CF$kuamut_bounds,
      deforest_path = bind_auto_manual_deforest,
      mgmt_strat_path = DefineStrata$mgmt_strata_union,
      workers = 3
    )
  ),
  tar_target(
    calculate_bio_growth_rate_mr2,
    calculate_proj_growth(
      agb_rasters = acd_path_list[1:8], # including 2023
      aoi_path = CF$kuamut_bounds,
      deforest_path = bind_auto_manual_deforest,
      mgmt_strat_path = DefineStrata$mgmt_strata_union,
      workers = 3
    )
  ),

  # from 'R/uncertainty-analysis.R'
  tar_target(
    calculate_uncertainty,
    uncertainty_analysis(
      aoi = CF$kuamut_bounds,
      mv31 = MLF$MerchVol_31cm,
      mv40 = MLF$MerchVol_40cm,
      gr = list(calculate_bio_growth_rate, calculate_bio_growth_rate_mr2)
    )
  ),
  tar_target(
    create_pre_16_strata,
    split_strata(
      coupes_src = CF$kuamut_logging_coups,
      strata_src = DefineStrata$mgmt_strata_union
    )
  ),

  # from 'R/AnnualisedCalculations_VM10v1.3.R'
  tar_target(
    AnnualisedCalculations,
    annualised_calcs(
      Parcels_src = ParcelCalculations$ParcelStrata,
      mosaic_src = DefineStrata$MosaicCoupesGo,
      nfm_coupes_src = DefineStrata$NFMCoupesGo,
      deforest_calcs = deforest_emission_calcs,
      C_RG_planted = MosaicAreaGrowthRate$AverageGrowthRate,
      nfm_lfme = EstimateLeakage$nfm_LFME,
      mosaic_lfme = EstimateLeakage$mosaic_LFME,
      prj_ongo_gr = list(
        calculate_bio_growth_rate$growth_rates,
        calculate_bio_growth_rate_mr2$growth_rates
      ),
      mgmt_strata = create_pre_16_strata,
      mosaic_spcf = sp_conserv_factor$spc_mosaic,
      nfm_spcf = sp_conserv_factor$spc_nfm,
      mr_periods = list(c(2016, 2021), c(2022, 2023))
    )
  ),

  # from "R/KuamutTree2Plot.R"
  tar_target(
    plot_level_growth_rate,
    forest_plot_growth_rates(TreeToPlotMetrics)
  ),
  tar_target(
    fgb_to_shp,
    spat_exp_to_shp(list(
      kuamut_bounds = CF$kuamut_bounds,
      kuamut_deforest = bind_auto_manual_deforest
    ))
  )
)

LEAF_MAPS <- list(
  tar_target(
    site_leaf,
    tm_site_leaf(CF$kuamut_bounds)
  ),
  tar_target(
    DEM_leaf,
    tm_raster_leaf(MergeDEM)
  ),
  tar_target(
    slope_leaf,
    tm_flat_areas_leaf(
      slope_ras = SlopeAnalysis$slope_raster,
      flat_areas = SlopeAnalysis$flat_areas
    )
  ),
  tar_target(
    strata_leaf,
    tm_strata_leaf(
      aoi = CF$kuamut_bounds,
      nfm_coupes = DefineStrata$NFMCoupesGo,
      mosaic_plant = DefineStrata$MosaicPLANT,
      mosaic_noplant = DefineStrata$MosaicUNplanted,
      nogo = DefineStrata$NoGo,
      flat = SlopeAnalysis$flat_areas
    )
  ),
  tar_target(
    plot_locs_leaf,
    tm_plot_locs_leaf(
      aoi = CF$kuamut_bounds,
      plot_locs = SlopeCorrection$slp_correct_spat
    )
  ),
  tar_target(
    slp_correct_leaf,
    tm_slope_cor_leaf(
      slope_ras = SlopeCorrection$detail_slp_ras,
      plot_locs = SlopeCorrection$slp_correct_spat
    )
  ),
  tar_target(
    deforest_zones_leaf,
    deforest_leaf(
      aoi = CF$kuamut_bounds,
      deforest_zones = bind_auto_manual_deforest
    )
  ),
  tar_target(
    volume_claim_leaf,
    claimed_areas(
      aoi = CF$kuamut_bounds,
      nfm_vol = TimberHarvestPlanNFM$nfm_vol31_spv,
      mos_unplant_vol = TimberHarvestPlanMosaic$mos_unplant_spv,
      mos_plant_vol = TimberHarvestPlanMosaic$mos_plant_spv,
      nogo = DefineStrata$NoGo
    )
  ),
  tar_target(
    reference_area_leaf,
    ref_area_leaf(
      ref_areas = CF$reference_areas_src,
      aoi = CF$kuamut_bounds
    )
  )
)

GGPLOT_TARGETS <- list(
  tar_target(
    plot_acd_viz,
    plot_acd_raincloud(TreeToPlotMetrics$plot_level_spat)
  ),
  tar_target(
    bio_growth_rate_ggplot,
    fix_effects_plot(
      calculate_bio_growth_rate$fix_effs_long,
      calculate_bio_growth_rate$growth_stats,
      fig_path = set_path(
        "FixEff-GrowthRate-multi-mr1.png",
        parent = "data_out/MR_figures"
      )
    )
  ),
  tar_target(
    bio_growth_rate_ggplot_mr2,
    fix_effects_plot(
      calculate_bio_growth_rate_mr2$fix_effs_long,
      calculate_bio_growth_rate_mr2$growth_stats,
      fig_path = set_path(
        "FixEff-GrowthRate-multi-mr2.png",
        parent = "data_out/MR_figures"
      )
    )
  ),
  tar_target(
    timber_harvest_rate_eval,
    evaluate_timber_harvest_rate(
      combine_timber_production,
      TimberHarvestPlanMosaic$`31cmDBH`
    )
  )
)

CREATE_GH_RELEASE <- tar_skip(
  github_release,
  create_github_release(
    DATA_VERSION, ML_VERSION_mr1, ML_VERSION_mr2, results_version,
    tag = make_tag(release_type),
    body = gh_release_body,
    after_target = c(AnnualisedCalculations, build_book, build_book_minimal)
  ),
  skip = !make_gh_release
)

BUILD_REFS <- tar_force(build_refs,
  command = {
    pkg_bib <- "src/packages.bib"
    journ_bib <- "src/add.bib"
    ref_bib <- "src/references.bib"
    knitr::write_bib(
      c(
        .packages(), "knitr", "quarto"
      ),
      pkg_bib
    )

    .l <- readLines(pkg_bib)

    file_conn <- file(ref_bib)
    writeLines(.l, file_conn)
    close(file_conn)

    .pl <- readLines(journ_bib)
    write(.pl, ref_bib, append = TRUE)
  },
  force = TRUE
)

BUILD_BOOK <- tarchetypes::tar_quarto(build_book, profile = "code")

BUILD_BOOK_MINIMAL <- tarchetypes::tar_quarto(build_book_minimal,
  profile = "minimal"
)

1.2 R Package Dependencies

This workflow relies on a number of R packages to perform various tasks and analyses. These packages have been carefully selected to ensure the project’s functionality and reproducibility. The following packages were used:

targets, tarchetypes, renv, checkmate, rlang, here, tidyr, measurements, purrr, furrr, stringr, units, readr, data.table, piggyback, fs, glue, usethis, BIOMASS, rmapshaper, sf, exactextractr, terra, ezwarp, demsrc, Rmisc, viridis, ggplot2, scico, ggpmisc, ggtext, ggdist, patchwork, DT, visNetwork, httpgd, gt, kableExtra, tmap, leafem, quarto, downlit, languageserver, remotes, optparse, rstudioapi, zip, lme4, broom, dplyr

For full citations to these packages, please see the Chapter 23