Commit a4a43a41 by Jan Wijffels

Connect the IRCEL calculation with the app

parent 3b0e41b0
...@@ -207,11 +207,11 @@ CitizenAir <- R6::R6Class("CitizenAir", ...@@ -207,11 +207,11 @@ CitizenAir <- R6::R6Class("CitizenAir",
} }
station station
}, },
getCombinedTimeseries = function(data = self$timeseries, type = c("rawdata", "minute", "hour", "day", "week", "month"), limit = FALSE){ getCombinedTimeseries = function(data = self$timeseries, type = c("rawdata", "min", "hour", "day", "week", "month"), limit = FALSE){
type <- match.arg(type) type <- match.arg(type)
x <- data x <- data
x <- setDT(x) x <- setDT(x)
if(type == "minute"){ if(type == "min"){
x <- x[, list(value = mean(value, na.rm=TRUE)), x <- x[, list(value = mean(value, na.rm=TRUE)),
by = list(type, timeseries_id, timeseries_label, phenomena_id, phenomena, time = fasttime::fastPOSIXct(format(time, "%Y-%m-%d %H:%M"), tz = "UTC"))] by = list(type, timeseries_id, timeseries_label, phenomena_id, phenomena, time = fasttime::fastPOSIXct(format(time, "%Y-%m-%d %H:%M"), tz = "UTC"))]
}else if(type == "hour"){ }else if(type == "hour"){
...@@ -235,7 +235,7 @@ CitizenAir <- R6::R6Class("CitizenAir", ...@@ -235,7 +235,7 @@ CitizenAir <- R6::R6Class("CitizenAir",
} }
x x
}, },
getCitizenData = function(type = c("meta", "timeseries", "phenomena", "fileinfo")){ getCitizenData = function(type = c("meta", "timeseries", "phenomena", "fileinfo", "timerange")){
type <- match.arg(type) type <- match.arg(type)
if(type == "meta"){ if(type == "meta"){
if(is.na(self$data$filename)){ if(is.na(self$data$filename)){
...@@ -276,6 +276,21 @@ CitizenAir <- R6::R6Class("CitizenAir", ...@@ -276,6 +276,21 @@ CitizenAir <- R6::R6Class("CitizenAir",
}else{ }else{
x <- data.frame(name = self$data$filename, file = self$data$file, stringsAsFactors = FALSE) x <- data.frame(name = self$data$filename, file = self$data$file, stringsAsFactors = FALSE)
} }
}else if(type == "timerange"){
x <- self$getCitizenData("timeseries")
x <- setDT(x)
startstop <- list()
if(nrow(x) > 0){
startstop$valid <- TRUE
startstop$start <- min(x$time)
startstop$end <- max(x$time)
}else{
startstop$valid <- FALSE
startstop$start <- as.POSIXct(character())
startstop$end <- as.POSIXct(character())
}
x <- startstop
x
} }
x x
} }
......
No preview for this file type
...@@ -87,32 +87,6 @@ reactive({ ...@@ -87,32 +87,6 @@ reactive({
}) })
``` ```
```{r}
## DEV APP SETTINGS
reactive <- reactiveValues()
plots <- reactiveValues()
mydata <- data.frame(date = as.POSIXct(seq.Date(Sys.Date(), Sys.Date()-999, by = "-1 day"), tz = "UTC"),
REF = rnorm(1000),
SENSOR = rnorm(1000))
linear_model <- lm(mydata$SENSOR ~ mydata$REF)
PI <- data.frame(predict(linear_model, interval="predict", level = 0.95))
colnames(PI) <- c("OLSQ", "lower PI", "upper PI")
CI <- data.frame(predict(linear_model, interval="confidence", level = 0.95))
colnames(CI) <- c("OLSQ", "lower CI", "upper CI")
mydata <- cbind(mydata, "OLSQ" = CI$OLSQ,"lower_CI" = CI$`lower CI`, "upper_CI"=CI$`upper CI`, "lower_PI" = PI$`lower PI`, "upper_PI" = PI$`upper PI`)
mydata$e <- mydata$OLSQ - mydata$SENSOR
mydata <- mydata[order(mydata[,1]), ]
summary <- data.frame(do.call(cbind, lapply(mydata, summary)))[-7,]
summary$date <- as.character(as.POSIXct(summary$date, origin="1970-01-01"))
summary <- rbind(summary,"# Valid"=colSums(!is.na(mydata)))
summary <- rbind(summary,"# NA's"=colSums(is.na(mydata)))
variables <- c("min", "25th Perc", "median", "mean","75th Perc", "max","# Valid", "# NA's")
summary <- cbind("variable"=variables,summary)
reactive[["linear_model"]] <- linear_model
reactive[["regData"]]<- mydata
reactive[["summaryRaw"]] <- summary
```
# Data {data-icon="fa-cloud-upload"} # Data {data-icon="fa-cloud-upload"}
```{r child="page01-introduction.Rmd", eval=TRUE} ```{r child="page01-introduction.Rmd", eval=TRUE}
......
...@@ -23,6 +23,7 @@ plot_timeseries <- function(x, ...){ ...@@ -23,6 +23,7 @@ plot_timeseries <- function(x, ...){
} }
``` ```
#### Officieel meetstation #### Officieel meetstation
```{r} ```{r}
...@@ -47,24 +48,33 @@ inputPanel( ...@@ -47,24 +48,33 @@ inputPanel(
ui ui
}), }),
radioButtons(inputId = "uiInput_aggregationlevel", label = "Vergelijk", radioButtons(inputId = "uiInput_aggregationlevel", label = "Vergelijk",
choices = c("Per minuut" = "minute", choices = c("Per minuut" = "min",
"Per uur" = "hour", "Per uur" = "hour",
"Per dag" = "day", "Per dag" = "day",
"Per maand" = "month", "Per maand" = "month",
"Op basis van ruwe data" = "rawdata"), "Op basis van ruwe data" = "rawdata"),
selected = "rawdata"), selected = "rawdata"),
#dateRangeInput(inputId = "uiInput_period", label = "Verander de periode", start = Sys.Date()-365, end = Sys.Date()), dateRangeInput(inputId = "uiInput_period", label = "Verander de periode", start = Sys.Date()-365, end = Sys.Date()),
actionButton(inputId = "uiInput_refresh", label = "Herbereken grafiek") actionButton(inputId = "uiInput_refresh", label = "Herbereken grafiek")
) )
observe({
userdata <- citizenair_userdata()
timerange <- userdata$appdata$getCitizenData("timerange")
if(timerange$valid){
updateDateRangeInput(session, inputId = "uiInput_period", start = as.Date(timerange$start), end = as.Date(timerange$end))
}
})
``` ```
```{r} ```{r}
observeEvent(input$uiInput_refresh, { observeEvent(input$uiInput_refresh, {
showNotification("Herberekening wordt uitgevoerd") showNotification("Herberekening wordt uitgevoerd")
input$uiInput_phenomena #input$uiInput_phenomena
input$uiInput_aggregationlevel #input$uiInput_aggregationlevel
input$uiInput_sheet #input$uiInput_sheet
input$uiInput_period #input$uiInput_period
#input$uiOutput_select_station_marker_click #input$uiOutput_select_station_marker_click
x <- APPDATA$getCombinedTimeseries(type = "rawdata") x <- APPDATA$getCombinedTimeseries(type = "rawdata")
x <- subset(x, type %in% "official" | (type %in% "citizen" & timeseries_id %in% input$uiInput_sheet)) x <- subset(x, type %in% "official" | (type %in% "citizen" & timeseries_id %in% input$uiInput_sheet))
...@@ -72,6 +82,66 @@ observeEvent(input$uiInput_refresh, { ...@@ -72,6 +82,66 @@ observeEvent(input$uiInput_refresh, {
x <- APPDATA$getCombinedTimeseries(x, type = input$uiInput_aggregationlevel, limit = TRUE) x <- APPDATA$getCombinedTimeseries(x, type = input$uiInput_aggregationlevel, limit = TRUE)
output$uiOutput_timeseries <- renderDygraph( output$uiOutput_timeseries <- renderDygraph(
plot_timeseries(x)) plot_timeseries(x))
mydata <- subset(x, phenomena_id %in% head(input$uiInput_phenomena, 1))
if(all(c("citizen", "official") %in% mydata$type) & input$uiInput_aggregationlevel %in% c("min", "hour", "day", "week", "month")){
## DATE/REF/SENSOR
mydata <- x
mydata$group <- ifelse(mydata$type %in% "official", "REF", "SENSOR")
mydata$date <- mydata$time
mydata <- dcast.data.table(data = mydata, formula = date ~ group, fun.aggregate = mean, na.rm=TRUE, value.var = "value")
#mydata <- data.frame(date = as.POSIXct(seq.Date(Sys.Date(), Sys.Date()-999, by = "-1 day"), tz = "UTC"),
# REF = rnorm(1000),
# SENSOR = rnorm(1000))
linear_model <- lm(mydata$SENSOR ~ mydata$REF)
PI <- data.frame(predict(linear_model, interval="predict", level = 0.95))
colnames(PI) <- c("OLSQ", "lower PI", "upper PI")
CI <- data.frame(predict(linear_model, interval="confidence", level = 0.95))
colnames(CI) <- c("OLSQ", "lower CI", "upper CI")
mydata <- cbind(mydata, "OLSQ" = CI$OLSQ,
"lower_CI" = CI$`lower CI`,
"upper_CI"=CI$`upper CI`,
"lower_PI" = PI$`lower PI`,
"upper_PI" = PI$`upper PI`)
mydata$e <- mydata$OLSQ - mydata$SENSOR
mydata <- mydata[order(mydata[,1]), ]
summary <- data.frame(do.call(cbind, lapply(mydata, summary)))[-7,]
summary$date <- as.character(as.POSIXct(summary$date, origin="1970-01-01"))
summary <- rbind(summary,"# Valid"=colSums(!is.na(mydata)))
summary <- rbind(summary,"# NA's"=colSums(is.na(mydata)))
variables <- c("min", "25th Perc", "median", "mean","75th Perc", "max","# Valid", "# NA's")
summary <- cbind("variable"=variables,summary)
# outputs of dataframe & summaries
output$data <- renderTable({
mydata <- timeAverage(mydata, avg.time=input$uiInput_aggregationlevel, data.thresh=85, statistic="mean")
mydata$date <- as.character(as.POSIXct(mydata$date, origin="1970-01-01"))
mydata
})
output$summaryRaw <- renderTable(summary)
# Scatterplot of linear regression + CI & PI
output$uiOutput_comparison <- renderPlot({
mydata <- mydata
linear_model <- linear_model
plot_PI_CI <- ggplot(data = mydata, aes(x = REF, y = SENSOR)) +
geom_ribbon(aes(ymin = lower_CI, ymax=upper_CI, fill = "CI"),alpha = 0.6) +
geom_ribbon(aes(ymin = lower_PI, ymax=upper_PI, fill = "PI"), alpha = 0.4) +
geom_point() +
geom_line(aes(y = OLSQ,col = "OLSQ")) +
geom_line(aes(y = REF,col = "y=x")) +
labs(title = paste("Scatterplot with regression estimates: intercept =",
format(summary(linear_model)$coefficients[1,1], digits = 3),
",slope =",format(summary(linear_model)$coefficients[2,1], digits = 3),
",R² =",format(summary(linear_model)$r.squared, digits = 3)),
x = "Reference value",
y = "Sensor") +
scale_fill_manual(name="",values=c("CI"="blue","PI"= "grey")) +
scale_colour_manual(name="", values = c("y=x"="red","OLSQ"="black"))
#plots[["scatter"]] <- plot_PI_CI
plot_PI_CI
})
}
}) })
``` ```
...@@ -85,6 +155,8 @@ observeEvent(input$uiInput_selectstation, { ...@@ -85,6 +155,8 @@ observeEvent(input$uiInput_selectstation, {
if(length(result$station) > 0){ if(length(result$station) > 0){
showNotification("Haalt data af vanaf de IRCEL webservice") showNotification("Haalt data af vanaf de IRCEL webservice")
} }
period <- lubridate::as.interval(as.POSIXct(input$uiInput_period[1] - 1), as.POSIXct(input$uiInput_period[1] + 1))
#official <- try(result$userdata$appdata$fetch_timeseries(result$station, time_span = period))
official <- try(result$userdata$appdata$fetch_timeseries(result$station)) official <- try(result$userdata$appdata$fetch_timeseries(result$station))
if(inherits(official, "try-error")){ if(inherits(official, "try-error")){
showNotification(sprintf("Er loopt iets fout met het ophalen van de data vanaf de IRCEL webservice: %s", error_message(official)), showNotification(sprintf("Er loopt iets fout met het ophalen van de data vanaf de IRCEL webservice: %s", error_message(official)),
...@@ -160,7 +232,7 @@ output$uiOutput_select_station <- renderLeaflet({ ...@@ -160,7 +232,7 @@ output$uiOutput_select_station <- renderLeaflet({
fillCol( fillCol(
tags$div( tags$div(
tags$h4("De Vlaamse Milieu Maatschappij heeft meetresultaten die aan de strengste kwaliteitseisen voldoen. Hieronder tonen tags$h4("De Vlaamse Milieu Maatschappij heeft meetresultaten die aan de strengste kwaliteitseisen voldoen. Hieronder tonen
we een eenvoudige statistische vergelijking tussen onze meetgegevens en die van jouw. Verander ook de VMM-meetplaats eens, een andere meetplaats komt misschien beter me de jouwe overeen!"), we een eenvoudige statistische vergelijking tussen onze meetgegevens en die van jouw. Verander ook de VMM-meetplaats eens, een andere meetplaats komt misschien beter met de jouwe overeen!"),
tags$h5("Onderstaande grafiek toont alle meetwaardes van uw toestel uitgezet tov alle meetwaardes van een toestel van de Vlaamse Milieu Maatschappij")), tags$h5("Onderstaande grafiek toont alle meetwaardes van uw toestel uitgezet tov alle meetwaardes van een toestel van de Vlaamse Milieu Maatschappij")),
fillRow(dropdownButton( fillRow(dropdownButton(
tags$h3("Selectie dagen van de week"), tags$h3("Selectie dagen van de week"),
...@@ -168,8 +240,8 @@ we een eenvoudige statistische vergelijking tussen onze meetgegevens en die van ...@@ -168,8 +240,8 @@ we een eenvoudige statistische vergelijking tussen onze meetgegevens en die van
choices = c("Alle dagen van de week" = "all", "Enkel weekdagen" = "weekdays", "Enkel weekend" = "weekend"), choices = c("Alle dagen van de week" = "all", "Enkel weekdagen" = "weekdays", "Enkel weekend" = "weekend"),
selected = "all"), selected = "all"),
tags$h3("Input voor de grafiek"), #tags$h3("Input voor de grafiek"),
sliderInput(inputId = "uiInput_slider", label = "TODO", value = 25, min = 5, max = 100, step = 1L, round = TRUE), #sliderInput(inputId = "uiInput_slider", label = "TODO", value = 25, min = 5, max = 100, step = 1L, round = TRUE),
circle = TRUE, status = "info", icon = icon("gear"), width = "300px", size = "sm", circle = TRUE, status = "info", icon = icon("gear"), width = "300px", size = "sm",
tooltip = tooltipOptions(title = "Verdere opties")), tooltip = tooltipOptions(title = "Verdere opties")),
plotOutput("uiOutput_comparison", height = "100%"), plotOutput("uiOutput_comparison", height = "100%"),
...@@ -177,29 +249,6 @@ we een eenvoudige statistische vergelijking tussen onze meetgegevens en die van ...@@ -177,29 +249,6 @@ we een eenvoudige statistische vergelijking tussen onze meetgegevens en die van
), ),
flex = c(NA, 2) flex = c(NA, 2)
) )
# Scatterplot of linear regression + CI & PI
output$uiOutput_comparison <- renderPlot({
req(c(reactive$regData,reactive$linear_model))
mydata <- reactive$regData
linear_model <- reactive$linear_model
plot_PI_CI <- ggplot(data = mydata, aes(x = REF, y = SENSOR)) +
geom_ribbon(aes(ymin = lower_CI, ymax=upper_CI, fill = "CI"),alpha = 0.6) +
geom_ribbon(aes(ymin = lower_PI, ymax=upper_PI, fill = "PI"), alpha = 0.4) +
geom_point() +
geom_line(aes(y = OLSQ,col = "OLSQ")) +
geom_line(aes(y = REF,col = "y=x")) +
labs(title = paste("Scatterplot with regression estimates: intercept =",
format(summary(linear_model)$coefficients[1,1], digits = 3),
",slope =",format(summary(linear_model)$coefficients[2,1], digits = 3),
",R² =",format(summary(linear_model)$r.squared, digits = 3)),
x = "Reference value",
y = "Sensor") +
scale_fill_manual(name="",values=c("CI"="blue","PI"= "grey")) +
scale_colour_manual(name="", values = c("y=x"="red","OLSQ"="black"))
plots[["scatter"]] <- plot_PI_CI
plot_PI_CI
})
``` ```
### Gevorderd overzicht van de vergelijking ### Gevorderd overzicht van de vergelijking
...@@ -213,18 +262,7 @@ fillCol( ...@@ -213,18 +262,7 @@ fillCol(
) )
``` ```
```{r}
# outputs of dataframe & summaries
output$data <- renderTable({
req(reactive[["regData"]])
mydata <- reactive[["regData"]]
avg.time <- "day"
mydata <- timeAverage(mydata, avg.time=avg.time, data.thresh=85,statistic="mean")
mydata$date <- as.character(as.POSIXct(mydata$date, origin="1970-01-01"))
mydata})
output$summaryRaw <- renderTable({reactive[["summaryRaw"]]})
output$summaryValid <- renderTable({reactive[["summaryValid"]]})
```
## Inputs {data-width=150} ## Inputs {data-width=150}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment