Packages & setup

# 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"))
d_vars    <- calcofi4r::get_variables()
d_cruises <- calcofi4r::get_cruises()
datatable(d_cruises)
# 1. choose cruise
(cruiseid <- d_cruises$cruiseid[1])
## [1] "2020-01-05-C-33RL"
# 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)
## 
##   60 63.3 66.7   70 73.3 76.7   80 81.7 81.8 83.3 85.4 86.7 86.8 88.5   90 91.7 
##    5    6    6    6    6    8    8    1    1   11    1   12    1    1   14    1 
## 93.3 93.4 
##   15    1
# 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)
## 
##   60 63.3 66.7   70 73.3 76.7   80 83.3 86.7   90 93.3 
##    5    6    6    6    6    8    8   11   12   14   15
# 2. choose station line
(sta_line <- casts$sta_line[1])
## [1] 60
# 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

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
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

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()

# rever to positive depth
d <- d %>% 
    mutate(depthm = depthm * -1)
# View(d)
skim(d)
Data summary
Name d
Number of rows 99
Number of columns 11
_______________________
Column type frequency:
character 1
Date 1
numeric 8
________________________
Group variables sta_line

Variable type: character

skim_variable sta_line n_missing complete_rate min max empty n_unique whitespace
sta_id 60 0 1 11 11 0 5 0

Variable type: Date

skim_variable sta_line n_missing complete_rate min max median n_unique
date 60 0 1 2020-01-25 2020-01-26 2020-01-25 2

Variable type: numeric

skim_variable sta_line n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
cast_count 60 0 1 35275.28 1.31 35273.00 35274.00 35275.00 35276.00 35277.00 ▃▇▇▇▇
sta_offshore 60 0 1 73.13 12.58 53.00 60.00 70.00 80.00 90.00 ▇▁▆▆▆
longitude 60 0 1 -124.55 0.91 -125.77 -125.05 -124.33 -123.60 -123.10 ▆▆▆▁▇
latitude 60 0 1 37.18 0.42 36.61 36.95 37.28 37.61 37.85 ▆▆▆▁▇
day 60 0 1 20.23 0.42 20.00 20.00 20.00 20.00 21.00 ▇▁▁▁▂
sta_line_n 60 0 1 5.00 0.00 5.00 5.00 5.00 5.00 5.00 ▁▁▇▁▁
depthm 60 0 1 157.24 154.88 0.00 30.00 100.00 250.00 516.00 ▇▃▂▁▁
t_degc 60 0 1 9.73 2.55 4.76 7.94 9.67 12.29 13.20 ▃▃▅▂▇ 
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")
# 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: beautiful colormaps for oceanography
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

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)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -4590   -4236   -3968   -3737   -3567     -70
##                   36.576923077 36.6502262444118 36.7235294118235
## -125.8                   -4482            -4497            -4502
## -125.728947368421        -4506            -4495            -4458
## -125.657894736842        -4502            -4502            -4356
## -125.586842105263        -4485            -4498            -4158
## -125.515789473684        -4491            -4461            -4113
##                   36.7968325792353 36.8701357466471
## -125.8                       -4462            -4452
## -125.728947368421            -4309            -4340
## -125.657894736842            -3586            -3430
## -125.586842105263            -3101            -2992
## -125.515789473684            -3509            -3692
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
## # A tibble: 2 × 3
##   sta_offshore   lon   lat
##          <dbl> <dbl> <dbl>
## 1           53 -123.  37.8
## 2           90 -126.  36.6
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
##          lon      lat    dist.km depth
## 1  -125.8000 36.57692   0.000000 -4482
## 2  -125.7289 36.60972   7.316763 -4506
## 3  -125.6579 36.64251  14.631186 -4502
## 4  -125.5868 36.67530  21.943267 -4498
## 5  -125.5158 36.70810  29.253002 -4113
## 6  -125.4447 36.74089  36.560388 -3975
## 7  -125.3737 36.77368  43.865422 -3661
## 8  -125.3026 36.80648  51.168102 -4068
## 9  -125.2316 36.83927  58.468423 -4008
## 10 -125.1605 36.87206  65.766383 -4049
## 11 -125.0895 36.90486  73.061978 -4333
## 12 -125.0184 36.93765  80.355206 -4339
## 13 -124.9474 36.97045  87.646063 -4280
## 14 -124.8763 37.00324  94.934547 -4245
## 15 -124.8053 37.03603 102.220653 -4229
## 16 -124.7342 37.06883 109.504379 -4178
## 17 -124.6632 37.10162 116.785722 -4107
## 18 -124.5921 37.13441 124.064678 -3954
## 19 -124.5211 37.16721 131.341245 -3857
## 20 -124.4500 37.20000 138.615418 -3938
## 21 -124.3789 37.23279 145.887196 -3959
## 22 -124.3079 37.26559 153.156575 -3980
## 23 -124.2368 37.29838 160.423551 -3956
## 24 -124.1658 37.33117 167.688122 -3922
## 25 -124.0947 37.36397 174.950284 -3820
## 26 -124.0237 37.39676 182.210034 -3764
## 27 -123.9526 37.42955 189.467370 -3735
## 28 -123.8816 37.46235 196.722287 -3636
## 29 -123.8105 37.49514 203.974782 -3600
## 30 -123.7395 37.52794 211.224853 -3509
## 31 -123.6684 37.56073 218.472496 -3420
## 32 -123.5974 37.59352 225.717708 -3364
## 33 -123.5263 37.62632 232.960486 -3287
## 34 -123.4553 37.65911 240.200826 -3006
## 35 -123.3842 37.69190 247.438725 -2655
## 36 -123.3132 37.72470 254.674181 -1894
## 37 -123.2421 37.75749 261.907190 -1113
## 38 -123.1711 37.79028 269.137748  -110
## 39 -123.1000 37.82308 276.365853   -70
plotProfile(trsect)

calcofi4r::plot_transect()

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)
}