Merge pull request #32 from KWB-R/rescale_and_distribute

Rescale and distribute

Author: hsonne

DESCRIPTION

@@ -19,6 +19,7 @@ RoxygenNote: 7.3.1
 Suggests: 
     covr,
     kwb.rect,
+    plumber,
     testthat (>= 3.0.0)
 Imports: 
     dplyr,

NAMESPACE

@@ -21,6 +21,8 @@ export(prepare_berlin_inputs)
 export(read_column_info)
 export(real_evapo_transpiration)
 export(run_rabimo)
+export(run_rabimo_with_measures)
+export(test_api)
 export(yearly_height_to_volume_flow)
 importFrom(dplyr,left_join)
 importFrom(kwb.utils,catAndRun)

R/distribute_measures.R

@@ -162,7 +162,8 @@ get_swale_connection_table <- function(unpaved_area_table, target, total_area)
   sealed_areas <- select_columns(unpaved_area_table, "sealed_area")
   corrected_scas <- select_columns(unpaved_area_table, "corr_sca")
   shares <- select_columns(unpaved_area_table, "corr_sca_frac")
-  to_swale_areas <- select_columns(unpaved_area_table, "corr_sca")
+  to_swale_areas <- select_columns(unpaved_area_table, "corr_sca") #???
+  zeros <- rep(0, nrow(unpaved_area_table))
 
   sealed_mean <- sum(sealed_areas) / total_area
 
@@ -177,14 +178,19 @@ get_swale_connection_table <- function(unpaved_area_table, target, total_area)
     shares > target
   }
 
-  ref_areas <- rep(0, length(sealed_areas))
+  ref_areas <- zeros
+  deltas <- zeros
+
   ref_areas[consider] <- if (to_increase) {
     sealed_areas[consider] - corrected_scas[consider]
   } else {
     to_swale_areas[consider]
   }
 
-  deltas <- get_weight(ref_areas) * total_diff_area
+  #deltas <- get_weight(ref_areas) * total_diff_area
+  if ((tot_area <- sum(ref_areas)) > 0) {
+    deltas <- ref_areas / tot_area * total_diff_area
+  }
 
   to_swale_areas <- corrected_scas + deltas
 

R/distribute_shares.R

@@ -18,15 +18,15 @@ distribute_shares <- function(
     length(target_value) == 1L
   )
 
-  total_diff_area <- sum(target_value * base_areas) - sum(partial_areas)
+  total_diff_area <- round(sum(target_value * base_areas) - sum(partial_areas),2)
 
   to_increase <- total_diff_area > 0
   shares <- ifelse(base_areas > 0, partial_areas / base_areas, 0)
 
   consider <- if (to_increase) {
-    shares < target_value
+    round(shares,2) < round(target_value,2)
   } else {
-    shares > target_value
+    round(shares,2) > round(target_value,2)
   }
 
   ref_areas <- rep(0, n_areas)

R/get_measure_stats.R

@@ -60,9 +60,9 @@ get_measure_stats <- function(blocks, reference_system = 2)
   }
 
   list(
-    green_roof = c(mean = mean_green_roof, max = max_green_roof),
-    unpaved = c(mean = mean_unpvd, max = max_unpvd),
-    to_swale = c(mean = mean_sca, max = max_sca)
+    green_roof = list(mean = mean_green_roof, max = max_green_roof),
+    unpaved = list(mean = mean_unpvd, max = max_unpvd),
+    to_swale = list(mean = mean_sca, max = max_sca)
   )
 }
 

R/rescale_target_values.R

@@ -0,0 +1,23 @@
+# rescale_target_values --------------------------------------------------------
+rescale_target_values <- function(new_targets, blocks)
+{
+  total_area <- sum(get_main_area(blocks))
+  total_roof_area <- sum(get_roof_area(blocks))
+  total_sealed_area <- sum(get_sealed_area(blocks))
+
+  green_roof_new <- select_elements(new_targets, "green_roof")
+  unpaved_new <- select_elements(new_targets, "unpaved")
+  to_swale_new <- select_elements(new_targets, "to_swale")
+
+  green_roof <- green_roof_new * total_area / total_roof_area
+  unpaved <- unpaved_new
+  to_swale <- to_swale_new * total_area / total_sealed_area
+
+  max_allowed_unpaved <- (total_area - total_roof_area)/total_area
+
+  stopifnot(green_roof <= 1)
+  stopifnot(unpaved <= max_allowed_unpaved)
+  stopifnot(to_swale <= 1)
+
+  list(green_roof = green_roof, unpaved = unpaved, to_swale = to_swale)
+}

R/run_rabimo.R

@@ -256,14 +256,22 @@ run_rabimo <- function(data, config, controls = define_controls())
   output_format <- control("output_format")
 
   result_data <- if (output_format == "abimo") {
+
     # Provide the same columns as Abimo does
     rename_columns(result_data_raw, name_mapping)
+
   } else if (output_format == "rabimo") {
-    remove_columns(result_data_raw, pattern = "_flow") %>%
-      remove_columns("total_runoff") %>%
-      move_columns_to_front(c("code", "total_area")) %>%
-      dplyr::rename_with(~gsub("total_","",.))
+
+    data.frame(
+      code = result_data_raw$code,
+      area = result_data_raw$total_area,
+      runoff = result_data_raw$total_surface_runoff,
+      infiltr = result_data_raw$total_infiltration,
+      evapor = result_data_raw$total_evaporation
+    )
+
   } else {
+
     clean_stop("controls$output_format must be either 'abimo' or 'rabimo'.")
   }
 

R/run_rabimo_with_measures.R

@@ -0,0 +1,30 @@
+#' Distribute Rainwater Management Measures and run R-Abimo
+#'
+#' @param blocks data frame of selected blocks (same columns as in
+#'   \code{\link{rabimo_inputs_2020}$data})
+#' @param measures list with elements \code{green_roof}, \code{unpaved},
+#'   \code{to_swale} representing the target percentages of the total areas
+#'   corresponding to each measure
+#' @param config configuration object, default:
+#'   \code{\link{rabimo_inputs_2020}$config}
+#' @export
+run_rabimo_with_measures <- function(
+    blocks,
+    measures,
+    config = kwb.rabimo::rabimo_inputs_2020$config
+)
+{
+  rescaled_targets <- rescale_target_values(
+    new_targets = measures,
+    blocks = blocks
+  )
+
+  run_rabimo(
+    distribute_measures(blocks = blocks, targets = rescaled_targets),
+    config = config)
+}
+
+# config_file = kwb.abimo::default_config()
+# config <- kwb.abimo:::read_config(file = safe_path(config_file))
+# new_config <- abimo_config_to_config(config, add_class_5 = TRUE)
+# config = prepare_config(new_config)

R/test_api.R

@@ -0,0 +1,8 @@
+#' Test the API to main functions using the plumber package and swagger
+#'
+#' @export
+test_api <- function()
+{
+  plumber_file <- system.file("scripts/plumber.R", package = "kwb.rabimo")
+  plumber::pr_run(plumber::pr(plumber_file))
+}

R/wabila.R

@@ -1,33 +1,41 @@
 # calculate_wabila_roof --------------------------------------------------------
-calculate_wabila_roof <- function(area, retention_height)
-{
+calculate_wabila_roof <- function(area, retention_height, return_code = TRUE) {
   # get climatic parameters from area
-  P <- area$prec_yr
-  ET_p <- area$epot_yr
+  P <- area["prec_yr"]
+  ET_p <- area["epot_yr"]
+  code <- area["code"]
 
   # runoff component
-  a <- 0.9115 + 0.00007063 * P - 0.000007498 * ET_p - 0.2063 *
-    log(retention_height + 1)
+  a <- 0.9115 + 0.00007063 * P - 0.000007498 * ET_p - 0.2063 * log(retention_height + 1)
 
   # infiltration component
-  g <- 0
+  g <- rep(0, length(P))  # if P is a vector, g will be a vector of zeros with the same length
 
-  #evaporation component
+  # evaporation component
   v <- 1 - a
 
-  c(a, g, v)
+  result <- setNames(c(a, g, v),
+                     c("surface_runoff", "infiltration", "evaporation"))
+
+  if (return_code){
+    return(cbind(code, result))
+  }
+
+  result
+
 }
 
 # calculate_wabila_green_roof --------------------------------------------------
-calculate_wabila_green_roof <- function(area, height, kf, w_diff)
+calculate_wabila_green_roof <- function(area, height, kf, w_diff, return_code = TRUE)
 {
   # h: depth, in cm
   # kf: permeability of soil, in mm/h
   # w_diff: difference between water holding capacity and wilting point, unitless
 
   # get climatic parameters from area
-  P <- area$prec_yr
-  ET_p <- area$epot_yr
+  P <- area["prec_yr"]
+  ET_p <- area["epot_yr"]
+  code <- area["code"]
 
   # runoff component
   a = -2.182 + 0.4293 * log(P) - 0.0001092 * P +  236.1/ET_p +
@@ -40,7 +48,14 @@ calculate_wabila_green_roof <- function(area, height, kf, w_diff)
   #evaporation component
   v <- 1 - a
 
-  c(a, g, v)
+  result <- setNames(c(a, g, v),
+                     c("surface_runoff", "infiltration", "evaporation"))
+
+  if (return_code){
+    return(cbind(code, result))
+  }
+
+  result
 }
 
 # calculate_wabila_paved -------------------------------------------------------
@@ -131,7 +146,7 @@ estimate_swale_area <- function(kf)
 
 
 # calculate_delta_mod ----------------------------------------------------------
-calculate_delta_mod <- function(results_mod_1, results_mod_2,
+calculate_delta_mod_v0 <- function(results_mod_1, results_mod_2,
                                 has_codes = TRUE,
                                 var_names = c("surface_runoff",
                                               "infiltration",
@@ -171,3 +186,45 @@ calculate_delta_mod <- function(results_mod_1, results_mod_2,
   return(delta_mod)
 
 }
+
+calculate_delta_mod <- function(results_mod_1, results_mod_2,
+                                has_codes = TRUE,
+                                col_patterns = c("runoff|infiltration|evaporation"),
+                                codes_name = "code"){
+
+  # force result objects to be data.frames
+  if(is.vector(results_mod_1)){
+    results_mod_1 <- as.data.frame(t(results_mod_1))
+    result_mod_2 <- as.data.frame(t(results_mod_2))
+  } else {
+    results_mod_1 <- as.data.frame(results_mod_1)
+    results_mod_2 <- as.data.frame(results_mod_2)
+  }
+
+  stopifnot(nrow(results_mod_1) == nrow(results_mod_2))
+
+  var_names_1 <- grep(pattern = col_patterns, x = names(results_mod_1), value = TRUE)
+  var_names_2 <- grep(pattern = col_patterns, x = names(results_mod_2), value = TRUE)
+
+  precipitation_1 <- rowSums(results_mod_1[var_names_1])
+  precipitation_2 <- rowSums(results_mod_2[var_names_2])
+  tolerance <- 1e-6
+  stopifnot(all(abs(precipitation_1 - precipitation_2) < tolerance))
+
+  delta_mod <- data.frame(
+    delta_mod = rowSums(
+      abs(results_mod_1[var_names_1] - results_mod_2[var_names_2])
+    ) * 0.5 / precipitation_1
+  )
+
+  if(has_codes){
+    stopifnot(
+      identical(results_mod_1[[codes_name]], results_mod_2[[codes_name]])
+    )
+    codes <- results_mod_1[[codes_name]]
+    return(cbind(code = codes, delta_mod = delta_mod))
+  }
+
+  delta_mod
+
+}