R/export_meteo.R

In aemon-j/LakeEnsemblR: Run Ensemble of Lake Models

#' Export LakeEnsemblR standardized input to model specific driver format
#'
#' Export driver files for each model
#'
#' @param config_file filepath; to LakeEnsemblr yaml master config file
#' @param model vector; model to export driving data. Options include
#'    c('GOTM', 'GLM', 'Simstrat', 'FLake')
#' @param folder filepath; to folder which contains the model folders generated by export_config()
#'
#' @examples
#' \dontrun{
#' export_meteo(model = c('GOTM', 'GLM', 'Simstrat', 'FLake'),
#'              meteo_file = 'LakeEnsemblR_meteo_standard.csv')
#' }
#' @importFrom gotmtools get_yaml_value calc_cc input_yaml calc_in_lwr
#' @importFrom glmtools read_nml set_nml write_nml
#' @importFrom zoo na.approx
#' @importFrom lubridate floor_date seconds
#'
#' @export
export_meteo <- function(config_file, model = c("GOTM", "GLM", "Simstrat", "FLake", "MyLake"),
                         folder = "."){

  # check model input
  model <- check_models(model)

  # It's advisable to set timezone to GMT in order to avoid errors when reading time
  original_tz  <-  Sys.getenv("TZ")
  Sys.setenv(TZ = "GMT")

  # Set working directory
  oldwd <- getwd()

  # this way if the function exits for any reason, success or failure, these are reset:
  on.exit({
    setwd(oldwd)
    Sys.setenv(TZ = original_tz)
  })


  yaml  <-  file.path(folder, config_file)

  meteo_file <- get_yaml_value(file = yaml, label = "meteo", key = "file")
  # Check if file exists
  if(!file.exists(meteo_file)){
    stop(meteo_file, " does not exist. Check filepath in ", config_file)
  }

  met_timestep <- get_meteo_time_step(file.path(folder, meteo_file))

  ### Import data
  message("Loading met data...", paste0("[", Sys.time(), "]"))

  met <- read.csv(file.path(folder, meteo_file), stringsAsFactors = FALSE)
  message("Finished loading met data!", paste0("[", Sys.time(), "]"))

  met[, 1] <- as.POSIXct(met[, 1])
  # Check time step
  tstep <- diff(as.numeric(met[, 1]))

  if((mean(tstep) - 86400) / 86400 < -0.05) {
    subdaily <- TRUE
  }else{
    subdaily <- FALSE
  }

  ### Naming conventions standard input
  # test if names are right
  chck_met <- sapply(list(colnames(met)), function(x) x %in% met_var_dic$standard_name)
  if(any(!chck_met)) {
    stop(paste0("Colnames of meteo file are not in standard notation!\n",
                "Colnames: ", paste0(colnames(met)[!chck_met], collapse = ", "),
                ifelse(sum(!chck_met)>1, " are", " is")," wrong.\n",
                "They should be one of: \n", paste0(met_var_dic$standard_name,
                                                    collapse = "\n")))
  }


  # FLake
  #####
  if("FLake" %in% model){

    fla_met <- format_met(met = met, model = "FLake", config_file = config_file)

    # Met output file name
    met_outfile <- "all_meteo_file.dat"
    met_outfpath <- file.path(folder, "FLake", met_outfile)


    # Write meteo file, potentially with the scaling factors in the config_file
    # Using create_scaling_factors in the helpers.R script
    scale_param <- create_scaling_factors(config_file, "FLake", folder)
    scale_met(fla_met, pars = scale_param, model = "FLake", out_file = met_outfpath)


    # Input values to nml
    nml_file <- file.path(folder, get_yaml_value(config_file, "config_files", "FLake"))

    input_nml(nml_file, "SIMULATION_PARAMS", "time_step_number", nrow(fla_met))
    input_nml(nml_file, "METEO", "meteofile", paste0("'", met_outfile, "'"))

    message("FLake: Created file ", file.path(folder, "FLake", met_outfile))

  }

  # GLM
  #####
  if("GLM" %in% model){
    glm_met <- format_met(met = met, model = "GLM", config_file = config_file)

    met_outfile <- file.path("GLM", "meteo_file.csv")

    # Write meteo file, potentially with the scaling factors in the config_file
    # Using create_scaling_factors in the helpers.R script
    scale_param <- create_scaling_factors(config_file, "GLM", folder)
    scale_met(glm_met, pars = scale_param, model = "GLM", out_file = met_outfile)

    # Input to nml file
    nml_path <- file.path(folder, get_yaml_value(config_file, "config_files", "GLM"))
    nml <- glmtools::read_nml(nml_path)

    nml_list <- list("subdaily" = subdaily, "lw_type" = "LW_IN", "meteo_fl" = "meteo_file.csv")
    nml <- glmtools::set_nml(nml, arg_list = nml_list)

    glmtools::write_nml(nml, nml_path)
    message("GLM: Created file ", file.path(folder, "GLM", "meteo_file.csv"))

  }

  ## GOTM
  if("GOTM" %in% model){

    yaml <- file.path(folder, get_yaml_value(config_file, "config_files", "GOTM"))

    met_outfile <- "meteo_file.dat"

    met_outfpath <- file.path(folder, "GOTM", met_outfile)

    got_met <- format_met(met, model = "GOTM", config_file = config_file)

    # Avoid bug where GOTM can crash if last date of met file == last date of simulation
    if(got_met[nrow(got_met), 1] == get_yaml_value(config_file, "time", "stop")){
      last_line <- got_met[nrow(got_met),]
      new_last_date <- format(as.POSIXct(last_line[, 1]) + seconds(met_timestep),
                              "%Y-%m-%d %H:%M:%S")
      last_line[1, 1] <- new_last_date
      got_met <- rbind(got_met, last_line)

      warning("Last date of met file equals last date of simulation. This could cause GOTM to crash ",
              "and therefore one extra time step has been added to the GOTM met file.")
    }

    # Write meteo file, potentially with the scaling factors in the config_file
    # Using create_scaling_factors in the helpers.R script
    scale_param <- create_scaling_factors(config_file, "GOTM", folder)
    scale_met(got_met, pars = scale_param, model = "GOTM", out_file = met_outfpath)

    # Format gotm.yaml file
    ## Set gotm.yaml met config - helper function
    set_met_config_yaml(met = met_outfpath, yaml_file = yaml)

    message("GOTM: Created file ", file.path(folder, "GOTM", met_outfile))

  }

  ## Simstrat
  if("Simstrat" %in% model){

    met_outfile <- "meteo_file.dat"
    par_file <- file.path(folder, get_yaml_value(config_file, "config_files", "Simstrat"))

    met_outfpath <- file.path(folder, "Simstrat", met_outfile)

    sim_met <- format_met(met = met, model = "Simstrat", config_file = config_file)

    # Write meteo file, potentially with the scaling factors in the config_file
    # Using create_scaling_factors in the helpers.R script
    scale_param <- create_scaling_factors(config_file, "Simstrat", folder)
    scale_met(sim_met, pars = scale_param, model = "Simstrat", out_file = met_outfpath)

    ### Write the name of the Simstrat meteo file in the par file
    input_json(file = par_file, label = "Input", key = "Forcing", "\"meteo_file.dat\"")

    message("Simstrat: Created file ", file.path(folder, "Simstrat", met_outfile))
  }

  ## MyLake
  if("MyLake" %in% model){

    met_outfile <- "meteo_file.dat"
    met_outfpath <- file.path(folder, "MyLake", met_outfile)

    # If met_timestep is not 24 hours, MyLake would crash
    # If met_timestep is lower than 24 hours, met is averaged to 24 hours
    if(met_timestep < 86400){
      warning("Meteo time step less than daily; averaging met file for MyLake simulation.")
      met_temp <- aggregate(met,
                            by = list(lubridate::floor_date(met$datetime, unit = "days")),
                            FUN = mean)
      met_temp$datetime <- NULL
      colnames(met_temp)[1] <- "datetime"
    }else if(met_timestep > 86400){
      stop("MyLake cannot be run with meteo forcing time steps larger than 1 day.")
    }else{
      met_temp <- met
    }

    mylake_met <- format_met(met = met_temp, model = "MyLake", config_file = config_file)

    # Write meteo file, potentially with the scaling factors in the config_file
    # Using create_scaling_factors in the helpers.R script
    scale_param <- create_scaling_factors(config_file, "MyLake", folder)
    scale_met(mylake_met, pars = scale_param, model = "MyLake", out_file = met_outfpath)

    message("MyLake: Created file ", file.path(folder, "MyLake", met_outfile))
  }

  # Set the timezone back to the original
  Sys.setenv(TZ = original_tz)

  message("export_meteo complete!")

}