---
title: "Wrapper"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Wrapper}
%\VignetteEncoding{UTF-8}
%\VignetteEngine{knitr::rmarkdown}
editor_options:
chunk_output_type: console
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
# Installation
```{r setup, eval = FALSE}
### Optionally: specify GitHub Personal Access Token (GITHUB_PAT)
### See here why this might be important for you:
### https://kwb-r.github.io/kwb.pkgbuild/articles/install.html#set-your-github_pat
# Sys.setenv(GITHUB_PAT = "mysecret_access_token")
# Install package "remotes" from CRAN
if (! require("remotes")) {
install.packages("remotes", repos = "https://cloud.r-project.org")
}
# Install KWB package 'kwb.heatsine' from GitHub
remotes::install_github("KWB-R/kwb.heatsine")
```
# Data prerequisites
The sinus fit optimisation works for daily time series of a surface water
and a groundwater monitoring point, which follow the following requirements:
## File naming
The time series must be provided in a `.csv` (comma separated value) file, which
follows the naming convention:
temperature_``_``.csv
where `` needs to be exactly `groundwater` or `surface-water`. In case any
other value for `` is chosen the program will not work (e.g `SURFACE
WATER`). The value of the `` will be used for labeling the
monitoring points in the plots and can be almost any value but must not contain
`_` (underscore) or `" "` (space).
## Data structure
The data structure of each `.csv` file is very simple as it contains only the
two columns `date` (formatted: `YYYY-MM-DD`, e.g. `2020-06-23`) and `value`
(temperature data in degree Celsius) as shown exemplary below:
```{r eval = FALSE}
date,value
2020-06-23,12.3
2020-06-24,12.3
```
Please make sure that both columns `date` and `value` are all lowercase and
correctly spelled as the program will not work if this is not the case.
# Data preparation
## Import
The R package includes the following temporal times series data for testing:
```{r}
list.files(kwb.heatsine::extdata_file())
```
In this the example the files `temperature_surface-water_Txxsxx-mxxxxsxxx.csv`
and `temperature_groundwater_Txxxx3.csv` are used for testing the functionality
and imported by using the code below:
```{r}
## Load the R package
library(kwb.heatsine)
load_temp <- function(base_name) {
kwb.heatsine::load_temperature_from_csv(
kwb.heatsine::extdata_file(base_name)
)
}
data_sw <- load_temp("temperature_surface-water_Txxsxx-mxxxxsxxx.csv")
data_gw <- load_temp("temperature_groundwater_Txxxx3.csv")
```
## Visualise
Plot time series for whole time period (moveover on data points to select)
### Surface water
```{r}
kwb.heatsine::plot_temperature_interactive(data_sw)
```
### Groundwater
```{r}
kwb.heatsine::plot_temperature_interactive(data_gw)
```
## Select
Reduce time period to a full sinus period which is definde by the user input
(`date_end` is optional, if not given it is set to 'date_start' + 365.25 days).
```{r}
data_sw_selected <- kwb.heatsine::select_timeperiod(
data_sw,
date_start = "2015-10-10",
date_end = "2016-10-14"
)
data_gw_selected <- kwb.heatsine::select_timeperiod(
data_gw,
date_start = "2015-12-28",
date_end = "2016-12-26"
)
```
Plot selected datasets:
### Surface Water (selected)
```{r}
kwb.heatsine::plot_temperature_interactive(df = data_sw_selected)
```
### Grundwater (selected)
```{r}
kwb.heatsine::plot_temperature_interactive(df = data_gw_selected)
```
# Sinus optimisation
## Background
According to [Cross Validated (2013)](https://stats.stackexchange.com/questions/77543/how-do-i-get-the-amplitude-and-phase-for-sine-wave-from-lm-summary)
the sinus fit can not only be written as `y = y0 + amplitude * sin(x + phase_shift)`, but also
in a linear form `y = y0 + alpha * sin(x) + beta * cos(x)`. The latter approach is used
within this R package as it enables the calculation of both `amplitude` (`a`) and
`phase_shift` (`x0`) by using R's build-in `stats::lm()` function as exemplary shown below:
```{r eval = TRUE}
### Generate synthetic temperature data following sine curve
period_length <- 365 # assumed period length
day_number <- seq_len(period_length) # daily sequence from 1 to period_length
temperature <- 26.9188 + 2.123641 * sin(2 * pi * day_number / period_length + 0.6049163)
x <- 2 * pi * day_number / period_length
### Print first data points:
tibble::tibble(day_number = day_number, temperature = temperature, x = x)
### Now check the if both formulas yield (almost) equal results
fit.lm2 <- lm(temperature ~ sin(x) + cos(x))
coeffs <- coef(fit.lm2)
y0 <- coeffs[1L] # y0
alpha <- coeffs[2L] # alpha = sin(phi)
beta <- coeffs[3L] # beta = cos(phi)
a <- sqrt(alpha ^ 2 + beta ^ 2) # amplitude
x0 <- atan2(beta, alpha) # phase t
### Print parameters of sinus fit
tibble::tibble(
y0 = y0,
alpha = alpha,
beta = beta,
a = a,
x0 = x0
)
x <- 2 * pi * day_number / period_length
par(mfrow = c(1, 2))
### Formula option 1: y = y0 + amplitude * sin(x + phase)
plot(
x = day_number,
y = y0 + a * sin(x + x0),
type = "l",
lty = 1,
lwd = 3,
col = "gray",
main = "Overplotted Graphs",
xlab = "day_number",
ylab = "temperature"
)
### Formula option 2: y = y0 + alpha * sin(x) + beta * cos(x)
lines(
x = day_number,
y = y0 + alpha * sin(x) + beta * cos(x),
lwd = 2,
lty = 2,
col = "red",
)
### Check numeric differences between both approaches
plot(
x = day_number,
y = y0 + a * sin(x + x0) - (y0 + alpha * sin(x) + beta * cos(x)),
lwd = 3,
col = "gray",
main = "Difference",
xlab = "day_number",
ylab = "temperature"
)
```
As shown above the differences between both sinus calculation approaches are mostly
zero but sometimes show a minimal offset (which is close to R's floating point
accuracy of 10^`r .Machine$double.ulp.digits * log10(2)` , i.e. `r 10^(.Machine$double.ulp.digits * log10(2))`
for the machine the calculation was carried out). As a consequence the second
approach can and will be used for sinus fitting within this R package.
## Run
Sinus fit operation is performed by minimizing the RMSE error between observed
and simulated temperature curves.
```{r}
# Generate data frame with simulated and observed values
predictions <- kwb.heatsine::run_optimisation(data_sw_selected,
data_gw_selected,
retardation_factor = 1.8
)
```
## Analyse
### Interactive Plot
```{r}
kwb.heatsine::plot_prediction_interactive(predictions)
```
### Table: Traveltimes
```{r}
knitr::kable(predictions$traveltimes)
```
### Table: Goodness of Fit
While only RMSE is used for optimising the sinus fit, a number of goodness-of-fit
parameters is also calculated and shown in the table below. For details on each
of these parameters the use is referred to the function `hydroGOF::gof()` which is
used for calculating all these values.
```{r}
knitr::kable(predictions$gof)
```
### Table: Optimisation Parameters
```{r}
knitr::kable(predictions$paras)
```
### Table: Data
```{r}
# only first entries
knitr::kable(head(predictions$data))
```
## Export
Export results to `csv`:
```{r eval = FALSE}
write_to_csv <- function(df, file_base) {
readr::write_csv(df, path = file.path(
kwb.utils::createDirectory("csv"), paste0(file_base, ".csv")
))
}
write_to_csv(predictions$data, "sinus-fit_predictions")
write_to_csv(predictions$paras, "sinus-fit_parameters")
write_to_csv(predictions$gof, "sinus-fit_goodness-of-fit")
write_to_csv(predictions$traveltimes, "sinus-fit_traveltimes")
write_to_csv(predictions$residuals, "sinus-fit_residuals")
## Export plots:
save_to_html <- function(p, file_base) {
htmlwidgets::saveWidget(p, paste0(file_base, ".html"))
}
withr::with_dir(new = kwb.utils::createDirectory("plots"), code = {
save_to_html(
plot_temperature_interactive(data_sw),
"temperature_surface-water_time-series_full"
)
save_to_html(
plot_temperature_interactive(data_sw_selected),
"temperature_surface-water_time-series_selected"
)
save_to_html(
plot_temperature_interactive(data_gw),
"temperature_groundwater_time-series_full"
)
save_to_html(
plot_temperature_interactive(data_gw_selected),
"temperature_groundwater_time-series_selected"
)
save_to_html(
plot_prediction_interactive(predictions),
"temperature_prediction"
)
# save_to_html(
# plot_residuals_interactive(prediction_df, binwidth = 0.5),
# "temperature_prediction_residuals"
# )
})
```
# Session Info
## Plattform
```{r plattform, eval = TRUE, echo= FALSE}
environment <- sessioninfo::session_info()
knitr::kable(tibble::enframe(unlist(environment$platform)))
```
## Packages
```{r packages, eval = TRUE, echo= FALSE}
environment$packages
```
### Pandoc
```{r pandoc, eval = TRUE, echo= FALSE}
if(rmarkdown::pandoc_available()) {
data.frame(pandoc_directory = rmarkdown::pandoc_exec(),
pandoc_version = as.character(rmarkdown::pandoc_version()),
stringsAsFactors = FALSE)
} else {
print("No PANDOC installed!")
}