Finnmaid+GETM.Rmd

---
title: "Finnmaid and GETM"
author: "Jens Daniel Müller und Lara Sophie Burchardt"
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(ncdf4)
library(vroom)
library(lubridate)
library(here)
library(seacarb)
library(oce)
library(patchwork)
library(zoo)
```

```{r subsetting_criteria}

# route
select_route <- "E"

# latitude limits
low_lat <- 57.3
high_lat <- 57.5

#depth range to subset GETM 3d files
d1_shallow <- 0
d1_deep <- 30

# date limits
start_date <- "2018-06-20"
end_date <- "2018-08-25"

fixed_values <- 
  read_csv(here::here("Data/_summarized_data_files", "BloomSail_fixed_values.csv"))


```

# GETM Data preparation

## Salinity and temperature profiles

Mean salinity and temperature profiles within the BloomSail area were extracted from daily GETM transects beneath the Finnmaid track.

```{r read_getm_3d_sal_tem, eval=FALSE}

nc <- nc_open(here::here("data/GETM", "Finnmaid.E.3d.2018.nc"))

lat <- ncvar_get(nc, "latc")

time_units <- nc$dim$time$units %>%     #we read the time unit from the netcdf file to calibrate the time 
    substr(start = 15, stop = 33) %>%   #calculation, we take the relevant information from the string
    ymd_hms()                           # and transform it to the right format
t <- time_units + ncvar_get(nc, "time") # read time vector
rm(time_units)

d <- ncvar_get(nc, "zax") # read depths vector

for (var_3d in c("salt", "temp")) {
  
array <- ncvar_get(nc, var_3d) # store the data in a 3-dimensional array
#dim(array) # should be 3d with dimensions: 544 coordinates, 51 depths, and number of days of month

fillvalue <- ncatt_get(nc, var_3d, "_FillValue")

# Working with the data
array[array == fillvalue$value] <- NA

    for (i in seq(1,length(t),1)){
      
      # i <- 3
      array_slice <- array[, , i] # slices data from one day
      
      array_slice_df <- as.data.frame(t(array_slice))
      array_slice_df <- as_tibble(array_slice_df)
      
      gt_3d_part <- array_slice_df %>%
        set_names(as.character(lat)) %>%
        mutate(dep = -d) %>%
        gather("lat", "value", 1:length(lat)) %>%
        mutate(lat = as.numeric(lat)) %>%
        filter(lat > low_lat, lat < high_lat,
               dep >= d1_shallow, dep <= d1_deep) %>%
        #summarise_all("mean") %>%
        mutate(var = var_3d,
               date_time=t[i]) %>% 
        dplyr::select(date_time, dep, value, var) #%>% 
        #filter(date_time >= start_date, date_time <= end_date)

      
      if (exists("gt_3d")) {
        gt_3d <- bind_rows(gt_3d, gt_3d_part)
        } else {gt_3d <- gt_3d_part}
      
  rm(array_slice, array_slice_df, gt_3d_part)
      
    }
rm(array, fillvalue)

}

nc_close(nc)
rm(nc)

gt_3d_jun_aug <- gt_3d %>% 
  filter(date_time >= start_date & date_time <= end_date) %>% 
  group_by(dep,var,date_time ) %>% 
  summarise_all(list(value=~mean(.,na.rm=TRUE))) %>%
  ungroup()

gt_3d_jun_aug %>% 
  vroom_write((here::here("data/_summarized_data_files", file = "gt_3d_jun_aug_2018.csv")))


rm(gt_3d, gt_3d_jun_aug, d1_deep, d1_shallow, i, lat, d, t, var_3d)

```


## Mixed layer depth

Regional mean mixed layer depth estimates based on sewater density and windspeed within the BloomSail area were extracted from 3h GETM surface data along the Finnmaid track.

```{r read_getm_2d_mld_wind, eval=FALSE}

nc_2d <- nc_open(here("data/GETM", "Finnmaid.E.2d.2018.nc"))
#print(nc_2d)

lat <- ncvar_get(nc_2d, "latc")

time_units <- nc_2d$dim$time$units %>%     #we read the time unit from the netcdf file to calibrate the time 
    substr(start = 15, stop = 33) %>%   #calculation, we take the relevant information from the string
    ymd_hms()                           # and transform it to the right format
t <- time_units + ncvar_get(nc_2d, "time") # read time vector
rm(time_units)

for (var in names(nc_2d$var)[c(3,4,6:12)]) {
  
#var <- "mld_rho"

array <- ncvar_get(nc_2d, var) # store the data in a 3-dimensional array
fillvalue <- ncatt_get(nc_2d, var, "_FillValue")
array[array == fillvalue$value] <- NA

array <- as.data.frame(t(array), xy=TRUE)
array <- as_tibble(array)
      
  gt_2d_part <- array %>%
  set_names(as.character(lat)) %>%
  mutate(date_time = t) %>%
  filter(date_time >= start_date & date_time <= end_date) %>% 
  gather("lat", "value", 1:length(lat)) %>%
  mutate(lat = as.numeric(lat)) %>%
  filter(lat > low_lat, lat<high_lat) %>%
  select(-lat) %>% 
  group_by(date_time) %>% 
  summarise_all(list(value=~mean(.,na.rm=TRUE))) %>%
  ungroup() %>% 
  mutate(var = var)
     
  if (exists("gt_2d")) {
    gt_2d <- bind_rows(gt_2d, gt_2d_part)
    } else {gt_2d <- gt_2d_part} 

rm(array, fillvalue, gt_2d_part)

}

nc_close(nc_2d)
rm(nc_2d)

gt_2d <- gt_2d %>% 
  mutate(value = round(value, 3)) %>% 
  pivot_wider(values_from = value, names_from = var) %>% 
  mutate(U_10 = round(sqrt(u10^2 + v10^2), 3)) %>% 
  select(-c(u10, v10))


gt_2d %>% 
  vroom_write((here::here("data/_summarized_data_files", file = "gt_2d_jun_aug_2018.csv")))

rm(t, var, gt_2d, lat)

```

# Hydrography

## Read summary files

```{r read_gt_summary_files}

gt_3d_jun_aug <- 
  read_tsv((here::here("data/_summarized_data_files", file = "gt_3d_jun_aug_2018.csv")))

gt_2d_jun_aug <- 
  read_tsv((here::here("data/_summarized_data_files", file = "gt_2d_jun_aug_2018.csv")))

gt_3d_jun_aug <- gt_3d_jun_aug %>% 
  pivot_wider(values_from = value, names_from = var) %>% 
  mutate(rho = swSigma(salinity = salt, temperature = temp, pressure = dep/10))

gt_2d_jun_aug_daily <- gt_2d_jun_aug %>% 
  mutate(day = yday(date_time)) %>% 
  group_by(day) %>% 
  summarise_all(list(~mean(.,na.rm=TRUE))) %>%
  ungroup() %>% 
  select(-day)

```


## Hovmoeller Plots

```{r plot_hovmoeller_gt_salt_temp_rho}

p_sal <- gt_3d_jun_aug %>% 
  ggplot()+
  geom_raster(aes(date_time, dep, fill=salt))+
  geom_vline(data=fixed_values, aes(xintercept = start))+
  geom_vline(data=fixed_values, aes(xintercept = end))+
  scale_fill_viridis_c(name="Salinity ", direction = -1)+
  scale_y_reverse()+
  coord_cartesian(expand = 0)+
  labs(y="Depth [m]")+
  theme_bw()+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())+
  geom_line(data= gt_2d_jun_aug_daily,
            aes(x = date_time, y = mld_rho), color = "white")+
  scale_color_discrete(name = "Legend", labels = c("MLD Rho"))


p_tem <- gt_3d_jun_aug %>% 
  ggplot()+
  geom_raster(aes(date_time, dep, fill=temp))+
  geom_vline(data=fixed_values, aes(xintercept = start))+
  geom_vline(data=fixed_values, aes(xintercept = end))+
  scale_fill_viridis_c(name="Temperature (°C)", option = "B")+
  scale_y_reverse()+
  coord_cartesian(expand = 0)+
  labs(y="Depth [m]", x="")+
  theme_bw()+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())+
  geom_line(data= gt_2d_jun_aug_daily,
            aes(x = date_time, y = mld_rho), color = "white")+
  scale_color_discrete(name = "Legend", labels = c("MLD Rho"))

p_rho <- gt_3d_jun_aug %>% 
  ggplot()+
  geom_raster(aes(date_time, dep, fill=rho))+
  geom_vline(data=fixed_values, aes(xintercept = start))+
  geom_vline(data=fixed_values, aes(xintercept = end))+
  scale_fill_viridis_c(name="d Rho (kg/m^3)", direction = -1)+
  scale_y_reverse()+
  coord_cartesian(expand = 0)+
  labs(y="Depth [m]", x="")+
  theme_bw()+
  theme(axis.title.x = element_blank())+
  geom_line(data= gt_2d_jun_aug_daily,
            aes(x = date_time, y = mld_rho), color = "white")+
  scale_color_discrete(name = "Legend", labels = c("MLD Rho"))

p_sal / p_tem / p_rho

rm(p_sal, p_tem, p_rho)

```

## Profiles

```{r plot_profiles_gt_salt_temp_rho}

gt_3d_jun_aug_long <- gt_3d_jun_aug %>% 
  pivot_longer(3:5, values_to = "value", names_to = "parameter")
#
lab_dates <- pretty(gt_3d_jun_aug_long$date_time)

gt_3d_jun_aug_long %>% 
  ggplot(aes(value, dep,
             col=as.numeric(date_time),
             group=as.factor(date_time)))+
  geom_path()+
  scale_y_reverse(expand = c(0,0))+
  scale_color_viridis_c(breaks = as.numeric(lab_dates),
                        labels = lab_dates,
                        name="Date")+
  theme_bw()+
  facet_grid(~parameter, scales = "free_x")

```


# Metrology

## Windspeeds

```{r plot_winspeed}

gt_2d_jun_aug %>% 
  ggplot()+
  geom_rect(data = fixed_values, aes(xmin=start, xmax=end, ymin=-Inf, ymax=Inf), alpha=0.2)+
  geom_line(aes(x= date_time, y = U_10, col="3-hourly"))+
  geom_line(data = gt_2d_jun_aug_daily,
            aes(x= date_time, y = U_10, col="Daily mean"))+
  labs(y="U (m/s)", x = "Date")+
  scale_color_brewer(palette = "Set1", name="", direction = -1)+
  theme_bw()

```



# Finnmaid

## Read data

Finnmaid data, including reconstructed data during LICOS operation failure.

```{r read_Finnmaid}

df <-
 read_csv(here::here("Data/_summarized_data_files",
                      "Finnmaid.csv"))

df <- df %>% 
  filter(Area == "BS",
         date > start_date,
         date < end_date) %>% 
  select(-c(Lon, Lat, patm, Teq, xCO2, route, Area)) %>% 
  mutate(ID = as.factor(ID)) %>% 
  rename(tem=Tem,
         sal=Sal)

```

## C~T~ calculation

Calculate based on fixed A~T~ and salinity mean values. 

```{r CT_calculation}

df <- df %>% 
  mutate(CT = carb(24, var1=pCO2, var2=1720*1e-6,
                   S=6.9, T=tem, k1k2="m10", kf="dg", ks="d",
                   gas="insitu")[,16]*1e6)

```

## Timeseries

```{r timeseries_of_regional_averages}

FM_ts_mean <- df %>% 
  arrange(date) %>% 
  group_by(ID, sensor) %>% 
  summarise_all(list(~mean(.)), na.rm=TRUE) %>%
  ungroup() %>% 
  pivot_longer(4:8, values_to = "mean", names_to = "parameter")

FM_ts_sd <- df %>% 
  arrange(date) %>% 
  group_by(ID, sensor) %>% 
  summarise_all(list(~sd(.)), na.rm=TRUE) %>%
  ungroup() %>% 
  pivot_longer(4:8, values_to = "sd", names_to = "parameter") %>% 
  select(-date)

FM_ts <- inner_join(FM_ts_mean, FM_ts_sd)

rm(FM_ts_mean, FM_ts_sd)

```


```{r Plot_Finnmaid_crossings, fig.asp=1.5}

FM_ts %>% 
  ggplot()+
  geom_rect(data = fixed_values, aes(xmin=start, xmax=end, ymin=-Inf, ymax=Inf), alpha=0.2)+
  geom_path(aes(x=date, y=mean, ymax=mean+sd, ymin=mean-sd))+
  geom_ribbon(aes(x=date, y=mean, ymax=mean+sd, ymin=mean-sd), alpha=0.5)+
  geom_point(aes(x=date, y=mean, ymax=mean+sd, ymin=mean-sd, col=sensor))+
  facet_grid(parameter~., scales = "free_y")+
  scale_color_brewer(palette = "Set1", direction=-1)+
  theme_bw()

```


# C~T~ vs temperature

## SST time series

```{r plot_FM_gt_temp_timeseries}

FM_ts %>% 
  filter(parameter == "tem") %>% 
  ggplot()+
  geom_rect(data = fixed_values, aes(xmin=start, xmax=end, ymin=-Inf, ymax=Inf), alpha=0.2)+
  geom_path(aes(x=date, y=mean, ymax=mean+sd, ymin=mean-sd))+
  geom_ribbon(aes(x=date, y=mean, ymax=mean+sd, ymin=mean-sd), alpha=0.5)+
  geom_point(aes(x=date, y=mean, col="Finnmaid"))+
  geom_path(data = gt_2d_jun_aug, aes(x=date_time, y=SST, col="GETM, 3h"))+
  geom_path(data = gt_2d_jun_aug_daily, aes(x=date_time, y=SST, col="GETM, daily"))+
  scale_color_brewer(palette = "Set1", name="")+
  labs(x="", y="SST (°C)")+
  theme_bw()

```

## dC~T~ vs dSST

As primary production (negative changes in C~T~) and increase in seawater temperature have a common driver (light), the relation between both changes was investigated.

The following analysis is restricted to the BloomSail period. 

```{r CT_vs_temperature}

CT_tem <- FM_ts %>%
  filter(date > fixed_values$start,
         date < fixed_values$end) %>% 
  filter(parameter %in% c("tem", "CT")) %>% 
  select(date, parameter, mean) %>% 
  pivot_wider(values_from = mean, names_from = parameter) %>% 
  set_names(c("date_time", "FM", "CT"))

CT_tem <- full_join(CT_tem,
                    gt_2d_jun_aug_daily %>% select(date_time, gt = SST))

CT_tem <- CT_tem %>% 
  arrange(date_time) %>%
  mutate(gt = na.approx(gt)) %>% 
  drop_na()

CT_tem <- CT_tem %>% 
  pivot_longer(cols = c(FM, gt), values_to = "SST", names_to = "obs")

CT_tem <- CT_tem %>% 
  group_by(obs) %>% 
  arrange(date_time) %>% 
  mutate(dCT = CT - lag(CT),
         dSST = SST - lag(SST),
         sign = if_else(dCT < 0, "neg", "pos")) %>% 
  ungroup()


lab_dates <- pretty(CT_tem$date_time)

CT_tem %>% 
  ggplot()+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0)+
  geom_smooth(aes(dSST, dCT), method = "lm", se=FALSE)+
  geom_point(aes(dSST, dCT, fill=as.numeric(date_time)), shape=21)+
  scale_fill_viridis_c(breaks = as.numeric(lab_dates),
                        labels = lab_dates,
                        name="Date")+
  guides(fill = guide_legend(override.aes=list(shape=21)))+
  facet_wrap(~obs)


```

# NCP

## MLD approach

Use dCT/dtem from Finnmaid, MLD from GETM

```{r NCP_MLD}

NCP_MLD <- FM_ts %>%
  filter(date > fixed_values$start,
         date < fixed_values$end) %>% 
  filter(parameter %in% c("tem", "CT")) %>% 
  select(date, parameter, mean) %>% 
  pivot_wider(values_from = mean, names_from = parameter) %>% 
  set_names(c("date_time", "SST", "CT"))

NCP_MLD <- full_join(NCP_MLD,
                     gt_2d_jun_aug_daily %>% select(-c(SSS, SST, U_10)))

NCP_MLD <- NCP_MLD %>% 
  pivot_longer(cols = 4:8, values_to = "MLD", names_to = "Parameter") %>% 
  arrange(date_time) %>%
  group_by(Parameter) %>% 
  mutate(MLD = na.approx(MLD)) %>%
  ungroup()

NCP_MLD %>% 
  ggplot()+
  geom_rect(data = fixed_values, aes(xmin=start, xmax=end, ymin=-Inf, ymax=Inf), alpha=0.2)+
  geom_line(aes(date_time, MLD, col=Parameter))+
  geom_point(data = NCP_MLD %>% drop_na, aes(date_time, MLD, col=Parameter))+
  scale_color_brewer(palette = "Set1")+
  labs(x="", y="MLD (m)")+
  theme_bw()


NCP_MLD <- NCP_MLD %>% 
  drop_na() %>% 
  group_by(Parameter) %>% 
  arrange(date_time) %>% 
  mutate(dCT = CT - lag(CT),
         dSST = SST - lag(SST),
         dNCP = dCT * MLD / 1000,
         NCP_cum = cumsum(replace_na(dNCP, 0))) %>% 
  ungroup()


```


### Incremental and cumulative timeseries

Total incremental and cumulative C~T~ changes inbetween cruise dates were calculated.


```{r plot_integrated_NCP_timeseries, fig.asp=1.5}

p_iNCP <- NCP_MLD %>% 
  ggplot(aes(date_time, dNCP, fill= Parameter))+
  geom_hline(yintercept = 0)+
  geom_col(col="black", position = "dodge")+
  scale_y_continuous(breaks = seq(-100, 100, 0.2))+
  labs(y="integrated, directional CT changes [mol/m2]", x="date")

p_iNCPcum <-  NCP_MLD %>% 
  ggplot(aes(date_time,  NCP_cum, 
             fill= Parameter))+
  geom_hline(yintercept = -0.4)+
  geom_area(position = "identity", col="black")+
  geom_hline(yintercept = 0)+
  scale_y_continuous(breaks = seq(-100, 100, 0.2))+
  facet_grid(Parameter~., scales = "free_y", space = "free_y")+
  theme(strip.background = element_blank(),
        strip.text = element_blank())+
  labs(y="integrated, cumulative, directional CT changes [mol/m2]", x="date")


(p_iNCP / p_iNCPcum)+
  plot_layout(guides = 'collect', heights = c(1, 2))

rm(p_iNCP, p_iNCPcum)
```


## delta T appraoch

```{r reconctruction_CT_incremental_changes_profiles}

CT_tem <- CT_tem %>% 
  filter(obs == "FM") %>% 
  mutate(factor = dCT/dSST,
         factor = if_else(is.na(factor), 0, factor)) %>% 
  select(date_time, factor)

CT_tem <- expand_grid(CT_tem, dep = unique(gt_3d_jun_aug$dep))

NCP_dT <- full_join(gt_3d_jun_aug %>% select(date_time, dep, temp), 
                    CT_tem)

NCP_dT <- NCP_dT %>% 
  arrange(dep, date_time) %>% 
  group_by(dep) %>% 
  mutate(temp = na.approx(temp)) %>% 
  ungroup() %>% 
  drop_na()

NCP_dT <- NCP_dT %>% 
  group_by(dep) %>% 
  arrange(date_time) %>% 
  mutate(dtem = temp - lag(temp)) %>% 
  ungroup() %>% 
  mutate(diff_value = dtem * factor) %>% 
  select(-factor)

```

The reconstructed incremental changes are added up to derive cummulative C~T~ changes throughout the water column.

```{r reconctruction_CT_cumulative_changes_profiles}

NCP_dT <- NCP_dT %>% 
  group_by(dep) %>% 
  arrange(date_time) %>% 
  mutate(diff_time  = as.numeric(date_time - lag(date_time)),
         diff_value_daily = diff_value / diff_time,
         cum_value = replace_na(diff_value, 0))

```

## Profiles of incremental changes

Changes of seawater parameters at each depth were reconstructed from one cruise day to the next and divided by the number of days inbetween.

```{r incremental_changes_profiles}

NCP_dT %>% 
  arrange(dep) %>% 
  ggplot(aes(diff_value_daily, dep, col=as.factor(date_time)))+
  geom_vline(xintercept = 0)+
  geom_point()+
  geom_path()+
  scale_y_reverse()+
  #scale_color_viridis_d()+
  #facet_wrap(~parameter, scales = "free_x")+
  labs(x="Change of value inbetween cruises per day")


```



## Profiles of cumulative changes

Cumulative changes of seawater parameters were calculated at each depth relative to the first cruise day on July 5.

```{r cumulative_changes_profiles}

NCP_dT %>% 
  arrange(dep) %>% 
  ggplot(aes(cum_value, dep, col=as.factor(date_time)))+
  geom_vline(xintercept = 0)+
  geom_point()+
  geom_path()+
  scale_y_reverse()+
  #scale_color_viridis_d()+
  labs(x="Cumulative change of value")

```

Reconstructed cumulative positive and negative changes of seawater parameters were calculated separately at each depth relative to the first cruise day on July 5.

```{r cumulative_directional_changes_profiles}

NCP_dT <- NCP_dT %>% 
  mutate(sign = if_else(diff_value < 0, "neg", "pos")) %>% 
  group_by(dep, sign) %>%
  arrange(date_time) %>%
  mutate(cum_value_sign = replace_na(diff_value, 0)) %>% 
  ungroup()

NCP_dT %>% 
  arrange(dep) %>% 
  ggplot(aes(cum_value_sign, dep, col=as.factor(date_time)))+
  geom_vline(xintercept = 0)+
  geom_point()+
  geom_path()+
  scale_y_reverse()+
  scale_color_viridis_d()+
  scale_fill_viridis_d()+
  facet_wrap(~sign, scales = "free_x", ncol=4)+
  labs(x="Cumulative directional change of value")

```