Ingest Primary Productivity

Code

librarian::shelf(
  cmocean, DT, fs, glue, here, htmltools, janitor, jsonlite, leaflet, 
  leaflet.extras, mapview, rerddap, R.utils, sf, stringr, terra, units,
  quiet = T)

source(here("libs/db.R")) # define: database connection (con)

yr               <- 2023  # vgpm-viirs most recent complete year
url_tar          <- glue("http://orca.science.oregonstate.edu/data/1x2/monthly/vgpm.r2022.v.chl.v.sst/hdf/vgpm.v.{yr}.tar")
file_tar         <- glue("/Users/bbest/My Drive/projects/msens/data/raw/oregonstate.edu/vgpm.v.{yr}.tar")
dir_out          <- "/Users/bbest/My Drive/projects/msens/data/derived"
dir_raw          <- "/Users/bbest/My Drive/projects/msens/data/raw/boem.gov"
vgpm_tif         <- glue("{dir_out}/vgpm-viirs_{yr}-daily-avg.tif")
oa_new_geo       <- glue("{dir_raw}/OCS_Area_Polygons.geojson")
pa_new_geo       <- glue("{dir_raw}/OCS_Proposed_Planning_Area_Polygons.geojson")
pa_oa_geo        <- glue("{dir_out}/pa_oa_2025.geojson")
verbose          <- F

1 Vertically Generalized Production Model (VGPM), VIIRS 2023 | OregonState

Source:

Ocean Productivity from Oregon State
- Standara VGPM: 1080 by 2160 Monthly HDF files from VIIRS R2022
- Frequently Asked Questions

Calculated sum of 2023 monthly net primary production (NPP) from VIIRS data.

Community guidance for developing this website was to provide a single productivity product as a Standard product. For this, we initially chose the original Vertically Generalized Production Model (VGPM) (Behrenfeld and Falkowski, 1997a) as the standard algorithm. The VGPM is a “chlorophyll-based” model that estimate net primary production from chlorophyll using a temperature-dependent description of chlorophyll-specific photosynthetic efficiency. For the VGPM, net primary production is a function of chlorophyll, available light, and the photosynthetic efficiency.

Standard products are based on chlorophyll, temperature and PAR data from SeaWiFS, MODIS and VIIRS satellites, along with estimates of euphotic zone depth from a model developed by Morel and Berthon (1989) and based on chlorophyll concentration.

Units: mg C / m² / day

Citation: Behrenfeld and Falkowski, 1997a

Code

if (!file.exists(vgpm_tif)) {
  if(!file.exists(file_tar))
    download.file(url_tar, file_tar)
  
  dir_untar <- path_ext_remove(file_tar)
  dir.create(dir_untar, showWarnings = F)
  untar(file_tar, exdir = dir_untar)
  
  # Get list of compressed files
  files_gz <- dir_ls(dir_untar, glob = "*.gz")
  
  dir_unzip <- file.path(dir_untar, "unzipped")
  dir.create(dir_unzip, showWarnings = F)
  
  # Unzip all files
  for (file_gz in files_gz) { # file_gz = files_gz[1]
    file_hdf <- file.path(dir_unzip, sub("\\.gz$", "", basename(file_gz)))
    
    if (!file.exists(file_hdf))
      gunzip(file_gz, destname = file_hdf)
  }
  
  # Get list of unzipped HDF files
  files_hdf <- dir_ls(dir_unzip, glob = "*.hdf") # length(files_hdf) = 12
  
  # Function to read VGPM data from HDF file
  read_vgpm_hdf <- function(file_hdf) { # file_hdf <- files_hdf[1]
    # First, query the HDF file to see the subdatasets
    sds_info <- terra::describe(file_hdf)
    if (verbose)
      print(sds_info)
    
    tryCatch({
      r <- terra::rast(file_hdf, noflip=T)
      
      # Standard VGPM product is typically 1/6 degree global grid (-180 to 180, -90 to 90)
      ext(r) <- c(-180, 180, -90, 90)
      crs(r) <- "EPSG:4326"
      r[r==-9999] <- NA
      # plot(r)
      
      return(r)
    }, error = function(e) {
      warning("Error reading", file_hdf, ":", e$message, "\n")
      return(NULL)
    })
  }
  
  # Initialize a list to store all monthly rasters
  monthly_rasters <- list()
  
  # Read all monthly files
  for (i in seq_along(files_hdf)) {
    if (verbose)
      message("Processing file", i, "of", length(files_hdf), ":", files_hdf[i], "\n")
    r <- read_vgpm_hdf(files_hdf[i])
    
    if (!is.null(r)) {
      monthly_rasters[[i]] <- r
    }
  }
  
  # Create a SpatRaster collection from all the monthly rasters
  
  s_vgpm <- rast(monthly_rasters)
  
  r <- app(s_vgpm, mean, na.rm = TRUE)
  names(r) <- "npp_daily_avg"
  
  writeRaster(r, vgpm_tif, overwrite = TRUE)
  
  message(glue(
    "Annual products created and saved.
    - daily avg for {yr}: {vgpm_tif}"))
}

r <- rast(vgpm_tif) # |> 
  # rotate(left=F) # -180 to 180 -> 0 to 360
  # shift(dx = 0.0001) # problem with -180 extent

mapView(
  r, maxpixels = 2332800, col.regions = cmocean("algae"), 
  layer.name = "raster: daily avg NPP") |> 
  slot("map") |> 
  addFullscreenControl() |> 
  addControl(tags$div(HTML("units: mg C / m<sup>2</sup> / day")), position = "topright")

2 NEW Proposed Planning Areas

BOEM internal:

Code

set_geom_ctr_area <- function(x){
  x |>
    st_set_geometry("geom") |>
    mutate(
      ctr      = st_centroid(geom),
      ctr_lon  = ctr |> st_coordinates() %>% .[,"X"],
      ctr_lat  = ctr |> st_coordinates() %>% .[,"Y"],
      area_km2 = st_area(geom) |>
        set_units(km^2) |>
        as.numeric() ) |>
    select(-ctr)
}

drop_ctr_area <- function(x){
  x |>
    select(-ctr_lon, -ctr_lat, -area_km2)
}

if (!file.exists(pa_oa_geo)){
  pa0 <- st_read(con, "ply_planareas_s05") |> 
    st_drop_geometry() |> 
    select(planarea_key, planarea_name) |> 
    mutate(
      planarea_name = str_replace(planarea_name, "Gulf of Mexico", "Gulf of America"),
      planarea_key = str_replace(planarea_key, "pa", ""),
      planarea_key = case_match(
        planarea_key, 
        "CGM" ~ "CGA",
        "EGM" ~ "EGA",
        "WGM" ~ "WGA",
        .default = planarea_key))
  
  pa <- read_sf(pa_new_geo) |> 
    clean_names() |> 
    rename(
      region_name  = region,
      planarea_key = planning_area) |>
    mutate(
      region_name = str_replace(region_name, "Gulf of Mexico", "Gulf of America"),
      region_key = case_match(
        region_name,
        "Alaska"          ~ "AK",
        "Atlantic"        ~ "AT",
        "Gulf of America" ~ "GA",
        "Pacific"         ~ "PA"),
      planarea_key = case_match(
        planarea_key, 
        "CG" ~ "CGA",
        "EG" ~ "EGA",
        "WG" ~ "WGA",
        .default = planarea_key)) |> 
    left_join(
      pa0, by = "planarea_key") |> 
    mutate(
      planarea_name = case_match(
        planarea_key,
        "HAR" ~ "High Arctic",
        .default = planarea_name)) # TODO: planarea_name for HAR is ...?
  
  oa <- read_sf(oa_new_geo) |> 
    clean_names() |> 
    rename(
      region_name   = region,
      planarea_name = area_name) |> 
    filter(!planarea_name %in% c("Alaska", "Pacific Coast", "Atlantic Coast", "Gulf of Mexico")) |> 
    mutate(
      region_key = case_match(
        region_name,
        "Atlantic" ~ "AT",
        "Pacific"  ~ "PA"),
      planarea_key = case_match(
        planarea_name,
        "American Samoa"                                        ~ "AMS",
        "Guam and Commonwealth of the Northern Mariana Islands" ~ "GUA",
        "Hawaiian Islands and Midway Atoll"                     ~ "HAW",
        "Howland and Baker Islands"                             ~ "HOW",
        "Jarvis Island"                                         ~ "JAR",
        "Johnston Atoll"                                        ~ "JOH",
        "Palmyra Atoll and Kingman Reef"                        ~ "PAL",
        "Puerto Rico and U.S. Virgin Islands"                   ~ "PUR",
        "Wake Island"                                           ~ "WAK"))
  
  pa_oa <- bind_rows(
    pa,
    oa) |>
    select(region_key, region_name, planarea_key, planarea_name) |>
    st_make_valid() |> 
    set_geom_ctr_area() |> 
    st_shift_longitude() |> 
    arrange(region_key, planarea_key)
  
  write_sf(pa_oa, pa_oa_geo)
}

pa_oa <- read_sf(pa_oa_geo)

mapView(pa_oa)

Code

pa_oa |> 
  st_drop_geometry() |> 
  datatable()

3 VGPM in Proposed Planning Areas

Extract daily average net primary production (NPP) values across all cells within the proposed planning areas.

Code

r <- rast(vgpm_tif)
r_r <- r |> 
  rotate(left=F) # -180 to 180 -> 0 to 360
  # shift(dx = 0.0001) # problem with -180 extent

pa_oa_npp <- pa_oa |>
  mutate(
    npp_cellavg = zonal(r_r, pa_oa |> vect(), fun = "mean", na.rm = T) |>
      pull(npp_daily_avg))

(mapView(
  r, maxpixels = 2332800, 
  col.regions = cmocean("algae"), 
  layer.name = "raster") + 
    mapView(
      pa_oa_npp, zcol = "npp_cellavg", 
      col.regions = cmocean("algae"), 
      layer.name = "planning areas")) |> 
  slot("map") |> 
  addFullscreenControl() |> 
  addControl(tags$div(HTML("units: mg C / m<sup>2</sup> / day")), position = "topright")

Code

pa_oa_npp |> 
  st_drop_geometry() |> 
  datatable()

Code

# Latest
# 
# -   [ERDDAP - search "productivity_viirs_snpp"](https://coastwatch.pfeg.noaa.gov/erddap/search/index.html?page=1&itemsPerPage=1000&searchFor=productivity_viirs_snpp)
# 
# -   monthly has lots of cloud cover -- see [WMS map](https://coastwatch.pfeg.noaa.gov/erddap/wms/productivity_viirs_snpp_monthly/index.html)
# 
# -   [Arctic Ocean Primary Productivity: The Response of Marine Algae to Climate Warming and Sea Ice Decline - NOAA Arctic](https://arctic.noaa.gov/report-card/report-card-2023/arctic-ocean-primary-productivity-the-response-of-marine-algae-to-climate-warming-and-sea-ice-decline-2023/#methods_data) ![](https://arctic.noaa.gov/wp-content/uploads/2023/11/ARC23-PrimaryProductivity-frey-Fig3-878x1024.png)

dir_nc <- here("data/productivity")
dates  <- c("2023-01-16", "2023-12-16")

# https://tile.marinesensitivity.org/public.ply_boemrgns.html
p_br <- read_sf(con, "ply_boemrgns")

# chunk east and west of antimeridian for extraction ----
bb_w <- st_bbox(c(xmin = 0, xmax = 180, ymin = -90, ymax = 90)) |> 
  st_as_sfc() |> 
  st_as_sf(crs = 4326)

bb_e <- st_bbox(c(xmin = -180, xmax = 0, ymin = -90, ymax = 90)) |> 
  st_as_sfc() |> 
  st_as_sf(crs = 4326)

p_br <- bind_rows(
  p_br |> 
    st_intersection(bb_w) |> 
    mutate(
      chunk = "w"),
  p_br |> 
    st_difference(bb_w) |> 
    mutate(
      chunk = "e")) |> 
  relocate(chunk, .after = boemrgn_key) |> 
  arrange(boemrgn_key, chunk)

leaflet() |> 
  addProviderTiles(providers$Esri.OceanBasemap) |> 
  addPolygons(data = st_shift_longitude(p_br))

p_br |> 
  st_drop_geometry() |> 
  datatable()

Code

# iterate over polygon chunks ----
for (i in 1:nrow(p_br)){ # p = p_br[3,]
  
  p  <- slice(p_br, i)
  bb <- st_bbox(p)
  
  nc <- here(glue(
    "data/productivity/{p$boemrgn_key}{p$chunk}_{paste(dates, collapse = '_')}.nc"))
  
  if (file.exists(nc)){
    message(glue("{i}/{nrow(p_br)} skip: {basename(nc)} already exists"))
    next()
  }
  
  message(glue("{i}/{nrow(p_br)} fetch: {basename(nc)}  ~ {Sys.time()}"))
  o <- griddap(
    "productivity_viirs_snpp_monthly",
    url = "https://coastwatch.pfeg.noaa.gov/erddap/",
    time      = dates,
    longitude = c(max(bb$xmin, -179.9792), min(bb$xmax, 179.9792)), 
    # TODO: compare with actual min/max from fetching
    latitude  = c(bb$ymin, bb$ymax),
    fields    = "productivity",
    store     = disk(path = dir_nc))
  file.rename(o$summary$filename, nc)
}
message(glue("{nrow(p_br)}/{nrow(p_br)} done: {basename(nc)}  ~ {Sys.time()}"))

boemrgn_key = "brAK"
chunk       = "e"
p <- p_br |> 
  filter(
    boemrgn_key == !!boemrgn_key,
    chunk       == !!chunk)
nc <- here(glue("data/productivity/{boemrgn_key}{chunk}_{paste(dates, collapse = '_')}.nc"))
r <- rast(nc) |> 
  flip() |> 
  shift(dx = 0.0001) # problem with -180 extent

leaflet() |> 
  addProviderTiles(providers$Esri.OceanBasemap) |> 
  addRasterImage(r[[6]]) |> 
  addPolygons(data = p)

r_cefi <- rast("/Users/bbest/Downloads/myplot.3189115.1740009451.8932478.nc")
plot(r_cefi[[1]])
terra::plot()

# extent      : 0.5, 172.5, 0.7455357, 172.2545  (xmin, xmax, ymin, ymax)
r_cefi_r <- r_cefi |> 
  rotate()
r_cefi_rf <- r_cefi_r |> 
  flip()
# dimensions  : 337, 172, 4  (nrow, ncol, nlyr)
# resolution  : 1, 0.5089286  (x, y)
# extent      : -85.5, 86.5, 0.7455357, 172.2545  (xmin, xmax, ymin, ymax)
plet(
  r_cefi_rf[[1]],
  tiles = "Esri.OceanBasemap")

Code

librarian::shelf(
  dplyr, DT, ggplot2, glue, here, jsonlite, leaflet, mapview, ncdf4, readr, reticulate, 
  sf, stringr, terra, tidyr,
  quiet = T)

user        <- "bbest1"
pass_txt    <- "~/My Drive/private/data.marine.copernicus.eu_bbest1-password.txt"
dir_cm      <- here("data/copernicus")
results_csv <- here("data/copernicus.csv")

# do once: create virtual enviroment and install copernicusmarine Python module
# virtualenv_create(envname = "CopernicusMarine")
# virtualenv_install(envname = "CopernicusMarine", packages = c("copernicusmarine"))

# use virtualenv and reticulate::import copernicusmarine Python module
use_virtualenv(virtualenv = "CopernicusMarine", required = TRUE)
cmt <- import("copernicusmarine")

# login
pass <- readLines(pass_txt)
loggedin <- cmt$login(user, pass, skip_if_user_logged_in = T) # py_help(cmt$login)
# writes to /Users/bbest/.copernicusmarine/.copernicusmarine-credentials

datasets <- list(
  # Product identifier: GLOBAL_ANALYSISFORECAST_BGC_001_028
  # Product name: Global Ocean Biogeochemistry Analysis and Forecast
  "cmems_mod_glo_bgc-bio_anfc_0.25deg_P1M-m" = list(   # Primary production and O2, monthly
    vars      = list("nppv", "o2"),
    date_rng  = c("2024-01-01","2024-12-01"),
    depth_rng =  c(0.5, 5727.9) )) # get all depths

datasets |> 
  toJSON(auto_unbox = T, pretty=T)

Code

dir_cm <- here("data/productivity/copernicus")
dir.create(dir_cm, showWarnings = F, recursive = F)

for (i in 1:nrow(p_br)){ # i = 1
  
  p  <- slice(p_br, i)
  bb <- p |> 
    # st_buffer(9.25*1000) |> # buffer by a pixel width 1/12° ~ 9.25 km at equator 
    # but -180 goes global width
    st_bbox()
  
  message(glue("{i}/{nrow(p_br)}: {p$boemrgn_key}"))
  
  for (j in 1:length(datasets)){ # j=1
    
    ds_id <- names(datasets)[j]
    ds <- datasets[[j]]
    
    nc <- here(glue(
      "{dir_cm}/{ds_id}_{paste(ds$date_rng, collapse = '_')}_{p$boemrgn_key}{p$chunk}.nc"))
    
    if (file.exists(nc)){
      message(glue("{j}/{length(datasets)} skip: {basename(nc)} already exists"))
      next()
    }
    
    message(glue("  {j}/{length(datasets)}: fetch {basename(nc)}  ~ {Sys.time()}"))
    
    # py_help(cmt$subset)
    out <- cmt$subset(
      dataset_id            = ds_id,
      # dataset_version     = "202309",
      variables             = ds$vars,
      minimum_longitude     = bb$xmin,
      maximum_longitude     = bb$xmax,
      minimum_latitude      = bb$ymin,
      maximum_latitude      = bb$ymax,
      start_datetime        = paste0(ds$date_rng[1], "T00:00:00"),
      end_datetime          = paste0(ds$date_rng[2], "T00:00:00"),
      minimum_depth         = ds$depth_rng[1],
      maximum_depth         = ds$depth_rng[2],
      output_filename       = nc,
      force_download        = T)
  } }
message(glue("{nrow(p_br)}/{nrow(p_br)} done: {basename(nc)}  ~ {Sys.time()}"))

Code

o <- nc_open(nc)
o

r <- rast(nc) |>
  crop(bb_e, snap = "in") # problem with -180 extent
r

i <- 1
cat(glue("
  index: {i}
  name: {names(r)[i]}
  time: {time(r[[i]])}
  varname: {varnames(r[[i]])}
  longname: {longnames(r[[i]])}"))

# Plot the raster with the beomrgn
plet(
  r[[i]], 
  main  = glue("{longnames(r[[i]])}\n{names(r)[i]}\n{time(r[[i]])}"),
  # main  = glue("{longnames(r[[i]])}\nthetao_depth=0.5\n{time(r[[i]])}"), 
  tiles = "Esri.OceanBasemap") |> 
  addPolygons(
    data        = p,
    fillOpacity = 0.2)

Code

for (i in 1:nrow(p_br)){ # i = 1
  
  p  <- slice(p_br, i)
   
  j = 1
  ds_id <- names(datasets)[j]
  ds <- datasets[[j]]
  
  nc <- here(glue(
    "{dir_cm}/{ds_id}_{paste(ds$date_rng, collapse = '_')}_{p$boemrgn_key}{p$chunk}.nc"))
  
  r <- rast(nc)
  
  times <- time(r["nppv_"]) |> as.character()
  for (m in 1:12){ # m = 1
    
    tif <- glue("{dir_cm}/nppv_{p$boemrgn_key}{p$chunk}_{sprintf('2024-%02d', m)}.tif")
    if (file.exists(tif))
      next()
    
    r["nppv_"] |> 
      subset(times == sprintf("2024-%02d-01", m)) |> 
      sum(na.rm = T) |> 
      writeRaster(tif)
  }
  
  tif <- glue("{dir_cm}/nppv_{p$boemrgn_key}{p$chunk}_2024_monthly-avg.tif")
  if (file.exists(tif))
    next()
  
  rast(list.files(dir_cm, glue("nppv_{p$boemrgn_key}{p$chunk}_2024-.*\\.tif$"), full.names = T)) |> 
    mean(na.rm = T) |> 
    writeRaster(tif)
}


r <- rast(tif)  # |>
  # crop(bb_e, snap = "in") # problem with -180 extent

plet(
  r, 
  main  = basename(tif),
  tiles = "Esri.OceanBasemap") |> 
  addPolygons(
    data        = p,
    fillOpacity = 0.2)

--- title: "Ingest Primary Productivity" editor: visual editor_options: chunk_output_type: console --- ```{r} #| label: setup librarian::shelf( cmocean, DT, fs, glue, here, htmltools, janitor, jsonlite, leaflet, leaflet.extras, mapview, rerddap, R.utils, sf, stringr, terra, units, quiet = T) source(here("libs/db.R")) # define: database connection (con) yr <- 2023 # vgpm-viirs most recent complete year url_tar <- glue("http://orca.science.oregonstate.edu/data/1x2/monthly/vgpm.r2022.v.chl.v.sst/hdf/vgpm.v.{yr}.tar") file_tar <- glue("/Users/bbest/My Drive/projects/msens/data/raw/oregonstate.edu/vgpm.v.{yr}.tar") dir_out <- "/Users/bbest/My Drive/projects/msens/data/derived" dir_raw <- "/Users/bbest/My Drive/projects/msens/data/raw/boem.gov" vgpm_tif <- glue("{dir_out}/vgpm-viirs_{yr}-daily-avg.tif") oa_new_geo <- glue("{dir_raw}/OCS_Area_Polygons.geojson") pa_new_geo <- glue("{dir_raw}/OCS_Proposed_Planning_Area_Polygons.geojson") pa_oa_geo <- glue("{dir_out}/pa_oa_2025.geojson") verbose <- F ``` ## Vertically Generalized Production Model (VGPM), VIIRS `r yr` | OregonState Source: - [Ocean Productivity from Oregon State](http://sites.science.oregonstate.edu/ocean.productivity) - [Standara VGPM: 1080 by 2160 Monthly HDF files from VIIRS R2022](http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.v.chl.v.sst.php) - [Frequently Asked Questions](http://orca.science.oregonstate.edu/faq01.php) Calculated sum of `r yr` monthly net primary production (NPP) from VIIRS data. > Community guidance for developing this website was to provide a single productivity product as a Standard product. For this, we initially chose the original Vertically Generalized Production Model (VGPM) (Behrenfeld and Falkowski, 1997a) as the standard algorithm. The VGPM is a "chlorophyll-based" model that estimate net primary production from chlorophyll using a temperature-dependent description of chlorophyll-specific photosynthetic efficiency. For the VGPM, net primary production is a function of chlorophyll, available light, and the photosynthetic efficiency. > > Standard products are based on chlorophyll, temperature and PAR data from SeaWiFS, MODIS and VIIRS satellites, along with estimates of euphotic zone depth from a model developed by Morel and Berthon (1989) and based on chlorophyll concentration. Units: mg C / m^2^ / day Citation: [Behrenfeld and Falkowski, 1997a](http://science.oregonstate.edu/ocean.productivity/references.php#Behrenfeld.1997a) ```{r} #| label: get_vgpm if (!file.exists(vgpm_tif)) { if(!file.exists(file_tar)) download.file(url_tar, file_tar) dir_untar <- path_ext_remove(file_tar) dir.create(dir_untar, showWarnings = F) untar(file_tar, exdir = dir_untar) # Get list of compressed files files_gz <- dir_ls(dir_untar, glob = "*.gz") dir_unzip <- file.path(dir_untar, "unzipped") dir.create(dir_unzip, showWarnings = F) # Unzip all files for (file_gz in files_gz) { # file_gz = files_gz[1] file_hdf <- file.path(dir_unzip, sub("\\.gz$", "", basename(file_gz))) if (!file.exists(file_hdf)) gunzip(file_gz, destname = file_hdf) } # Get list of unzipped HDF files files_hdf <- dir_ls(dir_unzip, glob = "*.hdf") # length(files_hdf) = 12 # Function to read VGPM data from HDF file read_vgpm_hdf <- function(file_hdf) { # file_hdf <- files_hdf[1] # First, query the HDF file to see the subdatasets sds_info <- terra::describe(file_hdf) if (verbose) print(sds_info) tryCatch({ r <- terra::rast(file_hdf, noflip=T) # Standard VGPM product is typically 1/6 degree global grid (-180 to 180, -90 to 90) ext(r) <- c(-180, 180, -90, 90) crs(r) <- "EPSG:4326" r[r==-9999] <- NA # plot(r) return(r) }, error = function(e) { warning("Error reading", file_hdf, ":", e$message, "\n") return(NULL) }) } # Initialize a list to store all monthly rasters monthly_rasters <- list() # Read all monthly files for (i in seq_along(files_hdf)) { if (verbose) message("Processing file", i, "of", length(files_hdf), ":", files_hdf[i], "\n") r <- read_vgpm_hdf(files_hdf[i]) if (!is.null(r)) { monthly_rasters[[i]] <- r } } # Create a SpatRaster collection from all the monthly rasters s_vgpm <- rast(monthly_rasters) r <- app(s_vgpm, mean, na.rm = TRUE) names(r) <- "npp_daily_avg" writeRaster(r, vgpm_tif, overwrite = TRUE) message(glue( "Annual products created and saved. - daily avg for {yr}: {vgpm_tif}")) } r <- rast(vgpm_tif) # |> # rotate(left=F) # -180 to 180 -> 0 to 360 # shift(dx = 0.0001) # problem with -180 extent mapView( r, maxpixels = 2332800, col.regions = cmocean("algae"), layer.name = "raster: daily avg NPP") |> slot("map") |> addFullscreenControl() |> addControl(tags$div(HTML("units: mg C / m<sup>2</sup> / day")), position = "topright") ``` ## NEW Proposed Planning Areas BOEM internal: - [Marine Ecoregions](https://boem.maps.arcgis.com/home/item.html?id=b00a15d7e4224ef98e496542d9a9d25f) - [OCS Areas](https://boem.maps.arcgis.com/home/item.html?id=8dbc5e97fc3346f58466c2de7b2fc870) - [Proposed Planning Areas](https://boem.maps.arcgis.com/home/item.html?id=b833761e456f4cbfb227a4b5a9632b89) ```{r} #| label: get_pa set_geom_ctr_area <- function(x){ x |> st_set_geometry("geom") |> mutate( ctr = st_centroid(geom), ctr_lon = ctr |> st_coordinates() %>% .[,"X"], ctr_lat = ctr |> st_coordinates() %>% .[,"Y"], area_km2 = st_area(geom) |> set_units(km^2) |> as.numeric() ) |> select(-ctr) } drop_ctr_area <- function(x){ x |> select(-ctr_lon, -ctr_lat, -area_km2) } if (!file.exists(pa_oa_geo)){ pa0 <- st_read(con, "ply_planareas_s05") |> st_drop_geometry() |> select(planarea_key, planarea_name) |> mutate( planarea_name = str_replace(planarea_name, "Gulf of Mexico", "Gulf of America"), planarea_key = str_replace(planarea_key, "pa", ""), planarea_key = case_match( planarea_key, "CGM" ~ "CGA", "EGM" ~ "EGA", "WGM" ~ "WGA", .default = planarea_key)) pa <- read_sf(pa_new_geo) |> clean_names() |> rename( region_name = region, planarea_key = planning_area) |> mutate( region_name = str_replace(region_name, "Gulf of Mexico", "Gulf of America"), region_key = case_match( region_name, "Alaska" ~ "AK", "Atlantic" ~ "AT", "Gulf of America" ~ "GA", "Pacific" ~ "PA"), planarea_key = case_match( planarea_key, "CG" ~ "CGA", "EG" ~ "EGA", "WG" ~ "WGA", .default = planarea_key)) |> left_join( pa0, by = "planarea_key") |> mutate( planarea_name = case_match( planarea_key, "HAR" ~ "High Arctic", .default = planarea_name)) # TODO: planarea_name for HAR is ...? oa <- read_sf(oa_new_geo) |> clean_names() |> rename( region_name = region, planarea_name = area_name) |> filter(!planarea_name %in% c("Alaska", "Pacific Coast", "Atlantic Coast", "Gulf of Mexico")) |> mutate( region_key = case_match( region_name, "Atlantic" ~ "AT", "Pacific" ~ "PA"), planarea_key = case_match( planarea_name, "American Samoa" ~ "AMS", "Guam and Commonwealth of the Northern Mariana Islands" ~ "GUA", "Hawaiian Islands and Midway Atoll" ~ "HAW", "Howland and Baker Islands" ~ "HOW", "Jarvis Island" ~ "JAR", "Johnston Atoll" ~ "JOH", "Palmyra Atoll and Kingman Reef" ~ "PAL", "Puerto Rico and U.S. Virgin Islands" ~ "PUR", "Wake Island" ~ "WAK")) pa_oa <- bind_rows( pa, oa) |> select(region_key, region_name, planarea_key, planarea_name) |> st_make_valid() |> set_geom_ctr_area() |> st_shift_longitude() |> arrange(region_key, planarea_key) write_sf(pa_oa, pa_oa_geo) } pa_oa <- read_sf(pa_oa_geo) mapView(pa_oa) pa_oa |> st_drop_geometry() |> datatable() ``` ## VGPM in Proposed Planning Areas Extract daily average net primary production (NPP) values across all cells within the proposed planning areas. ```{r} #| label: get_pa_vgpm r <- rast(vgpm_tif) r_r <- r |> rotate(left=F) # -180 to 180 -> 0 to 360 # shift(dx = 0.0001) # problem with -180 extent pa_oa_npp <- pa_oa |> mutate( npp_cellavg = zonal(r_r, pa_oa |> vect(), fun = "mean", na.rm = T) |> pull(npp_daily_avg)) (mapView( r, maxpixels = 2332800, col.regions = cmocean("algae"), layer.name = "raster") + mapView( pa_oa_npp, zcol = "npp_cellavg", col.regions = cmocean("algae"), layer.name = "planning areas")) |> slot("map") |> addFullscreenControl() |> addControl(tags$div(HTML("units: mg C / m<sup>2</sup> / day")), position = "topright") pa_oa_npp |> st_drop_geometry() |> datatable() ``` ```{r} #| label: old_br_east_west #| eval: false # Latest # # - [ERDDAP - search "productivity_viirs_snpp"](https://coastwatch.pfeg.noaa.gov/erddap/search/index.html?page=1&itemsPerPage=1000&searchFor=productivity_viirs_snpp) # # - monthly has lots of cloud cover -- see [WMS map](https://coastwatch.pfeg.noaa.gov/erddap/wms/productivity_viirs_snpp_monthly/index.html) # # - [Arctic Ocean Primary Productivity: The Response of Marine Algae to Climate Warming and Sea Ice Decline - NOAA Arctic](https://arctic.noaa.gov/report-card/report-card-2023/arctic-ocean-primary-productivity-the-response-of-marine-algae-to-climate-warming-and-sea-ice-decline-2023/#methods_data) ![](https://arctic.noaa.gov/wp-content/uploads/2023/11/ARC23-PrimaryProductivity-frey-Fig3-878x1024.png) dir_nc <- here("data/productivity") dates <- c("2023-01-16", "2023-12-16") # https://tile.marinesensitivity.org/public.ply_boemrgns.html p_br <- read_sf(con, "ply_boemrgns") # chunk east and west of antimeridian for extraction ---- bb_w <- st_bbox(c(xmin = 0, xmax = 180, ymin = -90, ymax = 90)) |> st_as_sfc() |> st_as_sf(crs = 4326) bb_e <- st_bbox(c(xmin = -180, xmax = 0, ymin = -90, ymax = 90)) |> st_as_sfc() |> st_as_sf(crs = 4326) p_br <- bind_rows( p_br |> st_intersection(bb_w) |> mutate( chunk = "w"), p_br |> st_difference(bb_w) |> mutate( chunk = "e")) |> relocate(chunk, .after = boemrgn_key) |> arrange(boemrgn_key, chunk) leaflet() |> addProviderTiles(providers$Esri.OceanBasemap) |> addPolygons(data = st_shift_longitude(p_br)) p_br |> st_drop_geometry() |> datatable() ``` ```{r} #| label: old_erddap_fetch #| eval: false # iterate over polygon chunks ---- for (i in 1:nrow(p_br)){ # p = p_br[3,] p <- slice(p_br, i) bb <- st_bbox(p) nc <- here(glue( "data/productivity/{p$boemrgn_key}{p$chunk}_{paste(dates, collapse = '_')}.nc")) if (file.exists(nc)){ message(glue("{i}/{nrow(p_br)} skip: {basename(nc)} already exists")) next() } message(glue("{i}/{nrow(p_br)} fetch: {basename(nc)} ~ {Sys.time()}")) o <- griddap( "productivity_viirs_snpp_monthly", url = "https://coastwatch.pfeg.noaa.gov/erddap/", time = dates, longitude = c(max(bb$xmin, -179.9792), min(bb$xmax, 179.9792)), # TODO: compare with actual min/max from fetching latitude = c(bb$ymin, bb$ymax), fields = "productivity", store = disk(path = dir_nc)) file.rename(o$summary$filename, nc) } message(glue("{nrow(p_br)}/{nrow(p_br)} done: {basename(nc)} ~ {Sys.time()}")) boemrgn_key = "brAK" chunk = "e" p <- p_br |> filter( boemrgn_key == !!boemrgn_key, chunk == !!chunk) nc <- here(glue("data/productivity/{boemrgn_key}{chunk}_{paste(dates, collapse = '_')}.nc")) r <- rast(nc) |> flip() |> shift(dx = 0.0001) # problem with -180 extent leaflet() |> addProviderTiles(providers$Esri.OceanBasemap) |> addRasterImage(r[[6]]) |> addPolygons(data = p) r_cefi <- rast("/Users/bbest/Downloads/myplot.3189115.1740009451.8932478.nc") plot(r_cefi[[1]]) terra::plot() # extent : 0.5, 172.5, 0.7455357, 172.2545 (xmin, xmax, ymin, ymax) r_cefi_r <- r_cefi |> rotate() r_cefi_rf <- r_cefi_r |> flip() # dimensions : 337, 172, 4 (nrow, ncol, nlyr) # resolution : 1, 0.5089286 (x, y) # extent : -85.5, 86.5, 0.7455357, 172.2545 (xmin, xmax, ymin, ymax) plet( r_cefi_rf[[1]], tiles = "Esri.OceanBasemap") ```  ```{r} #| label: old_copernicus_setup #| eval: false librarian::shelf( dplyr, DT, ggplot2, glue, here, jsonlite, leaflet, mapview, ncdf4, readr, reticulate, sf, stringr, terra, tidyr, quiet = T) user <- "bbest1" pass_txt <- "~/My Drive/private/data.marine.copernicus.eu_bbest1-password.txt" dir_cm <- here("data/copernicus") results_csv <- here("data/copernicus.csv") # do once: create virtual enviroment and install copernicusmarine Python module # virtualenv_create(envname = "CopernicusMarine") # virtualenv_install(envname = "CopernicusMarine", packages = c("copernicusmarine")) # use virtualenv and reticulate::import copernicusmarine Python module use_virtualenv(virtualenv = "CopernicusMarine", required = TRUE) cmt <- import("copernicusmarine") # login pass <- readLines(pass_txt) loggedin <- cmt$login(user, pass, skip_if_user_logged_in = T) # py_help(cmt$login) # writes to /Users/bbest/.copernicusmarine/.copernicusmarine-credentials datasets <- list( # Product identifier: GLOBAL_ANALYSISFORECAST_BGC_001_028 # Product name: Global Ocean Biogeochemistry Analysis and Forecast "cmems_mod_glo_bgc-bio_anfc_0.25deg_P1M-m" = list( # Primary production and O2, monthly vars = list("nppv", "o2"), date_rng = c("2024-01-01","2024-12-01"), depth_rng = c(0.5, 5727.9) )) # get all depths datasets |> toJSON(auto_unbox = T, pretty=T) ``` ```{r} #| label: old_copernicus_fetch #| eval: false dir_cm <- here("data/productivity/copernicus") dir.create(dir_cm, showWarnings = F, recursive = F) for (i in 1:nrow(p_br)){ # i = 1 p <- slice(p_br, i) bb <- p |> # st_buffer(9.25*1000) |> # buffer by a pixel width 1/12° ~ 9.25 km at equator # but -180 goes global width st_bbox() message(glue("{i}/{nrow(p_br)}: {p$boemrgn_key}")) for (j in 1:length(datasets)){ # j=1 ds_id <- names(datasets)[j] ds <- datasets[[j]] nc <- here(glue( "{dir_cm}/{ds_id}_{paste(ds$date_rng, collapse = '_')}_{p$boemrgn_key}{p$chunk}.nc")) if (file.exists(nc)){ message(glue("{j}/{length(datasets)} skip: {basename(nc)} already exists")) next() } message(glue(" {j}/{length(datasets)}: fetch {basename(nc)} ~ {Sys.time()}")) # py_help(cmt$subset) out <- cmt$subset( dataset_id = ds_id, # dataset_version = "202309", variables = ds$vars, minimum_longitude = bb$xmin, maximum_longitude = bb$xmax, minimum_latitude = bb$ymin, maximum_latitude = bb$ymax, start_datetime = paste0(ds$date_rng[1], "T00:00:00"), end_datetime = paste0(ds$date_rng[2], "T00:00:00"), minimum_depth = ds$depth_rng[1], maximum_depth = ds$depth_rng[2], output_filename = nc, force_download = T) } } message(glue("{nrow(p_br)}/{nrow(p_br)} done: {basename(nc)} ~ {Sys.time()}")) ``` ```{r} #| label: old_copernicus_nc #| eval: false o <- nc_open(nc) o r <- rast(nc) |> crop(bb_e, snap = "in") # problem with -180 extent r i <- 1 cat(glue(" index: {i} name: {names(r)[i]} time: {time(r[[i]])} varname: {varnames(r[[i]])} longname: {longnames(r[[i]])}")) # Plot the raster with the beomrgn plet( r[[i]], main = glue("{longnames(r[[i]])}\n{names(r)[i]}\n{time(r[[i]])}"), # main = glue("{longnames(r[[i]])}\nthetao_depth=0.5\n{time(r[[i]])}"), tiles = "Esri.OceanBasemap") |> addPolygons( data = p, fillOpacity = 0.2) ``` ```{r} #| label: old_copernicus_summarize #| eval: false for (i in 1:nrow(p_br)){ # i = 1 p <- slice(p_br, i) j = 1 ds_id <- names(datasets)[j] ds <- datasets[[j]] nc <- here(glue( "{dir_cm}/{ds_id}_{paste(ds$date_rng, collapse = '_')}_{p$boemrgn_key}{p$chunk}.nc")) r <- rast(nc) times <- time(r["nppv_"]) |> as.character() for (m in 1:12){ # m = 1 tif <- glue("{dir_cm}/nppv_{p$boemrgn_key}{p$chunk}_{sprintf('2024-%02d', m)}.tif") if (file.exists(tif)) next() r["nppv_"] |> subset(times == sprintf("2024-%02d-01", m)) |> sum(na.rm = T) |> writeRaster(tif) } tif <- glue("{dir_cm}/nppv_{p$boemrgn_key}{p$chunk}_2024_monthly-avg.tif") if (file.exists(tif)) next() rast(list.files(dir_cm, glue("nppv_{p$boemrgn_key}{p$chunk}_2024-.*\\.tif$"), full.names = T)) |> mean(na.rm = T) |> writeRaster(tif) } r <- rast(tif) # |> # crop(bb_e, snap = "in") # problem with -180 extent plet( r, main = basename(tif), tiles = "Esri.OceanBasemap") |> addPolygons( data = p, fillOpacity = 0.2) ```