In afsc-gap-products/coldpool: AFSC/RACE Groundfish Assessment Program EBS and NBS temperature products

1.0 Introduction

This document provides code and documentation related to regional model-based versions of the cold pool indices and interpolated bottom temperature rasters. These model-based indices are used in a variety of research and management contexts to complement the survey-based indices. This example demonstrates how to use the utilities in this package to develop these products.

This script was developed with the ROMS Bering10K hindcast in mind; see Kearney et al. (2020) for more information regarding this specific model. However, the code presented here is model agnostic, and can in theory be applied to any output that covers the eastern Bering Sea shelf region.

The calculation of these survey-based cold pool metrics is part of the annual hindcast update procedure that is run every summer. This script is dependent on two files that are generated as part of that update process prior to running:

• survey-replicated data file, .csv format. See "Survey-replicated data", below.
• cold pool indices file, .nc format. See "Annual indices", below.

For details regarding how those files are built, see the beringnpz/Bering10KPostprocessing repo (specifically, the surveyreplicatedbtemp.m, buildcoldpoolncfile.m, and addjul1btemp.m functions). The current coldpool-based code should be run after that process is complete in order to add survey and survey-replicated indices to the cold pool indices file.

While this script is written similar to a vignette, we note that the current paths to these files will not be accessible to most users. Therefore, the script is meant to be used by the regional modeling team to run periodic updates, and otherwise is provided here as documentation only.

1.1 Setup

The following packages are used in this script, and should be loaded prior to running:

• coldpool
• akfgmaps
• lubridate
• ncdf4

We start with some general setup that loads these packages and sets output figure resolution and the map projection to use for our calculations.

library(coldpool)

develflag = TRUE
if (develflag) {
  devtools::load_all()
} else {
  library(coldpool)
}

library(akgfmaps)
library(lubridate)
library(ncdf4)

# Global options, mirroring primary 0_update_cold_pool_index.Rmd
fig_res <- 600
proj_crs <- coldpool:::ebs_proj_crs

All model-related files (both input and output) will be placed in a folder outside this repo, alongside other model output. Here, we choose the Level 3 folder associated with the operational hindcast simulation of the Bering10K ROMS model (see the Bering10K Dataset Docs for more information). These details should be modified as needed for other users, computing locations, and/or models.

# Path to local mount of UW HPC folder mox.hyak.uw.edu:/gscratch/bumblereem/
moxdir = "~/Documents/mox_bumblereem" # Note: computer-specific, change as needed

# Specific simulation
simname = "B10K-K20_CORECFS" # weekly-archived "operational" hindcast
#simname = "B10K-K20nobio_CORECFS_daily" # daily-archived variant

# Model's Level 3 folder path
lev3path <- file.path(moxdir, "roms_for_public", simname, "Level3")

1.2 Survey-replicated data

We can perform a rough comparison of model output and the survey data by simply extracting the model bottom temperature field on a specific date of each year, for example, a mid-summer date like July 1. However, the groundfish survey temperature data are collected over a 2-3 month period rather than on a single day, so the rasters built by this package likewise reflect that non-synoptic nature, with regions in the southeast section typically more representative of early summer and the northern regions reflecting late summer conditions. To provide a true one:one comparison between model and observations, we need to build model-derived rasters that reflect a similar pattern.

The first step of this process is to build a survey-replicated dataset. We do this by sampling the model at the closest grid cell and archived output time step to each point in the survey dataset. The data is saved in a file that is identical to the [colpool]/data/index_hauls_temperature_data.csv file but with an additional column added for model bottom temperature. The file is generated external to this script.

srepcsv_path <- file.path(lev3path, sprintf("survey_replicates_%s.csv", simname))

1.3 Masks

The primary cold pool index calculations from this package rely on predefined masks from the akgfmaps R package. For our model-based indices, we use the southeastern Bering Sea (SEBS) mask from that package as well as a few variants: - The SEBS region without the northwest strata (i.e. strata 82 and 90). These strata were added to the sampling a few years later than most of the EBS, and for some projects we ignore that region to maintain data continuity. - The SEBS region (with and without the NW parts) with the shelf break set based on model bathymetry rather than real-world bathymetry. The shelf slope in the Bering10K ROMS model is smoothed in order to avoid numerical errors in the horizontal pressure gradient that are characteristic of sigma-coordinate models like ROMS in areas of steep topography. This results in the model having a slightly narrower shelf than the real world. This variant allows for better comparison between model and observations because it eliminates the area where we expect the two to disagree due to mismatched bottom depth.

# The default SEBS mask

akSEBS <- akgfmaps::get_base_layers(select.region = "bs.south", set.crs = proj_crs)

# Read polygon defining the Bering10K shelf (i.e. <200m depth)

modelshelfpolyfile = here::here("inst", "extdata", "B10K_lte200m_polygon.shp")

b10k_lte_200 <- sf::st_read(modelshelfpolyfile, quiet = TRUE) |>
  st_set_crs("+proj=longlat") |>
  sf::st_transform(crs = sf::st_crs(proj_crs)) 

# SEBS without northwest strata (i.e. strata ID <= 62)

akSEBSnonw <- akSEBS$survey.strata |>
  dplyr::filter(Stratum <= 62) |>
  dplyr::group_by(SURVEY) |>
  dplyr::summarise()

# Build 4 masks: SEBS w/o northwest, 
#                SEBS
#                SEBS w/o northwest and trimmed to model shelf
#                SEBS trimmed to model shelf
# (Note: order reflects that in the cold pool index .nc file)

mymask <- list(akSEBSnonw, 
               "sebs",
               st_intersection(akSEBSnonw, b10k_lte_200), 
               st_intersection(akSEBS$survey.area, b10k_lte_200))

maskname <- c("SEBS-noNW", "SEBS", "SEBS-noNW-modelshelf", "SEBS-modelshelf")

1.4 Cold pool index netCDF file

Model-based indices will be saved to a prexisting netCDF file.

cpoolnc_path <- file.path(lev3path, sprintf("%s_coldpool.nc", simname))

2.0 Generate bottom temperature rasters

With our datasets in place, we first generate raster images with bottom temperature (from both survey and survey-replicated model) interpolated to 5-km grid. These rasters will form the base for many of the following calculations. Note that the survey-derived rasters should be identical to those created as part of the annual coldpool updates.

# Years: Start with all years in file

temperature_df <- read.csv(file = srepcsv_path,
                           stringsAsFactors = FALSE)

names(temperature_df) <- tolower(names(temperature_df))
year_vec <- sort(unique(temperature_df$year))

# Interpolation parameters

cell_resolution = 5000 # m, i.e. 5km
methods = "Ste" # Stein's Matern kriging

# Create one interpolated raster per year per mask per data source

outpath = file.path(moxdir, "roms_for_public", simname, "Level3", "coldpool_tifs");

for (ii in 1:length(year_vec)) {

  for (im in 1:length(mymask)) {

    # File-naming convention

    outfilefun = function(mask, year, method, variable) 
      file.path(outpath, 
                sprintf("%s_%s_%s_%d.tif", maskname[im], method, variable, year))  

    print(sprintf("Creating rasters: mask #%d, %d", im, year_vec[ii]))

    # Create rasters

    if (!file.exists(outfilefun(NULL, year_vec[ii], methods, "gear_temperature"))) {

      # ... for survey data

      interpolate_variable(dat = dplyr::filter(temperature_df, year == year_vec[ii]),
                                   dat.year = year_vec[ii],
                                   in.crs = "EPSG:4326",
                                   interpolation.crs = proj_crs,
                                   cell.resolution = cell_resolution,
                                   lon.col = "longitude",
                                   lat.col = "latitude",
                                   var.col = "gear_temperature",
                                   nm = Inf,
                                   outputfilefun = outfilefun,
                                   select.region = mymask[[im]],
                                   methods = methods)
    }

    if (!file.exists(outfilefun(NULL, year_vec[ii], methods, "model_bottom_temp"))) {

      # ... for model data

      interpolate_variable(dat = dplyr::filter(temperature_df, year == year_vec[ii]),
                                   dat.year = year_vec[ii],
                                   in.crs = "EPSG:4326",
                                   interpolation.crs = proj_crs,
                                   cell.resolution = cell_resolution,
                                   lon.col = "longitude",
                                   lat.col = "latitude",
                                   var.col = "model_bottom_temp",
                                   nm = Inf,
                                   outputfilefun = outfilefun,
                                   select.region = mymask[[im]],
                                   methods = methods)
    }
  }
}

3.0 Annual indices

Metrics of mean bottom temperature and cold pool fractions are calculated for both the survey-based and model-based rasters. The resulting survey-based time series using the SEBS mask should be identical to the one included within the coldpool package.

# Open netCDF file

nc <- nc_open(cpoolnc_path, write=TRUE)

# Read relevant coordinate data

thresh    <- ncvar_get(nc, "threshold")
tdaysince <- ncvar_get(nc, "time")
tyear <- year(tdaysince + as.Date("1900-01-01"))

# Read existing data

btfile <- ncvar_get(nc, "average_bottom_temp", start=c(1,2,1))
cpfile <- ncvar_get(nc, "cold_pool_index", start=c(1,1,2,1))

# Loop over new rasters and add data to file as needed

for (ii in 1:length(year_vec)) {

  yidx = which(year_vec[ii] == tyear)

  filemissingdata <- all(is.na(btfile[,,yidx])) | all(is.na(cpfile[,,,yidx]))

  if (!filemissingdata & (any(is.na(btfile[,,yidx])) | any(is.na(cpfile[,,,yidx])))) {
    warning(sprintf("Data for year %d includes some missing; check this!", year_vec[ii]))
  }

  if (filemissingdata) {

    for (im in 1:length(mymask)) {

      print(sprintf("Saving data to file: %d mask %d:", year_vec[ii], im));

      # Raster files for this year/mask

      svyraster = terra::rast(file.path(outpath, 
                                      sprintf("%s_%s_%s_%d.tif", 
                                              maskname[im], 
                                              "Ste", 
                                              "gear_temperature",  
                                              year_vec[ii])))  
      modraster = terra::rast(file.path(outpath, 
                                      sprintf("%s_%s_%s_%d.tif", 
                                              maskname[im], 
                                              "Ste", 
                                              "model_bottom_temp", 
                                              year_vec[ii]))) 

      # Calculate average bottom temperature and save to file

      btsvy <- mean(terra::values(svyraster), na.rm = TRUE)
      btmod <- mean(terra::values(modraster), na.rm = TRUE)

      ncvar_put(nc, "average_bottom_temp", btsvy, start=c(im,2,yidx), count=c(1,1,1))
      ncvar_put(nc, "average_bottom_temp", btmod, start=c(im,3,yidx), count=c(1,1,1))

      # Calculate and save cold pool indices

      totarea <- svyraster |>
        cpa_from_raster(raster_units = "m", temperature_threshold = 100)

      for (ith in 1:length(thresh)) {

        cpsvy <- (svyraster |> 
                    cpa_from_raster(raster_units = "m", temperature_threshold = thresh[ith]))/totarea
        cpmod <- (modraster |> 
                    cpa_from_raster(raster_units = "m", temperature_threshold = thresh[ith]))/totarea

        ncvar_put(nc, "cold_pool_index", cpsvy, start=c(ith,im,2,yidx), count=c(1,1,1,1))
        ncvar_put(nc, "cold_pool_index", cpmod, start=c(ith,im,3,yidx), count=c(1,1,1,1))

      }
    }

    # Update history attribute

    hisold = ncatt_get(nc, 0, "history")
    hisnew = paste(format(Sys.time(), "%a %b %d %H:%M:%S %Y:"),
                   "Survey and survey-replicated data added for year", 
                   year_vec[ii])
    ncatt_put(nc, 0, "history", paste(hisnew, hisold$value, sep="\n"))

  }
}

afsc-gap-products/coldpool documentation built on Sept. 14, 2024, 7:40 p.m.