01_post-processing.Rmd

# Post-processing of raw eddy-covariance data

To derive turbulent fluxes and other turbulence characteristics from eddy-covariance measurements, the raw data have to be post-processed first. In the following, we will have a look a the required steps -- starting with the theoretical foundation of the method and the necessary physical and technical considerations for the raw data processing. The upcoming analyses (e.g., [02_basic-turbulence-diagnostics.ipynb](https://github.com/noctiluc3nt/ec_analyze/blob/main/notebooks/02_basic-turbulence-diagnostics.ipynb)) build upon this and allow for a more detailed investigation (e.g., [03_quadrant-analysis.ipynb](https://github.com/noctiluc3nt/ec_analyze/blob/main/notebooks/03_quadrant-analysis.ipynb)).


## Eddy-covariance method
How do we derive a flux from the sonic measurements? The eddy-covariance method is based on the Reynolds decomposition, which is applied to every measured quantity $x$ through
$$
x = \overline{x} + x' \quad \textrm{with} \quad \overline{x}:= \frac{1}{t_s} \int_0^{t_s} x \: dt.
$$
Therein, $\overline{x}$ represents the time average and $x'$ the deviation from it. In principle, one would have to consider the ensemble average (i.e., all possible paths in the phase space), but this is replaced by the time average assuming ergodicity. Now, we take a look at the flux, which is defined as the average product of two quantities $\overline{xw}$, and insert the Reynolds decomposition:
$$
    \overline{xw} = \overline{x}\:\overline{w} + \overline{x}\overline{w'} + \overline{x'}\overline{w} + \overline{x'w'} = \overline{x}\:\overline{w} + \overline{x'w'}
$$
The second and third term vanish because of the Reynolds postulates (i.e. $\overline{x'} = 0$). If additionally $\overline{w}=0$, the flux can be represented by the correlation, i.e.,  $\overline{xw} = \overline{x'w'}$. As you see, this procedure involves a lot of assumptions: For  $\overline{x'} = 0$ to be true, we assume that the flow is steady and for $\overline{w}=0$ you either have a homogeneous and flat surface or you need to rotate the sonic coordiante system (see the correction methods below). <br>
Usually, $w$ represents the vertical wind speed and $x$ can be an arbitrary scalar quantity. For the momentum flux, we use $x=u$ (streamwise velocity), for the sensible heat flux $x=T$ (temperature or potential temperature), for the latent heat flux $x=q$ (specific humidity) and for trace gas fluxes $x=c$ (the concentration or mixing ratio of the respective trace gas).

## Raw data handling and corrections

For applying the described eddy-covariance method to the measurements, we have to decide on an averaging time $t_s$,  have to perform quality control of the raw data and have to check the fulfillment of the made assumptions. 

### Choosing a suitable averaging time
As we have seen above, for calculating the fluxes with the eddy-covariance method, a suitable averaging time ($t_s$) should be chosen. For this, several things need to be considered: Generally, one can say that the longer the averaging time, the less likely it is that the conditions are steady (which is an assumption in the EC method), but short averaging times lose low frequency contributions. So usually, an averaging time of 30 minutes is chosen. However, this averaging time can still be too long if the turbulence is very intermittent, which is the case for very stable conditions. Consequently, some studies use the 30 minutes average just for unstable stratification and an averaging time of 1 to 10 minutes for stable conditions. It is recommendable to first study the characteristics of the sampled raw data (e.g., based on spectra, see [04_spectra.ipynb](https://github.com/noctiluc3nt/ec_analyze/blob/main/notebooks/04_spectra.ipynb))) and then choose an averaging time.<br>


### Overview of correction procedures applied to the raw data
- **Despiking** (`despiking`): Removing of spikes in the raw data. For this usually three methods are applied: (1) based on pre-defined thresholds (i.e., an expected temperature or wind speed interval), (2) median deviation (MAD) test, i.e., all measurements lie within a pre-defined range around the median, and (3) based on skewness (3. moment) and kurtosis (4. moment), i.e., the skewness and kurtosis of the considered interval do not exceed pre-defined values. Details see e.g. @Mauder2013.
- **Lag-time correction** (`shift2maxccf`): Lag-time can occur if several loggers are used or the inlet tube of the gas analyzer is very long, such that there is a time offset between the different measured variables. This can be corrected by calculating the maximum cross-correlation and shift the time series according to the lag difference, which is done by the function `shift2maxccf`.
- **Linear detrending** (`pracma::detrending`): Linear detrending by substracting the mean value (however, linear detrending is not consistent with Reynolds decomposition, so it should be used/interpreted with caution, see e.g. @Baldocchi2003). However, linear detrending should be applied before calculating higher order moments (as skewness and kurtosis) or spectra/structure functions.
- **Rotation** (`rotate_double` or `rotate_planar`): Rotation of the sonic coordinate system. Since the orientation of the sonic is arbitrary, the measurements have to be interpreted in a suitable coordinate system. For this, the natural coordinate system is used, which means that it follows the streamlines and the coordinates are then given by the streamwise velocity component $u$, the crosswise velocity component $v$ and the vertical velocity component $w$. To determine the orientation, two main approaches exist: Double rotation (`rotate_double`) alignes the sonic coordinate system with the streamlines for every averaging interval. Hence, this method can be applied in near-real time. Planar fit rotation (`rotate_planar`, @Wilczak2001) aligns the sonic coordinate system with the mean streamlines under the conditions that the mean vertical velocity $\overline{w}$ vanishes. For this, a longer time series (usually the whole measurement campaign) is required and thus near-real time processing is not possible. Generally, it is recommended to use double rotation in relatively simple topography (or for rapidly changing surface conditions, e.g. over growing gras), but planar fit in complex (micro-)topogaphy. The angles around which the sonic coordinate system was rotated also allow to estimate the quality of the measurements (the smaller the angle of rotation, the better and important is that they do not flip signs, in particular for close-to-zero fluxes under very stable stratification).
- **Quality flagging** (`flag_stationarity`,`flag_w`,`flag_distortion`,`flag_most`, e.g. @Foken1996): Several quality flags (in a 0-1-2 system with 0 (good quality), 1 (medium quality), and 2 (recommended to discard)) can be applied. Common flags are: a stationarity flag (`flag_stationarity`), that tests the agreement with the stationarity assumption of the eddy-covariance method, a vertical velocity flag (`flag_w`), which checks that the remaining vertical velocity after the rotation is small, a flow distortion flag (`flag_distortion`), which removes the pre-defined wind directions that are possibly affected/blocked by the mast, an integral turbulence characteristics flag (`flag_itc`), which tests for the agreement with Monin-Obukhov similarity theory. However, these flags have to be applied purpose-oriented and interpreted with caution. Sometimes particularly "poorly flagged" measurements are interesting for investigating turbulence under challenging (and therefore interesting) conditions, e.g. under very stable stratification.
- **SND (and cross-wind) correction** (`SNDcorrection`, @Schotanus1983): Converts the buoyancy flux to sensible heat flux. Since most sonics measure the sound propagation speed, the measured temperature (the so called sonic temperature $T_s$) is similar to the virtual temperature and thus the resulting flux $\overline{w'T_s'}$ represents the buoyancy flux, which is about 10-20 % larger than the sensible heat flux $\overline{w'T'}$. Note, the used constants in the method depend on the measurement device (default for CSAT3).
- **WPL correction** (`WPLcorrection`, @Webb1980): Converts volume- to mass-related quantities for trace gas fluxes. This correction only applies to trace gas concentrations (water vapor, carbon dioxide, methane, ...) measured with a gas analyzer. @Ibrom2007 suggested an adaption for closed-path systems acounting for desynchonized dilution of water vapor and carbon dioxide flux.
- **Unit conversion (density)** (`ppt2rho`): Converts ppt to density. Closed-path gas analyzers measure trace gas concentrations in parts-per-... (ppt: ... thousand, ppm: ... million), which has to be converted to a density based on the respective molar mass.
- **Unit conversion (fluxes)** (`cov2sh`, `cov2lh`,`cov2cf`): Converts covariances to the respective flux in W/m$^2$, i.e., cov(w,T) to sensible heat flux, cov(w,q) to latent heat flux and cov(w,c) to CO$_2$ flux.
- **Gap-filling** (`gapfilling`): Reproduction of timeseries by filling gaps. Gap-filling is necessary if budgets are calculated or a continuous timeseries is desired. Depending on the type of gaps (random gaps, systematic gaps, duration of gaps and available auxillary measurements) different physical and statistical methods for gap-filling exist (see review @Moffat2007). In meteorology, particulary for studying turbulent fluxes under stable conditions statistically, gap-filling is usually not used to avoid introducing "artificial measurements" created under the assumptions of the used gap-filling method.

There are more (and partly controversial) correction methods, details are discussed in @Foken2017. Further correction methods are required when spectra are considered (spectral corrections) or budgets are calculated (e.g. more sophisticated  gap-filling of missing values, to avoid a bias due to systematically discarded data).


## An example post-processing routine
Now, we create an example post-processing routine. For this, we use 3 days of raw eddy-covariance measurements (CSAT3/Campbell Scientific and Li-7200/LI-COR closed-path IR gas analyzer) with a 10 Hz sampling frequency from our station in Finse, Hardangervidda, Norway. The data set is in the folder `data/ec-data_10Hz_raw` and can be downloaded directly or the whole github repository can be cloned with `git clone git@github.com:noctiluc3nt/ec_analyze.git`.


```R
#loading Reddy package
#install.packages("../src/Reddy_0.0.0.9000.tar.gz",repos=NULL,source=TRUE,quiet=TRUE)
library(Reddy)
```


```R
#ec data files
dir_in="../data/ec-data_10Hz_raw" 
files=list.files(dir_in,full.names=TRUE)
nf=length(files)
```

Each given file contains 30 minutes of measurements, such that we just need to average over one file to get the 30 minutes averages and fluxes. In the notebook [04_multiresolution-decomposition.ipynb](https://github.com/noctiluc3nt/ec_analyze/blob/main/notebooks/04_multiresolution-decomposition.ipynb) we will have a look at a method that allows to determine suitable averaging times more accurately. In our example routine, the data is first despiked, than the wind is rotated and the sonic measurements are averaged. For the gas analyzer measurements, the unit is converted to a density. Then the turbulence intensities (standard deviations) and the fluxes (covariances) are calculated. For the sensible heat flux, we need to apply the SND correction additionally, the WPL correction can be applied to the latent heat flux. <br>
Here, the temperature is given directly in $^\circ$C, however, in the direct output of the LI-COR systems (.ghg-files) usually the speed of sound (sos in m/s) is stored, which can be easily converted to temperature with the function `sos2Ts`.


```R
#allocate output
var_out=c("time","u_mean","v_mean","w_mean","ws_mean","wd_mean","T_mean","h2o_mean","co2_mean",
          "u_sd","v_sd","w_sd","T_sd","h2o_sd","co2_sd",
          "cov_uv","cov_uw","cov_vw","cov_wT","cov_h2ow","cov_co2w","cov_wT_snd","cov_rhoH2Ow_wpl",
          "rot_angle1","rot_angle2","flag_stationarity","flag_w","flag_distortion")
nv=length(var_out)
dat=data.frame(array(NA,dim=c(nf,nv)))
colnames(dat)=var_out
```


```R
#postprocessing per file (30 mins)
for (i in 1:nf) {
    tmp=read.table(files[i],sep=",",header=T)
    dat$time[i]=tmp$X[1]
    #despiking
    tmp$T_degC=despiking(tmp$T_degC,-40,30)
    tmp$u_m.s=despiking(tmp$u_m.s,-25,25)
    tmp$v_m.s=despiking(tmp$v_m.s,-25,25)
    tmp$w_m.s=despiking(tmp$w_m.s,-5,5)
    #wind before rotation
    dat$ws_mean[i]=sqrt(mean(tmp$u_m.s,na.rm=T)^2+mean(tmp$v_m.s,na.rm=T)^2)
    dat$wd_mean[i]=(180+180/pi*atan2(mean(tmp$v_m.s,na.rm=T),mean(tmp$u_m.s,na.rm=T)))%%360
    #rotation
    rot_wind=rotate_double(tmp$u_m.s,tmp$v_m.s,tmp$w_m.s)
    tmp$u_m.s=rot_wind$u
    tmp$v_m.s=rot_wind$v
    tmp$w_m.s=rot_wind$w
    dat$rot_angle1[i]=rot_wind$theta
    dat$rot_angle2[i]=rot_wind$phi
    #averaging
    dat$u_mean[i]=mean(tmp$u_m.s,na.rm=T)
    dat$v_mean[i]=mean(tmp$v_m.s,na.rm=T)
    dat$w_mean[i]=mean(tmp$w_m.s,na.rm=T)
    dat$T_mean[i]=mean(tmp$T_degC,na.rm=T)
    #unit conversion for gases
    tmp$rhoH2O=ppt2rho(tmp$H2O_ppt,dat$T_mean[i]+273.15,87000)
    tmp$rhoCO2=ppt2rho(tmp$CO2_ppm/1000,dat$T_mean[i]+273.15,87000,gas="CO2")
    dat$h2o_mean[i]=mean(tmp$rhoH2O,na.rm=T)
    dat$co2_mean[i]=mean(tmp$rhoCO2,na.rm=T)
    #sds
    dat$u_sd[i]=sd(tmp$u_m.s,na.rm=T)
    dat$v_sd[i]=sd(tmp$v_m.s,na.rm=T)
    dat$w_sd[i]=sd(tmp$w_m.s,na.rm=T)
    dat$T_sd[i]=sd(tmp$T_degC,na.rm=T)
    dat$h2o_sd[i]=sd(tmp$rhoH2O,na.rm=T)
    dat$co2_sd[i]=sd(tmp$rhoCO2,na.rm=T)
    #covs
    dat$cov_uw[i]=cov(tmp$u_m.s,tmp$w_m.s,use="pairwise.complete.obs")
    dat$cov_uv[i]=cov(tmp$u_m.s,tmp$v_m.s,use="pairwise.complete.obs")
    dat$cov_vw[i]=cov(tmp$v_m.s,tmp$w_m.s,use="pairwise.complete.obs")
    dat$cov_wT[i]=cov(tmp$T_degC,tmp$w_m.s,use="pairwise.complete.obs")
    dat$cov_h2ow[i]=cov(tmp$rhoH2O,tmp$w_m.s,use="pairwise.complete.obs")
    dat$cov_co2w[i]=cov(tmp$rhoCO2,tmp$w_m.s,use="pairwise.complete.obs")
    #SND correction
    dat$cov_wT_snd[i]=SNDcorrection(dat$T_mean[i],dat$u_mean[i],dat$v_mean[i],dat$w_mean[i],dat$cov_uw[i],dat$cov_vw[i],dat$cov_wTs[i])
    #flagging
    dat$flag_stationarity[i]=flag_stationarity(tmp$w_m.s,tmp$T_degC,nsub=3000)
    dat$flag_w[i]=flag_w(dat$w_mean[i])
    dat$flag_distortion[i]=flag_distortion(tmp$u_m.s,tmp$v_m.s,dir_blocked=c(340,360))
}
```

The calculated covariances have to be converted to the fluxes (unit W/m$^2$), for example with the functions `cov2sh`and `cov2lh`:


```R
#calculate fluxes
dat$sh=cov2sh(dat$cov_wT) #dat$sh=cov2sh(dat$cov_wT_snd)
dat$lh=cov2lh(dat$cov_h2ow)
```

The output now contains the 30 minutes averages, standard deviations and covariances of the measured quantities, which will be used in the following notebooks for some more detailed analysis.


```R
#look at output
head(dat)
```


<table class="dataframe">
<caption>A data.frame: 6 × 30</caption>
<thead>
	<tr><th></th><th scope=col>time</th><th scope=col>u_mean</th><th scope=col>v_mean</th><th scope=col>w_mean</th><th scope=col>ws_mean</th><th scope=col>wd_mean</th><th scope=col>T_mean</th><th scope=col>h2o_mean</th><th scope=col>co2_mean</th><th scope=col>u_sd</th><th scope=col>⋯</th><th scope=col>cov_co2w</th><th scope=col>cov_wT_snd</th><th scope=col>cov_rhoH2Ow_wpl</th><th scope=col>rot_angle1</th><th scope=col>rot_angle2</th><th scope=col>flag_stationarity</th><th scope=col>flag_w</th><th scope=col>flag_distortion</th><th scope=col>sh</th><th scope=col>lh</th></tr>
	<tr><th></th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>⋯</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;lgl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>
</thead>
<tbody>
	<tr><th scope=row>1</th><td>2018-07-20 08:30:00</td><td>2.872084</td><td>0.3645577</td><td>0.1206889</td><td>2.873170</td><td>48.83370</td><td>15.86938</td><td>0.007254132</td><td>0.0006009614</td><td>1.062143</td><td>⋯</td><td>-9.861714e-08</td><td> 0.007904404</td><td>NA</td><td>228.8337</td><td>0.3484284</td><td>0</td><td>1</td><td>NA</td><td> 80.49175</td><td> 66.62497</td></tr>
	<tr><th scope=row>2</th><td>2018-07-20 09:00:00</td><td>2.864793</td><td>0.3848472</td><td>0.1330429</td><td>2.864538</td><td>46.72741</td><td>16.55190</td><td>0.007584267</td><td>0.0005978138</td><td>1.105122</td><td>⋯</td><td>-1.434213e-07</td><td> 0.007210505</td><td>NA</td><td>226.7274</td><td>0.7745417</td><td>0</td><td>1</td><td>NA</td><td>122.26788</td><td>107.00304</td></tr>
	<tr><th scope=row>3</th><td>2018-07-20 09:30:00</td><td>3.996526</td><td>0.4702586</td><td>0.1505031</td><td>4.002522</td><td>63.38033</td><td>17.05704</td><td>0.007472065</td><td>0.0005954894</td><td>1.409495</td><td>⋯</td><td>-1.228453e-07</td><td> 0.007451433</td><td>NA</td><td>243.3803</td><td>0.2888957</td><td>0</td><td>2</td><td>NA</td><td>123.00260</td><td> 84.88036</td></tr>
	<tr><th scope=row>4</th><td>2018-07-20 10:00:00</td><td>4.998016</td><td>0.6313088</td><td>0.1684323</td><td>4.997530</td><td>66.68404</td><td>17.60447</td><td>0.006762097</td><td>0.0005953881</td><td>1.289326</td><td>⋯</td><td>-9.420972e-08</td><td>-0.006557332</td><td>NA</td><td>246.6840</td><td>0.3483120</td><td>0</td><td>2</td><td>NA</td><td>123.44803</td><td> 83.31937</td></tr>
	<tr><th scope=row>5</th><td>2018-07-20 10:30:00</td><td>4.879095</td><td>0.5967435</td><td>0.1828934</td><td>4.880696</td><td>64.57826</td><td>18.08994</td><td>0.005410862</td><td>0.0005948927</td><td>1.469949</td><td>⋯</td><td>-1.257620e-07</td><td>-0.001824407</td><td>NA</td><td>244.5783</td><td>0.6531922</td><td>0</td><td>2</td><td>NA</td><td>145.37751</td><td>138.66441</td></tr>
	<tr><th scope=row>6</th><td>2018-07-20 11:00:00</td><td>5.225037</td><td>0.6579993</td><td>0.1693293</td><td>5.223984</td><td>50.44251</td><td>18.24207</td><td>0.004422858</td><td>0.0005954766</td><td>1.319227</td><td>⋯</td><td>-8.127438e-08</td><td>-0.004796992</td><td>NA</td><td>230.4425</td><td>0.4471074</td><td>0</td><td>2</td><td>NA</td><td>116.60730</td><td> 94.20120</td></tr>
</tbody>
</table>




```R
saveRDS(dat,file="../data/ec-data_30min_processed/processed_data_example.rds")
```

Additionally, the data can be gap-filled with `gapfilling` or averaged to other customized averaging intervals with `accumulate_timeseries`. The function `accumulate_timeseries` is based on the RcppRoll package, which provides fast and efficient "rolling" functions utilizing a C++ interface. <br> 
Generally, it should be noted that different sensors may require different correction methods, see e.g. @Foken2017 for a discussion, so a detailed knowledge of the instrumentation is required.