explore_depth.Rmd

---
title: "Depth Profile Plot Explore"
author: "Ben Best"
date: "`r Sys.Date()`"
output: 
  html_document:
    toc: true
    toc_float: true
    code_folding: show
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = F, message = F)
```

## Packages & setup

```{r}
# packages
if (!require("librarian")){
  install.packages("librarian")
  library(librarian)
}
librarian::shelf(
  calcofi/calcofi4r, 
  # marmap,
  cmocean, dplyr, DT, glue, here, interp,
  mapview, plotly, sf, skimr, tidyr)
mapviewOptions(fgb = F)

source(here("../apps/libs/db.R"))
```

```{r}
d_vars    <- calcofi4r::get_variables()
d_cruises <- calcofi4r::get_cruises()
datatable(d_cruises)

# 1. choose cruise
(cruiseid <- d_cruises$cruiseid[1])

# get casts, filtering by cruise
casts <- tbl(con, "ctd_casts") %>% 
  filter(cruiseid == !!cruiseid) %>% 
  select(cast_count, sta_id, date, longitude, latitude) %>% 
  collect() %>% 
  separate(
    sta_id, into = c("sta_line", "sta_offshore"), 
    sep = " ", convert = T, remove = F) %>% 
  mutate(
    day = difftime(date, min(date), units="days") %>% 
      as.integer()) %>% 
  st_as_sf(
    coords = c("longitude", "latitude"),
    crs = 4326, remove = F)

mapview(casts, zcol="sta_line")
# mapview(casts, zcol="sta_offshore")
# mapview(casts, zcol="day")

table(casts$sta_line)

# remove station lines with only one cast
casts <- casts %>% 
  group_by(sta_line) %>% 
  mutate(sta_line_n = n()) %>% 
  filter(sta_line_n > 1)

table(casts$sta_line)

# 2. choose station line
(sta_line <- casts$sta_line[1])

# filter by station line
casts <- casts %>% 
  filter(
    sta_line == !!sta_line)

bottles <- dbGetQuery(
  con,
  glue(
    "SELECT cast_count, depthm, t_degc 
    FROM ctd_bottles
    WHERE cast_count IN ({paste(casts$cast_count, collapse=',')})"))
d <- casts %>%
  inner_join(
    bottles,
    by="cast_count") %>% 
  st_drop_geometry() %>% 
  # select(sta_offshore, sta_line, depthm, t_degc) %>% 
  arrange(sta_offshore, sta_line, depthm)
datatable(d)
```

## plotly

```{r}
fig <- plot_ly(
  type='isosurface',
  x = d$sta_offshore,
  y = d$sta_line,
  z = d$depthm,
  value = d$t_degc,
  isomin = min(d$t_degc),
  isomax = max(d$t_degc))
fig
```

```{r}
library(plotly)

fig <- plot_ly(
  type='isosurface',
  x = c(0,0,0,0,1,1,1,1),
  y = c(1,0,1,0,1,0,1,0),
  z = c(1,1,0,0,1,1,0,0),
  value = c(1,2,3,4,5,6,7,8),
  isomin=2,
  isomax=6
  )

fig
```

## metR

* [Visualization tools • metR](https://eliocamp.github.io/metR/articles/Visualization-tools.html)

```{r}
shelf(metR)

d <- d %>% 
    mutate(depthm = depthm * -1)
  
ggplot(d, aes(sta_offshore, depthm, z = t_degc)) +
  geom_contour_fill(na.fill = TRUE)

ggplot(d, aes(sta_offshore, depthm)) +
  geom_contour_fill(aes(z = t_degc), kriging = T)

ggplot(d, aes(sta_offshore, depthm)) +
  geom_contour_fill(aes(z = t_degc), kriging = T) + 
  geom_point(size = 0.2)
```

## `ggplot2::geom_contour_filled()`

```{r}
# rever to positive depth
d <- d %>% 
    mutate(depthm = depthm * -1)
# View(d)
skim(d)

d %>% 
  ggplot(aes(x = sta_offshore, y = depthm)) +
  geom_contour_filled(aes(z = t_degc), bins = 10)+
  geom_point(size = 0.2) +
  scale_y_reverse() + 
  coord_cartesian(expand = F) +
  labs(x = "Offshore [id]", y = "Depth [m]")+
  theme_bw() %+%
  theme(
    panel.background  = element_rect(fill = "grey90"),
    panel.grid.major  = element_line(linetype = 3, colour = "grey60"),
    axis.text         = element_text(colour = 1, size = 10),
    axis.title        = element_text(colour = 1, size = 12),
    legend.background = element_blank(),
    legend.key        = element_blank(),
    legend.position   = "right") # +
  # metR::scale_x_longitude(ticks = 0.005,position = "bottom")
```

* [2D contours of a 3D surface — geom_contour • ggplot2](https://ggplot2.tidyverse.org/reference/geom_contour.html)


```{r}
# Irregular data
# Use a dataset from the interp package
data(franke, package = "interp")

origdata <- as.data.frame(interp::franke.data(1, 1, franke))

grid <- with(origdata, interp::interp(x, y, z))

griddf <- subset(
  data.frame(
    x = rep(grid$x, nrow(grid$z)),
    y = rep(grid$y, each = ncol(grid$z)),
    z = as.numeric(grid$z)),
  !is.na(z))

ggplot(griddf, aes(x, y, z = z)) +
  geom_contour_filled() +
  geom_point(data = origdata)
```

- [cmocean](https://cran.r-project.org/web/packages/cmocean/vignettes/cmocean.html): beautiful colormaps for oceanography

```{r}
g <- with(
  d, 
  interp::interp(
    sta_offshore, depthm, t_degc))

gd <- subset(
  data.frame(
    x = rep(g$x, nrow(g$z)),
    y = rep(g$y, each = ncol(g$z)),
    z = as.numeric(g$z)),
  !is.na(z))

nbins = 8
cmocean_pal = "thermal"

p <- ggplot(gd, aes(x, y, z = z)) +
  geom_contour_filled(bins=nbins) +
  coord_cartesian(expand = F) +
  scale_x_reverse() + 
  scale_y_reverse() + 
  labs(
    x     = "Offshore [id]", 
    y     = "Depth [m]",
    fill = "Temperature [C]") +
  scale_fill_manual(values = cmocean(cmocean_pal)(nbins)) +
  theme_bw()
p

# add points
p +
  geom_point(
    data = d %>%
      rename(
        x = sta_offshore,
        y = depthm,
        z = t_degc),
    size = 0.2)
```

## `marmap`: add bathymetry

```{r}
library(marmap)

b <- with(
  d,
  getNOAA.bathy(
    lon1 = min(longitude),
    lon2 = max(longitude),
    lat1 = min(latitude),
    lat2 = max(latitude),
    resolution=4))

# create color palettes
blues <- c("lightsteelblue4","lightsteelblue3","lightsteelblue2","lightsteelblue1")
greys <- c(grey(0.6),grey(0.93),grey(0.99))

plot(b, image = TRUE, land = TRUE, lwd = 0.03, bpal = list(c(0, max(b), greys), c(min(b), 0, blues)))
# Add coastline
plot(b, n = 1, lwd = 0.4, add = TRUE)

summary(b)


dr <- d %>% 
  filter(
    sta_offshore %in% c(max(sta_offshore), min(sta_offshore))) %>% 
  group_by(sta_offshore) %>% 
  summarize(
    lon = first(longitude),
    lat = first(latitude))

# Latitudinal  range:   36.58 to   37.82 ( 36.58 N to  37.82 N)
# Longitudinal range: -125.8  to -123.1  (125.8  W to 123.1  W)

dr

trsect <- get.transect(
    b,
    x1 = min(as.numeric(dimnames(b)[[1]])),
    x2 = max(as.numeric(dimnames(b)[[1]])),
    y1 = min(as.numeric(dimnames(b)[[2]])),
    y2 = max(as.numeric(dimnames(b)[[2]])),
    distance = T)
trsect

plotProfile(trsect)
```


## `calcofi4r::plot_transect()`

```{r, eval=F}

plot_transect <- function(x, y, z){
  g <- with(
  d, 
  interp::interp(
    sta_offshore, depthm, t_degc))

gd <- subset(
  data.frame(
    x = rep(g$x, nrow(g$z)),
    y = rep(g$y, each = ncol(g$z)),
    z = as.numeric(g$z)),
  !is.na(z))

nbins = 8
cmocean_pal = "thermal"

p <- ggplot(gd, aes(x, y, z = z)) +
  geom_contour_filled(bins=nbins) +
  coord_cartesian(expand = F) +
  scale_x_reverse() + 
  scale_y_reverse() + 
  labs(
    x     = "Offshore [id]", 
    y     = "Depth [m]",
    fill = "Temperature [C]") +
  scale_fill_manual(values = cmocean(cmocean_pal)(nbins)) +
  theme_bw()
p

# add points
p +
  geom_point(
    data = d %>%
      rename(
        x = sta_offshore,
        y = depthm,
        z = t_degc),
    size = 0.2)
}

```