Merge pull request #49 from afsc-gap-products/dev

Fix interpolate_variable

.gitignore

@@ -52,4 +52,10 @@ vignettes/*.pdf
 *.tiff
 /plots/cpchange
 /analysis/variability/
-/analysis/temperature_eof/
+/analysis/temperature_eof/
+
+/data/*.sbx
+/data/*.prj
+/data/*.shp.xml
+/data/*.cpg
+/data/*.sbn

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: coldpool
 Type: Package
 Title: AFSC/RACE Groundfish Assessment Program EBS and NBS temperature products
-Version: 3.0-2
+Version: 3.1-1
 Authors@R: c(person("Sean", "Rohan", email = "sean.rohan@noaa.gov", role = c("aut", "cre")),
     person("Lewis", "Barnett", email = "lewis.barnett@noaa.gov", role = c("aut", "ctb")),
     person("Kelly", "Kearney", role = c("ctb")),

R/data2raster.R

@@ -53,21 +53,23 @@ data2raster <- function(dat,
 { 
 
 
- if(is.character(bbox)) {
-   in_bbox <- NULL
- }
-
-  if(is.numeric(bbox)) {
-      bbox <- sf::st_bbox(st_buffer(mymask, 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
-      in_bbox <- bbox + c(xshift,yshift,xshift,yshift)
+  mymask <- parseregion(select.region, interpolation.crs)
+
+  if(is.character(bbox)) {
+    in_bbox <- NULL
+  }
+
+  if(is.null(bbox)) {
+    bbox <- sf::st_bbox(st_buffer(mymask, 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
+    in_bbox <- bbox + c(xshift,yshift,xshift,yshift)
   }
 
 
   for(ii in 1:length(fittypes)) {
-     fit_terra <- interpolate_variable(
+    fit_terra <- interpolate_variable(
       dat = dat,
       dat.year = dat.year,
       var.col = var.col,
@@ -93,8 +95,6 @@ data2raster <- function(dat,
     }
   }
 
-  mymask <- parseregion(select.region, interpolation.crs)
-
   out <- out |>
     akgfmaps::rasterize_and_mask(mymask)
 
@@ -366,6 +366,7 @@ calcindices_temp <- function (datafile,
                               bottomvar = "gear_temperature",
                               surfacevar = "surface_temperature") 
 {
+
   # Read temperature data from .csv file
 
   print("Reading temperature data...")
@@ -445,7 +446,7 @@ calcindices_temp <- function (datafile,
         mean(na.rm = TRUE)
 
       lt100_temp <- terra::mask(ras, lt100_strata, touches = FALSE)
-      bt_df$MEAN_BT_LT100M[i] <- terra::values(lt100_temp ) |> 
+      bt_df$MEAN_BT_LT100M[i] <- terra::values(lt100_temp) |> 
         mean(na.rm = TRUE)
 
       for (it in 1:nthresh) {

R/interpolate_variable.R

@@ -12,7 +12,7 @@
 #' @param cell.resolution Dimensions of interpolation grid cell in meters.
 #' @param nm Maximum number of cells to use for interpolation. 
 #' @param pre Prefix for file names in output (in development.)
-#' @param select.region Region for interpolation as a character string. Options = "ebs", "sebs", "nbs"
+#' @param select.region Region for interpolation as a character string OR an sf object. Character options = "ebs", "sebs", "nbs". If select.region is an sf object, must provide bbox.
 #' @param methods To use as a character vector. Valid choices: "NN", "IDW", "IDW4", "Exp", "Sph", "Bes", "Gau", "Cir", "Mat", "Ste", "Tps"
 #' @param bbox Bounding box for the interpolated grid.  Can be either a string 
 #'    ("sebs", "bs.south", "nebs", "ebs", "bs.north", "bs.all") corresponding to 
@@ -50,10 +50,6 @@ interpolate_variable <- function(dat,
     stopifnot("interpolate_variable: Can only choose one method if the raster is returned" = length(methods) == 1)
   }
 
-  if(!dir.exists(here::here("output", "raster", select.region, var.col))) {
-    dir.create(here::here("output", "raster", select.region, var.col), recursive = TRUE)
-  }
-
   names(dat)[which(names(dat) == var.col)] <- "var.col"
   names(dat)[which(names(dat) == lat.col)] <- "lat.col"
   names(dat)[which(names(dat) == lon.col)] <- "lon.col"
@@ -67,9 +63,29 @@ interpolate_variable <- function(dat,
       dplyr::filter(!is.na(var.col))
   }
 
-  # Load EBS survey exent for masking ----
-  region_layers <- akgfmaps::get_base_layers(select.region = select.region, 
-                                           set.crs = interpolation.crs)
+
+  # 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
+  } 
+
+  # 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 ----
 
@@ -212,12 +228,12 @@ interpolate_variable <- function(dat,
 
     # Write masked raster
     akgfmaps::rasterize_and_mask(sgrid = fit, 
-                                 amask = region_layers$survey.area) |>
+                                 amask = region_mask) |>
     coldpool::make_raster_file(filename = here::here("output", 
                                                      "raster", 
-                                                     select.region,
+                                                     fpath_subdir,
                                                      var.col, 
-                                                     paste0(select.region, "_", tolower(methods[ii]), "_", dat.year, "_", var.col, ".tif" )),
+                                                     paste0(fpath_subdir, "_", tolower(methods[ii]), "_", dat.year, "_", var.col, ".tif" )),
                                format = "GTiff",
                                overwrite = TRUE,
                                layer_name = dat.year)

analysis/roms_level3_coldpool_part2.R

@@ -1,21 +1,52 @@
 # roms_level3_coldpool_part2.R
 #
 # This script is part of the post-processing suite for the Bering10K ROMS 
-# hindcast simulation.  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)
+# 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)
-# source('extended_utils.R') Functions are now included in the package
 
-# Index file
+#----------------
+# 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")
 
-moxdir = "~/Documents/Research/Working/mox_bumblereem/"
-fname <- file.path(moxdir, "roms_for_public", "B10K-K20_CORECFS", "Level3", "B10K-K20_CORECFS_coldpool.nc")
-# fname = here::here("test.nc")
-ncin <- nc_open(fname)
+# 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
 
@@ -28,7 +59,7 @@ thresh <- ncvar_get(ncin, "threshold")
 
 nc_close(ncin)
 
-ismiss <- meanbtmp[1,2,] > 1e36 # Don't know how to query fill value in R...
+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]
@@ -43,7 +74,7 @@ cell_resolution <- 5000
 
 akSEBS <- akgfmaps::get_base_layers(select.region = "bs.south", set.crs = proj_crs)
 
-b10_lte_200 <- sf::st_read(here::here("data", "B10K_lte200m_polygon.shp"), quiet = TRUE) |>
+b10_lte_200 <- sf::st_read(modelshelfpolyfile, quiet = TRUE) |>
   st_set_crs("+proj=longlat") |>
   sf::st_transform(crs = sf::st_crs(proj_crs)) 
 
@@ -59,22 +90,20 @@ mymask <- list(akSEBSnonw,
 
 # Calculate indices using survey and survey-replicated data
 
-datafile <- here::here("data", "survey_replicates_B10K-K20_CORECFS.csv")
-
 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(datafile, mymask[[ii]],
+  svy  <- calcindices_temp(srepdatafile, mymask[[ii]],
                            select_years = select_years, 
                            bottomvar="gear_temperature", 
                            surfacevar=NULL, 
                            threshold=thresh)
 
-  srep <- calcindices_temp(datafile, mymask[[ii]],
+  srep <- calcindices_temp(srepdatafile, mymask[[ii]],
                            select_years = select_years, 
                            bottomvar="model_bottom_temp", 
                            surfacevar=NULL,
@@ -89,7 +118,7 @@ for (ii in 1:length(mymask)) {
 
 # Write to file
 
-ncout <- nc_open(fname, write=TRUE)
+ncout <- nc_open(hindcastcpoolfile, write=TRUE)
 
 ncvar_put(ncout, "average_bottom_temp", meanbtmp)
 ncvar_put(ncout, "cold_pool_index", cpindex)