diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index a1eacb98..a5a5ce41 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -57,8 +57,13 @@ jobs: key: ${{ runner.os }}-r-${{ matrix.config.r }}-3-${{ hashFiles('depends.Rds') }} restore-keys: ${{ runner.os }}-r-${{ matrix.config.r }}-3- + - name: Install system dependencies on Linux + if: runner.os == 'Linux' + run: | + sudo apt-get -y install cmake libssl-dev libgdal-dev gdal-bin libgeos-dev libproj-dev libsqlite3-dev libudunits2-dev + - name: Install system dependencies - if: runner.os == 'Linux (no, skip this!)' + if: runner.os == 'Linux (no, skip!)' env: RHUB_PLATFORM: linux-x86_64-ubuntu-gcc run: | diff --git a/.gitignore b/.gitignore index 4418d92d..6893b526 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ .Ruserdata docs kwb.rabimo.Rproj +inst/doc diff --git a/DESCRIPTION b/DESCRIPTION index 139faad1..5cfd679b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: kwb.rabimo Title: R Implementation of Water Balance Model Abimo -Version: 2.0.0 +Version: 2.1.0.9000 Authors@R: c( person("Hauke", "Sonnenberg", , "hauke.sonnenberg@kompetenz-wasser.de", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-9134-2871")), @@ -15,18 +15,21 @@ URL: https://github.com/KWB-R/kwb.rabimo BugReports: https://github.com/KWB-R/kwb.rabimo/issues Encoding: UTF-8 Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Suggests: ggplot2, jsonlite, + knitr, plumber, + rmarkdown, testthat (>= 3.0.0) Imports: dplyr, kwb.utils (>= 0.15.0), magrittr, parallel, - rlang + rlang, + sf Remotes: github::kwb-r/kwb.utils Config/testthat/edition: 3 @@ -34,3 +37,4 @@ Depends: R (>= 3.5.0) LazyData: true LazyDataCompression: xz +VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 9fb34966..3c0177f5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,34 +1,41 @@ -# Generated by roxygen2: do not edit by hand - -export("%>%") -export(calculate_delta_w) -export(data_to_natural) -export(define_controls) -export(generate_rabimo_area) -export(get_measure_stats) -export(get_soil_properties) -export(read_column_info) -export(real_evapo_transpiration) -export(run_rabimo) -export(run_rabimo_with_measures) -export(test_plumber_api) -export(triangle_of_fractions) -importFrom(kwb.utils,catAndRun) -importFrom(kwb.utils,createAccessor) -importFrom(kwb.utils,defaultIfNULL) -importFrom(kwb.utils,getAttribute) -importFrom(kwb.utils,printIf) -importFrom(kwb.utils,renameAndSelect) -importFrom(kwb.utils,renameColumns) -importFrom(kwb.utils,selectColumns) -importFrom(kwb.utils,selectElements) -importFrom(kwb.utils,stopFormatted) -importFrom(kwb.utils,stringList) -importFrom(magrittr,"%>%") -importFrom(parallel,detectCores) -importFrom(parallel,makeCluster) -importFrom(parallel,parLapply) -importFrom(parallel,stopCluster) -importFrom(rlang,.data) -importFrom(stats,approx) -importFrom(utils,globalVariables) +# Generated by roxygen2: do not edit by hand + +export("%>%") +export(calculate_delta_w) +export(crop_box) +export(data_to_natural) +export(define_controls) +export(generate_rabimo_area) +export(get_measure_stats) +export(get_soil_properties) +export(read_column_info) +export(real_evapo_transpiration) +export(run_rabimo) +export(run_rabimo_with_measures) +export(test_plumber_api) +export(triangle_of_fractions) +importFrom(kwb.utils,catAndRun) +importFrom(kwb.utils,createAccessor) +importFrom(kwb.utils,defaultIfNULL) +importFrom(kwb.utils,getAttribute) +importFrom(kwb.utils,printIf) +importFrom(kwb.utils,renameAndSelect) +importFrom(kwb.utils,renameColumns) +importFrom(kwb.utils,selectColumns) +importFrom(kwb.utils,selectElements) +importFrom(kwb.utils,stopFormatted) +importFrom(kwb.utils,stringList) +importFrom(magrittr,"%>%") +importFrom(parallel,detectCores) +importFrom(parallel,makeCluster) +importFrom(parallel,parLapply) +importFrom(parallel,stopCluster) +importFrom(rlang,.data) +importFrom(sf,st_as_sf) +importFrom(sf,st_as_sfc) +importFrom(sf,st_bbox) +importFrom(sf,st_crop) +importFrom(sf,st_drop_geometry) +importFrom(sf,st_sfc) +importFrom(stats,approx) +importFrom(utils,globalVariables) diff --git a/NEWS.md b/NEWS.md index a8c358a9..33c5d83c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,17 @@ +# kwb.rabimo 2.1.0 (2025-07-21) + +Contains a tutorial in the form of a vignette. To view the vignette, install +the package with + +`remotes::install_github("kwb-r/kwb.rabimo@dev", build_vignettes = TRUE)` + +and open the vignette with + +`vignette("tutorial", package = "kwb.rabimo")` + +You can read the tutorial also on our GitHub page: +https://kwb-r.github.io/kwb.rabimo/articles/tutorial.html + # kwb.rabimo 2.0.0 (2025-05-20) - contains new data for Berlin in which road areas do not belong to blocks diff --git a/R/calculate_delta_w.R b/R/calculate_delta_w.R index 895e993c..8dcdf035 100644 --- a/R/calculate_delta_w.R +++ b/R/calculate_delta_w.R @@ -25,7 +25,9 @@ calculate_delta_w <- function( ) { #kwb.utils::assignPackageObjects("kwb.rabimo") - #columns_water_balance = c("runoff", "infiltr", "evapor");column_code = "code" + #columns_water_balance=c("runoff","infiltr","evapor");column_code="code";digits=1L + + urban <- remove_geo_column_if_required(urban) columns <- c(column_code, columns_water_balance) data_urban <- select_columns(urban, columns) @@ -40,9 +42,22 @@ calculate_delta_w <- function( delta_ws <- rowSums(abs(joined_urban - joined_natural)) / rowSums(joined_natural) * 100 / 2 - data.frame( + delta_w <- data.frame( code = joined[[column_code]], delta_w = unname(round(delta_ws, digits)), stringsAsFactors = FALSE ) + + if (is.null(geometry <- attr(urban, "geometry"))) { + delta_w + } else { + restore_geo_column_if_required( + delta_w, + # unfortunately, the [] selection removes the attribute "sf_column" + geometry = structure( + geometry[match(delta_w$code, urban$code)], + sf_column = attr(geometry, "sf_column") + ) + ) + } } diff --git a/R/data_to_natural.R b/R/data_to_natural.R index 02d22b1f..88dfc773 100644 --- a/R/data_to_natural.R +++ b/R/data_to_natural.R @@ -25,6 +25,8 @@ data_to_natural <- function(data, type = "undeveloped", veg_class = 50) # data <- kwb.rabimo::rabimo_inputs_2020$data; type = "undeveloped" # data <- kwb.rabimo::rabimo_inputs_2025$data; type = "undeveloped" + data <- remove_geo_column_if_required(data) + # Check whether data look as expected stop_on_invalid_data(data) @@ -45,13 +47,18 @@ data_to_natural <- function(data, type = "undeveloped", veg_class = 50) } else if (type == "horticultural") { "horticultural" } else { - stop("please provide a known natural scenario type: undeveloped, horticultural or forested") + clean_stop( + 'Please provide a known natural scenario type: "undeveloped", ', + '"horticultural" or "forested"' + ) } } - check_or_convert_data_types( + data <- check_or_convert_data_types( data, types = get_expected_data_type(), convert = TRUE ) + + restore_geo_column_if_required(data) } diff --git a/R/generate_rabimo_area.R b/R/generate_rabimo_area.R index 32736517..80898d40 100644 --- a/R/generate_rabimo_area.R +++ b/R/generate_rabimo_area.R @@ -1,6 +1,6 @@ # generate_rabimo_area --------------------------------------------------------- -#' Generate an area in R-ABIMO format with default values +#' Generate an area in R-Abimo format with default values #' #' All default values can be overridden by entering new key-value pairs. #' diff --git a/R/rabimo_inputs_2020.R b/R/rabimo_inputs_2020.R index 2b980c86..3dea9944 100644 --- a/R/rabimo_inputs_2020.R +++ b/R/rabimo_inputs_2020.R @@ -17,7 +17,7 @@ #' \item{`prec_s`}{Long-term average of annual precipitation within summer months (May to October) in mm (integer)} #' \item{`epot_yr`}{Long-term average of annual potential evapotranspiration in mm (integer)} #' \item{`epot_s`}{Long-term average of annual potential evapotranspiration within summer months (May to October) in mm (integer)} -#' \item{`district`}{Number of Berlin "Bezirk" (district) in which the block area is located (character)} +#' \item{`district`}{Number of Berlin "Bezirk" (district) in which the block area is located (character). This column is Berlin-specific and optional, i.e. not required by the model.} #' \item{`total_area`}{Total block area in square metres (numeric)} #' \item{`main_frac`}{Fraction of the total area that is NOT considered as "road" area (numeric value between 0.0 and 1.0). This value should be 0.0 if roads are modelled separately, i.e. as block areas on their own.} #' \item{`roof`}{Fraction of the total area that is considered as "roof" area (numeric value between 0.0 and 1.0)} @@ -47,7 +47,7 @@ #' \item{`block_type`}{Block type identifier of the form "usage-type-id_block-type-id_usage-type-description_block-type-description" (character)} #' } #' -#' Note 1: The sum of surface class fractions `srf1_pvd`, `srf1_pvd`, `srf1_pvd`, `srf1_pvd`, `srf1_pvd` should be 1.0 within each block area. +#' Note 1: The sum of surface class fractions `srf1_pvd`, `srf2_pvd`, `srf3_pvd`, `srf4_pvd`, `srf5_pvd` should be 1.0 within each block area. #' #' Note 2: The fields with suffix "_r" are all zero because rows are modelled as their own blocks. In an earlier version of the dataset, roads were modelled as parts of the block area. #' diff --git a/R/rabimo_inputs_2025.R b/R/rabimo_inputs_2025.R index d6e9b33a..b313a7c8 100644 --- a/R/rabimo_inputs_2025.R +++ b/R/rabimo_inputs_2025.R @@ -9,18 +9,18 @@ #' @format ## `rabimo_inputs_2025` #' A list containing two elements: #' \describe{ -#' \item{data}{a data frame with the input data in r-abimo format ...(number of vars)} -#' \item{config}{a list object with configuration data} +#' \item{data}{a data frame with the input data in R-Abimo format (see below)} +#' \item{config}{a list object with configuration data (see below)} #' } #' @format ## `rabimo_inputs_2025$data` -#' A data.frame with 58531 observations of 25 variables: +#' A data.frame with 58531 observations of 26 variables: #' \describe{ #' \item{`code`}{Unique block area identifier (character)} #' \item{`prec_yr`}{Long-term average of annual precipitation in mm (integer)} #' \item{`prec_s`}{Long-term average of annual precipitation within summer months (May to October) in mm (integer)} #' \item{`epot_yr`}{Long-term average of annual potential evapotranspiration in mm (integer)} #' \item{`epot_s`}{Long-term average of annual potential evapotranspiration within summer months (May to October) in mm (integer)} -#' \item{`district`}{Number of Berlin "Bezirk" (district) in which the block area is located (character)} +#' \item{`district`}{Number of Berlin "Bezirk" (district) in which the block area is located (character). This column is Berlin-specific and optional, i.e. not required by the model.} #' \item{`total_area`}{Total block area in square metres (numeric)} #' \item{`roof`}{Fraction of the total area that is considered as "roof" area (numeric value between 0.0 and 1.0)} #' \item{`green_roof`}{Fraction of the roof area that belongs to green roofs (numeric value between 0.0 and 1.0). A value of 1.0 means that all roofs in the block area are green roofs.} @@ -32,21 +32,40 @@ #' \item{`srf3_pvd`}{Fraction of the paved area that belongs to surface class 3 (numeric value between 0.0 and 1.0, see note 1 below)} #' \item{`srf4_pvd`}{Fraction of the paved area that belongs to surface class 4 (numeric value between 0.0 and 1.0, see note 1 below)} #' \item{`srf5_pvd`}{Fraction of the paved area that belongs to surface class 5 (numeric value between 0.0 and 1.0, see note 1 below)} -#' \item{`to_swale`}{Fraction of sealed area (roof area + paved area) that is connected to an infiltration swale (numeric)} +#' \item{`to_swale`}{Fraction of sealed area (roof area + paved area) that is connected to an infiltration swale (numeric value between 0.0 and 1.0)} #' \item{`gw_dist`}{Distance between groundwater table and surface in metres (numeric)} #' \item{`ufc30`}{field capacity in 30 cm depth (numeric)} #' \item{`ufc150`}{field capacity in 150 cm depth (numeric)} #' \item{`land_type`}{land type, one of `forested`, `horticultural`, `urban`, `vegetationless`, `waterbody` (character)} -#' \item{`veg_class`}{vegetation class index (numeric), derived from an analysis tree volumes} +#' \item{`veg_class`}{vegetation class index (numeric), derived from an analysis of tree volumes} #' \item{`irrigation`}{irrigation in mm per year (integer)} -#' \item{`block_type`}{Block type identifier of the form "usage-type-id_block-type-id_usage-type-description_block-type-description" (character)} +#' \item{`block_type`}{Block type identifier of the form "usage-type-id_block-type-id_usage-type-description_block-type-description" (character). This column is Berlin-specific and optional, i.e. not required by the model.} +#' \item{`Shape`}{List structure containing geometry information on the different block areas. This column is optional. If provided, it will be appended to the model output so that model results can be plotted in the form of maps.} #' } +#' +#' Note 1: The sum of surface class fractions `srf1_pvd`, `srf2_pvd`, `srf3_pvd`, `srf4_pvd`, `srf5_pvd` should be 1.0 within each block area. +#' #' @format ## `rabimo_inputs_2025$config` #' A list with 3 named elements: #' \describe{ -#' \item{runoff_factors}{Runoff factors, vector of numeric with names `roof`, `surface1`, `surface2`, `surface3`, `surface4`, `surface5`} -#' \item{bagrov_values}{Bagrov values for sealed surfaces, vector of numeric with names `roof`, `green_roof`, `surface1`, `surface2`, `surface3`, `surface4`, `surface5`} -#' \item{swale}{Model parameter(s) related to the 'swale' measure, vector of numeric with currently one value, named `swale_evaporation_factor`} +#' \item{runoff_factors}{Runoff factors for roofs and five different surface +#' types, given as a vector of numeric with element names `roof`, `surface1`, +#' `surface2`, `surface3`, `surface4`, `surface5`. A runoff factor determines +#' the proportion of precipitation that, after subtraction of +#' evapotranspiration, becomes surface runoff from a paved area. The higher +#' the factor, the less permeable is the surface.} +#' \item{bagrov_values}{Bagrov values to calculate evapotranspiration from +#' paved surfaces, given as a vector of numeric with element names `roof`, +#' `green_roof`, `surface1`, `surface2`, `surface3`, `surface4`, `surface5`. +#' The higher the Bagrov value, the more evapotranspiration is generated by +#' the model. For a description of the evapotranspiration model and for a +#' figure that shows the influence of the Bagrov values (n) on the +#' evapotranspiration (in German), see \url{https://www.berlin.de/umweltatlas/wasser/wasserhaushalt/2001/methode/}} +#' \item{swale}{Model parameter(s) related to the 'swale' measure, given as a +#' vector of numeric with currently one value, named +#' `swale_evaporation_factor`. The swale evaporation factor determines which +#' fraction of the water going into a swale becomes evapotranspiration (the +#' rest becomes infiltration).} #' } #' @source #' @source diff --git a/R/run_rabimo.R b/R/run_rabimo.R index 186cc3ed..fca94ac6 100644 --- a/R/run_rabimo.R +++ b/R/run_rabimo.R @@ -1,436 +1,494 @@ -# run_rabimo ------------------------------------------------------------------- - -#' Run R-Abimo, the R-implementation of Water Balance Model Abimo -#' -#' @param data data frame similar to -#' \code{\link{rabimo_inputs_2025}$data} -#' @param config configuration object (list) similar to -#' \code{\link{rabimo_inputs_2025}$config} -#' @param controls list of settings that control how the function should behave. -#' Use \code{\link{define_controls}} to define such a list. The default is -#' the list returned by \code{define_controls()}. -#' @return data frame with columns as returned by Abimo -#' @export -#' @examples -#' # Get input data and config for Berlin (version 2020) -#' inputs_2020 <- kwb.rabimo::rabimo_inputs_2020 -#' -#' # Randomly select 1000 blocks (to reduce runtime) -#' data <- inputs_2020$data -#' data <- data[sample(seq_len(nrow(data)), size = 1000L), ] -#' -#' # Run R-Abimo -#' results_2020 <- kwb.rabimo::run_rabimo(data, inputs_2020$config) -#' -#' # Get input data and config for Berlin (version 2025) -#' inputs_2025 <- kwb.rabimo::rabimo_inputs_2025 -#' -#' # Randomly select 1000 blocks (to reduce runtime) -#' data <- inputs_2025$data -#' data <- data[sample(seq_len(nrow(data)), size = 1000L), ] -#' -#' # Run R-Abimo -#' results_2025 <- kwb.rabimo::run_rabimo(data, inputs_2025$config) -run_rabimo <- function(data, config, controls = define_controls()) -{ - # Provide functions and variables for debugging - # (Go to inst/scripts/test-rabimo.R to provide data and config for debugging) - if (FALSE) - { - kwb.utils::assignPackageObjects("kwb.rabimo") - inputs <- kwb.utils:::get_cached("rabimo_inputs_2020") - data <- inputs$data - config <- inputs$config - controls <- define_controls() - `%>%` <- magrittr::`%>%` - } - - # If road-area-specific columns are missing, create them - data <- handle_missing_columns(data) - - # Provide function to access the list of controls - control <- create_accessor(controls) - - # Check whether data and config have the expected structures - if (isTRUE(control("check"))) { - stop_on_invalid_data(data) - stop_on_invalid_config(config) - } - - # Get climate data - climate <- cat_and_run( - "Collecting climate related data", - get_climate(data) - ) - - # Create accessor functions to data columns and config elements - fetch_data <- create_accessor(data) - fetch_config <- create_accessor(config) - fetch_climate <- create_accessor(climate) - - # Prepare soil properties for all rows. They are required to calculate the - # actual evapotranspiration of unsealed areas. In the case of water bodies, - # all values are 0.0. (hsonne: really?) - soil_properties <- cat_and_run( - "Preparing soil property data for all block areas", - expr = get_soil_properties( - land_type = fetch_data("land_type"), - veg_class = fetch_data("veg_class"), - depth_to_water_table = fetch_data("gw_dist"), - field_capacity_30 = fetch_data("ufc30"), - field_capacity_150 = fetch_data("ufc150"), - dbg = FALSE - ) - ) - - # Pre-calculate all results of realEvapoTranspiration() - evaporation_sealed <- cat_and_run( - "Precalculating actual evapotranspirations for impervious areas", - expr = fetch_config("bagrov_values") %>% - lapply(function(x) { - real_evapo_transpiration( - potential_evaporation = fetch_climate("epot_yr"), - x_ratio = fetch_climate("x_ratio"), - bagrov_parameter = rep(x, nrow(data)), - use_abimo_algorithm = control("use_abimo_bagrov_solver") - ) - }) %>% - do.call(what = data.frame) - ) - - # Pre-calculate all results of actualEvaporationWaterbodyOrPervious() - evaporation_unsealed <- cat_and_run( - paste( - "Precalculating actual evapotranspirations for waterbodies or pervious", - "areas" - ), - actual_evaporation_waterbody_or_pervious( - usage_tuple = fetch_data(c("land_type", "veg_class", "irrigation")), - climate = climate, - soil_properties = soil_properties, - min_size_for_parallel = 100L, - #digits = 3L, - use_abimo_algorithm = control("use_abimo_bagrov_solver") - ) - ) - - runoff_all <- fetch_climate("prec_yr") - cbind( - evaporation_sealed, - unsealed = evaporation_unsealed - ) - - # Runoff for all sealed areas (including roofs) - - # Calculate roof related variables - - # total runoff of roof areas - # (total runoff, contains both surface runoff and infiltration components) - runoff_roof <- select_columns(runoff_all, "roof") - runoff_green_roof <- select_columns(runoff_all, "green_roof") - - # Provide runoff coefficients for impervious surfaces - runoff_factors <- fetch_config("runoff_factors") - - # actual runoff from roof surface (area based, with no infiltration) - runoff_roof_actual <- with(data, main_frac * - roof * (1 - green_roof) * swg_roof) * - runoff_factors[["roof"]] * runoff_roof - - # actual runoff from green roof surface (area based, with no infiltration) - runoff_green_roof_actual <- with(data, main_frac * - roof * green_roof * swg_roof) * - runoff_factors[["roof"]] * runoff_green_roof - - # actual infiltration from roof surface (area based, with no runoff) - infiltration_roof_actual <- with(data, main_frac * roof * - (1-green_roof) * (1-swg_roof)) * - runoff_roof - - # actual infiltration from green_roof surface (area based, with no runoff) - infiltration_green_roof_actual <- with(data, main_frac * roof * - green_roof * (1-swg_roof)) * - runoff_green_roof - - - # Calculate runoff for all surface classes at once - # (contains both surface runoff and infiltration components) - - # Identify active surface class columns in input data - surface_cols_no_rd <- matching_names(data, pattern_no_roads()) - surface_cols_rd <- matching_names(data, pattern_roads()) - digits <- gsub("\\D", "", surface_cols_no_rd) - surface_class_names <- paste0("surface",digits) - - # choose columns related to surface classes - runoff_sealed <- select_columns(runoff_all, surface_class_names) - # head(runoff_sealed) - - # Runoff from the actual partial areas that are sealed and connected - # (road and non-road) areas (for all surface classes at once) - - runoff_factor_matrix <- expand_to_matrix( - x = runoff_factors[surface_class_names], - nrow = nrow(data) - ) - - unbuilt_surface_fractions <- fetch_data(surface_cols_no_rd) - road_surface_fractions <- fetch_data(surface_cols_rd) - - # add an empty column in road_surface_fraction to match dimension if needed - if (!identical(length(surface_cols_no_rd), length(surface_cols_rd))) { - road_surface_fractions$srf5_pvd_r <- 0 - } - - runoff_sealed_actual <- runoff_sealed * ( - with(data, main_frac * pvd * swg_pvd) * unbuilt_surface_fractions + - with(data, road_frac * pvd_r * swg_pvd_r) * road_surface_fractions - ) * - runoff_factor_matrix - - # infiltration of sealed surfaces - # (road and non-road) areas (for all surface classes at once) - infiltration_sealed_actual <- runoff_sealed * ( - with(data, main_frac * pvd) * unbuilt_surface_fractions + - with(data, road_frac * pvd_r) * road_surface_fractions) - - runoff_sealed_actual - - # Total Runoff of unsealed surfaces (unsealedSurface_RUV) - runoff_unsealed <- fetch_climate("prec_yr") - as.numeric(evaporation_unsealed) # why as.numeric()? - - # Infiltration of road (unsealed areas) - infiltration_unsealed_roads <- - with(data, road_frac * (1 - pvd_r)) * - runoff_sealed[, ncol(runoff_sealed)] # last (less sealed) surface class - - fraction_unsealed <- with( - data, - ifelse(control("reproduce_abimo_error"), 1, main_frac) * (1 - (roof + pvd)) - ) - - infiltration_unsealed_surfaces <- fraction_unsealed * runoff_unsealed - - # Calculate runoff 'ROW' for entire block area (FLGES + STR_FLGES) (mm/a) - total_surface_runoff <- ( - runoff_roof_actual + runoff_green_roof_actual + - #orig.: runoff_unsealed_roads <- was set to zero in the master branch - rowSums(runoff_sealed_actual)) - - # Calculate infiltration rate 'RI' for entire block partial area (mm/a) - total_infiltration <- - (infiltration_roof_actual + - infiltration_green_roof_actual + - infiltration_unsealed_surfaces + - infiltration_unsealed_roads + - rowSums(infiltration_sealed_actual)) - - # Correct Surface Runoff and Infiltration if area has an infiltration swale - swale_delta <- total_surface_runoff * (fetch_data("to_swale")) - total_surface_runoff <- total_surface_runoff - swale_delta - total_infiltration <- total_infiltration + - swale_delta * (1 - fetch_config("swale")[["swale_evaporation_factor"]]) - - # Calculate "total system losses" 'R' due to runoff and infiltration - # for entire block partial area - total_runoff <- total_surface_runoff + total_infiltration - - # Calculate evaporation 'VERDUNST' by subtracting 'R', the sum of - # runoff and infiltration from precipitation of entire year, - # multiplied by precipitation correction factor - total_evaporation <- climate[["prec_yr"]] - total_runoff - - # Provide total area for calculation of "flows" - total_area <- fetch_data("total_area") - - # Calculate volume 'rowvol' from runoff (qcm/s) - surface_runoff_flow <- yearly_height_to_volume_flow( - total_surface_runoff, total_area - ) - - # Calculate volume 'rivol' from infiltration rate (qcm/s) - infiltration_flow <- yearly_height_to_volume_flow( - total_infiltration, total_area - ) - - # Calculate volume of "system losses" 'rvol' due to surface runoff and - # infiltration - total_runoff_flow <- surface_runoff_flow + infiltration_flow - - # Provide mapping between local variable names and ABIMO-output columns - name_mapping <- list( - code = "CODE", - total_runoff = "R", - total_surface_runoff = "ROW", - total_infiltration = "RI", - total_runoff_flow = "RVOL", - surface_runoff_flow = "ROWVOL", - infiltration_flow = "RIVOL", - total_area = "FLAECHE", - total_evaporation = "VERDUNSTUN" - ) - - # Compose result data frame. Use mget() to get the result vectors from the - # local environment and put them into the data frame - result_data_raw <- cbind( - fetch_data("code", drop = FALSE), - mget(names(name_mapping)[-1L]) - ) - - 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") { - - 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'.") - } - - # Round all columns to three digits (skip first column: "CODE") - result_data[-1L] <- lapply(result_data[-1L], round, 3L) - - if (isFALSE(control("intermediates"))) { - return(result_data) - } - - # Return intermediate results as attributes - structure( - result_data, - data = data, - intermediates = list( - climate = climate, - soil_properties = soil_properties, - evaporation_sealed = evaporation_sealed, - evaporation_unsealed = evaporation_unsealed, - roof = list( - evaporation_roof = evaporation_sealed[["roof"]], - runoff_roof = runoff_roof, - runoff_roof_actual = runoff_roof_actual, - infiltration_roof_actual = infiltration_roof_actual - ), - surface = list( - evaporation_sealed = evaporation_sealed[, -1L], - runoff_sealed = runoff_sealed, - runoff_sealed_actual = runoff_sealed_actual, - infiltration_sealed_actual = infiltration_sealed_actual - ), - runoff = cbind( - unsealedSurface = runoff_unsealed, - unsealedRoads = runoff_sealed[, ncol(runoff_sealed)], - sealed = runoff_sealed - ), - fraction_unsealed = fraction_unsealed - ) - ) -} - -# handle_missing_columns ------------------------------------------------------- -handle_missing_columns <- function(data) -{ - road_specific_columns <- c( - "road_frac", "pvd_r", "swg_pvd_r", - "srf1_pvd_r", "srf2_pvd_r", "srf3_pvd_r", "srf4_pvd_r" - ) - - missing_road_columns <- setdiff(road_specific_columns, names(data)) - - if (length(missing_road_columns)) { - for (column in missing_road_columns) { - data[[column]] <- 0 - } - } - - if (! "main_frac" %in% names(data)) { - data$main_frac <- 1 - } - - data -} - -# get_climate: provides climate relevant input data ---------------------------- -get_climate <- function(input) -{ - climate <- select_columns(input, c("prec_yr", "prec_s", "epot_yr", "epot_s")) - - climate[["x_ratio"]] <- climate[["prec_yr"]] / climate[["epot_yr"]] - - climate -} - -# yearly_height_to_volume_flow ------------------------------------------------- - -#' Convert Yearly Height (mm) to Volume Flow (unit?) -#' -#' @param height height in mm -#' @param area area in square metres -#' @keywords internal -yearly_height_to_volume_flow <- function(height, area) -{ - height * 3.171 * area / 100000.0 -} - -#' Define List of "Controls" -#' -#' Define a list of settings that control how the main function -#' \code{\link{run_rabimo}} should behave. -#' -#' @param check logical indicating whether the check functions are executed. -#' Default: \code{TRUE}. -#' @param use_abimo_bagrov_solver logical indicating whether or not to use the -#' (fast!) algorithm implemented in Abimo to solve the Bagrov equations. -#' Default: \code{TRUE}. -#' @param reproduce_abimo_error logical indicating whether or not to reproduce -#' the error that is contained in Abimo (missing area fraction factor). -#' Default: \code{FALSE}. -#' @param output_format one of "abimo" (upper case columns: CODE, R, ROW, RI, -#' RVOL, ROWVOL, RIVOL, FLAECHE, VERDUNSTUN), "rabimo" (lower case columns: -#' code, surface_runoff, infiltration, evaporation). -#' @param intermediates logical indicating whether the intermediate results are -#' returned as attributes. Default: \code{FALSE}. -#' @returns list with the arguments of this function as list elements -#' @export -#' @examples -#' inputs <- kwb.rabimo::rabimo_inputs_2020 -#' test_data <- inputs$data[sample(seq_len(nrow(inputs$data)), size = 1000L), ] -#' controls_default <- define_controls() -#' controls_no_check <- define_controls(check = FALSE) -#' controls_no_solver <- define_controls(use_abimo_bagrov_solver = FALSE) -#' system.time(result_default <- kwb.rabimo::run_rabimo( -#' test_data, inputs$config, controls_default -#' )) -#' system.time(result_no_check <- kwb.rabimo::run_rabimo( -#' test_data, inputs$config, controls_no_check -#' )) -#' identical(result_default, result_no_check) -#' \dontrun{ -#' system.time(result_no_solver <- kwb.rabimo::run_rabimo( -#' test_data, inputs$config, controls_no_solver -#' )) -#' } -define_controls <- function( - check = TRUE, - use_abimo_bagrov_solver = TRUE, - reproduce_abimo_error = FALSE, - output_format = "rabimo", - intermediates = FALSE -) -{ - list( - check = check, - use_abimo_bagrov_solver = use_abimo_bagrov_solver, - reproduce_abimo_error = reproduce_abimo_error, - output_format = output_format, - intermediates = intermediates - ) -} +# run_rabimo ------------------------------------------------------------------- + +#' Run R-Abimo, the R-implementation of Water Balance Model Abimo +#' +#' @param data data frame similar to +#' \code{\link{rabimo_inputs_2025}$data} +#' @param config configuration object (list) similar to +#' \code{\link{rabimo_inputs_2025}$config} +#' @param controls list of settings that control how the function should behave. +#' Use \code{\link{define_controls}} to define such a list. The default is +#' the list returned by \code{define_controls()}. +#' @param silent logical indicating whether to suppress console outputs, +#' the default is \code{FALSE} +#' @return data frame with columns as returned by Abimo +#' @importFrom sf st_as_sf st_drop_geometry st_sfc +#' @export +#' @examples +#' # Get input data and config for Berlin (version 2020) +#' inputs_2020 <- kwb.rabimo::rabimo_inputs_2020 +#' +#' # Randomly select 1000 blocks (to reduce runtime) +#' data <- inputs_2020$data +#' data <- data[sample(seq_len(nrow(data)), size = 1000L), ] +#' +#' # Run R-Abimo +#' results_2020 <- kwb.rabimo::run_rabimo(data, inputs_2020$config) +#' +#' # Get input data and config for Berlin (version 2025) +#' inputs_2025 <- kwb.rabimo::rabimo_inputs_2025 +#' +#' # Crop a box (to reduce runtime) +#' data <- crop_box(inputs_2025$data) +#' +#' # Run R-Abimo +#' results_2025 <- kwb.rabimo::run_rabimo(data, inputs_2025$config) +#' +#' plot(results_2025[, -1L]) +run_rabimo <- function( + data, config, controls = define_controls(), silent = FALSE +) +{ + # Provide functions and variables for debugging + # (Go to inst/scripts/test-rabimo.R to provide data and config for debugging) + if (FALSE) + { + kwb.utils::assignPackageObjects("kwb.rabimo") + data <- kwb.rabimo::rabimo_inputs_2025$data + config <- kwb.rabimo::rabimo_inputs_2025$config + controls <- define_controls() + `%>%` <- magrittr::`%>%` + } + + data <- remove_geo_column_if_required(data) + + # Save geometry data that may have stored in attribute "geometry" + geometry <- attr(data, "geometry") + + # If road-area-specific columns are missing, create them + data <- handle_missing_columns(data) + + # Provide function to access the list of controls + control <- create_accessor(controls) + + # Check whether data and config have the expected structures + if (isTRUE(control("check"))) { + stop_on_invalid_data(data) + stop_on_invalid_config(config) + } + + # Get climate data + climate <- cat_and_run( + dbg = !silent, + "Collecting climate related data", + get_climate(data) + ) + + # Create access functions to data columns and config elements + fetch_data <- create_accessor(data) + fetch_config <- create_accessor(config) + fetch_climate <- create_accessor(climate) + + # Prepare soil properties for all rows. They are required to calculate the + # actual evapotranspiration of unsealed areas. In the case of water bodies, + # all values are 0.0. (hsonne: really?) + soil_properties <- cat_and_run( + dbg = !silent, + "Preparing soil property data for all block areas", + expr = get_soil_properties( + land_type = fetch_data("land_type"), + veg_class = fetch_data("veg_class"), + depth_to_water_table = fetch_data("gw_dist"), + field_capacity_30 = fetch_data("ufc30"), + field_capacity_150 = fetch_data("ufc150"), + dbg = FALSE + ) + ) + + # Precalculate actual evapotranspirations for impervious areas + evaporation_sealed <- cat_and_run( + dbg = !silent, + "Precalculating actual evapotranspirations for impervious areas", + expr = fetch_config("bagrov_values") %>% + lapply(function(x) { + real_evapo_transpiration( + potential_evaporation = fetch_climate("epot_yr"), + x_ratio = fetch_climate("x_ratio"), + bagrov_parameter = rep(x, nrow(data)), + use_abimo_algorithm = control("use_abimo_bagrov_solver") + ) + }) %>% + do.call(what = data.frame) + ) + + # Precalculate actual evapotranspirations for waterbodies or pervious areas + evaporation_unsealed <- cat_and_run( + dbg = !silent, + paste( + "Precalculating actual evapotranspirations for waterbodies or pervious", + "areas" + ), + actual_evaporation_waterbody_or_pervious( + usage_tuple = fetch_data(c("land_type", "veg_class", "irrigation")), + climate = climate, + soil_properties = soil_properties, + min_size_for_parallel = 100L, + #digits = 3L, + use_abimo_algorithm = control("use_abimo_bagrov_solver") + ) + ) + + runoff_all <- fetch_climate("prec_yr") - cbind( + evaporation_sealed, + unsealed = evaporation_unsealed + ) + + # Runoff for all sealed areas (including roofs) + + # Calculate roof related variables + + # total runoff of roof areas + # (total runoff, contains both surface runoff and infiltration components) + runoff_roof <- select_columns(runoff_all, "roof") + runoff_green_roof <- select_columns(runoff_all, "green_roof") + + # Provide runoff coefficients for impervious surfaces + runoff_factors <- fetch_config("runoff_factors") + + # actual runoff from roof surface (area based, with no infiltration) + runoff_roof_actual <- with( + data, + main_frac * roof * (1 - green_roof) * swg_roof + ) * runoff_factors[["roof"]] * runoff_roof + + # actual runoff from green roof surface (area based, with no infiltration) + runoff_green_roof_actual <- with( + data, + main_frac * roof * green_roof * swg_roof + ) * runoff_factors[["roof"]] * runoff_green_roof + + # actual infiltration from roof surface (area based, with no runoff) + infiltration_roof_actual <- with( + data, main_frac * roof * (1-green_roof) * (1-swg_roof) + ) * runoff_roof + + # actual infiltration from green_roof surface (area based, with no runoff) + infiltration_green_roof_actual <- with( + data, + main_frac * roof * green_roof * (1-swg_roof) + ) * runoff_green_roof + + # Calculate runoff for all surface classes at once + # (contains both surface runoff and infiltration components) + + # Identify active surface class columns in input data + surface_cols_no_rd <- matching_names(data, pattern_no_roads()) + surface_cols_rd <- matching_names(data, pattern_roads()) + digits <- gsub("\\D", "", surface_cols_no_rd) + surface_class_names <- paste0("surface",digits) + + # choose columns related to surface classes + runoff_sealed <- select_columns(runoff_all, surface_class_names) + # head(runoff_sealed) + + # Runoff from the actual partial areas that are sealed and connected + # (road and non-road) areas (for all surface classes at once) + + runoff_factor_matrix <- expand_to_matrix( + x = runoff_factors[surface_class_names], + nrow = nrow(data) + ) + + unbuilt_surface_fractions <- fetch_data(surface_cols_no_rd) + road_surface_fractions <- fetch_data(surface_cols_rd) + + # add an empty column in road_surface_fraction to match dimension if needed + if (!identical(length(surface_cols_no_rd), length(surface_cols_rd))) { + road_surface_fractions$srf5_pvd_r <- 0 + } + + runoff_sealed_actual <- runoff_sealed * ( + with(data, main_frac * pvd * swg_pvd) * unbuilt_surface_fractions + + with(data, road_frac * pvd_r * swg_pvd_r) * road_surface_fractions + ) * + runoff_factor_matrix + + # infiltration of sealed surfaces + # (road and non-road) areas (for all surface classes at once) + infiltration_sealed_actual <- runoff_sealed * ( + with(data, main_frac * pvd) * unbuilt_surface_fractions + + with(data, road_frac * pvd_r) * road_surface_fractions) - + runoff_sealed_actual + + # Total Runoff of unsealed surfaces (unsealedSurface_RUV) + runoff_unsealed <- fetch_climate("prec_yr") - as.numeric(evaporation_unsealed) # why as.numeric()? + + # Infiltration of road (unsealed areas) + infiltration_unsealed_roads <- + with(data, road_frac * (1 - pvd_r)) * + runoff_sealed[, ncol(runoff_sealed)] # last (less sealed) surface class + + fraction_unsealed <- with( + data, + ifelse(control("reproduce_abimo_error"), 1, main_frac) * (1 - (roof + pvd)) + ) + + infiltration_unsealed_surfaces <- fraction_unsealed * runoff_unsealed + + # Calculate runoff 'ROW' for entire block area (FLGES + STR_FLGES) (mm/a) + total_surface_runoff <- ( + runoff_roof_actual + runoff_green_roof_actual + + #orig.: runoff_unsealed_roads <- was set to zero in the master branch + rowSums(runoff_sealed_actual)) + + # Calculate infiltration rate 'RI' for entire block partial area (mm/a) + total_infiltration <- + (infiltration_roof_actual + + infiltration_green_roof_actual + + infiltration_unsealed_surfaces + + infiltration_unsealed_roads + + rowSums(infiltration_sealed_actual)) + + # Correct Surface Runoff and Infiltration if area has an infiltration swale + swale_delta <- total_surface_runoff * (fetch_data("to_swale")) + total_surface_runoff <- total_surface_runoff - swale_delta + total_infiltration <- total_infiltration + + swale_delta * (1 - fetch_config("swale")[["swale_evaporation_factor"]]) + + # Calculate "total system losses" 'R' due to runoff and infiltration + # for entire block partial area + total_runoff <- total_surface_runoff + total_infiltration + + # Calculate evaporation 'VERDUNST' by subtracting 'R', the sum of + # runoff and infiltration from precipitation of entire year, + # multiplied by precipitation correction factor + total_evaporation <- climate[["prec_yr"]] - total_runoff + + # Provide total area for calculation of "flows" + total_area <- fetch_data("total_area") + + # Calculate volume 'rowvol' from runoff (qcm/s) + surface_runoff_flow <- yearly_height_to_volume_flow( + total_surface_runoff, total_area + ) + + # Calculate volume 'rivol' from infiltration rate (qcm/s) + infiltration_flow <- yearly_height_to_volume_flow( + total_infiltration, total_area + ) + + # Calculate volume of "system losses" 'rvol' due to surface runoff and + # infiltration + total_runoff_flow <- surface_runoff_flow + infiltration_flow + + # Provide mapping between local variable names and ABIMO-output columns + name_mapping <- list( + code = "CODE", + total_runoff = "R", + total_surface_runoff = "ROW", + total_infiltration = "RI", + total_runoff_flow = "RVOL", + surface_runoff_flow = "ROWVOL", + infiltration_flow = "RIVOL", + total_area = "FLAECHE", + total_evaporation = "VERDUNSTUN" + ) + + # Compose result data frame. Use mget() to get the result vectors from the + # local environment and put them into the data frame + result_data_raw <- cbind( + fetch_data("code", drop = FALSE), + mget(names(name_mapping)[-1L]) + ) + + 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") { + + 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'.") + } + + # Round all columns to three digits (skip first column: "code") + result_data[-1L] <- lapply(result_data[-1L], round, 3L) + + result_data <- restore_geo_column_if_required( + result_data, geometry = geometry + ) + + if (isFALSE(control("intermediates"))) { + return(result_data) + } + + # Return intermediate results as attributes + structure( + result_data, + data = data, + intermediates = list( + climate = climate, + soil_properties = soil_properties, + evaporation_sealed = evaporation_sealed, + evaporation_unsealed = evaporation_unsealed, + roof = list( + evaporation_roof = evaporation_sealed[["roof"]], + runoff_roof = runoff_roof, + runoff_roof_actual = runoff_roof_actual, + infiltration_roof_actual = infiltration_roof_actual + ), + surface = list( + evaporation_sealed = evaporation_sealed[, -1L], + runoff_sealed = runoff_sealed, + runoff_sealed_actual = runoff_sealed_actual, + infiltration_sealed_actual = infiltration_sealed_actual + ), + runoff = cbind( + unsealedSurface = runoff_unsealed, + unsealedRoads = runoff_sealed[, ncol(runoff_sealed)], + sealed = runoff_sealed + ), + fraction_unsealed = fraction_unsealed + ) + ) +} + +# handle_missing_columns ------------------------------------------------------- +handle_missing_columns <- function(data) +{ + road_specific_columns <- c( + "road_frac", "pvd_r", "swg_pvd_r", + "srf1_pvd_r", "srf2_pvd_r", "srf3_pvd_r", "srf4_pvd_r" + ) + + missing_road_columns <- setdiff(road_specific_columns, names(data)) + + if (length(missing_road_columns)) { + for (column in missing_road_columns) { + data[[column]] <- 0 + } + } + + if (! "main_frac" %in% names(data)) { + data$main_frac <- 1 + } + + data +} + +# get_climate: provides climate relevant input data ---------------------------- +get_climate <- function(input) +{ + climate <- select_columns(input, c("prec_yr", "prec_s", "epot_yr", "epot_s")) + + climate[["x_ratio"]] <- climate[["prec_yr"]] / climate[["epot_yr"]] + + climate +} + +# yearly_height_to_volume_flow ------------------------------------------------- + +#' Convert Yearly Height (mm) to Volume Flow (unit?) +#' +#' @param height height in mm +#' @param area area in square metres +#' @keywords internal +yearly_height_to_volume_flow <- function(height, area) +{ + height * 3.171 * area / 100000.0 +} + +#' Define List of "Controls" +#' +#' Define a list of settings that control how the main function +#' \code{\link{run_rabimo}} should behave. +#' +#' @param check logical indicating whether the check functions are executed. +#' Default: \code{TRUE}. +#' @param use_abimo_bagrov_solver logical indicating whether or not to use the +#' (fast!) algorithm implemented in Abimo to solve the Bagrov equations. +#' Default: \code{TRUE}. +#' @param reproduce_abimo_error logical indicating whether or not to reproduce +#' the error that is contained in Abimo (missing area fraction factor). +#' Default: \code{FALSE}. +#' @param output_format one of "abimo" (upper case columns: CODE, R, ROW, RI, +#' RVOL, ROWVOL, RIVOL, FLAECHE, VERDUNSTUN), "rabimo" (lower case columns: +#' code, surface_runoff, infiltration, evaporation). +#' @param intermediates logical indicating whether the intermediate results are +#' returned as attributes. Default: \code{FALSE}. +#' @returns list with the arguments of this function as list elements +#' @export +#' @examples +#' inputs <- kwb.rabimo::rabimo_inputs_2020 +#' test_data <- inputs$data[sample(seq_len(nrow(inputs$data)), size = 1000L), ] +#' controls_default <- define_controls() +#' controls_no_check <- define_controls(check = FALSE) +#' controls_no_solver <- define_controls(use_abimo_bagrov_solver = FALSE) +#' system.time(result_default <- kwb.rabimo::run_rabimo( +#' test_data, inputs$config, controls_default +#' )) +#' system.time(result_no_check <- kwb.rabimo::run_rabimo( +#' test_data, inputs$config, controls_no_check +#' )) +#' identical(result_default, result_no_check) +#' \dontrun{ +#' system.time(result_no_solver <- kwb.rabimo::run_rabimo( +#' test_data, inputs$config, controls_no_solver +#' )) +#' } +define_controls <- function( + check = TRUE, + use_abimo_bagrov_solver = TRUE, + reproduce_abimo_error = FALSE, + output_format = "rabimo", + intermediates = FALSE +) +{ + list( + check = check, + use_abimo_bagrov_solver = use_abimo_bagrov_solver, + reproduce_abimo_error = reproduce_abimo_error, + output_format = output_format, + intermediates = intermediates + ) +} + +#' Crop a box out of a shape +#' +#' @param x sf object +#' @param xoffset x-offset as fraction of original width (0..1) +#' @param yoffset y-offset as fraction of original height (0..1) +#' @param xscale new width as fraction of original width (0..1) +#' @param yscale new height as fraction of original height (0..1) +#' @importFrom sf st_as_sfc st_bbox st_crop +#' @export +crop_box <- function(x, xoffset = 0.45, yoffset = 0.45, xscale = 0.1, yscale = 0.1) +{ + stopifnot(inherits(x, "sf")) + sf::st_crop(x, sf::st_as_sfc(scale_bbox( + bbox = sf::st_bbox(x), xoffset, yoffset, xscale, yscale + ))) +} + +scale_bbox <- function(bbox, xoffset = 0.45, yoffset = 0.45, xscale = 0.1, yscale = 0.1) +{ + xmin <- bbox[["xmin"]] + xmax <- bbox[["xmax"]] + ymin <- bbox[["ymin"]] + ymax <- bbox[["ymax"]] + width <- xmax - xmin + height <- ymax - ymin + xmin_new <- xmin + xoffset * width + ymin_new <- ymin + yoffset * height + sf::st_bbox( + c( + xmin = xmin_new, + xmax = xmin_new + xscale * width, + ymin = ymin_new, + ymax = ymin_new + yscale * height + ), + crs = sf::st_crs(bbox) + ) +} diff --git a/R/stop_on_invalid_data.R b/R/stop_on_invalid_data.R index 0263471c..7de7ce08 100644 --- a/R/stop_on_invalid_data.R +++ b/R/stop_on_invalid_data.R @@ -1,129 +1,144 @@ -# stop_on_invalid_data --------------------------------------------------------- -#' @importFrom rlang .data -stop_on_invalid_data <- function(data) -{ - # Read information on column names and types - column_info <- read_column_info() - - # Helper function to lookup column names matching a property value - columns_with <- function(property, value) { - fetch <- create_accessor(column_info) - fetch("rabimo_berlin")[fetch(property) == value] - } - - # Helper function to check values in columns - check_columns <- function(data, columns, check, msg) { - for (column in columns) { - stopifnot(is.function(check)) - failed <- !check(select_columns(data, column)) - if (any(failed)) { - stop_formatted(msg, column, sum(failed)) - } - } - } - - # Stop if any required column is missing - missing <- setdiff(columns_with("type", "required"), names(data)) - - if (length(missing)) { - info <- dplyr::filter(column_info, .data[["rabimo_berlin"]] %in% missing) - clean_stop("There are missing columns:\n", paste(collapse = "\n", sprintf( - "- %s (%s)", info$rabimo_berlin, info$meaning - ))) - } - - # Stop if a column does not have the expected data type - check_or_convert_data_types( - data = data, - types = get_expected_data_type(columns = names(data)), - convert = FALSE - ) - - # Do not accept any NA - check_columns( - data = data, - columns = names(data) %>% - intersect(columns_with("data_type", "numeric")) %>% - intersect(columns_with("type", "required")), - check = function(x) !is.na(x), - msg = paste( - "Column '%s' must not contain missing values (NA, found %d times).", - "Please give a value (may be 0) in each row." - ) - ) - - # Check precipitation and evapotranspiration for negative values - check_columns( - data = data, - columns = c("prec_yr", "prec_s", "epot_yr", "epot_s"), - check = function(x) x >= 0, - msg = paste( - "There are negative values in column '%s' (%d-times).", - "Please make sure that all values are greater than or equal to zero." - ) - ) - - # Check fractions - check_columns( - data = data, - columns = columns_with("unit", "0..1") %>% - intersect(names(data)), - check = function(x) in_range(x, 0, 1), - msg = paste( - "Not all values in column '%s' are between 0 and 1 as expected", - "(%d failures)." - ) - ) - - check_sum_up_to_1_or_0(data, matching_names(data, pattern_no_roads())) - - if (length(columns <- matching_names(data, pattern_roads()))) { - check_sum_up_to_1_or_0(data, columns) - } -} - -# get_expected_data_type ------------------------------------------------------- -get_expected_data_type <- function(columns = NULL) -{ - columns_to_named_vector <- function(data, key_column, value_column) { - select_columns(data, value_column) %>% - stats::setNames(select_columns(data, key_column)) - } - - type_info <- read_column_info() %>% - columns_to_named_vector( - key_column = "rabimo_berlin", - value_column = "data_type" - ) - - if (is.null(columns)) { - return(type_info) - } - - type_info[intersect(names(type_info), columns)] -} - -# check_sum_up_to_1_or_0 ------------------------------------------------------- -check_sum_up_to_1_or_0 <- function(data, columns, tolerance = 0.005) -{ - # Helper function to check for equality allowing a tolerance - equals <- function(a, b) abs(a - b) <= tolerance - - sums <- rowSums(select_columns(data, columns)) - ok <- equals(sums, 0) | equals(sums, 1) - - if (all(ok)) { - return() - } - - cat("(First) invalid rows:\n") - - select_columns(data, c("code", columns))[!ok, ] %>% - utils::head() %>% - print() - - stop_formatted(string_list(columns), tolerance, x = paste( - "The sum of columns %s is not 1 or 0 in each row as expected", - "(see above). The tolerance was: %f" - )) -} +# stop_on_invalid_data --------------------------------------------------------- +#' @importFrom rlang .data +#' @importFrom kwb.utils stopFormatted +stop_on_invalid_data <- function(data) +{ + # Read information on column names and types + column_info <- read_column_info() + + # Helper function to lookup column names matching a property value + columns_with <- function(property, value) { + fetch <- create_accessor(column_info) + fetch("rabimo_berlin")[fetch(property) == value] + } + + # Helper function to check values in columns + check_columns <- function(data, columns, check, msg) { + for (column in columns) { + stopifnot(is.function(check)) + failed <- !check(select_columns(data, column)) + if (any(failed)) { + kwb.utils::stopFormatted(msg, column, sum(failed)) + } + } + } + + # Stop if any required column is missing + missing <- setdiff(columns_with("type", "required"), names(data)) + + if (length(missing)) { + info <- dplyr::filter(column_info, .data[["rabimo_berlin"]] %in% missing) + clean_stop("There are missing columns:\n", paste(collapse = "\n", sprintf( + "- %s (%s)", info$rabimo_berlin, info$meaning + ))) + } + + # Stop if a column does not have the expected data type + check_or_convert_data_types( + data = data, + types = get_expected_data_type(columns = names(data)), + convert = FALSE + ) + + # Do not accept any NA + check_columns( + data = data, + columns = names(data) %>% + intersect(columns_with("data_type", "numeric")) %>% + intersect(columns_with("type", "required")), + check = function(x) !is.na(x), + msg = paste( + "Column '%s' must not contain missing values (NA, found %d times).", + "Please give a value (may be 0) in each row." + ) + ) + + # Check precipitation and evapotranspiration for negative values + check_columns( + data = data, + columns = c("prec_yr", "prec_s", "epot_yr", "epot_s"), + check = function(x) x >= 0, + msg = paste( + "There are negative values in column '%s' (%d-times).", + "Please make sure that all values are greater than or equal to zero." + ) + ) + + # Check fractions + check_columns( + data = data, + columns = columns_with("unit", "0..1") %>% + intersect(names(data)), + check = function(x) in_range(x, 0, 1), + msg = paste( + "Not all values in column '%s' are between 0 and 1 as expected", + "(%d failures)." + ) + ) + + check_sum_up_to_1_or_0(data, matching_names(data, pattern_no_roads())) + + if (length(columns <- matching_names(data, pattern_roads()))) { + check_sum_up_to_1_or_0(data, columns) + } +} + +# get_expected_data_type ------------------------------------------------------- +get_expected_data_type <- function(columns = NULL) +{ + columns_to_named_vector <- function(data, key_column, value_column) { + select_columns(data, value_column) %>% + stats::setNames(select_columns(data, key_column)) + } + + type_info <- read_column_info() %>% + columns_to_named_vector( + key_column = "rabimo_berlin", + value_column = "data_type" + ) + + if (is.null(columns)) { + return(type_info) + } + + type_info[intersect(names(type_info), columns)] +} + +# check_sum_up_to_1_or_0 ------------------------------------------------------- +#' @importFrom kwb.utils stopFormatted stringList +check_sum_up_to_1_or_0 <- function(data, columns, tolerance = 0.005) +{ + select_columns <- kwb.utils::selectColumns + + # Helper function to check for equality allowing a tolerance + equals <- function(a, b) abs(a - b) <= tolerance + + column_data <- select_columns(data, columns, drop = FALSE) + + # Check for non-numeric columns + is_numeric <- sapply(column_data, is.numeric) + if (any(!is_numeric)) { + clean_stop( + "There are non-numeric columns in check_sum_up_to_1_or_0(): ", + kwb.utils::stringList(columns[!is_numeric]) + ) + } + + sums <- rowSums(column_data) + ok <- equals(sums, 0) | equals(sums, 1) + + if (all(ok)) { + return() + } + + cat("(First) invalid rows:\n") + + select_columns(data, c("code", columns))[!ok, ] %>% + utils::head() %>% + print() + + kwb.utils::stopFormatted(kwb.utils::stringList(columns), tolerance, x = paste( + "The sum of columns %s is not 1 or 0 in each row as expected", + "(see above). The tolerance was: %f" + )) +} diff --git a/R/utils.R b/R/utils.R index 7dbccee6..0f8419da 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1,184 +1,212 @@ -# cat_and_run ------------------------------------------------------------------ -#' @importFrom kwb.utils catAndRun -cat_and_run <- kwb.utils::catAndRun - -# clean_stop ------------------------------------------------------------------- -clean_stop <- function(...) -{ - stop(..., call. = FALSE) -} - -# check_or_convert_data_types -------------------------------------------------- -check_or_convert_data_types <- function( - data, types, convert = FALSE, dbg = TRUE -) -{ - columns <- intersect(names(data), names(types)) - - for (column in columns) { - #column <- columns[2L] - new_type <- types[[column]] - old_type <- class(data[[column]])[1L] - if (old_type != new_type) { - if (!convert) { - stop_formatted( - "Column '%s' (%s) does not have the expected data type (%s).", - column, old_type, new_type - ) - } - cat_and_run( - sprintf("Converting %s from %s to %s", column, old_type, new_type), - dbg = dbg, - data[[column]] <- do.call(paste0("as.", new_type), list(data[[column]])) - ) - } - } - - data -} - -# create_accessor -------------------------------------------------------------- -#' @importFrom kwb.utils createAccessor -create_accessor <- kwb.utils::createAccessor - -# default_if_null -------------------------------------------------------------- -#' @importFrom kwb.utils defaultIfNULL -default_if_null <- kwb.utils::defaultIfNULL - -# expand_to_matrix ------------------------------------------------------------- -expand_to_matrix <- function(x, nrow = NULL, ncol = NULL) -{ - if (is.null(nrow) && is.null(ncol) || !is.null(nrow) && !is.null(ncol)) { - clean_stop( - "Either nrow or ncol must be given but not both at the same time." - ) - } - - if (!is.null(nrow)) { - return(matrix(rep(x, nrow), nrow = nrow, byrow = TRUE)) - } - - if (!is.null(ncol)) { - return(matrix(rep(x, ncol), ncol = ncol, byrow = FALSE)) - } -} - -# get_attribute ---------------------------------------------------------------- -#' @importFrom kwb.utils getAttribute -get_attribute <- kwb.utils::getAttribute - -# helpers_index ---------------------------------------------------------------- -helpers_index <- function(x, values, epsilon = 0.0001, dbg = FALSE) -{ - if (length(x) > 1L) { - return(sapply(x, helpers_index, values, epsilon, dbg)) - } - - stopifnot(length(x) == 1L) - - indices <- which(x <= values + epsilon) - index <- ifelse(length(indices), min(indices), length(values)) - 1L - - print_if(dbg, x) - print_if(dbg, values) - print_if(dbg, indices) - print_if(dbg, index) - - index -} - -# int Calculation::index(float wert, float *feld, int anz) -# { -# int i; -# float eps = 0.0001; -# for (i = 0; i < anz; i++) -# if (wert <= feld[i] + eps) return(i); -# return(anz - 1); -# } - -# in_range --------------------------------------------------------------------- -in_range <- function(x, a, b, tolerance = 0.005) -{ - x + tolerance >= a & x - tolerance <= b -} - -# interpolate ------------------------------------------------------------------ -interpolate <- function(x, y, xout) -{ - yout <- rep(NA_real_, length(xout)) - - nx <- length(x) - - yout[xout <= x[1L]] <- y[1L] - yout[xout >= x[nx]] <- y[nx] - - if (any(is_na <- is.na(yout))) { - yout[is_na] <- sapply(xout[is_na], function(xi) { - i <- which(xi <= x[-1L])[1L] + 1L - (y[i - 1L] + y[i]) / 2 - }) - } - - yout -} - -interpolate_cpp <- function(xi, x, y) -{ - n <- length(x) - stopifnot(n == length(y)) - - if (xi <= x[1L]) { - return(y[1L]) - } - - if (xi >= x[n]) { - return(y[n]) - } - - for (i in seq_len(n)) { - print(i) - if (xi <= x[i + 1L]) { - print(y[i]) - print(y[i+1]) - print ((y[i] + y[i+1]) / 2.0) - return ((y[i] + y[i+1]) / 2.0) - } - } - - return(0.0) -} - -# matching_names --------------------------------------------------------------- -matching_names <- function(data, pattern) -{ - grep(pattern, names(data), value = TRUE) -} - -# print_if --------------------------------------------------------------------- -#' @importFrom kwb.utils printIf -print_if <- kwb.utils::printIf - -# rename_and_select ------------------------------------------------------------ -#' @importFrom kwb.utils renameAndSelect -rename_and_select <- kwb.utils::renameAndSelect - -# rename_columns --------------------------------------------------------------- -#' @importFrom kwb.utils renameColumns -rename_columns <- kwb.utils::renameColumns - -# select_columns --------------------------------------------------------------- -#' @importFrom kwb.utils selectColumns -select_columns <- kwb.utils::selectColumns - -# select_elements -------------------------------------------------------------- -#' @importFrom kwb.utils selectElements -select_elements <- kwb.utils::selectElements - -# stop_formatted --------------------------------------------------------------- -#' @importFrom kwb.utils stopFormatted -stop_formatted <- kwb.utils::stopFormatted - -# string_list ------------------------------------------------------------------ -#' @importFrom kwb.utils stringList -string_list <- kwb.utils::stringList - +# cat_and_run ------------------------------------------------------------------ +#' @importFrom kwb.utils catAndRun +cat_and_run <- kwb.utils::catAndRun + +# clean_stop ------------------------------------------------------------------- +clean_stop <- function(...) +{ + stop(..., call. = FALSE) +} + +# check_or_convert_data_types -------------------------------------------------- +#' @importFrom kwb.utils stopFormatted +check_or_convert_data_types <- function( + data, types, convert = FALSE, dbg = TRUE +) +{ + columns <- intersect(names(data), names(types)) + + for (column in columns) { + #column <- columns[2L] + new_type <- types[[column]] + old_type <- class(data[[column]])[1L] + if (old_type != new_type) { + if (!convert) { + kwb.utils::stopFormatted( + "Column '%s' (%s) does not have the expected data type (%s).", + column, old_type, new_type + ) + } + cat_and_run( + sprintf("Converting %s from %s to %s", column, old_type, new_type), + dbg = dbg, + data[[column]] <- do.call(paste0("as.", new_type), list(data[[column]])) + ) + } + } + + data +} + +# create_accessor -------------------------------------------------------------- +#' @importFrom kwb.utils createAccessor +create_accessor <- kwb.utils::createAccessor + +# default_if_null -------------------------------------------------------------- +#' @importFrom kwb.utils defaultIfNULL +default_if_null <- kwb.utils::defaultIfNULL + +# expand_to_matrix ------------------------------------------------------------- +expand_to_matrix <- function(x, nrow = NULL, ncol = NULL) +{ + if (is.null(nrow) && is.null(ncol) || !is.null(nrow) && !is.null(ncol)) { + clean_stop( + "Either nrow or ncol must be given but not both at the same time." + ) + } + + if (!is.null(nrow)) { + return(matrix(rep(x, nrow), nrow = nrow, byrow = TRUE)) + } + + if (!is.null(ncol)) { + return(matrix(rep(x, ncol), ncol = ncol, byrow = FALSE)) + } +} + +# get_attribute ---------------------------------------------------------------- +#' @importFrom kwb.utils getAttribute +get_attribute <- kwb.utils::getAttribute + +# helpers_index ---------------------------------------------------------------- +helpers_index <- function(x, values, epsilon = 0.0001, dbg = FALSE) +{ + if (length(x) > 1L) { + return(sapply(x, helpers_index, values, epsilon, dbg)) + } + + stopifnot(length(x) == 1L) + + indices <- which(x <= values + epsilon) + index <- ifelse(length(indices), min(indices), length(values)) - 1L + + print_if(dbg, x) + print_if(dbg, values) + print_if(dbg, indices) + print_if(dbg, index) + + index +} + +# int Calculation::index(float wert, float *feld, int anz) +# { +# int i; +# float eps = 0.0001; +# for (i = 0; i < anz; i++) +# if (wert <= feld[i] + eps) return(i); +# return(anz - 1); +# } + +# in_range --------------------------------------------------------------------- +in_range <- function(x, a, b, tolerance = 0.005) +{ + x + tolerance >= a & x - tolerance <= b +} + +# interpolate ------------------------------------------------------------------ +interpolate <- function(x, y, xout) +{ + yout <- rep(NA_real_, length(xout)) + + nx <- length(x) + + yout[xout <= x[1L]] <- y[1L] + yout[xout >= x[nx]] <- y[nx] + + if (any(is_na <- is.na(yout))) { + yout[is_na] <- sapply(xout[is_na], function(xi) { + i <- which(xi <= x[-1L])[1L] + 1L + (y[i - 1L] + y[i]) / 2 + }) + } + + yout +} + +interpolate_cpp <- function(xi, x, y) +{ + n <- length(x) + stopifnot(n == length(y)) + + if (xi <= x[1L]) { + return(y[1L]) + } + + if (xi >= x[n]) { + return(y[n]) + } + + for (i in seq_len(n)) { + print(i) + if (xi <= x[i + 1L]) { + print(y[i]) + print(y[i+1]) + print ((y[i] + y[i+1]) / 2.0) + return ((y[i] + y[i+1]) / 2.0) + } + } + + return(0.0) +} + +# matching_names --------------------------------------------------------------- +matching_names <- function(data, pattern) +{ + grep(pattern, names(data), value = TRUE) +} + +# print_if --------------------------------------------------------------------- +#' @importFrom kwb.utils printIf +print_if <- kwb.utils::printIf + +# remove_geo_column_if_required ------------------------------------------------ +remove_geo_column_if_required <- function(data) +{ + if (inherits(data, "sf")) { + if (is.null(sf_column <- attr(data, "sf_column"))) { + clean_stop("Missing attribute 'sf_column' in data.") + } + # save the geometry in attribute "geometry" and the original name of the + # geometry column in its attribute "sf_column" + structure( + sf::st_drop_geometry(data), + geometry = structure(sf::st_sfc(data[[sf_column]]), sf_column = sf_column) + ) + } else { + data + } +} + +# rename_and_select ------------------------------------------------------------ +#' @importFrom kwb.utils renameAndSelect +rename_and_select <- kwb.utils::renameAndSelect + +# rename_columns --------------------------------------------------------------- +#' @importFrom kwb.utils renameColumns +rename_columns <- kwb.utils::renameColumns + +# restore_geo_column_if_required ----------------------------------------------- +#' @importFrom sf st_as_sf +restore_geo_column_if_required <- function( + data, geometry = attr(data, "geometry") +) +{ + # Is there something stored in attribute "geometry"? + if (is.null(geometry)) { + data + } else { + # remove attribute "geometry" + attr(data, "geometry") <- NULL + # use original name for geometry column and remove attribute "sf_column" + data[[attr(geometry, "sf_column")]] <- structure(geometry, sf_column = NULL) + sf::st_as_sf(data) + } +} + +# select_columns --------------------------------------------------------------- +#' @importFrom kwb.utils selectColumns +select_columns <- kwb.utils::selectColumns + +# select_elements -------------------------------------------------------------- +#' @importFrom kwb.utils selectElements +select_elements <- kwb.utils::selectElements diff --git a/README.md b/README.md index ebb4b27a..b68cb666 100644 --- a/README.md +++ b/README.md @@ -7,9 +7,14 @@ # kwb.rabimo -The code in this package has been transferred from the C++ -code of ABIMO 3.3: Water Balance Model for Urban Areas -(https://github.com/KWB-R/abimo/). +R-implementation of a simple water balance model for urban areas, + +- based on "Wasserhaushaltsmodell Berlin ABIMO 3.2" (see Documentation below) and +- further developed by [KWB](https://kompetenz-wasser.de) + within [BMBF](https://www.bmbf.de/EN/Home/home_node.html)-funded + research project [AMAREX](https://amarex-projekt.de/en). + +For our Tutorial, click [here](https://kwb-r.github.io/kwb.rabimo/dev/articles/tutorial.html) (see also Documentation below). ## Installation @@ -18,21 +23,15 @@ code of ABIMO 3.3: Water Balance Model for Urban Areas install.packages("remotes", repos = "https://cloud.r-project.org") # Install package "kwb.rabimo" (latest "release") from GitHub -remotes::install_github("KWB-R/kwb.rabimo") - -# Install package "kwb.rabimo" (development version) from GitHub -remotes::install_github("KWB-R/kwb.rabimo@dev") +remotes::install_github("KWB-R/kwb.rabimo", build_vignettes = TRUE) ``` ## Basic Usage ### Provide input data and configuration -Compared to the original C++ version of Abimo we have modified the structures -of input data, output data and configuration. - -For the German city of Berlin, we provide data in the new structures in the -package: +For Berlin, the capital of Germany, we provide input data and model parameters +in the package: ```r # Load Berlin data in the original Abimo format @@ -74,6 +73,12 @@ kwb.rabimo::calculate_delta_w( ## Documentation -Release: [https://kwb-r.github.io/kwb.rabimo](https://kwb-r.github.io/kwb.rabimo) +### R-package kwb.rabimo + +- Package Home: https://kwb-r.github.io/kwb.rabimo +- Tutorial: https://kwb-r.github.io/kwb.rabimo/dev/articles/tutorial.html + +### Original software "Wasserhaushaltsmodell Berlin ABIMO 3.2" -Development: [https://kwb-r.github.io/kwb.rabimo/dev](https://kwb-r.github.io/kwb.rabimo/dev) +- Source code (C++): https://github.com/umweltatlas/abimo, +- User manual (in German): https://www.berlin.de/umweltatlas/_assets/literatur/goedecke_et_al_abimo2019_doku.pdf diff --git a/_pkgdown.yml b/_pkgdown.yml index 822fa52d..81a47a30 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -15,6 +15,7 @@ template: primary: '#007aff' border-radius: 0.5rem btn-border-radius: 0.25rem + math-rendering: mathjax development: mode: auto diff --git a/index.md b/index.md deleted file mode 100644 index 1a7daf0b..00000000 --- a/index.md +++ /dev/null @@ -1,30 +0,0 @@ -[![R-CMD-check](https://github.com/KWB-R/kwb.rabimo/workflows/R-CMD-check/badge.svg)](https://github.com/KWB-R/kwb.rabimo/actions?query=workflow%3AR-CMD-check) -[![pkgdown](https://github.com/KWB-R/kwb.rabimo/workflows/pkgdown/badge.svg)](https://github.com/KWB-R/kwb.rabimo/actions?query=workflow%3Apkgdown) -[![codecov](https://codecov.io/github/KWB-R/kwb.rabimo/branch/main/graphs/badge.svg)](https://codecov.io/github/KWB-R/kwb.rabimo) -[![Project Status](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://www.tidyverse.org/lifecycle/#experimental) -[![CRAN_Status_Badge](https://www.r-pkg.org/badges/version/kwb.rabimo)]() -[![R-Universe_Status_Badge](https://kwb-r.r-universe.dev/badges/kwb.rabimo)](https://kwb-r.r-universe.dev/) - -The code in this package has been transferred from the C++ -code of ABIMO 3.3: Water Balance Model for Urban Areas -(https://github.com/KWB-R/abimo/). - -## Installation - -For details on how to install KWB-R packages checkout our [installation tutorial](https://kwb-r.github.io/kwb.pkgbuild/articles/install.html). - -```r -### Optionally: specify GitHub Personal Access Token (GITHUB_PAT) -### See here why this might be important for you: -### https://kwb-r.github.io/kwb.pkgbuild/articles/install.html#set-your-github_pat - -# Sys.setenv(GITHUB_PAT = "mysecret_access_token") - -# Install package "remotes" from CRAN -if (! require("remotes")) { - install.packages("remotes", repos = "https://cloud.r-project.org") -} - -# Install KWB package 'kwb.rabimo' from GitHub -remotes::install_github("KWB-R/kwb.rabimo") -``` diff --git a/man/crop_box.Rd b/man/crop_box.Rd new file mode 100644 index 00000000..7baedf48 --- /dev/null +++ b/man/crop_box.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_rabimo.R +\name{crop_box} +\alias{crop_box} +\title{Crop a box out of a shape} +\usage{ +crop_box(x, xoffset = 0.45, yoffset = 0.45, xscale = 0.1, yscale = 0.1) +} +\arguments{ +\item{x}{sf object} + +\item{xoffset}{x-offset as fraction of original width (0..1)} + +\item{yoffset}{y-offset as fraction of original height (0..1)} + +\item{xscale}{new width as fraction of original width (0..1)} + +\item{yscale}{new height as fraction of original height (0..1)} +} +\description{ +Crop a box out of a shape +} diff --git a/man/generate_rabimo_area.Rd b/man/generate_rabimo_area.Rd index 8ed72d5b..c378d581 100644 --- a/man/generate_rabimo_area.Rd +++ b/man/generate_rabimo_area.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/generate_rabimo_area.R \name{generate_rabimo_area} \alias{generate_rabimo_area} -\title{Generate an area in R-ABIMO format with default values} +\title{Generate an area in R-Abimo format with default values} \usage{ generate_rabimo_area(code, ..., column_info = read_column_info()) } diff --git a/man/rabimo_inputs_2020.Rd b/man/rabimo_inputs_2020.Rd index 516bd63e..1c5e22c0 100644 --- a/man/rabimo_inputs_2020.Rd +++ b/man/rabimo_inputs_2020.Rd @@ -23,7 +23,7 @@ A data.frame with 58531 observations of 33 variables: \item{\code{prec_s}}{Long-term average of annual precipitation within summer months (May to October) in mm (integer)} \item{\code{epot_yr}}{Long-term average of annual potential evapotranspiration in mm (integer)} \item{\code{epot_s}}{Long-term average of annual potential evapotranspiration within summer months (May to October) in mm (integer)} -\item{\code{district}}{Number of Berlin "Bezirk" (district) in which the block area is located (character)} +\item{\code{district}}{Number of Berlin "Bezirk" (district) in which the block area is located (character). This column is Berlin-specific and optional, i.e. not required by the model.} \item{\code{total_area}}{Total block area in square metres (numeric)} \item{\code{main_frac}}{Fraction of the total area that is NOT considered as "road" area (numeric value between 0.0 and 1.0). This value should be 0.0 if roads are modelled separately, i.e. as block areas on their own.} \item{\code{roof}}{Fraction of the total area that is considered as "roof" area (numeric value between 0.0 and 1.0)} @@ -53,7 +53,7 @@ A data.frame with 58531 observations of 33 variables: \item{\code{block_type}}{Block type identifier of the form "usage-type-id_block-type-id_usage-type-description_block-type-description" (character)} } -Note 1: The sum of surface class fractions \code{srf1_pvd}, \code{srf1_pvd}, \code{srf1_pvd}, \code{srf1_pvd}, \code{srf1_pvd} should be 1.0 within each block area. +Note 1: The sum of surface class fractions \code{srf1_pvd}, \code{srf2_pvd}, \code{srf3_pvd}, \code{srf4_pvd}, \code{srf5_pvd} should be 1.0 within each block area. Note 2: The fields with suffix "_r" are all zero because rows are modelled as their own blocks. In an earlier version of the dataset, roads were modelled as parts of the block area. } diff --git a/man/rabimo_inputs_2025.Rd b/man/rabimo_inputs_2025.Rd index b4ed0b5b..005c5f43 100644 --- a/man/rabimo_inputs_2025.Rd +++ b/man/rabimo_inputs_2025.Rd @@ -9,21 +9,21 @@ A list containing two elements: \describe{ -\item{data}{a data frame with the input data in r-abimo format ...(number of vars)} -\item{config}{a list object with configuration data} +\item{data}{a data frame with the input data in R-Abimo format (see below)} +\item{config}{a list object with configuration data (see below)} } } \subsection{\code{rabimo_inputs_2025$data}}{ -A data.frame with 58531 observations of 25 variables: +A data.frame with 58531 observations of 26 variables: \describe{ \item{\code{code}}{Unique block area identifier (character)} \item{\code{prec_yr}}{Long-term average of annual precipitation in mm (integer)} \item{\code{prec_s}}{Long-term average of annual precipitation within summer months (May to October) in mm (integer)} \item{\code{epot_yr}}{Long-term average of annual potential evapotranspiration in mm (integer)} \item{\code{epot_s}}{Long-term average of annual potential evapotranspiration within summer months (May to October) in mm (integer)} -\item{\code{district}}{Number of Berlin "Bezirk" (district) in which the block area is located (character)} +\item{\code{district}}{Number of Berlin "Bezirk" (district) in which the block area is located (character). This column is Berlin-specific and optional, i.e. not required by the model.} \item{\code{total_area}}{Total block area in square metres (numeric)} \item{\code{roof}}{Fraction of the total area that is considered as "roof" area (numeric value between 0.0 and 1.0)} \item{\code{green_roof}}{Fraction of the roof area that belongs to green roofs (numeric value between 0.0 and 1.0). A value of 1.0 means that all roofs in the block area are green roofs.} @@ -35,24 +35,42 @@ A data.frame with 58531 observations of 25 variables: \item{\code{srf3_pvd}}{Fraction of the paved area that belongs to surface class 3 (numeric value between 0.0 and 1.0, see note 1 below)} \item{\code{srf4_pvd}}{Fraction of the paved area that belongs to surface class 4 (numeric value between 0.0 and 1.0, see note 1 below)} \item{\code{srf5_pvd}}{Fraction of the paved area that belongs to surface class 5 (numeric value between 0.0 and 1.0, see note 1 below)} -\item{\code{to_swale}}{Fraction of sealed area (roof area + paved area) that is connected to an infiltration swale (numeric)} +\item{\code{to_swale}}{Fraction of sealed area (roof area + paved area) that is connected to an infiltration swale (numeric value between 0.0 and 1.0)} \item{\code{gw_dist}}{Distance between groundwater table and surface in metres (numeric)} \item{\code{ufc30}}{field capacity in 30 cm depth (numeric)} \item{\code{ufc150}}{field capacity in 150 cm depth (numeric)} \item{\code{land_type}}{land type, one of \code{forested}, \code{horticultural}, \code{urban}, \code{vegetationless}, \code{waterbody} (character)} -\item{\code{veg_class}}{vegetation class index (numeric), derived from an analysis tree volumes} +\item{\code{veg_class}}{vegetation class index (numeric), derived from an analysis of tree volumes} \item{\code{irrigation}}{irrigation in mm per year (integer)} -\item{\code{block_type}}{Block type identifier of the form "usage-type-id_block-type-id_usage-type-description_block-type-description" (character)} +\item{\code{block_type}}{Block type identifier of the form "usage-type-id_block-type-id_usage-type-description_block-type-description" (character). This column is Berlin-specific and optional, i.e. not required by the model.} +\item{\code{Shape}}{List structure containing geometry information on the different block areas. This column is optional. If provided, it will be appended to the model output so that model results can be plotted in the form of maps.} } + +Note 1: The sum of surface class fractions \code{srf1_pvd}, \code{srf2_pvd}, \code{srf3_pvd}, \code{srf4_pvd}, \code{srf5_pvd} should be 1.0 within each block area. } \subsection{\code{rabimo_inputs_2025$config}}{ A list with 3 named elements: \describe{ -\item{runoff_factors}{Runoff factors, vector of numeric with names \code{roof}, \code{surface1}, \code{surface2}, \code{surface3}, \code{surface4}, \code{surface5}} -\item{bagrov_values}{Bagrov values for sealed surfaces, vector of numeric with names \code{roof}, \code{green_roof}, \code{surface1}, \code{surface2}, \code{surface3}, \code{surface4}, \code{surface5}} -\item{swale}{Model parameter(s) related to the 'swale' measure, vector of numeric with currently one value, named \code{swale_evaporation_factor}} +\item{runoff_factors}{Runoff factors for roofs and five different surface +types, given as a vector of numeric with element names \code{roof}, \code{surface1}, +\code{surface2}, \code{surface3}, \code{surface4}, \code{surface5}. A runoff factor determines +the proportion of precipitation that, after subtraction of +evapotranspiration, becomes surface runoff from a paved area. The higher +the factor, the less permeable is the surface.} +\item{bagrov_values}{Bagrov values to calculate evapotranspiration from +paved surfaces, given as a vector of numeric with element names \code{roof}, +\code{green_roof}, \code{surface1}, \code{surface2}, \code{surface3}, \code{surface4}, \code{surface5}. +The higher the Bagrov value, the more evapotranspiration is generated by +the model. For a description of the evapotranspiration model and for a +figure that shows the influence of the Bagrov values (n) on the +evapotranspiration (in German), see \url{https://www.berlin.de/umweltatlas/wasser/wasserhaushalt/2001/methode/}} +\item{swale}{Model parameter(s) related to the 'swale' measure, given as a +vector of numeric with currently one value, named +\code{swale_evaporation_factor}. The swale evaporation factor determines which +fraction of the water going into a swale becomes evapotranspiration (the +rest becomes infiltration).} } } } diff --git a/man/run_rabimo.Rd b/man/run_rabimo.Rd index d917f854..0cc32180 100644 --- a/man/run_rabimo.Rd +++ b/man/run_rabimo.Rd @@ -1,46 +1,50 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/run_rabimo.R -\name{run_rabimo} -\alias{run_rabimo} -\title{Run R-Abimo, the R-implementation of Water Balance Model Abimo} -\usage{ -run_rabimo(data, config, controls = define_controls()) -} -\arguments{ -\item{data}{data frame similar to -\code{\link{rabimo_inputs_2025}$data}} - -\item{config}{configuration object (list) similar to -\code{\link{rabimo_inputs_2025}$config}} - -\item{controls}{list of settings that control how the function should behave. -Use \code{\link{define_controls}} to define such a list. The default is -the list returned by \code{define_controls()}.} -} -\value{ -data frame with columns as returned by Abimo -} -\description{ -Run R-Abimo, the R-implementation of Water Balance Model Abimo -} -\examples{ -# Get input data and config for Berlin (version 2020) -inputs_2020 <- kwb.rabimo::rabimo_inputs_2020 - -# Randomly select 1000 blocks (to reduce runtime) -data <- inputs_2020$data -data <- data[sample(seq_len(nrow(data)), size = 1000L), ] - -# Run R-Abimo -results_2020 <- kwb.rabimo::run_rabimo(data, inputs_2020$config) - -# Get input data and config for Berlin (version 2025) -inputs_2025 <- kwb.rabimo::rabimo_inputs_2025 - -# Randomly select 1000 blocks (to reduce runtime) -data <- inputs_2025$data -data <- data[sample(seq_len(nrow(data)), size = 1000L), ] - -# Run R-Abimo -results_2025 <- kwb.rabimo::run_rabimo(data, inputs_2025$config) -} +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_rabimo.R +\name{run_rabimo} +\alias{run_rabimo} +\title{Run R-Abimo, the R-implementation of Water Balance Model Abimo} +\usage{ +run_rabimo(data, config, controls = define_controls(), silent = FALSE) +} +\arguments{ +\item{data}{data frame similar to +\code{\link{rabimo_inputs_2025}$data}} + +\item{config}{configuration object (list) similar to +\code{\link{rabimo_inputs_2025}$config}} + +\item{controls}{list of settings that control how the function should behave. +Use \code{\link{define_controls}} to define such a list. The default is +the list returned by \code{define_controls()}.} + +\item{silent}{logical indicating whether to suppress console outputs, +the default is \code{FALSE}} +} +\value{ +data frame with columns as returned by Abimo +} +\description{ +Run R-Abimo, the R-implementation of Water Balance Model Abimo +} +\examples{ +# Get input data and config for Berlin (version 2020) +inputs_2020 <- kwb.rabimo::rabimo_inputs_2020 + +# Randomly select 1000 blocks (to reduce runtime) +data <- inputs_2020$data +data <- data[sample(seq_len(nrow(data)), size = 1000L), ] + +# Run R-Abimo +results_2020 <- kwb.rabimo::run_rabimo(data, inputs_2020$config) + +# Get input data and config for Berlin (version 2025) +inputs_2025 <- kwb.rabimo::rabimo_inputs_2025 + +# Crop a box (to reduce runtime) +data <- crop_box(inputs_2025$data) + +# Run R-Abimo +results_2025 <- kwb.rabimo::run_rabimo(data, inputs_2025$config) + +plot(results_2025[, -1L]) +} diff --git a/tests/testthat/test-function-check_sum_up_to_1_or_0.R b/tests/testthat/test-function-check_sum_up_to_1_or_0.R index 8c1b6727..a953aca8 100644 --- a/tests/testthat/test-function-check_sum_up_to_1_or_0.R +++ b/tests/testthat/test-function-check_sum_up_to_1_or_0.R @@ -1,27 +1,33 @@ -test_that("check_sum_up_to_1_or_0() works", { - - f <- kwb.rabimo:::check_sum_up_to_1_or_0 - - expect_error(f()) - - expect_output(expect_error(f( - data.frame( - code = c("btf1", "btf2"), - a = 1:2, - b = 2:3, - c = 3:4 - ), - columns = c("a", "c") - ))) - - expect_null(f( - data.frame( - code = paste0("code_", 1:3), - a = c(0, 1, 0), - b = c(0, 1, 0), - c = c(0, 0, 1) - ), - columns = c("a", "c") - )) - -}) +#library(testthat) +test_that("check_sum_up_to_1_or_0() works", { + + f <- kwb.rabimo:::check_sum_up_to_1_or_0 + + expect_error(f()) + + expect_error( + f(data = data.frame(a = "a"), columns = "a"), + "There are non-numeric columns" + ) + + expect_output(expect_error(f( + data.frame( + code = c("btf1", "btf2"), + a = 1:2, + b = 2:3, + c = 3:4 + ), + columns = c("a", "c") + ))) + + expect_null(f( + data.frame( + code = paste0("code_", 1:3), + a = c(0, 1, 0), + b = c(0, 1, 0), + c = c(0, 0, 1) + ), + columns = c("a", "c") + )) + +}) diff --git a/tests/testthat/test-function-run_rabimo.R b/tests/testthat/test-function-run_rabimo.R index 302cceb1..c3db81d9 100644 --- a/tests/testthat/test-function-run_rabimo.R +++ b/tests/testthat/test-function-run_rabimo.R @@ -68,12 +68,34 @@ test_that("run_rabimo() works", { surface4 = -5, surface5 = -6 ), - swale = list(swale_evaporation_factor = 1) + swale = list( + swale_evaporation_factor = 1 + ) ) - expect_output(result <- f(data, config, controls = define_controls())) + expect_output( + result_1 <- f(data, config, controls = define_controls()) + ) + expect_silent( + result_2 <- f(data, config, controls = define_controls(), silent = TRUE) + ) - expect_s3_class(result, "data.frame") - expect_true(nrow(result) == nrow(data)) + expect_s3_class(result_1, "data.frame") + expect_true(nrow(result_1) == nrow(data)) + expect_identical(result_1, result_2) +}) + +test_that("run_rabimo() keeps the row order", { + inputs <- kwb.rabimo::rabimo_inputs_2020 + data <- inputs$data[sample(nrow(inputs$data), 10L), ] + expect_output(result <- kwb.rabimo::run_rabimo(data, config = inputs$config)) + expect_identical(data$code, result$code) +}) +test_that("run_rabimo() keeps geometry if data inherits from 'sf'", { + inputs <- kwb.rabimo::rabimo_inputs_2025 + data <- inputs$data[sample(nrow(inputs$data), 10L), ] + expect_true("sf" %in% class(data)) + expect_output(result <- kwb.rabimo::run_rabimo(data, config = inputs$config)) + expect_true("sf" %in% class(result)) }) diff --git a/tests/testthat/test-function-run_rabimo_with_measures.R b/tests/testthat/test-function-run_rabimo_with_measures.R index 2498dbfa..32a2d71f 100644 --- a/tests/testthat/test-function-run_rabimo_with_measures.R +++ b/tests/testthat/test-function-run_rabimo_with_measures.R @@ -5,49 +5,44 @@ test_that("run_rabimo_with_measures() works", { expect_error(f()) - data_2020 <- kwb.rabimo::rabimo_inputs_2020$data - data_2025 <- kwb.rabimo::rabimo_inputs_2025$data - - for (data in list(data_2020, data_2025)) { - #data <- data_2025 - + test_me <- function(data) { blocks <- data[sample(seq_len(nrow(data)), 10L), ] - stats <- kwb.rabimo:::get_measure_stats(blocks) - safety_factor <- 0.999 - + measures_max <- list( green_roof = safety_factor * stats$green_roof$max, unpaved = safety_factor * stats$unpaved$max, to_swale = safety_factor * stats$to_swale$max ) - + measures_too_big_1 <- list( green_roof = measures_max$green_roof + 0.01, unpaved = measures_max$unpaved, to_swale = measures_max$to_swale ) - + measures_too_big_2 <- list( green_roof = measures_max$green_roof, unpaved = measures_max$unpaved + 0.01, to_swale = measures_max$to_swale ) - + measures_too_big_3 <- list( green_roof = measures_max$green_roof, unpaved = measures_max$unpaved, to_swale = measures_max$to_swale + 0.01 ) - + expect_output(result <- f(blocks, measures = measures_max)) expect_true(all(result$surface_runoff == 0)) - + expect_error(f(blocks, measures = measures_too_big_1)) expect_error(f(blocks, measures = measures_too_big_2)) expect_error(f(blocks, measures = measures_too_big_3)) - } - + + test_me(data = kwb.rabimo::rabimo_inputs_2020$data) + test_me(data = kwb.rabimo::rabimo_inputs_2025$data) + }) diff --git a/vignettes/.gitignore b/vignettes/.gitignore new file mode 100644 index 00000000..097b2416 --- /dev/null +++ b/vignettes/.gitignore @@ -0,0 +1,2 @@ +*.html +*.R diff --git a/vignettes/tutorial.Rmd b/vignettes/tutorial.Rmd new file mode 100644 index 00000000..02c06f78 --- /dev/null +++ b/vignettes/tutorial.Rmd @@ -0,0 +1,475 @@ +--- +title: "How to Simulate Water Balance with R-Abimo" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{How to Simulate Water Balance with R-Abimo} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +FIGURE_WIDTH <- 5 +``` + +## Background + +This package provides an R-implementation and extension of the simple water +balance model Berlin ABIMO 3.2, which was originally developed by the Federal +Institute of Hydrology ([Bundesanstalt für Gewässerkunde](https://www.bafg.de)) +for rural areas and later adapted for urban areas, namely for Berlin, Germany. +In May 2022, the source code of the model was published on the developer +platform GitHub ([https://github.com/umweltatlas/abimo](https://github.com/umweltatlas/abimo)). + +During the research project [AMAREX](https://amarex-projekt.de/en) (an acronym +for the German translation of "adaptation of stormwater management to extreme +events"), funded by the former German Federal Ministry of Education and Research +([Bundesministerium für Bildung und Forschung -- BMBF](https://www.bmbf.de/EN/Home/home_node.html)), we, the package authors from [Kompetenzzentrum Wasser Berlin gGmbH (KWB)](https://kompetenz-wasser.de/de) +started to work on the original program code, written in the programming language C++ ([https://github.com/KWB-R/abimo](https://github.com/KWB-R/abimo/)). +We then decided to transfer the model to the +[programming language R](https://www.r-project.org/), to rename it to +R-Abimo (see e.g. [here](https://amarex-projekt.de/en/news/abimo-r-abimo)) +and to publish it in the form of this R package [kwb.rabimo](https://github.com/KWB-R/kwb.rabimo). +Compared to the original model, R-Abimo is more generic (i.e. can be more easily +adapted to other cities than Berlin) and it contains some additional features: + +- simulation of stormwater management measures (green roofs, swales), +- "conversion" of urban areas to "natural" areas, +- calculation of delta-W, an indicator for the deviation between the urban (status quo) state and an assumed "natural" state. + +## Prerequisites + +To use the package, you need to have R installed in a version >= 4.3.1. +You can download the current version of R from [here](https://cran.r-project.org/bin/windows/base/). + +## Installation + +In order to install kwb.rabimo directly from our GitHub account [KWB-R](https://github.com/kwb-r/), we recommend to install the R package +remotes first: + +```{r install_remotes, eval = FALSE} +# Install package "remotes" from CRAN +install.packages("remotes") +``` + +You can then install kwb.rabimo in the latest "official" version: + +```{r install_rabimo_release, eval = FALSE} +# Install package "kwb.rabimo" (latest "release") from GitHub +remotes::install_github("KWB-R/kwb.rabimo", build_vignettes = TRUE) +``` + +By setting `build_vignettes = TRUE` you make sure that this tutorial vignette +is installed together with the package. + +## Usage + +### Provide input data and configuration for Berlin + +Compared to the original C++ version of Abimo we have modified the structures +of input data, output data and configuration. +For Berlin, Germany, we provide data in the new structures in the +package: + +```{r load_berlin_data} +# Load Berlin data in the original Abimo format +abimo_inputs <- kwb.rabimo::rabimo_inputs_2025 +``` + +The object `abimo_inputs` is a list with two elements: + +- `abimo_inputs$data` is a data frame containing the actual input data. Each row +represents a block area and each column represents a block area's property. +- `abimo_inputs$config` is a list of parameters that configure the calculation +of runoff and evapotranspiration, respectively, and the infiltration swale model. + +Please refer to the help page of `rabimo_inputs_2025` for further information +on the different input columns and model parameters. + +To open the help page, run + +```{r help_input, eval = FALSE} +?kwb.rabimo::rabimo_inputs_2025 +``` + +In addition to the short descriptions given on the help page you may refer + +- to the [description of the original ABIMO model approaches given by Berlin's + Senate Department for Urban Development, Building and Housing (in German)](https://www.berlin.de/umweltatlas/wasser/wasserhaushalt/2017/methode/) or +- to the [original ABIMO user manual (in German)](https://www.berlin.de/umweltatlas/_assets/literatur/goedecke_et_al_abimo2019_doku.pdf). + +**Note:** We provide also an object `rabimo_inputs_2020`, that refers to an +older version of the Berlin data set. It can be used in almost the same way as +`rabimo_inputs_2025`. However, this old version does not contain geographical +information and we do not cover it within this tutorial. + +### Inspect the input data + +You may inspect the first rows (and only the most relevant columns) of the input +data with + +```{r view_input_data} +head(as.data.frame(abimo_inputs$data)[, 1:24]) +``` + +and you may print the whole configuration object with + +```{r view_input_config} +abimo_inputs$config +``` + +### Visualise the input data + +Since we provide the Berlin dataset together with geographical information +we can plot the spatial distribution of a variable (e.g. the annual +precipitation) in the form of a map: + +```{r, fig.width = FIGURE_WIDTH} +# Provide a subset of the data representing a zoom into the centre of Berlin +berlin_zoom <- kwb.rabimo::crop_box(abimo_inputs$data, + xoffset = 0.35, + yoffset = 0.5, + xscale = 0.07, + yscale = 0.07 +) + +# Plot annual precipitation +plot(berlin_zoom[, "prec_yr"], main = "Annual precipitation in mm") +``` + +### Run R-Abimo for the status quo (urban state) + +To run R-Abimo, call the main function `run_rabimo()` by passing both the input +data and the configuration object to the function: + + +```{r run_rabimo} +# Run R-Abimo +water_balance_urban <- kwb.rabimo::run_rabimo( + data = berlin_zoom, + config = abimo_inputs$config +) + +# Have a look at the first lines of the result data frame +head(sf::st_drop_geometry(water_balance_urban)) +``` +The result data frame contains two columns that are copied from the input data: + +- `code` - unique identifier of each block area, +- `area` - the total block area in square metres in column `area + +as well as three columns containing the model outputs, namely the water balance +components: + +- `runoff` - surface runoff in mm, +- `infiltr` - infiltration in mm, and +- `evapor` - evaporation in mm. + +If the model inputs contained geographical information (as in our case) that +information is restored in the output so that the model results can be plotted +in terms of maps: + +```{r plot_rabimo_result_runoff, fig.width = FIGURE_WIDTH} +# Plot model output "runoff" +plot(water_balance_urban[, "runoff"], main = "Annual runoff in mm") +``` + +```{r plot_rabimo_result_infiltration, fig.width = FIGURE_WIDTH} +# Plot model output "infiltration" +plot(water_balance_urban[, "infiltr"], main = "Annual infiltration in mm") +``` + +```{r plot_rabimo_result_evaporation, fig.width = FIGURE_WIDTH} +# Plot model output "evaporation" +plot(water_balance_urban[, "evapor"], main = "Annual evaporation in mm") +``` + +### Generate your own block areas + +The package provides a function `generate_rabimo_area()` to generate block areas +in the input format that R-Abimo requires. You need to pass at least the unique +identifier(s) (`code`) for the block area(s). If you leave all the other +arguments blank, all the block properties (i.e. columns of the data frame) are +filled with default values. +However, you can override the default values by passing arguments that are named +according to the names of the input columns: + +```{r generate_block_areas} +# Generate artificial block areas, using default values and differing only in +# the fractions of the areas that refer to roofs +codes <- paste0("area-", 1:5) +art_blocks <- kwb.rabimo::generate_rabimo_area( + code = codes, + roof = seq(0, 1, length.out = length(codes)) +) + +# Have a look at the data +art_blocks + +# Run R-Abimo on the block areas +art_water_balance <- kwb.rabimo::run_rabimo(art_blocks, config = abimo_inputs$config) + +# How does the roof area influence the runoff? +plot(art_blocks$roof, art_water_balance$runoff) +``` + +You may use this function to do some sensitivity analysis of the R-Abimo model. + +### Compare the urban state with a natural reference state + +In the research project [AMAREX](https://amarex-projekt.de/de) we propose to +use the deviation of the (urban) water balance from its assumed "natural" state as an +indicator for the vulnerability of an urban area to climate-related effects +(such as heat, flooding, negative impacts on the surface water quality). +As a measure of deviation we introduce $\Delta W$ that is calculated as follows: + +$$\Delta W = \frac{1}{2}\space(|e_{nat}-e_{urb}|+|i_{nat}-i_{urb}|+|r_{nat}-r_{urb}|)\space\frac{100\%}{P}$$ + +with + +- $e_{nat}$ and $e_{urb}$ = evaporation for the natural and urban state, respectively, in mm, +- $i_{nat}$ and $i_{urb}$ = infiltration for the natural and urban state, respectively, in mm, +- $r_{nat}$ and $r_{urb}$ = runoff for the natural and urban state, respectively, in mm, +- $P$ = annual precipitation in mm. + +The higher the value of $\Delta W$ the higher the deviation from the natural +water balance, i.e. the area is less similar to a natural area. +The lower the value of $\Delta W$ the lower is the deviation from the natural +water balance, i.e. the more similar the urban area is to a natural area. + +Our hypothesis is that a lower value of $\Delta W$ indicates a better +preparedness to climate change effects (e.g. because increased evaporation as +it occurs in a natural state scenario leads to a reduction of heat). + +### Run R-Abimo for a natural state scenario + +In the package we provide a function `data_to_natural()` that converts urban +block areas to corresponding "natural" block areas. The function converts +all roof areas and unbuilt paved areas to pervious areas and assumes an overall +vegetation as it can be found in a park (a mixture of meadows, bushes and trees). + +Let's use this function to simulate the water balance for the "natural" equivalent +of our urban area in the Berlin city centre. In all following calls to `run_rabimo()` +we will set `silent = TRUE` so that console outputs are suppressed. + +```{r calculate_natural_state} +# Convert urban state to "natural" state +berlin_zoom_natural <- kwb.rabimo::data_to_natural(berlin_zoom) + +# Calculate the water balance for the "natural" state +water_balance_natural <- kwb.rabimo::run_rabimo( + data = berlin_zoom_natural, + config = abimo_inputs$config, + silent = TRUE +) +``` + +### Calculate and plot "Delta-W" + +Calculate the deviation from the natural state in percent: + +```{r} +delta_w <- kwb.rabimo::calculate_delta_w( + urban = water_balance_urban, + natural = water_balance_natural +) +``` + +The function `claculate_delta_w` returns a data frame with two columns: + +- `code` - the code of the areas, +- `delta_w` - the value of $\Delta W$ in percent. + +```{r} +# Show the first rows of the delta_w data frame (without geometry) +head(sf::st_drop_geometry(delta_w)) +``` + +If the water balance results were provided together with their geographical +information, the spatial distribution of $\Delta W$ values can be visualised in +a map. + +In order to do so, we first define a helper function for plotting Delta-W with a +nice colour palette: + +```{r, fig.width = FIGURE_WIDTH} +# Define function to plot Delta-W +plot_delta_w <- function(data, main) { + palette <- colorRampPalette(c("white", "#a9e0ff", '#7b7bbc', "#350005"))(15) + breaks <- seq(0, 100, length.out = length(palette) + 1) + plot(data[, "delta_w"], main = main, pal = palette, breaks = breaks) +} +``` + +Now let's use the function to plot the Delta-W values in a map: + +```{r, fig.width=FIGURE_WIDTH} +# Plot deviation from natural water balance +plot_delta_w(delta_w, main = "Deviation from natural water balance in %") +``` + +### Simulate rainwater management measures + +R-Abimo currently allows to simulate three rainwater management measures: + +- green roofs, +- unsealing of paved areas, +- infiltration swales. + +These rainwater management measures are shortly demonstrated in the sections +below. In addition to that, the user may also modify the values in the other +columns, if desired, in order to calculate the effect these variables have on +the water balance. + +For example, + +- the composition of different pavement types is controlled by the variables + `srf1_pvd`, `srf2_pvd`, `srf3_pvd`, `srf4_pvd`, and `srf5_pvd` (their sum + represents the total paved area and must be equal to 1.0). Each class is + assigned a pair of model parameters, runoff factor and Bagrov value, + respectively, that are defined in the configuration object (see above). +- changing the value of `roof`, can be useful to simulate the effect of new + buildings or demolitions. +- the parameter `veg_class` can be adjusted to test the impact of vegetation on + the water balance, keeping in mind that a veg_class of 50 corresponds to a + typical urban park with trees, while a value of 75 corresponds to a forested + area. + +#### Simulate green roofs + +Green roofs increase the evaporation and reduce the runoff. In order to be able +to simulate the effect of green roofs in R-Abimo, we have introduced a column +`green_roof` to the input data. The values given in this column refer to the +fractions of the roof areas that refer to green roof areas. + +Let's see how the deviation from the natural state decreases if we introduce more green roofs to our test area in the city centre. We assume that in each block the fraction of roof areas that belong to green roofs is at least 50 percent: + +```{r} +# Make a copy of the original input data +zoom_green_roof <- berlin_zoom + +# Target green roof fraction +green_roof_target <- 0.5 + +# Which blocks need more green roofs? +needs_green_roofs <- zoom_green_roof$green_roof < green_roof_target + +# Set the target green roof fraction +zoom_green_roof$green_roof[needs_green_roofs] <- green_roof_target + +# Simulate the water balance +water_balance_green_roof <- kwb.rabimo::run_rabimo( + data = zoom_green_roof, + config = abimo_inputs$config, + silent = TRUE +) + +# Calculate Delta-W +delta_w_green_roof <- kwb.rabimo::calculate_delta_w( + urban = water_balance_green_roof, + natural = water_balance_natural +) +``` + +```{r, fig.width = FIGURE_WIDTH} +# Plot Delta-W for the status quo +plot_delta_w(delta_w, main = "Delta-W (status quo)") +``` + +```{r, fig.width = FIGURE_WIDTH} +# Plot Delta-W for the green roof scenario +plot_delta_w(delta_w_green_roof, main = "Delta-W (>= 50 % green roofs)") +``` + +The colours become brighter referring to lower Delta-W values i.e. +more "natural" water balances. + +#### Simulate the unsealing of paved areas + +The column `pvd` refers to the fractions of the total areas that are paved. In +order to simulate that paved areas are unsealed we decrease the fractions of +paved areas by a constant factor, let's say we reduce the fractions by 50 +percent: + +```{r} +# Make a copy of the original input data +zoom_unsealed <- berlin_zoom + +# Reduce the fraction of paved areas by 50 percent +zoom_unsealed$pvd <- zoom_unsealed$pvd * 0.5 + +# Simulate the water balance +water_balance_unsealed <- kwb.rabimo::run_rabimo( + data = zoom_unsealed, + config = abimo_inputs$config, + silent = TRUE +) + +# Calculate Delta-W +delta_w_unsealed <- kwb.rabimo::calculate_delta_w( + urban = water_balance_unsealed, + natural = water_balance_natural +) +``` + +```{r, fig.width = FIGURE_WIDTH} +# Plot Delta-W for the status quo +plot_delta_w(delta_w, main = "Delta-W (status quo)") +``` + +```{r, fig.width = FIGURE_WIDTH} +# Plot Delta-W for the green roof scenario +plot_delta_w(delta_w_unsealed, main = "Delta-W (pavement reduced by 50 %)") +``` + +We see an effect similar to that of the green roof scenario: colours become +brighter referring to smaller Delta-W values. + +#### Simulate infiltration swales + +In order to simulate infiltration swales we introduced a column `to_swale` to +the input data. The values in this column indicate the fraction of the total +sealed area (i.e. roof area + unbuilt paved area) that is connected to an +infiltration swale. The infiltration model is quite simple: We assume that the +infiltration swale works perfectly, i.e. a small part of the inflow is converted to evaporation and the remaining part of the inflow is converted to +infiltration. The configuration object contains the `swale_evaporation_factor` +that determines the fraction of the incoming inflow that is converted to evaporation. + +Let's (probably unrealistically) assume that in each block we could connect 50 percent of the sealed area to an infiltration swale. How would that change +Delta-W? + +```{r} +# Make a copy of the original input data +zoom_swale <- berlin_zoom + +# Connect 50 percent of the sealed (= roof + unbuilt paved) areas to swales +zoom_swale$to_swale <- 0.5 + +# Simulate the water balance +water_balance_swale <- kwb.rabimo::run_rabimo( + data = zoom_swale, + config = abimo_inputs$config, + silent = TRUE +) + +# Calculate Delta-W +delta_w_swale <- kwb.rabimo::calculate_delta_w( + urban = water_balance_swale, + natural = water_balance_natural +) +``` + +```{r, fig.width = FIGURE_WIDTH} +# Plot Delta-W for the status quo +plot_delta_w(delta_w, main = "Delta-W (status quo)") +``` + +```{r, fig.width = FIGURE_WIDTH} +# Plot Delta-W for the green roof scenario +plot_delta_w(delta_w_swale, main = "Delta-W (50% of sealed connected to swale)") +``` + +We let it as an exercise to the user to assess which of the three measures has the biggest impact!