---
title: "Tidal analysis"
author: "Evert de Froe"
date: "1/8/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = TRUE)
list.of.packages <- c("plot3D", "OceanView", "reshape2","RnetCDF", "oce", "lubridate", "stringi", "ggplot2", "patchwork", "factoextra", "ggfortify", "cluster", "xts", "dplyr", "Rfast")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
#if(length(new.packages)) install.packages(new.packages)
lapply(list.of.packages, require, character.only = TRUE)



```

## Aim of document

In this file I want to retrieve the slack and high tide for the four CTD sampling locations.

First I load the data, which are in .out format.


```{r load data}
wwidt <- c(4,7,2,8,5,19,3,7,3,7,3,7,4,6,2,8)



# Bank
bank_k1 <- read.fwf("Data/Tides_second/station6k1.out", widths = wwidt)[-c(1,3,5,7,9,11,13,15)]
names(bank_k1) <- c("lat", "lon", "time", "k1_u_tran", "k1_v_tran", "k1_u", "k1_v", "depth")

bank_m2 <- read.fwf("Data/Tides_second/station6m2.out", widths = wwidt)[-c(1,3,5,7,9,11,13,15)]
names(bank_m2) <- c("lat", "lon", "time", "m2_u_tran", "m2_v_tran", "m2_u", "m2_v", "depth")

# Upper slope
usl_k1 <- read.fwf("Data/Tides_second/station5k1.out", widths = wwidt)[-c(1,3,5,7,9,11,13,15)]
names(usl_k1) <-c("lat", "lon", "time", "k1_u_tran", "k1_v_tran", "k1_u", "k1_v", "depth")

usl_m2 <- read.fwf("Data/Tides_second/station5m2.out", widths = wwidt)[-c(1,3,5,7,9,11,13,15)]
names(usl_m2) <- c("lat", "lon", "time", "m2_u_tran", "m2_v_tran", "m2_u", "m2_v", "depth")

# Lower slope
lsl_k1 <- read.fwf("Data/Tides_second/station4k1.out", widths = wwidt)[-c(1,3,5,7,9,11,13,15)]
names(lsl_k1) <- c("lat", "lon", "time", "k1_u_tran", "k1_v_tran", "k1_u", "k1_v", "depth")

lsl_m2 <- read.fwf("Data/Tides_second/station4m2.out", widths = wwidt)[-c(1,3,5,7,9,11,13,15)]
names(lsl_m2) <- c("lat", "lon", "time", "m2_u_tran", "m2_v_tran", "m2_u", "m2_v", "depth")


# Oreo
oreo_k1 <- read.fwf("Data/Tides_second/station2k1.out", widths = wwidt)[-c(1,3,5,7,9,11,13,15)]
names(oreo_k1) <- c("lat", "lon", "time", "k1_u_tran", "k1_v_tran", "k1_u", "k1_v", "depth")

oreo_m2 <- read.fwf("Data/Tides_second/station2m2.out", widths = wwidt)[-c(1,3,5,7,9,11,13,15)]
names(oreo_m2) <- c("lat", "lon", "time", "m2_u_tran", "m2_v_tran", "m2_u", "m2_v", "depth")

```

# Combine tides 

Now I combine the K1 and M2 tide in one dataset


```{r combine_tides}
bank <- inner_join(bank_k1, bank_m2, by = c("lat", "lon", "time", "depth"))
usl <- inner_join(usl_k1, usl_m2, by = c("lat", "lon", "time", "depth"))
lsl <- inner_join(lsl_k1, lsl_m2, by = c("lat", "lon", "time", "depth"))
oreo <- inner_join(oreo_k1, oreo_m2, by = c("lat", "lon", "time", "depth"))


bank$time <- as.POSIXct(bank$time,format='%m.%d.%Y %H:%M:%S', tz = "UTC")
usl$time <- as.POSIXct(usl$time,format='%m.%d.%Y %H:%M:%S', tz = "UTC")
lsl$time <- as.POSIXct(lsl$time,format='%m.%d.%Y %H:%M:%S', tz = "UTC")
oreo$time <- as.POSIXct(oreo$time,format='%m.%d.%Y %H:%M:%S', tz = "UTC")


```

# Calculate tidal speed

The following plot shows the tidal constituents, the combined tidal currents, and a stickplot indicating the speed and direction of the tidal flow. 


```{r calc_tidal_speed, fig.height=12, fig.width=12}
layout(matrix(c(1:16), ncol = 4, byrow = FALSE))


# layout(c(1:12), ncol = 4, byrow = FALSE)
# Bank
bank$k1_speed = sqrt(bank$k1_u^2  + bank$k1_v^2)
bank$m2_speed = sqrt(bank$m2_u^2  + bank$m2_v^2)

bank$tot_u = bank$k1_u + bank$m2_u
bank$tot_v = bank$k1_v + bank$m2_v
bank$tot_speed = sqrt(bank$tot_u^2 + bank$tot_v^2)


# Upper slope
# Upper slope
usl$k1_speed = sqrt(usl$k1_u^2  + usl$k1_v^2)
usl$m2_speed = sqrt(usl$m2_u^2  + usl$m2_v^2)

usl$tot_u = usl$k1_u + usl$m2_u
usl$tot_v = usl$k1_v + usl$m2_v
usl$tot_speed = sqrt(usl$tot_u^2 + usl$tot_v^2)

# lower slope
lsl$k1_speed = sqrt(lsl$k1_u^2  + lsl$k1_v^2)
lsl$m2_speed = sqrt(lsl$m2_u^2  + lsl$m2_v^2)

lsl$tot_u = lsl$k1_u + lsl$m2_u
lsl$tot_v = lsl$k1_v + lsl$m2_v
lsl$tot_speed = sqrt(lsl$tot_u^2 + lsl$tot_v^2)

# Oreo mound 
oreo$k1_speed = sqrt(oreo$k1_u^2  + oreo$k1_v^2)
oreo$m2_speed = sqrt(oreo$m2_u^2  + oreo$m2_v^2)

oreo$tot_u = oreo$k1_u + oreo$m2_u
oreo$tot_v = oreo$k1_v + oreo$m2_v
oreo$tot_speed = sqrt(oreo$tot_u^2 + oreo$tot_v^2)

yyscale <- max(c(bank$tot_speed, usl$tot_speed, lsl$tot_speed, oreo$tot_speed))
# K1
with(bank, plot(time,k1_speed, xlab = "time", ylab = "u (cm/s)", main = "bank K1", ylim = c(-7,7)), col = "black")
with(bank, lines(time,k1_u, xlab = "time", col = "blue"))
with(bank, lines(time,k1_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))

# M2
with(bank, plot(time,m2_speed, xlab = "time", ylab = "u (cm/s)", main = "bank M2", ylim = c(-7,7)), col = "black")
with(bank, lines(time,m2_u, xlab = "time", col = "blue"))
with(bank, lines(time,m2_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))

# Total
with(bank, plot(time,tot_speed, xlab = "time", ylab = "u (cm/s)", main = "bank combined", ylim = c(-15,15)), col = "black")
with(bank, lines(time,tot_u, xlab = "time", col = "blue"))
with(bank, lines(time,tot_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))
plotSticks(x = bank$time, y = 0, u = bank$tot_u, v = bank$tot_v, yscale = yyscale, main = "bank stickplot" )




# K1
with(usl, plot(time,k1_speed, xlab = "time", ylab = "u (cm/s)", main = "upper slope K1", ylim = c(-7,7)), col = "black")
with(usl, lines(time,k1_u, xlab = "time", col = "blue"))
with(usl, lines(time,k1_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))

with(usl, plot(time,m2_speed, xlab = "time", ylab = "u (cm/s)", main = "upper slope M2", ylim = c(-7,7)), col = "black")
with(usl, lines(time,m2_u, xlab = "time", col = "blue"))
with(usl, lines(time,m2_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))

with(usl, plot(time,tot_speed, xlab = "time", ylab = "u (cm/s)", main = "upper slope combined", ylim = c(-15,15)), col = "black")
with(usl, lines(time,tot_u, xlab = "time", col = "blue"))
with(usl, lines(time,tot_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))
plotSticks(x = usl$time, y = 0, u = usl$tot_u, v = usl$tot_v, yscale = yyscale, main = "upper slope stickplot")


# Lowe slope


# K1
with(lsl, plot(time,k1_speed, xlab = "time", ylab = "u (cm/s)", main = "lower slope K1", ylim = c(-7,7)), col = "black")
with(lsl, lines(time,k1_u, xlab = "time", col = "blue"))
with(lsl, lines(time,k1_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))

with(lsl, plot(time,m2_speed, xlab = "time", ylab = "u (cm/s)", main = "lower slope M2", ylim = c(-7,7)), col = "black")
with(lsl, lines(time,m2_u, xlab = "time", col = "blue"))
with(lsl, lines(time,m2_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))

with(lsl, plot(time,tot_speed, xlab = "time", ylab = "u (cm/s)", main = "lower slope combined",  ylim = c(-15,15)), col = "black")
with(lsl, lines(time,tot_u, xlab = "time", col = "blue"))
with(lsl, lines(time,tot_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))
plotSticks(x = lsl$time, y = 0, u = lsl$tot_u, v = lsl$tot_v, yscale = yyscale, main = "lower slope stickplot")



# K1
with(oreo, plot(time,k1_speed, xlab = "time", ylab = "u (cm/s)", main = "oreo K1", ylim = c(-7,7)), col = "black")
with(oreo, lines(time,k1_u, xlab = "time", col = "blue"))
with(oreo, lines(time,k1_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))

with(oreo, plot(time,m2_speed, xlab = "time", ylab = "u (cm/s)", main = "oreo M2", ylim = c(-7,7)), col = "black")
with(oreo, lines(time,m2_u, xlab = "time", col = "blue"))
with(oreo, lines(time,m2_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))

with(oreo, plot(time,tot_speed, xlab = "time", ylab = "u (cm/s)", main = "Oreo combined", ylim =   c(-15,15)), col = "black")
with(oreo, lines(time,tot_u, xlab = "time", col = "blue"))
with(oreo, lines(time,tot_v, xlab = "time", col = "red"))
legend("topright", legend = c("speed", "u", "v"), lty = 1, col = c("black", "blue", "red"))
plotSticks(x = oreo$time, y = 0, u = oreo$tot_u, v = oreo$tot_v, yscale =yyscale, main = "oreo stickplot")

```


# Plot stickplots seperately


```{r stickplots, fig.height = 12, fig.width=4, eval = FALSE}
tiff("figures/fig2_tides.tiff", width = 4, height = 12, unit = "in", res = 300)   
par(mfrow = c(4,1))
  plotSticks(x = bank$time, y = 0, u = bank$tot_u, v = bank$tot_v, yscale = yyscale)
  plotSticks(x = usl$time, y = 0, u = usl$tot_u, v = usl$tot_v, yscale = yyscale)
  plotSticks(x = lsl$time, y = 0, u = lsl$tot_u, v = lsl$tot_v, yscale = yyscale)
  plotSticks(x = oreo$time, y = 0, u = oreo$tot_u, v = oreo$tot_v, yscale =yyscale)
dev.off( )


```


# Retrieve direction of tidal current

One thing to notice: bank, upper slope, lower slope, all follow quite well a clockwise rotation of the tidal direction, but Oreo mound has a swing in between. 

```{r get direction}
# Bank
bank <- bank %>%
  as_tibble() %>% 
  mutate(flowdir = (90 - (180 / pi) *  atan2(tot_v,tot_u))%%360)
# Upper slope
usl <- usl %>%
  as_tibble() %>% 
  mutate(flowdir = (90 - (180 / pi) *  atan2(tot_v,tot_u))%%360)
# Lower slope
lsl <- lsl %>%
  as_tibble() %>% 
  mutate(flowdir = (90 - (180 / pi) *  atan2(tot_v,tot_u))%%360)
# Oreo
oreo <- oreo %>%
  as_tibble() %>% 
  mutate(flowdir = (90 - (180 / pi) *  atan2(tot_v,tot_u))%%360)

ggplot(data = bank, aes(x = time, y = flowdir))+
  geom_point()
ggplot(data = usl,aes(x = time, y = flowdir))+
  geom_point()
ggplot(data = lsl,aes(x = time, y = flowdir))+
  geom_point()
ggplot(data = oreo,aes(x = time, y = flowdir))+
  geom_point()

```

# Extract times of slack and high tidal speeds

For the bank, there is a prolonged period of time with low tidal current speeds. 


```{r extract slack/high velo times}
# Bank

sl1 <- bank[nth(bank$tot_speed,5, index.return = TRUE ),]$time # slack end
sl2 <- bank[nth(bank$tot_speed, 1, index.return = TRUE),]$time # slack begind

hi1 <- bank[nth(bank$tot_speed,1, index.return = TRUE, descending = TRUE),]$time # slack end
hi2 <- bank[nth(bank[1:20,]$tot_speed,1, index.return = TRUE, descending = TRUE),]$time
hi3 <- bank[bank$tot_speed == nth(bank[90:nrow(bank),]$tot_speed,k =  1,descending = TRUE),]$time
ggplot(data = bank, aes(x = time, y = tot_speed))+
  geom_point()+
  geom_vline(xintercept = sl1, linetype = "dotted")+
  geom_vline(xintercept = sl2, linetype = "dotted")+
  geom_vline(xintercept = hi1, linetype = "dotted", col = "red")+
  geom_vline(xintercept = hi2, linetype = "dotted", col = "red")+
  geom_vline(xintercept = hi3, linetype = "dotted", col = "red")

bank_times <- as.data.frame(matrix(data = NA, nrow = 5, ncol = 3 ))
names(bank_times)   <- c("area", "descript", "time")
bank_times$area     <- "bank"
bank_times$descript <- c("S_W", "S_W", "H_N", "H_S", "H_N")
bank_times$time     <- c(sl1, sl2, hi1, hi2, hi3)


# upper slope
sl1 <- usl[1,]$time # slack first
sl2 <- usl[nth(usl$tot_speed,2, index.return = TRUE),]$time # slack 2
sl3 <- usl[usl$tot_speed == nth(usl[50:90,]$tot_speed,k =  1),]$time
sl4 <- usl[97,]$time # slack end

hi1 <- usl[usl$tot_speed == nth(usl[1:30,]$tot_speed,k =  1,descending = TRUE),]$time
hi2 <- usl[nth(usl$tot_speed,1, index.return = TRUE, descending = TRUE),]$time # slack end
hi3 <- usl[usl$tot_speed == nth(usl[1:70,]$tot_speed,k =  1, descending = TRUE),]$time


ggplot(data = usl, aes(x = time, y = tot_speed))+
  geom_point()+
  geom_vline(xintercept = c(sl1,sl2, sl3, sl4), linetype = "dotted")+
  geom_vline(xintercept = c(hi1, hi2, hi3), linetype = "dotted", col = "red")
  


usl_times <- as.data.frame(matrix(data = NA, nrow =7, ncol = 3 ))
names(usl_times)   <- c("area", "descript", "time")
usl_times$area     <- "upper_slope"
usl_times$descript <- c("S_NW", "S_NW", "S_NE", "S_SW","H_NW", "H_S","H_N")
usl_times$time     <- c(sl1, sl2,sl3, sl4, hi1, hi2, hi3)

# Lower slope
sl1 <- lsl[nth(lsl$tot_speed, 2, index.return = TRUE),]$time # slack west
sl2 <- lsl[lsl$tot_speed == nth(lsl[30:50,]$tot_speed,k =  1),]$time
sl3 <- lsl[nth(lsl$tot_speed, 1, index.return = TRUE),]$time # slack begind

# high velocity
hi1 <- lsl[lsl$tot_speed == nth(lsl[1:40,]$tot_speed,k =  1,descending = TRUE),]$time
hi2 <- lsl[lsl$tot_speed == nth(lsl$tot_speed,k =  1,descending = TRUE),]$time
hi3 <- lsl[lsl$tot_speed == nth(lsl[80:nrow(lsl),]$tot_speed,k =  1,descending = TRUE),]$time

ggplot(data = lsl, aes(x = time, y = tot_speed))+
  geom_point()+
  geom_vline(xintercept = c(sl1, sl2, sl3), linetype = "dotted")+
  geom_vline(xintercept = c(hi1, hi2, hi3), linetype = "dotted", col = "red")

lsl_times <- as.data.frame(matrix(data = NA, nrow =6, ncol = 3 ))
names(lsl_times)   <- c("area", "descript", "time")
lsl_times$area     <- "lower_slope"
lsl_times$descript <- c("S_W", "S_E", "S_W", "H_N","H_S", "H_N")
lsl_times$time     <- c(sl1, sl2,sl3, hi1, hi2, hi3)

# Oreo mound
## Slack
sl1 <- oreo[oreo$tot_speed == nth(oreo[1:20,]$tot_speed,k =  1),]$time
sl2 <- oreo[nth(oreo$tot_speed, 1, index.return = TRUE),]$time # slack west
sl3 <- oreo[oreo$tot_speed == nth(oreo[80:nrow(oreo),]$tot_speed,k =  1),]$time 

# High velo
hi1 <- oreo[oreo$tot_speed == nth(oreo[1:20,]$tot_speed,k =  1,descending = TRUE),]$time
hi2 <- oreo[oreo$tot_speed == nth(oreo[20:40,]$tot_speed,k =  1,descending = TRUE),]$time
hi3 <- oreo[oreo$tot_speed == nth(oreo[60:80,]$tot_speed,k =  1,descending = TRUE),]$time
hi4 <- oreo[oreo$tot_speed == nth(oreo[80:nrow(oreo),]$tot_speed,k =  1,descending = TRUE),]$time


ggplot(data = oreo, aes(x = time, y = tot_speed))+
  geom_point()+
  geom_vline(xintercept = c(sl1, sl2, sl3), linetype = "dotted")+
  geom_vline(xintercept = c(hi1, hi2, hi3, hi4), linetype = "dotted", col = "red")

oreo_times <- as.data.frame(matrix(data = NA, nrow =7, ncol = 3 ))
names(oreo_times)   <- c("area", "descript", "time")
oreo_times$area     <- "oreo"
oreo_times$descript <- c("S_E", "S_SW", "S_W", "H_N","H_S", "H_SW", "H_N2")
oreo_times$time     <- c(sl1, sl2,sl3, hi1, hi2, hi3, hi4)


tidal_times <- rbind(bank_times, usl_times, lsl_times, oreo_times)


```


```{r write dataset, eval = FALSE}

writexl::write_xlsx(tidal_times, "Data/data_cleaned/tidal_times.xlsx")
tides <- list(bank, usl, lsl, oreo)
save(file = "Data/CTD/tides.rda", tides)


```


