R package interface to the NWSRFS (National Weather Service River Forecast System) Fortran hydrologic models: SAC-SMA, SNOW-17, Unit Hydrograph, Lag-K, Chanloss, and Consuse.
Provides both low-level Fortran wrappers and a high-level orchestration layer that reads NWRFC autocalibration directories, auto-detects model components, and runs the full model chain. Also includes an AdjustQ module for preprocessing upstream flow.
Installation
Using pixi (recommended for development — provides R, gfortran, and test dependencies):
See the top-level README for pixi setup details.
Or from R:
# install.packages("devtools")
devtools::install_github("NOAA-NWRFC/nwsrfs-hydro-models", subdir = "nwsrfs_r")Or from a local clone:
Requires gfortran for compiling the Fortran source.
Quick Start
load_example() runs the full model chain on bundled data and returns an nwsrfs_run object:
library(nwsrfsr)
# NRKW1: Nooksack River — 2 SAC-SMA/SNOW17 zones + 3 Lag-K upstream tribs
run = load_example("NRKW1")
print(run)
#> NWSRFS Model Run
#> Zones: NRKW1-1, NRKW1-2
#> Timestep: 6 hours
#> Length: 62820 timesteps
#> Components:
#> SAC-SMA/SNOW17/UH: TRUE
#> Lag-K: TRUE
#> Chanloss: FALSE
#> Consuse: FALSE
# SFLN2: Salmon Falls Creek — 2 zones + CU zone, chanloss, consumptive use
run2 = load_example("SFLN2")
print(run2)
#> NWSRFS Model Run
#> Zones: SFLN2-1, SFLN2-2, SFLN2-CU
#> Timestep: 6 hours
#> Length: 62828 timesteps
#> Components:
#> SAC-SMA/SNOW17/UH: TRUE
#> Lag-K: FALSE
#> Chanloss: TRUE
#> Consuse: TRUEThe $sim element is the simulated streamflow vector (cfs, 6-hour timestep):
summary(run$sim)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 561.6 1894.5 2888.4 3590.3 4302.7 52279.9
plot(run$sim, type = "l", ylab = "Flow (cfs)", main = "NRKW1 Simulated Flow")You can also run from an NWRFC autocalibration directory on disk via nwsrfs_run():
# example, assuming you have downloaded the sameple data: https://doi.org/10.5281/zenodo.18829935
run = nwsrfs_run("nwsrfs-ac-sample-runs/1zone/FSSO3")Sample autocalibration data for 5 basins is available from the autocalibration repository and archived on Zenodo (DOI: 10.5281/zenodo.18829935).
Exploring Bundled Data
The package ships example data for two stations. The parameter tables have columns p_name (unique ID), name, type (model component), zone, and value:
head(nrkw1_pars)
#> p_name name type zone value
#> 1 init_co_MFNW1 init_co lagk MFNW1 247.4027145
#> 2 init_if_MFNW1 init_if lagk MFNW1 111.5698487
#> 3 init_of_MFNW1 init_of lagk MFNW1 391.9661841
#> 4 init_stor_MFNW1 init_stor lagk MFNW1 113.7500227
#> 5 ktbl_a_MFNW1 ktbl_a lagk MFNW1 -0.7160900
#> 6 ktbl_b_MFNW1 ktbl_b lagk MFNW1 -0.1514212Forcing data is a named list of data.frames, one per zone, with 6-hour timesteps:
names(nrkw1_forcing)
#> [1] "NRKW1-1" "NRKW1-2"
head(nrkw1_forcing[["NRKW1-1"]])
#> year month day hour map_mm mat_degc ptps
#> 1 1979 10 1 0 0 9.278 0
#> 2 1979 10 1 6 0 12.759 0
#> 3 1979 10 1 12 0 17.602 0
#> 4 1979 10 1 18 0 10.296 0
#> 5 1979 10 2 0 0 7.491 0
#> 6 1979 10 2 6 0 10.907 0Upstream flow data (NRKW1 only) is a named list of data.frames for each Lag-K reach:
Running Individual Model Components
The model chain can be run step-by-step using the low-level wrappers. This follows the same sequence the high-level load_example() uses internally: FA -> SAC-SMA/SNOW17 -> UH -> Lag-K.
library(nwsrfsr)
dt_hours = 6
# 1. Apply forcing adjustments (monthly climatological corrections + PET)
forcing_adj = fa_nwrfc(dt_hours, nrkw1_forcing, nrkw1_pars)
head(forcing_adj[["NRKW1-1"]])
#> year month day hour map_mm mat_degc ptps pet_mm etd_mm
#> 1 1979 10 1 0 0 8.063764 0 0.5458438 0.5108825
#> 2 1979 10 1 6 0 11.544764 0 0.5458438 0.5108825
#> 3 1979 10 1 12 0 16.387764 0 0.5458438 0.5108825
#> 4 1979 10 1 18 0 9.081764 0 0.5458438 0.5108825
#> 5 1979 10 2 0 0 6.275488 0 0.6351838 0.5918346
#> 6 1979 10 2 6 0 9.691488 0 0.6351838 0.5918346
# 2. SAC-SMA / SNOW-17 — returns total channel inflow (mm) per zone
tci = sac_snow(dt_hours, forcing_adj, nrkw1_pars)
dim(tci)
#> [1] 62820 2
head(tci)
#> [,1] [,2]
#> [1,] 0.2462325 0.9752563
#> [2,] 0.2488699 0.9984781
#> [3,] 0.2476989 0.9733128
#> [4,] 0.2465621 0.9494683
#> [5,] 0.2454517 0.9268678
#> [6,] 0.2443759 0.9054463
# 3. Unit Hydrograph routing — converts TCI to streamflow (cfs)
local_flow = uh(dt_hours, tci, nrkw1_pars,
sum_zones = TRUE, start_of_timestep = TRUE, backfill = TRUE)
head(local_flow)
#> [1] 485.3301 485.3301 629.9145 661.2546 661.6406 653.3936
# 4. Lag-K routing of upstream tributaries
lagk_flow = lagk(dt_hours, nrkw1_upflow, nrkw1_pars, sum_routes = TRUE)
head(lagk_flow)
#> [1] 887.7466 979.0988 672.8225 674.2631 635.2833 633.9531
# 5. Combine local + routed upstream flow
sim = local_flow + lagk_flow
head(sim)
#> [1] 1373.077 1464.429 1302.737 1335.518 1296.924 1287.347See the API reference for each function: fa_nwrfc(), sac_snow(), uh(), lagk(), chanloss(), consuse()
AdjustQ Preprocessing
AdjustQ replicates the CHPS FEWS AdjustQ tool, creating 6-hour upstream flow time series from daily and instantaneous observations, optionally blended with simulation output.
The quickest way to run it is with adjustq_load_example(), which uses bundled NRKW1 data:
# With simulation (runs load_example("NRKW1") internally)
aq = adjustq_load_example(sim = TRUE)
head(aq)
#> datetime flow_cfs
#> 1 1979-10-01 00:00:00 1089.111
#> 2 1979-10-01 06:00:00 1161.571
#> 3 1979-10-01 12:00:00 1033.318
#> 4 1979-10-01 18:00:00 1059.320
#> 5 1979-10-02 00:00:00 1032.968
#> 6 1979-10-02 06:00:00 1029.565
# Without simulation (observation-only merge)
aq2 = adjustq_load_example(sim = FALSE)To call adjustq() manually with your own data:
run = load_example("NRKW1")
# Build a sim data.frame with year/month/day/hour/flow_cfs columns
sim_df = data.frame(
run$forcings[[1]][, c("year", "month", "day", "hour")],
flow_cfs = run$sim
)
result = adjustq(
daily_flow = nrkw1_daily_flow,
inst_flow = nrkw1_inst_flow,
sim = sim_df
)Bundled Datasets
| Object | Description |
|---|---|
nrkw1_pars |
NRKW1 optimal parameters (366 rows) |
nrkw1_forcing |
NRKW1 6-hour forcing, 2 zones |
nrkw1_upflow |
NRKW1 upstream flow, 3 Lag-K reaches (MFNW1, NFNW1, NSSW1) |
nrkw1_daily_flow |
NRKW1 observed daily streamflow |
nrkw1_inst_flow |
NRKW1 observed 6-hour instantaneous flow |
sfln2_pars |
SFLN2 optimal parameters (514 rows) |
sfln2_forcing |
SFLN2 6-hour forcing, 3 zones (incl. CU) |
sfln2_daily_flow |
SFLN2 observed daily streamflow |
sfln2_inst_flow |
SFLN2 observed 6-hour instantaneous flow |
NRKW1 — Nooksack River at North Cedarville, WA (USGS-12210700). 2-zone SAC-SMA/SNOW17/UH model with 3-reach Lag-K routing, no chanloss or consuse.
SFLN2 — Salmon Falls Creek NR San Jacinto, NV (USGS-13105000). 2-zone SAC-SMA/SNOW17/UH model with chanloss and consumptive use, no upstream flow.
Citation
citation("nwsrfsr")This includes the software citation and the accepted 2026 JAWRA paper reference. Repository-level citation metadata is in /CITATION.cff (repo root).
Testing & Building
pixi run test-r # run R tests (testthat)
pixi run check-r # R CMD check (no manual)
pixi run build-r # build R source tarball
pixi run check-r-cran # CRAN submission check on built tarballOr from R:
devtools::test("nwsrfs_r")Documentation
Build locally with pixi:
pixi run build-docs-r # build pkgdown site
pixi run build-docs-r-readme # render README.qmd to HTML
pixi run build-docs # build all docs (Python + R + R README)Or from R:
pkgdown::build_site(pkg = "nwsrfs_r", install = FALSE, preview = FALSE)