Automate R script

README.md

@@ -4,7 +4,7 @@
 # Regionalization and Formulation Testing and Selection (RaFTS)
 
 **Description**:  
-The formulation-selector tool, aka Regionalization and Formulation Testing & Selection (RaFTS), is under development. For more information, see the [Wiki](https://github.com/NOAA-OWP/formulation-selector/wiki). 
+The formulation-selector tool, aka Regionalation and Formulation Testing & Selection (RaFTS), is under development. For more information, see the [Wiki](https://github.com/NOAA-OWP/formulation-selector/wiki). 
 
 As NOAA OWP builds the model-agnostic NextGen framework, the hydrologic modeling community will need to know how to optimally select model formulations and estimate parameter values across ungauged catchments. This problem becomes intractable when considering the unique combinations of current and future model formulations combined with the innumerable possible parameter combinations across the continent. To simplify the model selection problem, we apply an analytical tool that predicts hydrologic formulation performance (Bolotin et al., 2022, Liu et al., 2022) using community-generated data. The regionalization and formulation testing and selection (RaFTS) tool readily predicts how models might perform across catchments based on catchment attributes. This decision support tool is designed such that as the hydrologic modeling community generates more results, better decisions can be made on where formulations would be best suited.
 

pkg/proc.attr.hydfab/R/proc_attr_grabber.R

@@ -105,16 +105,16 @@ retrieve_attr_exst <- function(comids, vars, dir_db_attrs, bucket_conn=NA){
 }
 
 
-proc_attr_std_hfsub_name <- function(comid,custom_name='', ext='gpkg'){
+proc_attr_std_hfsub_name <- function(comid,custom_name='', fileext='gpkg'){
   #' @title Standardidze hydrofabric subsetter's local filename
   #' @description Internal function that ensures consistent filename
   #' @param comid the USGS common identifier, generated by nhdplusTools
   #' @param custom_name Desired custom name following 'hydrofab_'
-  #' @param ext file extension of the hydrofrabric data. Default 'gpkg'
+  #' @param fileext file extension of the hydrofrabric data. Default 'gpkg'
 
   hfsub_fn <- base::gsub(pattern = paste0(custom_name,"__"),
                          replacement = "_",
-                         base::paste0('hydrofab_',custom_name,'_',comid,'.',ext))
+                         base::paste0('hydrofab_',custom_name,'_',comid,'.',fileext))
   return(hfsub_fn)
 }
 
@@ -185,7 +185,7 @@ proc_attr_usgs_nhd <- function(comid,usgs_vars){
 }
 
 
-proc_attr_hf <- function(comid, dir_db_hydfab,custom_name="{lyrs}_",ext = 'gpkg',
+proc_attr_hf <- function(comid, dir_db_hydfab,custom_name="{lyrs}_",fileext = 'gpkg',
                          lyrs=c('divides','network')[2],
                          hf_cat_sel=TRUE, overwrite=FALSE){
 
@@ -195,7 +195,7 @@ proc_attr_hf <- function(comid, dir_db_hydfab,custom_name="{lyrs}_",ext = 'gpkg'
   #' @param comid character class. The common identifier USGS location code for a surface water feature.
   #' @param dir_db_hydfab character class. Local directory path for storing hydrofabric data
   #' @param custom_name character class. A custom name to insert into hydrofabric file. Default \code{glue("{lyrs}_")}
-  #' @param ext character class. file extension of hydrofabric file. Default 'gpkg'
+  #' @param fileext character class. file extension of hydrofabric file. Default 'gpkg'
   #' @param lyrs character class. The layer name(s) of interest from hydrofabric. Default 'network'.
   #' @param hf_cat_sel boolean. TRUE for a total catchment characterization specific to a single comid, FALSE (or anything else) for all subcatchments
   #' @param overwrite boolean. Overwrite local data when pulling from hydrofabric s3 bucket? Default FALSE.
@@ -204,7 +204,7 @@ proc_attr_hf <- function(comid, dir_db_hydfab,custom_name="{lyrs}_",ext = 'gpkg'
   # Build the hydfab filepath
   name_file <- proc.attr.hydfab:::proc_attr_std_hfsub_name(comid=comid,
                                    custom_name=glue::glue('{lyrs}_'),
-                                   ext=ext)
+                                   fileext=fileext)
   fp_cat <- base::file.path(dir_db_hydfab, name_file)
 
   if(!base::dir.exists(dir_db_hydfab)){
@@ -225,7 +225,7 @@ proc_attr_hf <- function(comid, dir_db_hydfab,custom_name="{lyrs}_",ext = 'gpkg'
 
   # Read the hydrofabric file gpkg for each layer
   hfab_ls <- list()
-  if (ext == 'gpkg') {
+  if (fileext == 'gpkg') {
     # Define layers
     layers <- sf::st_layers(dsn = fp_cat)
     for (lyr in layers$name){
@@ -469,7 +469,7 @@ proc_attr_gageids <- function(gage_ids,featureSource,featureID,Retr_Params,
       # Retrieve the variables corresponding to datasets of interest & update database
       loc_attrs <- try(proc.attr.hydfab::proc_attr_wrap(comid=comid,
                                                     Retr_Params=Retr_Params,
-                                                    lyrs='network',overwrite=FALSE))
+                                                    lyrs=lyrs,overwrite=FALSE))
       if("try-error" %in% class(loc_attrs)){
         message(glue::glue("Skipping gage_id {gage_id} corresponding to comid {comid}"))
       }

pkg/proc.attr.hydfab/flow/fs_attrs_grab.R

@@ -12,6 +12,7 @@
 #' term 'gage_id' is used as a variable in glue syntax to create featureID
 #' @seealso [fs_proc] A python package that processes input data for the
 #' formulation-selector
+#' @usage Rscript fs_attrs_grab.R "/path/to/attribute_config.yaml"
 
 # Changelog / Contributions
 #   2024-07-24 Originally created, GL
@@ -21,61 +22,81 @@ library(yaml)
 library(ncdf4)
 library(proc.attr.hydfab)
 library(glue)
+
 # TODO is AWS_NO_SIGN_REQUEST necessary??
 # Sys.setenv(AWS_NO_SIGN_REQUEST="YES")
 
-# TODO create config yaml.
-
-# TODO read in config yaml, must populate NA for items that are empty.
-
-# Define input directory:
-# TODO change this to reading the standardized metadata, not the generated data
+# Define command line argument
+cmd_args <- commandArgs("trailingOnly" = TRUE)
 
-# raw_config <- yaml::read_yaml("/Users/guylitt/git/formulation-selector/scripts/eval_ingest/xssa/xssa_attr_config.yaml")
-raw_config <- yaml::read_yaml("/Users/guylitt/git/formulation-selector/scripts/eval_ingest/SI/SI_attr_config.yaml")
+if(base::length(cmd_args)!=1){
+  warning("Unexpected to have more than one argument in Rscript fs_attrs_grab.R /path/to/attribute_config.yaml.")
+}
 
+# Read in config file, e.g.  "~/git/formulation-selector/scripts/eval_ingest/SI/SI_attr_config.yaml"
+path_attr_config <- cmd_args[1] # "~/git/formulation-selector/scripts/eval_ingest/xssa/xssa_attr_config.yaml"
+raw_config <- yaml::read_yaml(path_attr_config)
 
-datasets <- ds <- raw_config$formulation_metadata[[grep("datasets",raw_config$formulation_metadata)]]$datasets#c("juliemai-xSSA",'all')[1] # A listing of datasets to grab attributes. Dataset names match what is inside dir_std_base.  'all' processes all datasets inside dir_std_base.
-#ds_nc_filenames <- c('juliemai-xSSA_Raven_blended.nc','*.nc')[1]
-
+# A listing of datasets to grab attributes. Dataset names match what is inside dir_std_base.  'all' processes all datasets inside dir_std_base.
+datasets <- raw_config$formulation_metadata[[grep("datasets",
+        raw_config$formulation_metadata)]]$datasets #c("juliemai-xSSA",'all')[1]
 
+# Define directory paths from the config file
 home_dir <- Sys.getenv("HOME")
-dir_base <- file.path(home_dir,'noaa','regionalization','data')
-
-dir_std_base <- file.path(dir_base,"input","user_data_std") # The location of standardized data generated by fs_proc python package
-dir_db_hydfab <- file.path(dir_base,'input','hydrofabric') # The local dir where hydrofabric data are stored to limit s3 connections
-dir_db_attrs <- file.path(dir_base,'input','attributes') # The parent dir where each comid's attribute parquet file is stored in the subdirectory 'comid/', and each dataset's aggregated parquet attributes are stored in the subdirectory '/{dataset_name}
-
-
-s3_base <- "s3://lynker-spatial/tabular-resources" # s3 path containing hydrofabric-formatted attribute datasets
-s3_bucket <- 'lynker-spatial' # s3 bucket containing hydrofabric data
-
-s3_path_hydatl <- glue::glue('{s3_base}/hydroATLAS/hydroatlas_vars.parquet') # path to hydroatlas data formatted for hydrofabric
+dir_base <- glue::glue(base::unlist(raw_config$file_io)[['dir_base']])#file.path(home_dir,'noaa','regionalization','data')
+dir_std_base <- glue::glue(base::unlist(raw_config$file_io)[['dir_std_base']]) #file.path(dir_base,"input","user_data_std") # The location of standardized data generated by fs_proc python package
+dir_db_hydfab <- glue::glue(base::unlist(raw_config$file_io)[['dir_db_hydfab']]) # file.path(dir_base,'input','hydrofabric') # The local dir where hydrofabric data are stored to limit s3 connections
+dir_db_attrs <- glue::glue(base::unlist(raw_config$file_io)[['dir_db_attrs']])  # file.path(dir_base,'input','attributes') # The parent dir where each comid's attribute parquet file is stored in the subdirectory 'comid/', and each dataset's aggregated parquet attributes are stored in the subdirectory '/{dataset_name}
+
+# Read s3 connection details
+s3_base <- base::unlist(raw_config$hydfab_config)[['s3_base']]#s3://lynker-spatial/tabular-resources" # s3 path containing hydrofabric-formatted attribute datasets
+s3_bucket <- base::unlist(raw_config$hydfab_config)[['s3_bucket']] #'lynker-spatial' # s3 bucket containing hydrofabric data
+
+# s3 path to hydroatlas data formatted for hydrofabric
+if ("s3_path_hydatl" %in% names(base::unlist(raw_config$attr_select))){
+  s3_path_hydatl <- glue::glue(base::unlist(raw_config$attr_select)[['s3_path_hydatl']])  # glue::glue('{s3_base}/hydroATLAS/hydroatlas_vars.parquet')
+} else {
+  s3_path_hydatl <- NULL
+}
 
 # Additional config options
-hf_cat_sel <- c("total","all")[1] # total: interested in the single location's aggregated catchment data; all: all subcatchments of interest
-ext <- 'gpkg'
-attr_sources <- c("hydroatlas","usgs") # "streamcat",
-# TODO communicate to user that these are standardized variable names
-ha_vars <- c('pet_mm_s01', 'cly_pc_sav', 'cly_pc_uav','cly_pc_sav','ari_ix_sav') # hydroatlas variables
-sc_vars <- c() # TODO look up variables. May need to select datasets first
-usgs_vars <- c('TOT_TWI','TOT_PRSNOW','TOT_POPDENS90','TOT_EWT','TOT_RECHG','TOT_PPT7100_ANN','TOT_AET','TOT_PET','TOT_SILTAVE','TOT_BASIN_AREA','TOT_BASIN_SLOPE','TOT_ELEV_MEAN','TOT_ELEV_MAX','TOT_Intensity','TOT_Wet','TOT_Dry' ) # list of variables retrievable using nhdplusTools::get_characteristics_metadata()
+hf_cat_sel <-  base::unlist(raw_config$hydfab_config)[['hf_cat_sel']] #c("total","all")[1] # total: interested in the single location's aggregated catchment data; all: all subcatchments of interest
+ext <- base::unlist(raw_config$hydfab_config)[['ext']] # 'gpkg'
+
+#-----------------------------------------------------
+# Variable listings:
+names_attr_sel <- base::unlist(base::lapply(raw_config$attr_select,
+                                            function(x) base::names(x)))
+
+# Transform into single named list of lists rather than nested sublists
+idxs_vars <- base::grep("_vars", names_attr_sel)
+var_names <- names_attr_sel[idxs_vars]
+sub_attr_sel <- base::lapply(idxs_vars, function(i)
+  raw_config$attr_select[[i]][[1]])
+base::names(sub_attr_sel) <- var_names
+
+# Subset to only those non-null variables:
+sub_attr_sel <- sub_attr_sel[base::unlist(base::lapply(sub_attr_sel,
+                          function(x) base::any(!base::is.null(unlist(x)))))]
+var_names_sub <- names(sub_attr_sel)
 #-----------------------------------------------------
+message(glue::glue("Attribute dataset sources include the following:\n
+  {paste0(var_names_sub,collapse='\n')}"))
 
-# TODO generate this listing structure based on what is provided in yaml config
-#   & accounting for empty entries
+message(glue::glue("Attribute variables to be acquired include :\n
+  {paste0(sub_attr_sel,collapse='\n')}"))
 
-Retr_Params <- list(paths = list(# Note that if a path is provided, ensure the
+Retr_Params <- base::list(paths = base::list(
+  # Note that if a path is provided, ensure the
   # name includes 'path'. Same for directory having variable name with 'dir'
                         dir_db_hydfab=dir_db_hydfab,
                         dir_db_attrs=dir_db_attrs,
                         s3_path_hydatl = s3_path_hydatl,
                         dir_std_base = dir_std_base),
-                   vars = list(usgs_vars = usgs_vars,
-                               ha_vars = ha_vars,
-                               sc_vars = sc_vars),
+                   vars = sub_attr_sel,
                    datasets = datasets
        )
+# PROCESS ATTRIBUTES
 
 ls_comids <- proc.attr.hydfab:::grab_attrs_datasets_fs_wrap(Retr_Params,overwrite = TRUE)
 

pkg/proc.attr.hydfab/inst/extdata/xssa_attr_config_all_vars_avail.yaml

@@ -16,7 +16,7 @@ file_io: # May define {home_dir} for python's '{home_dir}/string_path'.format(ho
   - 'save_type': 'netcdf' #  Required. Use 'csv' to create a directory structure & save multiple files. May also save as hierarchical files 'netcdf' or 'zarr', or if 'csv' chosen, a directory structure is created
   - 'save_loc': 'local' #  Required.  Use 'local' for saving to a local path via dir_save. Future work will create an approach for 'aws' or other cloud saving methods
   - 'dir_base' : '{home_dir}/noaa/regionalization/data/input' # Required. The save location of standardized output
-  - 'dir_std_base' : '{dir_base}/user_data_std' # Required. The location of standardized data generated by fs_proc python package
+  - 'dir_std_base' : '{dir_base}/user_data_std' # Required. The location of standardized data generated by proc_fs python package
   - 'dir_db_hydfab' : '{dir_base}/hydrofabric' # Required. The local dir where hydrofabric data are stored (limits the total s3 connections)
   - 'dir_db_attrs' : '{dir_base}/attributes' # Required. The parent dir where each comid's attribute parquet file is stored in the subdirectory 'comid/', and each dataset's aggregated parquet attributes are stored in the subdirectory '/{dataset_name}
 formulation_metadata:

pkg/proc.attr.hydfab/inst/extdata/xssa_attr_config_missing_vars.yaml

@@ -16,7 +16,7 @@ file_io: # May define {home_dir} for python's '{home_dir}/string_path'.format(ho
   - 'save_type': 'netcdf' #  Required. Use 'csv' to create a directory structure & save multiple files. May also save as hierarchical files 'netcdf' or 'zarr', or if 'csv' chosen, a directory structure is created
   - 'save_loc': 'local' #  Required.  Use 'local' for saving to a local path via dir_save. Future work will create an approach for 'aws' or other cloud saving methods
   - 'dir_base' : '{home_dir}/noaa/regionalization/data/input' # Required. The save location of standardized output
-  - 'dir_std_base' : '{dir_base}/user_data_std' # Required. The location of standardized data generated by fs_proc python package
+  - 'dir_std_base' : '{dir_base}/user_data_std' # Required. The location of standardized data generated by proc_fs python package
   - 'dir_db_hydfab' : '{dir_base}/hydrofabric' # Required. The local dir where hydrofabric data are stored (limits the total s3 connections)
   - 'dir_db_attrs' : '{dir_base}/attributes' # Required. The parent dir where each comid's attribute parquet file is stored in the subdirectory 'comid/', and each dataset's aggregated parquet attributes are stored in the subdirectory '/{dataset_name}
 formulation_metadata:

pkg/proc.attr.hydfab/tests/testthat/test_proc_attr_grabber.R

@@ -292,7 +292,7 @@ testthat::test_that("proc_attr_usgs_nhd", {
 
 testthat::test_that("proc_attr_hf not a comid",{
  testthat::expect_error(proc.attr.hydfab::proc_attr_hf(comid="13Notacomid14", dir_db_hydfab,
-                                                       custom_name="{lyrs}_",ext = 'gpkg',
+                                                       custom_name="{lyrs}_",fileext = 'gpkg',
                                                        lyrs=c('divides','network')[2],
                                                        hf_cat_sel=TRUE, overwrite=FALSE))
 })