Explore GEBCO Bathymetry

1 GEBCO - General Bathymetric Chart of the Oceans

Download:

• GEBCO Gridded Bathymetry Data
  GEBCO_2022_sub_ice_topo.nc

Code

# packages ----
librarian::shelf(
  fs, glue, leaflet, leaflet.extras,
  marinesensitivity/msens,
  sf, terra, yaml,
  quiet = T)
# devtools::load_all("~/Github/MarineSensitivity/msens")

g_nc     <- "/Users/bbest/big/gebco_2022_sub_ice_topo/GEBCO_2022_sub_ice_topo.nc"
dir_data <- "/Users/bbest/My Drive/projects/msens/data"
g_tif    <- glue("{dir_data}/derived/gebco_depth.tif")

if (!file.exists(g_tif)){

  # prep shelfs polygon ----
  p <- msens::ply_shlfs
  b <- st_bbox(p)
  v <- vect(p)
  
  # read GEBCO netcdf ----
  r_g <- rast(g_nc) |>            # extent: -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    rotate(r_g, left=F) |>        # extent:    0, 360, -90, 90
    crop(r_gr, b) |>              # crop to bounding box
    mask(r_grc, v, touches = T)   # mask to shelfs polygon
  
  # convert elevation to depth ----
  r_depth <- r_g * -1
  names(r_depth) <- "depth"
  crs(r_depth) <- "epsg:4326"
  
  # write tif as COG ----
  writeRaster(
    r_depth, 
    g_tif, 
    overwrite = T,
    datatype  = "INT2S",
    gdal      = c(
      "TILED=YES",
      "COMPRESS=DEFLATE"))
}

# break up raster to either side of dateline ----

rast_to_cog3857eastwest <- function(
    r,  # raster or path to tif; must be True: crs(r) == crs("epsg:4326")
    cog_pfx       = "/Users/bbest/My Drive/projects/msens/data/derived/gebco_depth",
    cog_east_sfx  = "_cog3857e",
    cog_west_sfx  = "_cog3857w",
    datatype      = "INT2S", 
    method        = "bilinear"){
  
  if (is.character(r))
    r <- rast(r) # for analysis
  
  stopifnot(crs(r) == crs("epsg:4326"))
  
  # output each {h}eading
  for (h in c("east","west")){ # h = "west"
    message(glue("shift, crop, project {h} of antimeridian"))
    
    if (h == "east"){
      sfx <- cog_east_sfx
      r_h <- shift(r, -360) |> 
        crop(ext(-180,180,-90,90), snap="in") |>
        project("epsg:3857", method = method)
    }
    if (h == "west"){
      sfx <- cog_west_sfx
      r_h <- r |> 
        crop(ext(-180,180,-90,90), snap="in") |>
        project("epsg:3857", method = method)
    }
    
    tmp <- tempfile(fileext = ".tif")
    cog <- glue("{cog_pfx}{sfx}.tif")
    
    r_h |> 
      writeRaster(
        tmp, 
        overwrite = T,
        datatype  = datatype,
        gdal      = c(
          "TILED=YES",
          "COMPRESS=DEFLATE"))
    
    system(glue("rio cogeo create '{tmp}' '{cog}'"))
    unlink(tmp)
  }
  
  list(
    range       = global(r, "range", na.rm=T),
    bbox        = ext(r) |> st_bbox(),
    tif4326     = basename(glue("{cog_pfx}.tif")),
    cog3857east = basename(glue("{cog_pfx}{cog_east_sfx}.tif")),
    cog3857west = basename(glue("{cog_pfx}{cog_west_sfx}.tif"))) |> 
    write_yaml(
      glue("{cog_pfx}.yml"))
}

cog_pfx <- glue("{dir_data}/derived/gebco_depth")
tifs_ew <- glue("{cog_pfx}_cog3857{c('e','w')}.tif")
if(!all(file.exists(tifs_ew))){
  r <- rast(g_tif) 
  rast_to_cog3857eastwest(r, cog_pfx)
}

cog_base       = "https://file.marinesensitivity.org/cog/env/"
cog_palette    = "spectral_r"
cog_method     = "average"
cog_opacity    = 0.9
lgnd_palette   = "Spectral"
lgnd_palette_r = TRUE
title          = "Depth (m)"

m <- read_yaml(glue("{cog_base}/gebco_depth.yml"))

cog_range <- m$range |> as.numeric()
e_cog_url <- glue("{cog_base}/{m$cog3857east}")
w_cog_url <- glue("{cog_base}/{m$cog3857west}")

tile_opts <- glue(
  "resampling_method={cog_method}&rescale={paste(cog_range, collapse=',')}&return_mask=true&colormap_name={cog_palette}")

e_tile_url <- glue(
  "https://api.cogeo.xyz/cog/tiles/WebMercatorQuad/{{z}}/{{x}}/{{y}}@2x?url={e_cog_url}&{tile_opts}")
w_tile_url <- glue(
  "https://api.cogeo.xyz/cog/tiles/WebMercatorQuad/{{z}}/{{x}}/{{y}}@2x?url={w_cog_url}&{tile_opts}")

bb <- m$bbox

msens::ms_basemap() |>
  # West
  addTiles(
    urlTemplate = w_tile_url,
    options     = tileOptions(
      opacity   = cog_opacity)) |>
  # East
  addTiles(
    urlTemplate = e_tile_url,
    options     = tileOptions(
      opacity   = cog_opacity)) |>
  fitBounds(bb[[1]], bb[[2]], bb[[3]], bb[[4]]) |>
  addLegend(
    pal    = colorNumeric(lgnd_palette, cog_range[1]:cog_range[2], reverse = lgnd_palette_r),
    values = c(cog_range[1], cog_range[2]),
    title  = title) |> 
  leaflet.extras::addFullscreenControl()

2 VGPM

TODO:

• align VGPM with GEBCO
• convert to COG

Code

rast_to_cog3857eastwest(
  r       = glue("{dir_data}/derived/vgpm_2021.tif"), 
  cog_pfx = glue("{dir_data}/derived/vgpm_2021"),
  datatype = "FLT4S")