feat(heatmap): add objective curve

Refs #14, #16

R/calcRiskHeatMap.R

@@ -2,26 +2,33 @@ calcRiskHeatMap <- function(con, ruleset, station, level, lake) {
   iObjective <- which(objectives$station == station & objectives$level == level)
   cdfs <- tbl(con, "CDFs")
   ruleset <- as.character(ruleset)
-  cdf <- cdfs %>% 
-    filter(id_ruleset == {{ ruleset }}, 
+  cdf <- cdfs %>%
+    filter(id_ruleset == {{ ruleset }},
            id_objective == {{ iObjective }},
-           id_lake == {{ lake }}) %>% 
+           id_lake == {{ lake }}) %>%
     collect()
   cdf$id_cal_day <- as.factor(cdf$id_cal_day)
-  discretStorage <- seq(0, rvgest::lakes[lake, "max"] - rvgest::lakes[lake, "min"], 1) * 1E6
-  listDiscretProbs <- lapply(unique(cdf$id_cal_day), function(x) {
-    dataLength <- length(which(cdf$id_cal_day == x))
+  cdf$V <- cdf$V / 1E6 # -> hm3
+  if(!objectives$flood[iObjective]) cdf$prob <- 1 - cdf$prob
+  discretStorage <- seq(0, rvgest::lakes[lake, "max"] - rvgest::lakes[lake, "min"], 1)
+  caldays <- unique(cdf$id_cal_day)
+  listDiscretProbs <- lapply(seq.int(length(caldays)), function(i) {
+    calday <- caldays[i]
+    dataLength <- length(which(cdf$id_cal_day == calday))
     if(dataLength > 1) {
-      l <- approx(cdf$V[cdf$id_cal_day == x], cdf$prob[cdf$id_cal_day == x], discretStorage, rule = 2)
-      data.frame(day = x, V = l$x, prob = l$y)
+      l <- approx(cdf$V[cdf$id_cal_day == calday], cdf$prob[cdf$id_cal_day == calday], discretStorage, rule = 2)
+      data.frame(julian = i, day = calday, V = l$x, prob = round(l$y, digits = 1))
     } else if(dataLength == 1){
-      data.frame(day = x, V = discretStorage, prob = 0)
+      data.frame(julian = i, day = calday, V = discretStorage, prob = 0)
     } else {
       warning("No data for lake", lake, " at ", station, " for objective ", level)
     }
   })
   # if(length(listDiscretProbs) == 365) {
-    return(do.call(rbind, listDiscretProbs))
+    df <- do.call(rbind, listDiscretProbs)
+    df$V <- round(df$V + rvgest::lakes[lake, "min"])
+    df$julian <- 1:nrow(df)
+    return(df)
   # } else {
   #   return(NULL)
   # }

R/getObjectiveStorage.R

@@ -3,22 +3,26 @@
 #' @param date [Date] (default current date)
 #'
 #' @return Named [numeric] vector with objective storage for each reservoir
-#' @export
 #'
 #' @examples
 #' getObjectiveStorage()
 getObjectiveStorage <- function(date = Sys.Date()) {
+  date_julian <- which(format(date, "%m%d") == objectiveStorageCurves$day)
+  objectiveStorageCurves[date_julian, -1]
+}
+
+getObjectiveStoragePivots <- function() {
   pivots <- list(
-    data.frame(mmdd = c('1101', '0201', '0301', '0701'),
-               S = c(24.5, 84., 130., 170.)),
-    data.frame(mmdd = c('1108', '1208', '0208', '0408', '0708'),
-               S = c(18.7, 40., 95., 175., 207.8)),
     data.frame(mmdd = c('1101', '1231', '0131', '0701'),
                S = c(12.3, 34., 50., 80.)),
+    data.frame(mmdd = c('1108', '1208', '0208', '0408', '0708'),
+               S = c(18.7, 40., 95., 175., 207.8)),
+    data.frame(mmdd = c('1101', '0201', '0301', '0701'),
+               S = c(24.5, 84., 130., 170.)),
     data.frame(mmdd = c('1101', '0101', '0201', '0301', '0701'),
                S = c(25, 100, 170, 260, 350))
   )
-  names(pivots) <- c("AUBE", "SEINE", "YONNE","MARNE")
+  names(pivots) <- rvgest::lakes$name
   caldays <- format(seq(as.Date("2001-01-01"), as.Date("2001-12-31"), 1), "%m%d")
   pivots <- lapply(pivots, function(x) {
     # sort array by date
@@ -31,9 +35,28 @@ getObjectiveStorage <- function(date = Sys.Date()) {
     x$julian[nrow(x)] <- x$julian[nrow(x)] + 365
     x
   })
-  date_julian <- which(format(date, "%m%d") == caldays)
-  unlist(lapply(pivots, function(x) {
-    approx(x$julian, x$S, date_julian)$y
-  }))
 }
 
+objectiveStoragePivots <- getObjectiveStoragePivots()
+
+getObjectiveStorageCurve <- function(lake) {
+  caldays <- format(seq(as.Date("2001-01-01"), as.Date("2001-12-31"), 1), "%m%d")
+  data.frame(day = as.factor(format(seq(as.Date("2001-01-01"), as.Date("2001-12-31"), 1), "%m%d")),
+             S = sapply(seq.int(365), function(i) {
+               approx(objectiveStoragePivots[[lake]]$julian, objectiveStoragePivots[[lake]]$S, i)$y
+             }))
+}
+
+getObjectiveStorageCurves <- function() {
+  l <- lapply(rvgest::lakes$name, function(lake) {
+    df <- getObjectiveStorageCurve(lake)
+    names(df)[2] <- lake
+    df
+    })
+  df <- do.call(cbind, l)
+  df <- df[, c(1, 2 * seq.int(nrow(rvgest::lakes)))]
+  df$julian = seq(365)
+  return(df)
+}
+
+objectiveStorageCurves <- getObjectiveStorageCurves()

R/mod_instant_risk_overview.R

@@ -4,9 +4,9 @@
 #'
 #' @param id,input,output,session Internal parameters for {shiny}.
 #'
-#' @noRd 
+#' @noRd
 #'
-#' @importFrom shiny NS tagList 
+#' @importFrom shiny NS tagList
 mod_instant_risk_overview_ui <- function(id){
   ns <- NS(id)
   choices <- seq(length(rvgest::rulesets$rules))
@@ -30,7 +30,7 @@ mod_instant_risk_overview_ui <- function(id){
       lapply(seq.int(nrow(rvgest::lakes)), function(i) {
         column(width = 3,
           sliderInput(ns(paste0("V", i)),
-                      label = rvgest::lakes$name[i],
+                      label = paste(rvgest::lakes$name[i], "lake (hm3)"),
                       value = round(getObjectiveStorage()[rvgest::lakes$name[i]]),
                       min = rvgest::lakes$min[i],
                       max = rvgest::lakes$max[i]))
@@ -39,17 +39,17 @@ mod_instant_risk_overview_ui <- function(id){
     # tableOutput(ns("table"))
     plotOutput(ns("plot"))
   )
-  
+
 }
-    
+
 #' instant_risk_overview Server Function
 #'
-#' @noRd 
+#' @noRd
 #' @import ggplot2
 mod_instant_risk_overview_server <- function(id, con) {
   moduleServer(id, function(input, output, session){
     ns <- session$ns
-    
+
     df <- reactive({
       storages <- c(input$V1, input$V2, input$V3, input$V4)
       names(storages) <- lakes$name
@@ -64,10 +64,10 @@ mod_instant_risk_overview_server <- function(id, con) {
         coord_flip())
   })
 }
-    
+
 ## To be copied in the UI
 # mod_instant_risk_overview_ui("instant_risk_overview_ui_1")
-    
+
 ## To be copied in the server
 # mod_instant_risk_overview_server("instant_risk_overview_ui_1")
- 
+

R/mod_one_objective_focus.R

@@ -4,9 +4,9 @@
 #'
 #' @param id,input,output,session Internal parameters for {shiny}.
 #'
-#' @noRd 
+#' @noRd
 #'
-#' @importFrom shiny NS tagList 
+#' @importFrom shiny NS tagList
 #' @importFrom rvgest objectives
 mod_one_objective_focus_ui <- function(id){
   ns <- NS(id)
@@ -50,10 +50,10 @@ mod_one_objective_focus_ui <- function(id){
     uiOutput(ns("plots"))
   )
 }
-    
+
 #' one_objective_focus Server Functions
 #'
-#' @noRd 
+#' @noRd
 mod_one_objective_focus_server <- function(id, con){
   moduleServer( id, function(input, output, session){
     ns <- session$ns
@@ -66,18 +66,21 @@ mod_one_objective_focus_server <- function(id, con){
         dfDiscretProbs <- calcRiskHeatMap(con, input$ruleset, objective$station, objective$level, lake)
         if(!is.null(dfDiscretProbs)) {
           breaks <- sapply(iBreaks, function(x) levels(dfDiscretProbs$day)[x])
-          p <- ggplot(dfDiscretProbs, aes(x = day, y = V / 1e6)) + 
+          objCurve <- data.frame(day = objectiveStorageCurves$day, V = round(objectiveStorageCurves[[lake]]))
+          objCurve <- objCurve[seq.int(1, 365,3),]
+          p <- ggplot(dfDiscretProbs, aes(x = day, y = V)) +
             geom_tile(aes(fill = prob)) +
             scale_fill_continuous(low = "green", high = "red", name = "Failure probability") +
+            geom_point(data = objCurve, color = "black", size = 1, shape = 21, fill = "grey") +
             ggtitle(paste("Lake", lake)) +
             scale_x_discrete("Calendar days", breaks = breaks, labels = labels) +
             scale_y_continuous(name="Reservoir storage (hm3)")
           output[[paste("plot", lake, sep = "_")]] <- renderPlot({p},width = 640, height = 280)
         }
       })
-      
+
     })
-    
+
     output$plots <- renderUI({
       objective <- rvgest::objectives[rvgest::objectives$station == input$station & rvgest::objectives$level == input$level,]
       plot_output_list <- lapply(objective$lakes[[1]]$name, function(lake) {
@@ -88,9 +91,9 @@ mod_one_objective_focus_server <- function(id, con){
      })
   })
 }
-    
+
 ## To be copied in the UI
 # mod_one_objective_focus_ui("one_objective_focus_ui_1")
-    
+
 ## To be copied in the server
 # mod_one_objective_focus_server("one_objective_focus_ui_1")

irmara.Rproj

@@ -12,6 +12,9 @@ Encoding: WINDOWS-1252
 RnwWeave: Sweave
 LaTeX: pdfLaTeX
 
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
 BuildType: Package
 PackageUseDevtools: Yes
 PackageInstallArgs: --no-multiarch --with-keep.source

vignettes/debugging/HeatMap.Rmd

+---
+title: "Untitled"
+author: "David Dorchies"
+date: "24/04/2021"
+output: html_document
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+
+# Setup
+
+```{r cars}
+library(irmara)
+library(ggplot2)
+con <- DBI::dbConnect(RSQLite::SQLite(), dbname = irmara:::app_sys("cdf", "Qgen_5000y_unbiased_median.sqlite"))
+objectiveStorageCurves <- irmara:::objectiveStorageCurves
+```
+
+#
+
+```{r}
+ruleset <- 1
+iObjective <- 1
+objective <- rvgest::objectives[iObjective,]
+lake <- objective$lakes[[1]]$name[1]
+```
+
+```{r}
+objCurve <- data.frame(day = objectiveStorageCurves$day, V = round(objectiveStorageCurves[[lake]]))
+objCurve <- objCurve[seq.int(1, 365,3),]
+head(objCurve)
+```
+
+
+```{r}
+indexMonths <- c(1,31,28,31,30,31,30,31,31,30,31,30)
+labels <- unlist(strsplit("JFMAMJJASOND", NULL))
+iBreaks <- sapply(seq.int(12), function(i) sum(indexMonths[1:i]))
+```
+
+```{r}
+dfDiscretProbs <- irmara:::calcRiskHeatMap(con, ruleset, objective$station, objective$level, lake)
+breaks <- sapply(iBreaks, function(x) levels(dfDiscretProbs$day)[x])
+head(dfDiscretProbs)
+```
+
+```{r}
+p <- ggplot(dfDiscretProbs, aes(x = day, y = V)) +
+            geom_tile(aes(fill = prob)) +
+            scale_fill_continuous(low = "green", high = "red", name = "Failure probability") +
+            geom_point(data = objCurve, color = "black", size = 1, shape = 21, fill = "grey") +
+            ggtitle(paste("Lake", lake)) +
+            scale_x_discrete("Calendar days", breaks = breaks, labels = labels) +
+            scale_y_continuous(name="Reservoir storage (hm3)")
+p
+```
+