---
author: "Thomas Wutzler"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    keep_md: true
vignette: >
  %\VignetteEngine{knitr::knitr}
  %\VignetteIndexEntry{DEGeb example: UStar estimation}
  %\usepackage[UTF-8]{inputenc}
---

```{r, include = FALSE}
# do not execute on CRAN: 
# https://stackoverflow.com/questions/28961431/computationally-heavy-r-vignettes
is_check <- ("CheckExEnv" %in% search()) || any(c("_R_CHECK_TIMINGS_",
             "_R_CHECK_LICENSE_") %in% names(Sys.getenv()))
knitr::opts_chunk$set(eval = !is_check)
```

```{r setup, include = FALSE}
library(knitr)
# rmarkdown::render("vignettes/DEGebExample.Rmd")
opts_chunk$set(out.extra = 'style = "display:block; margin: auto"'
    #, fig.align = "center"
    , fig.width = 4.3, fig.height = 3.2, dev.args = list(pointsize = 10)
    , message = FALSE
    )
knit_hooks$set(spar = function(before, options, envir) {
    if (before) {
        par( las = 1 )                   #also y axis labels horizontal
        par(mar = c(2.0,3.3,0,0) + 0.3 )  #margins
        par(tck = 0.02 )                          #axe-tick length inside plots             
        par(mgp = c(1.1,0.2,0) )  #positioning of axis title, axis labels, axis
     }
})
# genVigs("DEGebExample")
```

Crop example demonstrating multiple years and user defined uStar-Seasons
========================================================================

```{r results = 'hide'}
#isDevelopMode <- TRUE
if (!exists("isDevelopMode")) library(REddyProc)
set.seed(0815)      # for reproducible results
```

First, the data is loaded. This example uses data that has been downloaded 
from http://www.europe-fluxdata.eu
and preprocessed by `fLoadEuroFlux16`, where the DateTime Column has been created, 
and the variables renamed to the BGC-convention (e.g. Tair instead of Ta).

```{r}
data(DEGebExample)
summary(DEGebExample)
```


VPD was not given with the original dataset and is calculated from Tair and rH.
```{r}
DEGebExample$VPD <- fCalcVPDfromRHandTair(DEGebExample$rH, DEGebExample$Tair)
EProc <- sEddyProc$new('DE-Geb', DEGebExample, c('NEE','Rg','Tair','VPD', 'Ustar'))
EProc$sSetLocationInfo(LatDeg = 51.1, LongDeg = 10.9, TimeZoneHour = 1)  #Location of Gebesee
```

Defining Seasons with different surface friction conditions 
-----------------------------------------------------------

The site is a crop site. The harvesting times are visible as sharp edges in the plots of NEE.

The micrometeorological conditions differ between the different cropping periods,
because the friction at the surface differs.
Also not the cropping periods do not correspond very well to seasons.
Hence, for the estimation of uStar-Thresholds, we apply a user-defined splitting 
of uStar-seasons. With function `usCreateSeasonFactorYdayYear` we provide the starting
points of the seasons.  

Note that, here, the seasons are not constrained within one calendaryear. 
There are other variants of a user-specified season that do respect calendaryear boundaries, 
or that let seasons start at the same day each year. 

```{r DEGeb_estUStar1a, spar = TRUE, fig.width = 10}
seasonStarts <- as.data.frame( do.call( rbind, list(
	  c(70,2004)
          ,c(210,2004)
          ,c(320,2004)
          ,c(70,2005)
          ,c(180,2005)
          ,c(320,2005)
          ,c(120,2006)
          ,c(305,2006) 		
)))
seasonFactor <- usCreateSeasonFactorYdayYear(
  DEGebExample$DateTime - 15*60, starts = seasonStarts)
plot( NEE ~ DateTime, DEGebExample )
seasonStartsDate <- fConvertTimeToPosix( data.frame(Year = seasonStarts[,2]
	, DoY = seasonStarts[,1], Hour = 0.25), 'YDH'
	, Year = "Year", Day = "DoY", Hour = "Hour")
abline( v = seasonStartsDate$DateTime)
```

The user-specific seasoning is provided to the gap-filling by the argument `seasonFactor`. 
```{r DEGeb_estUStar1b, cache = TRUE}
(uStarTh <- EProc$sEstUstarThold(seasonFactor = seasonFactor))
#EProc$useSeaonsalUStarThresholds()
# estimation can be inspected by plotting the saturation of NEE with UStar 
# for the temperatures of one season
#EProc$sPlotNEEVersusUStarForSeason( levels(uStarTh$season)[2] )
```

Note that there is an estimate for each season. Further, an annual estimate is 
obtained by taking the maximum across the seasons, and the overall estimate is 
the mean across the years.

By default the gap-filling uses annually aggregated estimates of uStar-Threshold.
This usually works for sites with continuous vegetation cover.
For the crop-site of this example, we will use a different threshold for each of 
the defined seasons, by calling `sEddyProc_useSeaonsalUStarThresholds` 
before gapfilling.

```{r DEGeb_estUStar1c}
EProc$useSeaonsalUStarThresholds()
EProc$sGetUstarScenarios()
```	

```{r DEGeb_gapFill1, cache = TRUE, eval = FALSE, output = 'hide', message = FALSE}
EProc$sMDSGapFillAfterUstar('NEE', FillAll = FALSE, isVerbose = FALSE)
```		


Uncertainty introduced by the uStar Threshold estimate: bootstrap  
-----------------------------------------------------------------

With a lower estimate of uStar threshold, more records with lower NEE are kept in 
the dataset instead of marked as gaps. Therefore annual estimate of NEE will decrease
with lower uStar Threshold. Also the partitioning of the net-flux to GPP and Reco is 
sensitive to inclusion of data at dawn period with conditions of low uStar.

In order to quantify this uncertainty, a lower, median and upper estimates of uStar 
are obtained from a bootstrapped sample of half-hourly NEE measurements.
The Gap-Filling and computation of derived quantities such as GPP are then repeated for
different estimates of the uStar Threshold.


```{r DEGeb_estUStarBoot1, cache = TRUE, message = FALSE}
EProc <- sEddyProc$new('DE-Geb', DEGebExample, c('NEE','Rg','Tair','VPD', 'Ustar'))
EProc$sSetLocationInfo(LatDeg = 51.1, LongDeg = 10.9, TimeZoneHour = 1)  #Location of Gebesee
# here, because of processing time only 30 samples instead of 100, and 10% and 90% 
# percentile instead of default 5%,50%, and 95% with 100 samples
EProc$sEstimateUstarScenarios( 
  seasonFactor = seasonFactor, nSample = 30L, probs = c(0.1,0.9))
#(uStarScens <- usGetSeasonalSeasonUStarMap(
#  EProc$sGetEstimatedUstarThresholdDistribution()
#))
#EProc$sSetUstarScenarios(uStarScens)
EProc$useSeaonsalUStarThresholds()
EProc$sGetUstarScenarios()
```	

After setting the uStar scenarios, they will be used by method 
`sMDSGapFillUStarScens`.
```{r DEGeb_gapFillBoot1, cache = TRUE, message = FALSE, output = 'hide'}
EProc$sMDSGapFillUStarScens('NEE', FillAll = FALSE)
```

Additional output columns are produced for each uStar quantile. 

```{r}
grep("^NEE.*_f$", colnames( EProc$sExportResults()), value = TRUE )
```		

Several methods provide processing steps for all scenarios.
In addition, method `sApplyUStarScen` calls a user-specified function repeatedly
with a modified argument `suffix`. 
For example, the flux partitioning across all threshold scenarios can be
invoked by the following code to produce several GPP estimates.

```{r DEGeb_fluxPart1, cache = TRUE, message = FALSE}
EProc$sMDSGapFill('Tair', FillAll = FALSE)
EProc$sApplyUStarScen( EProc$sMRFluxPartition )
#grep("U10", colnames(EProc$sExportResults()), value = TRUE) 	
grep("^GPP.*_f$", colnames( EProc$sExportResults()), value = TRUE )
```

Using change point detection instead of moving point method for UStar Threshold estimation
------------------------------------------------------------------------------------------

The package also provides another method of estimating the point where NEE saturates 
with increasing uStar.
With the ChangePointDetection (CPT) method, the data is not binned by classes of 
uStar but the changepoint
is estimated based on the entire subset within one seasons and one temperature class.
The user invokes this method by specifying argument 
`ctrlUstarEst.l = usControlUstarEst(isUsingCPTSeveralT = TRUE)`
to `EstUstarThold` or `sEstUstarThresholdDistribution`.

The CPT method is usually yields higher thresholds and marks more data as Gap. 
  
```{r DEGeb_estUStarCPT, cache = TRUE}
EProc <- sEddyProc$new(
  'DE-Geb', DEGebExample, c('NEE','Rg','Tair','VPD', 'Ustar'))
resUStar <- EProc$sEstUstarThold(
					ctrlUstarEst = usControlUstarEst(isUsingCPTSeveralT = TRUE)
					, seasonFactor = seasonFactor
			)
#(uStarThCP <- usGetSeasonalSeasonUStarMap(resUStar))
EProc$useSeaonsalUStarThresholds()
EProc$sGetUstarScenarios()
```


		




