response_time.Rmd

---
title: "Response time correction of Contros HydroC pCO~2~ data"
author: "Jens Daniel Müller"
date:  "`r format(Sys.time(), '%d %B, %Y')`"
output: 
  workflowr::wflow_html:
    number_sections: true
    toc_depth: 3
    toc_float:
      collapsed: false
editor_options:
  chunk_output_type: console
---

```{r global_options, include = FALSE}
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
```

```{r packages}
library(tidyverse)
library(seacarb)
library(broom)
library(lubridate)
library(tibbletime)
library(patchwork)

```


```{r ggplot_theme, include = FALSE}
theme_set(theme_bw())
```


# List of relevant parameters

Response time determination

- Fit interval length
- tau residual threshold
- Mean vs T-dependent tau

Response time correction

- correction algorithm (Bittig vs. Fiedler, Miloshevich, Fietzek)
- smoothing

Quality assesment of response time correction

- Down-Up-difference vs reference value
- Depth interval width for offset calculation
- Max depth for Down-Up-difference
- NA criterion for included Down-Up-difference



# Response time determination

## HydroC sensor settings

The sensor was first run with a low power pump (1W). Later - and for most parts of the expedition - with a stronger (8W) pump. Pumps were switched between recordings (data file: SD_datafile_20180718_170417CO2-0618-001.txt):  

* 2018-07-17;13:08:34
* 2018-07-17;13:08:35

Logging frequency for all measurement modes (Measure, Zero, Flush) was increased in two steps, It was:  

10 sec for all recordings including SD_datafile_20180714_073641CO2-0618-001.txt  

2 sec after change in SD_datafile_20180717_121052CO2-0618-001.txt at:  

* 2018-07-14;07:52:02
* 2018-07-14;07:52:12
* 2018-07-14;07:52:14

1 sec after change in SD_datafile_20180718_170417CO2-0618-001 at:  

* 2018-07-17;12:27:25
* 2018-07-17;12:27:27
* 2018-07-17;12:27:28


```{r pump_switch_date}

pump_switch <- ymd_hms("2018-07-17;13:08:34")

```


## Model fitting

Response times were determined by fitting a nonlinear least-squares model with the `nls` function as described [here](http://douglas-watson.github.io/post/2018-09_exponential_curve_fitting/) by Douglas Watson.

* Flush period length: variable
* Flush period restricted to equilibration phase, avoiding initial gas mixing effects occuring at the start of each Flush period
* only completed Flush periods (duration > 500 sec) included


```{r read_prepare_tm_data}

tm <- read_csv(here::here("data/_merged_data_files/merging_interpolation",
                          "tm.csv"),
               col_types = cols(ID = col_character(),
                                pCO2_analog = col_double(),
                                pCO2_corr = col_double(),
                                Zero = col_factor(),
                                Flush = col_factor(),
                                Zero_counter = col_integer(),
                                deployment = col_integer(),
                                duration = col_double(),
                                mixing = col_character(),
                                lat = col_double(),
                                lon = col_double()))

tm <- tm %>%
  select(date_time, ID, dep, tem, Flush, pCO2_corr, Zero_counter, duration, mixing)

tm_flush <- tm %>%
  filter(Flush == 1, mixing == "equilibration")

rm(tm)

tm_flush <- tm_flush %>% 
  group_by(Zero_counter) %>% 
  mutate(duration = duration- min(duration)) %>% 
  ungroup()

```

### Example plot

An example plot for a `nls` model fitted to pCO~2~ observations during a Flush phase is shown below.

```{r response_determination_example_plot, fig.cap="Example response time determination by non-linear least squares fit to the pCO~2~ recovery signal after zeroing. The vertical line indicates the determined response time tau. The horizontal line indicates 63% of the difference between start and final fitted pCO~2~."}

i_Zero_counter <- 51

tm_flush_counter <- tm_flush %>%
  filter(Zero_counter == i_Zero_counter,
         duration <= 300)

fit <- tm_flush_counter %>%
  nls(pCO2_corr ~ SSasymp(duration, yf, y0, log_alpha), data = .)

tau <- as.numeric(exp(-tidy(fit)[3,2]))
pCO2_corr_end <- as.numeric(tidy(fit)[1,2])
pCO2_corr_start <- as.numeric(tidy(fit)[2,2])
pCO2_corr_delta = pCO2_corr_end - pCO2_corr_start
resid_abs_mean <- mean(abs(resid(fit)))

augment(fit) %>%
  ggplot(aes(duration, pCO2_corr))+
  geom_point(shape=21)+
  geom_line(aes(y = .fitted), col="red")+
  geom_vline(xintercept = tau)+
  geom_hline(yintercept = pCO2_corr_start + 0.63 *(pCO2_corr_delta))+
  labs(y=expression(pCO[2]~(µatm)), x="t (sec)")

ggsave(here::here("output/Plots/Figures_publication/appendix", "tau_fit_example.pdf"),
       width = 180, height = 90, dpi = 300, units = "mm")

rm(tm_flush_counter, fit, i_Zero_counter, tau,
   pCO2_corr_delta, pCO2_corr_end, pCO2_corr_start,
   resid_abs_mean)

```

### Flush intervals

```{r set_flush_duration_intervals}

duration_intervals <- seq(150,500,50)

```


As there was some speculation about the dependence of determined response times ($\tau$) on the chosen duration of the fit interval, the response time $\tau$ was determined for all zeroings and for total durations of:  

`r duration_intervals` secs

### Residual limits

```{r set_residual_limit}

resid_lim <- 1

```

A mean absolute residual threshold of >`r resid_lim`% of final pCO~2~ was defined.

### Fitting algorithm


```{r fit_plot_all_response_times, include=TRUE, eval=FALSE}

pdf(file=here::here("output/Plots/response_time",
                    "tau_determination_pCO2_corr_flushperiods_nls.pdf"),
    onefile = TRUE, width = 7, height = 4)

for (i in unique(tm_flush$Zero_counter)) {
  for (max_duration in duration_intervals) {

    tm_flush_counter <- tm_flush %>%
      filter(Zero_counter == i, duration <= max_duration)
    
    fit <- 
      try(
      tm_flush_counter %>%
          nls(pCO2_corr ~ SSasymp(duration, yf, y0, log_alpha), data = .),
      TRUE)
    
    if (class(fit) == "nls"){
    
      tau <- as.numeric(exp(-tidy(fit)[3,2]))
      pCO2_corr_end <- as.numeric(tidy(fit)[1,2])
      pCO2_corr_start <- as.numeric(tidy(fit)[2,2])
      pCO2_corr_delta = pCO2_corr_end - pCO2_corr_start
      resid_abs_mean <- mean(abs(resid(fit))/pCO2_corr_end)*100
      
      temp <- as_tibble(bind_cols(Zero_counter = i,
                                  duration = max_duration,
                                  date_time = mean(tm_flush_counter$date_time),
                                  dep = mean(tm_flush_counter$dep),
                                  tem = mean(tm_flush_counter$tem),
                                  pCO2_corr = pCO2_corr_end,
                                  tau = tau,
                                  resid = resid_abs_mean))
      
      if (exists("tau_values")){tau_values <- bind_rows(tau_values, temp)}
        else {tau_values <- temp}
      
      if (resid_abs_mean > resid_lim){warn <- "orange"}
        else {warn <- "black"}
      
      print(
      augment(fit) %>%
        ggplot(aes(duration, pCO2_corr))+
        geom_point(col = warn)+
        geom_line(aes(y = .fitted))+
        geom_vline(xintercept = tau)+
        geom_hline(yintercept = pCO2_corr_start + 0.63 *(pCO2_corr_delta))+
        labs(y=expression(pCO[2]~(µatm)), x="Duration of Flush period (s)",
             title = paste("Zero_counter: ", i,
                           "Tau: ", round(tau,1),
                           "Mean absolute residual (%): ", round(resid_abs_mean, 2)))+
        xlim(0,600)
      )
      
    }
    
    
    else {
    
      temp <- as_tibble(bind_cols(Zero_counter = i,
                                  duration = max_duration,
                                  date_time = mean(tm_flush_counter$date_time),
                                  dep = mean(tm_flush_counter$dep),
                                  tem = mean(tm_flush_counter$tem),
                                  pCO2_corr = pCO2_corr_end,
                                  tau = NaN,
                                  resid = NaN))
      
      if (exists("tau_values")){tau_values <- bind_rows(tau_values, temp)}
        else {tau_values <- temp}
      
      print(
      tm_flush_counter %>%
        ggplot(aes(duration, pCO2_corr))+
        geom_point(col="red")+
        labs(y=expression(pCO[2]~(µatm)), x="Duration of Flush period (s)",
             title = paste("Zero_counter: ", i,
                           "nls model failed"))+
        xlim(0,600)
      )
      
    }
  }
}

dev.off()

rm(tm_flush_counter, fit, i, tau, pCO2_corr_delta, pCO2_corr_end, pCO2_corr_start, temp, max_duration, resid_abs_mean, warn)

tau_values %>% 
  write_csv(here::here("data/_merged_data_files/response_time",
                       "tau_values.csv"))

rm(tau_values, tm_flush)

```

A pdf with plots of all individual response time fits can be accessed [here](https://github.com/jens-daniel-mueller/BloomSail/tree/master/output/Plots/response_time/tau_determination_pCO2_corr_flushperiods_nls.pdf){target="_blank"}. In this pdf, response time fits that exceed the residual criterion (Mean absolute residual >`r resid_lim`% of final pCO~2~) are printed in orange. Data from flush periods without succesful fit are printed red.


```{r apply_residual_thresholds_for_tau, results='hide'}

tau_values <-
 read_csv(here::here("data/_merged_data_files/response_time",
                       "tau_values.csv"))

# define periods of different pumps used
max_Zero_counter <- max(unique(tau_values[tau_values$date_time < pump_switch,]$Zero_counter))

tau_values <- tau_values %>% 
  mutate(pump_power = if_else(Zero_counter <= max_Zero_counter, "1W", "8W"))


# subset determined tau values by residual threshold
tau_resid <- tau_values %>% 
  group_by(Zero_counter) %>% 
  mutate(resid_max = max(resid, na.rm = TRUE)) %>% 
  filter(resid_max < resid_lim) %>% 
  select(-resid_max) %>% 
  ungroup()

tau_resid_out <- tau_values %>% 
  group_by(Zero_counter) %>% 
  mutate(resid_max = max(resid, na.rm = TRUE)) %>% 
  filter(resid_max > resid_lim) %>% 
  select(-resid_max) %>% 
  ungroup()

# Flush periods where model failure occured
tau_values %>% 
  filter(is.na(resid)) %>% 
  group_by(Zero_counter) %>% 
  summarise(n()) %>% 
  ungroup()

# Flush periods removed due to residual criterion
tau_resid_out %>% 
  group_by(Zero_counter) %>% 
  summarise(n()) %>% 
  ungroup()

# mean tau of first RT determination
tau_resid %>% 
  filter(Zero_counter == 2) %>% 
  summarise(tau = mean(tau))

# mean tau of all RT determinations before pump switch, except first
tau_resid %>% 
  filter(Zero_counter != 2, Zero_counter <= 20) %>% 
  summarise(tau = mean(tau))

tau_resid <- tau_resid %>% 
  filter(Zero_counter != 2)

# calculate some metrics for the subsetting
n_Zero_counters <- tau_values %>% 
  group_by(Zero_counter) %>% 
  n_groups()
n_duration_intervals <- length(duration_intervals)
n_tau_max <- n_Zero_counters * length(duration_intervals)
n_tau_total <- nrow(tau_values %>% filter(!is.na(resid)))
n_tau_resid <- nrow(tau_resid)

```


## Outcome

### General considerations

Estimated $\tau$ values were only taken into account when stable environmental pCO~2~ levels were present. Absence of stable environmental pCO~2~ was assumed when the mean absolute fit residual was above `r resid_lim` % of the final equilibrium pCO~2~. If one model fit (irrespective the chosen fit interval length) of a particular flush period did not match that criterion, the flush period was ignored entirely. Usually, fits with higher duration did not meet this criterion. For some unexplained reason the first $\tau$ determination resulted in values about twice as high as all other flush periods and was therefore removed as an outlier.

Metrics to characterize the fitting procedure include the number of:

* Flush periods: `r n_Zero_counters`
* Duration intervals: `r n_duration_intervals`
* Exercised response time fits: `r n_tau_max`
* Succesful response times determinations: `r n_tau_total` (`r round(100*n_tau_total/n_tau_max, 1)`)%
* $\tau$'s after removing groups of fits with high absolute fit residual: `r n_tau_resid` (`r round(100*n_tau_resid/n_tau_max, 1)` %)

It should be noted that all failed model fits occured in flush periods where the residual criterion was not meet by at least one other fit (i.e. fitting only failed under unstable conditions).

`r rm(n_Zero_counters, n_duration_intervals,n_tau_max, n_tau_total, n_tau_resid, n_tau_resid_min_max, resid_abs_mean, tau_resid_out)`

```{r tau_fit_residuals, fig.cap="Histogram of residuals from fit displayed for the investigate durations of the fit interval. Vertical line represents the chosen threshold."}

tau_values %>% 
  ggplot(aes(resid))+
  geom_histogram()+
  facet_wrap(~duration, labeller = label_both)+
  geom_vline(xintercept = resid_lim)+
  labs(x=expression(Mean~absolute~residuals~("%"~of~equilibrium~pCO[2])))

```

### Fit interval length

No clear dependence of $\tau$ on the length of the flushing period was found.

```{r tau_duration_individual, fig.cap="Determined tau values as a function of the fit interval duration, displayed individually for each flush period.", fig.asp=2.5}

tau_resid %>% 
  group_by(Zero_counter) %>% 
  mutate(d_tau = tau - mean(tau)) %>% 
  ggplot(aes(duration, d_tau))+
  geom_hline(yintercept = 0)+
  geom_smooth()+
  geom_point()+
  facet_wrap(~Zero_counter, ncol = 4, labeller = label_both)+
  labs(x="Duration (sec)", y="Deviation from mean tau (sec)")

```


```{r tau_duration_mean, fig.cap="Determined tau values as a function of the fit interval duration, pooled for all flush period."}

tau_resid %>% 
  group_by(Zero_counter) %>% 
  mutate(d_tau = tau - mean(tau)) %>% 
  ggplot(aes(duration, d_tau, group=duration))+
  geom_violin()+
  geom_point()+
  labs(x="Duration (sec)", y="Deviation from mean tau (sec)")+
  facet_wrap(~pump_power)

```

```{r subset_best_duration}

duration_min_tau_sd <- tau_resid %>%
  group_by(Zero_counter) %>%
  mutate(d_tau = tau - mean(tau)) %>%
  ungroup() %>%
  group_by(duration) %>%
  summarise(d_tau_sd = sd(d_tau, na.rm = TRUE)) %>%
  ungroup() %>% 
  slice(which.min(d_tau_sd)) %>%
  select(duration) %>%
  pull()

tau_resid <- tau_resid %>% 
  filter(duration == duration_min_tau_sd) %>% 
  select(-duration)

rm(duration_min_tau_sd)

```



### Time series

```{r tau_vs_date, fig.cap="Tau for all Zeroings with color representing water depth. Red lines represent linear regression trends for tau determined in surface waters (<10m)."}

ggplot()+
  geom_smooth(data = tau_resid %>% filter(dep < 10),
              aes(date_time, tau, linetype = pump_power),
              method = "lm",
              se = FALSE,
              col="red")+
  geom_point(data = tau_resid,
             aes(date_time, tau, col=dep, shape=pump_power))+
  scale_color_viridis_c(name="Depth (m)")+
  labs(y="Tau (sec)", x="Date")+
  ylim(0,NA)

```



### Temperature dependence

A temperature dependence of determined response times $\tau$ was found, with similar slopes but different intercepts for both pumps used.

```{r tau_vs_tem, fig.cap="Tau as a function of temperature for all zeroings determined with low power (left) and strong (right) pump. Color represents the water depth."}
tau_resid %>% 
  ggplot(aes(tem, tau, col=dep))+
  geom_smooth(method = "lm")+
  geom_point()+
  scale_color_viridis_c(name="Depth (m)")+
  labs(y="Tau (sec)", x="Temperature (deg C)")+
  facet_wrap(~pump_power, labeller = label_both)

```

For the response times determined near the surface (<10m, restricted temperature range), no clear temperature dependence of $\tau$ was detected. 

```{r tau_vs_tem_surface, fig.cap="Surface tau (<10m) as a function of temperature for all zeroings determined with low power (left) and strong (right) pump. Color represents the water depth."}
tau_resid %>% 
  filter(dep < 10) %>% 
  ggplot(aes(tem, tau, col=dep))+
  geom_smooth(method = "lm")+
  geom_point()+
  scale_color_viridis_c(name="Depth (m)")+
  labs(y="Tau (sec)", x="Temperature (deg C)")+
  facet_wrap(~pump_power, labeller = label_both)

```


## Final tau

Finally, the mean response times are:

```{r mean_response_time}

tau_mean <- tau_resid %>% 
  group_by(pump_power) %>% 
  summarise(tau = mean(tau, na.rm = TRUE))

tau_mean

```

But we can also fit response times as a function of water temperature:

```{r linear_model_response_time}

tau_fit <- tau_resid %>% 
  group_by(pump_power) %>% 
  do(fit = lm(tau ~ tem, data = .)) %>% 
  tidy(fit) %>% 
  select(pump_power, term, estimate) %>% 
  spread(term, estimate)

tau_fit

tau_fit %>% write_csv(here::here("data/_merged_data_files/response_time",
                                "tau_fit.csv"))

rm(list=setdiff(ls(), c("tau_resid", "tau_fit", "pump_switch")))

```

Both response time estimated (constant mean vs T-dependent) will be applied to correct the recorded pCO~2~ profiles.

# Response time correction

## Data preparation

Following tasks were performed to prepare data for the response time correction:

* Select only profiles
* Assign deployment periods with 1W- and 8W- pump


```{r response_correction_data_preparation}

tm <- read_csv(here::here("data/_merged_data_files/merging_interpolation",
                          "tm.csv"),
               col_types = cols(ID = col_character(),
                                pCO2_analog = col_double(),
                                pCO2_corr = col_double(),
                                Zero = col_factor(),
                                Flush = col_factor(),
                                Zero_counter = col_integer(),
                                deployment = col_integer(),
                                duration = col_double(),
                                mixing = col_character(),
                                lat = col_double(),
                                lon = col_double()))

# extract relevant parts

tm <- tm %>%
  select(date_time, ID, type, station, dep, sal, tem,
         Zero, Flush, pCO2_corr, deployment, Zero_counter)

tm <- tm %>% 
  filter(type == "P")

tm <- tm %>%
  group_by(ID, station) %>% 
  mutate(duration = as.numeric(date_time - min(date_time)),
         pump_power = if_else(date_time < pump_switch, "1W", "8W")) %>%
  arrange(date_time)

```

* Include manually derived meta-information about the profiling status

```{r include_profile_meta_data}
# Load profile meta data 

meta <- read_csv(here::here("Data/TinaV/Sensor",
                            "Sensor_meta.csv"),
                 col_types = cols(ID = col_character()))


# Merge data and meta information 

tm <- full_join(tm, meta)
rm(meta)


# creating descriptive variables ------------------------------------------
tm <- tm %>% 
  mutate(phase = "standby",
         phase = if_else(duration >= start &
                           duration < down & 
                           !is.na(down) &
                           !is.na(start),
                         "down", phase),
         phase = if_else(duration >= down  & 
                           duration < lift & 
                           !is.na(lift) & 
                           !is.na(down ),   
                         "low",  phase),
         phase = if_else(duration >= lift  & 
                           duration < up   & 
                           !is.na(up  ) & 
                           !is.na(lift  ),  
                         "mid",  phase),
         phase = if_else(duration >= up    & 
                           duration < end  & 
                           !is.na(end ) & 
                           !is.na(up   ),   
                         "up",   phase))

tm <- tm %>% 
  select(-c(start, down, lift, up, end, comment))

tm <- tm %>% 
  filter(Zero == 0, Flush == 0)


```

* Subset reference pCO~2~ recordings at the end of equilibration periods executed at constant depth

```{r equi_pCO2_reference_points}

tm_pCO2_equi <- tm %>% 
  filter(phase %in% c("mid")) %>% 
  group_by(ID, station) %>% 
  top_n(5, row_number()) %>% 
  summarise(date_time = mean(date_time),
            duration = mean(duration),
            pCO2_corr = mean(pCO2_corr, na.rm = TRUE),
            dep = mean(dep, na.rm = TRUE)) %>% 
  ungroup()


tm_pCO2_equi %>% 
 write_csv(here::here("data/_merged_data_files/response_time",
                      "tm_pCO2_equi.csv"))
```

* Select profile for example plots

```{r example_profile}

example_ID <- "180802"
example_station <- "P03"

```


```{r plot_profiling_timeseries_example, fig.cap="Example timeseries of profiling depth and pCO~2~. Colors represent manually assigned profiling phases. The black points represent reference data collected at the end of the mid equilibration period."}

cast_dep <- tm %>% 
  pivot_longer(c(dep, pCO2_corr), names_to = "parameter", values_to = "value")

cast_dep_equi <- tm_pCO2_equi %>% 
  pivot_longer(c(dep, pCO2_corr), names_to = "parameter", values_to = "value")

max_duration <- round(max(cast_dep$duration)/1000,0)*1000


i_ID <- example_ID
i_station <- example_station

cast_dep_equi_sub <- cast_dep_equi %>%
  filter(ID == i_ID,
         station == i_station)


cast_dep %>%
  filter(ID == i_ID,
         station == i_station) %>%
  ggplot(aes(duration, value, col=phase))+
    geom_point(size=0.5)+
    geom_point(data = cast_dep_equi_sub, aes(duration, value), col="black")+
    scale_y_reverse()+
    scale_x_continuous(breaks = seq(0,6000,500))+
    labs(title = str_c("Date: ",i_ID," | Station: ",i_station))+
    facet_grid(parameter~., scales = "free_y")
      

rm(cast_dep, cast_dep_equi, cast_dep_equi_sub, i_station, i_ID, max_duration)

```

A pdf with all timeseries plots of profiling depth and pCO~2~ can be accessed [here](https://github.com/jens-daniel-mueller/BloomSail/tree/master/output/Plots/response_time/time_series_depth_pCO2_corr_by_profile.pdf){target="_blank"}.


```{r plot_profiling_timeseries, eval=FALSE}

cast_dep <- tm %>% 
  pivot_longer(c(dep, pCO2_corr), names_to = "parameter", values_to = "value")

cast_dep_equi <- tm_pCO2_equi %>% 
  pivot_longer(c(dep, pCO2_corr), names_to = "parameter", values_to = "value")

max_duration <- round(max(cast_dep$duration)/1000,0)*1000

pdf(file=here::here("output/Plots/response_time",
                    "time_series_depth_pCO2_corr_by_profile.pdf"),
    onefile = TRUE, width = 7, height = 4)

for(i_ID in unique(cast_dep$ID)){
  for(i_station in unique(cast_dep$station)){

    if (nrow(cast_dep %>% filter(ID == i_ID, station == i_station)) > 0){
      
        cast_dep_equi_sub <- cast_dep_equi %>%
        filter(ID == i_ID,
               station == i_station)

      print(
      
        cast_dep %>%
        filter(ID == i_ID,
               station == i_station) %>%
        ggplot(aes(duration, value, col=phase))+
        geom_point(size=0.5)+
        geom_point(data = cast_dep_equi_sub, aes(duration, value), col="black")+
        scale_y_reverse()+
        scale_x_continuous(breaks = seq(0,6000,500))+
        labs(title = str_c("Date: ",i_ID," | Station: ",i_station))+
        facet_grid(parameter~., scales = "free_y")+
        theme_bw()
      
      )

    }

  }
}

dev.off()

rm(cast_dep, cast_dep_equi, cast_dep_equi_sub, i_station, i_ID, max_duration)

```


## Reponse time correction

The executed response time correction featured the following aspects:

* Correction according to Bittig et al. (2018, supplement)
* RT: Constant mean vs. T-dependent response times applied (both independently quantified for 1W- and 8W-pump)
* tau_factor: Factor ranging from 0.8 - 1.6 applied to determined tau values
* Post-smoothing: 30 sec running mean (eg across 15 observations at 2 sec measurement frequency) 

```{r RT_correction, eval=FALSE}

# Response time correction approach after Bittig --------------------------

RT_corr <- function(c1, c0, dt, tau) {
  ( 1 / ( 2* (( 1+(2*tau/dt) )^(-1) ))) * (c1 - (1-(2* (( 1+(2*tau/dt) )^(-1) ))) * c0)
}


# Assign mean response time (tau) values ----------------------------------------------

# df_mean <- full_join(df, tau_mean)
# 
# df_mean <- df_mean %>% 
#   mutate(RT = "constant")

# Assign T-dependent response time (tau) values ----------------------------------------------

tau_fit <- tau_fit %>% 
  rename(tau_intercept = `(Intercept)`, tau_slope=tem)

tm_fit <- full_join(tm, tau_fit)

tm_fit <- tm_fit %>% 
  mutate(tau = tau_intercept + tau_slope *tem) %>% 
  select(-tau_intercept, -tau_slope)

# tm_fit <- tm_fit %>% 
#   mutate(RT = "T-dependent")


# Merge data sets with constand and T-dependent tau

tm <- tm_fit
#tm <- bind_rows(tm_fit, tm_mean)

rm(tm_fit)

# p1 <- tm %>% 
#   ggplot(aes(tau, dep, col=pump_power))+
#   geom_point()+
#   scale_y_reverse()
# 
# p2 <- tm %>% 
#   ggplot(aes(tem, dep, col=pump_power))+
#   geom_point()+
#   scale_y_reverse()
# 
# p1 | p2
# 
# rm(p1, p2)


# Prepare data set for RT correction --------------------------------------

tm <- tm %>% 
  #group_by(RT) %>% 
  arrange(date_time) %>% 
  mutate(dt = as.numeric(as.character(date_time - lag(date_time))))
 # ungroup()

# measurement frequency

freq <- tm %>% 
  filter(dt < 13) %>% 
  group_by(ID) %>% 
  summarise(dt_mean = round(mean(dt, na.rm = TRUE),0))

tm <- full_join(tm, freq)

# tau factors

tm <- expand_grid(tm, tau_factor = seq(0.8, 1.6, 0.2))

tm <- tm %>%
  mutate(tau_test = tau*tau_factor)


# Apply RT correction to entire data set

for(i_ID in unique(tm$ID)){
  
#i_ID <- "180716"
  
  freq_sub <- freq %>% filter(ID == i_ID) %>% pull(dt_mean)
  
  window <- 30 / freq_sub
  rolling_mean   <- rollify(~mean(.x, na.rm = TRUE), window = window)

  tm_sub <- tm %>%
    filter(ID == i_ID) %>% 
    group_by(station, tau_factor) %>% 
    mutate(pCO2_RT = RT_corr(pCO2_corr, lag(pCO2_corr), dt, tau_test),
           pCO2_RT = if_else(pCO2_RT %in% c(Inf, -Inf), NaN, pCO2_RT),
           window = window,
           pCO2 = rolling_mean(pCO2_RT)
           #pCO2_RT_median = rolling_median(pCO2_RT)
           ) %>% 
    ungroup()
  
  # time shift RT corrected data
  shift <- as.integer(as.character(window/2))
  
  tm_sub <- tm_sub %>%
    group_by(station, tau_factor) %>% 
    mutate(pCO2 = lead(pCO2, shift)) %>% 
    ungroup()
  
  
  if (exists("tm_RT")){tm_RT <- bind_rows(tm_RT, tm_sub)}
  else{tm_RT <- tm_sub}
  
  rm(tm_sub, freq_sub, rolling_mean, shift, window)
  
}

tm_RT %>% 
 write_csv(here::here("data/_merged_data_files/response_time",
                      "tm_RT_profiles_by_tau_factor.csv"))

rm(i_ID, freq, RT_corr, tm_RT, tau_fit, tau_resid)

```

As an alternative, the high resolution RT correction can be run.

```{r RT_correction_high_res, eval=FALSE}

# Response time correction approach after Bittig --------------------------

RT_corr <- function(c1, c0, dt, tau) {
  ( 1 / ( 2* (( 1+(2*tau/dt) )^(-1) ))) * (c1 - (1-(2* (( 1+(2*tau/dt) )^(-1) ))) * c0)
}


# Assign mean response time (tau) values ----------------------------------------------

# df_mean <- full_join(df, tau_mean)
# 
# df_mean <- df_mean %>% 
#   mutate(RT = "constant")

# Assign T-dependent response time (tau) values ----------------------------------------------

tau_fit <- tau_fit %>% 
  rename(tau_intercept = `(Intercept)`, tau_slope=tem)

tm_fit <- full_join(tm, tau_fit)

tm_fit <- tm_fit %>% 
  mutate(tau = tau_intercept + tau_slope *tem) %>% 
  select(-tau_intercept, -tau_slope)

# tm_fit <- tm_fit %>% 
#   mutate(RT = "T-dependent")


# Merge data sets with constand and T-dependent tau

tm <- tm_fit
#tm <- bind_rows(tm_fit, tm_mean)

rm(tm_fit)

# p1 <- tm %>% 
#   ggplot(aes(tau, dep, col=pump_power))+
#   geom_point()+
#   scale_y_reverse()
# 
# p2 <- tm %>% 
#   ggplot(aes(tem, dep, col=pump_power))+
#   geom_point()+
#   scale_y_reverse()
# 
# p1 | p2
# 
# rm(p1, p2)


# Prepare data set for RT correction --------------------------------------

tm <- tm %>% 
  #group_by(RT) %>% 
  arrange(date_time) %>% 
  mutate(dt = as.numeric(as.character(date_time - lag(date_time))))
 # ungroup()

# measurement frequency

freq <- tm %>% 
  filter(dt < 13) %>% 
  group_by(ID) %>% 
  summarise(dt_mean = round(mean(dt, na.rm = TRUE),0))

tm <- full_join(tm, freq)

# tau factors

tm <- expand_grid(tm, tau_factor = seq(0.8, 1.6, 0.05))

tm <- tm %>%
  mutate(tau_test = tau*tau_factor)


# Apply RT correction to entire data set

for(i_ID in unique(tm$ID)){
  
#i_ID <- "180716"
  
  freq_sub <- freq %>% filter(ID == i_ID) %>% pull(dt_mean)
  
  window <- 30 / freq_sub
  rolling_mean   <- rollify(~mean(.x, na.rm = TRUE), window = window)

  tm_sub <- tm %>%
    filter(ID == i_ID) %>% 
    group_by(station, tau_factor) %>% 
    mutate(pCO2_RT = RT_corr(pCO2_corr, lag(pCO2_corr), dt, tau_test),
           pCO2_RT = if_else(pCO2_RT %in% c(Inf, -Inf), NaN, pCO2_RT),
           window = window,
           pCO2 = rolling_mean(pCO2_RT)
           #pCO2_RT_median = rolling_median(pCO2_RT)
           ) %>% 
    ungroup()
  
  # time shift RT corrected data
  shift <- as.integer(as.character(window/2))
  
  tm_sub <- tm_sub %>%
    group_by(station, tau_factor) %>% 
    mutate(pCO2 = lead(pCO2, shift)) %>% 
    ungroup()
  
  
  if (exists("tm_RT")){tm_RT <- bind_rows(tm_RT, tm_sub)}
  else{tm_RT <- tm_sub}
  
  rm(tm_sub, freq_sub, rolling_mean, shift, window)
  
}

tm_RT %>% 
 write_csv(here::here("data/_merged_data_files/response_time",
                      "tm_RT_profiles_by_tau_factor_high_res.csv"))

rm(i_ID, freq, RT_corr, tm_RT, tau_fit, tau_resid, tm)

```




```{r plot_optimized_RT_corrected_profiles_example, fig.cap="Example plot of response time corrected and raw pCO~2~ profiles. Panels highlight the effect of constant vs T-dependent tau estimates (columns) and the optimization by applying a constant factor (rows).", fig.asp = 2.3}

tm <- read_csv(here::here("data/_merged_data_files/response_time",
                          "tm_RT_profiles_by_tau_factor.csv"),
               col_types = cols(ID = col_character(),
                                Zero = col_factor(),
                                Flush = col_factor(),
                                p_type = col_factor(),
                                Zero_counter = col_integer(),
                                deployment = col_integer()))

i_ID <- example_ID
i_station <- example_station

  equi_cast <- tm_pCO2_equi %>%
    filter(ID == i_ID,
           station == i_station)

  tm %>%
    filter(ID == i_ID,
           station == i_station,
           phase %in% c("up", "down")) %>%
    ggplot()+
    geom_path(aes(pCO2_corr, dep, linetype = phase, col="raw"))+
    geom_path(aes(pCO2, dep, linetype = phase, col="corrected"))+
    geom_point(data = equi_cast, aes(pCO2_corr, dep))+
    scale_y_reverse()+
    coord_cartesian(ylim = c(35,0),
                    xlim = c(70, 270))+
    scale_color_brewer(palette = "Set1", name="")+
    labs(title = str_c("Date: ",i_ID," | Station: ",i_station))+
    facet_grid(tau_factor~., labeller = label_both)

  
rm(equi_cast)

```

A pdf with all timeseries plots of profiling depth and pCO~2~ can be accessed [here](https://github.com/jens-daniel-mueller/BloomSail/tree/master/output/Plots/response_time/RT_correction_pCO2_profiles.pdf){target="_blank"}



```{r plot_optimized_RT_corrected_profiles, eval=FALSE}

pdf(file=here::here("output/Plots/response_time",
    "profiles_pCO2.pdf"), onefile = TRUE, width = 7, height = 11)

for(i_ID in unique(tm$ID)){
  for(i_station in unique(tm$station)){

    if (nrow(tm %>% filter(ID == i_ID, station == i_station)) > 0){
      
      equi_cast <- tm_pCO2_equi %>%
        filter(ID == i_ID,
               station == i_station)
      
      print(

  tm %>%
    filter(ID == i_ID,
           station == i_station,
           phase %in% c("up", "down")) %>%
    ggplot()+
    geom_path(aes(pCO2_corr, dep, linetype = phase, col="raw"))+
    geom_path(aes(pCO2, dep, linetype = phase, col="corrected"))+
    geom_point(data = equi_cast, aes(pCO2_corr, dep))+
    scale_y_reverse()+
    coord_cartesian(ylim = c(35,0),
                    xlim = c(70, 270))+
    scale_color_brewer(palette = "Set1", name="")+
    labs(title = str_c("Date: ",i_ID," | Station: ",i_station))+
    theme_bw()+
    facet_grid(tau_factor~., labeller = label_both)
      
      )


    }

  }
}

dev.off()

rm(equi_cast, i_ID, i_station)

```


# Diagnosis

In the following, the success of the response time correction is assessed through the offset between the **downcast** and:

1. Upcast
2. pCO~2~ reference value (recorded after equilibration period during upcast)

The offset comparison requires to discretize the depth recording. **Depth intervals of 1m** were chosen.

First, we analyse all profiles individually. Later we'll merge the information across profiles and come up with a single metric to quantive the quality of the response time correction

## Down- vs upcast


```{r compute_pCO2_offset_profiles, eval=FALSE}

# pCO2 offset up - down cast
tm_grid <- tm %>% 
  filter(phase %in% c("down", "up")) %>% 
  mutate(dep_grid = as.numeric(as.character( cut(dep, seq(0,40,1), seq(0.5,39.5,1)))),
         tau_factor = as.factor(tau_factor)) %>% 
  select(ID, station, tau_factor, p_type, dep_grid, phase, pCO2_corr, pCO2) %>% 
  group_by(ID, station, tau_factor, p_type, dep_grid, phase) %>%
  summarise_all("mean", na.rm = TRUE) %>% 
  ungroup() %>% 
  pivot_longer(cols = c(pCO2_corr, pCO2), names_to = "correction") %>% 
  pivot_wider(names_from = phase, values_from = value) %>% 
  mutate(pCO2_delta = up - down,
         pCO2_up_down_average = (down + up)/2,
         pCO2_delta_rel = 100 * pCO2_delta / pCO2_up_down_average)


tm_high_res <- read_csv(here::here("data/_merged_data_files/response_time",
                          "tm_RT_profiles_by_tau_factor_high_res.csv"),
               col_types = cols(ID = col_character(),
                                Zero = col_factor(),
                                Flush = col_factor(),
                                p_type = col_factor(),
                                Zero_counter = col_integer(),
                                deployment = col_integer()))


# pCO2 offset up - down cast
tm_grid_high_res <- tm_high_res %>% 
  filter(phase %in% c("down", "up")) %>% 
  mutate(dep_grid = as.numeric(as.character( cut(dep, seq(0,40,1), seq(0.5,39.5,1)))),
         tau_factor = as.factor(tau_factor)) %>% 
  select(ID, station, tau_factor, p_type, dep_grid, phase, pCO2_corr, pCO2) %>% 
  group_by(ID, station, tau_factor, p_type, dep_grid, phase) %>%
  summarise_all("mean", na.rm = TRUE) %>% 
  ungroup() %>% 
  pivot_longer(cols = c(pCO2_corr, pCO2), names_to = "correction") %>% 
  pivot_wider(names_from = phase, values_from = value) %>% 
  mutate(pCO2_delta = up - down,
         pCO2_up_down_average = (down + up)/2,
         pCO2_delta_rel = 100 * pCO2_delta / pCO2_up_down_average)



# descritize depth recordings of equilibrated pCO2 values

tm_pCO2_equi_grid <- tm_pCO2_equi %>%
  mutate(dep_grid = as.numeric(as.character(cut(dep, seq(0,40,1), seq(0.5,39.5,1))))) %>% 
  select(ID, station, dep_grid, pCO2_equi=pCO2_corr)


tm_grid %>% 
 write_csv(here::here("data/_merged_data_files/response_time",
                      "tm_grid.csv"))

tm_grid_high_res %>% 
 write_csv(here::here("data/_merged_data_files/response_time",
                      "tm_grid_high_res.csv"))

tm_pCO2_equi_grid %>% 
 write_csv(here::here("data/_merged_data_files/response_time",
                      "tm_pCO2_equi_grid.csv"))


rm(tm_grid, tm_pCO2_equi_grid, tm_grid_high_res)


```

```{r plot_corrected_discretized_pCO2_profiles_example, fig.cap="Example plot of discretized, response time corrected and raw pCO~2~ profiles. Panels highlight the effect of constant vs T-dependent tau estimates (columns) and the optimization by applying a constant factor (rows). The black point indicates the reference pCO~2~ value.", fig.asp = 2.3}


tm_grid <-
 read_csv(here::here("data/_merged_data_files/response_time",
                      "tm_grid.csv"),
               col_types = cols(ID = col_character()))

tm_grid_high_res <-
 read_csv(here::here("data/_merged_data_files/response_time",
                      "tm_grid_high_res.csv"),
               col_types = cols(ID = col_character()))

tm_pCO2_equi_grid <-
 read_csv(here::here("data/_merged_data_files/response_time",
                      "tm_pCO2_equi_grid.csv"),
          col_types = cols(ID = col_character()))



i_ID <- example_ID
i_station <- example_station

tm_pCO2_equi_grid_sub <- tm_pCO2_equi_grid %>%
  filter(ID == i_ID,
         station == i_station)


      tm_grid %>%
        filter(ID == i_ID,
               station == i_station) %>%
        arrange(dep_grid) %>%
        ggplot()+
        geom_path(aes(down, dep_grid, col=correction, linetype="down"))+
        geom_path(aes(up, dep_grid, col=correction, linetype ="up"))+
        geom_point(data = tm_pCO2_equi_grid_sub, aes(pCO2_equi, dep_grid))+
        scale_y_reverse()+
        coord_cartesian(ylim = c(35,0),
                    xlim = c(70, 270))+
        scale_linetype(name="cast")+
        scale_color_brewer(palette = "Set1", direction = -1)+
        labs(y="Depth [m]", x=expression(pCO[2]~(µatm)),
             title = str_c("Date: ",i_ID," | Station: ",i_station))+
        facet_grid(tau_factor~., labeller = label_both)
      
rm(tm_pCO2_equi_grid_sub, i_ID, i_station)

```


```{r plot_corrected_discretized_pCO2_profiles_appendix, include=FALSE}

i_ID <- example_ID
i_station <- example_station
i_tau_factor <- 1.25

tm_grid_high_res %>%
  filter(ID == i_ID,
         station == i_station,
         tau_factor == i_tau_factor) %>%
  arrange(dep_grid) %>%
  ggplot()+
  geom_path(aes(down, dep_grid, col=correction, linetype="down"))+
  geom_path(aes(up, dep_grid, col=correction, linetype ="up"))+
  scale_y_reverse()+
  scale_linetype(labels = c("upcast", "downcast"))+
  scale_color_brewer(palette = "Set1", direction = -1, labels = c("raw", "corrected"))+
  labs(y="Depth (m)", x=expression(pCO[2]~(µatm)))+
  theme(legend.title = element_blank(),
        legend.position = c(0.8,0.8),
        legend.background = element_rect(colour = "black"))

ggsave(here::here("output/Plots/Figures_publication/appendix", "pCO2_RT_corrected_profile_example.pdf"),
       width = 100, height = 140, dpi = 300, units = "mm")
      
rm(i_ID, i_station, i_tau_factor)

```





A pdf with all discretized pCO~2~ profiles can be assessed [here](https://github.com/jens-daniel-mueller/BloomSail/tree/master/output/Plots/response_time/profiles_pCO2_grid.pdf){target="_blank"}


```{r plot_corrected_discretized_pCO2_profiles, eval=FALSE}

pdf(file=here::here("output/Plots/response_time",
    "profiles_pCO2_grid.pdf"), onefile = TRUE, width = 7, height = 11)

for(i_ID in unique(tm_grid$ID)){
  for(i_station in unique(tm_grid$station)){

    if (nrow(tm_grid %>% filter(ID == i_ID, station == i_station)) > 0){
      
      
        tm_pCO2_equi_grid_sub <- tm_pCO2_equi_grid %>%
        filter(ID == i_ID,
               station == i_station)

      print(
      
      tm_grid %>%
        filter(ID == i_ID,
               station == i_station) %>%
        arrange(dep_grid) %>%
        ggplot()+
        geom_path(aes(down, dep_grid, col=correction, linetype="down"))+
        geom_path(aes(up, dep_grid, col=correction, linetype ="up"))+
        geom_point(data = tm_pCO2_equi_grid_sub, aes(pCO2_equi, dep_grid))+
        scale_y_reverse()+
        coord_cartesian(ylim = c(35,0),
                        xlim = c(70, 270))+
        scale_linetype(name="cast")+
        scale_color_brewer(palette = "Set1", direction = -1)+
        labs(y="Depth [m]", title = str_c("Date: ",i_ID," | Station: ",i_station))+
        theme_bw()+
        facet_grid(tau_factor~., labeller = label_both)
      
      )
      
      rm(tm_pCO2_equi_grid_sub)
    }

  }
}

dev.off()

```


```{r plot_delta_pCO2_absolute_profiles_example, fig.cap="Example plot of absolute pCO~2~ offset profiles. Panels highlight the effect of constant vs T-dependent tau estimates. Colour indicates the optimization by applying a constant factor to tau. Vertical red lines mark an arbitray 10 µatm pCO~2~ threshold."}

i_ID <- example_ID
i_station <- example_station


tm_grid %>%
  filter(ID == i_ID,
         station == i_station,
         correction == "pCO2") %>%
  arrange(dep_grid) %>%
  ggplot(aes(pCO2_delta, dep_grid, col=as.factor(tau_factor)))+
  geom_path()+
  geom_point()+
  scale_y_reverse(breaks=seq(0,40,2))+
  scale_color_discrete(name="tau factor")+
  labs(x=expression(Delta~pCO[2]~(µatm)),
       y="Depth [m]",
       title = str_c("Date: ",i_ID," | Station: ",i_station))+
  geom_vline(xintercept = 0)+
  geom_vline(xintercept = c(-10,10), col="red")

rm(i_ID, i_station)

```


```{r plot_delta_pCO2_absolute_profiles, eval=FALSE}

pdf(file=here::here("output/Plots/response_time",
    "profiles_pCO2_delta_grid.pdf"), onefile = TRUE, width = 7, height = 7)

for(i_ID in unique(tm_grid$ID)){
  for(i_station in unique(tm_grid$station)){

    if (nrow(tm_grid %>% filter(ID == i_ID, station == i_station)) > 0){

      print(
      
      tm_grid %>%
        filter(ID == i_ID,
               station == i_station,
               correction == "pCO2") %>%
        arrange(dep_grid) %>%
        ggplot(aes(pCO2_delta, dep_grid, col=as.factor(tau_factor)))+
        geom_path()+
        geom_point()+
        scale_y_reverse(breaks=seq(0,40,2))+
        scale_color_discrete(name="tau factor")+
        labs(x = "delta pCO2 [µatm]", y="Depth [m]", title = str_c("Date: ",i_ID," | Station: ",i_station))+
        geom_vline(xintercept = 0)+
        geom_vline(xintercept = c(-10,10), col="red")+
        theme_bw()
      
      )
      
    }

  }
}

dev.off()

```


A pdf with all absolute pCO~2~ offset profiles can be assessed [here](https://github.com/jens-daniel-mueller/BloomSail/tree/master/output/Plots/response_time/profiles_pCO2_delta_grid.pdf){target="_blank"}.



```{r plot_delta_pCO2_relative_profiles_example, fig.cap="Example plot of relative offset pCO~2~ profiles. Panels highlight the effect of constant vs T-dependent tau estimates. Colour indicates the optimization by applying a constant factor to tau. Vertical red lines mark an arbitray 10% threshold."}

i_ID <- example_ID
i_station <- example_station

      tm_grid %>%
        filter(ID == i_ID,
               station == i_station,
               correction == "pCO2") %>%
        arrange(dep_grid) %>%
        ggplot(aes(pCO2_delta_rel, dep_grid, col=as.factor(tau_factor)))+
        geom_path()+
        geom_point()+
        scale_y_reverse(breaks=seq(0,40,2))+
        scale_color_discrete(name="tau factor")+
        labs(x=expression(Delta~pCO[2]~("%"~of~absolute~value)), y="Depth [m]",
                          title = str_c("Date: ",i_ID," | Station: ",i_station))+
        geom_vline(xintercept = 0)+
        geom_vline(xintercept = c(-10,10), col="red")
      
```

A pdf with all relative pCO~2~ offset profiles can be assessed [here](https://github.com/jens-daniel-mueller/BloomSail/tree/master/output/Plots/response_time/RT_correction_delta_pCO2_relative_profiles_discrete.pdf){target="_blank"}.


```{r plot_delta_pCO2_relative_profiles, eval=FALSE}

pdf(file=here::here("output/Plots/response_time",
    "profiles_pCO2_delta_rel_grid.pdf"), onefile = TRUE, width = 7, height = 7)

for(i_ID in unique(tm_grid$ID)){
  for(i_station in unique(tm_grid$station)){

    if (nrow(tm_grid %>% filter(ID == i_ID, station == i_station)) > 0){

      print(
      
      tm_grid %>%
        filter(ID == i_ID,
               station == i_station,
               correction == "pCO2") %>%
        arrange(dep_grid) %>%
        ggplot(aes(pCO2_delta_rel, dep_grid, col=as.factor(tau_factor)))+
        geom_path()+
        geom_point()+
        scale_y_reverse(breaks=seq(0,40,2))+
        scale_color_discrete(name="tau factor")+
        labs(x = "delta pCO2 [% of absolute value]",
             y="Depth [m]",
             title = str_c("Date: ",i_ID," | Station: ",i_station))+
        geom_vline(xintercept = 0)+
        geom_vline(xintercept = c(-10,10), col="red")+
        theme_bw()
      
      )
      
    }

  }
}

dev.off()

```

## Downcast vs reference value

```{r compute_pCO2_offset_reference_value}

tm_equi_delta <- full_join(tm_grid, tm_pCO2_equi_grid) %>% 
  filter(!is.na(pCO2_equi)) %>% 
  mutate(pCO2_delta_equi = down - pCO2_equi,
         pCO2_delta_equi_rel = 100 * pCO2_delta_equi / pCO2_equi)

```

```{r plot_pCO2_offset_reference_value, fig.cap="Offset between pCO~2~ downcast and upcast reference value as a function of absolute pCO~2~. (Conditions: T-dependent tau, Factor = 1.2."}

tm_equi_delta %>% 
  filter(tau_factor == 1.2, correction=="pCO2") %>% 
  ggplot(aes(pCO2_equi, pCO2_delta_equi))+
  geom_hline(yintercept = 0)+
  geom_point()+
  labs(x=expression(Reference~pCO[2]~(µatm)), y=expression(Delta~pCO[2]~from~reference~(µatm)))

```


```{r plot_pCO2_offset_reference_value_violin, fig.cap="Offset between pCO~2~ downcast and upcast reference value. Panels highlight the effect of constant vs T-dependent tau estimates. Colour distinguish raw and corrected offsets."}

tm_equi_delta %>% 
  ggplot(aes(as.factor(tau_factor), pCO2_delta_equi, fill=correction))+
  geom_hline(yintercept = 0)+
  geom_violin()+
  #labs(y=expression(Delta~pCO[2]~from~reference~(µatm)), x="Tau factor")+
  scale_fill_brewer(palette = "Set1")

```


## Summary metrics

In order to decide, which conditions resulted in the best response correction the **mean absoulte and relative pCO~2~ offset across all profiles** was calculated for:

* the offset between downcast and reference value
* the offset between downcast and upcast
* constant and T-dependent tau
* applied tau factors

Summary statistics were restricted to complete shallow profiles (not more than 2 observations missing from 1m depth intervals, maximum depth 20m). 



```{r identify_max_dep_pCO2_within_profiles}

tm_grid_stat <- tm_grid %>% 
  filter(correction == "pCO2") %>% 
  #drop_na() %>% 
  group_by(ID, station) %>% 
  summarise(dep_max = max(dep_grid),
            pCO2_max = max(down)) %>% 
  ungroup()

tm_grid_stat %>% 
  ggplot(aes(dep_max, pCO2_max))+
  geom_point()

tm_grid <- full_join(tm_grid, tm_grid_stat)
tm_grid_high_res <- full_join(tm_grid_high_res, tm_grid_stat)

meta <- read_csv(here::here("Data/TinaV/Sensor",
                            "Sensor_meta.csv"),
                 col_types = cols(ID = col_character()))

tm_grid_stat <- full_join(meta, tm_grid_stat)

tm_grid_stat  %>% 
 write_csv(here::here("data/_merged_data_files/response_time",
                      "tm_grid_stat_meta.csv"))

rm(tm_grid_stat, meta)

```


### Define dep and pCO2 thresholds

```{r dep_pCO2_threshold}

dep_max_lim <- 30
pCO2_max_lim <- 300
dep_lim <- 20

```

Summary metric are restricted to profiles that did not exceed a maximum sampling depth of `r dep_max_lim` m and pCO2 of `r pCO2_max_lim` µatm.

```{r calculate_performance_metrics}

tm_grid_shallow <- tm_grid %>% 
  filter(dep_max <= dep_max_lim,
         pCO2_max <= pCO2_max_lim,
         dep_grid <= dep_lim) %>% 
  group_by(ID, station, tau_factor, correction) %>% 
  mutate(nr_na = sum(is.na(pCO2_delta))) %>% 
  ungroup() %>% 
  filter(nr_na <= 2)

tm_grid_shallow_high_res <- tm_grid_high_res %>% 
  filter(dep_max <= dep_max_lim,
         pCO2_max <= pCO2_max_lim,
         dep_grid <= dep_lim) %>% 
  group_by(ID, station, tau_factor, correction) %>% 
  mutate(nr_na = sum(is.na(pCO2_delta))) %>% 
  ungroup() %>% 
  filter(nr_na <= 2)


tm_grid_shallow_sum <- tm_grid_shallow %>% 
  mutate(pCO2_delta_abs = abs(pCO2_delta),
         pCO2_delta_rel_abs = abs(pCO2_delta_rel)) %>% 
  group_by(tau_factor, dep_grid, correction) %>% 
  summarise(mean         = mean(pCO2_delta, na.rm = TRUE),
            sd           = sd(pCO2_delta, na.rm = TRUE),
            mean_abs     = mean(pCO2_delta_abs, na.rm = TRUE),
            mean_rel     = mean(pCO2_delta_rel, na.rm = TRUE),
            sd_rel       = sd(pCO2_delta_rel, na.rm = TRUE),
            mean_rel_abs = mean(pCO2_delta_rel_abs, na.rm = TRUE)) %>% 
  ungroup() %>% 
  pivot_longer(cols = sd:mean_rel_abs,
               names_to = "estimate", values_to = "pCO2_delta")

tm_grid_shallow_sum_high_res <- tm_grid_shallow_high_res %>% 
  mutate(pCO2_delta_abs = abs(pCO2_delta),
         pCO2_delta_rel_abs = abs(pCO2_delta_rel)) %>% 
  group_by(tau_factor, dep_grid, correction) %>% 
  summarise(mean         = mean(pCO2_delta, na.rm = TRUE),
            sd           = sd(pCO2_delta, na.rm = TRUE),
            mean_abs     = mean(pCO2_delta_abs, na.rm = TRUE),
            mean_rel     = mean(pCO2_delta_rel, na.rm = TRUE),
            sd_rel       = sd(pCO2_delta_rel, na.rm = TRUE),
            mean_rel_abs = mean(pCO2_delta_rel_abs, na.rm = TRUE)) %>% 
  ungroup() %>% 
  pivot_longer(cols = sd:mean_rel_abs,
               names_to = "estimate", values_to = "pCO2_delta")


tm_equi_delta_sum <- tm_equi_delta %>%
  mutate(pCO2_delta_equi_abs = abs(pCO2_delta_equi),
         pCO2_delta_equi_rel_abs = abs(pCO2_delta_equi_rel)) %>%
  group_by(correction, tau_factor) %>%
  summarise(mean         = mean(pCO2_delta_equi, na.rm = TRUE),
            mean_abs     = mean(pCO2_delta_equi_abs, na.rm = TRUE),
            mean_rel = mean(pCO2_delta_equi_rel, na.rm = TRUE),
            mean_rel_abs = mean(pCO2_delta_equi_rel_abs, na.rm = TRUE)) %>%
  ungroup() %>%
  pivot_longer(cols = mean:mean_rel_abs,
               names_to = "estimate", values_to = "dpCO2")

```


```{r plot_cast_delta_pCO2, fig.cap="Offset between up- and downcast. Panel columns: Constant and T-dependent tau. Panel rows from top to bottom: Mean of absolute offset, mean of relative absolute offset, standard deviation of offset, standard deviation of relative offset.", fig.asp=2}

tm_grid_shallow_sum %>% 
  filter(correction == "pCO2",
         estimate %in% c("mean_abs", "mean_rel_abs", "sd", "sd_rel")) %>% 
  ggplot()+
  geom_vline(xintercept = 0)+
  geom_hline(yintercept = 20)+
  geom_vline(xintercept = c(10), col="red")+
  geom_path(aes(pCO2_delta, dep_grid, col=as.factor(tau_factor)))+
  scale_y_reverse()+
  scale_color_discrete(name="Tau factor")+
  labs(x=expression(Delta~pCO[2]~(µatm)), y="Depth intervals (1m)")+
  facet_grid(estimate~.)

```

```{r plot_cast_delta_pCO2_example, include=FALSE}

tm_grid_shallow_sum_high_res %>% 
  filter(correction == "pCO2",
         estimate %in% c("mean_abs"),
         tau_factor == 1.25) %>% 
  ggplot()+
  geom_path(aes(pCO2_delta, dep_grid))+
  scale_y_reverse()+
  scale_color_discrete(name="Tau factor")+
  labs(x=expression(Mean~"|"~Delta~pCO[2]~"|"~(µatm)), y="Depth (m)")+
  xlim(0,NA)

ggsave(here::here("output/Plots/Figures_publication/appendix",
                  "pCO2_RT_corrected_profile_mean_abs_offset.pdf"),
       width = 100, height = 140, dpi = 300, units = "mm")

```


```{r plot_tau_optimum_cast_delta_pCO2, fig.cap="Mean offset between up- and downcast for all profiles up to 20m water depth. The lines between discrete tau factors result from the same analysis performed with high resolution of the tau factor. Left Panel: Mean absolute offset (µatm). Right panel: Mean relative offset (% of absolute value)."}


tm_grid_shallow_sum_mean_high_res <- tm_grid_shallow_sum_high_res %>% 
  group_by(estimate, correction, tau_factor) %>% 
  summarise(pCO2_delta_mean = mean(pCO2_delta)) %>% 
  ungroup()


ggplot()+
  geom_line(data = tm_grid_shallow_sum_mean_high_res %>% filter(estimate %in% c("mean_abs", "mean_rel_abs")),
             aes(tau_factor, pCO2_delta_mean, linetype=correction, shape=correction))+
  ylim(0,NA)+
  facet_wrap(~estimate)

```

Below we determined the tau factor that corresponds to lowest absolute mean offset.

```{r cast_tau_with_minimum_offfset}

pCO2_delta_mean_min <- tm_grid_shallow_sum_mean_high_res %>%
  filter(correction == "pCO2", estimate %in% c("mean_abs")) %>% 
  slice(which.min(pCO2_delta_mean)) %>% 
  select(pCO2_delta_mean) %>% 
  pull() %>% 
  round(2)

tau_factor_min <- tm_grid_shallow_sum_mean_high_res %>%
  filter(correction == "pCO2", estimate %in% c("mean_abs")) %>% 
  slice(which.min(pCO2_delta_mean)) %>% 
  select(tau_factor) %>% 
  pull()

```

A tau factor of `r tau_factor_min` resulted in the lowest absolute mean offset (`r pCO2_delta_mean_min` µatm) and will be used for response time correction.

Likewise, we analyse the offset from the pCO~2~ reference value:


```{r tau_optimum_reference_pCO2, fig.cap="Mean pCO~2~ offset from reference values as a function of the factor applied to tau. The lines between discrete tau factors result from the same analysis performed with high resolution of the tau factor. Left Panel: Mean absolute offset (µatm). Right panel: Mean relative offset (% of absolute value)."}

# tm_equi_delta_sum_highres <-
#  read_csv(here::here("data/_merged_data_files",
#                       "X_BloomSail_CTD_HydroC_profiles_RT_reference-offset_highres_taufactor_mean.csv")) %>% 
#   pivot_longer(cols = 4:7, names_to = "estimate", values_to = "dpCO2") %>% 
#   filter(estimate %in% c("mean_abs", "mean_rel_abs"))

tm_equi_delta_sum %>%
  filter(estimate %in% c("mean_abs", "mean_rel_abs")) %>% 
  ggplot(aes(tau_factor, dpCO2, linetype=correction, shape=correction))+
  # geom_line(data = tm_equi_delta_sum_highres,
  #           aes(tau_factor, dpCO2))+
  geom_point()+
  geom_hline(yintercept = 0)+
  labs(x="Tau factor", y=expression(Mean~Delta~pCO[2]))+
  facet_wrap(~estimate)


```

## Summary plots


```{r plot_profiling_timeseries_all, eval=FALSE}

i_tau_factor <- "1.2"

cast_dep <- tm %>% 
  filter(tau_factor == i_tau_factor) %>% 
  pivot_longer(c(dep, pCO2_corr, pCO2),
               names_to = "parameter", values_to = "value")

cast_dep_equi <- tm_pCO2_equi %>% 
  pivot_longer(c(dep, pCO2_corr), names_to = "parameter", values_to = "value")

tm_sub <- tm %>% 
  filter(tau_factor == i_tau_factor)

tm_grid_sub <- tm_grid %>% 
  filter(tau_factor == i_tau_factor)

max_duration <- round(max(cast_dep$duration)/1000,0)*1000

pdf(file=here::here("output/Plots/response_time",
                    "all_plots.pdf"),
    onefile = TRUE, width = 7, height = 10)

for(i_ID in unique(tm$ID)){
  for(i_station in unique(tm$station)){

    # i_ID <- unique(cast_dep$ID)[1]
    # i_station <- unique(cast_dep$station)[1]
    i_ID
    i_station
    
    if (nrow(cast_dep %>% filter(ID == i_ID, station == i_station)) > 0){
          if (nrow(tm %>% filter(ID == i_ID, station == i_station)) > 0){
     if (nrow(tm_grid %>% filter(ID == i_ID, station == i_station)) > 0){
      
      cast_dep_equi_sub <- cast_dep_equi %>%
        filter(ID == i_ID,
               station == i_station)

      p_time_series <- cast_dep %>%
        filter(ID == i_ID,
               station == i_station) %>%
        ggplot(aes(duration, value, col=phase))+
        geom_point(size=0.5)+
        geom_point(data = cast_dep_equi_sub, aes(duration, value), col="black")+
        scale_y_reverse()+
        scale_x_continuous(breaks = seq(0,6000,500))+
        labs(title = str_c("ID: ",i_ID," | Station: ",i_station))+
        facet_grid(parameter~., scales = "free_y")+
        theme_bw()

      p_profile <- tm_sub %>%
        filter(ID == i_ID,
               station == i_station,
               phase %in% c("up", "down")) %>%
        ggplot()+
        geom_path(aes(pCO2_corr, dep, linetype = phase, col="raw"))+
        geom_path(aes(pCO2, dep, linetype = phase, col="corrected"))+
        #geom_point(data = equi_cast, aes(pCO2_corr, dep))+
        scale_y_reverse()+
        coord_cartesian(ylim = c(25,0), xlim = c(70,250))+
        scale_color_brewer(palette = "Set1", name="", guide = FALSE)+
        scale_linetype(guide = FALSE)+
        labs(y="Depth [m]", x="pCO2",
             title = "full res")

      
      tm_pCO2_equi_grid_sub <- tm_pCO2_equi_grid %>%
        filter(ID == i_ID,
               station == i_station)

      p_profile_grid <- tm_grid_sub %>%
        filter(ID == i_ID,
               station == i_station) %>%
        arrange(dep_grid) %>%
        ggplot()+
        geom_path(aes(down, dep_grid, col=correction, linetype="down"))+
        geom_path(aes(up, dep_grid, col=correction, linetype ="up"))+
        geom_point(data = tm_pCO2_equi_grid_sub, aes(pCO2_equi, dep_grid))+
        scale_y_reverse()+
        coord_cartesian(ylim = c(25,0), xlim = c(70,250))+
        scale_linetype(name="cast")+
        scale_color_brewer(palette = "Set1", direction = -1)+
        labs(x="pCO2",
             title = "1m grid")+
        theme(axis.title.y = element_blank(),
              axis.text.y = element_blank())


      p_delta_abs <- tm_grid %>%
        filter(ID == i_ID,
               station == i_station,
               correction == "pCO2") %>%
        arrange(dep_grid) %>%
        ggplot(aes(pCO2_delta, dep_grid, col=as.factor(tau_factor)))+
        geom_path()+
        geom_point()+
        scale_y_reverse(breaks=seq(0,40,2))+
        scale_color_discrete(name="tau factor", guide=FALSE)+
        labs(x = "delta pCO2 [µatm]", y="Depth [m]")+
        geom_vline(xintercept = 0)+
        geom_vline(xintercept = c(-10,10), col="red")
      
      p_delta_rel <- tm_grid %>%
        filter(ID == i_ID,
               station == i_station,
               correction == "pCO2") %>%
        arrange(dep_grid) %>%
        ggplot(aes(pCO2_delta_rel, dep_grid, col=as.factor(tau_factor)))+
        geom_path()+
        geom_point()+
        scale_y_reverse(breaks=seq(0,40,2))+
        scale_color_discrete(name="tau factor")+
        labs(x = "delta pCO2 [% of absolute value]", y="Depth [m]")+
        geom_vline(xintercept = 0)+
        geom_vline(xintercept = c(-10,10), col="red")+
        theme(axis.title.y = element_blank(),
              axis.text.y = element_blank())
      
      print(
  
  p_time_series / (p_profile | p_profile_grid) / (p_delta_abs | p_delta_rel)
  
)
    rm(p_time_series, p_profile, p_profile_grid, p_delta_abs, p_delta_rel)

     }}}



  }
}

dev.off()

rm(cast_dep, cast_dep_equi, cast_dep_equi_sub, i_station, i_ID,
   max_duration, tm_pCO2_equi_grid_sub)


```



# Conclusion

* Taking the temperature dependence of tau into account resulted in a slightly better agreement between up- and downcast, as well as downcast and reference value
* For all quality metrics we find improved agreement for slightly positive tau factor ranging from 1.04 - 1.24
* As the focus is on pCO~2~ changes in the upper 10 m of the water column, we will apply the tau factor 1.24 leading to the lowest relative offset between up- and downcast

# Correct entire data set

Finally, the response time correction was applied to the full data set (not only profile data) based on the optimum parameterization determined above.


```{r RT_correction_all, eval=FALSE}


# Response time correction approach after Bittig --------------------------
RT_corr <- function(c1, c0, dt, tau) {
  ( 1 / ( 2* (( 1+(2*tau/dt) )^(-1) ))) * (c1 - (1-(2* (( 1+(2*tau/dt) )^(-1) ))) * c0)
}



tm <- read_csv(here::here("data/_merged_data_files/merging_interpolation",
                          "tm.csv"),
               col_types = cols(ID = col_character(),
                                pCO2_analog = col_double(),
                                pCO2_corr = col_double(),
                                Zero = col_factor(),
                                Flush = col_factor(),
                                Zero_counter = col_integer(),
                                deployment = col_integer(),
                                duration = col_double(),
                                mixing = col_character(),
                                lat = col_double(),
                                lon = col_double()))

tm <- tm %>%
  group_by(ID, station) %>% 
  mutate(duration = as.numeric(date_time - min(date_time)),
         pump_power = if_else(date_time < ymd_hms("2018-07-17;13:08:34"), "1W", "8W")) %>%
  arrange(date_time)


tau_fit <- read_csv(here::here("data/_merged_data_files/response_time",
                              "tau_fit.csv"))



# Assign T-dependent response time (tau) values ----------------------------------------------

tau_fit <- tau_fit %>% 
  rename(tau_intercept = `(Intercept)`, tau_slope=tem)

tm <- full_join(tm, tau_fit)


#Assign tau
tm <- tm %>% 
  mutate(tau = tau_intercept + tau_slope *tem) %>% 
  select(-tau_intercept, -tau_slope)


# Prepare data set for RT correction --------------------------------------

tm <- tm %>% 
  group_by(ID, station) %>% 
  arrange(date_time) %>% 
  mutate(dt = as.numeric(as.character(date_time - lag(date_time)))) %>% 
  ungroup()

# measurement frequency

freq <- tm %>% 
  filter(dt < 13) %>% 
  group_by(ID) %>% 
  summarise(dt_mean = round(mean(dt, na.rm = TRUE),0))

tm <- full_join(tm, freq)

# tau factors
tm <- expand_grid(tm, tau_factor = tau_factor_min)

tm <- tm %>%
  mutate(tau_test = tau*tau_factor)


# Apply RT correction to entire data set

for(i_ID in unique(tm$ID)){
  
#i_ID <- "180716"
  
  freq_sub <- freq %>% filter(ID == i_ID) %>% pull(dt_mean)
  
  window <- 30 / freq_sub
  rolling_mean   <- rollify(~mean(.x, na.rm = TRUE), window = window)

  tm_sub <- tm %>%
    filter(ID == i_ID) %>% 
    group_by(station) %>% 
    mutate(pCO2_RT = RT_corr(pCO2_corr, lag(pCO2_corr), dt, tau_test),
           pCO2_RT = if_else(pCO2_RT %in% c(Inf, -Inf), NaN, pCO2_RT),
           window = window,
           pCO2 = rolling_mean(pCO2_RT)
           ) %>% 
    ungroup()
  
  # time shift RT corrected data
  shift <- as.integer(as.character(window/2))
  
  tm_sub <- tm_sub %>%
    group_by(station) %>% 
    mutate(pCO2 = lead(pCO2, shift)) %>% 
    ungroup()
  
  
  if (exists("tm_corr")){tm_corr <- bind_rows(tm_corr, tm_sub)}
  else{tm_corr <- tm_sub}
  
  rm(tm_sub, freq_sub, rolling_mean, shift, window)
  
}


rm(RT_corr, i_ID, freq)

tm_corr <- tm_corr %>% 
  select(-c(tau, tau_factor, tau_test, window))

tm_corr %>% 
 write_csv(here::here("data/_merged_data_files/response_time",
                      "tm_RT_all.csv"))
rm(tm, tm_corr)

```


# Sensitivity considerations

A change in DIC of 1 µmol kg^-1^ corresponds to a change in pCO~2~ of around 1 µatm, in the Central Baltic Sea at a pCO~2~ of around 100 µatm (summertime conditions).

```{r sensitivity_estimate, fig.cap="pCO~2~ sensitivity to changes in DIC."}

df <- data.frame(cbind(
  (c(1720)),
  (c(7))))

Tem <- seq(5,25,5)
pCO2<-seq(50,500,20)

df<-merge(df, Tem)
names(df) <- c("AT", "S", "Tem")  

df<-merge(df, pCO2)
names(df) <- c("AT", "S", "Tem", "pCO2")  

df$AT<-df$AT*1e-6

df$DIC<-carb(flag=24, var1=df$pCO2, var2=df$AT, S=df$S, T=df$Tem, k1k2="m10", kf="dg", pHscale="T")[,16]
df$pCO2.corr<-carb(flag=15, var1=df$AT, var2=df$DIC, S=df$S, T=df$Tem, k1k2="m10", kf="dg", pHscale="T")[,9]

df$pCO2.2<-df$pCO2.corr + 25
df$DIC.2<-carb(flag=24, var1=df$pCO2.2, var2=df$AT, S=df$S, T=df$Tem, k1k2="m10", kf="dg", pHscale="T")[,16]


df$ratio<-(df$pCO2.2-df$pCO2.corr)/(df$DIC.2*1e6-df$DIC*1e6)

df %>% 
  ggplot(aes(pCO2, ratio, col=as.factor(Tem)))+
  geom_line()+
  scale_color_viridis_d(option = "C",name="Tem [°C]")+
  labs(x=expression(pCO[2]~(µatm)), y=expression(Delta~pCO[2]~"/"~Delta~DIC~(µatm~µmol^{-1}~kg)))+
  scale_y_continuous(limits = c(0,8), breaks = seq(0,10,1))

rm(df, Tem, pCO2)

```




# Open tasks / questions

- Clean workspace for plotting code after RT correction