Packages & setup

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

source(here("../apps/libs/db.R"))
calcofi4r::stations %>% 
  st_drop_geometry() %>% 
  skimr::skim()
Data summary
Name Piped data
Number of rows 2634
Number of columns 9
_______________________
Column type frequency:
character 1
logical 4
numeric 4
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
sta_id 0 1 11 11 0 2634 0

Variable type: logical

skim_variable n_missing complete_rate mean count
is_offshore 0 1 0.42 FAL: 1534, TRU: 1100
is_cce 0 1 0.04 FAL: 2521, TRU: 113
is_ccelter 0 1 0.03 FAL: 2568, TRU: 66
is_sccoos 0 1 0.00 FAL: 2625, TRU: 9

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
sta_id_line 0 1 92.01 31.45 1.00 73.30 93.00 113.30 176.70 ▁▃▇▅▁
sta_id_station 0 1 72.87 67.35 0.00 39.00 56.00 80.00 531.00 ▇▁▁▁▁
lon 0 1 -120.74 6.71 -153.16 -123.96 -119.87 -116.33 -107.13 ▁▁▂▇▃
lat 0 1 31.64 5.28 17.41 28.52 31.60 34.49 47.94 ▁▅▇▂▁ 
nrow(calcofi4r::stations) # 2,634
## [1] 2634
mapview(calcofi4r::stations, zcol="sta_id_line")

Weber & Moore (2013)

Source: Appendix of Weber & Moore (2013) Corrected conversion algorithms for the CalCOFI station grid and their implementation in several computer languages. California Cooperative Oceanic Fisheries Investigations Reports

deg2rad <- function(deg) deg * pi / 180
rad2deg <- function(rad) rad * 180 / pi
inverse_mercator <- function(mercatorlat, iterations = 3){ 
  
  approxlat <- mercatorlat
  
  iterlatitude <- function(mercatorlat, approxlat){
    approxlat <- 2 * (atan(exp(deg2rad(mercatorlat) + 0.00676866 * sin(deg2rad(approxlat)))) * 180 / pi - 45)
    approxlat
  }
  
  for (i in 1:iterations)
    approxlat <  iterlatitude(mercatorlat, approxlat)
  
  approxlat
}

to_mercator <- function(latitude) { 
  y <- rad2deg(log(tan(deg2rad(45 + latitude / 2))) - 0.00676866 * sin(deg2rad(latitude)))
  y
}

station_to_latlon <- function(x, roundlines = true) {
  
  if (length(x) == 2 & class(x) != 'matrix'){
    x <- matrix(x, 1, 2)
  }
  line <- x[, 1]
  station <- x[, 2]
  
  reflatitude <- 34.15 - 0.2 * (line - 80) * cos(deg2rad(30))
  latitude <- reflatitude - (station - 60) * sin(deg2rad(30)) / 15
  l1 <- (to_mercator(latitude) - to_mercator(34.15)) * tan(deg2rad(30))
  l2 <- (to_mercator(reflatitude) - to_mercator(latitude)) / (cos(deg2rad(30)) * sin(deg2rad(30)))
  longitude <- -1 * (l1 + l2 + 121.15)
  cbind(lon = longitude, lat = latitude)
}
# https://proj.org/operations/projections/calcofi.html # -121.15 34.15   80.00   60.00
station_to_latlon(c(80.0, 60.0))
##          lon   lat
## [1,] -121.15 34.15
#          lon   lat
# [1,] -121.15 34.15

lonlat_to_station <- function(x){
  # x = c(-121.15, 34.15)
  if (length(x) == 2 & class(x) != 'matrix'){
    x <- matrix(x, 1, 2)
  }
  longitude <- x[, 1]
  latitude  <- x[, 2]
  # longitude <- -121.15 
  # latitude  <- 34.15
  
  # assume we're in the western hemispere
  longitude[longitude > 180] <- -1 * (longitude[longitude > 180] - 360)
  longitude[longitude < 0]   <- longitude[longitude < 0] * -1
  
  l1 <- (to_mercator(latitude) - to_mercator(34.15)) * tan(deg2rad(30))
  l2 <- longitude - l1 - 121.15
  mercreflatitude <- l2 * cos(deg2rad(30)) * sin(deg2rad(30)) + to_mercator(latitude)
  reflatitude     <- inverse_mercator(mercreflatitude)
  
  line    <- 80 - (reflatitude - 34.15) * 5 / cos(deg2rad(30))
  station <- 60 + (reflatitude - latitude) * 15 / sin(deg2rad(30))
  
  cbind(line = line, station = station)
}
# https://proj.org/operations/projections/calcofi.html # 80.0 60.0   -121.15 34.15
lonlat_to_station(c(-121.15, 34.15))
##          line station
## [1,] 68.42567 120.142
# 68.42567 120.142 # Doh! Wrong answer

proj-bin

For Debian Linux, install proj binary:

sudo apt-get install proj-bin
which proj
## Reading package lists...
## Building dependency tree...
## Reading state information...
## proj-bin is already the newest version (6.3.1-1).
## 0 upgraded, 0 newly installed, 0 to remove and 13 not upgraded.
## /usr/bin/proj

Define functions:

lonlat_to_stationid <- function(lon, lat){
  # convert station ID to lon, lat using the proj library
  
  system(glue("echo {lon} {lat} | proj +proj=calcofi +epsg=4326 -f '%05.1f'"), intern=T) %>%
    stringr::str_replace("\t", " ")
}
lonlat_to_stationid(-121.15, 34.15) # "080.0 060.0"
## [1] "080.0 060.0"
stationid_to_lonlat <- function(stationid){
  # using 5th decimal place, a la CCE_Stations.txt
  
  system(glue("echo {stationid} | proj +proj=calcofi +epsg=4326  -I -d 5"), intern=T) %>%
    stringr::str_replace("\t", " ")
}
stationid_to_lonlat("080.0 060.0") # "-121.15000 34.15000"
## [1] "-121.15000 34.15000"
# A: offshore/1, line/1
sta_a.1.1 <- calcofi4r::stations %>% 
  st_drop_geometry() %>% 
  mutate(
    x = round(sta_id_station),
    y = round(sta_id_line)) %>% 
  group_by(x, y) %>% 
  summarize(
    n = n(),
    .groups = "drop") %>% 
  mutate(
    z = glue("{y} {x}")) %>%
  mutate(
    lonlat = map_chr(z, stationid_to_lonlat)) %>% 
  separate(lonlat, c("lon", "lat"), sep = " ", convert=T) %>%
  st_as_sf(
    coords = c("lon", "lat"), crs=4326, remove = F) %>% 
  rename(
    n_stations = n,
    sta_id = z)
nrow(sta_a.1.1) # 2,270
## [1] 2270
mapview(sta_a.1.1, zcol="y")
# B: offshore/5, line/2
sta_b.5.2 <- calcofi4r::stations %>% 
  st_drop_geometry() %>% 
  mutate(
    x = round(sta_id_station/5) * 5,
    y = round(sta_id_line/2) * 2) %>% 
  group_by(x, y) %>% 
  summarize(
    n = n(),
    .groups = "drop") %>% 
  mutate(
    z = glue("{y} {x}")) %>%
  mutate(
    lonlat = map_chr(z, stationid_to_lonlat)) %>% 
  separate(lonlat, c("lon", "lat"), sep = " ", convert=T) %>%
  st_as_sf(
    coords = c("lon", "lat"), crs=4326, remove = F) %>% 
  rename(
    n_stations = n,
    sta_id = z)
nrow(sta_b.5.2) # 1,128
## [1] 1128
mapview(sta_b.5.2, zcol="y")
# C: offshore/10, line/3
sta_c.10.3 <- calcofi4r::stations %>% 
  st_drop_geometry() %>% 
  mutate(
    x = round(sta_id_station/10) * 10,
    y = round(sta_id_line/3) * 3) %>% 
  group_by(x, y) %>% 
  summarize(
    n = n(),
    .groups = "drop") %>% 
  mutate(
    z = glue("{y} {x}")) %>%
  mutate(
    lonlat = map_chr(z, stationid_to_lonlat)) %>% 
  separate(lonlat, c("lon", "lat"), sep = " ", convert=T) %>%
  st_as_sf(
    coords = c("lon", "lat"), crs=4326, remove = F) %>% 
  rename(
    # sta_offshore = x,
    # sta_alongshore = y,
    n_stations = n,
    sta_id = z)
nrow(sta_c.10.3) # 702
## [1] 702
mapview(sta_c.10.3, zcol="y")
# D: offshore/10, line/4
sta_d.10.4 <- calcofi4r::stations %>% 
  st_drop_geometry() %>% 
  mutate(
    x = round(sta_id_station/10) * 10,
    y = round(sta_id_line/4) * 4) %>% 
  group_by(x, y) %>% 
  summarize(
    n = n(),
    .groups = "drop") %>% 
  mutate(
    z = glue("{y} {x}")) %>%
  mutate(
    lonlat = map_chr(z, stationid_to_lonlat)) %>% 
  separate(lonlat, c("lon", "lat"), sep = " ", convert=T) %>%
  st_as_sf(
    coords = c("lon", "lat"), crs=4326, remove = F) %>% 
  rename(
    # sta_offshore = x,
    # sta_alongshore = y,
    n_stations = n,
    sta_id = z)
nrow(sta_d.10.4) # 611
## [1] 611
mapview(sta_d.10.4, zcol="y")
hull <- st_convex_hull(st_union(calcofi4r::stations)) %>% 
  st_transform(leaflet:::epsg3857) %>% 
  st_buffer(1000)

sta_d.10.4_v <- sta_d.10.4 %>% 
  st_transform(leaflet:::epsg3857) %>% 
  st_union() %>% 
  st_voronoi(hull)

sta_d.10.4_vp <- sta_d.10.4_v %>% 
  st_collection_extract(
  type = "POLYGON")
sta_d.10.4_vpi <- st_intersection(sta_d.10.4_vp, hull)
sta_d.10.4_vpig <- sta_d.10.4_vpi %>% 
  st_transform(4326)

mapview(sta_d.10.4_vpi) +
  mapview(sta_d.10.4, zcol="y")
librarian::shelf(
  rnaturalearth, rnaturalearthdata)
# devtools::install_github("ropensci/rnaturalearthhires")

land_l <- ne_countries(scale = "large", returnclass = "sf") %>% 
  st_combine() %>% 
  st_make_valid()
# land_m <- ne_countries(scale = "medium", returnclass = "sf")
# land_s <- ne_countries(scale = "small", returnclass = "sf")

# sta_d.10.4_vpigl_0 <- sta_d.10.4_vpigl
sta_d.10.4_vpigl <- st_difference(sta_d.10.4_vpig, land_l)
mapview(sta_d.10.4_vpigl, alpha.regions = 0.5) +
  mapview(sta_d.10.4, zcol="y", cex=2)