---
title: "Rosalia: an experimental research site to study hydrological processes in a forest catchment - data repository "
author: 
 - name: Josef Fürst, Hans Peter Nachtnebel
   affiliation: Institute for Hydrology and Water Management (HyWa), University of Natural Resources and Life Sciences Vienna (BOKU), Vienna, Austria.
   email: josef.fuerst@boku.ac.at
 - name: Josef Gasch
   affiliation: Forest Demonstration Centre, University of Natural Resources and Life Sciences Vienna  -  (BOKU), Vienna, Austria.
 - name: Reinhard Nolz, Michael Paul Stockinger, Christine Stumpp
   affiliation: Institute for Soil Physics and Rural Water Management (SoPhy), University of Natural Resources and Life Sciences Vienna (BOKU), Vienna Austria
 - name: Karsten Schulz
   affiliation: Institute for Hydrology and Water Management (HyWa), University of Natural Resources and Life Sciences Vienna (BOKU), Vienna, Austria.
   

date: "`r format(Sys.time(), '%B %d, %Y')`"
output:
  bookdown::html_document2: 
    df_print: paged
    code_folding: hide
    fig_width: 9
    number_sections: yes
    toc: yes
    toc_float: yes
  html_document:
    df_print: paged
    code_folding: hide
    fig_width: 9
    number_sections: yes
    toc: yes
    toc_float: yes
  html_notebook:
    code_folding: hide
    fig_width: 9
    number_sections: yes
    toc: yes
    toc_float: yes
always_allow_html: yes
editor_options: 
  chunk_output_type: inline
---
```{r, setup, include=FALSE}
# set this option in the first code chunk in the document
knitr::opts_chunk$set(cache = FALSE)
```

# Purpose of this document

This is a HTML document created by R Markdown which contains the technical description of the datasets related to **Fürst et al. (2020) Rosalia: an experimental research site to study hydrological processes in a forest catchment (submitted in August 2020 to Earth System Science Data (ESSD)**. To read in depth about the Rosalia site and the monitoring network implemented there, please see this publication.

This document has two purposes:

1. to demonstrate basic options of using the Rosalia time series data in R. For this purpose, the database was exported to a simple SQLite database, which is accessed by tools in the RSQLite package. Accessing the data by RSQLite is demonstrated in this document, and since SQLite is available for many languages and software tools, this can serve as a template for use in other environments.The R code chunks to perform a selection, computation or visualisation of the data can be turned on and off in the created HTML page.

2. to visually present the most important data and to comment on some of their characteristics  

Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>. This document generates dynamic graphics using the`dygraphs` package and therefore output is only possible to HTML.

If you opened the RosaliaESSD.Rmd file in RStudio, clicking the **Knit** button will generate a HTML document that includes both content as well as the output of any embedded R code chunks within the document. In the HTML file, the code chunks are hidden by default. By clicking the **Code** button next to a text paragraph, code chunks can be unfolded. The **Code** button right to the title turns on display of all code chunks in the document.

# Packages used

The essential R packages required for this document are DBI (DBMS interfaces), dygraphs (interactive timeseries graphs) and xts (time series). Some more packages are used as well, but could be replaced by others.

```{r message=FALSE, warning = FALSE}
library(tidyverse)
library(DBI)
library(xts)
library(dygraphs)
library(lubridate)
library(timeDate)
library(knitr)
library(reshape2)
library(grid)
library(gridExtra)
library(rgdal)
library(maptools)
library(raster)
library(readxl)


# A few constants
database = paste(getwd(),"/RosaliaESSD.db",sep="") # SQLite file to use
ph = 200 # height of precip plots
pw = 800 # width of precip plots
qh = 300 # height of Q plots
qw = 800 # width of Q plots
rswf = c(0, 0.3) # Plotting range for soil water content
```

# Dataset structure

The main part of the dataset is the time series database implemented in SQLite. Additionally, there is a folder which contains the shapefiles of the spatial datasets, and an MS Excel file wich contains the isotopic data.

For the time series, we are using a version of the database in file `r database` for this document. Please replace the name of your database file, if necessary. The database has a simple structure, consisting of 2 data tables only: The `catalog` is basically a table of contents, with column `sid` being the identifier of the time series. All time series observations are stored in table `timeseries` with columns `sid, times, val`. An auxiliary table `CatalogFields`describes the fields in the catalog. Read the code chunks below to learn how to access the database and perform simple SQL queries.

```{r tscatalog, message=FALSE, warning = FALSE}
con = dbConnect(RSQLite::SQLite(), dbname=database)
fullcatalog = dbGetQuery( con,'select * from catalog' )
knitr::kable(fullcatalog, caption = "Full catalog of available time series")
```

## Explanation of the catalog

The description of the fields in the `catalog` table is included in table `CatalogFields` and is defined as follows:
```{r message=FALSE, warning = FALSE}
ptimeseries = dbGetQuery(con, "select * from CatalogFields;")
kable(ptimeseries, caption = "Explanation of catalog fields")
```

## Timeseries table
The timeseries table contains all data (all sites and variables in 1 table). A specific timeseries is identified by its `sid` field. Examples are given below. Time in table timeseries is stored as `unixepoch` time (seconds since 1.January 1970 00:00:00 GMT). It is converted to an xts-compatible time by `dt1 = timeDate(times, format = "%Y-%m-%d %H:%M", FinCenter="UTC")`.



```{r message=FALSE, warning = FALSE}
ptimeseries = dbGetQuery(con, "select sid, times, val from timeseries limit 5;")
kable(ptimeseries, caption = "First 5 records of table timeseries")
```

# Discharge data

## Read the data

Reading the data from the database is explained in some detail related to the discharge data. These methods are applicable analogously for all the other time series as well.

Typical selections require 2 steps:

1. Select the timeseries from the `catalog`. The catalog can be printed as a table (here, by `kable`). Next, `cat1` is a query that retrieves a dataframe of the selected time series in the database. Adjust the SQL query to your requirements (in the following example, we select only discharge data with a quality flag of `ED1`, i.e., data have been cleaned from obvious errors and artefacts)! The result is a subset of the catalog (not yet the time series data!).


```{r message=FALSE, warning = FALSE}
sql1 = "select * from catalog where variable = 'DISCHARGE' and quality = 'ED1';"
cat1 = dbGetQuery( con, sql1)
kable(cat1, caption = paste("Selected time series by SQL statement:", sql1 ))
```

2. Select the observations for the selected `sid` in `cat1` from table `timeseries`. As a result, all the selected time series are put into a multivariate xts series `tsall`. In the example, the site name is used for the column names (assumes, for the plots with one common y-axis to be useful, that the same variable is selected for different sites). **Remark:** times of the data are recorded in UTC. To avoid possible confusion by "on-the-fly" conversions in some of the time series packages, UTC is explicitly specified wherever meaningful.

```{r message=FALSE, warning = FALSE}
tsall = xts(tzone="UTC")
for (sid in cat1$sid) {
  sel1 = paste("select times,val from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$times, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$val,order.by = dt1,tzone="UTC")
  tsall = merge(tsall, ts1,tzone="UTC")
  colnames(tsall)[ncol(tsall)] = cat1$site[cat1$sid == sid]
}

```
2. (optional): Select a time window only, e.g. after 2017-01-01

```{r message=FALSE, warning = FALSE}
tsallp = xts(tzone="UTC")
begin = as.numeric(as.POSIXct("2017-01-01",tz="UTC"))
for (sid in cat1$sid) {
  sel1 = paste("select times,val from timeseries where sid =",sid,"and times >=", begin,"; ")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$times, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$val,order.by = dt1,tzone="UTC")
  tsallp = merge(tsallp, ts1,tzone="UTC")
  colnames(tsallp)[ncol(tsallp)] = cat1$site[cat1$sid == sid]
}
dbDisconnect(con) # close database connection when all data are retrieved
```


## Aggregate time series for daily, weekly and monthly means

Here are a few examples how time series can be aggregated. We can use the `apply.monthly` and related functions in the `xts` package to do these operations efficiently on a whole `xts` object with a single function call. The results have timestamps at end-of-interval. Therefore, we use `floor_data` to shift the timestamps to begin of interval, to produce more intuitive graphs with `dygraph`.

```{r message=FALSE, warning = FALSE}
tsam = apply.monthly(tsall,mean)
time(tsam)= floor_date(time(tsam),unit="month") # shift time to begin of interval
tsaw = apply.weekly(tsall,mean)
time(tsaw)= floor_date(time(tsaw),unit="week")
tsad = apply.daily(tsall,mean)
time(tsad)= floor_date(time(tsad),unit="day")
```

## Plot graphs

`dygraph` is a very handy and powerful, interactive time series visualisation tool. It allows zooming and panning by mouse interaction, even in the HTML output. Multiple `dygraphs` can be put into a group and will zoom synchronously on the time axis. It is used in this document to provide an easy tool to explore the time series data.

### Plot original 10 min data with dygraph

```{r message=FALSE, warning = FALSE}
col1=c("black","red","green","purple")

dygraph(tsall, main = "Q hydrographs 10min", width = qw, height = qh, group ="Q") %>%
  dyAxis("y", label = "Q[l/s]") %>%
  dyRangeSelector() %>%
  dyOptions(colors=col1) %>%
dyHighlight(highlightCircleSize = 4, 
            highlightSeriesBackgroundAlpha = 0.4,
            hideOnMouseOut = TRUE)
```

**Remark:** On June 13th, 2018 at 16:00 there was an extreme peak of 582 l/s at Q3 due to a heavy storm event. During this event, large amounts of sediments were mobilised which have also blocked the culvert at Q2. Data are not reliable during this event.

### Plot original 10 min data in log scale

Because the watershed area of Q4 with 223 ha is `r round(223/9,1)` times larger than that of Q1 with 9 ha, a display of log scaled discharge supports easier comparison of the hydrographs.

```{r message=FALSE, warning = FALSE}
col1=c("black","red","green","purple")

dygraph(tsall, main = "Q hydrographs 10min - log scale", width = qw, height = qh, group ="Q") %>%
  dyAxis("y", label = "Q[l/s]") %>%
  dyRangeSelector() %>%
  dyOptions(colors=col1,logscale=TRUE) %>%
dyHighlight(highlightCircleSize = 4, 
            highlightSeriesBackgroundAlpha = 0.4,
            hideOnMouseOut = TRUE)
```

### Daily time series as step plot

```{r}
dygraph(tsad,main = "Q hydrographs daily",  width = qw, height = qh, group="Q")  %>%
  dyAxis("y", label = "Q[l/s]") %>%
  dyOptions(stepPlot=TRUE, colors=col1)
```

```{r}
dygraph(tsad,main = "Q hydrographs daily log scale ",  width = qw, height = qh, group="Q")  %>%
  dyAxis("y", label = "Q[l/s]") %>%
  dyOptions(stepPlot=TRUE, colors=col1, logscale = TRUE)
```

### Weekly time series as step plot
```{r message=FALSE, warning = FALSE}
dygraph(tsaw,main = "Q hydrographs weekly", width = qw, height = qh, group="Q")  %>%
  dyAxis("y", label = "Q[l/s]") %>%
  dyOptions(stepPlot=TRUE, colors=col1)

```

**Remark:** Interval means are not calculated if the time series contains at least one `NaN` in the interval. Data gaps therefore appear to be longer than they are.

# Precipitation data

## Read the data

As a result, all the selected time series are put into a multivariate xts series `psall`. Site name is set as column names (to get unique names, it is assumed that the same variable is selected for different sites). 

```{r message=FALSE, warning = FALSE}
con = dbConnect(RSQLite::SQLite(), dbname=database)
plist = c('K1','K2','K3', 'Q1', 'Q2', 'Q3', 'Q4')
cat1 = dbGetQuery( con,"select * from catalog where variable = 'PRECIPITATION' and quality = 'ED1' and (site in ('K1', 'K2', 'K3', 'Q1', 'Q2', 'Q3', 'Q4'));" )
kable(cat1, caption="Overview of precipitation timeseries")
psall = xts(tzone="UTC")
for (sid in cat1$sid) {
  sel1 = paste("select times,val from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$times, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$val,order.by = dt1,tzone="UTC")
  psall = merge(psall, ts1,tzone="UTC")
  colnames(psall)[ncol(psall)] = cat1$site[cat1$sid == sid]
}
dbDisconnect(con)
```

## Aggregate time series: yearly, monthly, daily sums

```{r message=FALSE, warning = FALSE}
options(knitr.kable.NA = '-')
psam = apply.monthly(psall,colSums)
time(psam)= floor_date(time(psam),unit="month") # shift to begin of interval
Month=as.character(as.yearmon(date(psam)))
psad = apply.daily(psall,colSums)
time(psad)= floor_date(time(psad),unit="day") # shift to begin of interval
psay = apply.yearly(psall,colSums)
Year = year(date(psay))
kable(cbind(Year,psay)[1:nrow(psay)-1,], caption="Mean annual precipitation (mm)")
```


```{r message=FALSE, warning = FALSE}
psamt=as.data.frame(psam)
kable(cbind(Month,psamt)[1:nrow(psamt)-1,],  row.names=FALSE,caption="Mean monthly precipitation (mm)")

```

**Remark:** Sites Q1 - Q4 use tipping bucket devices which are not working during the frost period and therefore do not have complete records for any year. K3 was installed in August 2018. 

### Plot daily time series as step plot

```{r message=FALSE, warning = FALSE, fig.height=2.5}
res = lapply(1:(length(plist)), function(i) dygraph(psad[,i],
    main = paste(plist[i], " daily precipitation"), 
    width = pw, height = ph, group="PD")  %>%
    dyAxis("y", valueRange = c(0, 60), label = "P[mm]") %>%
    dySeries(plist[i], stepPlot=TRUE, fillGraph = TRUE, color="blue"))
htmltools::tagList(res)
```


### Plot monthly time series as step plot

```{r message=FALSE, warning = FALSE, fig.height=2.5}
res = lapply(1:(length(plist)), function(i) dygraph(psam[,i],
    main = paste(plist[i], " monthly precipitation"), 
    width = pw, height = ph, group="PM")  %>%
    dyAxis("y", valueRange = c(0, 200), label = "P[mm]") %>%
    dySeries(plist[i], stepPlot=TRUE, fillGraph = TRUE, color="blue"))
htmltools::tagList(res)
```



**Remark:** Interval sums are not calculated if the time series contains at least one `NaN` in the interval. Data gaps therefore appear to be longer than they are.

# Soil water profiles

Hydraprobe II sensors are installed at depths 10, 20, 40, (60) cm. Sensor addresses are 4, 5, 6, (7). (I.e., Q2S0.04 refers to the sensor at depth 10 cm, site Q2S0, etc.). Soil moisture is recorded as "soil moisture water fraction" (swf) and soil temperature in °C.

## Soil moisture

```{r message=FALSE, warning = FALSE}
depth4 = c(0.1,0.2,0.4,0.6)
depth3 = c(0.1,0.2,0.4)
con = dbConnect(RSQLite::SQLite(), dbname=database)
swlistq1s0 = c('Q1S0.04','Q1S0.05','Q1S0.06','Q1S0.07')
swlistq2s0 = c('Q2S0.04','Q2S0.05','Q2S0.06','Q2S0.07')
swlistq2s1 = c('Q2S1.04','Q2S1.05','Q2S1.06')
swlistq2s2 = c('Q2S2.04','Q2S2.05','Q2S2.06')
selq1s0 = "select * from catalog where variable = 'SOIL MOISTURE WATER FRACTION' and quality = 'ED1' and (site in ('Q1S0.04', 'Q1S0.05', 'Q1S0.06', 'Q1S0.07'));"
selq2s0 = "select * from catalog where variable = 'SOIL MOISTURE WATER FRACTION' and quality = 'ED1' and (site in ('Q2S0.04', 'Q2S0.05', 'Q2S0.06', 'Q2S0.07'));"
selq2s1 = "select * from catalog where variable = 'SOIL MOISTURE WATER FRACTION' and quality = 'ED1' and (site in ('Q2S1.04', 'Q2S1.05', 'Q2S1.06'));"
selq2s2 = "select * from catalog where variable = 'SOIL MOISTURE WATER FRACTION' and quality = 'ED1' and (site in ('Q2S2.04', 'Q2S2.05', 'Q2S2.06'));"
catq1s0 = dbGetQuery( con, selq1s0)
catq1s0$depth = depth4
kable(catq1s0,caption="Soil moisture series at site Q2S0")
```


```{r message=FALSE, warning = FALSE}
swallq1s0 = xts(tzone="UTC")
for (sid in catq1s0$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  swallq1s0 = merge(swallq1s0, ts1,tzone="UTC")
  colnames(swallq1s0)[ncol(swallq1s0)] = catq1s0$site[catq1s0$sid == sid]
}
catq2s0 = dbGetQuery( con, selq2s0)
catq2s0$depth = depth4
kable(catq2s0,caption="Soil moisture series at site Q2S0")
```


```{r message=FALSE, warning = FALSE}
swallq2s0 = xts(tzone="UTC")
for (sid in catq2s0$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  swallq2s0 = merge(swallq2s0, ts1,tzone="UTC")
  colnames(swallq2s0)[ncol(swallq2s0)] = catq2s0$site[catq2s0$sid == sid]
}
catq2s1 = dbGetQuery( con, selq2s1)
catq2s1$depth = depth3
kable(catq2s1,caption="Soil moisture series at site Q2S1")
```


```{r message=FALSE, warning = FALSE}
swallq2s1 = xts(tzone="UTC")
for (sid in catq2s1$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  swallq2s1 = merge(swallq2s1, ts1,tzone="UTC")
  colnames(swallq2s1)[ncol(swallq2s1)] = catq2s1$site[catq2s1$sid == sid]
}
catq2s2 = dbGetQuery( con, selq2s2)
catq2s2$depth = depth3
kable(catq2s2,caption="Soil moisture series at site Q2S2")
swallq2s2 = xts(tzone="UTC")
for (sid in catq2s2$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  swallq2s2 = merge(swallq2s2, ts1,tzone="UTC")
  colnames(swallq2s2)[ncol(swallq2s2)] = catq2s2$site[catq2s2$sid == sid]
}
dbDisconnect(con)
```
## Plot full resolution data (example: Q2S0)

Soil moisture data are quite voluminous if plotted in full 10-min resolution. Only data for Q2S0 are displayed here as an example. Other sites can be accessed analogously.

```{r message=FALSE, warning=FALSE}
#' Create plots.
invisible(Sys.setlocale("LC_TIME", "English"))
col1 = c("black", "red", "green", "purple")

dygraph(swallq2s0,main = "10-min soil moisture at site Q2S0", width = qw, height = qh, group="SWF")  %>%
  dySeries("Q2S0.04", label = "10 cm") %>%
  dySeries("Q2S0.05", label = "20 cm") %>%
  dySeries("Q2S0.06", label = "40 cm") %>%
  dySeries("Q2S0.07", label = "60 cm") %>%
  dyAxis("y", label = "swf[-]") %>%
  dyOptions(stepPlot=TRUE, colors=col1)
```


## Aggregate time series for daily means

```{r message=FALSE, warning = FALSE}
swadq1s0 = apply.daily(swallq1s0,mean)
time(swadq1s0)= floor_date(time(swadq1s0),unit="day") # shift to begin of interval
swadq2s0 = apply.daily(swallq2s0,mean)
time(swadq2s0)= floor_date(time(swadq2s0),unit="day") # shift to begin of interval
swadq2s1 = apply.daily(swallq2s1,mean)
time(swadq2s1)= floor_date(time(swadq2s1),unit="day") # shift to begin of interval
swadq2s2 = apply.daily(swallq2s2,mean)
time(swadq2s2)= floor_date(time(swadq2s2),unit="day") # shift to begin of interval
```

### Plot daily data

```{r message=FALSE, warning = FALSE}

#' Create plots.
invisible(Sys.setlocale("LC_TIME", "English"))
col1 = c("black", "red", "green", "purple")
res = list(
dygraph(swadq1s0,main = "Daily soil moisture at site Q1S0", width = qw, height = qh, group="SWF")  %>%
  dySeries("Q1S0.04", label = "10 cm") %>%
  dySeries("Q1S0.05", label = "20 cm") %>%
  dySeries("Q1S0.06", label = "40 cm") %>%
  dySeries("Q1S0.07", label = "60 cm") %>%
  dyAxis("y", label = "swf[-]") %>%
  dyOptions(stepPlot=TRUE, colors=col1),
dygraph(swadq2s0,main = "Daily soil moisture at site Q2S0", width = qw, height = qh, group="SWF")  %>%
  dySeries("Q2S0.04", label = "10 cm") %>%
  dySeries("Q2S0.05", label = "20 cm") %>%
  dySeries("Q2S0.06", label = "40 cm") %>%
  dySeries("Q2S0.07", label = "60 cm") %>%
  dyAxis("y", label = "swf[-]") %>%
  dyOptions(stepPlot=TRUE, colors=col1),
dygraph(swadq2s1,main = "Daily soil moisture at site Q2S1", width = qw, height = qh, group="SWF")  %>%
  dySeries("Q2S1.04", label = "10 cm") %>%
  dySeries("Q2S1.05", label = "20 cm") %>%
  dySeries("Q2S1.06", label = "40 cm") %>%
  dyAxis("y", label = "swf[-]") %>%
  dyOptions(stepPlot=TRUE, colors=col1),
dygraph(swadq2s2,main = "Daily soil moisture at site Q2S2", width = qw, height = qh, group="SWF")  %>%
  dySeries("Q2S2.04", label = "10 cm") %>%
  dySeries("Q2S2.05", label = "20 cm") %>%
  dySeries("Q2S2.06", label = "40 cm") %>%
  dyRangeSelector() %>%
  dyAxis("y", label = "swf[-]") %>%
  dyOptions(stepPlot=TRUE, colors=col1))
htmltools::tagList(res)
```

## Soil temperature

```{r message=FALSE, warning = FALSE}
depth4 = c(0.1,0.2,0.4,0.6)
depth3 = c(0.1,0.2,0.4)
con = dbConnect(RSQLite::SQLite(), dbname=database)
swlistq1s0 = c('Q1S0.04','Q1S0.05','Q1S0.06','Q1S0.07')
swlistq2s0 = c('Q2S0.04','Q2S0.05','Q2S0.06','Q2S0.07')
swlistq2s1 = c('Q2S1.04','Q2S1.05','Q2S1.06')
swlistq2s2 = c('Q2S2.04','Q2S2.05','Q2S2.06')
selq1s0 = "select * from catalog where variable = 'SOIL TEMPERATURE' and quality = 'ED1' and (site in ('Q1S0.04', 'Q1S0.05', 'Q1S0.06', 'Q1S0.07'));"
selq2s0 = "select * from catalog where variable = 'SOIL TEMPERATURE' and quality = 'ED1' and (site in ('Q2S0.04', 'Q2S0.05', 'Q2S0.06', 'Q2S0.07'));"
selq2s1 = "select * from catalog where variable = 'SOIL TEMPERATURE' and quality = 'ED1' and (site in ('Q2S1.04', 'Q2S1.05', 'Q2S1.06'));"
selq2s2 = "select * from catalog where variable = 'SOIL TEMPERATURE' and quality = 'ED1' and (site in ('Q2S2.04', 'Q2S2.05', 'Q2S2.06'));"
catq1s0 = dbGetQuery( con, selq1s0)
catq1s0$depth = depth4
kable(catq1s0,caption="Soil temperature series at site Q1S0")
```


```{r message=FALSE, warning = FALSE}
stallq1s0 = xts(tzone="UTC")
for (sid in catq1s0$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  stallq1s0 = merge(stallq1s0, ts1,tzone="UTC")
  colnames(stallq1s0)[ncol(stallq1s0)] = catq1s0$site[catq1s0$sid == sid]
}
catq2s0 = dbGetQuery( con, selq2s0)
catq2s0$depth = depth4
kable(catq2s0,caption="Soil temperature series at site Q2S0")
```


```{r message=FALSE, warning = FALSE}
stallq2s0 = xts(tzone="UTC")
for (sid in catq2s0$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  stallq2s0 = merge(stallq2s0, ts1,tzone="UTC")
  colnames(stallq2s0)[ncol(stallq2s0)] = catq2s0$site[catq2s0$sid == sid]
}
catq2s1 = dbGetQuery( con, selq2s1)
catq2s1$depth = depth3
kable(catq2s1,caption="Soil temperature series at site Q2S1")
```


```{r message=FALSE, warning = FALSE}
stallq2s1 = xts(tzone="UTC")
for (sid in catq2s1$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  stallq2s1 = merge(stallq2s1, ts1,tzone="UTC")
  colnames(stallq2s1)[ncol(stallq2s1)] = catq2s1$site[catq2s1$sid == sid]
}
catq2s2 = dbGetQuery( con, selq2s2)
catq2s2$depth = depth3
kable(catq2s2,caption="Soil temperature series at site Q2S2")
stallq2s2 = xts(tzone="UTC")
for (sid in catq2s2$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  stallq2s2 = merge(stallq2s2, ts1,tzone="UTC")
  colnames(stallq2s2)[ncol(stallq2s2)] = catq2s2$site[catq2s2$sid == sid]
}
dbDisconnect(con)
```
## Aggregate time series for daily means

```{r message=FALSE, warning = FALSE}
stadq1s0 = apply.daily(stallq1s0,mean)
time(stadq1s0)= floor_date(time(stadq1s0),unit="day") # shift to begin of interval
stadq2s0 = apply.daily(stallq2s0,mean)
time(stadq2s0)= floor_date(time(stadq2s0),unit="day") # shift to begin of interval
stadq2s1 = apply.daily(stallq2s1,mean)
time(stadq2s1)= floor_date(time(stadq2s1),unit="day") # shift to begin of interval
stadq2s2 = apply.daily(stallq2s2,mean)
time(stadq2s2)= floor_date(time(stadq2s2),unit="day") # shift to begin of interval
```

### Plot daily data

```{r message=FALSE, warning = FALSE}

#' Create plots.
invisible(Sys.setlocale("LC_TIME", "English"))
col1 = c("black", "red", "green", "purple")
res = list(
dygraph(stadq1s0,main = "Daily soil temperature at site Q1S0", width = qw, height = qh, group="T")  %>%
  dySeries("Q1S0.04", label = "10 cm") %>%
  dySeries("Q1S0.05", label = "20 cm") %>%
  dySeries("Q1S0.06", label = "40 cm") %>%
  dySeries("Q1S0.07", label = "60 cm") %>%
  dyAxis("y", label = "Soil temperature [°C]") %>%
  dyOptions(stepPlot=TRUE, colors=col1),
dygraph(stadq2s0,main = "Daily soil temperature at site Q2S0", width = qw, height = qh, group="T")  %>%
  dySeries("Q2S0.04", label = "10 cm") %>%
  dySeries("Q2S0.05", label = "20 cm") %>%
  dySeries("Q2S0.06", label = "40 cm") %>%
  dySeries("Q2S0.07", label = "60 cm") %>%
  dyAxis("y", label = "Soil temperature [°C]") %>%
  dyOptions(stepPlot=TRUE, colors=col1),
dygraph(stadq2s1,main = "Daily soil temperature at site Q2S1", width = qw, height = qh, group="T")  %>%
  dySeries("Q2S1.04", label = "10 cm") %>%
  dySeries("Q2S1.05", label = "20 cm") %>%
  dySeries("Q2S1.06", label = "40 cm") %>%
  dyAxis("y", label = "Soil temperature [°C]") %>%
  dyOptions(stepPlot=TRUE, colors=col1),
dygraph(stadq2s2,main = "Daily soil temperature at site Q2S2", width = qw, height = qh, group="T")  %>%
  dySeries("Q2S2.04", label = "10 cm") %>%
  dySeries("Q2S2.05", label = "20 cm") %>%
  dySeries("Q2S2.06", label = "40 cm") %>%
  dyRangeSelector() %>%
  dyAxis("y", label = "Soil temperature [°C]") %>%
  dyOptions(stepPlot=TRUE, colors=col1))
htmltools::tagList(res)
```

# Air temperature, relative humidity and water temperature

Air temperature and relative humidity are recorded at sites K1, K2, Q1 - Q4 using an Adcon TR1 sensor. Water temperature is measured by the conductivity sensors C4E at sites Q1 and Q2, and by the Condu:lyte sensor at Q4.

## Air temperature

```{r}
con = dbConnect(RSQLite::SQLite(), dbname=database)
selatemp = "select * from catalog where variable = 'TEMPERATURE' and quality = 'ED1' ;"
catatemp = dbGetQuery( con, selatemp)
kable(catatemp,caption="Air temperature series")
tempall = xts(tzone="UTC")
col6 = c("black", "red", "green", "purple","orange","cyan")
for (sid in catatemp$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  tempall = merge(tempall, ts1,tzone="UTC")
  colnames(tempall)[ncol(tempall)] = catatemp$site[catatemp$sid == sid]
}
# Aggregate hourly to reduce graph complexity
ep = endpoints(tempall, "hours")
temphr = period.apply(tempall,ep,mean)
time(temphr)= floor_date(time(temphr),unit="hour") # shift to begin of interval

dygraph(temphr,main = "Air temperature (hourly)", width = qw, height = qh, group="T")  %>%
  dyRangeSelector() %>%
  dyAxis("y", label = "Air Temperature[°C]") %>%
  dyOptions(colors=col6)

dbDisconnect(con)
```

## Relative humidity

```{r}
con = dbConnect(RSQLite::SQLite(), dbname=database)
selrh = "select * from catalog where variable = 'RELATIVE HUMIDITY' and quality = 'ED1' ;"
catrh = dbGetQuery( con, selrh)
kable(catrh,caption="Relative Humidity series")
rhall = xts(tzone="UTC")
col6 = c("black", "red", "green", "purple","orange","cyan")
for (sid in catrh$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  rhall = merge(rhall, ts1,tzone="UTC")
  colnames(rhall)[ncol(rhall)] = catrh$site[catrh$sid == sid]
}
# Aggregate hourly to reduce graph complexity
ep = endpoints(rhall, "hours")
rhhr = period.apply(rhall,ep,mean)
time(rhhr)= floor_date(time(rhhr),unit="hour") # shift to begin of interval
dygraph(rhhr,main = "Relative humidity", width = qw, height = qh, group="T")  %>%
  dyRangeSelector() %>%
  dyAxis("y", label = "rH [%]") %>%
  dyOptions(colors=col6)

dbDisconnect(con)
```

## Conductivity

```{r}
con = dbConnect(RSQLite::SQLite(), dbname=database)
selcon = "select * from catalog where variable = 'CONDUCTIVITY' and quality = 'ED1' ;"
catcon = dbGetQuery( con, selcon)
kable(catcon,caption="Conductivity series")
conall = xts(tzone="UTC")
col6 = c("black", "red", "green", "purple","orange","cyan")
for (sid in catcon$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  conall = merge(conall, ts1,tzone="UTC")
  colnames(conall)[ncol(conall)] = catcon$site[catcon$sid == sid]
}
dygraph(conall,main = "Conductivity", width = qw, height = qh, group="T")  %>%
  dyRangeSelector() %>%
  dyAxis("y", label = "Conductivity [µS/cm]") %>%
  dyOptions(colors=col6)

dbDisconnect(con)
```
Conductivity data should be used with caution. Unfortunately, the devices installed at Q1 and Q2 (Ponsel C4E) tended to be rather sensitive to small surface contamination, leading to significant offsets of the recorded conductivity. In both creeks, the baseflow conductivity is quite close to 120 µS/cm. However, the relative changes of conductivity during rainfall-runoff events are very plausible and the data are trusted to be still useful. At site Q4, a different type of instrument (s::can Condulyser) is installed, which up to now delivered reliable data. The considerably higher level of baseflow conductivity of about 200 µS/cm has been verified and can be explained by the much larger watershed (140 ha, versus 27 ha) as well as by changes in the geological conditions.

## Water temperature

```{r message=TRUE, warning=TRUE}
con = dbConnect(RSQLite::SQLite(), dbname=database)
selwt = "select * from catalog where variable = 'WATER TEMPERATURE' and quality = 'ED1' ;"
catwt = dbGetQuery( con, selwt)
kable(catwt,caption="Water temperature series")
wtall = xts(tzone="UTC")
col6 = c("black", "red", "green", "purple","orange","cyan")
for (sid in catwt$sid) {
  sel1 = paste("select times as Zeit,val as Wert from timeseries where sid =",sid,";")
  ts1 = dbGetQuery(con,sel1)
  dt1 = timeDate(ts1$Zeit, format = "%Y-%m-%d %H:%M", FinCenter="UTC")
  ts1 = as.xts(ts1$Wert,order.by = dt1,tzone="UTC")
  wtall = merge(wtall, ts1,tzone="UTC")
  colnames(wtall)[ncol(wtall)] = catwt$site[catwt$sid == sid]
}
# Aggregate hourly to reduce graph complexity
ep = endpoints(wtall, "hours")
wthr = period.apply(wtall,ep,mean)
time(wthr)= floor_date(time(wthr),unit="hour") # shift to begin of interval

dygraph(wthr,main = "Water temperature (hourly)", width = qw, height = qh,group="T")  %>%
  dyRangeSelector() %>%
  dyAxis("y", label = "Temperature [°C]") %>%
  dyOptions(colors=col6)

dbDisconnect(con)
```
# Isotopic data

The stable isotopes of water as an environmental tracer have found versatile use in catchment science in, e.g., estimating water transit times or hydrograph separation studies. They are a powerful tool, complementing the hydrometeorological network of the Rosalia. However, data acquisition was started as late as in 2019. At discharge site Q4 river and precipitation samples are collected since June and October 2019. To date, the precipitation data have been collected as bi-weekly bulk samples and are compared to the daily river water grab samples. Comparison of precipitation with river isotope values already shows the response of the discharge to the precipitation input tracer signal.

Due to the shorter time series and the required laboratory analysis for every single sample, the data are provided in the format of an MS Excel sheet. Isotopic data are stored in a separate sheet for each sampling point and series. Each series has a site ID as described below.

```{r}
isoxlsx = paste(getwd(),"/Isotope_ESSD.xlsx",sep="")
Explanation = read_xlsx(isoxlsx,sheet="Explanation",skip=2)
kable(Explanation,caption="Description of isotopic data file content")
K1PIso = read_xlsx(isoxlsx,sheet="K1PIso")
Q4PIso = read_xlsx(isoxlsx,sheet="Q4PIso")
Q1GrabIso = read_xlsx(isoxlsx,sheet="Q1GrabIso")
Q2GrabIso = read_xlsx(isoxlsx,sheet="Q2GrabIso")
Q4DailyIso = read_xlsx(isoxlsx,sheet="Q4DailyIso")

tsiso = xts(tzone="UTC")
ts1 = as.xts(Q4DailyIso$d18O, order.by= Q4DailyIso$`Date (UTC)`,tzone="UTC")
tsiso = merge(tsiso, ts1,tzone="UTC")
colnames(tsiso)[ncol(tsiso)] = "d18O"
ts1 = as.xts(Q4DailyIso$d2H, order.by= Q4DailyIso$`Date (UTC)`,tzone="UTC")
#tsiso = merge(tsiso, ts1,tzone="UTC") # not displayed in plot
#colnames(tsiso)[ncol(tsiso)] = "d2H"
ts1 = as.xts(Q4DailyIso$`d-excess`, order.by= Q4DailyIso$`Date (UTC)`,tzone="UTC")
#tsiso = merge(tsiso, ts1,tzone="UTC") # not displayed in plot
#colnames(tsiso)[ncol(tsiso)] = "d-excess"
ts1 = as.xts(K1PIso$d18O, order.by= K1PIso$Datefrom,tzone="UTC")
tsiso = merge(tsiso, ts1,tzone="UTC")
colnames(tsiso)[ncol(tsiso)] = 'd18O_P_K1'
ts1 = as.xts(Q4PIso$d18O, order.by= Q4PIso$Datefrom,tzone="UTC")
tsiso = merge(tsiso, ts1,tzone="UTC")
colnames(tsiso)[ncol(tsiso)] = "d18O_P_Q4"
dygraph(tsiso,main = "Isotopic data Q4", width = qw, height = qh, group="I")  %>%
  dyRangeSelector() %>%
  dyAxis("y", label = "River d18O[‰]", valueRange = c(-13,-7)) %>%
  dyAxis("y2",label = "Precip d18O[‰]", valueRange = c(-16,-4), independentTicks = T)%>%
  dySeries("d18O_P_Q4",stepPlot = T, axis="y2") %>%
  dySeries("d18O_P_K1",stepPlot = T, axis="y2") %>%
  dyOptions(connectSeparatedPoints = TRUE)
  
```


# GIS datasets

To facilitate working with the data, a few essential GIS datasets are provided: Watershed boundaries for the 4 stream gauges Q1 - Q4, a shapefile of the monitoring sites, the streams (result of ground survey) and a 10x10 m DEM.

```{r fig.height=10, fig.width=10, message=F, warning=F}
#ogrDrivers()
dsn = paste(getwd(),"/gis",sep="")
#ogrListLayers(dsn)
ws14 <- readOGR( 
  dsn= dsn, 
  layer="WatershedsQ14",  # Watersheds of gauges Q1 - Q4
  verbose=FALSE
)
gew <- readOGR( 
  dsn= dsn, 
  layer="Streams2018", # Streams
  verbose=FALSE
)
ms <- readOGR( 
  dsn= dsn, 
  layer="SitesESSD",  # Monitoring sites
  verbose=FALSE
)
dem = raster(paste(dsn,"/RosaliaDEM10m.tif",sep=""))
plot(dem, ext=extent(-4200,-2000, 284500,286800), axes=T, box=F, main="Map of sites - Overview",
     legend.args = list(text = 'Height (m a.s.l.)'))
plot(ws14, add=TRUE)
gewlwd = c(1,1.5,2,2.5)[gew$Breite]
plot(gew, col = "blue", lwd = gewlwd, add=TRUE)
plot(ms, pch=ms$Symbol, add=TRUE, cex = 1.5)
pointLabel(coordinates(ms),labels=ms$Site)
legend("bottomright",
       legend = c("Stream gauge","Rain gauge","Soil water","Relais", "Base station"),
       pch=c(17,18,16,5,15),
       bty = "n", # turn off border
       cex = 1) # adjust legend font size

#Zoom in detail around Q2
plot(dem, ext=extent(-2800,-2650, 284950,285100), axes=T, box=F, main="Map of sites - Detail Q2",
     legend.args = list(text = 'Height (m a.s.l.)'))
plot(ws14, add=TRUE)
gewlwd = c(1,1.5,2,2.5)[gew$Breite]
plot(gew, col = "blue", lwd = gewlwd, add=TRUE)

plot(ms, pch=ms$Symbol, add=TRUE, cex = 1.5)


pointLabel(coordinates(ms),labels=ms$Site)
legend("bottomright",
       legend = c("Stream gauge","Rain gauge","Soil water","Relais", "Base station"),
       pch=c(17,18,16,5,15),
       bty = "n", # turn off border
       cex = 1) # adjust legend font size


```


