librarian:: shelf (
calcofi/ calcofi4r,
dplyr, glue, here, jsonlite, leaflet, lubridate, mapview, purrr, sf, terra,
tibble, tidyr)
cm = "/opt/homebrew/Caskroom/mambaforge/base/envs/copernicusmarine/bin/copernicusmarine"
user = "bbest1"
# pass = readLines("~/My Drive/private/data.marine.copernicus.eu_bbest1-password.txt")
pass = "temporary"
verbose = F
# CalCOFI bounding box
b <- st_union (cc_grid_zones) |> st_bbox ()
# b
# xmin ymin xmax ymax
# -135.23008 18.42757 -105.77692 49.23891
m_json <- here (glue ("data/copernicusmarine/{dataset_id}.json" ))
cmd <- glue ("{cm} describe --include-datasets -c {dataset_id}" )
if (! file.exists (m_json))
system (cmd, intern = TRUE ) |> writeLines (m_json)
m_o <- read_json (m_json)
# listviewer::jsonedit(m_o)
get_coord = function (o, v){
if (! v %in% names (o))
return (NA )
if (length (o[[v]]) == 0 )
return (NA )
o[[v]] |>
as.character ()
}
coord_date <- function (df, col){
df |>
filter (coord_id == "time" ) |>
pull (!! col) |>
(\(x) as.double (x) / 1000 )() |>
as_datetime (origin = "1970-01-01" , tz = "UTC" ) |>
as_date ()
}
m_vars <- m_o |>
pluck ("products" , 1 , "datasets" , 1 , "versions" , 1 , "parts" , 1 , "services" ) |>
keep (\(x) x$ service_type$ service_name == "arco-geo-series" ) |>
pluck (1 , "variables" )
d_vars_coords <- m_vars %>%
tibble (
short_name = map_chr (., "short_name" ),
standard_name = map_chr (., "standard_name" ),
units = map_chr (., "units" ),
coord = map (., "coordinates" ) ) |>
select (- ` . ` ) |>
unnest (coord) |>
mutate (
coord_id = map_chr (coord, "coordinates_id" ),
coord_units = map_chr (coord, "units" ),
coord_min_val = map_chr (coord, get_coord, "minimum_value" ),
coord_max_val = map_chr (coord, get_coord, "maximum_value" ),
coord_step = map_chr (coord, get_coord, "step" ),
coord_chunk_length = map_chr (coord, get_coord, "chunking_length" ),
coord_chunk_type = map_chr (coord, get_coord, "chunk_type" ),
coord_chunk_ref = map_chr (coord, get_coord, "chunk_reference_coordinate" ),
coord_chunk_geo = map_chr (coord, get_coord, "chunk_geometric_factor" ),
coord_values = map (coord, "values" ),
coord_values_n = map_int (coord_values, length))
#
# d_vars_coords |> select(short_name, coord_id) |> table()
# coord_id
# short_name depth latitude longitude time
# bottomT 1 1 1 1
# mlotst 1 1 1 1
# siconc 1 1 1 1
# sithick 1 1 1 1
# so 1 1 1 1
# thetao 1 1 1 1
# uo 1 1 1 1
# usi 1 1 1 1
# vo 1 1 1 1
# vsi 1 1 1 1
# zos 1 1 1 1
#
# d_vars_coords |> select(coord_id, coord_values_n) |> table()
# coord_values_n
# coord_id 0 50
# depth 0 11
# latitude 11 0
# longitude 11 0
# time 11 0
d_vars <- d_vars_coords |>
nest (.by = c (short_name, standard_name, units)) |>
mutate (
date_min = map_dbl (data, coord_date, "coord_min_val" ) |>
as.Date (),
date_max = map_dbl (data, coord_date, "coord_max_val" ) |>
as.Date ())
depths <- d_vars |>
filter (
short_name == "thetao" ) |>
pull (data) |>
(\(x) x[[1 ]])() |>
filter (
coord_id == "depth" ) |>
pull (coord_values) |>
unlist ()
# depths <- obj$products[[1]]$datasets[[1]]$versions[[1]]$parts[[1]]$services[[2]]$variables[[1]]$coordinates[[4]][["values"]] |>
# unlist()
depth <- depths[length (depths)] * - 1
# output filename
nc <- here (glue ("data/copernicusmarine/{dataset_id}__{variable}__{date}.nc" ))
# https://help.marine.copernicus.eu/en/articles/7970514-copernicus-marine-toolbox-installation
# mamba update --name copernicusmarine copernicusmarine --yes
# copernicusmarine, version 1.2.2
cmd <- glue ("
{cm} subset -i {dataset_id} \\
-x {b$xmin} -X {b$xmax} \\
-y {b$ymin} -Y {b$ymax} \\
-t {date} -T {date + months(1) - days(1)} \\
-z {depth - 0.01} -Z {depth + - 0.01} \\
-v {variable} \\
-o {dirname(nc)} -f {basename(nc)} \\
--username {user} --password '{pass}' \\
--force-download" )
if (verbose)
print (cmd)
if (! file.exists (nc))
system (cmd, intern = TRUE )
r <- rast (nc)
# dimensions : 369, 353, 1 (nrow, ncol, nlyr)
# resolution : 0.08333334, 0.08333334 (x, y)
# extent : -135.2083, -105.7917, 18.45833, 49.20833 (xmin, xmax, ymin, ymax)
r_3857 <- leaflet:: projectRasterForLeaflet (r, method = "bilinear" )
# dimensions : 403, 313, 1 (nrow, ncol, nlyr)
# resolution : 10458.57, 10458.57 (x, y)
plet (r_3857, tiles = "Esri.WorldImagery" )