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:
+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.
+The following packages are used in this script, and should be loaded +prior to running:
+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)library(akgfmaps)
+library(lubridate)
+library(ncdf4)
+
+# Global options, mirroring primary 0_update_cold_pool_index.Rmd
+fig_res <- 600
+proj_crs <- coldpool:::ebs_proj_crsAll 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")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))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")Model-based indices will be saved to a prexisting netCDF file.
+cpoolnc_path <- file.path(lev3path, sprintf("%s_coldpool.nc", simname))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)
+ }
+ }
+}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)
+
+# 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])
+ ncattput(nc, 0, "history", paste(hisnew, hisold$value, sep="\n"))
+
+ }
+}