From fc31ee30a006dbdc1b01977dcdfbf405e1ad55bb Mon Sep 17 00:00:00 2001 From: "Kelly.Kearney" Date: Thu, 8 Feb 2024 11:49:59 -0800 Subject: [PATCH 1/6] Updated interpolate_variable to support more flexible masking - Added outputfilefun input option to custimize raster output path files - Added support for sf object region masks Both of these options allow for indices related to the ROMS model-based cold pool indices. Changes should not conflict with (and all defaults support) the primary survey data workflow. --- R/interpolate_variable.R | 133 +++++++++++++++++++++++++-------------- 1 file changed, 87 insertions(+), 46 deletions(-) diff --git a/R/interpolate_variable.R b/R/interpolate_variable.R index e7a4723..ab8a8d0 100644 --- a/R/interpolate_variable.R +++ b/R/interpolate_variable.R @@ -19,6 +19,7 @@ #' specific bounds, or a 4-element named vector with values xmin, xmax, ymin, #' and ymax. #' @param return.raster Should a raster be returned? +#' @param outputfilefun Function of the form f(mask, year, method, var) that returns a string holding the full path name where the output raster should be stored #' @return Returns a data frame containing cold pool areas estimated by interpolation methods, for cutoffs at zero degrees and two degrees C. If argument write.to.file = TRUE, also writes a GeoTIFF raster to output directory. #' @export @@ -36,7 +37,12 @@ interpolate_variable <- function(dat, select.region = "sebs", methods = c("NN", "IDW", "IDW4", "Exp", "Sph", "Bes", "Gau", "Cir", "Mat", "Ste", "Tps"), bbox = NULL, - return_raster = FALSE) { + return_raster = FALSE, + outputfilefun = NULL) { + + #--------------------- + # Input checks + #--------------------- if(!all(tolower(methods) %in% c("nn", "idw", "idw4", "exp", "sph", "bes", "gau", "cir", "mat", "ste", "tps" ))) { @@ -64,11 +70,12 @@ interpolate_variable <- function(dat, } + #--------------------- + # Raster setup + #--------------------- + # When select.region is a character vector if(is.character(select.region)) { - - fpath_subdir <- select.region - # Load survey extent for masking ---- region_mask <- akgfmaps::get_base_layers(select.region = select.region, set.crs = interpolation.crs)$survey.area @@ -76,59 +83,66 @@ interpolate_variable <- function(dat, # When select.region is an sf object if(is.data.frame(select.region)) { - - fpath_subdir <- ifelse(is.na(pre), "var", pre) - region_mask <- select.region - - } - - if(!dir.exists(here::here("output", "raster", fpath_subdir, var.col))) { - dir.create(here::here("output", "raster", fpath_subdir, var.col), recursive = TRUE) } # Set dimensions for raster cells ---- if(is.null(bbox)) { + if (is.character(select.region)) { # Default interpolation bounds for named regions - if(select.region %in% c("sebs", "bs.south")) { - # Make raster for interpolation ---- - n_dim <- floor(abs(-1625000 - -35000))/cell.resolution - - interp_raster <- terra::rast(xmin = -1625000, - xmax = -35000, - ymin = 379500, - ymax = 1969500, - nrow = n_dim, - ncol = n_dim, - crs = interpolation.crs) + if(select.region %in% c("sebs", "bs.south")) { + # Make raster for interpolation ---- + n_dim <- floor(abs(-1625000 - -35000))/cell.resolution - } else { + interp_raster <- terra::rast(xmin = -1625000, + xmax = -35000, + ymin = 379500, + ymax = 1969500, + nrow = n_dim, + ncol = n_dim, + crs = interpolation.crs) - if(select.region %in% c("bs.north", "nbs")) {extrap.box = c(xmn = -179.5, xmx = -157, ymn = 50, ymx = 68)} + } else { - if(select.region %in% c("bs.all", "ebs")) {extrap.box = c(xmn = -179.5, xmx = -157, ymn = 50, ymx = 68)} + if(select.region %in% c("bs.north", "nbs")) {extrap.box = c(xmn = -179.5, xmx = -157, ymn = 50, ymx = 68)} - plot.boundary <- akgfmaps::transform_data_frame_crs(data.frame(x = c(extrap.box['xmn'], extrap.box['xmx']), - y = c(extrap.box['ymn'], extrap.box['ymx'])), - out.crs = interpolation.crs) + if(select.region %in% c("bs.all", "ebs")) {extrap.box = c(xmn = -179.5, xmx = -157, ymn = 50, ymx = 68)} + plot.boundary <- akgfmaps::transform_data_frame_crs(data.frame(x = c(extrap.box['xmn'], extrap.box['xmx']), + y = c(extrap.box['ymn'], extrap.box['ymx'])), + out.crs = interpolation.crs) - n_dim <- floor(abs(plot.boundary$x[1] - plot.boundary$x[2]))/cell.resolution - interp_raster <- terra::rast(xmin = plot.boundary$x[1], - xmax = plot.boundary$x[2], - ymin = plot.boundary$y[1], - ymax = plot.boundary$y[2], - nrow = n_dim, - ncol = n_dim, - crs = interpolation.crs) + n_dim <- floor(abs(plot.boundary$x[1] - plot.boundary$x[2]))/cell.resolution + + interp_raster <- terra::rast(xmin = plot.boundary$x[1], + xmax = plot.boundary$x[2], + ymin = plot.boundary$y[1], + ymax = plot.boundary$y[2], + nrow = n_dim, + ncol = n_dim, + crs = interpolation.crs) + + } + } else { # Default interpolation bounds for sf object regions + bbox <- sf::st_bbox(st_buffer(region_mask, cell.resolution*5)) + bbox <- round(bbox/cell_resolution)*cell_resolution # round to align rows/cols + xshift <- ((-1625000/cell_resolution) - floor(-1625000/cell_resolution))*cell_resolution + yshift <- ((379500/cell_resolution) - floor(379500/cell_resolution))*cell_resolution + bbox <- bbox + c(xshift,yshift,xshift,yshift) + interp_raster <- terra::rast(xmin = bbox["xmin"], + xmax = bbox["xmax"], + ymin = bbox["ymin"], + ymax = bbox["ymax"], + nrows = floor((bbox["ymax"]-bbox["ymin"])/cell.resolution), + ncols = floor((bbox["xmax"]-bbox["xmin"])/cell.resolution), + crs = interpolation.crs) } - } else { + } else { # User-customized interpolation bounds - # Section to maintain backwards compatibility for data2raster interp_raster <- terra::rast(xmin = bbox["xmin"], xmax = bbox["xmax"], ymin = bbox["ymin"], @@ -136,7 +150,6 @@ interpolate_variable <- function(dat, nrows = floor((bbox["ymax"]-bbox["ymin"])/cell.resolution), ncols = floor((bbox["xmax"]-bbox["xmin"])/cell.resolution), crs = interpolation.crs) - } @@ -151,6 +164,10 @@ interpolate_variable <- function(dat, sf::st_transform(crs = interpolation.crs) + #--------------------- + # Interpolate + #--------------------- + for(ii in 1:length(methods)) { if(tolower(methods[ii]) == "nn") { @@ -226,14 +243,38 @@ interpolate_variable <- function(dat, return(fit) } - # Write masked raster + #--------------------- + # Write raster to file + #--------------------- + + # Output folder and filename setup + + if (is.character(select.region)) { + maskname <- select.region + } else if (is.data.frame(select.region)) { + maskname <- ifelse(is.na(pre), "var", pre) + } + + if (is.null(outputfilefun)) { + # Default location: [coldpool]/output/raster/[mask]/[var]/[mask]_[method]_[year]_[var].tif + outfolder <- here::here("output", "raster", maskname, var.col) + outfile <- paste0(maskname, "_", tolower(methods[ii]), "_", dat.year, "_", var.col, ".tif" ) + } else { + # User-customized location + outfile <- outputfilefun(maskname, dat.year, methods[ii], var.col) + outfolder <- dirname(outfile) + outfile <- basename(outfile) + } + + if (!dir.exists(outfolder)) { + dir.create(outfolder, recursive = TRUE) + } + + # Rasterize and save to geotiff + akgfmaps::rasterize_and_mask(sgrid = fit, amask = region_mask) |> - coldpool::make_raster_file(filename = here::here("output", - "raster", - fpath_subdir, - var.col, - paste0(fpath_subdir, "_", tolower(methods[ii]), "_", dat.year, "_", var.col, ".tif" )), + coldpool::make_raster_file(filename = file.path(outfolder, outfile), format = "GTiff", overwrite = TRUE, layer_name = dat.year) From ada75587268cfad8bbfe639f1a75514729be049c Mon Sep 17 00:00:00 2001 From: "Kelly.Kearney" Date: Thu, 8 Feb 2024 11:54:34 -0800 Subject: [PATCH 2/6] Added regional model workflow example (WIP) --- analysis/cold_pool_for_regional_models.Rmd | 208 +++++++++++++++++++++ 1 file changed, 208 insertions(+) create mode 100644 analysis/cold_pool_for_regional_models.Rmd diff --git a/analysis/cold_pool_for_regional_models.Rmd b/analysis/cold_pool_for_regional_models.Rmd new file mode 100644 index 0000000..538fa07 --- /dev/null +++ b/analysis/cold_pool_for_regional_models.Rmd @@ -0,0 +1,208 @@ +--- +title: "Bottom temperature and Cold Pool Index for regional models" +author: +- affiliation: + description: + email: Kelly.Kearney@noaa.gov + name: Kelly Kearney +output: html_document +fontsize: 12pt +addr: + l1: 7600 Sand Point Way NE + l2: NMFS RACE Division, Groundfish Assessment Program + l3: Seattle, WA 98115 +--- + +# Introduction + +This document provides code and documentation related to regional model-based versions of the cold pool indices and interpolated bottom temperature rasters. + +TODO: more detail once finished + +## Setup + +```{r setup, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} +library(coldpool) +library(akgfmaps) + +# Global options, mirroring primary 0_update_cold_pool_index.Rmd +fig_res <- 600 +proj_crs <- coldpool:::ebs_proj_crs +``` + +## Survey-replicated data + +We can perform a rough comparison of the model output and the survey data by simply extracting the model bottom temperature field on a specific date of each year. + +TODO: describe survey-replicated data, and how I build it. + +```{r sim_paths, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} +# Path to mox-hyak /gscratch/bumblereem remote mount +# Note: this is computer-specific, and will need to be changed by any user +# trying to run this. + +moxdir = "~/Documents/mox_bumblereem" + +# Survey-replicated model data, in .csv table format + +#simname = "B10K-K20_CORECFS" +simname = "B10K-K20nobio_CORECFS_daily" + +srepcsv_path <- file.path(moxdir, "roms_for_public", simname, "Level3", sprintf("survey_replicates_%s.csv", simname)) + +``` +## 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. + +```{r masks, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} +# 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("data", "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") + +``` + +## Cold pool index netCDF file + +TODO: Describe the coldpool.nc file + +```{r netcdfpath, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} +cpoolnc_path <- file.path(moxdir, "roms_for_public", simname, "Level3", + sprintf("%s_coldpool.nc", simname)) +``` + +# Generate bottom temperature rasters + +These rasters will form the base for many of the following calculations. + +```{r raster_calculations, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE, results='hide'} + +# 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) + } + } +} + +``` +# Calculate cold pool and mean bottom temperature 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. + +The resulting metrics are saved to an existing netCDF file that also includes corresponding metrics based on the model values on July 1 of each simulation year (those metrics use the same masks as in this example but do not involve any interpolation; instead, metrics are averaged across the model's native grid cells.) + +```{r} + +system(paste('ncdump -h', cpoolnc_path)) + +# for(i in 1:length(bottom_temp_files)) { +# bt_raster <- terra::rast(bottom_temp_files[i]) +# bt_df$YEAR[i] <- as.numeric(gsub("[^0-9.-]", "", names(bt_raster))) # Extract year +# bt_df$AREA_LTE2_KM2[i] <- bt_raster |> +# cpa_from_raster(raster_units = "m", temperature_threshold = 2) +# bt_df$AREA_LTE1_KM2[i] <- bt_raster |> +# cpa_from_raster(raster_units = "m", temperature_threshold = 1) +# bt_df$AREA_LTE0_KM2[i] <- bt_raster |> +# cpa_from_raster(raster_units = "m", temperature_threshold = 0) +# bt_df$AREA_LTEMINUS1_KM2[i] <- bt_raster |> +# cpa_from_raster(raster_units = "m", temperature_threshold = -1) +# bt_df$MEAN_GEAR_TEMPERATURE[i] <- mean(terra::values(bt_raster), na.rm = TRUE) +# lt100_temp <- terra::mask(bt_raster, +# lt100_strata, +# touches = FALSE) +# bt_df$MEAN_BT_LT100M[i] <- mean(terra::values(lt100_temp), na.rm = TRUE) +``` +} + + + From f8c711e14f5fe44fbcf62120d0e891cf4e9b4c5a Mon Sep 17 00:00:00 2001 From: Kelly Kearney Date: Fri, 9 Feb 2024 14:31:54 -0800 Subject: [PATCH 3/6] Moved B10K_lte200m_polygon shapefile to proper folder --- inst/extdata/B10K_lte200m_polygon.dbf | Bin 0 -> 85 bytes inst/extdata/B10K_lte200m_polygon.shp | Bin 0 -> 27204 bytes inst/extdata/B10K_lte200m_polygon.shx | Bin 0 -> 108 bytes 3 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 inst/extdata/B10K_lte200m_polygon.dbf create mode 100644 inst/extdata/B10K_lte200m_polygon.shp create mode 100644 inst/extdata/B10K_lte200m_polygon.shx diff --git a/inst/extdata/B10K_lte200m_polygon.dbf b/inst/extdata/B10K_lte200m_polygon.dbf new file mode 100644 index 0000000000000000000000000000000000000000..6b6f261aa0044d94e62a1f037fdfe4f7e4430df3 GIT binary patch literal 85 zcmZRs;$mlHU|?`$5C)Q%z$Y;&H3uT>45Eb4l<+Dzr50u8r5hS}$LEx!#v2(JpgM^Xmf+U-y_m~m1ZFE(DBc{k8Fso-qxWU(n>Q^lr^ zO$(boHo!*0W{J%n8v~mMHXm$$*rKqd8GB}&TTKIB|GD^EG6cA@SyCnn=}*UA9r0O1 z0KV7F8bcg3$oczt3ma*V43^<1q-!^Pw$8);|NZ@c)5_zU#Q!t?j;jUaepLc|Uw>;* zU4RCxZyQ5;J=e|J^^uzoSDF`!`x$YZlu@0K1BAkwFm}cnq`P0tWVinh^EK)fmB;L(!Px7+VfT57aPNp%sm}o#q%swU?U4R*Lqhog z(vR8`Pi$Qzz>SR;zPe%m>uJp&&mjH3eP1S17?h6RJN;29IwwJdvX6x;B1tsR@{>EQ zjC9wHzR*KRpRQ<}oag!f&M&>7HoKvQ1|x6pq%KdAKq)-o;-w}Us9pN-{?`NvF34ST z6+&7n^F{vj00|bSBn;E5X>jnhRmH*}3GTVM(B+=ffXsENX=#WA?{4S6&woaP#Xik# zbMHvt>JVC|@fh{TubeX7MS?4~?fex*G-x$>I??~01O`9&`OT2_yjsfH-AsZ|!bbTI zw`frC_*{L(D-xs+s$6AWqe1w9WaO0=5)`R?UiUAT25awpbx*29+W3r5@)a6zTDxSp zRgqwh&-3$Er0Zox2W%gbU}N;ux|)kLDB^EPm3u@2m6_RKSET>$*k-o&774ztJn71b z^z)CWMps=S!L8({Z?9cIc`TSdu0uNOwN1Jy(xo5XzIhx^g2?vwUq{c;;9{cdF3$uK z{OD?WaUq-rLs{R7)+drcBG3Dl9@34UJil^BkU)U!LK_UEL5)BH_g$pr6B7kek+v_~ z*O~55f~4(@G6#^p>%8j9StbcSf7QPH_aqIDeAIUKaU{Wu(EH8zj?;kGvNo9aFbR$+ zR*ktMEnT?yWih&9VtiNA- zl#8@(vY)bt{D64-rBxt&_~8(!eBmO-&)vkj6}5BmFj6XWxYpB1os&(Ds?rpjd;= zHHGw}k+_9*#x#($iWWB;CBlbel~#cCiH7|=tQjJx9EKcG0~+KXin?5m^s@>lVO^yC zDbXK3tVH>*-CW{`beCAK&fFRjsBioic^YXD-+fnnBMD@=gL2DpJd4hm*cutMmon9k z&qxP!JsRp)B|(ZeSzpiu^^saaT+k%}`KWcT&VOiHR;{=R+TZUS`%nT6{*;94lpACJ z(C^0zkyd%udZyEY1T7n`PyDe!`xdYs@J8BuTS3$oOB#&RfB%U(g7VoKOmjp!*lmMw zjtAQR$<}b#{93iNIG)IF3X4(^W0A;D7B z=Iv#kG%zCf>;KCpfren*@)FX!k}OsU7m^^_kR6bQ{ll&H+h!o$e|NKCH_}2kcDUUx zBSABL$#nyR2CIF{9LnktFH}kOab7g|TpjJa9ce`}Ja0$(c2vZuUI&g}Jim6mH{##v zSla{rB)GnCp#3P)%Hun4T8@(-FYcDfrxP@o7A7TEd?&%aJ28{0{)i{*>s=b=NN_|o zI{fY_8bljy=-7gEO!4)c?GZFM|GHn{%q+@tEXkxO2II%!N!vFwB;fh9f`c~}?ct2X z>-Q_k@PKpO!}vrRaBbP=D}(ggW%B`%G#coH&Hl0HBtwE|o=DLp8vIZSXIuTj@rL+X z&5?fZ?paD0A%TBr0Ie<)@lWA`*y(R1DCsI^8|KkqZZG+8;a8;NkCzl8O|3aD^57S) z|4C({WIhep+g4^j;UGiSx-YAvkp5|!u<%Tf41LYbJc4&&Lr$L1MHG8|JX4v9co@=_DetP&Z-$d{Z%%FtfUChBq_9g`QG zaP7kRf+>mdnFB#1HGZF%AtfNEg!r_39DdE=NZ$ z?IDBjD2LaZ7W`gQ^Xr5@8I0BQpBN(@_E7ux4>dB}jasUI`U>s6s=4O392r(c`DUmi zEkX3Nn^{W+HCGp}TeUQh4JhX1k|2W~f9C~xq=Ujb%UAM};q<==_Z*Bj$72f&ITuOr zzr`P4GtwEEW*YRkEPeX(mjnl-sh=|6puadcG6v=me~2=E!bl$tN#yNE`}yzq_22!E z{ty*XX~+1nUnF&4H5vAWCT%HvkMTm?+4kisGVm1C1gu9og{<8367B21zn7P*XjM~p zhe)j$rze}q(3q?(#-ZU3J%1W=ZfKEVEl<_f^E&SEAk!w3b0-;;nqRp-)^>+hgM$v@ z%4GQO@BjPvH-XKs^!B&|pWS0%qyDtT=1N@j-NCXJ9*3+W!@pZ2wF6r25GSSHtc>>Y z-{0@y+~S-ahs!FSAI<(V5~%r6?$d9J!cj^|Ji@c|i3 zzg*n$_=Y<)2b?MpDJH|snN-W6n>b#wnp)3uGRS`LG*7wb4vXExm8s=qh@$c})fKtJ z(95!e4h>{*n_I)RaNiw@w(r*6_?8S1d5x+} zR!LC6PUF}@%_8b2-9%DkEd^+oUq00R=MKGFK7M^7O97t7;=mQF=wQCrS8_>)0`=aF zf+4HvU=ebRo~ufM@S}xp2L~Ia*r3CCkOF*H zW#8H;(?K@#poWwQ1$MQKQa5j>!^UvlkE%u#xGEr%QmuybTON8-L8ic+sy44QZ93#J zuN|r;P@qIAXqd5w4z-11$rU)Br2*eBo&9u(pmayx$MLSY2wsWbM~CsDi(BP)Q=lb* zlxTC14q3*QwE{X6IMP=uCPkqG=yNGdYEpm|b4Bix4IM7L87I`MQ{cHj*TbB{bdXT7 z-+yTv1sDeoR3CPu!y-FZ&{2f~IReMkW!&jdx@%wlCv6HCq#if!^PofP^_TlHaJ<=g zk&f$LbU4Vm^-j_d=kt=7^Y)>Gq2O_P8m`|kUHG~FaXO?`mQIP-QQ-97Oa*>FI#@py z@H$VUfOFu9tGhhuz~O#jvYSBx2QJaotG($^@(t%aRUbCi2z9?-aP(!p%Nrjq=D#Gpd}WfLrc# zygrVXoVKfZ9VZps|5z2vqW;_-AH7qsk_z6lDZ!&u$+{H_sGYGPE_#mlz8W<-abe|kA*1gNkrfwU(`pAJ{Dm2=ALlt8{R@SKHrcN>G8!!JpZI_Dl@i9dHll`)5%9DoO+O;WVc@B146Eb?NJS zHR%vySXd{gNCmsz&8o9I=}@gNZtS2+h4eDQG^ILb1H#NyG>4 z;{5tOR0u6kJWCCt!+G^&$!Q}heC5ok&OA+rl{32tCi`$bzgck^5p{n%jggFEQ{i5N zqjFmU9lF9L&+PT2Lcuq;O-G{X;5;t*Ys!lX(apQ&J#jsUz4T5ApFsKVpH>LR^>Fup zD&HPVg(EK1C(Ds259TnV?-Uijo!os*IEoHRh3SC{XQ*JJq4h;R2JwrexK}cU3jRC& z&fiO={lxc5Gb_i`MrbAnie(s6mRJi?p?NnwF9cIhE ztvT#Lg^G+j?~dQ4!%y25O(Sf3WDdh+>qIxzCznT6x`{$kJ{ONXa^=DU*7UhWXPe`T&AHt) zN(C)B4ufO(y_JP*mJ?8+*u_>>2*3BstI4P6q5hw`H)!DZW@%<@vt3m9+<2{;BFX^c zqwW3%=wAa@gRg!QWqJrwC_QKmH40v1c++i2ut={_7>)AvGD7^Y7xf<~?Z&IW!-iiT? z!cxPwwN$VnK2Z*KV1Vi~o{o#6R7m`|`j)E;12h$6I<3}Fq2Zgc!9ON%owkb!{tU2<>k^vYiuf;R<6#uY0E->%=Z1G+yde~8??29f{CDc-HBp}I zVZHnrKL*SntXkNN`Z~Ph*XE)C2CP40HY0`lBpiR(!U|%*NhVjI5b9f+%lnE-2m{`P zD~By3e#t(vR9cJk<;;s;BcZ)WizUT(o<#jyrCc#UdwL<`yX|r?1Gq;Iu2pkDf6&{> z?;Og2z+DV}X|&fp4sO;FXK}o#o`fd!mxuQ}Ci4LH|(m<`z7a z#{iAXi5uOS=ns;o6R+H6z}e5;=NHib9ZKa#;_fm)VAyu>BjUqP5B~b_`wS3W?V3}L z_>re5u`}oa1A?_Ut1h2JJWATV`{ZK=Y&d1O-!ByLIUsr{sEh&eGtx?=NTe^`&n&ED zz{(*9llP~ophOSJ@2+P6hhOuv!(mi--6|(N+QfjuO>%uZ&SLx?b+G=|!hmgDzpYjw zJ{I+kX*AX`ApA3RDFX4+NYC#;{0jzTh?dMPBv66l-rr2p6SUuME2Cq{Xy1hEn!Lpf zNXptbY;Fy6%{LvCB$omSDKU`qIo}p#YvWK`{3&U6Ki4541A#a-c6y^DL zAj>WW@q_1nz{oQy5NEa*_MBz_E%-;X?h7g`KDS@6MSHU-&@_ywLjONon(_wi^*Ob& zt)ZF<8d9YJ#%u;0ppFFc)KMX*-1?e31N~Dyqg%0(3S0w`>wBCTAf)+N&iN(g3-^xv zHnU?u3h`1{Y7OQmH%oLMQyGx^Is1DpeouVK9F;OdeEX>DKHY$Lw`bGHAyc$JZvtZ# z&L^$E*8c_i|5>NAr&VjIVE;iUT22@3i|gez^Li>Ad;LoBv^oQR423ngH)4LHpljWO zc*JD~lc$=gP#!oH9IMIzX>Iuk&o@-)@pxJCL7oAfcMh)#@1R2W)ar-PN(`uXqA`5B zsUVZRrT4WA1M}%F!>}GIyjc}~Cm8Xr@X}>3?*S_K`VmMDSCAL@`o}t1f{l(aqGZ-%?RwYvBFkj#oBEO%+ ze5B!6uHj!QggCu7W&WUpPW*j;Ee;DXiTs=4H%bRJPTw1wI4z*-_=AMtK{|Z8U3hgB zmj!TLb1qnk@yu>lap(fC1@zL&#W;E~elKgp^zmCj;ahX%IgF2+W7t;Zf)?O089Va* z6&>q=9KFfH7U1>5G>P6y2YY?)GQ-&3t*T|b2sf`K<}M@q<1SUAiPkRa1$8tXE5N(lV!}0I=yPb4r06z zJ=v1FM1_wJ&Doad4-un}zIiTCA*!2qw;cM5ob>l8!+FeK4``}8w~JbuJK}{{V%o_#8n`E%T}0OU5);6E8Lb&eZ+vhP1)nw zeVC6$GGsqJLws7hY^U)N{f{m8h+e~h_fOqMD$w8fbf=`pni&vFuFEm)LOj@h_z<&$ z0S6mpPd9hsdh^yE9{j)nIxQf{tR44nO~WcSgADkX(9~Jeiuio0;#tiI;zN@xd3OsH zE_2-x+&;~K^o!?K-f2R-8~-4cJj(#G1=pN7;)`cR=FrqXj9)wNs)r(;ST;KyFfOlV6pA@`Nv@)+W& zN{FR&umTgl)NQHNdxHM?{q0nlDieh7T>f(I5#}q)r`&(;U;=;f0Idb_ne#_{u<~vu zyp>lr6S|N2-dpc&Hoychzi^3gmkKv+as$H-FyXuPZt_{g`?SE9hx5ox2r{xcS%&d| z$2y^_%7zJlU-OEOI2@p8yyETVspKGKt<;r_#wlIFgW2`BVwI09U#@SrH8;RqiSyf&RR zZaaehvQ}@Kl?W4L>}z-xG5;yxe)vyX2>BUXt>|QUTKV00zgn`h^FJ}&-y}UiV&P$mIUC$Gh z8cnF+aZ&k#2HIQPHz_t3U_M~GL2(7zn>W>I?7&_sG{19RYSCbVo0&3CH0G113MWG2 zx1;<>2$M zDsWSrP9|di5+{?wQU(0}V#~-reI{^jvw1opONA$6=DLvvOt|vxx>1ob6`tp)-qJE* z!s$X8+kWg{M16L25bf)zo!KSRQw?4+nL z;jFLy%9#l>jEc@cNh&ybr|eX4V?v;_w7xp_?=C3~t70)>X=B)aGwgr$-`AJ2UQ8HU z;j6od>mfzIefr#w36ghilh$DWpB@dyN+;2NTk`8GaXq}YFS;s2nb3GsaqKnrzn-pe zSK|y5hHEYj1mXI>PaUXhie`ez;je2Vu>Y(V%s3@5!RU0dU=;Ru@|f8>k&67C4e2>2 zgY_E^X6DIECe|h6Po3F-_C0i4aOoNouD>t*9f|ee+=ia6JOo)y*FTTAA?ZJ7Iw#q{$6f*d-O_Xqd{HM6?WhoPkZ66f& zt5PB5^f#t+ITIfJoajEOL4}rEO5ew-nBb6Ww($t=xBgD`^(56YVfX59=|54ws=e*> z4X>DR#5r@k1M4%FwvXmFwld)(ulPMTZS=prvdYwVsQ()vY=Hewq-m9`{=kIisiH!! zbdi6}EfN>{@q6Y-*j2Qju*VJ{@RbR@$&cUp?7;mqf7fNo5EEX_4V?bE73*cIePq)= zF@aTiV7DX6v-fVyo8cZNh^@5ZZ2gbCucg&^zh%NMEeZ8Gld;P?Ne`Qdc8B$EED|MScwC=u2EHu%18+ zKiD+Pgpf_IKjtq`KyH&ITM_vwb`YIBJ57P?5}{1L5ya1_3I7da6zI}l-|^-P6QnLw z(ic8aK;wAuwB{fa?2dfP=;sACHp%^B%`fyu{MD$AqHLaO)zhzkcUn@>BL; z-t8p?q|DZQ5W@a5a)MEM&nPhP=9N(v+VAEk;>y!l-;E9Quu;SQKgdV+=Uu0OF6S|x zx*OP^PyY)wg974Wo>oWmm_Xb+C$t!k>l5Nq(pfK)vIuq!htvjB@ zP~gqwrlTiPnDFM|ulKJKC@|8cTeA5A6DHjpIy+M-;L+G0(sd5U<58>Vyn^zkSJa*k zXTn!c@q|}ckAB+uDONd%3FkGQ8a%s6fvnpPJp23+FFh}E5HhhoEVe5hmKOyaBXy_`u)leeq|1m41*#=A6C<(zZG|_nv_lkF z$(XRnM@3TtPR#sT=m!4PZUt`Y%3HvUN`kfSLKlOOcFO2sh zyb(#-TTq_dkscoqKd&U_^5jWV;BOk~X)Y1*?T1RLi#*oP)A!gVm@#2blXqzYp6BSF zeIb0>lnD*@ZQE4&DG;-c((S(={W-bdOWZOU`eTCc_+WfbyuWer+b=SlcIylC+RKE# z;$u=qc%J1a;Bhex`;Ww}eU&mqhHw9#9=L@4c^c^I6J2E3=AuaU(PP4k`{aWecpi4% zSU|)d;|J5sHD=-s87gyk{kV+d^_F^HH-1Hi0QQrZ9PGa}xCTMtz9WH+?huPkBwL5 zC$N83>Ck~~m&wp6EWD)sKh^Uyz{J>DFLn4lY|nynj5hC(|>nMFq? zwEHfU&`yxy@sx|_k{!l>*>SyECK))d$I*TsVS;ZwwQJm!3^}$#MZdAXOu2pq(UuHC zE%A+ivA=9kY2$4I88poHNAO!Qf$bGg8f-y^YrRkOC8SmmHg;m<%V}IGM7nkaCrxpU;q;|k8$v{XfKldd*MC2-7#CSEryZGoQ3ASI8-!_PNA?aPrm? z?dMGR`89LOcLf=8>*~G=K4rp(-sxZ5cz=X3cdqfp3nnP`ceOv7!uv$o!L1$`kLxdg z^>r8}f%x@)(fuzG|E=~I{_4Z^>fh@hMSPE%iWa%_9@mp1l3iVc@hz-+zZ%|OGy7t= z*R6>OtNq_oTJU~M5o|uO7yC=tTYPsf$NOl*b0?&*f55Ou!eS8#;yu!KFIA%dNp%<` z<&)s%gYjoWxZX31TR7#?QT}@Gi2>DUKN|)E>+pWhFs)H(bt4llX@`7pJx_wGL+2XC zaJ_PO$R}pQNf5@%ed9zk@+b3Zb8-*~E`Am+b!b3;iRs#->`j7aAD;RdV}H9^6`wL! z5(x6Db?-p_9^A35Jz$IXcjO1Fi=HuIy@t$8z7q*LjUxE+%9t=NRJMM*CDO{deex(T zDq&*ZTLKBDaw5fZksq4tI*${0Un-^jK`@l z{TfEgn5n}=z{LpPF;6Z|c@lTZvUzw2o_Ft7-0LtU@ zi#vZtF&_=&*jaXl1XAiRPWerty_H>jqIe4LKY4G9nfQ(V^)sFlcLMKo-Ijgh_80dH zm;6_Zdg1*=wdYn{%ZN{=?VZ3N!KLG>VkWCtaH@0j>*S-T&k-BrGdwJ~p!jg29p|44 zQ#d#&$N~qA;dDuN64;+%N7#z7AfSO2V2Ag?!Wx5rJ(pww*Ui=o$tdslq@wa*aTXY7 zSbT9L;(bfUL64pdENBhvE>A{zY7SaE=*Y4_(P6Q>%76qXu4FI#l4e2qjZx{uop?V{ zuvFh%o(1;Xr1_t1BjLRO@2CVN77P{bHyPC>fxesA#jY(Zkp5jqU%Q_KR@#hHAqp&r z^YsavGDCiL#?b1OST|eN9`l$VN!+kpo^s4sUVf(W-K5ta9OpYymuVEQPxjpLG6

zSs!JAAms4eM0+#+ZDIZOFn*sBT747cG1tDC7-G+YVCffX*U`R5A{7)R4za*sl46yE z`YN+Z>w1XuW%A1CT+$UclZ-15<10t*7gevyI{NwAgQRC$jD z3%u`-pYW8y-^Fi7d?um&To5yL5+%W(CU#zeISUq5r^d{=Nf4Og+}eoay|^Uw(0~K) zM{PVh?oVXF83p^wZFt{uZo6!?GRnhc;a|bJF(N!Xsc-(#gaw~3Rr@c0B7)?UVT=II z&tms^&%P(Zm5wzlY7Vf#be&*Pe+v;d{T-DJHD-Z!{qTp{mqe&JA2p=6j|IY|>*sEk z5+S|eHY}okBF@{dQrxHVckr*k7=_MT7$9(4Veo z5Bb?=71XkbAot{Ly@omqzI_R=<47jLVflCk&TTB%@pIthn{!04e1CrBd$g|`Z*(C$ zmPfdj{dbT%9N|wnFyC0*%s>SSil=i7>l+h zLdrw#c|lPY6r2^Vv2i4V>?q%ZSs@l|nGf?8q!Pi!W&F(-ew_bvigmRy-bY>iv$$q8 z3uGW~EP_CU%1=)?OINZ$;gN2^4j_WV%JSku4)ots8}IdL5kW!oY+u|0=2!JzJjwb* zpuH6?KR3sOiMoLq9StG`N$nnu{Db&EsJFXE5&7G)nSSmk=KptJzYShbgsf$rZE=WS zkFtkbC54EvK1^IEWt<7mB$~H9;v~X6>3n!L;@v&ZieJQ80tBnuhn@Yxgq3Ow#luSk z_$+hMh}n;Lnevw-kc$YW>0)Ptdzmn+pL(8qo&dXb?0+3b{2j8d7f)IyfNSM->IuYW zN6wD=^}h)aG#7co2Jw4(eOp}mUjmG~kS2~JzP}pb^OBw-z~51gHVce52gc5@lmFm- z{C8?@ei(1e?i9rcj}hQ#Z;*f`#=C9O?v=sQ1b8@M{wxsVT@RCE(eO|a?%##f&fjPx zfNK8L`E_@(KCriFlSl;tc$0@@sJF4+N4@0~Tta{|1(BrlH!*+imNp-}K>*Wu&a+MV zOo+cPIA(H%04t7{JuAWdM3i<>FD{n=dx|bk{?5dFi6(J$G!@5Za&6MjV1kBSud!(? z0V1PSEfiC+-V$@NG&GR_qn{20`Xu9i=a;_byGR1!#gm7&7nq>qBtOR)M}Uh8>2Geu zV*OJ3oVZOGu0K;|W8zuNCk}4Z2s%rE^J5k2Uoii=AX+tF7f1k?y;AATAsCOMV^{;& zUyS9Jtr(2?VJY9(WcHv{gb}QOqHF$eDnOtIVA(MpMz^o z*qmfStZi!TXEOqf7?VV3$bUukN6WVT1elk3=({Q$^PBv0rqmI}l(dyVo@o z^VQD7CV{@z1o*m{`{?#GChT7yuP=3!03=0`BfC)l;qzIH0eiIHjuQQeY^+xfc(q2l z6X2xMJ+lGKZ_j`F*7nka0A{grifPz?U!<=F{UiY7H0Ha(7J);z;|ZYq zdc^D#_7@pdmR-0=fbA_5-aPEzb1lN`;#C58N`I?(jQQ}MZ1qlq3&UX3Vr_n#8_pEtaj`^SXKxV`R0!*#( zx%UbERh6Xtbwd=&GsD*CEBf=-wY%HzqJ5c~%3Uf&e?N3~=!E8Z0toGY`#To%PoI+< z&mX4|pkeI8)v!7y*dBg)NHv=PVOO&>tT2CFPo6qcR6qdbDS`K@SnujrCLGjwg!*-< zAK(23?bDDU`}iIK)=aC8n7qUKiO7S?J8uy{ZG)e?VmIQa;}hq*4+vl$KP{u)$Akzq zZgsi41X$m}c%n0e_7SFkD-r$MQ2#=rHR9b!e`T-?`nS$Uy1wlfl#kKhA2BGu@~dq- z^uOW$|EGBCYz_ewmXzZ6Af5{3+ugQLNBM3l@l-;*)qUpuv^t*v)c!l zcZaC)ClII)Wm&lq;?$7 zXO-UI66UK?@)~z~3D7Pvui1z3Wq--!?8*`BuW0x78|J5Pzq_nCI7ooadLr#3n4fOm zoAN$$9PuT$_}352PbaFbL)S0?#$*O`ruZ=bx@FyWt_STm>`~biHw&(>;O*9ILwxZ2 zF8fn}1@DshYP=9%PkYol{1V1|=WYqDq7(foiMw+ai4!h(?S7x$e05uo%AfBzme7A)^Ca669i?b3;kuZMQ9K=*`V z+)a#se+hF3PGEkyV$tziw;&PLHSOZQq|1WKO`rXFFn&foJ*4s!^Ua`Kf`&H6*I3uS z2lp{w(X(~7P~S>~Hw}B==4r4%{n()$3YtXNSns^5VJ8ciZ1Yvtdx=1}`aI*r}V}9pxMzHlg<|8NeFLI4xKAE;M z{qhv%E8O1?aBabSGv`JejSKS`{cg?`I+(xqaOoNiyAffc*IFs)2n*8P_eB2k#(bbk z{e7hi3t}cS%l`!tfwx#D@HL$UI|Ys$?C{0+Hdjk5*0WK+p{;{mp+wldDE>py2iH65 zSeP4t`ONw9X@6f9@Gn|t<(whH!L@L{*`Ec`qFg>17cjqyDyowYVnO@qVwvR2_&&{e zHdyoo3wV?+3w7jRzOuFEwdY9|EN#eAZN823A5VSN9L$2nxA?mBfO(jZ;EPpNWtw0I&JdSU}l*JTt482+6088%AWXAeKAk zQNTAMNNA0OEoI?;hr_~;^b_A3xlJU69wShUaq@ z1|1a`Eyv&Qm57_@;C^jDd-&`N7L18}_FDwpPglB^D^}xvs3q61823*$7krZ&>rfs8 z)xJk@|K#Gu;WOHd`E&gys|GR&N{^S7q_m>GDB5OalS-|~h9^+uQV90N!&UVGB#^jj88wi>=(JdM9MY?!cVV}Wbn z^xd2od>=%s)6aj!f-F-#`}t6;PsBKK1=Qns3Tw5K&yhgQ%eBP#B@4Fq={Ly7i~M!kJgit~!g@^X%tjjWcV!z#PR3gjj2(Zl zJMTT>o2*pxEY@ozv+{O(d}P64g}ZAKu%5GCus(cOKMO=JH-yO!VSVJCgec!<7N{Pb z7rZ}?^{<+=nf|XVsH>`fvU`>UZSItvvk&#yV$Z!^tzs z%KIZ18+3me*6wy915w^8u9t@mk^^~DR~clepZj1iB)|sF3BSn|zGSfea4uzh4I2hr zho3r!kb!Zo>HNZ4Hq1E76P`wqVTJaZyVMPAc-?VOvMY=XhRj!K?i<-4bImP$b1c?> z<$ZR9%d?@0|6$gNi)47ZBQPRcnGI_>NvsNdFS7S1JE=;Q4bN5J(~1J*UpDopxCR>z zKbcJD$t8nCuVN})n+>UE{Bu0%WT-2-U7xg@4O*Lww(Y}uvD4t;ipsrgFs(J1J#iP` zPX$zj{o2QdO#20k>!oCPcgaFzn<*RG^V!;MRalSx?Hd+CU_)!Et0c2Fs_r~RQJMS11QB&@$Kt+_UD$OeXE>&uyFT>r(64a&xBc$DQLx;cdmYmP|VWSOzy`mLsn zHF5a<>w`)_0htYNl|n!NIzxtOuk{H-R&3x&SULM9hztw$JLe?q*kC^RdTAQp`%p6W zIg^~&aAa5OM9nEOB=Kz(3U_6Lk&tI;k{{OdQ$Jlhc!UkN-=|~*vhe-W%`!#l-LABM88{U7?HJqxJzx9hXPwfH1|btoAY7aCfNblG6M z-)bfX*E1hYonmRCJkS2Hi^uhTeDvByb~_s`%nbURM}Exslz-6bk!8b)L(b2OX!t&AL(=XWW8?do z)|bYv_+Dpdr+=^*8$!l&;>NI^FMp%ud$IC;NILOfxue8m~Pk_~Vs=dc{s*Ht4mZ#+T&UCHIEYc`DU>2iZZI?%s|9^|He z=)wA8Oj-02FWP_cTjAvntdEiM_bK9hG+)x;jju@X`1pK3Kdwie>jO3F9SI`Dw|CA; zu%Y+Ukl|QA)^EC-Auo6y}CVpNKFKgp%~AAUtXF6XosVtq6@=tl7}*6)aS-yOP!^|j4W zJ6e%n0uPvlXON(wK05jv@>|8q&$)*6aO=-5AytBG_~fy@P9%#2z6_n#F4T{C{e~T? z7*BE^uTy=D{Ek6GY6`xGy+DcgY?NZ-dxPUr*V0JPwC}I#-A!!xlGb0c0_SVKcUdJ> zfel-jYt4?}dY-;8^0~2v4I_2OC7iK7PkGAKma+}?;XC^J;XUMs`SRph#IuG--;-R( zuWFTi{zb$yHf>$p?0pix+cqzZLp)2k=dlmUQJz}TtD+Fk?t1d}o~# zZ_aG6J!N(62%gsrm&zorq_N?Sy~j}K0tp;{ScVC)(4WjQIrrfC(%LiEXVtyf(Aem< z!fFZsz5ru`rLHgfXZtVj zlMv1Zu|rX>>+yZ9!A;(h7m;kJzNYc4k&g`Ps5<1QVTdPUF`Jde@IC9Zn=bc**bw+T z>3*mb8D7a=9v=!t{4o*~ic}(lF!!PK&OkPd_8rUa-cAO*UigU*fA69yob1yj!vQ6s zD$_G){{)4OzJ2JwJkc9>gs|b#hBtQF=42Q-YS1ll67hO!@JkfN+lLXqZ`t5_|MiI6 zjyp_-?KT4MwUHk?9=2)^;*;2k9&3@~INme0JOwW@oLRp@$jYA$A2L|2s)(l-cwV2$ z!1eyESN7VDcrRtF{;CrBF;d>o*^BQrHD^xx2L_|R%Jn<-A^uhKURTOMd3Z|Q`5a3p zLpgkrdwdr0=B)`s;|RVdo$^%q7sG~+NevR4nPiCUXMQ8Z<9OO@x39-|Z##Uijh4a& zx^hru3C4pL`J!TFX>4fWS=7x$eE!;@FTN%d^+?CsLZh5Y0O8Y1b(F6W`m`y z*=a@q#@hvIna&OLr{gr2L{B`A_MiK1bcYQSPtCN1P+oigOu3uiNB`O%rW);t_)i*B zAw6Qlm5BhtUzE4ap&Wg}QtAe!2@VV0KiSb=F4Ewz0`?U?8fy(A%h^zr#!ARTdt2o0od5U&<5k&^jHfyoQskvRe8csNj;ys@y^9QX4fQ$3 z_`9+6v8`=rZ)GYVe6t$k-O^9GF2Mc!1eKzaJ&W?s0{Um^DT%MVplcH^vLj2f9uba?QGaJ{pxb&4l+cym>N5G zu_2h`aA;%;=3oE!-=|3_zUHScONP#nZ}0o?cP3Mc9FFlzFx$p5tb+|C!$Ztwj9hRReV7DR2>_}r>jpTPT}7Z2>r;NS%vwM zQn8;FzPDCd`FAh*ImU+z{2!T0c0!|INk9KV(C3ZMlF3;@hPX zQ-}T{^nZo}$rka?Ua3!H4(%z?@rT&(cM^~c)^L`iJ*jJ$cdo>Ewca9>efk>OSIDH` z#%>b${1EwfDGTwx<-C5*00|5Sww`xM$M{;Bv%0v81Pb$tu{%=HzLNzHSKxl&3K{<# zMIsya%~wvUW4xrA3L36IkN6xy`>KWU^yjg)@>`-Xzx0r8&h8{ZqE>1kX9VJTL5|lI z+;3cKSJxRtJaH&4dXR$qk!D)4=>*~pHM8;0L<0$obW?QOPhfnxmo?W{L4q5kmguL) z*zmLNXJOMr#6Q`E9XCDMprU+QqU1ISCPY~t#SAuxD<{19TY~3X5xr)4ZYb~im0~@& zNKlPs?R3QBf2*RT9%tiu;pq9<6h}5Bd6Z3`zeob^LjLE;wrmj6*_Cl57S9*$qpDM_ zFh03wNo@~DeZ`#EeHHV~8O>iQPSL2Jqz%V!n6W|QO_F*{BAzd9{g_jP@u+c(&}EMM zx3*Sk-sc#v9_-n`;KlRPm;MoVUSmA_v~*^q<`~Aa1BA$0Bh0s&=LW9h-#htbQ@!oU z0mOsNcbCJD;(1`~*3-8LY*_J$+*rXzeZ3%0rD1&9CjM+!qAS`9_qbXF#-HYW3%w3_ zzGt8Cd?69z+n8UX6JU{~4pq@swB^$~@YM^W{?!PCE9I}v~ z+TNb&(48bG{7Jm}81q-}MB8jb75w`zwkCf&aKG?nIghxs5%Z;8Cb84VkIZTd4aLnQ z;5WHEDr&%nUoZ97_Tb+`c`moXZnqvAR*CBJUcmD>+Svn5*4k_+m8L}8*ns1;OBtTt zf&S`xH{-t?!ZZ_mipTv7&|EI7sji;&$-~S^+$UKxnMQI=zl1k*3`JCejC5jMdDyk=G zqznxr%2Yy$RHP^kG$^Z9vR{9dllS!?aJ_g;ISwXW}d zUmJcw(>m||yZ#ayyd!+%9E*$>9(HQQ`l~gsGhQ_HUp;30wdL4vGK(c*_8eo;h{lP??e+XVF7Yvq}HK`gqlva2bp89erxVwNiOt&%Pi zG4=Q0yDE)Kg1k|`cSy^;UO>^Du66Uk|Kz7xEVv6jLLtog{a#P-wd=3COkWAeKTrGj zb7vOTvt5h_&jl2+#_8TG#UCOr=WeI;WMFh@Y4k;D}?1(6ccHqK0QT1PrM%u zZ<~bUQ{bEm=v%Q9V`Jl_SY#;AySghuKnGa9<8y^rWWO%Q=0vQ3ev0ItkdPaU; zPqbWj^w$jZr_h?&?%@wTXSd;8_4iDgSif}ITmkfwoo53StC{pLZNrU04)n*r*Hwx= zVN&F&b7w2og9r29G6wEIPxf1v@X1p^8+`RtXO=Q) zt%Cj?tajAkK9jcXiXKPrtrCmk_GdGxz#+~|t?pyfgZIHxN{bIqNHK^072wuh4nA1(NBG%6 zQvsc4|G1Nq1b+04{Rwm$#|y(y>rHb&x7DK zi4O8Lb&~{i*JeD&8ocJ|qYTFm1p&3l`l?L^ACha&O#dN`=k>{l+zyQIK%vj$AxZE9 zQT>1fjPJnl38zKH@cbBhZe9%@TGge^mKK8kmy=|!oC@Cew~)BnD4&>fPu(gy3f|DA z(yj~M9X(Dr-yM8&%6#FL<>1}39{Xzz?u7pt+d$kHJR_^@w4W8m|0M3*BTev(c{493 zW#atg6f4pdz|$40vaCmO{v!>$mBhd^Ql%U+%)rN{XPR;Qr355w;Z?nT59W`CZBnZ| z&aWk_DQ6pa#+z3e&nF7VbX#XdT>z8xxR>SMARcFz&i)fw0&_4@~*t!Xr zRQRWEwlL0DL4I&G4grl$E4p26_jUbGZX2k+r|hiB{J_|92f zYa_six+W#W)bN=U6*uI1_A{SuJ!*F5dorn5v1D0(JD+^&_7#ZY_;#mb#B0E}iUmF; z32sa}av(0d3H&QDNjPP;7x-9;iiP+)@YF+PL-{Va-`qbL@^#=*HGK{y;GOqPv?S;? zcxhP2gktc}8|E=0k>Gi!Oh8F;DCZVka6@VNYx8zZ;DQ#(^LIpK)s8v68oxgmJ) zmJlz&ZNz(}?qPE#_3aO;6e{G? zT>kivU##&w-t^pPd7V$4ej>qZOuz%z2g%#!gZJrw)xJTPKX%3;G2r)ga+Ryp7J}!R zg-tqk0r3i5i$emIg0CBdIljx{(?g?`mM`?sKR2+WF&*)_5#Fs6jd4GghoxXQ@*bO>YwFKJ)h@- zc*@>H&SX_4oi-o8@-!Rw`*`TbA&louxR!OaJ>m^JdTVMiew*V~JEU2Bn$g%Tdc}}Q zX|Cadc+^+2pJ*KkUOy}pE9B6t264`7wIo4GeUf-M6ZSd z-haNTQU$h#i2rY0a8C#Ii)XEAW0(O2|*UHwF1WYxZ>U zlkk7`kACI?b>s)_5L)Pw!6dawUEzI`@P4x|E&6$$N#kf_RDk-MqEtjb7ci+oX!2t% zW$@(-hUxo?m?XP4;ASf7$Ha7h-u8$|8hoCY@o)_#3M#vz}oy(#&?-VyNZoF`gB$hW%Toe;ba z>j%4Sd#b%rAr=Lgn?ey5|rn9J0f^*9J zJMy)}UnG`deRKNnhQiCeJX)cU&U|XXqU|RxJzVycM=A=LtF*D6s>v}Ks6>9?`bC~g z_M5V(ztQ)_;w~O-`+21(2J5%AVcMw=+L0gW*voEpV9{Q0@zrjv$e;I*`CeziqN~+s zTbth_pY!sMotm~R`VpRWnp@AKt|=>?aarKCu5sR?b;u{|5q;u}_$<{ozq_Jd^5|Q2 zm_!WLfB7475C4OFMpZEcnme&5SJd~s;5qVxMdfZbxPtF~x^r2&0)7umc`_vo7QH5^ zkn80W^#Xohz2?A4R~O5MPFjaAK!N>1~rLf1yq3Gs3pWKH$8Pa@w^ zIa_%-*57YCSwV3L7|)fJ$vPp3uai3@+8+bI0b{Fk(&1R2zmaHN9f|ApZR?^>2f@#G zj>;B>^T?CYVt*$Z@qtGSs=|&4Uj{%WTo;){*U^~C&; zmpM0F#G(wt(^0XLG|Q4<3oV z-Tl1a4~vf3j91NY{nszVYO1{m;=N79v?|cvDkqGWg#F;YVUc44&OD0yk{dEh8Tv!d z3T=&5JeoNrA~;aPfl7{zUbl6`{I8Q`|Iu=woQJw^w>a_Ww_D&S%MkvRF86oGF_F*j zyejOqu>)N)XKUSd!2BpT{Gee8{rhM8s~TJQYYaL_oW*{4*?Z&RVSD8NK97<6;o?Az zzgsWLTj6_xR4o`A>e(&fSjl8`r;xH-|@g34PkdyRpCS-etkpgr9-LxpfnuKd8<*m$Py@ z{3ty7Z(Z67e^Gx&f7dCPKTfweJ!tQd5^!jh5|3_{$_?p2Z&+S#KI(*g{Qb4MTfNYp zXOzZTE5{>`%pYE5(as!a#~$XK$g_zjmX z#6(=U2YwK?o%Qj21M=P5`cs!&f&LwKVEjSk@2jkgSga17UEK1w78Vy+v22tHN^!;(Xn~ng5DQn%fYU^U{GzPgO{%p#8lmuPy!G!@qQLuR>=fmxT7m zc3ZdLcpi_UULgO!UEeXh1bSN4PFtH~*2KuPwYMoBGLim?ec|VIYWYfU0 zDT6=Y7j&(cvswoFY>82Qc4`5a4xKD_cnH1ow*LdS;Oku4&@{oo+z~vY{!IAPJoqs= zI~R6C54KuVT0WAG`_cNNOW*~){Kg@llQ~>^8+E*2&X-LyCk`9!$c4X}TJDZm@Jywn zGfu6(4F91=Pn0CVyQ>-%qCUYtsVw_q=_T;YvlrBNJjurWiElkO=LqsgUKWjJU*saM z!sApX`1Avjnd*93TvC+HZlN^r;FMp}WzoJ|rg7pO@a?ityFIb#T*^;yD0IFAy><`p zOu;GmOC{F0*QI0p6D?yFoPfVoq=S&}SvI{iG2gl_iA%fpi=Fs*9K2agXSg|*OJU(w zEI;t%tJRUC-jQ7Tn7?J-d+_u!v!e3j2hsne!CY75KRMQi1gnH{sWEuy%%em$H6Ah9 z+Z4v7pn^#Y95TQYRP;mL_i!n_xoYk8EH;VMa3!j@;rMPwnF<(h{9SQU#j=bYB$z_JrS=HDj4;74&`S0Q-C|E>-H-vfJMx|4Z@OX&HC;FKs+r zs@4R)5WGHmBmCeVL@3Q?w%~pzXX{>Lx`;VrT?}wl3t-KAC#WAmrc8yJ?v}Q;$nKC4+Wsh2P?1 z(>sc>y6^+rx#npD{215hY>Iz1i%Y%A$x4;7$p33mt5lfIrBA-9+BXz%yrzGy9h!>y z5K%2U3C9a~_r^Sd_KuW1mts9`4n9SGpMlph^%zy_kP6ZA!}Xmg!9BRDheMLK^F5{8 z;NO||v%KdkhvxY*W_Y{-kBW1j5DWiFm7X&TbDkkTQ09DTa|eeWHa{M=evbH&)t(0p znmJ^5%qXeh1@hCZYkL}+;FqaC>4`=u_%=VlTeqG=H=le@yIA(``yrpzScl{Jx}W$~ z1pZNF`M2Rq4mkvJt_5H|bsRH#-3foqgYjA6`FU)r+*5F!$~pAdvEYm5RW`Y+8HWzT zpYz}$@zXw-zn+d>v)oEJbTOdU`x)XhLN~2-oB}^jPoJ1kQ9Lj4DrFCL7IDb$?NUuY zJf9DxHB2oEIiyJoZlxV%6F12yJo5^Nb}*0I^x*j!=i;I??jnaSopsGK!SmBy9XIH8 zmP6gM1HSCR^Vk?Iv+!;jhjyQR|JQXq|1Ul2bCXkXeJw4vw@0%PTqfH*kO)6mk3+(W zv4{t-*N=Rgfc6JHPg{ij*PnIl{>tA;_{9WU=}X@eGsFTKT|RDu1N z*>tI6QSo`0S}VuGq*y%8QAdDW+R`wW95{Y1q+ z4Z$C1_H&y{i46Kyu5tgvPxy5_9`B-_!l3`Fy`5gE|1PZmI|?uOKEA}DZEBpO5!g?- zY}u)7e}O^AG#J|#jf1{1lyliU^}olz6uJJ_emOp^ul4#_RmPyi7Okjb*zZIq2ed6M zVNh!9yny`JA<(Ersw-}bRk7xqgT<>5nZEt26dZByR5=~tuOU#abh)t?7t1p=~shZV4>`p_LxCGon4QfME%w5 z309f6|2w|_|2>DwSE<($qZric#0}`L@dY2Qb8SZ1-TI$}O0~WeuNKsm6~dtL{1xvE zP}YCCYN22Yg9H~EJRMOM)CVon^<|LG&*@fg@O=n_pa!Ln>lx%`{wne&%GZ;AcdDcQ zO`W>)NhmKrwKUKj{f#}oMN*dZ*U!Fme?x%Qa596gypOmP*yBrn@&0+!PcrDt(NE@K zC~ujkKHL_~AX^#Bj#1QKk?_x<5|qcbFL&UiTmSZ@`E91_W(qSY>k?PJ1K(qityYOM zKso!JwdiM*wcRJj4gSRZDZbRrLVcmtRh%xA$F?_Ewa;Mrb`C`foP#9s|La(~?5f=k z4#~_;^lr!huYOm=KP@46AGn@$ZUzsppDW&27tEpWf7xV-`LN+H^uOEhIusMG=mWnZ zR$$!YhxT;D>qR$zJ^rtil&Vg?uqTHB=^61n_=~dV9=58qr;9crGA3N^|9pRAZ2LD$ z_4*fMeSKlgJ6^9GgYxE2lUH%U`rTbD&`X^`73x71bFhAo8@Kd+sS1Nc%R*MzVtvmm zX_&c2ia}%BtI3RBt$ojN!vIkRVQmdkX?ZM@%LLYhWZ)wnq$zNf%?Op z>Bh2M3>w?Nm%6-9?;_|s!Zkw`GSKg4j5iXPKyP|v(R>x}qZgY$>Ve)f7|+=Di|6kwbAx{P zboK6qs_@sDnS>C|7Mh(H>-z#ezSY}eXQRH znR2}^QC7)L%DX!a`O&E-R<8ZvM-PuJ>5j*GzGo1|vhd&z!*U%b)HG!?`?;ldzG%Wr#hPy0>zr!P@a}|Ma&!JO}O1zoMoRM9e;M2UKn9LJtKcn7s`EN zW-(>4lL-5jhqO7^^S+yPizq%_2kAiHFXZ9 z+Be;|wP90#+QA3;GvIgsB}C`G6`QtI$`n{=b4co8=4M?>tgnP6_f^9$Utzj?Y&^y{ z)<3`6_g5ot2;Uz$D6h2v^XXbt`kdJ)E3ME-jhw}%C81uoQhV_|2WHddH?#jczsB|# zVY1?N;9L#`{$9J^9M@~As>Rk0O^k2Lko`-Zfkk>xMhmJkJvE%)JH~QzX3z_2mo(^K0qrW=d&T{BHq;pY4KF|TX Mb`((~1_qNi0NeBt+W-In literal 0 HcmV?d00001 From 4e10b0f35f6068f592701f97c3e8ec31278c8b0a Mon Sep 17 00:00:00 2001 From: Kelly Kearney Date: Mon, 12 Feb 2024 16:33:10 -0800 Subject: [PATCH 4/6] Completed update to analysis/cold_pool_for_regional_models.Rmd --- R/interpolate_variable.R | 2 +- analysis/cold_pool_for_regional_models.Rmd | 181 +++-- analysis/cold_pool_for_regional_models.html | 720 ++++++++++++++++++++ analysis/roms_level3_coldpool_part2.R | 136 ---- data/B10K_lte200m_polygon.dbf | Bin 85 -> 0 bytes data/B10K_lte200m_polygon.shp | Bin 27204 -> 0 bytes data/B10K_lte200m_polygon.shx | Bin 108 -> 0 bytes man/interpolate_variable.Rd | 5 +- 8 files changed, 862 insertions(+), 182 deletions(-) create mode 100644 analysis/cold_pool_for_regional_models.html delete mode 100644 analysis/roms_level3_coldpool_part2.R delete mode 100644 data/B10K_lte200m_polygon.dbf delete mode 100644 data/B10K_lte200m_polygon.shp delete mode 100644 data/B10K_lte200m_polygon.shx diff --git a/R/interpolate_variable.R b/R/interpolate_variable.R index ab8a8d0..706a31f 100644 --- a/R/interpolate_variable.R +++ b/R/interpolate_variable.R @@ -19,7 +19,7 @@ #' specific bounds, or a 4-element named vector with values xmin, xmax, ymin, #' and ymax. #' @param return.raster Should a raster be returned? -#' @param outputfilefun Function of the form f(mask, year, method, var) that returns a string holding the full path name where the output raster should be stored +#' @param outputfilefun Function of the form f(mask, year, method, var) that returns a string holding the full path name where the output raster should be stored. If not included, default is [coldpool_dir]/output/raster/[mask]/[var]/[mask]_[method]_[year]_[var].tif #' @return Returns a data frame containing cold pool areas estimated by interpolation methods, for cutoffs at zero degrees and two degrees C. If argument write.to.file = TRUE, also writes a GeoTIFF raster to output directory. #' @export diff --git a/analysis/cold_pool_for_regional_models.Rmd b/analysis/cold_pool_for_regional_models.Rmd index 538fa07..4466f0c 100644 --- a/analysis/cold_pool_for_regional_models.Rmd +++ b/analysis/cold_pool_for_regional_models.Rmd @@ -1,7 +1,7 @@ --- title: "Bottom temperature and Cold Pool Index for regional models" author: -- affiliation: +- affiliation: NOAA Alaska Fisheries Science Center description: email: Kelly.Kearney@noaa.gov name: Kelly Kearney @@ -15,41 +15,66 @@ addr: # Introduction -This document provides code and documentation related to regional model-based versions of the cold pool indices and interpolated bottom temperature rasters. +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. -TODO: more detail once finished +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 following files must be in place prior to running this example: + +- survey-replicated data file, .csv format. See "Survey-replicated data", below. +- cold pool indices file, .nc format. See "Annual indices", below. ## Setup -```{r setup, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} +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 sets output figure resolution and the map projection to use for our calculations. + +```{r setup_nodevel, eval=FALSE} library(coldpool) +``` +```{r setup_devel, include=FALSE} +devtools::load_all() +``` +```{r setup, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} 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 ``` -## Survey-replicated data - -We can perform a rough comparison of the model output and the survey data by simply extracting the model bottom temperature field on a specific date of each year. -TODO: describe survey-replicated data, and how I build it. +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](https://beringnpz.github.io/roms-bering-sea/B10K-dataset-docs/) for more information). These details should be modified as needed for other users, computing locations, and/or models. ```{r sim_paths, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} -# Path to mox-hyak /gscratch/bumblereem remote mount -# Note: this is computer-specific, and will need to be changed by any user -# trying to run this. -moxdir = "~/Documents/mox_bumblereem" +# 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 -# Survey-replicated model data, in .csv table format +# Model's Level 3 folder path +lev3path <- file.path(moxdir, "roms_for_public", simname, "Level3") +``` + +## Survey-replicated data -#simname = "B10K-K20_CORECFS" -simname = "B10K-K20nobio_CORECFS_daily" +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. -srepcsv_path <- file.path(moxdir, "roms_for_public", simname, "Level3", sprintf("survey_replicates_%s.csv", simname)) +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. +```{r csv_path, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} +srepcsv_path <- file.path(lev3path, sprintf("survey_replicates_%s.csv", simname)) ``` ## Masks @@ -64,7 +89,7 @@ akSEBS <- akgfmaps::get_base_layers(select.region = "bs.south", set.crs = proj_c # Read polygon defining the Bering10K shelf (i.e. <200m depth) -modelshelfpolyfile = here::here("data", "B10K_lte200m_polygon.shp") +modelshelfpolyfile = here::here("inst", "extdata", "B10K_lte200m_polygon.shp") b10k_lte_200 <- sf::st_read(modelshelfpolyfile, quiet = TRUE) |> st_set_crs("+proj=longlat") |> @@ -94,16 +119,15 @@ maskname <- c("SEBS-noNW", "SEBS", "SEBS-noNW-modelshelf", "SEBS-modelshelf") ## Cold pool index netCDF file -TODO: Describe the coldpool.nc file +Model-based indices will be saved to a prexisting netCDF file. ```{r netcdfpath, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} -cpoolnc_path <- file.path(moxdir, "roms_for_public", simname, "Level3", - sprintf("%s_coldpool.nc", simname)) +cpoolnc_path <- file.path(lev3path, sprintf("%s_coldpool.nc", simname)) ``` # Generate bottom temperature rasters -These rasters will form the base for many of the following calculations. +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. ```{r raster_calculations, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE, results='hide'} @@ -130,7 +154,9 @@ for (ii in 1:length(year_vec)) { # File-naming convention - outfilefun = function(mask, year, method, variable) file.path(outpath, sprintf("%s_%s_%s_%d.tif", maskname[im], method, variable, year)) + 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])) @@ -175,34 +201,101 @@ for (ii in 1:length(year_vec)) { } ``` -# Calculate cold pool and mean bottom temperature annual indices +# Annual indices: calculate and save cold pool and mean bottom temperature metrics to file 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. The resulting metrics are saved to an existing netCDF file that also includes corresponding metrics based on the model values on July 1 of each simulation year (those metrics use the same masks as in this example but do not involve any interpolation; instead, metrics are averaged across the model's native grid cells.) -```{r} - -system(paste('ncdump -h', cpoolnc_path)) - -# for(i in 1:length(bottom_temp_files)) { -# bt_raster <- terra::rast(bottom_temp_files[i]) -# bt_df$YEAR[i] <- as.numeric(gsub("[^0-9.-]", "", names(bt_raster))) # Extract year -# bt_df$AREA_LTE2_KM2[i] <- bt_raster |> -# cpa_from_raster(raster_units = "m", temperature_threshold = 2) -# bt_df$AREA_LTE1_KM2[i] <- bt_raster |> -# cpa_from_raster(raster_units = "m", temperature_threshold = 1) -# bt_df$AREA_LTE0_KM2[i] <- bt_raster |> -# cpa_from_raster(raster_units = "m", temperature_threshold = 0) -# bt_df$AREA_LTEMINUS1_KM2[i] <- bt_raster |> -# cpa_from_raster(raster_units = "m", temperature_threshold = -1) -# bt_df$MEAN_GEAR_TEMPERATURE[i] <- mean(terra::values(bt_raster), na.rm = TRUE) -# lt100_temp <- terra::mask(bt_raster, -# lt100_strata, -# touches = FALSE) -# bt_df$MEAN_BT_LT100M[i] <- mean(terra::values(lt100_temp), na.rm = TRUE) -``` +We include some quick checks in this section to make sure we don't overwrite data already in the file, and to flag if our most recent calculations diverge from any previously-saved values. + +```{r annualindices, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} + +# 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")) + + } } +``` diff --git a/analysis/cold_pool_for_regional_models.html b/analysis/cold_pool_for_regional_models.html new file mode 100644 index 0000000..ab98765 --- /dev/null +++ b/analysis/cold_pool_for_regional_models.html @@ -0,0 +1,720 @@ + + + + + + + + + + + + + +Bottom temperature and Cold Pool Index for regional models + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

+ + + + + + + +
+

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 following files must be in place prior to running this +example:

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

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 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_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")
+
+
+

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))
+
+
+

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")
+
+
+

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))
+
+
+
+

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)
+    }
+  }
+}
+
+
+

Annual indices: calculate and save cold pool and mean bottom +temperature metrics to file

+

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.

+

The resulting metrics are saved to an existing netCDF file that also +includes corresponding metrics based on the model values on July 1 of +each simulation year (those metrics use the same masks as in this +example but do not involve any interpolation; instead, metrics are +averaged across the model’s native grid cells.)

+

We include some quick checks in this section to make sure we don’t +overwrite data already in the file, and to flag if our most recent +calculations diverge from any previously-saved values.

+
# 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"))
+    
+  }
+}
+
+ + + + +
+ + + + + + + + + + + + + + + diff --git a/analysis/roms_level3_coldpool_part2.R b/analysis/roms_level3_coldpool_part2.R deleted file mode 100644 index 912e773..0000000 --- a/analysis/roms_level3_coldpool_part2.R +++ /dev/null @@ -1,136 +0,0 @@ -# roms_level3_coldpool_part2.R -# -# This script is part of the post-processing suite for the Bering10K ROMS -# hindcast simulation (see https://github.com/beringnpz/Bering10KPostprocessing). -# It calculates the survey-based and survey-replicated cold pool indices and -# adds those to the existing cold pool index file (created by -# roms_level3_coldpool_part1.m). -# -# Note: Because it interacts with some remotely-stored netCDF files (located on the -# UW mox-hyak HPC system and accessed via a very computer-specific path to a -# remote mount), this example will not run out of the box. - -library(coldpool) -library(ncdf4) -library(lubridate) - -#---------------- -# User setup -#---------------- - -# Cold pool index file: on UW mox-hyak at -# /gscratch/bumblereem/roms_for_public/B10K-K20_CORECFS/Level3/B10K-K20_CORECFS_coldpool.nc -# Path below is the remote mount path via K. Kearney's computer. -# I opt not to include a copy of this file in this repository to avoid problems -# associated with regularly-updated binary files in git repositories. - -moxdir = "~/Documents/mox_bumblereem/" # Path to mox-hyak /gscratch/bumblereem remote mount -hindcastcpoolfile <- file.path(moxdir, "roms_for_public", "B10K-K20_CORECFS", "Level3", "B10K-K20_CORECFS_coldpool.nc") - -# Survey-replicates data file. File format is identical to data/index_hauls_temperature_data.csv -# but with an extra column added for the model-based bottom temp values. This -# file is a copy of -#/gscratch/bumblereem/roms_for_public/B10K-K20_CORECFS/Level3/survey_replicates_B10K-K20_CORECFS.csv. - -srepdatafile <- here::here("data", "survey_replicates_B10K-K20_CORECFS.csv") - -# This shapefile indicates the portion of the Bering10K ROMS model shelf that is -# less than 200m (model shelf is narrower than real-world shelf due to -# numerically-necessary bathymetric smoothing). - -modelshelfpolyfile = here::here("data", "B10K_lte200m_polygon.shp") - -#---------------- -# Cold pool index -# calculation for -# model and obs -#---------------- - -ncin <- nc_open(hindcastcpoolfile) - -# Read existing data and parse years we need to calculate - -time <- ncvar_get(ncin,"time") -tdate <- lubridate::ymd("1900-01-01") + lubridate::days(time) - -meanbtmp <- ncvar_get(ncin,"average_bottom_temp") -cpindex <- ncvar_get(ncin, "cold_pool_index") -thresh <- ncvar_get(ncin, "threshold") - -nc_close(ncin) - -ismiss <- is.na(meanbtmp[1,2,]) | meanbtmp[1,2,] > 1e36 # Don't know how to query fill value in R... - -select_years <- year(tdate[ismiss]) -select_years <- select_years[select_years >= 1982 & select_years != 2020] -# select_years <- c(2005, 2017, 2018, 2019) - -# User-set variables - -proj_crs <- coldpool:::ebs_proj_crs -cell_resolution <- 5000 - -# Masks: w/ and w/o NW (82+90), w/ and w/o narrowed-to-model-200m-isobath - -akSEBS <- akgfmaps::get_base_layers(select.region = "bs.south", set.crs = proj_crs) - -b10_lte_200 <- sf::st_read(modelshelfpolyfile, quiet = TRUE) |> - st_set_crs("+proj=longlat") |> - sf::st_transform(crs = sf::st_crs(proj_crs)) - -akSEBSnonw <- akSEBS$survey.strata |> - dplyr::filter(Stratum <= 62) |> - dplyr::group_by(SURVEY) |> - dplyr::summarise() - -mymask <- list(akSEBSnonw, - akSEBS$survey.area, - st_intersection(akSEBSnonw, b10_lte_200), - st_intersection(akSEBS$survey.area, b10_lte_200)) - -# Calculate indices using survey and survey-replicated data - -cpname <- gsub("-", "m", sprintf("AREA_LTE_%.1f", thresh)) - -yrmask <- year(tdate) %in% select_years -for (ii in 1:length(mymask)) { - - totarea <- as.numeric(st_area(mymask[[ii]]))/1e6 # m^2 to km^2 - - svy <- calcindices_temp(srepdatafile, mymask[[ii]], - select_years = select_years, - bottomvar="gear_temperature", - surfacevar=NULL, - threshold=thresh) - - srep <- calcindices_temp(srepdatafile, mymask[[ii]], - select_years = select_years, - bottomvar="model_bottom_temp", - surfacevar=NULL, - threshold=thresh) - - meanbtmp[ii,2,yrmask] <- svy$MEAN_GEAR_TEMPERATURE - meanbtmp[ii,3,yrmask] <- srep$MEAN_GEAR_TEMPERATURE - - cpindex[,ii,2,yrmask] <- aperm(data.matrix(svy[,cpname]), c(2,1))/totarea - cpindex[,ii,3,yrmask] <- aperm(data.matrix(srep[,cpname]), c(2,1))/totarea -} - -# Write to file - -ncout <- nc_open(hindcastcpoolfile, write=TRUE) - -ncvar_put(ncout, "average_bottom_temp", meanbtmp) -ncvar_put(ncout, "cold_pool_index", cpindex) - -hisstr <- ncatt_get(ncout, 0, "history") -if (length(select_years) == 1) { - newhis <- paste0(date(), ": Survey and survey-replicated data added for year ", select_years) -} else { - newhis <- paste0(date(), ": Survey and survey-replicated data added for year ", min(select_years), "-", max(select_years)) -} -hisstr <- paste(newhis, hisstr$value, sep="\n") - -ncatt_put(ncout, 0, "history", hisstr) - -nc_close(ncout) diff --git a/data/B10K_lte200m_polygon.dbf b/data/B10K_lte200m_polygon.dbf deleted file mode 100644 index 6b6f261aa0044d94e62a1f037fdfe4f7e4430df3..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 85 zcmZRs;$mlHU|?`$5C)Q%z$Y;&H3uT>45Eb4l<+Dzr50u8r5hS}$LEx!#v2(JpgM^Xmf+U-y_m~m1ZFE(DBc{k8Fso-qxWU(n>Q^lr^ zO$(boHo!*0W{J%n8v~mMHXm$$*rKqd8GB}&TTKIB|GD^EG6cA@SyCnn=}*UA9r0O1 z0KV7F8bcg3$oczt3ma*V43^<1q-!^Pw$8);|NZ@c)5_zU#Q!t?j;jUaepLc|Uw>;* zU4RCxZyQ5;J=e|J^^uzoSDF`!`x$YZlu@0K1BAkwFm}cnq`P0tWVinh^EK)fmB;L(!Px7+VfT57aPNp%sm}o#q%swU?U4R*Lqhog z(vR8`Pi$Qzz>SR;zPe%m>uJp&&mjH3eP1S17?h6RJN;29IwwJdvX6x;B1tsR@{>EQ zjC9wHzR*KRpRQ<}oag!f&M&>7HoKvQ1|x6pq%KdAKq)-o;-w}Us9pN-{?`NvF34ST z6+&7n^F{vj00|bSBn;E5X>jnhRmH*}3GTVM(B+=ffXsENX=#WA?{4S6&woaP#Xik# zbMHvt>JVC|@fh{TubeX7MS?4~?fex*G-x$>I??~01O`9&`OT2_yjsfH-AsZ|!bbTI zw`frC_*{L(D-xs+s$6AWqe1w9WaO0=5)`R?UiUAT25awpbx*29+W3r5@)a6zTDxSp zRgqwh&-3$Er0Zox2W%gbU}N;ux|)kLDB^EPm3u@2m6_RKSET>$*k-o&774ztJn71b z^z)CWMps=S!L8({Z?9cIc`TSdu0uNOwN1Jy(xo5XzIhx^g2?vwUq{c;;9{cdF3$uK z{OD?WaUq-rLs{R7)+drcBG3Dl9@34UJil^BkU)U!LK_UEL5)BH_g$pr6B7kek+v_~ z*O~55f~4(@G6#^p>%8j9StbcSf7QPH_aqIDeAIUKaU{Wu(EH8zj?;kGvNo9aFbR$+ zR*ktMEnT?yWih&9VtiNA- zl#8@(vY)bt{D64-rBxt&_~8(!eBmO-&)vkj6}5BmFj6XWxYpB1os&(Ds?rpjd;= zHHGw}k+_9*#x#($iWWB;CBlbel~#cCiH7|=tQjJx9EKcG0~+KXin?5m^s@>lVO^yC zDbXK3tVH>*-CW{`beCAK&fFRjsBioic^YXD-+fnnBMD@=gL2DpJd4hm*cutMmon9k z&qxP!JsRp)B|(ZeSzpiu^^saaT+k%}`KWcT&VOiHR;{=R+TZUS`%nT6{*;94lpACJ z(C^0zkyd%udZyEY1T7n`PyDe!`xdYs@J8BuTS3$oOB#&RfB%U(g7VoKOmjp!*lmMw zjtAQR$<}b#{93iNIG)IF3X4(^W0A;D7B z=Iv#kG%zCf>;KCpfren*@)FX!k}OsU7m^^_kR6bQ{ll&H+h!o$e|NKCH_}2kcDUUx zBSABL$#nyR2CIF{9LnktFH}kOab7g|TpjJa9ce`}Ja0$(c2vZuUI&g}Jim6mH{##v zSla{rB)GnCp#3P)%Hun4T8@(-FYcDfrxP@o7A7TEd?&%aJ28{0{)i{*>s=b=NN_|o zI{fY_8bljy=-7gEO!4)c?GZFM|GHn{%q+@tEXkxO2II%!N!vFwB;fh9f`c~}?ct2X z>-Q_k@PKpO!}vrRaBbP=D}(ggW%B`%G#coH&Hl0HBtwE|o=DLp8vIZSXIuTj@rL+X z&5?fZ?paD0A%TBr0Ie<)@lWA`*y(R1DCsI^8|KkqZZG+8;a8;NkCzl8O|3aD^57S) z|4C({WIhep+g4^j;UGiSx-YAvkp5|!u<%Tf41LYbJc4&&Lr$L1MHG8|JX4v9co@=_DetP&Z-$d{Z%%FtfUChBq_9g`QG zaP7kRf+>mdnFB#1HGZF%AtfNEg!r_39DdE=NZ$ z?IDBjD2LaZ7W`gQ^Xr5@8I0BQpBN(@_E7ux4>dB}jasUI`U>s6s=4O392r(c`DUmi zEkX3Nn^{W+HCGp}TeUQh4JhX1k|2W~f9C~xq=Ujb%UAM};q<==_Z*Bj$72f&ITuOr zzr`P4GtwEEW*YRkEPeX(mjnl-sh=|6puadcG6v=me~2=E!bl$tN#yNE`}yzq_22!E z{ty*XX~+1nUnF&4H5vAWCT%HvkMTm?+4kisGVm1C1gu9og{<8367B21zn7P*XjM~p zhe)j$rze}q(3q?(#-ZU3J%1W=ZfKEVEl<_f^E&SEAk!w3b0-;;nqRp-)^>+hgM$v@ z%4GQO@BjPvH-XKs^!B&|pWS0%qyDtT=1N@j-NCXJ9*3+W!@pZ2wF6r25GSSHtc>>Y z-{0@y+~S-ahs!FSAI<(V5~%r6?$d9J!cj^|Ji@c|i3 zzg*n$_=Y<)2b?MpDJH|snN-W6n>b#wnp)3uGRS`LG*7wb4vXExm8s=qh@$c})fKtJ z(95!e4h>{*n_I)RaNiw@w(r*6_?8S1d5x+} zR!LC6PUF}@%_8b2-9%DkEd^+oUq00R=MKGFK7M^7O97t7;=mQF=wQCrS8_>)0`=aF zf+4HvU=ebRo~ufM@S}xp2L~Ia*r3CCkOF*H zW#8H;(?K@#poWwQ1$MQKQa5j>!^UvlkE%u#xGEr%QmuybTON8-L8ic+sy44QZ93#J zuN|r;P@qIAXqd5w4z-11$rU)Br2*eBo&9u(pmayx$MLSY2wsWbM~CsDi(BP)Q=lb* zlxTC14q3*QwE{X6IMP=uCPkqG=yNGdYEpm|b4Bix4IM7L87I`MQ{cHj*TbB{bdXT7 z-+yTv1sDeoR3CPu!y-FZ&{2f~IReMkW!&jdx@%wlCv6HCq#if!^PofP^_TlHaJ<=g zk&f$LbU4Vm^-j_d=kt=7^Y)>Gq2O_P8m`|kUHG~FaXO?`mQIP-QQ-97Oa*>FI#@py z@H$VUfOFu9tGhhuz~O#jvYSBx2QJaotG($^@(t%aRUbCi2z9?-aP(!p%Nrjq=D#Gpd}WfLrc# zygrVXoVKfZ9VZps|5z2vqW;_-AH7qsk_z6lDZ!&u$+{H_sGYGPE_#mlz8W<-abe|kA*1gNkrfwU(`pAJ{Dm2=ALlt8{R@SKHrcN>G8!!JpZI_Dl@i9dHll`)5%9DoO+O;WVc@B146Eb?NJS zHR%vySXd{gNCmsz&8o9I=}@gNZtS2+h4eDQG^ILb1H#NyG>4 z;{5tOR0u6kJWCCt!+G^&$!Q}heC5ok&OA+rl{32tCi`$bzgck^5p{n%jggFEQ{i5N zqjFmU9lF9L&+PT2Lcuq;O-G{X;5;t*Ys!lX(apQ&J#jsUz4T5ApFsKVpH>LR^>Fup zD&HPVg(EK1C(Ds259TnV?-Uijo!os*IEoHRh3SC{XQ*JJq4h;R2JwrexK}cU3jRC& z&fiO={lxc5Gb_i`MrbAnie(s6mRJi?p?NnwF9cIhE ztvT#Lg^G+j?~dQ4!%y25O(Sf3WDdh+>qIxzCznT6x`{$kJ{ONXa^=DU*7UhWXPe`T&AHt) zN(C)B4ufO(y_JP*mJ?8+*u_>>2*3BstI4P6q5hw`H)!DZW@%<@vt3m9+<2{;BFX^c zqwW3%=wAa@gRg!QWqJrwC_QKmH40v1c++i2ut={_7>)AvGD7^Y7xf<~?Z&IW!-iiT? z!cxPwwN$VnK2Z*KV1Vi~o{o#6R7m`|`j)E;12h$6I<3}Fq2Zgc!9ON%owkb!{tU2<>k^vYiuf;R<6#uY0E->%=Z1G+yde~8??29f{CDc-HBp}I zVZHnrKL*SntXkNN`Z~Ph*XE)C2CP40HY0`lBpiR(!U|%*NhVjI5b9f+%lnE-2m{`P zD~By3e#t(vR9cJk<;;s;BcZ)WizUT(o<#jyrCc#UdwL<`yX|r?1Gq;Iu2pkDf6&{> z?;Og2z+DV}X|&fp4sO;FXK}o#o`fd!mxuQ}Ci4LH|(m<`z7a z#{iAXi5uOS=ns;o6R+H6z}e5;=NHib9ZKa#;_fm)VAyu>BjUqP5B~b_`wS3W?V3}L z_>re5u`}oa1A?_Ut1h2JJWATV`{ZK=Y&d1O-!ByLIUsr{sEh&eGtx?=NTe^`&n&ED zz{(*9llP~ophOSJ@2+P6hhOuv!(mi--6|(N+QfjuO>%uZ&SLx?b+G=|!hmgDzpYjw zJ{I+kX*AX`ApA3RDFX4+NYC#;{0jzTh?dMPBv66l-rr2p6SUuME2Cq{Xy1hEn!Lpf zNXptbY;Fy6%{LvCB$omSDKU`qIo}p#YvWK`{3&U6Ki4541A#a-c6y^DL zAj>WW@q_1nz{oQy5NEa*_MBz_E%-;X?h7g`KDS@6MSHU-&@_ywLjONon(_wi^*Ob& zt)ZF<8d9YJ#%u;0ppFFc)KMX*-1?e31N~Dyqg%0(3S0w`>wBCTAf)+N&iN(g3-^xv zHnU?u3h`1{Y7OQmH%oLMQyGx^Is1DpeouVK9F;OdeEX>DKHY$Lw`bGHAyc$JZvtZ# z&L^$E*8c_i|5>NAr&VjIVE;iUT22@3i|gez^Li>Ad;LoBv^oQR423ngH)4LHpljWO zc*JD~lc$=gP#!oH9IMIzX>Iuk&o@-)@pxJCL7oAfcMh)#@1R2W)ar-PN(`uXqA`5B zsUVZRrT4WA1M}%F!>}GIyjc}~Cm8Xr@X}>3?*S_K`VmMDSCAL@`o}t1f{l(aqGZ-%?RwYvBFkj#oBEO%+ ze5B!6uHj!QggCu7W&WUpPW*j;Ee;DXiTs=4H%bRJPTw1wI4z*-_=AMtK{|Z8U3hgB zmj!TLb1qnk@yu>lap(fC1@zL&#W;E~elKgp^zmCj;ahX%IgF2+W7t;Zf)?O089Va* z6&>q=9KFfH7U1>5G>P6y2YY?)GQ-&3t*T|b2sf`K<}M@q<1SUAiPkRa1$8tXE5N(lV!}0I=yPb4r06z zJ=v1FM1_wJ&Doad4-un}zIiTCA*!2qw;cM5ob>l8!+FeK4``}8w~JbuJK}{{V%o_#8n`E%T}0OU5);6E8Lb&eZ+vhP1)nw zeVC6$GGsqJLws7hY^U)N{f{m8h+e~h_fOqMD$w8fbf=`pni&vFuFEm)LOj@h_z<&$ z0S6mpPd9hsdh^yE9{j)nIxQf{tR44nO~WcSgADkX(9~Jeiuio0;#tiI;zN@xd3OsH zE_2-x+&;~K^o!?K-f2R-8~-4cJj(#G1=pN7;)`cR=FrqXj9)wNs)r(;ST;KyFfOlV6pA@`Nv@)+W& zN{FR&umTgl)NQHNdxHM?{q0nlDieh7T>f(I5#}q)r`&(;U;=;f0Idb_ne#_{u<~vu zyp>lr6S|N2-dpc&Hoychzi^3gmkKv+as$H-FyXuPZt_{g`?SE9hx5ox2r{xcS%&d| z$2y^_%7zJlU-OEOI2@p8yyETVspKGKt<;r_#wlIFgW2`BVwI09U#@SrH8;RqiSyf&RR zZaaehvQ}@Kl?W4L>}z-xG5;yxe)vyX2>BUXt>|QUTKV00zgn`h^FJ}&-y}UiV&P$mIUC$Gh z8cnF+aZ&k#2HIQPHz_t3U_M~GL2(7zn>W>I?7&_sG{19RYSCbVo0&3CH0G113MWG2 zx1;<>2$M zDsWSrP9|di5+{?wQU(0}V#~-reI{^jvw1opONA$6=DLvvOt|vxx>1ob6`tp)-qJE* z!s$X8+kWg{M16L25bf)zo!KSRQw?4+nL z;jFLy%9#l>jEc@cNh&ybr|eX4V?v;_w7xp_?=C3~t70)>X=B)aGwgr$-`AJ2UQ8HU z;j6od>mfzIefr#w36ghilh$DWpB@dyN+;2NTk`8GaXq}YFS;s2nb3GsaqKnrzn-pe zSK|y5hHEYj1mXI>PaUXhie`ez;je2Vu>Y(V%s3@5!RU0dU=;Ru@|f8>k&67C4e2>2 zgY_E^X6DIECe|h6Po3F-_C0i4aOoNouD>t*9f|ee+=ia6JOo)y*FTTAA?ZJ7Iw#q{$6f*d-O_Xqd{HM6?WhoPkZ66f& zt5PB5^f#t+ITIfJoajEOL4}rEO5ew-nBb6Ww($t=xBgD`^(56YVfX59=|54ws=e*> z4X>DR#5r@k1M4%FwvXmFwld)(ulPMTZS=prvdYwVsQ()vY=Hewq-m9`{=kIisiH!! zbdi6}EfN>{@q6Y-*j2Qju*VJ{@RbR@$&cUp?7;mqf7fNo5EEX_4V?bE73*cIePq)= zF@aTiV7DX6v-fVyo8cZNh^@5ZZ2gbCucg&^zh%NMEeZ8Gld;P?Ne`Qdc8B$EED|MScwC=u2EHu%18+ zKiD+Pgpf_IKjtq`KyH&ITM_vwb`YIBJ57P?5}{1L5ya1_3I7da6zI}l-|^-P6QnLw z(ic8aK;wAuwB{fa?2dfP=;sACHp%^B%`fyu{MD$AqHLaO)zhzkcUn@>BL; z-t8p?q|DZQ5W@a5a)MEM&nPhP=9N(v+VAEk;>y!l-;E9Quu;SQKgdV+=Uu0OF6S|x zx*OP^PyY)wg974Wo>oWmm_Xb+C$t!k>l5Nq(pfK)vIuq!htvjB@ zP~gqwrlTiPnDFM|ulKJKC@|8cTeA5A6DHjpIy+M-;L+G0(sd5U<58>Vyn^zkSJa*k zXTn!c@q|}ckAB+uDONd%3FkGQ8a%s6fvnpPJp23+FFh}E5HhhoEVe5hmKOyaBXy_`u)leeq|1m41*#=A6C<(zZG|_nv_lkF z$(XRnM@3TtPR#sT=m!4PZUt`Y%3HvUN`kfSLKlOOcFO2sh zyb(#-TTq_dkscoqKd&U_^5jWV;BOk~X)Y1*?T1RLi#*oP)A!gVm@#2blXqzYp6BSF zeIb0>lnD*@ZQE4&DG;-c((S(={W-bdOWZOU`eTCc_+WfbyuWer+b=SlcIylC+RKE# z;$u=qc%J1a;Bhex`;Ww}eU&mqhHw9#9=L@4c^c^I6J2E3=AuaU(PP4k`{aWecpi4% zSU|)d;|J5sHD=-s87gyk{kV+d^_F^HH-1Hi0QQrZ9PGa}xCTMtz9WH+?huPkBwL5 zC$N83>Ck~~m&wp6EWD)sKh^Uyz{J>DFLn4lY|nynj5hC(|>nMFq? zwEHfU&`yxy@sx|_k{!l>*>SyECK))d$I*TsVS;ZwwQJm!3^}$#MZdAXOu2pq(UuHC zE%A+ivA=9kY2$4I88poHNAO!Qf$bGg8f-y^YrRkOC8SmmHg;m<%V}IGM7nkaCrxpU;q;|k8$v{XfKldd*MC2-7#CSEryZGoQ3ASI8-!_PNA?aPrm? z?dMGR`89LOcLf=8>*~G=K4rp(-sxZ5cz=X3cdqfp3nnP`ceOv7!uv$o!L1$`kLxdg z^>r8}f%x@)(fuzG|E=~I{_4Z^>fh@hMSPE%iWa%_9@mp1l3iVc@hz-+zZ%|OGy7t= z*R6>OtNq_oTJU~M5o|uO7yC=tTYPsf$NOl*b0?&*f55Ou!eS8#;yu!KFIA%dNp%<` z<&)s%gYjoWxZX31TR7#?QT}@Gi2>DUKN|)E>+pWhFs)H(bt4llX@`7pJx_wGL+2XC zaJ_PO$R}pQNf5@%ed9zk@+b3Zb8-*~E`Am+b!b3;iRs#->`j7aAD;RdV}H9^6`wL! z5(x6Db?-p_9^A35Jz$IXcjO1Fi=HuIy@t$8z7q*LjUxE+%9t=NRJMM*CDO{deex(T zDq&*ZTLKBDaw5fZksq4tI*${0Un-^jK`@l z{TfEgn5n}=z{LpPF;6Z|c@lTZvUzw2o_Ft7-0LtU@ zi#vZtF&_=&*jaXl1XAiRPWerty_H>jqIe4LKY4G9nfQ(V^)sFlcLMKo-Ijgh_80dH zm;6_Zdg1*=wdYn{%ZN{=?VZ3N!KLG>VkWCtaH@0j>*S-T&k-BrGdwJ~p!jg29p|44 zQ#d#&$N~qA;dDuN64;+%N7#z7AfSO2V2Ag?!Wx5rJ(pww*Ui=o$tdslq@wa*aTXY7 zSbT9L;(bfUL64pdENBhvE>A{zY7SaE=*Y4_(P6Q>%76qXu4FI#l4e2qjZx{uop?V{ zuvFh%o(1;Xr1_t1BjLRO@2CVN77P{bHyPC>fxesA#jY(Zkp5jqU%Q_KR@#hHAqp&r z^YsavGDCiL#?b1OST|eN9`l$VN!+kpo^s4sUVf(W-K5ta9OpYymuVEQPxjpLG6

zSs!JAAms4eM0+#+ZDIZOFn*sBT747cG1tDC7-G+YVCffX*U`R5A{7)R4za*sl46yE z`YN+Z>w1XuW%A1CT+$UclZ-15<10t*7gevyI{NwAgQRC$jD z3%u`-pYW8y-^Fi7d?um&To5yL5+%W(CU#zeISUq5r^d{=Nf4Og+}eoay|^Uw(0~K) zM{PVh?oVXF83p^wZFt{uZo6!?GRnhc;a|bJF(N!Xsc-(#gaw~3Rr@c0B7)?UVT=II z&tms^&%P(Zm5wzlY7Vf#be&*Pe+v;d{T-DJHD-Z!{qTp{mqe&JA2p=6j|IY|>*sEk z5+S|eHY}okBF@{dQrxHVckr*k7=_MT7$9(4Veo z5Bb?=71XkbAot{Ly@omqzI_R=<47jLVflCk&TTB%@pIthn{!04e1CrBd$g|`Z*(C$ zmPfdj{dbT%9N|wnFyC0*%s>SSil=i7>l+h zLdrw#c|lPY6r2^Vv2i4V>?q%ZSs@l|nGf?8q!Pi!W&F(-ew_bvigmRy-bY>iv$$q8 z3uGW~EP_CU%1=)?OINZ$;gN2^4j_WV%JSku4)ots8}IdL5kW!oY+u|0=2!JzJjwb* zpuH6?KR3sOiMoLq9StG`N$nnu{Db&EsJFXE5&7G)nSSmk=KptJzYShbgsf$rZE=WS zkFtkbC54EvK1^IEWt<7mB$~H9;v~X6>3n!L;@v&ZieJQ80tBnuhn@Yxgq3Ow#luSk z_$+hMh}n;Lnevw-kc$YW>0)Ptdzmn+pL(8qo&dXb?0+3b{2j8d7f)IyfNSM->IuYW zN6wD=^}h)aG#7co2Jw4(eOp}mUjmG~kS2~JzP}pb^OBw-z~51gHVce52gc5@lmFm- z{C8?@ei(1e?i9rcj}hQ#Z;*f`#=C9O?v=sQ1b8@M{wxsVT@RCE(eO|a?%##f&fjPx zfNK8L`E_@(KCriFlSl;tc$0@@sJF4+N4@0~Tta{|1(BrlH!*+imNp-}K>*Wu&a+MV zOo+cPIA(H%04t7{JuAWdM3i<>FD{n=dx|bk{?5dFi6(J$G!@5Za&6MjV1kBSud!(? z0V1PSEfiC+-V$@NG&GR_qn{20`Xu9i=a;_byGR1!#gm7&7nq>qBtOR)M}Uh8>2Geu zV*OJ3oVZOGu0K;|W8zuNCk}4Z2s%rE^J5k2Uoii=AX+tF7f1k?y;AATAsCOMV^{;& zUyS9Jtr(2?VJY9(WcHv{gb}QOqHF$eDnOtIVA(MpMz^o z*qmfStZi!TXEOqf7?VV3$bUukN6WVT1elk3=({Q$^PBv0rqmI}l(dyVo@o z^VQD7CV{@z1o*m{`{?#GChT7yuP=3!03=0`BfC)l;qzIH0eiIHjuQQeY^+xfc(q2l z6X2xMJ+lGKZ_j`F*7nka0A{grifPz?U!<=F{UiY7H0Ha(7J);z;|ZYq zdc^D#_7@pdmR-0=fbA_5-aPEzb1lN`;#C58N`I?(jQQ}MZ1qlq3&UX3Vr_n#8_pEtaj`^SXKxV`R0!*#( zx%UbERh6Xtbwd=&GsD*CEBf=-wY%HzqJ5c~%3Uf&e?N3~=!E8Z0toGY`#To%PoI+< z&mX4|pkeI8)v!7y*dBg)NHv=PVOO&>tT2CFPo6qcR6qdbDS`K@SnujrCLGjwg!*-< zAK(23?bDDU`}iIK)=aC8n7qUKiO7S?J8uy{ZG)e?VmIQa;}hq*4+vl$KP{u)$Akzq zZgsi41X$m}c%n0e_7SFkD-r$MQ2#=rHR9b!e`T-?`nS$Uy1wlfl#kKhA2BGu@~dq- z^uOW$|EGBCYz_ewmXzZ6Af5{3+ugQLNBM3l@l-;*)qUpuv^t*v)c!l zcZaC)ClII)Wm&lq;?$7 zXO-UI66UK?@)~z~3D7Pvui1z3Wq--!?8*`BuW0x78|J5Pzq_nCI7ooadLr#3n4fOm zoAN$$9PuT$_}352PbaFbL)S0?#$*O`ruZ=bx@FyWt_STm>`~biHw&(>;O*9ILwxZ2 zF8fn}1@DshYP=9%PkYol{1V1|=WYqDq7(foiMw+ai4!h(?S7x$e05uo%AfBzme7A)^Ca669i?b3;kuZMQ9K=*`V z+)a#se+hF3PGEkyV$tziw;&PLHSOZQq|1WKO`rXFFn&foJ*4s!^Ua`Kf`&H6*I3uS z2lp{w(X(~7P~S>~Hw}B==4r4%{n()$3YtXNSns^5VJ8ciZ1Yvtdx=1}`aI*r}V}9pxMzHlg<|8NeFLI4xKAE;M z{qhv%E8O1?aBabSGv`JejSKS`{cg?`I+(xqaOoNiyAffc*IFs)2n*8P_eB2k#(bbk z{e7hi3t}cS%l`!tfwx#D@HL$UI|Ys$?C{0+Hdjk5*0WK+p{;{mp+wldDE>py2iH65 zSeP4t`ONw9X@6f9@Gn|t<(whH!L@L{*`Ec`qFg>17cjqyDyowYVnO@qVwvR2_&&{e zHdyoo3wV?+3w7jRzOuFEwdY9|EN#eAZN823A5VSN9L$2nxA?mBfO(jZ;EPpNWtw0I&JdSU}l*JTt482+6088%AWXAeKAk zQNTAMNNA0OEoI?;hr_~;^b_A3xlJU69wShUaq@ z1|1a`Eyv&Qm57_@;C^jDd-&`N7L18}_FDwpPglB^D^}xvs3q61823*$7krZ&>rfs8 z)xJk@|K#Gu;WOHd`E&gys|GR&N{^S7q_m>GDB5OalS-|~h9^+uQV90N!&UVGB#^jj88wi>=(JdM9MY?!cVV}Wbn z^xd2od>=%s)6aj!f-F-#`}t6;PsBKK1=Qns3Tw5K&yhgQ%eBP#B@4Fq={Ly7i~M!kJgit~!g@^X%tjjWcV!z#PR3gjj2(Zl zJMTT>o2*pxEY@ozv+{O(d}P64g}ZAKu%5GCus(cOKMO=JH-yO!VSVJCgec!<7N{Pb z7rZ}?^{<+=nf|XVsH>`fvU`>UZSItvvk&#yV$Z!^tzs z%KIZ18+3me*6wy915w^8u9t@mk^^~DR~clepZj1iB)|sF3BSn|zGSfea4uzh4I2hr zho3r!kb!Zo>HNZ4Hq1E76P`wqVTJaZyVMPAc-?VOvMY=XhRj!K?i<-4bImP$b1c?> z<$ZR9%d?@0|6$gNi)47ZBQPRcnGI_>NvsNdFS7S1JE=;Q4bN5J(~1J*UpDopxCR>z zKbcJD$t8nCuVN})n+>UE{Bu0%WT-2-U7xg@4O*Lww(Y}uvD4t;ipsrgFs(J1J#iP` zPX$zj{o2QdO#20k>!oCPcgaFzn<*RG^V!;MRalSx?Hd+CU_)!Et0c2Fs_r~RQJMS11QB&@$Kt+_UD$OeXE>&uyFT>r(64a&xBc$DQLx;cdmYmP|VWSOzy`mLsn zHF5a<>w`)_0htYNl|n!NIzxtOuk{H-R&3x&SULM9hztw$JLe?q*kC^RdTAQp`%p6W zIg^~&aAa5OM9nEOB=Kz(3U_6Lk&tI;k{{OdQ$Jlhc!UkN-=|~*vhe-W%`!#l-LABM88{U7?HJqxJzx9hXPwfH1|btoAY7aCfNblG6M z-)bfX*E1hYonmRCJkS2Hi^uhTeDvByb~_s`%nbURM}Exslz-6bk!8b)L(b2OX!t&AL(=XWW8?do z)|bYv_+Dpdr+=^*8$!l&;>NI^FMp%ud$IC;NILOfxue8m~Pk_~Vs=dc{s*Ht4mZ#+T&UCHIEYc`DU>2iZZI?%s|9^|He z=)wA8Oj-02FWP_cTjAvntdEiM_bK9hG+)x;jju@X`1pK3Kdwie>jO3F9SI`Dw|CA; zu%Y+Ukl|QA)^EC-Auo6y}CVpNKFKgp%~AAUtXF6XosVtq6@=tl7}*6)aS-yOP!^|j4W zJ6e%n0uPvlXON(wK05jv@>|8q&$)*6aO=-5AytBG_~fy@P9%#2z6_n#F4T{C{e~T? z7*BE^uTy=D{Ek6GY6`xGy+DcgY?NZ-dxPUr*V0JPwC}I#-A!!xlGb0c0_SVKcUdJ> zfel-jYt4?}dY-;8^0~2v4I_2OC7iK7PkGAKma+}?;XC^J;XUMs`SRph#IuG--;-R( zuWFTi{zb$yHf>$p?0pix+cqzZLp)2k=dlmUQJz}TtD+Fk?t1d}o~# zZ_aG6J!N(62%gsrm&zorq_N?Sy~j}K0tp;{ScVC)(4WjQIrrfC(%LiEXVtyf(Aem< z!fFZsz5ru`rLHgfXZtVj zlMv1Zu|rX>>+yZ9!A;(h7m;kJzNYc4k&g`Ps5<1QVTdPUF`Jde@IC9Zn=bc**bw+T z>3*mb8D7a=9v=!t{4o*~ic}(lF!!PK&OkPd_8rUa-cAO*UigU*fA69yob1yj!vQ6s zD$_G){{)4OzJ2JwJkc9>gs|b#hBtQF=42Q-YS1ll67hO!@JkfN+lLXqZ`t5_|MiI6 zjyp_-?KT4MwUHk?9=2)^;*;2k9&3@~INme0JOwW@oLRp@$jYA$A2L|2s)(l-cwV2$ z!1eyESN7VDcrRtF{;CrBF;d>o*^BQrHD^xx2L_|R%Jn<-A^uhKURTOMd3Z|Q`5a3p zLpgkrdwdr0=B)`s;|RVdo$^%q7sG~+NevR4nPiCUXMQ8Z<9OO@x39-|Z##Uijh4a& zx^hru3C4pL`J!TFX>4fWS=7x$eE!;@FTN%d^+?CsLZh5Y0O8Y1b(F6W`m`y z*=a@q#@hvIna&OLr{gr2L{B`A_MiK1bcYQSPtCN1P+oigOu3uiNB`O%rW);t_)i*B zAw6Qlm5BhtUzE4ap&Wg}QtAe!2@VV0KiSb=F4Ewz0`?U?8fy(A%h^zr#!ARTdt2o0od5U&<5k&^jHfyoQskvRe8csNj;ys@y^9QX4fQ$3 z_`9+6v8`=rZ)GYVe6t$k-O^9GF2Mc!1eKzaJ&W?s0{Um^DT%MVplcH^vLj2f9uba?QGaJ{pxb&4l+cym>N5G zu_2h`aA;%;=3oE!-=|3_zUHScONP#nZ}0o?cP3Mc9FFlzFx$p5tb+|C!$Ztwj9hRReV7DR2>_}r>jpTPT}7Z2>r;NS%vwM zQn8;FzPDCd`FAh*ImU+z{2!T0c0!|INk9KV(C3ZMlF3;@hPX zQ-}T{^nZo}$rka?Ua3!H4(%z?@rT&(cM^~c)^L`iJ*jJ$cdo>Ewca9>efk>OSIDH` z#%>b${1EwfDGTwx<-C5*00|5Sww`xM$M{;Bv%0v81Pb$tu{%=HzLNzHSKxl&3K{<# zMIsya%~wvUW4xrA3L36IkN6xy`>KWU^yjg)@>`-Xzx0r8&h8{ZqE>1kX9VJTL5|lI z+;3cKSJxRtJaH&4dXR$qk!D)4=>*~pHM8;0L<0$obW?QOPhfnxmo?W{L4q5kmguL) z*zmLNXJOMr#6Q`E9XCDMprU+QqU1ISCPY~t#SAuxD<{19TY~3X5xr)4ZYb~im0~@& zNKlPs?R3QBf2*RT9%tiu;pq9<6h}5Bd6Z3`zeob^LjLE;wrmj6*_Cl57S9*$qpDM_ zFh03wNo@~DeZ`#EeHHV~8O>iQPSL2Jqz%V!n6W|QO_F*{BAzd9{g_jP@u+c(&}EMM zx3*Sk-sc#v9_-n`;KlRPm;MoVUSmA_v~*^q<`~Aa1BA$0Bh0s&=LW9h-#htbQ@!oU z0mOsNcbCJD;(1`~*3-8LY*_J$+*rXzeZ3%0rD1&9CjM+!qAS`9_qbXF#-HYW3%w3_ zzGt8Cd?69z+n8UX6JU{~4pq@swB^$~@YM^W{?!PCE9I}v~ z+TNb&(48bG{7Jm}81q-}MB8jb75w`zwkCf&aKG?nIghxs5%Z;8Cb84VkIZTd4aLnQ z;5WHEDr&%nUoZ97_Tb+`c`moXZnqvAR*CBJUcmD>+Svn5*4k_+m8L}8*ns1;OBtTt zf&S`xH{-t?!ZZ_mipTv7&|EI7sji;&$-~S^+$UKxnMQI=zl1k*3`JCejC5jMdDyk=G zqznxr%2Yy$RHP^kG$^Z9vR{9dllS!?aJ_g;ISwXW}d zUmJcw(>m||yZ#ayyd!+%9E*$>9(HQQ`l~gsGhQ_HUp;30wdL4vGK(c*_8eo;h{lP??e+XVF7Yvq}HK`gqlva2bp89erxVwNiOt&%Pi zG4=Q0yDE)Kg1k|`cSy^;UO>^Du66Uk|Kz7xEVv6jLLtog{a#P-wd=3COkWAeKTrGj zb7vOTvt5h_&jl2+#_8TG#UCOr=WeI;WMFh@Y4k;D}?1(6ccHqK0QT1PrM%u zZ<~bUQ{bEm=v%Q9V`Jl_SY#;AySghuKnGa9<8y^rWWO%Q=0vQ3ev0ItkdPaU; zPqbWj^w$jZr_h?&?%@wTXSd;8_4iDgSif}ITmkfwoo53StC{pLZNrU04)n*r*Hwx= zVN&F&b7w2og9r29G6wEIPxf1v@X1p^8+`RtXO=Q) zt%Cj?tajAkK9jcXiXKPrtrCmk_GdGxz#+~|t?pyfgZIHxN{bIqNHK^072wuh4nA1(NBG%6 zQvsc4|G1Nq1b+04{Rwm$#|y(y>rHb&x7DK zi4O8Lb&~{i*JeD&8ocJ|qYTFm1p&3l`l?L^ACha&O#dN`=k>{l+zyQIK%vj$AxZE9 zQT>1fjPJnl38zKH@cbBhZe9%@TGge^mKK8kmy=|!oC@Cew~)BnD4&>fPu(gy3f|DA z(yj~M9X(Dr-yM8&%6#FL<>1}39{Xzz?u7pt+d$kHJR_^@w4W8m|0M3*BTev(c{493 zW#atg6f4pdz|$40vaCmO{v!>$mBhd^Ql%U+%)rN{XPR;Qr355w;Z?nT59W`CZBnZ| z&aWk_DQ6pa#+z3e&nF7VbX#XdT>z8xxR>SMARcFz&i)fw0&_4@~*t!Xr zRQRWEwlL0DL4I&G4grl$E4p26_jUbGZX2k+r|hiB{J_|92f zYa_six+W#W)bN=U6*uI1_A{SuJ!*F5dorn5v1D0(JD+^&_7#ZY_;#mb#B0E}iUmF; z32sa}av(0d3H&QDNjPP;7x-9;iiP+)@YF+PL-{Va-`qbL@^#=*HGK{y;GOqPv?S;? zcxhP2gktc}8|E=0k>Gi!Oh8F;DCZVka6@VNYx8zZ;DQ#(^LIpK)s8v68oxgmJ) zmJlz&ZNz(}?qPE#_3aO;6e{G? zT>kivU##&w-t^pPd7V$4ej>qZOuz%z2g%#!gZJrw)xJTPKX%3;G2r)ga+Ryp7J}!R zg-tqk0r3i5i$emIg0CBdIljx{(?g?`mM`?sKR2+WF&*)_5#Fs6jd4GghoxXQ@*bO>YwFKJ)h@- zc*@>H&SX_4oi-o8@-!Rw`*`TbA&louxR!OaJ>m^JdTVMiew*V~JEU2Bn$g%Tdc}}Q zX|Cadc+^+2pJ*KkUOy}pE9B6t264`7wIo4GeUf-M6ZSd z-haNTQU$h#i2rY0a8C#Ii)XEAW0(O2|*UHwF1WYxZ>U zlkk7`kACI?b>s)_5L)Pw!6dawUEzI`@P4x|E&6$$N#kf_RDk-MqEtjb7ci+oX!2t% zW$@(-hUxo?m?XP4;ASf7$Ha7h-u8$|8hoCY@o)_#3M#vz}oy(#&?-VyNZoF`gB$hW%Toe;ba z>j%4Sd#b%rAr=Lgn?ey5|rn9J0f^*9J zJMy)}UnG`deRKNnhQiCeJX)cU&U|XXqU|RxJzVycM=A=LtF*D6s>v}Ks6>9?`bC~g z_M5V(ztQ)_;w~O-`+21(2J5%AVcMw=+L0gW*voEpV9{Q0@zrjv$e;I*`CeziqN~+s zTbth_pY!sMotm~R`VpRWnp@AKt|=>?aarKCu5sR?b;u{|5q;u}_$<{ozq_Jd^5|Q2 zm_!WLfB7475C4OFMpZEcnme&5SJd~s;5qVxMdfZbxPtF~x^r2&0)7umc`_vo7QH5^ zkn80W^#Xohz2?A4R~O5MPFjaAK!N>1~rLf1yq3Gs3pWKH$8Pa@w^ zIa_%-*57YCSwV3L7|)fJ$vPp3uai3@+8+bI0b{Fk(&1R2zmaHN9f|ApZR?^>2f@#G zj>;B>^T?CYVt*$Z@qtGSs=|&4Uj{%WTo;){*U^~C&; zmpM0F#G(wt(^0XLG|Q4<3oV z-Tl1a4~vf3j91NY{nszVYO1{m;=N79v?|cvDkqGWg#F;YVUc44&OD0yk{dEh8Tv!d z3T=&5JeoNrA~;aPfl7{zUbl6`{I8Q`|Iu=woQJw^w>a_Ww_D&S%MkvRF86oGF_F*j zyejOqu>)N)XKUSd!2BpT{Gee8{rhM8s~TJQYYaL_oW*{4*?Z&RVSD8NK97<6;o?Az zzgsWLTj6_xR4o`A>e(&fSjl8`r;xH-|@g34PkdyRpCS-etkpgr9-LxpfnuKd8<*m$Py@ z{3ty7Z(Z67e^Gx&f7dCPKTfweJ!tQd5^!jh5|3_{$_?p2Z&+S#KI(*g{Qb4MTfNYp zXOzZTE5{>`%pYE5(as!a#~$XK$g_zjmX z#6(=U2YwK?o%Qj21M=P5`cs!&f&LwKVEjSk@2jkgSga17UEK1w78Vy+v22tHN^!;(Xn~ng5DQn%fYU^U{GzPgO{%p#8lmuPy!G!@qQLuR>=fmxT7m zc3ZdLcpi_UULgO!UEeXh1bSN4PFtH~*2KuPwYMoBGLim?ec|VIYWYfU0 zDT6=Y7j&(cvswoFY>82Qc4`5a4xKD_cnH1ow*LdS;Oku4&@{oo+z~vY{!IAPJoqs= zI~R6C54KuVT0WAG`_cNNOW*~){Kg@llQ~>^8+E*2&X-LyCk`9!$c4X}TJDZm@Jywn zGfu6(4F91=Pn0CVyQ>-%qCUYtsVw_q=_T;YvlrBNJjurWiElkO=LqsgUKWjJU*saM z!sApX`1Avjnd*93TvC+HZlN^r;FMp}WzoJ|rg7pO@a?ityFIb#T*^;yD0IFAy><`p zOu;GmOC{F0*QI0p6D?yFoPfVoq=S&}SvI{iG2gl_iA%fpi=Fs*9K2agXSg|*OJU(w zEI;t%tJRUC-jQ7Tn7?J-d+_u!v!e3j2hsne!CY75KRMQi1gnH{sWEuy%%em$H6Ah9 z+Z4v7pn^#Y95TQYRP;mL_i!n_xoYk8EH;VMa3!j@;rMPwnF<(h{9SQU#j=bYB$z_JrS=HDj4;74&`S0Q-C|E>-H-vfJMx|4Z@OX&HC;FKs+r zs@4R)5WGHmBmCeVL@3Q?w%~pzXX{>Lx`;VrT?}wl3t-KAC#WAmrc8yJ?v}Q;$nKC4+Wsh2P?1 z(>sc>y6^+rx#npD{215hY>Iz1i%Y%A$x4;7$p33mt5lfIrBA-9+BXz%yrzGy9h!>y z5K%2U3C9a~_r^Sd_KuW1mts9`4n9SGpMlph^%zy_kP6ZA!}Xmg!9BRDheMLK^F5{8 z;NO||v%KdkhvxY*W_Y{-kBW1j5DWiFm7X&TbDkkTQ09DTa|eeWHa{M=evbH&)t(0p znmJ^5%qXeh1@hCZYkL}+;FqaC>4`=u_%=VlTeqG=H=le@yIA(``yrpzScl{Jx}W$~ z1pZNF`M2Rq4mkvJt_5H|bsRH#-3foqgYjA6`FU)r+*5F!$~pAdvEYm5RW`Y+8HWzT zpYz}$@zXw-zn+d>v)oEJbTOdU`x)XhLN~2-oB}^jPoJ1kQ9Lj4DrFCL7IDb$?NUuY zJf9DxHB2oEIiyJoZlxV%6F12yJo5^Nb}*0I^x*j!=i;I??jnaSopsGK!SmBy9XIH8 zmP6gM1HSCR^Vk?Iv+!;jhjyQR|JQXq|1Ul2bCXkXeJw4vw@0%PTqfH*kO)6mk3+(W zv4{t-*N=Rgfc6JHPg{ij*PnIl{>tA;_{9WU=}X@eGsFTKT|RDu1N z*>tI6QSo`0S}VuGq*y%8QAdDW+R`wW95{Y1q+ z4Z$C1_H&y{i46Kyu5tgvPxy5_9`B-_!l3`Fy`5gE|1PZmI|?uOKEA}DZEBpO5!g?- zY}u)7e}O^AG#J|#jf1{1lyliU^}olz6uJJ_emOp^ul4#_RmPyi7Okjb*zZIq2ed6M zVNh!9yny`JA<(Ersw-}bRk7xqgT<>5nZEt26dZByR5=~tuOU#abh)t?7t1p=~shZV4>`p_LxCGon4QfME%w5 z309f6|2w|_|2>DwSE<($qZric#0}`L@dY2Qb8SZ1-TI$}O0~WeuNKsm6~dtL{1xvE zP}YCCYN22Yg9H~EJRMOM)CVon^<|LG&*@fg@O=n_pa!Ln>lx%`{wne&%GZ;AcdDcQ zO`W>)NhmKrwKUKj{f#}oMN*dZ*U!Fme?x%Qa596gypOmP*yBrn@&0+!PcrDt(NE@K zC~ujkKHL_~AX^#Bj#1QKk?_x<5|qcbFL&UiTmSZ@`E91_W(qSY>k?PJ1K(qityYOM zKso!JwdiM*wcRJj4gSRZDZbRrLVcmtRh%xA$F?_Ewa;Mrb`C`foP#9s|La(~?5f=k z4#~_;^lr!huYOm=KP@46AGn@$ZUzsppDW&27tEpWf7xV-`LN+H^uOEhIusMG=mWnZ zR$$!YhxT;D>qR$zJ^rtil&Vg?uqTHB=^61n_=~dV9=58qr;9crGA3N^|9pRAZ2LD$ z_4*fMeSKlgJ6^9GgYxE2lUH%U`rTbD&`X^`73x71bFhAo8@Kd+sS1Nc%R*MzVtvmm zX_&c2ia}%BtI3RBt$ojN!vIkRVQmdkX?ZM@%LLYhWZ)wnq$zNf%?Op z>Bh2M3>w?Nm%6-9?;_|s!Zkw`GSKg4j5iXPKyP|v(R>x}qZgY$>Ve)f7|+=Di|6kwbAx{P zboK6qs_@sDnS>C|7Mh(H>-z#ezSY}eXQRH znR2}^QC7)L%DX!a`O&E-R<8ZvM-PuJ>5j*GzGo1|vhd&z!*U%b)HG!?`?;ldzG%Wr#hPy0>zr!P@a}|Ma&!JO}O1zoMoRM9e;M2UKn9LJtKcn7s`EN zW-(>4lL-5jhqO7^^S+yPizq%_2kAiHFXZ9 z+Be;|wP90#+QA3;GvIgsB}C`G6`QtI$`n{=b4co8=4M?>tgnP6_f^9$Utzj?Y&^y{ z)<3`6_g5ot2;Uz$D6h2v^XXbt`kdJ)E3ME-jhw}%C81uoQhV_|2WHddH?#jczsB|# zVY1?N;9L#`{$9J^9M@~As>Rk0O^k2Lko`-Zfkk>xMhmJkJvE%)JH~QzX3z_2mo(^K0qrW=d&T{BHq;pY4KF|TX Mb`((~1_qNi0NeBt+W-In diff --git a/man/interpolate_variable.Rd b/man/interpolate_variable.Rd index 5338bdc..758a789 100644 --- a/man/interpolate_variable.Rd +++ b/man/interpolate_variable.Rd @@ -20,7 +20,8 @@ interpolate_variable( methods = c("NN", "IDW", "IDW4", "Exp", "Sph", "Bes", "Gau", "Cir", "Mat", "Ste", "Tps"), bbox = NULL, - return_raster = FALSE + return_raster = FALSE, + outputfilefun = NULL ) } \arguments{ @@ -53,6 +54,8 @@ interpolate_variable( specific bounds, or a 4-element named vector with values xmin, xmax, ymin, and ymax.} +\item{outputfilefun}{Function of the form f(mask, year, method, var) that returns a string holding the full path name where the output raster should be stored. If not included, default is [coldpool_dir]/output/raster/[mask]/[var]/[mask]_[method]_[year]_[var].tif} + \item{return.raster}{Should a raster be returned?} } \value{ From 62b54302f1c792b66dcb7c3e477c5f6a7ff5676b Mon Sep 17 00:00:00 2001 From: Kelly Kearney Date: Tue, 13 Feb 2024 14:24:04 -0800 Subject: [PATCH 5/6] Minor text updates to cold_pool_for_regional_models.Rmd --- analysis/cold_pool_for_regional_models.Rmd | 36 +++++++++------- analysis/cold_pool_for_regional_models.html | 46 +++++++++++---------- 2 files changed, 46 insertions(+), 36 deletions(-) diff --git a/analysis/cold_pool_for_regional_models.Rmd b/analysis/cold_pool_for_regional_models.Rmd index 4466f0c..87e2a6d 100644 --- a/analysis/cold_pool_for_regional_models.Rmd +++ b/analysis/cold_pool_for_regional_models.Rmd @@ -13,18 +13,23 @@ addr: l3: Seattle, WA 98115 --- -# Introduction +# 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 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 following files must be in place prior to running this example: +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. -## Setup +For details regarding how those files are built, see the [beringnpz/Bering10KPostprocessing](https://github.com/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: @@ -33,13 +38,18 @@ The following packages are used in this script, and should be loaded prior to ru - lubridate - ncdf4 -We start with some general setup that sets output figure resolution and the map projection to use for our calculations. +We start with some general setup that loads these packages and sets output figure resolution and the map projection to use for our calculations. ```{r setup_nodevel, eval=FALSE} library(coldpool) ``` ```{r setup_devel, include=FALSE} -devtools::load_all() +develflag = TRUE +if (develflag) { + devtools::load_all() +} else { + library(coldpool) +} ``` ```{r setup, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} library(akgfmaps) @@ -67,7 +77,7 @@ simname = "B10K-K20_CORECFS" # weekly-archived "operational" hindcast lev3path <- file.path(moxdir, "roms_for_public", simname, "Level3") ``` -## Survey-replicated data +## 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. @@ -76,7 +86,7 @@ The first step of this process is to build a survey-replicated dataset. We do t ```{r csv_path, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} srepcsv_path <- file.path(lev3path, sprintf("survey_replicates_%s.csv", simname)) ``` -## Masks +## 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. @@ -117,7 +127,7 @@ maskname <- c("SEBS-noNW", "SEBS", "SEBS-noNW-modelshelf", "SEBS-modelshelf") ``` -## Cold pool index netCDF file +## 1.4 Cold pool index netCDF file Model-based indices will be saved to a prexisting netCDF file. @@ -125,7 +135,7 @@ Model-based indices will be saved to a prexisting netCDF file. cpoolnc_path <- file.path(lev3path, sprintf("%s_coldpool.nc", simname)) ``` -# Generate bottom temperature rasters +# 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. @@ -201,14 +211,10 @@ for (ii in 1:length(year_vec)) { } ``` -# Annual indices: calculate and save cold pool and mean bottom temperature metrics to file +# 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. -The resulting metrics are saved to an existing netCDF file that also includes corresponding metrics based on the model values on July 1 of each simulation year (those metrics use the same masks as in this example but do not involve any interpolation; instead, metrics are averaged across the model's native grid cells.) - -We include some quick checks in this section to make sure we don't overwrite data already in the file, and to flag if our most recent calculations diverge from any previously-saved values. - ```{r annualindices, include=TRUE, message=FALSE, warning=FALSE, echo = TRUE} # Open netCDF file diff --git a/analysis/cold_pool_for_regional_models.html b/analysis/cold_pool_for_regional_models.html index ab98765..ca5ed49 100644 --- a/analysis/cold_pool_for_regional_models.html +++ b/analysis/cold_pool_for_regional_models.html @@ -360,7 +360,7 @@

Kelly Kearney

-

Introduction

+

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 @@ -372,16 +372,28 @@

Introduction

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 following files must be in place prior to running this -example:

+

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.

-

Setup

+

1.1 Setup

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

    @@ -390,8 +402,9 @@

    Setup

  • lubridate
  • ncdf4
-

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

+

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)
@@ -417,7 +430,7 @@ 
Setup

 lev3path <- file.path(moxdir, "roms_for_public", simname, "Level3")
-

Survey-replicated data

+

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, @@ -438,7 +451,7 @@

Survey-replicated data

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

Masks

+

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 @@ -488,13 +501,13 @@

Masks

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

Cold pool index netCDF file

+

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))
-

Generate bottom temperature rasters

+

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 @@ -570,21 +583,12 @@

Generate bottom temperature rasters

} }
-
-

Annual indices: calculate and save cold pool and mean bottom -temperature metrics to file

+
+

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.

-

The resulting metrics are saved to an existing netCDF file that also -includes corresponding metrics based on the model values on July 1 of -each simulation year (those metrics use the same masks as in this -example but do not involve any interpolation; instead, metrics are -averaged across the model’s native grid cells.)

-

We include some quick checks in this section to make sure we don’t -overwrite data already in the file, and to flag if our most recent -calculations diverge from any previously-saved values.

# Open netCDF file
 
 nc <- nc_open(cpoolnc_path)

From 564b8d148b5c2a12c7a173adab4cc4a08c03f31c Mon Sep 17 00:00:00 2001
From: Kelly Kearney 
Date: Wed, 4 Sep 2024 15:16:02 -0700
Subject: [PATCH 6/6] Correct typos (missing write mode and function name typo)

---
 analysis/cold_pool_for_regional_models.Rmd | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/analysis/cold_pool_for_regional_models.Rmd b/analysis/cold_pool_for_regional_models.Rmd
index 87e2a6d..f9a71a3 100644
--- a/analysis/cold_pool_for_regional_models.Rmd
+++ b/analysis/cold_pool_for_regional_models.Rmd
@@ -219,7 +219,7 @@ Metrics of mean bottom temperature and cold pool fractions are calculated for bo
 
 # Open netCDF file
 
-nc <- nc_open(cpoolnc_path)
+nc <- nc_open(cpoolnc_path, write=TRUE)
 
 # Read relevant coordinate data
 
@@ -297,7 +297,7 @@ for (ii in 1:length(year_vec)) {
     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"))
+    ncatt_put(nc, 0, "history", paste(hisnew, hisold$value, sep="\n"))
     
   }
 }