Exploring turbidity fluctuations

This is an example workflow of how to perform basic turbidity analysis with loadflux package and tidyverse ecosystem.

library(dplyr)
library(purrr)
library(tidyr)
library(tsibble)
library(loadflux)

Hydrological Events

First, we need to demarcate hydrological events

data(djanturb)
df <- hydro_events(dataframe = djanturb,
                   q = discharge,
                   datetime = time,
                   window = 21)

df %>% 
  event_plot(q = discharge,
             datetime = time,
             he = he,
             ssc = ntu,
             y2label = "Turbidity")

Turbidity Index

Then we can calculate Turbidity Index TI for every hydrological event


TI_index <- df %>% 
  group_by(he) %>% 
  nest() %>% 
  mutate(TI = map_dbl(data, ~TI(.x, ntu, time))) %>% 
  select(-data) %>% 
  ungroup()

TI_index
#> # A tibble: 6 x 2
#>      he    TI
#>   <dbl> <dbl>
#> 1     1 0.193
#> 2     2 0.113
#> 3     3 0.147
#> 4     4 0.153
#> 5     5 0.108
#> 6     6 0.137

Hydrological events parameters

To summarize the hydrological events parameters an approach from features package can be used. For this purpose we need to transform our dataframe into tsibble object:

library(tsibble)

df_ts <- df %>% 
  as_tsibble(key = he,
             index = time)

df_ts
#> # A tsibble: 620 x 4 [10m] <Europe/Moscow>
#> # Key:       he [6]
#>       he time                discharge   ntu
#>    <dbl> <dttm>                  <dbl> <dbl>
#>  1     1 2016-06-16 09:30:00      1.23  751.
#>  2     1 2016-06-16 09:40:00      1.23  622.
#>  3     1 2016-06-16 10:00:00      1.23  964 
#>  4     1 2016-06-16 10:10:00      1.23  804.
#>  5     1 2016-06-16 10:30:00      1.25  927.
#>  6     1 2016-06-16 10:40:00      1.26  673.
#>  7     1 2016-06-16 11:00:00      1.27  631.
#>  8     1 2016-06-16 11:10:00      1.28  645.
#>  9     1 2016-06-16 11:30:00      1.28  665.
#> 10     1 2016-06-16 11:40:00      1.29  726.
#> # ... with 610 more rows

Then we can calculate start, end and length of the every hydrological event:

library(feasts)
#> Загрузка требуемого пакета: fabletools

df_ts %>% 
  features(time,
           feat_event)
#> # A tibble: 6 x 4
#>      he start               end                 length        
#>   <dbl> <dttm>              <dttm>              <drtn>        
#> 1     1 2016-06-16 09:30:00 2016-06-17 10:50:00 25.33333 hours
#> 2     2 2016-06-17 11:00:00 2016-06-18 09:50:00 22.83333 hours
#> 3     3 2016-06-18 10:00:00 2016-06-19 12:20:00 26.33333 hours
#> 4     4 2016-06-19 12:30:00 2016-06-20 11:00:00 22.50000 hours
#> 5     5 2016-06-20 12:00:00 2016-06-21 11:20:00 23.33333 hours
#> 6     6 2016-06-21 11:40:00 2016-06-22 10:20:00 22.66667 hours

Or with the help of brolgar and feasts packages we can calculate turbidity statistics, autocorrelation and spectral functions:

library(brolgar)
library(feasts)

df_ts %>% 
  features(ntu, feat_five_num)
#> # A tibble: 6 x 6
#>      he   min   q25   med   q75   max
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1     1 133.  193.   348.  493. 1018.
#> 2     2 126.  222.   306.  439. 1209.
#> 3     3  83.4 118.   154.  226.  520.
#> 4     4  62.4  86.2  114.  196.  451.
#> 5     5  51.8  63.3   76   106.  374.
#> 6     6  81.3 121.   164   262. 1364.

df_ts %>% 
  features(ntu, c(feat_spectral,
                  feat_acf))
#> # A tibble: 6 x 9
#>      he spectral_entropy  acf1 acf10 diff1_acf1 diff1_acf10 diff2_acf1
#>   <dbl>            <dbl> <dbl> <dbl>      <dbl>       <dbl>      <dbl>
#> 1     1            0.477 0.735  3.42     -0.433      0.382      -0.598
#> 2     2            0.787 0.782  1.90     -0.127      0.0851     -0.533
#> 3     3            0.792 0.656  2.70     -0.436      0.217      -0.648
#> 4     4            0.517 0.769  4.26     -0.268      0.272      -0.477
#> 5     5            0.846 0.618  1.44     -0.201      0.107      -0.455
#> 6     6            0.647 0.875  4.24     -0.221      0.132      -0.600
#> # ... with 2 more variables: diff2_acf10 <dbl>, season_acf1 <dbl>