README.Rmd

---
output:
  md_document:
    variant: markdown_github
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, echo = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "README-"
)
```

# openeddy

This package provides utilities for eddy covariance data handling, 
quality checking (similar to Mauder et al., 2013), processing, summarising 
and plotting. It aims to standardise the automated quality checking and make 
data processing reproducible.

## Installation

1. Install devtools package if not available yet.
```r
install.packages("devtools")
```
1. Install openeddy
```r
devtools::install_github("lsigut/openeddy")
```

## Extended Example

An extended example describing the intended eddy covariance data processing
workflow is available at:

https://github.com/lsigut/EC_workflow

## Short Example

Demonstration of selected generally applicable `openeddy` functions.

```{r}
library(openeddy)
library(REddyProc)
library(ggplot2)
```

Example data from `REddyProc` package do not include statistics extracted from
raw data or quality control (QC) flags. Notice that units are included. 

```{r}
data(Example_DETha98)
str(Example_DETha98[1:4])
DETha98 <- fConvertTimeToPosix(Example_DETha98, 'YDH', 
                               Year = 'Year', Day = 'DoY', Hour = 'Hour')
```

Certain variable names are expected by `openeddy`, thus DETha98 variable names need to be renamed. Net ecosystem exchange (NEE) is already filtered based on
QC information (qc_NEE) that is not included. Though flag 2 values can be recreated, respective NEE values are missing and flag 1 cannot be resolved.

```{r}
rename <- names(DETha98) %in% c("DateTime", "Rg")
names(DETha98)[rename] <- c("timestamp", "GR")
DETha98$qc_NEE <- ifelse(is.na(DETha98$NEE), NA, 0)
```

Application of three general filters is presented. 

1) Low covariance between vertical wind component and CO~2~ concentration 
((CO~2~)) can be caused by frozen ultrasonic anemometer or problems with (CO~2~) 
measurements. Such cases are flagged and saved to qc_NEE_lowcov. 
1) Runs with repeating values are a sign of malfunctioning equipment 
(qc_NEE_runs).
1) Spikes in low frequency data cause problems during gap-filling and should be
excluded. Since DETha98 was already quality checked, the amount of detected
spikes is limited. In order to correctly evaluate spikes, preliminary QC 
(qc_NEE_prelim) that combines available QC tests or filters should be produced
and used in `despikeLF`.

```{r}
DETha98$qc_NEE_lowcov <- 
  apply_thr(DETha98$NEE, c(-0.01, 0.01), "qc_NEE_lowcov", "between")
table(DETha98$qc_NEE_lowcov)
DETha98$qc_NEE_runs <- flag_runs(DETha98$NEE, "qc_NEE_runs")
table(DETha98$qc_NEE_runs)
DETha98$qc_NEE_prelim <- 
  combn_QC(DETha98, 
           c("qc_NEE", "qc_NEE_lowcov", "qc_NEE_runs"), 
           "qc_NEE_prelim", additive = FALSE, na.as = NA)
DETha98$qc_NEE_despikeLF <- 
  despikeLF(DETha98, "NEE", "qc_NEE_prelim", "qc_NEE_despikeLF", 
            light = NULL)
table(DETha98$qc_NEE_despikeLF)
```

The QC results can be summarized in tabular or graphical form using 
`summary_QC`. It is possible to summarize each filter independently or 
summarize the cummulative effect of applied filters. Note that the fraction of
flagged records in this example is negligible as DETha98 dataset was already
quality checked.

```{r}
summary_QC(DETha98, 
           c("qc_NEE", "qc_NEE_lowcov", "qc_NEE_runs", "qc_NEE_despikeLF"),
           na.as = c(NA, NA, NA, 0))
summary_QC(DETha98, 
           c("qc_NEE", "qc_NEE_lowcov", "qc_NEE_runs", "qc_NEE_despikeLF"),
           na.as = c(NA, NA, NA, 0), cumul = TRUE, plot = TRUE, flux = "NEE")
```

Although individual QC columns should be stored as they are useful to 
distinguish the reason why certain records were excluded, only the combined QC
column (qc_NEE_composite) is usually used in further data processing and 
analysis.

```{r}
DETha98$qc_NEE_composite <- 
  combn_QC(DETha98, 
           c("qc_NEE", "qc_NEE_lowcov", "qc_NEE_runs", "qc_NEE_despikeLF"), 
           "qc_NEE_composite", additive = FALSE, na.as = c(NA, NA, NA, 0))
```

Function `plot_eddy` is useful for visualization of the whole dataset including
flux values, its respective QC flags and the most important micrometeorological
parameters in monthly and weekly time resolution. Only a two week subset is 
presented here to limit the extent of output.

```{r}
DETha98[, c("P", "PAR", "Rn")] <- NA
(varnames <- varnames(DETha98))
(units <- openeddy::units(DETha98))
sub <- DETha98$DoY >= 29 & DETha98$DoY < 43
DETha98_sub <- DETha98[sub, ]
openeddy::units(DETha98) <- units
plot_eddy(DETha98_sub, "NEE", "qc_NEE", "qc_NEE_composite", skip = "monthly",
          light = "GR")
```

In addition to actual despiking, `despikeLF` can be used also for visualization
of the internally computed double-differenced time series in order to inspect
selected 13 days blocks. See section Plotting in `despikeLF` help file for
further description.

```{r}
despikeLF_plots <- 
  despikeLF(DETha98, "NEE", "qc_NEE_prelim", "qc_NEE_despikeLF", 
            light = NULL, plot = TRUE)$plots
despikeLF_plots$`iter 1`$all$`1998-01-27 - 1998-02-08`
```

## References

Publication describing openeddy is not yet available. When describing 
the proposed quality control scheme, please refer to it as similar to:

Mauder, M., Cuntz, M., Drüe, C., Graf, A., Rebmann, C., Schmid, H.P., 
Schmidt, M., Steinbrecher, R., 2013. A strategy for quality and uncertainty 
assessment of long-term eddy-covariance measurements. Agric. For. Meteorol. 
169, 122-135, https://doi.org/10.1016/j.agrformet.2012.09.006