# Wageningen Environmental Research
# waterwijzerlandbouw@wur.nl
#------------------------------------------------

# file KNMI:
# variable  |  alias  |  unit
# temp         T         [degree C]
# precip       RH        [mm]
# radiation    Q         [W m-2]
# max-temp     TX        [degree C]
# min-temp     TN        [degree C]
# rel-humidity U         [%]
# windspeed    FG        [m s-1]

# set start of program
TIMPRGSTART <- Sys.time()

# --- arguments ---

# set modus
test <- FALSE

# load arguments
if (!test) {
  command_args <- commandArgs(trailingOnly = TRUE)
  if (length(command_args) != 1) stop("wrong usage of script, arguments should be:\n- controlfile")
} else {
  setwd("c:/WaterwijzerLandbouw/WWL_regionaal")
  command_args <- NULL
  command_args[1] <- "control_meteo_projection.inp"
}

# set arguments
ctrl <- command_args[1]

# --- load libraries ---

source("tools/R/libraries_meteo.R")

message("\nProgram R-PROG started...\n")

# --- initial part of procedure ---

# load controlfile
text <- read_lines(file = ctrl, lazy = FALSE, num_threads = 1, progress = FALSE)

# set directories
dir_out <- get_dir(text = text, item = "DIROUT", create = FALSE)

# data from KNMI (https://klimaatscenarios-data.knmi.nl/tijdreeks)
file_station <- get_record(text = text, item = "STATION", split = TRUE)
file_replace <- get_record(text = text, item = "REPLACE", split = TRUE, item_exists = FALSE)

# period
tstart <- get_record(text = text, item = "TSTART", item_exists = FALSE)
tend <- get_record(text = text, item = "TEND", item_exists = FALSE)

# split met-files
split <- get_switch(text = text, item = "SPLIT", item_exists = FALSE, default = FALSE)

# intensity of precipitation (optional)
intensity <- get_switch(text = text, item = "INTENSITY", item_exists = FALSE, default = TRUE)
if (intensity) file_cofwet <- get_record(text = text, item = "COFWET")

# windspeed factor
f_windspeed <- get_record(text = text, item = "FWIND", item_exists = FALSE, default = 1.0)

# humidity factor
f_humidity <- get_record(text = text, item = "FHUM", item_exists = FALSE, default = 1.0)

# add startup period
startup <- get_switch(text = text, item = "STARTUP", item_exists = FALSE, default = TRUE)

# --- main part of procedure ---

# load coefficients to determine precipitation intensity
if (intensity) {
  db_cofwet <- read_csv(file = file_cofwet, col_types = "idd", lazy = FALSE, num_threads = 1, progress = FALSE) |>
    mutate(MM = str_pad(string = month, width = 2, pad = "0")) |>
    select(MM,a,b)
}

# loop over station
for (i_stn in seq_along(file_station)) {
  
  done <- FALSE
  n_rpl <- length(file_replace)
  i_rpl <- 0
  while (!done) {
    
    # set filename knmi
    load_data <- TRUE
    if (i_rpl == 0) {
      file_knmi <- file_station[i_stn]
    } else {
      file_knmi <- file_replace[i_rpl]
      if (file_knmi == file_station[i_stn]) load_data <- FALSE
    }
    
    # process data
    if (load_data) {
      
      # extract info knmi-data
      knmi <- path_ext_remove(path = path_file(path = file_knmi))
      stn_id <- str_split_i(string = knmi, pattern = "_", i = 1)
      station <- str_split_i(string = knmi, pattern = "_", i = 2)
      projection <- str_split_i(string = knmi, pattern = "_", i = 4)
      
      if (i_rpl == 0) {
        message(str_c("\nprocessing station ", station, " (", stn_id, ")"))
        station_info <- get_station_KNMI(station_id = stn_id, type = "precipitation")
        station_ref <- station
        projection_ref <- projection
        replace <- NULL
      } else {
        if (projection_ref != projection) stop("wrong projection selected for replacement")
        replace <- c(replace, station)
      }
      
      # load data
      db_raw <- read_csv(file = file_knmi, comment  = "#", col_types = "i-dddddd-", col_names = c("date", "RH", "Q", "TX", "TN", "U", "FG")) |>
        mutate(
          date = ymd(.data$date),
          YYYY = year(.data$date),
          MM = str_pad(string = as.character(month(.data$date)), width = 2, pad = "0"),
          DD = str_pad(string = day(.data$date), width = 2, pad = "0")
        )
      
      # default period
      if (is.null(tstart)) tstart <- min(db_raw$date)
      if (is.null(tend)) tend <- max(db_raw$date)
      
      # check selected period
      if (tstart < min(db_raw$date) || tend > max(db_raw$date)) stop("selected period 'TSTART' or 'TEND' outside range")
      
      # filter period
      db_raw <- db_raw |>
        filter(.data$date >= tstart, .data$date <= tend)
      
      # set missing values
      db_raw[db_raw < -50] <- NA_real_
      
      # set station
      db_raw <- db_raw |> mutate(STN = str_c("'", stn_id, "'"))
      
      # replace missing values
      if (i_rpl == 0) {
        db <- db_raw
        if (length(db[is.na(db)]) > 0) message(str_c("\nmissing observations found for station ", station, "(", stn_id, ")"))
      } else {
        for (i_col in 1:ncol(db)) {
          if (names(db)[i_col] %in% c("RH", "Q", "TN", "TX", "U", "FG")) {
            if (length(db[is.na(db[, i_col]), i_col]) > 0) {
              db[is.na(db[, i_col]), "STN"] <- str_c("'", stn_id, "'")
              db[is.na(db[, i_col]), i_col] <- db_raw[is.na(db[, i_col]), i_col]
            }
          }
        }
      }
    }
    
    # check missing data
    if (!any(is.na(db)) || i_rpl >= n_rpl) {
      done <- TRUE
    } else {
      i_rpl <- i_rpl + 1
    }
  }
  
  if (length(db[is.na(db)]) == 0 && i_rpl > 0) message("all missing observation(s) replaced")
  
  # convert radiation [W m-2] to [kJ m-2]
  db <- db |>
    mutate(Q = .data$Q * 86.4)
  
  # convert relative humidity UG [%] to actual vapour pressure HUM [kPa]; optional usage of humidity factor
  db <- db |> 
    mutate(HUM = actual_vapour_pressure(TN = .data$TN, TX = .data$TX, UG = .data$U) * f_humidity)
  
  # optional usage of windspeed factor
  db <- db |> 
    mutate(FG = .data$FG * f_windspeed)

  # add missing data
  db <- db |> mutate(EV24 = NA_real_, DR = NA_real_)
  
  # determine precipitation intensity
  if (intensity) {
    db <- left_join(x = db, y = db_cofwet, by = join_by("MM")) |>
      mutate(DR = if_else(.data$RH > 0, pmax(0.001, pmin(1.0, .data$RH / (24*(.data$a*.data$RH^.data$b)))), 0)) |>
      select(-"a", -"b")
  }
  
  # add start-up period
  if (startup) {
    db_startup <- db |>
      filter(.data$YYYY == 1991) |>
      mutate(YYYY = 1990)
    db <- rbind(db_startup, db)
  }
  
  # prepare output
  db <- db |>
    mutate(
      Q = if_else(is.na(.data$Q), NA_character_, formatC(x = .data$Q, format = "f", digits = 1)),
      TN = if_else(is.na(.data$TN), NA_character_, formatC(x = .data$TN, format = "f", digits = 1)),
      TX = if_else(is.na(.data$TX), NA_character_, formatC(x = .data$TX, format = "f", digits = 1)),
      HUM = if_else(is.na(.data$HUM), NA_character_, formatC(x = .data$HUM, format = "f", digits = 6)),
      FG = if_else(is.na(.data$FG), NA_character_, formatC(x = .data$FG, format = "f", digits = 1)),
      RH = if_else(is.na(.data$RH), NA_character_, formatC(x = .data$RH, format = "f", digits = 3)),
      EV24 = if_else(is.na(.data$EV24), NA_character_, formatC(x = .data$EV24, format = "f", digits = 1)),
      DR = if_else(is.na(.data$DR), NA_character_, formatC(x = .data$DR, format = "f", digits = 3))
    )
  
  # create header
  header <- c(
    "*------------------------------------------------------------------------------------",
    "* Source of data    : Royal Netherlands Meterological Institute (https://www.knmi.nl)",
    "* File content      : Meteo from KNMI; input file for SWAP 4.2 (https://swap.wur.nl)",
    "*",
    str_c("* Station           : ", str_to_upper(station_ref), " (", station_info$STATION_ID[1],")"),
    str_c("* Longitude         : ", ifelse(is.na(station_info$LON[1]), "unknown", station_info$LON[1])),
    str_c("* Latitude          : ", ifelse(is.na(station_info$LAT[1]), "unknown", station_info$LAT[1])),
    str_c("* Altitude          : ", ifelse(is.na(station_info$ALT[1]), "unknown", station_info$ALT[1])),
    str_c("* Projection        : ", projection),
    "*",
    "* Rad               : Solar radiation (kJ m-2)",
    "* Tmin              : Minimum air temperature (degree C)",
    "* Tmax              : Maximum air temperature (degree C)",
    str_c("* Hum               : Air humidity (kPa); humidity factor ", f_humidity, " (-)"),
    str_c("* Wind              : Average windspeed; 24h (m s-1); windspeed factor ", f_windspeed, " (-)"),
    "* Rain              : Rainfall (mm)",
    "* ETref             : Reference evapotranspiration; Makkink (mm)",
    str_c("* Wet               : Fraction of time with rainfall (-)", ifelse(!intensity, "", "; intensity based on daily rainfall"))
  )
  
  if (i_rpl > 0) {
    header <- c(header, "*")
    rpl_id <- db |> filter(.data$STN != str_c(",", station_info$STATION_ID[1], "'")) |> pull("STN") |> unique() |> str_replace_all(pattern = "'", replacement = "")
    for (i_rpl in seq_along(rpl_id)) {
      replace_info <- get_station_KNMI(station_id = rpl_id[i_rpl], type = "precipitation")
      header <- c(header, str_c("* Missing values are replaced by station ", replace_info$STATION_NAME[1], " (", replace_info$STATION_ID[1],")"))
    }
  }
  
  if (startup) {
    header <- c(header, "*", "* Start-up period used (duplicated 1991)")
  }
  
  header <- c(header,
              "*",
              "* File generated by : Wageningen Environmental Research",
              str_c("* File generated at : ",str_sub(string = Sys.time(), start = 1, end = 20)),
              "*------------------------------------------------------------------------------------",
              "Station,DD,MM,YYYY,Rad,Tmin,Tmax,Hum,Wind,Rain,ETref,Wet"
  )
  
  # create output directory
  dir_create(dir_out)
  
  # create meteo files
  if (!split) {
    
    # prepare output
    db <- db |> select("STN", "DD", "MM", "YYYY", "Q", "TN", "TX", "HUM", "FG", "RH", "EV24", "DR")
    
    # write data
    file <- str_c(dir_out, "/", str_pad(string = station_info$STATION_ID[1], width = 3, pad = "0"), "_", str_to_lower(string = str_c(str_sub(string = projection, start = 1, end = 4) , "_", str_sub(string = projection, start = 5, end = 6))), "_fw_fh.met")
    writeLines(text = header, con = file)
    write.table(x = db, file = file, na = "-99.9", sep = ",", dec = ".", append = TRUE, quote = FALSE, col.names = FALSE, row.names = FALSE)
    
    if (any(is.na(db |> select(-"DR", -"EV24")))) warning("incomplete data for file: ", file, "\n")
    
  } else {
    
    for (year in unique(db$YYYY)) {
      
      s_db <- db |> filter(YYYY == year)
      
      # prepare output
      s_db <- s_db |> select("STN", "DD", "MM", "YYYY", "Q", "TN", "TX", "HUM", "FG", "RH", "EV24", "DR")
      
      # write data
      file <- str_c(dir_out, "/", str_pad(string = station_info$STATION_ID[1], width = 3, pad = "0"), "_", str_to_lower(string = str_c(str_sub(string = projection, start = 1, end = 4) , "_", str_sub(string = projection, start = 5, end = 6))), "_fw_fh.", str_sub(string = year, start = 2, end = 4))
      writeLines(text = header, con = file)
      write.table(x = s_db, file = file, na = "-99.9", sep = ",", dec = ".", append = TRUE, quote = FALSE, col.names = FALSE, row.names = FALSE)
      
      if (any(is.na(s_db |> select(-"DR", -"EV24")))) warning("incomplete data for file: ", file, "\n")
      
    }
  }
}

# --- return of procedure ---

TIMPRGEND <- Sys.time()
TIMPRGCALC <- as.numeric(difftime(time1 = TIMPRGEND, time2 = TIMPRGSTART, units = "secs"))
message(paste0("\nProgram R-PROG successfully ended in ", floor(TIMPRGCALC / 60), " minutes and ", ceiling(TIMPRGCALC %% 60), " seconds"))
q(save = "no")
