#' Import Abiotic data
#' 
#' Function to import abiotic data using the SM data download toolbox webservices (https://www.scheldemonitor.be/dataproducts/en/download/)
#' Written by Wout Van Echelpoel (UGent)
#' Wout.VanEchelpoel@UGent.be - Sep 2021
#' for the T2021 ScheldeMonitor report 
#' https://github.com/scheldemonitor/T2021/blob/main/00_Algemeen/Scripts/Function_ImportData.R
#'
#' @param param parameter id or list of parameter ids
#' @param start start year
#' @param end end year
#'
#' @return a dataframe with the available data for the parameters and years selected
#' @export
#'
#' @examples 
#' # Define parameters (Scheldemonitor-code)
#' v.par.dutch <-c(9694,9695) # hoog- en laagwater in NAP+cm # Nederlandse codes
#' 
#' # download data
#' df.dat.dutch <- f.importAbioticData(v.par.dutch, start = 2010, end = 2023)
f.importAbioticData <- function(param, start = 2015, end = 2021){
  # Define time period
  v.years <- c(start:end)
  
  # Define fixed elements of URL
  v.url <- c('http://geo.vliz.be/geoserver/wfs/ows?service=WFS&version=1.1.0&request=GetFeature&typeName=Dataportal%3Aabiotic_observations&resultType=results&viewParams=where%3Aobs.context+%26%26+ARRAY%5B1%5D+AND+standardparameterid+IN+%28',
             '%29+AND+%28%28datetime_search+BETWEEN+%27',
             '-01-01%27+AND+%27',
             '-12-32%27+%29%29%3Bcontext%3A0001%3Bloggedin%3A1&propertyName=stationname%2Clongitude%2Clatitude%2Cdatetime%2Cdepth%2Cparametername%2Cvaluesign%2Cvalue%2Cunit%2Cdataprovider%2Cdatasettitle%2Clod%2Cloq%2Ccategory%2Cseason%2Cclassunit%2Cclass%2Caphiaid%2Cdateprecision%2Cdatafichetitle%2Cdataficheid%2Cstandardparameterid&outputFormat=csv')
  
  # Define parameter-specific section of URL
  n.par <- paste(param, collapse = '%5C%2C')
  
  # Define temporary list for output
  lst.dat <- list()
  
  # Import data for each year
  for (i in 1:length(v.years)) {
    message(paste0('Import for year ', v.years[i]))
    # Compose year-specific URL and import data
    s.url <- paste0(v.url[1], n.par, v.url[2], v.years[i], v.url[3], 
                    v.years[i], v.url[4])
    lst.dat[[i]] <- data.frame(st_read(s.url, stringsAsFactors = F, quiet = T))
    
    # Print feedback to user
    if (nrow(lst.dat[[i]]) == 0) { 
      message(paste0('---No data for year ', v.years[i], '---'))
    } else {
      message(paste0('---Data for year ', v.years[i], ' added---'))
    }
  }
  
  # Generate dataframe from list
  df.dat <- do.call(rbind, lst.dat)
  
  # Generate output
  return(df.dat)
}


#' Check the amount of records
#' 
#' Check if per station each year has sufficient data for realiable analyses
#' written by Matthijs Gensen (HKV, Schelde in Beeld)
#' for the T2021 ScheldeMonitor report
#' https://github.com/scheldemonitor/T2021/blob/main/02_Hydrodynamiek/Function_Checkamount.R
#'
#' @param data a dataframe resulting from f.importAbioticData()
#' @param n.min minimum number of records per year
#'
#' @return a dataframe with only the valid years per station
#' @export
#'
#' @examples
#' 
#' # Define the minumun number a station must have per year to be valid
#' n.min <- 0.75 * 705 # Check HW, LW sepeartely (cyclus equals 12h25min, 705 values per year)
#' 
#' df.HW.dutch <- f.checkamount(df.HW.dutch, n.min)
f.checkamount <- function(data ,n.min){
  # Remove years with insufficient data
  for (station in unique(data$stationname)){
    df.station <- data[which(data$stationname %in% station),]
    for (y in unique(df.station$year)){
      df.stationyear <- df.station[which(df.station$year %in% y),]
      if (nrow(df.stationyear)>= n.min){
        if(!exists("df.checked")) {
          df.checked <- df.stationyear
        } else {
          df.checked <- bind_rows(df.checked, df.stationyear)
        }
      } else {
        print(paste('Insufficient data for', station, 'in', as.character(y)))
      }
    } 
  }
  return(df.checked)
}

#' Helper function to calculate statistic
#' 
#' Written by Wout Van Echelpoel (UGent) for the T2021 ScheldeMonitor report
#' Wout.VanEchelpoel@UGent.be - Jan 2022 
#' https://github.com/scheldemonitor/T2021/blob/main/00_Algemeen/Scripts/Function_TemporalAggregating.R
#'
#' @param data a dataframe
#' @param stat the statistic to calculate. Options are 'mean', 'sd', 'se', 'median', 'min', 'max', '5perc', '10perc', '90perc', '99perc','sum'
#'
#' @return the statistic
#' @export
#'
#' @examples 
f.calculateStatistic <- function(data, stat = 'mean'){
  # Determine if requested statistic is included
  if(!tolower(stat) %in% c('mean', 'sd', 'se', 'median', 'min', 'max', '5perc', 
                           '10perc', '90perc', '99perc','sum')){
    stop(paste0('Options for statistic are: (1) mean, (2) sd, (3) se, ', 
                '(4) median, (5) min, (6) max, (7) 5perc, (8) 10perc, ', 
                '(9) 90perc, (10) 99perc and (11) sum'), call. = F)
  }
  
  # Determine requested statistic
  if(tolower(stat) == 'mean'){
    n.stat <- mean(data$value, na.rm = T)
  } else if(tolower(stat) == 'sd'){
    n.stat <- sd(data$value, na.rm = T)
  } else if(tolower(stat) == 'se'){
    n.stat <- sd(data$value, na.rm = T) / sqrt(nrow(data[!is.na(data$value)]))
  } else if(tolower(stat) == 'median'){
    n.stat <- median(data$value, na.rm = T)
  } else if(tolower(stat) == 'min'){
    n.stat <- min(data$value, na.rm = T)
  } else if(tolower(stat) == 'max'){
    n.stat <- max(data$value, na.rm = T)
  } else if(tolower(stat) == '5perc'){
    n.stat <- quantile(data$value, probs = 0.05, na.rm = T)
  } else if(tolower(stat) == '10perc'){
    n.stat <- quantile(data$value, probs = 0.10, na.rm = T)
  } else if(tolower(stat) == '90perc'){
    n.stat <- quantile(data$value, probs = 0.90, na.rm = T)
  } else if(tolower(stat) == '99perc'){
    n.stat <- quantile(data$value, probs = 0.99, na.rm = T)
  } else if(tolower(stat) == 'sum'){
    n.stat <- sum(data$value, na.rm = T)
  }
  # Return calculated statistic
  return(n.stat)
}


#' Define function for yearly aggregation per station
#' 
#' Written by Wout Van Echelpoel (UGent) for the T2021 ScheldeMonitor report
#' Wout.VanEchelpoel@UGent.be - Jan 2022 
#' https://github.com/scheldemonitor/T2021/blob/main/00_Algemeen/Scripts/Function_TemporalAggregating.R
#' 
#' @param data a dataframe. Select the columns to aggregate.
#' @param stat the statistic to calculate
#'
#' @return a dataframe  the statistic calculated
#' @export
#'
#' @examples 
#' df.mean <- f.yearlyPerStation(df.HW.dutch[c("year", "stationname", "value")], 'mean') #calculate the parameter per year and station
#' df.temp <- f.yearlyPerStation(df.HW.dutch[c("year", "stationname", "value")], '99perc')
f.yearlyPerStation <- function(data, stat = 'mean'){
  # Determine if all required columns are present
  if(length(which(c('year', 'stationname', 'value') %in% names(data))) != 3){
    stop('Required input: (1) year, (2) stationname and (3) value')
  } else {
    message(paste0('Calculate ', stat, ' per year and per station'))
  }
  
  # Remove stations without stationname
  message(paste0("---Removed ", sum(is.na(data$stationname)), 
                 " observations without stationname---"))
  data <- data[!is.na(data$stationname), ]
  
  # Split data conditional to station and year
  lst.dat <- split(data, list(data$stationname, data$year), drop = T)
  
  # Define output list
  lst.yr <- list()
  
  # For each subset, calculate statistic
  for (i in c(1:length(lst.dat))){
    # Add general information to output
    v.names <- c('latitude', 'longitude', 'stationname', 'parametername',
                 'unit', 'segment', 'zone', 'year')
    lst.yr[[i]] <- lst.dat[[i]][1, names(lst.dat[[i]]) %in% v.names]
    
    # Determine requested statistic and add to output
    lst.yr[[i]]$value <- f.calculateStatistic(lst.dat[[i]], stat = stat)
  }
  
  # Merge individual statistics and return sorted overview
  df.yr <- do.call('rbind', lst.yr)
  message(paste0('---Yearly ', stat, ' per station calculated---'))
  return(df.yr[order(df.yr$stationname, df.yr$year), ])
}