In ryanrbart/rhessysR: Tools for Running RHESSys from R

knitr::opts_chunk$set(echo = F, results = 'hide', eval = F)

RHESSysIOinR Input Functions

This vignette shows example uses of the IOin_X functions, which together construct input data objects used to run one or many RHESSys simulations. For each function see the R help using ?, e.g. ?IOin_input_rhessys for more detailed documentation on each argument.

Input functions:

• IOin_rhessys_input: Basic info for RHESSys command line call
• IOin_tec_std: Dataframe with basic tec events
• IOin_hdr: List of def files to be included in the hdr file
• IOin_def_pars_simple: Def file parameters to be modified
• IOin_cmd_pars: Command line parameter modifiers, being phased out of common use
• IOin_clim: Data required for a climate base station

Setup

Using rhessys/Testing as example input data. If you don't already have develop branch of RHESSys cloned, get that now too.

# to clone and checkout develop branch, a bit slow for some reason
# gert::git_clone(url = "https://github.com/RHESSys/RHESSys",
#                 path = "~/Repos/rhessys-devtest",
#                 branch = "develop")

# Set working dir to the testing folder
knitr::opts_knit$set(root.dir = "~/Repos/rhessys-devtest/Testing/")

Compile RHESSys

If needed, compile RHESSys.

library(RHESSysIOinR)
# This helper function should work for both Mac and WSL (Windows), and avoids dealing with obnoxious WSL paths on the command line
compile_rhessys(location = "~/Repos/rhessys-devtest/")

IOin_rhessys_input

This is the basic information neeeded to construct a RHESSys command line call. For full details see: RHESSys command line options.

input_rhessys = IOin_rhessys_input(
  version = "../rhessys/rhessys7.4",
  tec_file = "tecfiles/tec.test",
  world_file = "worldfiles/w8TC.world",
  world_hdr_prefix = "w8TC",
  flowtable = "flowtables/w8TC.flow",
  start = "1988 9 30 1",
  end = "2000 10 1 1",
  output_folder = "out/",
  output_prefix = "w8TC",
  commandline_options = c("-g")
)

# Along with paths to the other basic inputs, this is the simplest way to assemble a rhessys run, though is limited in its options
run_rhessys_single(
  input_rhessys = input_rhessys,
  hdr_files = "worldfiles/w8TC.hdr",
  output_filter = "tecfiles/testing_filter.yml"
)

IOin_tec_std

This assembles data for a tec file into a dataframe, automatically adding entries for print_daily_on and print_daily_growth_on based on the start date, and conditionally adding output_current_state at the end date. Start and end date can accept R Date formatted objects or RHESSys formatted date strings with hours (space delimited).

input_tec_data = IOin_tec_std(start = "1988 10 1 2",
                              end = "2000 9 30 24",
                              output_state = TRUE)

IOin_hdr

This function creates a list with data for a header file, which will be created via run_rhessys_single. The created hdr file will be named based on world_hdr_prefix from IOin_rhessys_input, and in a folder of the same name. Multiple header files with appended run IDs will be created there to reference different def files modified based on a function like IOin_def_pars_sobol.

input_hdr = IOin_hdr(
  basin = "defs/basin.def",
  hillslope = "defs/hill.def",
  zone = "defs/zone.def",
  soil = "defs/soil_sandyloam.def",
  landuse = "defs/lu_undev.def",
  stratum = "defs/veg_douglasfir.def",
  basestations = "clim/w8_base"
)

Output Filters

Though RHESSys still supports output without output filters, we are moving towards using output filters as a basic component of RHESSys, and so they are included here in the example of the most basic RHESSys run. Output filters require a output filter file in yaml format, containing any number of filters. run_rhessys_single() can write a filter file based on an input R data object (list). The functions below allow for reading of existing output filters, modification, and creation of new output filters.

# read an existing filter file
filter1 = read_output_filter(filter_in = "tecfiles/testing_filter.yml")
# equivalent:
# filter = IOin_output_filters(filter_in = "tecfiles/testing_filter.yml")

# create a new filter R object
filter2 = build_output_filter(
  timestep = "daily",
  output_format = "csv",
  output_path = "../Testing/out",
  output_filename = "basin_daily",
  spatial_level = "basin",
  spatial_ID = as.integer("1"),
  variables = c("patch.total_water_in", "patch.streamflow", "patch.evaporation")
)

# modify an existing filter, either from file or R object, and return R obj
# all options left null (the default) will use the existing value.
filter3 = modify_output_filter(
  filter_in = "tecfiles/testing_filter.yml",
  variables = c("patch.total_water_in", "patch.streamflow", "patch.evaporation")
)

# combine output filters and create object ready to input to run_rhessys_single/multi
filter4 = IOin_output_filters(filter2, filter_in = "tecfiles/testing_filter.yml", file_name = "test_output_filters")

# for just a single filter with named filter file
output_filter = IOin_output_filters(filter2, file_name = "tecfiles/test_output_filter.yml")

Run RHESSys

Together, these basic inputs (basic RHESSys info, tec file, header file, and output filter) can be used to run RHESSys for a single simulation with very simple options. In addition to running using the IOin_X functions to create data objects and later write the needed files, you can use run_rhessys_single by referencing existing tec, hdr, and output filter files as shown above. The data contained by the IOin_rhessys_input() output object is always required to assemble the command line call.

run_rhessys_single(
  input_rhessys = input_rhessys,
  hdr_files = input_hdr,
  tec_data = input_tec_data,
  output_filter = output_filter
)

IOin_def_pars

This function modifies def file parameters. This is done by creating modified def files with the altered parameters, and pointing to that new def file in the header file. This allows for many parameter sets to be generated through many def files and associated header files.

input_def_pars = IOin_def_pars_simple(
  list("defs/soil_sandyloam.def", "m", 0.04269528),
  list("defs/soil_sandyloam.def", "Ksat_0", 1954),
  list("defs/soil_sandyloam.def", "m_z", 0.1423176),
  list("defs/soil_sandyloam.def", "pore_size_index", 0.2337106),
  list("defs/soil_sandyloam.def", "psi_air_entry", 0.8651277),
  list("defs/soil_sandyloam.def", "sat_to_gw_coeff", 0),
  list("defs/hill.def", "gw_loss_coeff", 0)
)

run_rhessys_single(
  input_rhessys = input_rhessys,
  hdr_files = input_hdr,
  tec_data = input_tec_data,
  def_pars = input_def_pars,
  output_filter = output_filter,
  return_cmd = T
)

IOin_clim

This is just to generate a basestation on the fly, but can be useful when running across multiple climate series, and to ensure the paths in the base station file are correct. Future additions will support more complex modification of climate sequences and more.

input_clim = IOin_clim(
  base_station_id = 101,
  x_coordinate = 100.0,
  y_coordinate = 100.0,
  z_coordinate = 975,
  effective_lai = 3.5,
  screen_height = 160,
  daily_prefix = "clim/w8_daily"
)
# if set to existing file name, existing base station will be overwritten.
input_hdr$base_stations = "clim/w8_base_test"

run_rhessys_single(
  input_rhessys = input_rhessys,
  hdr_files = input_hdr,
  tec_data = input_tec_data,
  clim_base = input_clim,
  output_filter = output_filter
)

IOin_cmd_pars

Command line parameters (standard parameters) are multipliers on some commonly varied existing parameter values (found in def files). Because it has become easier to modify the def file values directly, we are encouraging users to move away from this method of parameter variation, but this function is nonetheless how you would input/modify command line parameters.

input_cmd_pars = IOin_cmd_pars(
  m = 0.355794,
  k = 651.390265,
  m_v = 0.355794,
  k_v = 651.390265,
  pa = 1.083102,
  po = 1.193924,
  gw1 = 0.116316,
  gw2 = 0.916922
)

Running Multiple Scenarios

The above input functions all generate RHESSys inputs for a single RHESSys run. With additional options, or different IOin_X functions, multiple RHESSys scenarios can be generated, and run via run_rhessys_multi().

The method for generating those multiple scenarios is flexible, as long as the output data format adheres to what is expected input for run_rhessys_multi(). See the code/documentation/comments on IOin_def_pars_simple() for more info on the data format.

Multiple parameter sets

We previously used IOin_def_pars_simple to change parameters for a single parameter set. IOin_def_pars_simple can also be passed multiple parameter sets, as long as the number of parameters for each set is equal, as in the following example.

input_def_pars_multiple = IOin_def_pars_simple(
  list("defs/soil_sandyloam.def", "m", c(0.04269528, 0.08539056, 0.1280858)),
  list("defs/soil_sandyloam.def", "Ksat_0", c(1954, 977, 488.5)),
  list("defs/soil_sandyloam.def", "m_z", c(0.1423176, 0.2846352, 0.4269528)),
  list("defs/soil_sandyloam.def", "pore_size_index", c(0.2337106, 0.4674212, 0.7011318)),
  list("defs/soil_sandyloam.def", "psi_air_entry", c(0.8651277, 1.730255, 2.595383)),
  list("defs/soil_sandyloam.def", "sat_to_gw_coeff", c(0, 0.2, 0.4)),
  list("defs/hill.def", "gw_loss_coeff", c(0, 0.2, 0.4))
)

run_rhessys_multi(
  input_rhessys = input_rhessys,
  hdr_files = input_hdr,
  tec_data = input_tec_data,
  def_pars = input_def_pars_multiple,
  output_filter = output_filter
)

Monte Carlo

To generate random samples across a range of parameters (i.e. monte-carlo simulation), the above approach can be modified to pass multiple random values to IOin_def_pars_simple. In the example below, runif() generates random samples over a uniform distribution between the minimum and maximum value in the parameter range. Other distributions, such as log-normal, may be substituted.

n=5
input_def_pars_monte_carlo = IOin_def_pars_simple(
  list("defs/soil_sandyloam.def", "m", runif(n, 0.1, 1)),
  list("defs/soil_sandyloam.def", "Ksat_0", runif(n, 0.1, 2000)),
  list("defs/soil_sandyloam.def", "m_z", runif(n, 0.1, 1)),
  list("defs/soil_sandyloam.def", "pore_size_index", runif(n, 0.01, 1)),
  list("defs/soil_sandyloam.def", "psi_air_entry",runif(n, 0.1, 3)),
  list("defs/soil_sandyloam.def", "sat_to_gw_coeff", runif(n, 0, .7)),
  list("defs/hill.def", "gw_loss_coeff", runif(n, 0, .7))
)

run_rhessys_multi(
  input_rhessys = input_rhessys,
  hdr_files = input_hdr,
  tec_data = input_tec_data,
  def_pars = input_def_pars_monte_carlo,
  output_filter = output_filter
)

IOin_def_pars_latin_hypercube

Many more sophisticated approaches for establishing parameter sets (e.g. latin-hypercube, sobol sensitivity) are generated by sampling a multidimensional distribution among all the parameters, rather than sampling the distribution of each individual parameter. For these approaches, IOin_def_pars_simple is too simple and a tailored IOin function will need to be used.

IOin_def_pars_latin_hypercube generates parameters according to [latin hypercube sampling] (https://en.wikipedia.org/wiki/Latin_hypercube_sampling). In this example, the third element of each list represents the parameter ranges to be passed to IOin_def_pars_latin_hypercube, with the three values in the vector representing the total number of parameter sets, the minimum value of parameter sampling range, the and maximum value of parameter sampling range.

n=5
input_def_pars_latin_hypercube = IOin_def_pars_latin_hypercube(
  list("defs/soil_sandyloam.def", "m", c(n, 0.1, 1)),
  list("defs/soil_sandyloam.def", "Ksat_0", c(n, 0.1, 2000)),
  list("defs/soil_sandyloam.def", "m_z", c(n, 0.1, 1)),
  list("defs/soil_sandyloam.def", "pore_size_index", c(n, 0.01, 1)),
  list("defs/soil_sandyloam.def", "psi_air_entry",c(n, 0.1, 3)),
  list("defs/soil_sandyloam.def", "sat_to_gw_coeff", c(n, 0, .7)),
  list("defs/hill.def", "gw_loss_coeff", c(n, 0, .7))
)

run_rhessys_multi(
  input_rhessys = input_rhessys,
  hdr_files = input_hdr,
  tec_data = input_tec_data,
  def_pars = input_def_pars_latin_hypercube,
  output_filter = output_filter
)

IOin_def_pars_sobol

This is still work in progress, but shows how Sobol parameter sampling can be used to generate RHESSys inputs.

# this gives the most flexibility to how users want to generate a distribution
n = 10
pars = list(
  list("defs/veg_douglasfir.def", "epc.waring_pa", runif(n, 0.1, 0.4)),
  list("defs/veg_douglasfir.def", "epc.proj_sla", runif(n, 2, 6)),
  list("defs/soil_sandyloam.def", "pore_size_index", runif(n, 0.2, 0.5)),
  list("defs/soil_sandyloam.def", "psi_air_entry", runif(n, 0.5, 8))
)

input_def_pars_sobol = IOin_def_pars_sobol(pars)

# For speed of the example...
input_def_pars_sobol = lapply(input_def_pars_sobol, FUN = function(x) {x[[3]] = x[[3]][1:5]; return(x)})

run_rhessys_multi(
  input_rhessys = input_rhessys,
  hdr_files = input_hdr,
  tec_data = input_tec_data,
  def_pars = input_def_pars_sobol,
  output_filter = output_filter
)

ryanrbart/rhessysR documentation built on March 29, 2024, 4:30 p.m.