Revamped sim function to tackle degradation Cstock at Cstock level an…

…d not E level: fct_combine_mcs_cstock() instead of fct_combine_mcs_E()
openforis · Dec 4, 2024 · 05d3b34 · 05d3b34
1 parent c5157de
commit 05d3b34
Show file tree

Hide file tree

Showing 15 changed files with 243 additions and 103 deletions.
diff --git a/.gitignore b/.gitignore
@@ -4,6 +4,7 @@
 vignettes/*.png
 vignettes/*.csv
 inst/extdata/example1_conf.xlsx
+inst/extdata/example2_conf.xlsx
 
 # History files
 .Rhistory

diff --git a/R/fct_calc_res.R b/R/fct_calc_res.R
@@ -43,6 +43,6 @@ fct_calc_res <- function(
       E_ME  = round((.data$E_ciupper - .data$E_cilower) / 2),
       E_U   = round(.data$E_ME / .data$E * 100),
     ) %>%
-    dplyr::select(!!col_id, .data$E, .data$E_U, .data$E_ME, .data$E_cilower, .data$E_ciupper)
+    dplyr::select(!!col_id, "E", "E_U", "E_ME", "E_cilower", "E_ciupper")
 
 }
diff --git a/R/fct_check_data2.R b/R/fct_check_data2.R
@@ -50,8 +50,8 @@
 #'   cat_racti  = c("DF", "DG", "EN", "EN_AF", "EN_RE"),
 #'   cat_ptype  = c("REF", "REF[0-9]", "MON", "MON[0-9]"),
 #'   cat_pdf    = c("normal", "beta"),
-#'   col_usr    = c("trunc_pdf", "n_iter", "ran_seed", "c_unit", "c_fraction", "dg_pool",
-#'                  "ad_annual", "conf_level"),
+#'   col_usr    = c("trunc_pdf", "n_iter", "ran_seed", "c_unit", "c_fraction", "c_fraction_se",
+#'                  "c_fraction_pdf", "dg_ext", "dg_pool", "ad_annual", "conf_level"),
 #'   col_time   = c("period_no", "year_start", "year_end", "period_type"),
 #'   col_ad     = c("trans_no",	"trans_id",	"trans_period",	"redd_activity", "lu_initial_id",
 #'                  "lu_initial",	"lu_final_id", "lu_final", "trans_area", "trans_se",
@@ -60,7 +60,7 @@
 #'                	"c_se",	"c_pdf",	"c_pdf_a",	"c_pdf_b",	"c_pdf_c")
 #' )
 #'
-#' app_checklist$cat_cpools_all <- c(app_checklist$cat_cpools, "RS", "DG_ratio")
+#' app_checklist$cat_cpools_all <- c(app_checklist$cat_cpools, "RS", "DG_ratio", "C_all")
 #'
 #' fct_check_data2(.ad = ad, .cs = cs, .usr = usr, .time = time, .checklist = app_checklist)
 #'
@@ -79,7 +79,7 @@ fct_check_data2 <- function(.usr, .time, .ad, .cs, .checklist){
   tmp <- list()
   out <- list()
 
-  ## Check 1. tables have at least the correct columns
+  ## Check 1. tables have at least the correct columns #########################
   tmp$cols_usr_ok  <- all(.checklist$col_usr %in% names(.usr))
   tmp$cols_time_ok <- all(.checklist$col_time %in% names(.time))
   tmp$cols_ad_ok   <- all(.checklist$col_ad %in% names(.ad))
@@ -100,7 +100,7 @@ fct_check_data2 <- function(.usr, .time, .ad, .cs, .checklist){
     out$cols_pb <- c(tmp$cols_usr_pb, tmp$cols_time_pb, tmp$cols_ad_pb, tmp$cols_cs_pb)
   }
 
-  ## Check 2. tables dimensions
+  ## Check 2. tables dimensions ################################################
   tmp$size_usr_ok  <- nrow(.usr) == 1 ## usr has only one row
   tmp$size_time_ok <- nrow(.time) >= 2 ## at least one ref and one monitoring
   tmp$size_ad_ok   <- nrow(.ad) >= 2 ## at least one lu transition for ref and monitoring
@@ -122,7 +122,7 @@ fct_check_data2 <- function(.usr, .time, .ad, .cs, .checklist){
   }
 
 
-  ## Check 3. data types are correct
+  ## Check 3. data types are correct ###########################################
   ## - usr tab
   tmp$datatypes_usr_ok <- all(
     is.logical(.usr$trunc_pdf),
@@ -197,7 +197,7 @@ fct_check_data2 <- function(.usr, .time, .ad, .cs, .checklist){
     out$datatypes_pb <- c(tmp$datatypes_usr_pb, tmp$datatypes_time_pb, tmp$datatypes_ad_pb, tmp$datatypes_cs_pb)
   }
 
-  ## Check 4. category variables are correct
+  ## Check 4. category variables are correct ###################################
   ## - usr
   ## Get usr$dg_pool as vector
   ## ex. c("AGB", "BGB", "DW") %in% c("AGB", "BGB", "DW", "LI", "SOC", "ALL")
@@ -207,7 +207,7 @@ fct_check_data2 <- function(.usr, .time, .ad, .cs, .checklist){
 
   tmp$cats_usr_ok <- all(
     .usr$c_unit %in% .checklist$cat_cunits,
-    all(dg_pool %in% .checklist$cat_cpools)
+    all(dg_pool %in% .checklist$cat_cpools_all)
   )
 
   ## - time
@@ -249,7 +249,8 @@ fct_check_data2 <- function(.usr, .time, .ad, .cs, .checklist){
   }
 
 
-  ## 5. Check Unique IDs
+  ## Check 5. Check Unique IDs #################################################
+
   tmp$ids_time_ok <- nrow(.time) == length(unique(.time$period_no))
   tmp$ids_ad_ok   <- nrow(.ad) == length(unique(.ad$trans_id))
   tmp$ids_cs_ok   <- nrow(.cs) == length(unique(.cs$c_id))
@@ -266,7 +267,8 @@ fct_check_data2 <- function(.usr, .time, .ad, .cs, .checklist){
     out$ids_pb <- c(tmp$ids_time_pb, tmp$ids_ad_pb, tmp$ids_cs_pb)
   }
 
-  ## 6. Check matching and logical interactions
+
+  ## Check 6. Check matching and logical interactions ##########################
 
   ## - Period matching exactly between tables
   tmp$match_period_ad_ok <- all(sort(unique(.ad$trans_period)) == sort(unique(.time$period_no)))
@@ -282,7 +284,7 @@ fct_check_data2 <- function(.usr, .time, .ad, .cs, .checklist){
   lu_ad <- sort(c(unique(.ad$lu_initial_id), unique(.ad$lu_final_id)))
   lu_cs <- sort(unique(.cs$lu_id))
 
-  tmp$match_lu_ok <- all(lu_ad == lu_cs)
+  tmp$match_lu_ok <- all(lu_ad %in% lu_cs)
 
   ## - At least one ref and one monitoring period
   nb_ref <- .time %>% dplyr::filter(stringr::str_detect(.data$period_type, pattern = "REF|REF[0-9]"))
@@ -292,7 +294,22 @@ fct_check_data2 <- function(.usr, .time, .ad, .cs, .checklist){
   tmp$match_mon_ok <- nrow(nb_mon) > 0
 
   ## - CF if DM
-  ## TBD should add carbon fraction as row in c_stocks instead of user_inputs
+  tmp$match_dm_ok <- (.usr$c_unit == "DM" & is.numeric(.usr$c_fraction)) | .usr$c_unit == "C"
+
+  ## - DEG ext working
+  if (!is.na(.usr$dg_ext) & "DG_ratio" %in% unique(.cs$c_pool)) {
+    dg_lu <- .cs |>
+      dplyr::filter(.data$c_pool == "DG_ratio") |>
+      dplyr::pull("lu_id") |>
+      stringr::str_remove(pattern = .usr$dg_ext)
+
+    tmp$match_dg_ok <- all(dg_lu %in% unique(.cs$lu_id))
+
+  } else if (is.na(.usr$dg_ext) & "DG_ratio" %in% unique(.cs$c_pool)) {
+    tmp$match_dg_ok <- FALSE
+  } else {
+    tmp$match_dg_ok <- FALSE
+  }
 
   ## - DG method: either (1) dg_ratio applied to all pools, (2) dg_ratio applied to some pools other kept intact (dg_expool == T), (3) diff in Cstocks.
   ## NOT IMPLEMENTED YET
@@ -303,7 +320,9 @@ fct_check_data2 <- function(.usr, .time, .ad, .cs, .checklist){
     tmp$match_period_cs_ok,
     tmp$match_lu_ok,
     tmp$match_ref_ok,
-    tmp$match_mon_ok
+    tmp$match_mon_ok,
+    tmp$match_dm_ok,
+    tmp$match_dg_ok
   )
 
   if (out$matches_ok) {
@@ -314,8 +333,9 @@ fct_check_data2 <- function(.usr, .time, .ad, .cs, .checklist){
     if (tmp$match_lu_ok)        tmp$match_lu_pb        <- NULL else tmp$match_lu_pb        <- "land use mismatch between AD and CS"
     if (tmp$match_ref_ok)       tmp$match_ref_pb       <- NULL else tmp$match_ref_pb       <- "missing REF in time_periods"
     if (tmp$match_mon_ok)       tmp$match_mon_pb       <- NULL else tmp$match_mon_pb       <- "missing MON in time_periods"
-
-    out$matches_pb <- c(tmp$match_period_ad_pb, tmp$match_period_cs_pb, tmp$match_lu_pb, tmp$match_ref_pb, tmp$match_mon_pb)
+    if (tmp$match_dm_ok)        tmp$match_dm_pb        <- NULL else tmp$match_dm_pb        <- "user_inputs c_units is 'DM' but c_fraction missing"
+    if (tmp$match_dg_ok)        tmp$match_dg_pb        <- NULL else tmp$match_dg_pb        <- "missing dg_ext for DG_ratio or mismatch with dg_ext and AD"
+    out$matches_pb <- c(tmp$match_period_ad_pb, tmp$match_period_cs_pb, tmp$match_lu_pb, tmp$match_ref_pb, tmp$match_mon_pb, tmp$match_dm_pb, tmp$match_dg_pb)
 
   }
 

diff --git a/R/fct_combine_mcs_E.R b/R/fct_combine_mcs_E.R
@@ -76,7 +76,8 @@ fct_combine_mcs_E <- function(.ad, .cs, .usr){
   mcs_trans <- purrr::map(vec_trans, function(x){
 
     ## !! FOR TESTING ONLY
-    # x = "T1_ev_wet_closed_dg_ev_wet_closed"
+    # x = "T1_H_H_deg" #"T1_P_Crop" #"T1_ev_wet_closed_dg_ev_wet_closed"
+    # print(x)
     ## !!
 
     ad_x   <- .ad %>% dplyr::filter(.data$trans_id == x)
@@ -127,44 +128,46 @@ fct_combine_mcs_E <- function(.ad, .cs, .usr){
 
     }
 
-    ## If degradation is ratio, using .usr$dg_pool to calculate C_all_f
-    if (redd_x == "DG" & !is.na(.usr$dg_pool)) {
-
-      dg_pool <- stringr::str_split(.usr$dg_pool, pattern = ",") |> purrr::map(stringr::str_trim) |> unlist()
-      dg_pool_i <- paste0(dg_pool, "_i")
-
-      combi <- combi %>%
-        dplyr::rowwise() %>%
-        dplyr::mutate(C_all_f = .data$DG_ratio_f * sum(!!!rlang::syms(dg_pool_i)) * 44/12) %>%
-        dplyr::ungroup()
-
-      ## If degradation has unaffected pools, we identify them by difference and add them to final C stock
-      c_pools <- .cs |>
-        dplyr::filter(.data$lu_id == ad_x$lu_initial_id) |>
-        dplyr::filter(!(is.na(.data$c_value) & is.na(.data$c_pdf_a))) |>
-        dplyr::pull("c_pool") |>
-        unique()
-      dg_expool <- paste0(setdiff(c_pools, dg_pool), "_i")
-
-      if (length(dg_expool) > 0) {
-        combi <- combi %>%
-          dplyr::rowwise() %>%
-          dplyr::mutate(C_all_f = .data$C_all_f + sum(!!!rlang::syms(dg_expool))) %>%
-          dplyr::ungroup()
-      }
-
-    }
+    # ## If degradation is ratio, using .usr$dg_pool to calculate C_all_f
+    # if (redd_x == "DG" & !is.na(.usr$dg_pool)) {
+    #
+    #   dg_pool <- stringr::str_split(.usr$dg_pool, pattern = ",") |> purrr::map(stringr::str_trim) |> unlist()
+    #   dg_pool_i <- paste0(dg_pool, "_i")
+    #
+    #   combi <- combi %>%
+    #     dplyr::rowwise() %>%
+    #     dplyr::mutate(C_all_f = .data$DG_ratio_f * sum(!!!rlang::syms(dg_pool_i))) %>%
+    #     dplyr::ungroup()
+    #
+    #   ## If degradation has unaffected pools, we identify them by difference and add them to final C stock
+    #   if (.usr$dg_pool != "C_all") {
+    #     c_pools <- .cs |>
+    #       dplyr::filter(.data$lu_id == ad_x$lu_initial_id) |>
+    #       dplyr::filter(!(is.na(.data$c_value) & is.na(.data$c_pdf_a))) |>
+    #       dplyr::pull("c_pool") |>
+    #       unique()
+    #     dg_expool <- paste0(setdiff(c_pools, dg_pool), "_i")
+    #
+    #     if (length(dg_expool) > 0) {
+    #       combi <- combi %>%
+    #         dplyr::rowwise() %>%
+    #         dplyr::mutate(C_all_f = .data$C_all_f + sum(!!!rlang::syms(dg_expool))) %>%
+    #         dplyr::ungroup()
+    #     }
+    #   }
+#
+#     }
 
     combi
 
   }) |> purrr::list_rbind()
   ## END LOOP
 
   ## Re-arrange columns and add EF and E (emissions at transition level)
-  mcs_trans %>%
+  tt <- mcs_trans %>%
     dplyr::mutate(
-      EF = .data$C_all_i - .data$C_all_f,
-      E_sim  = .data$AD * .data$EF
+      EF = round((.data$C_all_i - .data$C_all_f) * 44/12, 0),
+      E_sim  = round(.data$AD * .data$EF, 0)
     ) %>%
     # dplyr::mutate(dplyr::across(c(.data$E_sim, .data$AD, .data$EF, .data$C_all_i, .data$C_all_f))) |>
     dplyr::select(

diff --git a/R/fct_combine_mcs_cpools.R b/R/fct_combine_mcs_cpools.R
@@ -42,7 +42,7 @@
 fct_combine_mcs_cpools <- function(.c_sub, .usr){
 
   ## !! FOR TESTING ONLY
-  # .c_sub  <- cs |> filter(lu_id == "postdef_open") ## "dg_ev_wet_closed"
+  # .c_sub  <- cs |> dplyr::filter(lu_id == "P_deg") ## "dg_ev_wet_closed" ## "ALL_P_AGB"
   # .usr    <- usr
   ## !!
 
@@ -79,14 +79,13 @@ fct_combine_mcs_cpools <- function(.c_sub, .usr){
   ## End map()
 
   ## ADD CF if needed
-  if (c_check$has_CF) {
-    params <- .c_sub %>% dplyr::filter(.data$c_pool == "CF")
+  if (c_check$has_CF & !("DG_ratio" %in% names(SIMS))) {
     SIMS$CF <- fct_make_mcs(
       .n_iter = .usr$n_iter,
-      .pdf    = params$c_pdf,
-      .mean   = params$c_value,
-      .se     = params$c_se,
-      .params = c(params$c_pdf_a, params$c_pdf_b, params$c_pdf_c),
+      .pdf    = .usr$c_fraction_pdf,
+      .mean   = .usr$c_fraction,
+      .se     = .usr$c_fraction_se,
+      #.params = c(params$c_pdf_a, params$c_pdf_b, params$c_pdf_c),
       .trunc  = .usr$trunc_pdf
     )
   }

diff --git a/R/fct_combine_mcs_cstock.R b/R/fct_combine_mcs_cstock.R
@@ -57,6 +57,10 @@ fct_combine_mcs_cstock <- function(.ad, .cs, .usr){
 
   mcs_c <- purrr::pmap(combi, function(lu, period){
 
+    ## !!! FOR TESTING ONLY
+    # lu = "P_deg"
+    # period = "ALL"
+    ## !!!
     c_sub <- .cs %>%
       dplyr::filter(.data$lu_id == lu, .data$c_period == period) %>%
       dplyr::filter(!(is.na(.data$c_value) & is.na(.data$c_pdf_a)))
@@ -69,7 +73,54 @@ fct_combine_mcs_cstock <- function(.ad, .cs, .usr){
     dplyr::select("sim_no", "c_period", "lu_id", "C_all", "C_form", dplyr::everything())
 
   ## CHECK
-  # mcs_c |> dplyr::filter(.data$sim_no == 1)
+  # tt <- mcs_c |> dplyr::filter(.data$sim_no == 1)
+
+  if ("DG_ratio" %in% unique(mcs_c$C_form)) {
+
+    ## Get pools used for DG
+    if (.usr$dg_pool == "ALL") {
+      dg_pool <- "C_all"
+    } else {
+      dg_pool <- stringr::str_split(.usr$dg_pool, pattern = ",") |> purrr::map(stringr::str_trim) |> unlist()
+    }
+    dg_pool_intact <- paste0(dg_pool, "_intact")
+
+    ## Filter DG to modify formula and recalculate
+    mcs_dg <- mcs_c |>
+      dplyr::filter(.data$C_form == "DG_ratio") |>
+      dplyr::mutate(
+        lu_intact = stringr::str_remove(.data$lu_id, pattern = .usr$dg_ext)
+      )
+
+    mcs_join <- mcs_c |>
+      dplyr::filter(lu_id %in% unique(mcs_dg$lu_intact)) |>
+      dplyr::select("sim_no", lu_intact = "lu_id", !!!rlang::syms(dg_pool))
+
+    names(mcs_join)[!(names(mcs_join) %in% c("sim_no", "lu_intact"))] <- dg_pool_intact
+
+    mcs_dg2 <- mcs_dg |>
+      dplyr::left_join(mcs_join, by = c("sim_no", "lu_intact")) |>
+      dplyr::rowwise() |>
+      dplyr::mutate(
+        C_form = paste0(.data$C_form, " * (", paste0(dg_pool_intact, collapse = " + "), ")"),
+        C_all  = .data$C_all * sum(!!!rlang::syms(dg_pool_intact))
+      ) |>
+      dplyr::ungroup() |>
+      dplyr::select(-"lu_intact", -dplyr::all_of(dg_pool_intact))
+
+    mcs_c2 <- mcs_c |>
+      dplyr::filter(.data$C_form != "DG_ratio") |>
+      dplyr::bind_rows(mcs_dg2) |>
+      dplyr::distinct()
+
+  } else {
+    mcs_c2 <- mcs_c
+  }
+
+  mcs_c2
+
+  ## Check
+  # tt <- mcs_c2 |> dplyr::filter(.data$sim_no == 1)
 
 }
 

diff --git a/R/fct_make_formula.R b/R/fct_make_formula.R
@@ -35,8 +35,8 @@
 #' @export
 fct_make_formula <- function(.c_check, .c_unit){
 
-  c_eq <- c("(", "(",  "AGB", " + ", "BGB", ")", " * ", "CF", " + ", "DW", " + ", "LI", " + ", "SOC", ")", " * ", "44/12")
-  names(c_eq) <- c("all_(", "cf_(", "AGB", "plus_bgb", "BGB", "cf_)", "times_cf", "CF", "plus_dw", "DW", "plus_li", "LI", "plus_soc", "SOC", "all_)", "times_mol", "mol")
+  c_eq <- c("(",  "AGB", " + ", "BGB", ")", " * ", "CF", " + ", "DW", " + ", "LI", " + ", "SOC")
+  names(c_eq) <- c("cf_(", "AGB", "plus_bgb", "BGB", "cf_)", "times_cf", "CF", "plus_dw", "DW", "plus_li", "LI", "plus_soc", "SOC")
 
   ## Handle dg_ratio
   if (.c_check$has_DG) return("DG_ratio")
@@ -75,7 +75,7 @@ fct_make_formula <- function(.c_check, .c_unit){
   }
 
   if (.c_check$has_AL){
-    c_eq_out <- "ALL * 44/12"
+    c_eq_out <- "ALL"
   }
 
   ## Output