###########################################################
# Schelde water quality interface                         #
#                                                         #
# Auteurs: Willem Stolte                                  #
#                                                         #
#                                                         #
#                                                         #
# Datum : 2019-11-07                                      #
# Bedrijf: Deltares                                       #
# Licentie: GNU General Public License                    #
#                                                         #           
# Contact : Willem Stolte of Marcel Taal                  #
# Email : willem.stolte@deltares.nl                       #
#                                                         #
########################################################### 


##=PREPARATION========

##==legend=====

## COMMENTS
## ##      = comment
## #       = outcommented code
## #!#     = should be checked if location of script changes 
## #?#     = to do
## ##==    = block separator

## OBJECTS
## fnXXX   = filename
## dirXXX  = directory (path) / foldername
## oXXX    = loaded object
## dfXXX   = dataframe      
## dtXXX   = datatable 

##==install packages and open libraries ====

library(shiny)
library(tidyverse)
library(shinydashboard)
library(leaflet)
# library(plotly)
library(lubridate)
library(shinyjs)
library(shinycssloaders)
library(V8)
library(htmlwidgets)
library(htmltools)
library(mapview)
library(readxl)
library(scales)
# library(xtable)

##==settings=======

# options(shiny.display.mode="showcase")
# options(scipen = 999)  ## no scientific numbering
# options(shiny.reactlog = T)

##== load functions ====


##==set paths =======

#!# directories

#!# files

# fysChemDataFile <- "data/ow_algemeen.csv"
fysChemDataFile <- "data/ddl_ow_algemeen.csv" 

##== set variables =======

# misschien gebruiken voor type stof?
# variables = reactiveValues(soortgroep = NULL, soortgroepnaam = NULL, watertype = NULL)

##==load files=========

df <- read_delim(fysChemDataFile, delim = ";", guess_max = 50000) %>%
  filter(!is.na(value), value != 0)

# test
if("christmas" == "easter"){
  subdf <- df %>%
    filter(year  >= 1998,  year <= 2013) %>%
    filter(parameterlabel %in% "ZS mg/l") %>%
    filter(station %in% "Terneuzen boei 20")
  coefficients <- subdf %>% 
    group_by(year, month, parameter, station) %>% summarize(monthlyMean = mean(value, na.rm = T)) %>%
    group_by(year, parameter, station) %>% summarize(yearlyMean = mean(monthlyMean, na.rm = T)) %>%
    ungroup() %>% group_by(station, parameter) %>%
    do(trend = lm(yearlyMean ~ year, data = .)) %>%
    broom::tidy(trend) %>%
    filter(term == "year")
  
  View(coefficients)
  
  subdf %>%
    group_by(year, month, parameter, station) %>% summarize(monthlyMean = mean(value, na.rm = T)) %>%
    group_by(year, parameter, station) %>% summarize(yearlyMean = mean(monthlyMean, na.rm = T)) %>%
    ungroup() %>%
    nest(station, parameter) %>% 
    mutate(
      fit = map(data, ~ lm(yearlyMean ~ year, data = .x)),  # S3 list-col
      tidied = map(fit, broom::tidy)
    ) %>% 
    unnest(tidied) %>%
    filter(term == "year")
  
}


# "ZS mg/l"

##==make catalogue=========

substances = unique(df$parameterlabel)
stations = unique(df$station)
stationOrder = unname(unlist(unique(df[order(-df$lat + df$lng),]$station)))
df$station <- factor(df$station, levels = stationOrder)

##== set functions=====

## no functions yet

##==USER INTERFACE==================

header  <- dashboardHeader(title = "Westerschelde general water quality")

#=== sidebar menu ===========================
sidebar <- dashboardSidebar(
  sidebarMenu(id = "side_tabs",
              img(src = "deltares_logo.png", width = "175px"),
              menuItem(text = "Time series", tabName = "timeseries", icon = icon("chart-line")),
              menuItem(text = "Correlations", tabName = "correlations", icon = icon("project-diagram")),
              # menuItem(text = "Brush", tabName = "brush", icon = icon("paint-brush")),
              menuItem(text = "Info", tabName = "info", icon = icon("book-reader")),
              numericInput("startYear", label = "Start:", 2000, min = 1970, max = 2018),
              numericInput("endYear", label = "End:", 2018, min = 1970, max = 2018),
              checkboxInput("removeOutliers", label = "remove outliers?", value = T),
              wellPanel(
                downloadButton("dnld", label = "Download png")
              )
  )
)
## rows have a grid width of 12, so a box with width = 4, takes up one third of the space
## tops will be lined out, bottoms not
## heights are in pixels.. 

body    <- dashboardBody(
  
  tabItems(
    ##===time series tab==============
    tabItem(tabName = "timeseries",
            fluidRow(
              box(title = "parameter",
                  solidHeader = T,
                  status = "success",
                  collapsible = T,
                  collapsed = F,
                  width = 4,
                  selectInput("parameter", label = "Parameter:", selected = "NO3 Nnf mg/l", choices = substances, multiple = TRUE),
                  selectInput("location", label = "Location:", choices = stations, selected = 'Hansweert geul', multiple = TRUE),
                  checkboxInput("facetStation", label = "split stations", value = T)
              ),
              box(title = "lines",
                  solidHeader = T,
                  status = "success",
                  collapsible = T,
                  collapsed = F,
                  width = 4,
                  checkboxInput("trendline", label = "Trendline"),
                  checkboxInput("pointline", label = "Pointline"),
                  checkboxInput("loessline", label = "Loessline"),
                  conditionalPanel(
                    condition = "input.loessline",                  collapsed = F,
                    
                    sliderInput("lspan", label = NULL, min = 0, max = 1, step = 0.05, value = 0.7)
                  )
              ),
              box(title = "statistics",
                  solidHeader = T,
                  status = "success",
                  collapsible = T,
                  collapsed = F,
                  width = 4,
                  checkboxInput("percSummer", label = "10 & 90 percentile summer"),
                  checkboxInput("percWinter", label = "10 & 90 percentile winter"),
                  checkboxInput("boxplot", label = "Boxplot")
              ),
              tabBox(
                # status = "info",
                # collapsible = T,
                width = 12,
                # tabItems(
                tabPanel("plot",
                         # textOutput("height"),
                         uiOutput("ui_plot")
                         # plotOutput("facetplot", height = '200px')  #, width = '100%', height = '600px'
                ),
                tabPanel("coefficients",
                         tableOutput("coefTable")
                )
                # )
                
              )
            )
    ),
    
    ##== correlations tab ==============
    tabItem(tabName = "correlations",
            fluidRow(
              box(title = "Correlations",
                  solidHeader = T,
                  status = "success",
                  collapsible = F,
                  collapsed = F,
                  width = 2,
                  selectInput("paramx", label = "Parameter X:", selected = "SALNTT DIMSLS", choices = substances),
                  selectInput("paramy", label = "Parameter Y:", selected = "NO3 Nnf mg/l", choices = substances)
              ),
              box(
                status = "info",
                width = 8,
                plotOutput("scatterplot", width = '100%', height = )
              )
            )
    ),
    
    ##========= brush tab == in develepment =============
    # tabItem(tabName = "brush",
    #         fluidRow(
    #           box(title = "Brush",
    #               solidHeader = T,
    #               status = "success",
    #               collapsible = T,
    #               collapsed = F,
    #               width = 12,
    #               p("In development")
    #           )
    #         )
    # ),
    
    ##====achtergrondinformatie pagina ====
    tabItem(tabName = "info",
            fluidRow(
              box(title = "Uitleg",
                  solidHeader = T,
                  status = "success",
                  collapsible = T,
                  collapsed = F,
                  width = 12,
                  # h3("Handleiding AqMaD"),
                  # downloadButton("handleiding_aqmad", "Download handleiding voor AqMaD"),
                  includeMarkdown("include.md")
              ),
              box(title = "Download achtergrondinformatie",
                  solidHeader = T,
                  status = "success",
                  collapsible = T,
                  collapsed = F,
                  width = 12,
                  # h3("Handleiding AqMaD"),
                  # downloadButton("handleiding_aqmad", "Download handleiding voor AqMaD"),
                  p("De R-code voor deze tool is op aanvraag beschikbaar bij Willem Stolte (Deltares) via:"),
                  p("willem.stolte(at)deltares.nl")
              )
            )
    )
  )
)

ui <- dashboardPage(skin = "green", header, sidebar, body,  tags$head(tags$link(rel = "stylesheet", type = "text/css", href = "custom.css"), tags$style(".leaflet-top {z-index:999!important;}"))) ## tags$head etc is noodzakelijk om de zoombuttons van leaflet achter de dropdown box te krijgen. 
# ui <- dashboardPage(header, sidebar, body,  tags$head(tags$style(".leaflet-top {z-index:999!important;}"))) ## tags$head etc is noodzakelijk om de zoombuttons van leaflet achter de dropdown box te krijgen. 

##==SET SERVER=========


server <- function(input, output, session) {
  
  remove_outliers <- function(x, na.rm = TRUE, ...) {
    qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
    H <- 1.5 * IQR(x, na.rm = na.rm)
    y <- x
    y[x < (qnt[1] - H)] <- NA
    y[x > (qnt[2] + H)] <- NA
    y
  }
  
  plot_height <- reactive({
    max(200, 200*length(input$parameter) - 30*length(input$parameter))
  })
  
  ##==UI OPBOUW=========================================
  
  #!# When height is to be based on number of parameters 
  #?# Werkt niet zoals het hoort
  output$ui_plot <- renderUI({
    plotOutput("facetplot", height = paste0(plot_height(), "px"),
               width = "100%"
               )
  })
  
  output$height <- renderText(
    paste0(plot_height(), "px")
  )
  
  ##== OUTPUT =========================================
  
  subdf <- reactive({
    subdf <- df %>%
      filter(year  >= input$startYear,  year <= input$endYear) %>%
      filter(parameterlabel %in% input$parameter) %>%
      filter(station %in% input$location) %>% ungroup()
    if(input$removeOutliers) {
      subdf <- subdf %>%
        # group_by(parameterlabel, station) %>%
        mutate(value = remove_outliers(value))
    }
    subdf$season <- factor(subdf$season, levels = c('summer', 'winter'))
    return(subdf)
  })
  
  
  #== timeplot ============
  output$facetplot <- renderPlot({
    
    ylabel <- paste(unique(subdf()$parameterlabel), sep = ", ")
    
    if(nrow(subdf())==0) return(NULL)

    if(input$trendline){
      
      coefficients <- subdf() %>%
        group_by(year, month, parameter, station) %>% summarize(monthlyMean = mean(value, na.rm = T)) %>%
        group_by(year, parameter, station) %>% summarize(yearlyMean = mean(monthlyMean, na.rm = T)) %>%
        ungroup() %>%
        nest(-station, -parameter) %>% 
        mutate(
          fit = map(data, ~ lm(yearlyMean ~ year, data = .x)),  # S3 list-col
          tidied = map(fit, broom::tidy)
        ) %>% 
        unnest(tidied) %>% select(-data, -fit) %>%
        filter(term == "year" & p.value < 0.05)
      
      trenddata <- subdf() %>% 
        filter(station %in% coefficients$station & parameter %in% coefficients$parameter) %>%
        group_by(year, month, parameterlabel, station) %>% summarize(monthlyMean = mean(value, na.rm = T)) %>%
        group_by(year, parameterlabel, station) %>% summarize(yearlyMean = mean(monthlyMean, na.rm = T)) %>%
        ungroup() %>% filter(!is.na(station))
      
      if(length(trenddata$year) > 0){
        trenddata$datetime <- as_datetime((paste(as.character(trenddata$year), "07", "01")))
      }
      
    }

    # make base plot object
    pl <- ggplot(subdf(), aes(datetime, value))
    pl <- pl + geom_point(aes(color = station), alpha = 0.8)
    
    if(input$trendline){
      if(length(trenddata$year) >0){
        pl <- pl + 
          geom_smooth(data = trenddata, aes(datetime, yearlyMean, color = station), method='lm') #+
        # geom_text(data = coefficients, aes(as_datetime((paste(as.character(input$startYear), "12", "01"))), 0, label = paste("slope:", signif(estimate, 3), "|", "|", "p:", signif(p.value, 3))))
        #?# waarom komt deze bij alle parameters te staan? ook als maar een een
        #significante trend vertoont?
      }
    }
    
    if(input$loessline){
      # add loess line
      pl <- pl + geom_smooth(aes(color = station), method='loess', span = input$lspan, size = 1, n = 1000)
    }
    
    if(input$boxplot) {
      # add box and whiskers to plot
      pl <- pl + geom_boxplot(aes(x = as.POSIXct(paste0(year, '-07-01 00:00:00')),
                                  y = value,
                                  group = interaction(year, season), fill = season), 
                              alpha = 0.6,
                              outlier.size = 0)
    }
    
    if(input$percSummer){
      # calculate quantile statistics
      if(input$facetStation){
        subdf_year <- subdf() %>% 
          filter(season == 'summer') %>%
          group_by(year, season, parameterlabel, station) %>%
          summarize(
            perc90 = quantile(value, probs = 0.9, na.rm = TRUE),
            perc10 = quantile(value, probs = 0.1, na.rm = TRUE),
            value = mean(value))
      }
      
      if(!input$facetStation){
        subdf_year <- subdf() %>% 
          filter(season == 'summer') %>%
          group_by(year, season, parameterlabel) %>%
          summarize(
            perc90 = quantile(value, probs = 0.9, na.rm = TRUE),
            perc10 = quantile(value, probs = 0.1, na.rm = TRUE),
            value = mean(value))
      }
      pl <- pl + geom_ribbon(data = subdf_year, aes(as_datetime(paste(year, 7, 1)), ymin = perc10, ymax = perc90), fill = 'red', alpha = 0.2)
    }
    
    if(input$percWinter){
      # calculate quantile statistics
      if(input$facetStation){
        subdf_year <- subdf() %>% 
          filter(season == 'winter') %>%
          group_by(year, parameterlabel, station) %>%
          summarize(
            perc90 = quantile(value, probs = 0.9, na.rm = TRUE),
            perc10 = quantile(value, probs = 0.1, na.rm = TRUE),
            value = mean(value))
        # mutate(season <- factor(season, levels = c('summer', 'winter'))) %>%
      }
      
      if(!input$facetStation){
        subdf_year <- subdf() %>% 
          filter(season == 'winter') %>%
          group_by(year, season, parameterlabel) %>%
          summarize(
            perc90 = quantile(value, probs = 0.9, na.rm = TRUE),
            perc10 = quantile(value, probs = 0.1, na.rm = TRUE),
            value = mean(value))
        # mutate(season <- factor(season, levels = c('summer', 'winter'))) %>%
      }
      pl <- pl + geom_ribbon(data = subdf_year, aes(as_datetime(paste(year, 7, 1)), ymin = perc10, ymax = perc90), fill = 'blue', alpha = 0.2)
    }
    
    if(input$pointline) pl <- pl + geom_line(aes(color = station))
    # additional styling
    pl <- pl +
      theme_minimal() + 
      theme(
        # legend.position = "none",
        text = element_text(size = 16), 
        # legend.justification=c(1,1), legend.position = c(0.3,1), 
        panel.spacing.y = unit(0, "cm")#, 
        # panel.background = element_rect(fill = "transparent"),
        # panel.border = element_(color = "darkgrey")
      ) +
      scale_x_datetime(minor_breaks = date_breaks("1 year")) +
      xlab("year") #+ ylab(ylabel)
    
    if(input$facetStation) {pl <- pl + facet_grid(parameterlabel ~ station, scales = "free_y") + guides(color=FALSE)} else {
      pl <- pl + facet_grid(parameterlabel ~ ., scales = "free_y")
    }
    
    # return plot to ui
    print(pl)
    
  })
  
  
  #===coefficientsTable =================
  
  output$coefTable <- renderTable({
    # iris
    subdf() %>%
      group_by(year, month, parameter, station) %>% summarize(monthlyMean = mean(value, na.rm = T)) %>%
      group_by(year, parameter, station) %>% summarize(yearlyMean = mean(monthlyMean, na.rm = T)) %>%
      ungroup() %>%
      nest(-station, -parameter) %>%
      mutate(
        fit = map(data, ~ lm(yearlyMean ~ year, data = .x)),  # S3 list-col
        tidied = map(fit, broom::tidy)
      ) %>%
      unnest(tidied) %>% select(-data, -fit) %>%
      filter(term == "year" & p.value < 0.05)
  })
  
  
  
  
  #== scatterplot ============
  output$scatterplot <- renderPlot({
    
    tempdf <- df %>%
      filter(year  >= input$startYear,  year <= input$endYear) %>%
      filter(parameterlabel %in% c(input$paramx, input$paramy))
    
    if(input$removeOutliers) {
      tempdf <- tempdf %>%
        group_by(parameterlabel, station) %>%
        mutate(value = remove_outliers(value))
    }
    
    xydata <- tempdf %>%  
      group_by(station, parameter, year, datetime, season) %>%
      summarize(value = mean(value)) %>%
      spread(parameter, value)
    
    map = df %>% distinct(parameter, parameterlabel)
    
    PX = unname(unlist(map$parameter[map$parameterlabel==input$paramx]))
    PY = unname(unlist(map$parameter[map$parameterlabel==input$paramy]))
    
    sformula = paste(PX,  "~", PY)
    reg <- do.call("lm", list(formula = sformula, data = quote(xydata)))
    p.val = coef(summary(reg))[2,4]
    slope = coef(summary(reg))[2,1]
    intercept = coef(summary(reg))[1,1]
    
    
    # names(p.vals) = unique(xydata$station)
    
    # trenddata = xydata[xydata$station %in% names(p.vals)[p.vals < 0.05],]
    # if(length(trenddata$year) >0){
    #   trenddata$datetime <- as_datetime((paste(as.character(trenddata$year), "07", "01")))
    # }
    
    
    # if(weeklyAverage) {group_by(station, parameter, year, month, season) %>%}
    pl <- xydata %>%
      arrange(year) %>%
      ggplot(aes_string(PX, PY)) +
      geom_point(aes(color = station))
    # pl <- pl + geom_smooth(data = trenddata, aes(datetime, value, color = station), method='lm')
    if(p.val < 0.05) {
      pl = pl + geom_smooth(method = "lm")
      pl = pl + labs(caption =  paste("slope:", signif(slope, 3), "|", "intercept:", signif(intercept, 3), "|", "p:", signif(p.val, 3)))
    }
    pl + theme_light() + 
      xlab(sub("_","/" , input$paramx)) + ylab(sub("_","/" , input$paramy)) +
      theme(
        # legend.position = "bottom",
        text = element_text(size = 20)
      )
  })
  
  
  
  output$dnld = downloadHandler(
    filename = paste(input$location, cat(input$parameter, sep = "_"), '.png', sep=''),
    content = function(file) {
      ggsave(file, plot = last_plot(), device = "png") #, height = 5, width = 10
    }
  )
}

shinyApp(ui, server)