Downscaling AquaMaps

v01: blue whale, GEBCO SoCal

1 Overview

Goal: Downscale AquaMaps.org species distributions (Kaschner et al. 2023; Ready et al. 2010) from 0.5 decimal degrees to 15 arc seconds (111.11 km to 0.46 km at the equator), using the R package aquamapsdata and the the General Bathymetric Chart of the Oceans GEBCO.

We start with the “Blue Whale” (Balaenoptera musculus) and Southern California.

Later we’ll iterate over species and expand to global, which will require large raster handling techniques using Cloud-Optimized GeoTIFFS (COGs; see cogeo.org).

All code and files (except the large global GEBCO grid) are found in this repository:

• github.com/marinebon/aquamaps-downscaled

Code

# packages ----
if (!"librarian" %in% installed.packages())
  install.packages("librarian")
if (!"rcrypt" %in% installed.packages())
  devtools::install_bitbucket("bklamer/rcrypt") # dependency for aquamapsdata
librarian::shelf(
 bklamer/rcrypt,
 raquamaps/aquamapsdata,
 dplyr, ggplot2, glue, here, knitr, leaflet, 
 # TODO: migrate raster to terra
 # terra, 
 raster, rnaturalearth, sf, stringr, tidyr,
 quiet = T)
select = dplyr::select

# initial run-once step required to install remote db locally
# download_db(force = TRUE)

# aquamaps database ----
am_db <- default_db("sqlite")

# paths ----
dir_big         <- "/Users/bbest/big"
gebco_nc        <- glue("{dir_big}/gebco_2022_sub_ice_topo/GEBCO_2022_sub_ice_topo.nc")
gebco_socal_tif <- here("data/gebco_socal.tif")
land_socal_geo  <- here("data/land_socal.geojson")
bo_tif          <- here("data/bio-oracle.tif")

# custom functions ----
add_ocean_basemap <- function(m){
  # m: leaflet() map
  
  m |>
    # add base: blue bathymetry and light brown/green topography
    addProviderTiles(
      "Esri.OceanBasemap",
      options = providerTileOptions(
        variant = "Ocean/World_Ocean_Base")) |>
    # add reference: placename labels and borders
    addProviderTiles(
      "Esri.OceanBasemap",
      options = providerTileOptions(
        variant = "Ocean/World_Ocean_Reference"))
}

add_am_raster <- function(
    m, 
    r,
    title,
    cols = c("#FEB24C", "#FD8D3C", "#FC4E2A", "#E31A1C", "#B10026"),
    truncate_to_zero = T){
  # m: leaflet() map
  # r: raster 
  # TODO: migrate to terra::rast()
  
  # r     = r_gebco_bb
  # title = "GEBCO depth (m)"
  # cols  = RColorBrewer::brewer.pal(7, "Blues")
  
  r   <- leaflet::projectRasterForLeaflet(r, method = "bilinear")
  
  # truncate to 0 to prevent negative values 
  #   that were generated by projecting the raster 
  #   from geographic projection (decimal degrees) to Web Mercator (meters)
  if (truncate_to_zero){
    v <- values(r)
    v[v<0] <- 0
    values(r) <- v
  }
  
  pal <- leaflet::colorBin(
    cols, na.omit(unique(values(r))), 
    bins = length(cols), pretty = TRUE, na.color = "#00000000")
  
  e <- raster::extent(r) |> 
    sf::st_bbox() |> 
    st_as_sfc() |> 
    st_as_sf(crs=3857) |> 
    st_transform(4326) |> 
    st_bbox()
    
  m |> 
    leaflet::addRasterImage(
      r, project = F, colors = pal, opacity = 0.8) |> 
    addLegend(
      values = raster::values(r), 
      title = title, pal = pal) |> 
    leaflet::fitBounds(
      lng1 = e[["xmin"]], 
      lat1 = e[["ymin"]], 
      lng2 = e[["xmax"]], 
      lat2 = e[["ymax"]])
}

2 Species map (blue whale)

Code

# fuzzy search allows full text search operators AND, OR, NOT and +
# see https://www.sqlitetutorial.net/sqlite-full-text-search/
sp_term <- "blue whale"
key <- am_search_fuzzy(search_term = sp_term) |> 
  pull(key) # "ITS-Mam-180528"

# get the identifier for the species
r <- am_raster(key)

# show the native habitat map
m <- leaflet() |> 
  add_ocean_basemap() |> 
  add_am_raster(r, title = sp_term)
m 

Figure 1: Map of blue whale (Balaenoptera musculus) distribution from AquaMaps.

2.1 Zoom to SoCal

Notice the very large pixels, far bigger than useful for smaller planning purposes, such as for Sanctuaries or BOEM Wind Energy Areas.

Code

# Southern California
bbox <- c(-121, 32, -117, 35)
m |>
  fitBounds(lng1 = bbox[1], lat1 = bbox[2], lng2 = bbox[3], lat2 = bbox[4])

Figure 2: Map of blue whale (Balaenoptera musculus) distribution from AquaMaps zoomed into Southern California. Notice the very large pixels, far bigger than useful for smaller planning purposes, such as for Sanctuaries or BOEM Wind Energy Areas.

3 Environmental preferences

Here are the environmental preferences for the species in the database.

Code

sp_env <- am_hspen() |> 
  filter(SpeciesID == key) |> 
  head(1) |> 
  collect()

sp_env |> 
  mutate(across(everything(), as.character)) |> 
  pivot_longer(everything()) |> 
  kable()

Table 1: Table of blue whale (Balaenoptera musculus) environmental suitability parameters from Aquamaps.

name	value
SpeciesID	ITS-Mam-180528
Speccode	69007
LifeStage	adults
FAOAreas	18, 21, 27, 31, 34, 41, 47, 48, 51, 57, 58, 61, 67, 71, 77, 81, 87, 88
FAOComplete	NA
NMostLat	90
SMostLat	-90
WMostLong	-180
EMostLong	180
DepthYN	1
DepthMin	0
DepthPrefMin	1000
DepthPrefMax	4000
DepthMax	8000
MeanDepth	1
Pelagic	0
TempYN	1
TempMin	-1.8
TempPrefMin	-1.3
TempPrefMax	27.87
TempMax	32.07
SalinityYN	1
SalinityMin	3.58
SalinityPrefMin	32.57
SalinityPrefMax	35.49
SalinityMax	38.84
PrimProdYN	1
PrimProdMin	0.1
PrimProdPrefMin	1.4
PrimProdPrefMax	16.07
PrimProdMax	119.58
IceConYN	1
IceConMin	-0.88
IceConPrefMin	0
IceConPrefMax	0.49
IceConMax	0.96
OxyYN	0
OxyMin	1.1
OxyPrefMin	116.82
OxyPrefMax	275.01
OxyMax	408.48
LandDistYN	0
LandDistMin	0
LandDistPrefMin	17
LandDistPrefMax	733
LandDistMax	1740
Remark	FAO areas,bounding box and/or pelagic flag based on last review.
DateCreated	2019-08-07 00:00:00
DateModified	NA
expert_id	NA
DateExpert	NA
Layer	s
Rank	1
MapOpt	1
ExtnRuleYN	1
Reviewed	1

Now let’s convert all variables having {Var}YN == 1 into the relative environmental suitability rhomboids (Kaschner et al. 2006).

Code

var <- "Depth"

d_probs <- tribble(
  ~prob_name, ~prob_value,
  "Min"      , 0,
  "PrefMin"  , 1,
  "PrefMax"  , 1,
  "Max"      , 0)

vars_yes <- sp_env |> 
  select(ends_with("YN")) |> 
  pivot_longer(
    everything()) |> 
  filter(value == 1) |> 
  pull(name) |> 
  str_replace("YN","")

d <- sp_env |> 
  select(starts_with(vars_yes)) |>
  select(!ends_with("YN")) |> 
  pivot_longer(
    everything(),
    values_to = "var_value") |> 
  separate_wider_regex(
    name,
    c(var       = paste(vars_yes, collapse = "|"), # "",
      prob_name = paste(d_probs$prob_name, collapse = "|"))) |> 
  left_join(
    d_probs,
    by = "prob_name")

kable(d)

Table 2: Table environmental suitability parameters from Aquamaps that are applicable to blue whale (Balaenoptera musculus), i.e. {Var}YN == 1 in Table 1.

var	prob_name	var_value	prob_value
Depth	Min	0.00	0
Depth	PrefMin	1000.00	1
Depth	PrefMax	4000.00	1
Depth	Max	8000.00	0
Temp	Min	-1.80	0
Temp	PrefMin	-1.30	1
Temp	PrefMax	27.87	1
Temp	Max	32.07	0
Salinity	Min	3.58	0
Salinity	PrefMin	32.57	1
Salinity	PrefMax	35.49	1
Salinity	Max	38.84	0
PrimProd	Min	0.10	0
PrimProd	PrefMin	1.40	1
PrimProd	PrefMax	16.07	1
PrimProd	Max	119.58	0
IceCon	Min	-0.88	0
IceCon	PrefMin	0.00	1
IceCon	PrefMax	0.49	1
IceCon	Max	0.96	0

Code

g <- ggplot(d, aes(var_value, prob_value)) +
  geom_area() +
  theme_light() +
  facet_wrap(
    vars(var), 
    scales = "free") +
  labs(
    title    = sp_term,
    subtitle = "environmental envelope",
    x        = NULL,
    y        = "probability of presence")
g

Figure 3: Plots of environmental suitability parameters from Aquamaps that are applicable to blue whale (Balaenoptera musculus) from Table 2.

4 Depth (GEBCO) for SoCal

Code

# limit to bounding box for now
ply_bb <- extent(
  c(bbox[1], bbox[3], bbox[2], bbox[4])) |> 
  st_bbox() |> 
  st_as_sfc() |> 
  st_as_sf(crs = 4326)

# land
if (!file.exists(land_socal_geo)){
  ply_land <- ne_download(
    scale       = 10, # 110/50/10: high spatial resolution (10 m)
    type        = "land", 
    category    = "physical",
    returnclass = "sf")
  # plot(ply_land)
  ply_land_bb <- ply_land |> 
    st_intersection(ply_bb)
  # plot(ply_land_bb)
  write_sf(ply_land_bb, land_socal_geo)
}
ply_land_bb <- read_sf(land_socal_geo)
# plot(ply_land_bb)

if (!file.exists(gebco_socal_tif)){
  # read large GEBCO netcdf file outside repo
  r_gebco <- raster(gebco_nc)
  
  # crop to SoCal bounding box
  r_gebco_bb <- r_gebco |> 
    crop(ply_bb)
  
  # mask out land, ie > 0 
  r_gebco_bb <- r_gebco_bb |> 
    mask(r_gebco_bb <= 0, maskvalue = 0) * -1
  
  # write to TIF
  writeRaster(r_gebco_bb, gebco_socal_tif, overwrite = T)
}
r_gebco_bb <- raster(gebco_socal_tif)

m <- leaflet() |> 
  add_ocean_basemap() |> 
  add_am_raster(
    r                = r_gebco_bb, 
    title            = "GEBCO depth (m)", 
    cols             = RColorBrewer::brewer.pal(7, "Blues"))
m 

Figure 4: Map of depth from GEBCO zoomed into Southern California. Notice the much higher resolution compared to Figure 2.

5 Ramp depth with species preference

5.1 Create ramp_env() function

Code

ramp_env <- function(v, min, min_pref, max, max_pref){
  
  x <- c(min, min_pref, max, max_pref)
  y <- c(  0,        1,   1,        0)

  approx(
    x, y, 
    xout   = v, 
    yleft  = 0,
    yright = 0,
    rule   = 2,
    method = "linear")$y
}

p <- d |> 
  filter(var == !!var)
p <- setNames(p$var_value, p$prob_name) |> as.list()
# p
#   Min PrefMin PrefMax     Max 
#     0    1000    4000    8000

x_pref <- c(p$Min, p$PrefMin, p$PrefMax, p$Max)
y_pref <- c(  0,        1,   1,        0)
x_new <- seq(-200, 10000, by=100)
y_new <- ramp_env(x_new, p$Min, p$PrefMin, p$PrefMax, p$Max)

plot(
  x_pref,
  y_pref, 
  xlim = range(c(x_pref, x_new), na.rm=T), 
  ylim = range(c(y_pref, y_new), na.rm=T),
  xlab = var,
  ylab = sp_term)
points(x_new, y_new, col = 2, pch = "*")

Figure 5: Plot of original depth preferences for 4 points (black circles) and interpolated values (red asterisks) using new ramp_env() function.

5.2 Apply to SoCal

Apply the ramp_env() function to the SoCal depth using blue whale preferences.

Code

r_sp_depth_bb <- terra::app(
  x        = terra::rast(r_gebco_bb), 
  fun      = ramp_env, 
  min      = p$Min, 
  min_pref = p$PrefMin, 
  max_pref = p$PrefMax,
  max      = p$Max) |> 
  raster()

m <- leaflet() |> 
  add_ocean_basemap() |> 
  add_am_raster(
    r                = r_sp_depth_bb, 
    title            = glue("{sp_term}, {var} only"))
m 

Figure 6: Map of depth preference for r sp_term applied to SoCal depth with the ramp_env() function.

6 sdmpredictors

Code

librarian::shelf(
  sdmpredictors, skimr)

# exploring the marine datasets 
datasets <- list_datasets(terrestrial = FALSE, marine = TRUE)

kable(datasets)

	dataset_code	terrestrial	marine	url	description	citation
2	Bio-ORACLE	FALSE	TRUE	https://bio-oracle.org/	Bio-ORACLE is a set of GIS rasters providing geophysical, biotic and environmental data for surface and benthic marine realms at a spatial resolution 5 arcmin (9.2 km) in the ESRI ascii and tif format.	Tyberghein L., Verbruggen H., Pauly K., Troupin C., Mineur F. & De Clerck O. Bio-ORACLE: a global environmental dataset for marine species distribution modeling. Global Ecology and Biogeography. doi: 10.1111/j.1466-8238.2011.00656.x
3	MARSPEC	FALSE	TRUE	http://marspec.org/	MARSPEC is a set of high resolution climatic and geophysical GIS data layers for the world ocean. Seven geophysical variables were derived from the SRTM30_PLUS high resolution bathymetry dataset. These layers characterize the horizontal orientation (aspect), slope, and curvature of the seafloor and the distance from shore. Ten “bioclimatic” variables were derived from NOAA’s World Ocean Atlas and NASA’s MODIS satellite imagery and characterize the inter-annual means, extremes, and variances in sea surface temperature and salinity. These variables will be useful to those interested in the spatial ecology of marine shallow-water and surface-associated pelagic organisms across the globe. Note that, in contrary to the original MARSPEC, all layers have unscaled values.	Sbrocco, EJ and Barber, PH (2013) MARSPEC: Ocean climate layers for marine spatial ecology. Ecology 94: 979. doi: 10.1890/12-1358.1

Code

# exploring the marine layers 
layers <- list_layers(datasets)

# names(layers)

# skim(layers)

# table(layers$dataset_code)
# Bio-ORACLE    MARSPEC 
#        918         42 
# table(layers$primary_type)
#                                             ''                     GEBCO / EMODnet Bathymetry 
#                                              1                                              3 
#                            in situ measurement    in situ measurements, monthly climatologies 
#                                              7                                             17 
#                                          Model  Satellite (Aqua-MODIS), monthly climatologies 
#                                            831                                             31 
# Satellite (Aqua-MODIS), seasonal climatologies     Satellite (SeaWIFS), monthly climatologies 
#                                              2                                              4 
#                               Satellite (SRTM)        Satellite (Terra-MODIS), monthly images 
#                                              7                                              6 
#                              satellite imagery 
#                                             51 

# table(layers$units)
#               %         Celsius         degrees degrees Celcius      E/m^2/year Einstein/m_/day        fraction 
#               6              21               2              72              54               4              18 
#       g/m^3/day      kilometers               m             m/s            m^-1          meters          mg/m^3 
#              72               1              18              72               6               4              76 
#      micromol/L    micromol/m^3            ml/l         mol/m^3             PSS             psu         radians 
#               3             432               1               2              73              17               2 
#        unitless 
#               4

# layers |> 
  # select(
  #   layer_code, name, description, units, 
  #   primary_spatial_resolution, primary_source) |> 
  # filter(str_detect(description, regex("roductivity", ignore_case = T))) |> 
  # filter(str_detect(description, regex("Mean sea surface net primary productivity", ignore_case = T))) |> 
  # filter(str_detect(description, regex("Sea surface temperature (mean)", ignore_case = T))) |> 
  # filter(str_detect(description, regex("salinity", ignore_case = T))) |> 
  # filter(str_detect(description, regex("ice", ignore_case = T))) |> 
  # View()


# https://aquamaps.org/main/files/HCAF_v7.zip
# 
# - PrimProdMean: Original unit in gC·m-3·day-1, converted to mgC·m-3·day-1. If OceanArea=0 (i.e. land cell), non-zero values were set to null (53 cells). Bio-ORACLE grid only goes down to 78.5 S, 1,869 ocean cells tagged -9999.

# layers |> 
#   filter(layer_code == "BO22_ppmean_ss") |> 
#   mutate(across(everything(), as.character)) |> 
#   pivot_longer(everything()) |> 
#   kable()
# units: g/m^3/day    

# get Bio-Oracle (BO) latest version (22)
lyrs <- c(
  Temp     = "BO22_tempmean_ss",
  Salinity = "BO22_salinitymean_ss",
  PrimProd = "BO22_ppmean_ss",
  IceCon   = "BO22_icecovermean_ss")
  
# download to temporary directory
if (!file.exists(bo_tif)){

  stk <- load_layers(lyrs, datadir = here("data/bio-oracle/temp"))

  # class      : RasterStack 
  # dimensions : 2160, 4320, 9331200, 4  (nrow, ncol, ncell, nlayers)
  # resolution : 0.08333333, 0.08333333  (x, y)
  # extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
  # crs        : +proj=longlat +datum=WGS84 +no_defs 
  # names      : BO22_tempmean_ss, BO22_salinitymean_ss, BO22_ppmean_ss, BO22_icecovermean_ss 
  # min values :        -1.797733,             0.059304,       0.000060,             0.000000 
  # max values :        30.178629,            40.651300,       0.257881,             0.974730

  # get x arg for [ee.Image.interpolate()](https://developers.google.com/earth-engine/apidocs/ee-image-interpolate)
  # d |> 
  #   # filter(var == "Temp") |>       # -1.8, -1.3, 27.87, 32.07
  #   # filter(var == "Salinity") |>   # 3.58, 32.57, 35.49, 38.84
  #   # filter(var == "PrimProd") |>   # 0.1, 1.4, 16.07, 119.58
  #   filter(var == "IceCon") |>       # -0.88, 0, 0.49, 0.96
  #   pull(var_value) |> 
  #   paste(collapse=", ") |> 
  #   cat()
  
  raster::writeRaster(stk, bo_tif)
}

# librarian::shelf(terra)
# plet(rast(bo_tif)[[1]])

7 manual upload with retry

Code

librarian::shelf(googleCloudStorageR)

gcs_auth_json <- "/Users/bbest/My Drive/private/offhab-google-service-account_09e7228ac965.json"
gcs_bucket    <- "offhab_lyrs"
name          <- "bio-oracle.tif"

Sys.setenv(
  "GCS_DEFAULT_BUCKET" = gcs_bucket,
  "GCS_AUTH_FILE"      = gcs_auth_json)

u <- googleCloudStorageR::gcs_upload(
  file          = bo_tif,
  bucket        = "offhab_lyrs",
  name          = basename(bo_tif))
# ℹ 2023-07-30 18:38:57.609799 > Found resumeable upload URL:  https://www.googleapis.com/upload/storage/v1/b/offhab_lyrs/o/?uploadType=resumable&name=bio-oracle.tif&predefinedAcl=private&upload_id=ADPycduT4eJUJFQA3c-cdYj6U-3ljbVu6IItZxtQi3WyE68uSdBS_o6uuL1Pq8kJnFvr5jdZfJgcw9_ZPGmlZKzXC3C8Zw
# retry
u2 <- googleCloudStorageR::gcs_retry_upload(u)

# https://code.earthengine.google.com/?asset=projects/eq-am-fine/assets/sdmpredictors

gcs_to_gee <- function (
    gcs_name, 
    gee_name          = fs::path_ext_remove(gcs_name), 
    properties        = "", 
    gcs_bucket        = "offhab_lyrs", 
    gee_asset         = "projects/eq-am-fine/assets/sdmpredictors", 
    pyramiding_policy = "MEAN"
    # missing_data      = 255
    ){
  
    f_json <- tempfile(fileext = ".json")
    properties_json <- jsonlite::toJSON(properties, pretty = T,
        auto_unbox = T)
    # writeLines(glue::glue("{{\n       \"name\": \"{gee_asset}/{gee_name}\",\n       \"tilesets\":[{{\"sources\":[{{\"uris\":[\"gs://{gcs_bucket}/{gcs_name}\"]}}]}}],\n       \"pyramidingPolicy\":\"MEAN\",\n       \"properties\": {properties_json},\n       \"missing_data\":{{\"values\":[{missing_data}]}}\n    }}"), 
    writeLines(glue::glue("{{\n       \"name\": \"{gee_asset}/{gee_name}\",\n       \"tilesets\":[{{\"sources\":[{{\"uris\":[\"gs://{gcs_bucket}/{gcs_name}\"]}}]}}],\n       \"pyramidingPolicy\":\"MEAN\",\n       \"properties\": {properties_json} }}"), 
        f_json)
    cmd <- glue::glue("earthengine upload image --manifest '{f_json}'")
    message(cmd)
    system(cmd)
    return(glue("{gee_asset}/{gee_name}"))
}


gcs_to_gee(
  basename(bo_tif),
  properties = list(
    source  = "sdmpredictors"))

# readLines("/var/folders/sl/7s3zmk1129jcrgsn1c4hcs2r0000gn/T//Rtmp2X9EGE/file5883322f397b.json") |> cat()


# tif <- "/Users/bbest/Downloads/depth_m0000000000-0000000000.tif"
# tif <- "/Users/bbest/Downloads/depth_m0000023552-0000000000.tif"
# tif <- "/Users/bbest/Downloads/depth_m0000376832-0000000000.tif"
# r <- rast(tif)
# r
# terra::plet(r)

• problems here possibly in version 0.6-1 of stars > # version 0.6-1 > * [<-.stars_proxy() clones environment, so that after r[r > 100]<-NA we don’t get infinite recursion when realizing r

• General Bathymetric Chart of the Oceans (GEBCO) - awesome-gee-community-catalog

• Resampling and Reducing Resolution | Google Earth Engine | Google for Developers

Generate from ramps in GEE with:

• ee.Algorithms.If()
  • Computations using Images | Google Earth Engine | Google for Developers
  • javascript - Google Earth Engine: Apply scaling function to band values of image - Geographic Information Systems Stack Exchange
  • Applying Raster Calculation to Image Collection in Google Earth Engine - Geographic Information Systems Stack Exchange
• Image.unitScale(low, high)
  • google earth engine - Rescale NDVI [-1,1] to 0-255 using GEE - Geographic Information Systems Stack Exchange

8 Try COG

• GEE: export_depth

Code

# devtools::install_github("r-spatial/stars")
librarian::shelf(
  offshorewindhabitat/offhabr, glue, here, fs, stars, terra)

depth_tif <- glue("{path_ext_remove(gebco_nc)}_bathy.tif")
depth_cog <- glue("{path_ext_remove(gebco_nc)}_bathy_cog.tif")

# read large GEBCO netcdf file outside repo
r_gebco <- rast(gebco_nc)
# terra::writeRaster(r_gebco, depth_cog)

r_d190 <- read_stars(gebco_nc) |> 
  st_downsample(190)
r_d190[r_d190 > 0] = NA
plot(r_d190)

system.time({r <- rast(gebco_nc)}) # 10516.380 secs elapsed
system.time({r[r > 0] = NA})       # 5915.483  secs elapsed 
system.time({r <- r * -1})         # 190.878   secs elapsed
system.time({writeRaster(r, depth_tif, datatype = "INT2U")})

r <- rast(depth_tif)

system.time({
offhabr::write_rast(
  # r,
  depth_tif,
  depth_cog,
  datatype     = "INT2U",
  overwrite    = TRUE,
  method       = "average",
  threads      = "ALL_CPUS",
  epsg         = 4326,
  verbose      = T,
  web_optimize = F,
  use_gdal_cog = T)          # system.time(): 2183.913 secs elapsed
})

# install.packages("googleCloudStorageR")
librarian::shelf(googleCloudStorageR)
# ✔ Setting scopes to https://www.googleapis.com/auth/devstorage.full_control and https://www.googleapis.com/auth/cloud-platform
# ✔ Successfully auto-authenticated via /Users/bbest/My Drive/private/offhab-google-service-account_09e7228ac965.json
# ✔ Set default bucket name to 'offhab_lyrs'
upload_to_gcs(
  file          = depth_cog,
  name          = "depth_gebco_m.tif",
  make_public   = T,
  gcs_auth_json =
    "/Users/bbest/My Drive/private/offhab-google-service-account_09e7228ac965.json",
  gcs_bucket    = "offhab_lyrs") # 3 GB at 3 MBPS: 2 hrs 23 mins
# ℹ 2023-07-30 14:58:35.94178 > File size detected as  3.1 Gb
# ℹ 2023-07-30 14:58:37.120278 > Found resumeable upload URL:  https://www.googleapis.com/upload/storage/v1/b/offhab_lyrs/o/?uploadType=resumable&name=depth_gebco_m.tif&predefinedAcl=private&upload_id=ADPycdun6IdndMKBjpKcpL2KvQV145Q6pzsFt5F7-c7ZdRgtmzTH_-sSe9JfKUebpJiKZib-a6MCKYyGFzNh9hYABimF

# manual upload with retry
gcs_auth_json <- "/Users/bbest/My Drive/private/offhab-google-service-account_09e7228ac965.json"
gcs_bucket    <- "offhab_lyrs"
name          <- "depth_gebco_m.tif"
Sys.setenv(
  "GCS_DEFAULT_BUCKET" = gcs_bucket,
  "GCS_AUTH_FILE"      = gcs_auth_json)

u <- googleCloudStorageR::gcs_upload(
  file          = depth_cog,
  bucket        = "offhab_lyrs",
  name          = name)
Sys.time()

tif <- "/Users/bbest/Github/marinebon/aquamaps-downscaled/data/gebco_socal_cog.tif"
u <- googleCloudStorageR::gcs_upload(
  file          = tif,
  bucket        = "offhab_lyrs",
  name          = basename(tif))
# make publicly available
googleCloudStorageR::gcs_update_object_acl(
   basename(tif), entity_type = "allUsers")


# retry
u2 <- googleCloudStorageR::gcs_retry_upload(u)

# make publicly available
googleCloudStorageR::gcs_update_object_acl(
  name, entity_type = "allUsers")



Sys.time()

# TODO: upload depth_cog to GCStorage and try in leaflet

# OLD...
write_stars(st_g, depth_cog, type="Int16")

s <- read_stars(depth_cog)
plot(s)

st_gebco <- read_stars(
  depth_cog, 
  proxy=T)
plot(st_gebco < 0) # downsample set to 190

system.time({
st_gebco[st_gebco > 0] = NA  # system.time():  mins elapsed
})
plot(st_gebco, downsample=190) 
# Error in `[<-.stars`(x = x, i = i, value = value, downsample = downsample) : 
#   unused argument (downsample = downsample)
#   https://github.com/r-spatial/stars/issues/627
write_stars(st_gebco, depth_cog)


plot(st_gebco)

system.time({
st_gebco <- st_gebco * -1     # system.time(): 4.042083 mins elapsed
})

Code

librarian::shelf(
  gdalcubes, glue, here, fs, stars, terra)

gebco_socal_cog <-  glue("{path_ext_remove(gebco_socal_tif)}_cog.tif")
r <- read_stars(gebco_socal_tif)
d <- stars::st_dimensions(r)

# does not work:
# file_delete(gebco_socal_cog)
# write_stars(r, gebco_socal_cog, driver="COG", type="UInt16")

ic <- create_image_collection(
  files     = gebco_socal_tif, 
  date_time = as.Date("1970-01-01"))

cv <- cube_view(
  extent      = extent(ic,"EPSG:4326"), 
  srs         = "EPSG:4326", 
  dt          = "P1Y", 
  nx          = d$x$to, 
  ny          = d$y$to, 
  aggregation = "mean", 
  resampling  = "bilinear")

write_tif(
  raster_cube(ic, cv),
  dir              = here("data"),
  prefix           = "cube_",
  overviews        = T,
  COG              = T,
  rsmpl_overview   = "bilinear",
  creation_options = NULL,
  write_json_descr = T)


# cube_1970-01-01.tif

Code

librarian::shelf(
  r-spatial/leafem, leaflet, stars, viridisLite)

url <- "https://storage.googleapis.com/offhab_lyrs/cube_1970-01-01.tif"
r <- rast(glue("/vsicurl/{url}"))

pal = hcl.colors(256, "inferno")
brks = NULL

myCustomJSFunc = htmlwidgets::JS(
  "
    pixelValuesToColorFn = (raster, colorOptions) => {
      const cols = colorOptions.palette;
      var scale = chroma.scale(cols);

      if (colorOptions.breaks !== null) {
        scale = scale.classes(colorOptions.breaks);
      }
      var pixelFunc = values => {
        let clr = scale.domain([raster.mins, raster.maxs]);
        if (isNaN(values)) return colorOptions.naColor;
        if (values < 50) return chroma('pink').alpha(0.7).hex();
        return clr(values).hex();
      };
      return pixelFunc;
    };
  "
)

leaflet() |>
  addTiles() |>
  addCOG(
    url          = url, 
    layerId      = "band1",
    group        = "COG", 
    resolution   = 512, 
    opacity      = 0.9,
    autozoom     = TRUE,
    colorOptions = colorOptions(
      palette = viridisLite::inferno,
      breaks  = seq(1,5349, by=500),
      na.color = "transparent"))
    # colorOptions = leafem:::colorOptions(
    #   palette  = pal, #viridisLite::cividis
    #   # breaks   = brks,
    #   na.color = "transparent"),
    # pixelValuesToColorFn = myCustomJSFunc)

Code

librarian::shelf(offhabr)

base_opacity   = 0.5

cog_range      = c(1, 5349)
cog_method     = "average"
cog_palette    = "spectral_r"
cog_opacity    = 0.9
# lgnd_palette   = "Spectral"
# lgnd_palette_r = TRUE

cog_url   <- "https://storage.googleapis.com/offhab_lyrs/cube_1970-01-01.tif"
tile_opts <- glue(
  "resampling_method={cog_method}&rescale={paste(cog_range, collapse=',')}&return_mask=true&colormap_name={cog_palette}")
tile_url  <- glue(
  "https://api.cogeo.xyz/cog/tiles/WebMercatorQuad/{{z}}/{{x}}/{{y}}@2x?url={cog_url}&{tile_opts}")

# oh_map_cog_lyr("cube_1970-01-01")

oh_map() |>
  addTiles(
    urlTemplate=tile_url,
    options = tileOptions(
      opacity = cog_opacity))

• How do I set up CORS for my Google Cloud Storage Bucket?

echo '[{"origin": ["*"],"responseHeader": ["*"],"method": ["GET", "HEAD"],"maxAgeSeconds": 3600}]' > cors-config.json

echo '[{"origin": ["*"],"method": ["*"]}]' > cors-config.json

gsutil cors set cors-config.json gs://offhab_lyrs

gsutil cors get gs://offhab_lyrs

9 Next Steps

Goal: Downscale global AquaMaps with all env preferences

TODO:

• Gather finer resolution global data based on yes/no parameters (*YN):

  • Depth: GEBCO

  • Temp: BO22_tempmean_ss

  • Salinity: BO22_salinitymean_ss

  • PrimProd: BO22_ppmean_ss

  • IceCon: BO22_icecovermean_ss

  • Oxy:

  • LandDist:

  • ExtnRule:

    • Note: “BO22_” refers to Bio-Oracle version 2.2 from sdmpredictors

• Work out rest of individual species workflow

  • Apply ramp_env() to each environmental parameter applicable to the species (ie {var}YN ==1)
  • Average all environmental parameter grids for the species
  • Clip to NMostLat, SMostLat, WMostLong, EMostLong
  • Mask to FAOAreas
  • Figure out unknown fields: ExtnRuleYN, MapOpt, …

• Repeat for workflow for ALL species

• Render maps from anywhere

  • Store in COG with each species as a separate layer; or try individual layer COGs
  • Upload COG to Google Cloud Storage
  • Install TiTiler on MarineBON.app server
  • Render map layers from COG using TiTiler in new function(s) borrowing from offhabr functions like oh_map_cog_lyr()

• Build Shiny app

  • Render map of selected species from dropdown
  • Summarize species from drawn area
  • Summarize species from selected area from existing:
    • EEZ
    • LME
    • Sanctuary
    • BOEM Wind Energy Area
    • …

• Calculate biodiversity metrics

  • Richness
  • Abundance
  • Extinction Risk
  • Endemism
  • Foundation Species
  • …

10 Links from GEE

TODO:

• convert Bio-Oracle SDM predictors to same resolution
• experiment with mosaic from for simple export
  • google earth engine - How to split the world in four quadrants as ee.Image - Geographic Information Systems Stack Exchange
  • Google Earth Engine Tutorial: Split Image By Grid and Export to Google Drive - YouTube
  • ee.ImageCollection.mosaic | Google Earth Engine | Google for Developers
  • ee.ImageCollection.fromImages | Google Earth Engine | Google for Developers
  • Image mosaic/composite creation for Landsat and Sentinel-2 in Google Earth Engine - openMRV

Try:

• /Users/bbest/Github/offshorewindhabitat/scripts/offhab_gee1.ipynb
• geemap: see below
• zarr | IOOS
• stars
  • stars proxy objects
  • IPBES raster gdalcubes
  • Large data and cloud native | Ch. 9 Spatial Data Science
  • Accessing data from large online rasters with Cloud-Optimized-Geotiff, GDAL, and terra R package | Francisco Rodríguez-Sánchez
  • Cloud-based processing of satellite image collections in R using STAC, COGs, and on-demand data cubes
• viz with leafem: seemed slow
  • new addGeoRaster method · Issue #25 · r-spatial/leafem
  • Add Cloud Optimised Geotiff (COG) to a leaflet map. — addCOG • leafem
• gdalcubes again
  • w/ rstac
  • break up into smaller tiles
• 2.2. Generate a Regional Composite Through Spatial Tiling
• Compositing, Masking, and Mosaicking | GEE
  By combining the concepts of image collections, logical operators, masking and compositing, you can achieve interesting cartographic results. For example, suppose you want an image in which land pixels are displayed in true-color and all the other pixels are displayed in blue
• other env predictors
  • distance from shore
    • 30-m global shorelin on GEE
    • ee.Image.distance | Google Earth Engine | Google for Developers
    • ee.Image.fastDistanceTransform | Google Earth Engine | Google for Developers
  • SST by month
    • GEE: GSST 2002-2019
• model data products
  • ECCO
    • https://podaac.jpl.nasa.gov/dataset/ECCO_L4_TEMP_SALINITY_05DEG_MONTHLY_V4R4#
    • https://ecco-v4-python-tutorial.readthedocs.io/intro.html
• install titiler
  • https://github.com/developmentseed/titiler
• gee tricks
  • client vs server | GEE
  • GFW use of GEE
  • Jenks natural breaks w/ SLD
• geemap
  • 11 export image
  • 44 cog stac
  • 92 plotly - geemap
  • 95 create cog
  • 100 numpy to cog
  • 103 split control
  • plotlymap module add_cog_layer()
  • plotlymap module add_stac_layer(
  • foliumap module add_cog_layer()
  • foliumap add_stac_layer()
• indicator portal
  • Data Management Tutorials - IPBES ICT guide
    • Part 2 - Preparing and Mapping Data to IPBES Regions and Sub-regions - IPBES ICT guide
    • Part 11 - How to Document an Indicator - IPBES ICT guide
• Map of Life: indicators, portal setup
  • Map of Life | Map of Life
  • Map of Life - Indicators
  • Map of Life - Patterns
  • Discovery Potential

11 References

Kaschner, K., K. Kesner-Reyes, C. Garilao, J. Segschneider, J. Rius-Barile, T. Rees, and R. Froese. 2023. “AquaMaps: Predicted Range Maps for Aquatic Species. Retrieved from https://www.aquamaps.org.”

Kaschner, K., R. Watson, A. W. Trites, and D. Pauly. 2006. “Mapping World-Wide Distributions of Marine Mammal Species Using a Relative Environmental Suitability (RES) Model.” Marine Ecology Progress Series 316 (July): 285–310. https://doi.org/10.3354/meps316285.

Ready, Jonathan, Kristin Kaschner, Andy B. South, Paul D. Eastwood, Tony Rees, Josephine Rius, Eli Agbayani, Sven Kullander, and Rainer Froese. 2010. “Predicting the Distributions of Marine Organisms at the Global Scale.” Ecological Modelling 221 (3): 467–78. https://doi.org/10.1016/j.ecolmodel.2009.10.025.