Sensitivity

library(keys.lid)

### Download "EPA SWMM 5.1.015.7z" from KWB cloud
### "https://cloud.kompetenz-wasser.de/index.php/f/182136"
### and extract to "C:/_UserProg/EPA SWMM 5.1.015/"

### define paths
paths_list <- list(
  root_data = keys.lid::extdata_file(),
  root_swmm = "C:/_UserProg",
  swmm_version = "5.1.013",
  swmm_exe = "<root_swmm>/EPA SWMM <swmm_version>/swmm5.exe",
  lid_models = "<root_data>/sensitivity/models",
  lid_models_zone1 = "<lid_models>/zone1",
  params = "<root_data>/sensitivity/params",
  weather_data = "<root_data>/rawdata/data_weather_sponge_regions",
  sensitivity_results = "<root_data>/sensitivity/results"
)

paths <- kwb.utils::resolve(paths_list)


# go to SWMM GUI and create input file (*.inp)
# input rainfall and temperature data for the desired sponge city climate zone
# are selected there. these files are in paths$weather_data

# read swmm model
swmm_file <- swmmr::read_inp(x = file.path(paths$lid_models_zone1,
                                           'model_greenroof_zone1.inp'))

# read user-defined tables for parameter min and max values
params_min <- read.table(file.path(paths$params, 'params_greenroof_min.csv'),
                         sep = ';',
                         header = TRUE)
params_max <- read.table(file.path(paths$params, 'params_greenroof_max.csv'),
                         sep = ';',
                         header = TRUE)

# number of parameter combinations
l <- 10

# from swmm_file, find which LID parameters are being used (these change based on
# type of LID being modeled)
param_positions <- apply(
  X = params_min[, 2:ncol(params_min)],
  FUN  = function(x) which(!is.na(x)),
  MARGIN = 1)
active_rows <- which(sapply(X = param_positions, FUN = length) > 0)

# initialize output data.frame

# years of analysis
years <- c(2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019)

# make column names
colnamessens <- vector(mode = 'character', length = 0)
for(row in active_rows){
  for(col in param_positions[[row]]){
    colnamessens <- c(colnamessens,
                      paste(params_min$layer[row],
                            colnames(swmm_file$lid_controls)[2 + col],
                            sep = '_'))
  }
}

# create data.frame
sensitivity_results <- data.frame(matrix(
  data = NA,
  ncol = length(colnamessens) + length(years),
  nrow = l,
  dimnames = list(NULL, c(colnamessens,
                          paste('VRR', years, sep='_')))))

# run sensitivity analysis, input time series ideally should cover several years
for(i in seq_len(l)){

  cat('\nmodel run no.', i, '\n')

  # set parameters for simulation -> draw from uniform distribution with ranges
  # given by user in params_min and params_max
  for(row in active_rows){
    for(col in param_positions[[row]]){
      swmm_file$lid_controls[row, col+2] <-
        runif(n = 1,
              min = params_min[[row, col + 1]],
              max = params_max[[row, col + 1]])
    }
  }

  # if the random draw produced field capacity > porosity, set field capacity =
  # 0.99*porosity
  if('SOIL' %in% swmm_file$lid_controls$`Type/Layer`){
    porosity <-
      swmm_file$lid_controls$Par2[swmm_file$lid_controls$`Type/Layer` == 'SOIL']

    field_capacity <-
      swmm_file$lid_controls$Par3[swmm_file$lid_controls$`Type/Layer` == 'SOIL']

    if(field_capacity > porosity){
      swmm_file$lid_controls$Par3[swmm_file$lid_controls$`Type/Layer` == 'SOIL'] <-
        0.99*swmm_file$lid_controls$Par2[swmm_file$lid_controls$`Type/Layer` == 'SOIL']
    }
  }

  # write parameter values in the output data.frame
  for(row in active_rows){
    for(col in param_positions[[row]]){

      currcol <- paste(params_min$layer[row],
                       colnames(swmm_file$lid_controls)[2 + col],
                       sep = '_')

      sensitivity_results[i, currcol] <- swmm_file$lid_controls[row, col+2]
    }
  }

  # save the changed input file, overwriting the original file to avoid writing
  # thousands of copies (one per run)
  swmmr::write_inp(swmm_file, file.path(paths$lid_models_zone1, 'tmp.inp'))

  # run swimm with changed input file
  results <- swmmr::run_swmm(inp = file.path(paths$lid_models_zone1, 'tmp.inp'),
                             exec = paths$swmm_exe)
  swmm_run_has_errors <- !file.exists(file.path(paths$lid_models_zone1, 'tmp.out'))

  if(swmm_run_has_errors) {
      warning(sprintf("SWMM run %d did not complete successfully without errors.",
                      i))

  } else {
  # read out results for itype 3 (= system) and vIndex 4 (= runoff) and 1(= rainfall)
  results_runoff <- swmmr::read_out(results$out, iType = 3, vIndex = 4)
  results_rainfall_rate <- swmmr::read_out(results$out, iType = 3, vIndex = 1)

  # store results in data frame
  results <- data.frame(
    dateTime = zoo::index(results_runoff$system_variable$total_runoff),
    rainfall_rate = zoo::coredata(results_rainfall_rate$system_variable$total_rainfall),
    runoff = zoo::coredata(results_runoff$system_variable$total_runoff),
    years = format(zoo::index(results_runoff$system_variable$total_runoff),
                   format = '%Y'))

  # convert runoff from l/s to mm/s
  flow_units <- swmm_file$options[
    swmm_file$options$Option == 'FLOW_UNITS', 'Value'][[1]]
  if(flow_units == 'LPS'){
    lidarea <- swmm_file$subcatchments$Area
    results$runoff <- results$runoff/(1e4*lidarea)
  }

  # compute rainfall depth [mm] based on rainfall rate [mm/hour] and time
  # step of data [hours]

  # get time interval in hours
  dt <- as.numeric(strsplit(x = swmm_file$raingages$Interval,
                 split = ':')[[1]])
  dt <- dt[1] + dt[2]/60

  # rainfall depth = rainfall rate * time
  results$rainfall_depth <- results$rainfall_rate*dt

  # compute annual VRR (volume rainfall retention) for all analysis years
  # and add it to output data.frame
  years <- unique(results$years)
  vrr <- vector(mode = 'numeric', length = length(years))
  names(vrr) <- years
  for(j in seq_along(years)){
    yearj <- results[results$years == years[j], ]
    runoff_volume <- keys.lid::computeVol(data = yearj,
                                timeColumn = 'dateTime',
                                Qcolumn = 'runoff')
    rainfall_volume <- sum(yearj$rainfall_depth, na.rm = TRUE)

    vrr[j] <- 1 - runoff_volume/rainfall_volume
  }

  sensitivity_results[i,
                      (ncol(sensitivity_results) -
                         length(vrr) + 1):ncol(sensitivity_results)] <- vrr

  rm(results_rainfall_rate)
  fs::file_delete(file.path(paths$lid_models_zone1, 'tmp.out'))
  }

  rm(results)
  fs::file_delete(file.path(paths$lid_models_zone1, 'tmp.inp'))
  fs::file_delete(file.path(paths$lid_models_zone1, 'tmp.rpt'))
}

# write out results
write.table(sensitivity_results,
            file.path(paths$sensitivity_results, 'greenroof_zone1.txt'),
            sep = ';',
            row.names = FALSE,
            quote = FALSE)


# annual VRR vs weather vs. lid parameters
y <- apply(X = sensitivity_results[, (ncol(sensitivity_results) -
                                        12 + 1):ncol(sensitivity_results)],
           MARGIN = 1,
           FUN = mean)

X <- sensitivity_results[, 1:length(colnamessens)]
data_aov <- cbind(y, X)

aov(y ~ ., data = data_aov)
summary(aov(y ~ ., data = data_aov))