pschmidtwalter / LWFBrook90R

Run the LWF-Brook90 hydrological model within R.
https://pschmidtwalter.github.io/LWFBrook90R/
11 stars 7 forks source link

potential adjustment of run_LWFB90 #56

Closed rhabel closed 3 years ago

rhabel commented 3 years ago

Hi,

The chk_errors()-function in L297 (as well as L408ff) of run_LWFB90() stops the code while throwing an error-message https://github.com/pschmidtwalter/LWFBrook90R/blob/e113798157b92cf997fc07ac5667d87cc3ebf287/R/runLWFB90.R#L297

when, say, we're running thousands of points automatically with run_multisite_LWFB90, there will inevitably be points that create errors of some kind. In this case, it would be useful if stopping the whole modelling process through an error would be an optional choice. In a local version of the run_LWFB90-function I wrote a workaround like this:

    simout <- LWFBrook90R:::r_lwfbrook90(
      siteparam = data.frame(simyears[1],
                             as.integer(format(options_b90$startdate, "%j")),
                             param_b90$coords_y, param_b90$snowini, param_b90$gwatini,
                             options_b90$prec_interval),
      climveg = cbind(climate[, c("yr", "mo", "da","globrad","tmax","tmin",
                                  "vappres","windspeed","prec","mesfl")],
                      standprop_daily[, c("densef", "height", "lai", "sai", "age")]),
      precdat = precip,
      param = LWFBrook90R:::param_to_rlwfbrook90(param_b90, options_b90$imodel),
      pdur = param_b90$pdur,
      soil_materials = param_b90$soil_materials,
      soil_nodes = param_b90$soil_nodes[,c("layer","midpoint", "thick", "mat", "psiini", "rootden")],
      output_log = verbose,
      timelimit = timelimit)

    if (simout$error_code != 0L) {
      if (simout$error_code == 1L){outline <- "Simulation terminated abnormally: 'initial matrix psi > 0'"}
      if (simout$error_code == 2L){outline <- "Simulation initialazation failed: 'FWETK failed to determine wetness at KF'"}
      if (simout$error_code == 3L){outline <- "Simulation terminated abnormally: 'inconsistent dates in climate!'"}
      if (simout$error_code == 4L){outline <- "Simulation terminated abnormally: 'inconsistent dates in precipitation input!'"}
      if (simout$error_code == 5L){outline <- "Simulation terminated abnormally: 'wrong precipitation interval input!'"}
      if (simout$error_code == 6L){outline <- "Simulation terminated abnormally: 'negative soil water storage!'"}
      if (simout$error_code == 7L){outline <- "Simulation terminated abnormally: 'water storage exceeds water capacity!'"}

      write(paste0("Error in ID ", unique(soil$id_custom), ": ", outline),
               file = paste0(errorpath, "errors.txt"), append=TRUE)

      simres <- "error"

    }else{
      finishing_time <- Sys.time()
      ...

This way, all points get modelled and one can check which ones caused errors in the error-file. I know, this code-snippet is quite specific for my application (having an \code{id_custom} in the soil-dataframe plus \code{errorpath} in the additional arguments) but potentially the general optional-stopping-functionality might be useful to implement.

pschmidtwalter commented 3 years ago

Hi Raphael,

this is a classic case for an error-routine.... In your custom version, you could use tryCatch() to catch potential errors. This is also what the multi-run functions do in their foreach-loops: When a call of run_LWFB90() fails, the error instead of the simulation results is appended to the list of results. There, you can identify the site by the name of the list entry.

Here an example with run_multi_LWFB90(), but it works the same with run_multisite_LWB90().

library(LWFBrook90R)
data("slb1_meteo")
data("slb1_soil")

# Set up lists containing model control options and model parameters:
parms <- set_paramLWFB90()
opts <- set_optionsLWFB90(startdate = as.Date("2003-06-01"),
                          enddate = as.Date("2003-06-30"))

# Derive soil hydraulic properties from soil physical properties using pedotransfer functions
soil <- cbind(slb1_soil, hydpar_wessolek_tab(slb1_soil$texture))

#set up data.frame with variable parameters
n <- 2
set.seed(2021)
vary_parms <- data.frame(maxlai = runif(n,2,7))

#multi-run
b90.multi <- run_multi_LWFB90(paramvar = vary_parms,
                              param_b90 = parms,
                              options_b90 = opts,
                              climate = slb1_meteo,
                              soil = soil)
str(b90.multi)

# now the same simulation, but with psiini > 0 to provoque an error:
parms$psiini <- 0.1
b90.multi <- run_multi_LWFB90(paramvar = vary_parms,
                              param_b90 = parms,
                              options_b90 = opts,
                              climate = slb1_meteo,
                              soil = soil) 

b90.multi