In hydrosolutions/riversCentralAsia: Accompanying R Package to Applied Modeling of Hydrological Systems in Central Asia Course

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", 
  eval = nzchar(Sys.getenv("riversCentralAsia_vignette_eval"))
)
knitr::opts_knit$set(
  rmarkdown.html_vignette.check_title = FALSE)

Preparation of a GIS layer with hydrological response unit

Typically, a river basin is discretized into hydrological response units (HRU) that capture the hydrological properties of a zone within the river basin. Often, HRU are derived from digital elevation models (DEMs) to be able to represent the elevation-dependent timing of snow melt in a basin. This is especially suitable for mountainous basins where gravity flow determines the river runoff. However, in lower lying basins, geology or land use may dominate the hydrology and thus HRU may have to be derived from geological or land use maps.

The function gen_basinElevationBands allows the generation of a shape layer with elevation bands derived from a DEM. Note that you can play with the parameters of the function to adjust the shapes of the elevation bands.

To reproduce this vignette, please make sure you have this package as well as the example data package downloaded (see README{target="_blank"}).

library(tidyverse)
library(lubridate)
library(timetk)
library(riversCentralAsia)

gis_path <- "../../riversCentralAsia_ExampleData/16076_DEM.tif"
crs_project <- "+proj=longlat +datum=WGS84"  # EPSG:4326, latlon WGS84

# Parameter definition for the generation of the elevation bands
band_interval <- 500 # in meters. Note that normally you want to work with band intervals of 100 m to 200 m. To make the model less computationally demanding, we work with a coarser resolution of 500 m for this demo. 
holeSize_km2 <- .1 # cleaning holes smaller than that size
smoothFact <- 2 # level of band smoothing
demAggFact <- 2 # dem aggregation factor (carefully fine-tune this)
## Delineation
hru_shp <- gen_basinElevationBands(
  dem_PathN = gis_path, demAggFact = demAggFact, band_interval = band_interval, 
  holeSize_km2 = holeSize_km2, smoothFact = smoothFact)

# Control output
hru_shp %>% plot()

Now, only minor post-processing (like adding a name attribute) is further required and the HRU layer is ready for import to RS Minerve for use in the automated model creation. The open-source book Modeling of Hydrological Systems in Semi-Arid Central Asia, section on geospatial data preparation{target="_blank"} gives step-by-step instructions for this.

Extraction of climate data on hydrological response units

We recommend the use of daily precipitation and temperature data from the CHELSA data set (made available by WSL)[@karger_climatologies_2017; @karger_high-resolution_2020]. The data is downloaded from the CHELSA server, then extracted to each hydrological response unit and reformated to the RS MINERVE forcing input format using function gen_HRU_Climate_CSV_RSMinerve. Please note that the global CHELSA data set at 1km grid resolution requires several GB of free storage space on your computer.

This process is demonstrated in detail in the book Modelling of Hydrological Systems in Semi-Arid Central Asia, section Climate Data - Downscaling and Bias Correction{target="_blank"}. Below we show an example useage from this book chapter. This example is reproducible with the (larger) example data set from the book but not with the example data set in this package.

# Temperature data processing
hist_obs_tas <- riversCentralAsia::gen_HRU_Climate_CSV_RSMinerve(
  climate_files = <List of CHELSA temperature files>,
  catchmentName = <river_name>,
  temp_or_precip = "Temperature",  # To process precipitation data, type "Precipitation"
  elBands_shp = <hru_shp>,
  startY = <start year>,
  endY = <end year>,
  obs_frequency = <'hour', 'day', or 'month'>,  # Note that CHELSA data is available at daily or monthly frequency, not hourly.
  climate_data_type = <'hist_obs', 'hist_sim', or 'fut_sim'>,
  crs_in_use = <a projection code, for example '+proj=longlat +datum=WGS84'>)

An example using a single input file is given below.

# Assign names and elevations to the HRU
# Here we use demo elevations. Typically you would extract elevations from DEM. 
#   Z <- exactextractr::exact_extract(dem, hru_shp, 'mean')
hru_shp <- hru_shp |> 
  mutate(name = paste0("Atb_", layer), 
         Z = c(1000, 2000, 3000, 4000))  

# Temperature data processing
fut_sim_tas <- riversCentralAsia::gen_HRU_Climate_CSV_RSMinerve(
  climate_files = 
    "../../riversCentralAsia_ExampleData/tas_day_GFDL-ESM4_ssp126_r1i1p1f1_gr1.nc",
  catchmentName = "Atbashy",
  temp_or_precip = "Temperature",
  elBands_shp = hru_shp,  # from the code chunk above
  startY = 2012,
  endY = 2017,
  obs_frequency = "day", 
  climate_data_type = "fut_sim",
  crs_in_use = '+proj=longlat +datum=WGS84')

# This gives you the temperature forcing for each elevation band in a format that can be imported to RS MINERVE. 
fut_sim_tas
# A tibble: 2,200 x 5
#   Station             Atb_1             Atb_2             Atb_3             Atb_4            
#   <chr>               <chr>             <chr>             <chr>             <chr>            
# 1 Station             Atb_1             Atb_2             Atb_3             Atb_4            
# 2 X                   1121699.41259477  1137845.07398682  1152073.09136446  1162907.10435686 
# 3 Y                   4593793.26535944  4597636.57285844  4600746.13179295  4602406.26249851 
# 4 Z                   1000              2000              3000              4000             
# 5 Sensor              T                 T                 T                 T                
# 6 Category            Temperature       Temperature       Temperature       Temperature      
# 7 Unit                C                 C                 C                 C                
# 8 Interpolation       Linear            Linear            Linear            Linear           
# 9 01.01.2012 00:00:00 -14.048846244812  -14.0820846557617 -14.0855207443237 -14.0795421600342
#10 02.01.2012 00:00:00 -13.4489545822144 -13.4161062240601 -13.4127111434937 -13.4186210632324
# ... with 2,190 more rows
# i Use `print(n = ...)` to see more rows

# You now just need to write the table to csv and import it to RS MINERVE. 
# readr::write_csv(fut_sim_tas, <filename>, na = "NA", col_names = FALSE)

This produces a tibble with the daily average temperatures for each elevation band in the input elBands_shp layer.

Bias correction using quantile mapping

Also here, due to the size of the input files, we do not provide a reproducible example in this package but refer the interested reader to the book Modelling of Hydrological Systems in Semi-Arid Central Asia, section Climate Data - Downscaling and Bias Correction{target="_blank"} where reproducible examples of the downscaling of projected climate data are available. The bias correction is done using the R package qmap{target="_blank"} [@gudmundsson_qmap_2016].

The synthetic example below demonstrates how the wrapper for qmap in riversCentral Asia is used to downscale climate projections to each HRU of the basin.

hist_obs <- tibble::tribble(~Date, ~Basin, ~Pr,
                            "1979-01-01", "K_eb1", 0.1,
                            "1979-01-01", "K_eb2", 0.2,
                            "1979-01-01", "K_eb3", 0.3,
                            "1979-01-02", "K_eb1", 0.4,
                            "1979-01-02", "K_eb2", 0.5,
                            "1979-01-02", "K_eb3", 0.6) |>
  dplyr::mutate(Date = as.Date(Date))
hist_sim <- hist_obs |>
  dplyr::filter(Basin == "K_eb1") |>
  dplyr::select(-Basin) |>
  dplyr::mutate(Pr = Pr + 1, Model = "A")
hist_sim <- hist_sim |>
  dplyr::add_row(hist_sim |>
                   dplyr::mutate(Pr = Pr + 2, Model = "B"))
fut_sim <- hist_sim |>
  dplyr::mutate(Scenario = "a") |>
  dplyr::add_row(hist_sim |>
                   dplyr::mutate(Pr = Pr + 1, Scenario = "b"))
fut_sim <- fut_sim |>
  dplyr::mutate(Date = as.Date(Date) + 2) |> 
  dplyr::add_row(fut_sim |>
                   dplyr::mutate(Date = as.Date(Date) + 4))

results <- doQuantileMapping(hist_obs, hist_sim, fut_sim)
mapped_hist_sim <- results[[1]]
mapped_fut_sim <- results[[2]]

mapped_hist_sim <- mapped_hist_sim |> 
  mutate(Model = ifelse(Model == "ERA5-CHELSA", "baseline", Model))

# Visually inspect results
ggplot() + 
  geom_line(data = hist_sim, aes(Date, Pr, colour = Model), linetype = 2) + 
  geom_line(data = mapped_hist_sim, aes(Date, Pr, colour = Model)) + 
  geom_line(data = hist_obs, aes(Date, Pr), colour = "black", linetype = 3) + 
  geom_line(data = fut_sim, aes(Date, Pr, colour = Model), linetype = 2) + 
  geom_line(data = mapped_fut_sim, aes(Date, Pr, colour = Model)) + 
  facet_wrap(c("Basin", "Scenario")) + 
  theme_bw()

References

hydrosolutions/riversCentralAsia documentation built on Feb. 7, 2023, 4:50 p.m.