In KevinSee/QRFcapacity: Carrying Capacity with QRF
# setwd('reports')
# load knitr for markdown
library(knitr)
knitr::opts_chunk$set(echo=FALSE,
warning=FALSE,
# error = FALSE,
message=FALSE)
#options(tinytex.verbose = TRUE)
options(knitr.kable.NA = '-')
# options(knitr.table.format = "pandoc")
library(kableExtra)
# load needed packages
library(QRFcapacity)
library(tidyverse)
library(magrittr)
library(sf)
library(janitor)
library(ggpubr)
library(readxl)
# library(ggrepel)
# library(scales)
# library(ggspatial)
# theme_set(theme_bw())
theme_set(theme_pubr(base_size = 10))
knitr::write_bib(c("base",
"survey",
"sf",
"knitr",
"rmarkdown"),
file = 'packages.bib')
# Chinook
data("rch_200")
ira_chnk_list = list("East Fork Salmon" = rch_200 %>%
filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - East Fork Salmon River",
chnk),
"Valley Creek" = rch_200 %>%
filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Valley Creek",
chnk),
"Yankee Fork" = rch_200 %>%
filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Yankee Fork",
chnk),
"North Fork Salmon" = rch_200 %>%
filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - North Fork Salmon River",
GNIS_Name %in% c("Dahlonega Creek",
"Hughes Creek",
# "Moose Creek",
"North Fork Salmon River",
"Sheep Creek",
"Twin Creek"),
chnk),
"Lemhi" = rch_200 %>%
filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Lemhi River",
chnk,
HUC8_code == "17060204"),
"Pahsimeroi" = rch_200 %>%
filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Pahsimeroi River",
chnk),
"Panther Creek" = rch_200 %>%
filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Panther Creek",
GNIS_Name %in% c("Big Deer Creek",
"Blackbird Creek",
"Clear Creek",
"Moyer Creek",
"Musgrove Creek",
"Napias Creek",
"Panther Creek"),
chnk),
"Upper Salmon" = rch_200 %>%
filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Salmon River Upper Mainstem above Redfish Lake",
chnk))
wtsd_chnk_poly = ira_chnk_list %>%
map(.f = function(x) {
x %>%
st_union() %>%
st_convex_hull()
})
# get Chinook range from 200m reaches
chnk_wtsd_range = ira_chnk_list %>%
map(.f = function(x) {
x %>%
mutate(Species = "Chinook Salmon",
SciName = "Oncorhynchus tshawytscha") %>%
select(StreamName = GNIS_Name,
Species, SciName,
UseType = chnk_use,
ESU_DPS = chnk_ESU_DPS,
MPG = chnk_MPG,
NWR_POPID = chnk_NWR_POPID,
NWR_NAME = chnk_NWR_NAME)
})
chnk_domain = do.call(rbind, chnk_wtsd_range)
# Steelhead
ira_sthd_list = list("East Fork Salmon" = rch_200 %>%
filter(sthd_NWR_NAME == "Steelhead (Snake River Basin DPS) - East Fork Salmon River",
sthd),
"Valley Creek" = rch_200 %>%
filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Valley Creek",
sthd),
"Yankee Fork" = rch_200 %>%
filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Yankee Fork",
sthd),
"North Fork Salmon" = rch_200 %>%
filter(sthd_NWR_NAME == "Steelhead (Snake River Basin DPS) - North Fork Salmon River",
GNIS_Name %in% c("Dahlonega Creek",
"Hughes Creek",
# "Moose Creek",
"North Fork Salmon River",
"Sheep Creek",
"Twin Creek"),
sthd),
"Lemhi" = rch_200 %>%
filter(sthd_NWR_NAME == "Steelhead (Snake River Basin DPS) - Lemhi River",
sthd,
HUC8_code == "17060204"),
"Pahsimeroi" = rch_200 %>%
filter(sthd_NWR_NAME == "Steelhead (Snake River Basin DPS) - Pahsimeroi River",
GNIS_Name %in% c("Sulphur Creek",
"Pahsimeroi River",
"Patterson Creek"),
sthd),
"Panther Creek" = rch_200 %>%
filter(sthd_NWR_NAME == "Steelhead (Snake River Basin DPS) - Panther Creek",
GNIS_Name %in% c("Big Deer Creek",
"Blackbird Creek",
"Clear Creek",
"Moyer Creek",
"Musgrove Creek",
"Napias Creek",
"Panther Creek"),
sthd),
"Upper Salmon" = rch_200 %>%
filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Salmon River Upper Mainstem above Redfish Lake",
sthd))
wtsd_sthd_poly = ira_sthd_list %>%
map(.f = function(x) {
x %>%
st_union() %>%
st_convex_hull()
})
# get steelhead range from 200m reaches
sthd_wtsd_range = ira_sthd_list %>%
map(.f = function(x) {
x %>%
mutate(Species = "steelhead",
SciName = "Oncorhynchus mykiss") %>%
select(StreamName = GNIS_Name,
Species, SciName,
UseType = sthd_use,
ESU_DPS = sthd_ESU_DPS,
MPG = sthd_MPG,
NWR_POPID = sthd_NWR_POPID,
NWR_NAME = sthd_NWR_NAME)
})
sthd_domain = do.call(rbind, sthd_wtsd_range)
library(nhdplusTools)
hayden_comid = tibble(site = "Hayden Confluence",
lat = 44.869942,
lon = -113.626182) %>%
st_as_sf(coords = c("lon", "lat"),
crs = 4326) %>%
nhdplusTools::discover_nhdplus_id()
upper_lem = nhdplusTools::plot_nhdplus(outlets = list(hayden_comid),
streamorder = 2,
actually_plot = F)
upper_lem$basin = upper_lem$basin %>%
nngeo::st_remove_holes() %>%
st_transform(st_crs(chnk_domain))
cap_df = c('juv_summer',
# 'juv_summer_dash',
'redds',
"juv_winter") %>%
as.list() %>%
rlang::set_names() %>%
map_df(.id = "lifestage",
.f = function(x) {
rch_cap = st_read(paste0("../output/gpkg/Rch_Cap_RF_", x,
".gpkg"),
quiet = T)
ira_chnk_rch = ira_chnk_list %>%
map_df(.id = "watershed",
.f = function(x) {
tibble(chnk_rch = list(x %>%
inner_join(rch_cap %>%
select(-c(GNIS_Name:Watershed)) %>%
st_drop_geometry() %>%
as_tibble(),
by = "UniqueID")))
})
ira_sthd_rch = ira_sthd_list %>%
map_df(.id = "watershed",
.f = function(x) {
tibble(sthd_rch = list(x %>%
inner_join(rch_cap %>%
select(-c(GNIS_Name:Watershed)) %>%
st_drop_geometry() %>%
as_tibble(),
by = "UniqueID")))
})
if(x == "juv_summer") {
rch_pts = st_read("../output/gpkg/Sum_Juv_Capacity.gpkg",
quiet = T) %>%
st_transform(st_crs(rch_cap))
} else if(x == "juv_summer_dash") {
rch_pts = st_read("../output/gpkg/Sum_Juv_Capacity_DASH.gpkg",
quiet = T) %>%
st_transform(st_crs(rch_cap))
} else if(x == "redds") {
rch_pts = st_read("../output/gpkg/Redds_Capacity.gpkg",
quiet = T) %>%
st_transform(st_crs(rch_cap))
} else if(x == "juv_winter") {
rch_pts = st_read("../output/gpkg/Win_Juv_Capacity.gpkg",
quiet = T) %>%
st_transform(st_crs(rch_cap))
}
ira_chnk_pts = wtsd_chnk_poly %>%
map_df(.id = "watershed",
.f = function(x) {
tibble(chnk_pts = list(st_intersection(rch_pts, x)))
})
ira_sthd_pts = wtsd_sthd_poly %>%
map_df(.id = "watershed",
.f = function(x) {
tibble(sthd_pts = list(st_intersection(rch_pts, x)))
})
ira_chnk_rch %>%
full_join(ira_sthd_rch,
by = "watershed") %>%
full_join(ira_chnk_pts,
by = "watershed") %>%
full_join(ira_sthd_pts,
by = "watershed") %>%
return()
}) %>%
left_join(wtsd_chnk_poly %>%
map_df(.id = "watershed",
.f = function(x) {
tibble(wtsd_chnk_poly = list(x))
}),
by = "watershed") %>%
left_join(wtsd_sthd_poly %>%
map_df(.id = "watershed",
.f = function(x) {
tibble(wtsd_sthd_poly = list(x))
}),
by = "watershed")
cap_df %<>%
pivot_longer(cols = -(lifestage:watershed),
names_to = "name",
values_to = "value") %>%
mutate(Species = if_else(grepl('chnk', name), "Chinook",
if_else(grepl('sthd', name), "Steelhead",
NA_character_))) %>%
mutate(name = str_remove(name, "chnk_"),
name = str_remove(name, "sthd_")) %>%
pivot_wider(names_from = "name",
values_from = "value") %>%
left_join(tibble(Species = c("Chinook", "Steelhead"),
spp_domain = list(chnk_domain,
sthd_domain)),
by = "Species")
cap_est = cap_df %>%
nest(inputs = rch:spp_domain) %>%
mutate(cap_rch = map2(inputs,
Species,
.f = function(y,z) {
calc_watershed_cap(wtsd_polygon = y$wtsd_poly[[1]],
capacity_sf = y$rch[[1]],
capacity_name = if_else(z == "Chinook",
"chnk_per_m",
"sthd_per_m"),
capacity_se_name = if_else(z == "Chinook",
"chnk_per_m_se",
"sthd_per_m_se"),
# spp_range = y$spp_domain[[1]],
by_stream = F) %>%
rename(n = n_rchs)
})) %>%
mutate(cap_pts = map2(inputs,
Species,
.f = function(y,z) {
calc_watershed_cap(wtsd_polygon = y$wtsd_poly[[1]],
capacity_sf = y$pts[[1]],
capacity_name = if_else(z == "Chinook",
"chnk_per_m",
"sthd_per_m"),
capacity_se_name = if_else(z == "Chinook",
"chnk_per_m_se",
"sthd_per_m_se"),
spp_range = y$spp_domain[[1]],
by_stream = F) %>%
select(-area) %>%
rename(n = n_pts)
})) %>%
pivot_longer(cols = c(cap_rch,
cap_pts),
names_to = "source",
values_to = "cap") %>%
mutate(source = recode(source,
"cap_rch" = "Reaches",
"cap_pts" = "Points")) %>%
unnest(cols = cap) %>%
mutate(across(tot_length,
~ . / 1000))
cap_comp_df = cap_est %>%
select(-inputs)
params = read_excel('../../MRA_QRF/data/IRA/USIRA Tables (all Upper Salmon).xlsx',
"Escapement Data - Chinook",
range = "B2:C7") %>%
add_column(Species = 'Chinook',
.before = 0) %>%
bind_rows(read_excel('../../MRA_QRF/data/IRA/USIRA Tables (all Upper Salmon).xlsx',
"Escapement Data - Steelhead",
range = "B2:C7") %>%
add_column(Species = 'Steelhead',
.before = 0))
# set recovery goals
rec_goals = tibble(Watershed = c("Upper Salmon",
"Valley Creek",
"Yankee Fork",
"East Fork Salmon",
"Pahsimeroi",
"Lemhi",
"North Fork Salmon",
"Panther Creek"),
Species = c("Chinook"),
Mean = c(656,
506,
248,
283,
402,
347,
NA,
NA),
Max = c(1419,
739,
343,
343,
822,
718,
NA,
NA),
MAT = c(1000,
500,
500,
1000,
1000,
2000,
500,
750)) %>%
bind_rows(tibble(Watershed = c("Upper Salmon",
"Valley Creek",
"Yankee Fork",
"East Fork Salmon",
"Pahsimeroi",
"Lemhi",
"North Fork Salmon",
"Panther Creek"),
Species = c("Steelhead"),
Mean = c(92,
193,
95,
30,
1156,
337,
252,
449),
Max = c(154,
278,
213,
54,
1614,
417,
349,
650),
MAT = c(483,
247,
270,
1000,
1000,
1000,
500,
500))) %>%
mutate(`MAT+25%` = MAT * 1.25) %>%
pivot_longer(cols = -c(Watershed:Species),
names_to = "Scenario",
values_to = "Escapement") %>%
inner_join(params %>%
mutate(Parameter = recode(Parameter,
"Female Ratio" = "prop_fem",
"Redds/Female" = "redd_per_fem",
"Fecundity" = 'fecud',
'Egg:Parr' = "egg_to_parr",
"Parr:Presmolt" = "parr_to_presmolt")) %>%
pivot_wider(id_cols = "Species",
names_from = "Parameter",
values_from = "Value"),
by = "Species") %>%
mutate(Redds = Escapement * prop_fem * redd_per_fem,
Eggs = Redds * fecud,
Parr = Eggs * egg_to_parr,
Presmolts = Parr * parr_to_presmolt) %>%
select(-c(prop_fem:parr_to_presmolt)) %>%
pivot_longer(cols = c(Escapement, Redds:Presmolts),
names_to = "Lifestage",
values_to = "Abundance")
comp_df = cap_comp_df %>%
rename(Lifestage = lifestage,
Watershed = watershed) %>%
mutate(Lifestage = recode(Lifestage,
"juv_summer" = "Parr",
"redds" = "Redds",
"juv_winter" = "Presmolts")) %>%
mutate(Scenario = paste0("QRF_cap_", source)) %>%
rename(Abundance = tot_cap) %>%
select(any_of(names(rec_goals))) %>%
bind_rows(rec_goals) %>%
pivot_wider(names_from = "Scenario",
values_from = "Abundance") %>%
select(Species, Watershed, Lifestage, Mean:`MAT+25%`, starts_with("QRF")) %>%
# mutate(cap_def = `MAT+25%` - QRF_cap,
# rel_def = cap_def / QRF_cap) %>%
mutate(Lifestage = factor(Lifestage,
levels = c("Redds",
"Parr",
"Presmolts"))) %>%
arrange(Watershed, Lifestage) %>%
filter(!is.na(Lifestage)) %>%
pivot_longer(cols = Mean:QRF_cap_Points,
names_to = "scenario",
values_to = "value") %>%
mutate(scenario = factor(scenario,
levels = c("QRF_cap_Reaches",
"QRF_cap_Points",
"Mean",
"Max",
"MAT",
"MAT+25%")),
scenario = fct_rev(scenario))
bar_list_p = comp_df %>%
split(list(.$Species, .$Lifestage)) %>%
map(.f = function(df) {
df %>%
ggplot(aes(y = scenario,
x = value,
fill = scenario)) +
geom_col(position = "dodge") +
# scale_fill_brewer(palette = "Paired") +
scale_fill_manual(values = c("QRF_cap_Reaches" = "black",
"QRF_cap_Points" = "gray",
"Mean" = "lightgreen",
"Max" = "darkgreen",
"MAT" = "lightblue",
"MAT+25%" = "darkblue"),
name = "Abundance") +
facet_wrap(~ Watershed,
scales = "free_x") +
theme(legend.position = "bottom",
axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank()) +
labs(x = "Scenario",
title = unique(df$Lifestage))
})
Introduction
When the IRA report came out, our QRF model included an extrapolation model based on GRTS master sample points [@See2021]. Since then, we have updated the extrapolation model to one that is based on a linear network of 200m reaches. This linear network is also based on an updated NHDPlus network, which often increases the length of various streams as the updated version has finer resolution.
To examine the impact of this change in extrapolation models, we re-compiled capacity estimates at the watershed scale for all IRA watersheds, across 3 lifestages and both species. We computed capacity using the master sample points, and also using the 200m reaches. We used the updated NHDPlus network for both, so the total stream kilometers for each watershed is the same. For some small tributaries, there may not be a master sample point on that tributary within the range of that species. In such cases, we have no estimate of the carrying capacity of that tributary, so we set it to zero.
Methods
We used the same life-history parameters used in the IRA addendum, shown in Table \@ref(tab:param-table).
params %>%
kable(caption = "Life history parameters used to translate recovery goals into various life-stages.",
digits = 2) %>%
kable_styling()
Results
Figure \@ref(fig:all-comp-fig) shows the comparison between estimates from points and reaches, along with the 95% confidence intervals, broken down by species and lifestage. These results are shown in Table \@ref(tab:chnk-cap-tab) for Chinook, and Table \@ref(tab:sthd-cap-tab) for steelhead.
comp_p_list = cap_comp_df %>%
rename(Lifestage = lifestage,
Watershed = watershed,
Length = tot_length,
Capacity = tot_cap,
SE = tot_cap_se) %>%
mutate(Lifestage = recode(Lifestage,
"juv_summer" = "Parr",
"redds" = "Redds",
"juv_winter" = "Presmolts")) %>%
mutate(Lifestage = factor(Lifestage,
levels = c("Parr",
"Presmolts",
"Redds"))) %>%
select(-n) %>%
pivot_wider(names_from = "source",
values_from = c("Length", "Capacity", "SE"),
names_glue = "{source}_{.value}") %>%
split(list(.$Species, .$Lifestage)) %>%
map(.f = function(z) {
my_range = c(range(z$Reaches_Capacity),
range(z$Points_Capacity))
z %>%
ggplot(aes(x = Points_Capacity,
y = Reaches_Capacity,
color = Watershed)) +
geom_abline(linetype = 2) +
geom_errorbar(aes(ymin = Reaches_Capacity + Reaches_SE * qnorm(0.025),
ymax = Reaches_Capacity + Reaches_SE * qnorm(0.975)),
width = 0) +
geom_errorbarh(aes(xmin = Points_Capacity + Points_SE * qnorm(0.025),
xmax = Points_Capacity + Points_SE * qnorm(0.975)),
height = 0) +
geom_point(size = 3) +
coord_cartesian(xlim = range(my_range),
ylim = range(my_range)) +
# scale_x_continuous(limits = range(my_range)) +
# scale_y_continuous(limits = range(my_range)) +
labs(x = "Points",
y = "Reaches",
title = paste(unique(z$Species), unique(z$Lifestage), "Capacity"))
})
comp_p = ggarrange(plotlist = comp_p_list,
nrow = 3,
ncol = 2,
common.legend = T,
legend = "bottom")
comp_p
The contemporary mean and max of recent escapements (translated to other lifestages), as well as the recovery goal of mean abundance threshold (MAT) and a 25% buffer on top of that (MAT+25%) are shown in Figures \@ref(fig:chnk-parr-fig) - \@ref(fig:sthd-redd-fig), together with estimates of current capacity from master sample points and reaches.
Chinook
bar_list_p$Chinook.Parr
bar_list_p$Chinook.Presmolts
bar_list_p$Chinook.Redds
cap_comp_df %>%
filter(Species == "Chinook") %>%
rename(Lifestage = lifestage,
Watershed = watershed) %>%
mutate(Lifestage = recode(Lifestage,
"juv_summer" = "Parr",
"redds" = "Redds",
"juv_winter" = "Presmolts")) %>%
mutate(Lifestage = factor(Lifestage,
levels = c("Redds",
"Parr",
"Presmolts"))) %>%
mutate(tot_cap_cv = tot_cap_se / tot_cap,
Capacity = paste0(prettyNum(round(tot_cap), big.mark = ","), " (", round(tot_cap_cv, 2), ")")) %>%
select(Lifestage:source, tot_length, Capacity) %>%
rename(Length = tot_length) %>%
pivot_wider(names_from = "source",
values_from = c("Length", "Capacity"),
names_glue = "{source}_{.value}") %>%
kable(digits = 1,
caption = "Estimated capacities (CVs) for Chinook in various lifestages using the master sample points and linear network reaches for extrapolation. Total watershed length for each method is shown as well.") %>%
kable_styling(fixed_thead = T)
Steelhead
bar_list_p$Steelhead.Parr
bar_list_p$Steelhead.Presmolts
bar_list_p$Steelhead.Redds
cap_comp_df %>%
filter(Species == "Steelhead") %>%
rename(Lifestage = lifestage,
Watershed = watershed) %>%
mutate(Lifestage = recode(Lifestage,
"juv_summer" = "Parr",
"redds" = "Redds",
"juv_winter" = "Presmolts")) %>%
mutate(Lifestage = factor(Lifestage,
levels = c("Redds",
"Parr",
"Presmolts"))) %>%
mutate(tot_cap_cv = tot_cap_se / tot_cap,
Capacity = paste0(prettyNum(round(tot_cap), big.mark = ","), " (", round(tot_cap_cv, 2), ")")) %>%
select(Lifestage:source, tot_length, Capacity) %>%
rename(Length = tot_length) %>%
pivot_wider(names_from = "source",
values_from = c("Length", "Capacity"),
names_glue = "{source}_{.value}") %>%
kable(digits = 1,
caption = "Estimated capacities (CVs) for steelhead in various lifestages using the master sample points and linear network reaches for extrapolation. Total watershed length for each method is shown as well.") %>%
kable_styling(fixed_thead = T)
Lemhi Specific
We broke down the Lemhi available capacity into the Upper and Lower Lemhi (defined by the confluence of Hayden Creek).
upper_lem_poly = upper_lem$basin
lower_lem_poly = cap_est %>%
filter(watershed == "Lemhi") %>%
slice(1) %>%
pull(inputs) %>%
map("wtsd_poly") %>%
extract2(1) %>%
extract2(1) %>%
st_difference(upper_lem_poly)
lem_cap = cap_est %>%
filter(watershed == "Lemhi") %>%
select(lifestage:source) %>%
mutate(upp_lem_cap = map(inputs,
.f = function(y) {
calc_watershed_cap(wtsd_polygon = upper_lem_poly,
capacity_sf = y$rch[[1]],
spp_range = y$spp_domain[[1]],
by_stream = F) %>%
rename(n = n_rchs)
})) %>%
mutate(low_lem_cap = map(inputs,
.f = function(y) {
calc_watershed_cap(wtsd_polygon = lower_lem_poly,
capacity_sf = y$rch[[1]],
spp_range = y$spp_domain[[1]],
by_stream = F) %>%
rename(n = n_rchs)
})) %>%
pivot_longer(cols = ends_with("lem_cap"),
names_to = "area",
values_to = "cap") %>%
unnest(cols = cap) %>%
mutate(area = recode(area,
"upp_lem_cap" = "Upper Lemhi",
"low_lem_cap" = "Lower Lemhi")) %>%
select(-inputs) %>%
mutate(across(tot_length,
~. / 1000))
lem_comp_df = lem_cap %>%
rename(Watershed = watershed,
Lifestage = lifestage) %>%
mutate(Lifestage = recode(Lifestage,
"juv_summer" = "Parr",
"redds" = "Redds",
"juv_winter" = "Presmolts")) %>%
mutate(Scenario = paste0("QRF_cap_", source)) %>%
select(Lifestage:area,
Scenario,
Abundance = tot_cap) %>%
filter(source == "Points") %>%
bind_rows(rec_goals %>%
filter(Watershed == "Lemhi",
Lifestage %in% c("Parr", "Redds", "Presmolts"))) %>%
mutate(Lifestage = factor(Lifestage,
levels = c("Redds",
"Parr",
"Presmolts")),
Scenario = factor(Scenario,
levels = c("QRF_cap_Points",
"Mean",
"Max",
"MAT",
"MAT+25%"))) %>%
mutate(area = if_else(is.na(area), as.character(Scenario), area)) %>%
filter(Scenario %in% c("QRF_cap_Points",
"MAT+25%")) %>%
mutate(Scenario = fct_expand(Scenario, "Available Capacity (QRF)", "Required Capacity\n(MAT +25%)"),
Scenario = fct_recode(Scenario,
"Available Capacity (QRF)" = "QRF_cap_Points",
"Required Capacity\n(MAT +25%)" = "MAT+25%"),
Scenario = fct_rev(Scenario)) %>%
mutate(area = factor(area,
levels = c("MAT+25%",
"Upper Lemhi",
"Lower Lemhi"))) %>%
arrange(Species, Lifestage, area, Scenario)
lem_bar_p = lem_comp_df %>%
split(list(.$Species, .$Lifestage)) %>%
map(.f = function(df) {
df %>%
ggplot(aes(y = Scenario,
x = Abundance,
fill = area)) +
geom_col(position = "stack") +
scale_fill_manual(values = c("Lower Lemhi" = "darkblue",
"Upper Lemhi" = "lightblue",
"MAT+25%" = "black"),
breaks = c("MAT+25%", "Lower Lemhi", "Upper Lemhi"),
name = NULL) +
# scale_fill_brewer(palette = "Paired",
# name = "Scenario") +
theme(legend.position = "bottom",
# axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank()) +
scale_x_continuous(labels = function(x) format(x, scientific = F)) +
labs(x = "Abundance",
title = paste(unique(df$Species), unique(df$Lifestage)))
})
Chinook
lem_bar_p$Chinook.Parr
lem_bar_p$Chinook.Presmolts
lem_bar_p$Chinook.Redds
Steelhead
lem_bar_p$Steelhead.Parr
lem_bar_p$Steelhead.Presmolts
lem_bar_p$Steelhead.Redds
walk(names(lem_bar_p), function(.x) {
ggsave(
path = "figures/",
filename = paste0(.x, ".png"),
plot = lem_bar_p[[.x]],
width = 6.5,
height = 2
)
})
Discussion
In general, we noted the following patterns in comparing QRF estimates of capacity using master sample points versus a linear network:
- Chinook parr usually had greater capacity estimates with linear network (with a couple exceptions)
- Chinook presmolts had similar capacity estimates (some smaller, some larger)
- Chinook redds had much smaller capacity estimates with linear network
- Steelhead parr had similar or slightly greater capacity estimates with linear network
- Steelhead presmolts had very similar capacity estimates
- Steelhead redds had much smaller capacity estimates with linear network
We do not believe there is anything inherently incorrect in either extrapolation method. Differences in estimates may arise due to differences in available covariates between master sample points and 200m reaches (and therefore differences in the extrapolation models), but also because the master sample point version must "smooth" across a linear network from a series of points (i.e. taking the average fish/m capacity from points on a stream, and multiplying ny the length of that stream), whereas the 200m reaches are already on a linear network, and we are making extrapolation estimates at every single reach. For this reason, we would recommend using the results from the 200m reach extrapolation going forward.
One of the goals of this exercise was to make as close to apple-to-apples comparisons between estimates based on master sample points and estimates based on 200m reaches on a linear network. It needs to be noted that there are a few streams in the upper Salmon area that are not included in the 200m reach linear network (e.g. Big Springs, Little Springs, Pratt Creek and Lee Creek in the Lemhi, Patterson Side Channels in the Pahsimeroi, Yellowbelly Creek in the Upper Salmon, etc.). Therefore, the estimates for the master sample points presented here may not match those in the IRA report exactly. But because of that, the lengths in Tables \@ref(tab:chnk-cap-tab) and \@ref(tab:sthd-cap-tab) are nearly identical, so any differences in capacity estimates here are not simply due to different species ranges being used; we tried to put them on a level playing field in that regard. This does mean that the linear network could be improved, especially in the Lemhi (and less so in the Pahsimeroi) to allow the 200m reach extrapolation to include those missing streams. That will involve generating not just shapefiles of the stream lines, but either calculating directly or finding a proxy for the various covariates that need to be attached to each reach for the extrapolation model.
References
KevinSee/QRFcapacity documentation built on Feb. 27, 2023, 3:57 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# setwd('reports') # load knitr for markdown library(knitr) knitr::opts_chunk$set(echo=FALSE, warning=FALSE, # error = FALSE, message=FALSE) #options(tinytex.verbose = TRUE) options(knitr.kable.NA = '-') # options(knitr.table.format = "pandoc") library(kableExtra)
# load needed packages library(QRFcapacity) library(tidyverse) library(magrittr) library(sf) library(janitor) library(ggpubr) library(readxl) # library(ggrepel) # library(scales) # library(ggspatial) # theme_set(theme_bw()) theme_set(theme_pubr(base_size = 10))
knitr::write_bib(c("base", "survey", "sf", "knitr", "rmarkdown"), file = 'packages.bib')
# Chinook data("rch_200") ira_chnk_list = list("East Fork Salmon" = rch_200 %>% filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - East Fork Salmon River", chnk), "Valley Creek" = rch_200 %>% filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Valley Creek", chnk), "Yankee Fork" = rch_200 %>% filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Yankee Fork", chnk), "North Fork Salmon" = rch_200 %>% filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - North Fork Salmon River", GNIS_Name %in% c("Dahlonega Creek", "Hughes Creek", # "Moose Creek", "North Fork Salmon River", "Sheep Creek", "Twin Creek"), chnk), "Lemhi" = rch_200 %>% filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Lemhi River", chnk, HUC8_code == "17060204"), "Pahsimeroi" = rch_200 %>% filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Pahsimeroi River", chnk), "Panther Creek" = rch_200 %>% filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Panther Creek", GNIS_Name %in% c("Big Deer Creek", "Blackbird Creek", "Clear Creek", "Moyer Creek", "Musgrove Creek", "Napias Creek", "Panther Creek"), chnk), "Upper Salmon" = rch_200 %>% filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Salmon River Upper Mainstem above Redfish Lake", chnk)) wtsd_chnk_poly = ira_chnk_list %>% map(.f = function(x) { x %>% st_union() %>% st_convex_hull() }) # get Chinook range from 200m reaches chnk_wtsd_range = ira_chnk_list %>% map(.f = function(x) { x %>% mutate(Species = "Chinook Salmon", SciName = "Oncorhynchus tshawytscha") %>% select(StreamName = GNIS_Name, Species, SciName, UseType = chnk_use, ESU_DPS = chnk_ESU_DPS, MPG = chnk_MPG, NWR_POPID = chnk_NWR_POPID, NWR_NAME = chnk_NWR_NAME) }) chnk_domain = do.call(rbind, chnk_wtsd_range) # Steelhead ira_sthd_list = list("East Fork Salmon" = rch_200 %>% filter(sthd_NWR_NAME == "Steelhead (Snake River Basin DPS) - East Fork Salmon River", sthd), "Valley Creek" = rch_200 %>% filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Valley Creek", sthd), "Yankee Fork" = rch_200 %>% filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Yankee Fork", sthd), "North Fork Salmon" = rch_200 %>% filter(sthd_NWR_NAME == "Steelhead (Snake River Basin DPS) - North Fork Salmon River", GNIS_Name %in% c("Dahlonega Creek", "Hughes Creek", # "Moose Creek", "North Fork Salmon River", "Sheep Creek", "Twin Creek"), sthd), "Lemhi" = rch_200 %>% filter(sthd_NWR_NAME == "Steelhead (Snake River Basin DPS) - Lemhi River", sthd, HUC8_code == "17060204"), "Pahsimeroi" = rch_200 %>% filter(sthd_NWR_NAME == "Steelhead (Snake River Basin DPS) - Pahsimeroi River", GNIS_Name %in% c("Sulphur Creek", "Pahsimeroi River", "Patterson Creek"), sthd), "Panther Creek" = rch_200 %>% filter(sthd_NWR_NAME == "Steelhead (Snake River Basin DPS) - Panther Creek", GNIS_Name %in% c("Big Deer Creek", "Blackbird Creek", "Clear Creek", "Moyer Creek", "Musgrove Creek", "Napias Creek", "Panther Creek"), sthd), "Upper Salmon" = rch_200 %>% filter(chnk_NWR_NAME == "Chinook Salmon (Snake River Spring/Summer-run ESU) - Salmon River Upper Mainstem above Redfish Lake", sthd)) wtsd_sthd_poly = ira_sthd_list %>% map(.f = function(x) { x %>% st_union() %>% st_convex_hull() }) # get steelhead range from 200m reaches sthd_wtsd_range = ira_sthd_list %>% map(.f = function(x) { x %>% mutate(Species = "steelhead", SciName = "Oncorhynchus mykiss") %>% select(StreamName = GNIS_Name, Species, SciName, UseType = sthd_use, ESU_DPS = sthd_ESU_DPS, MPG = sthd_MPG, NWR_POPID = sthd_NWR_POPID, NWR_NAME = sthd_NWR_NAME) }) sthd_domain = do.call(rbind, sthd_wtsd_range)
library(nhdplusTools) hayden_comid = tibble(site = "Hayden Confluence", lat = 44.869942, lon = -113.626182) %>% st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% nhdplusTools::discover_nhdplus_id() upper_lem = nhdplusTools::plot_nhdplus(outlets = list(hayden_comid), streamorder = 2, actually_plot = F) upper_lem$basin = upper_lem$basin %>% nngeo::st_remove_holes() %>% st_transform(st_crs(chnk_domain))
cap_df = c('juv_summer', # 'juv_summer_dash', 'redds', "juv_winter") %>% as.list() %>% rlang::set_names() %>% map_df(.id = "lifestage", .f = function(x) { rch_cap = st_read(paste0("../output/gpkg/Rch_Cap_RF_", x, ".gpkg"), quiet = T) ira_chnk_rch = ira_chnk_list %>% map_df(.id = "watershed", .f = function(x) { tibble(chnk_rch = list(x %>% inner_join(rch_cap %>% select(-c(GNIS_Name:Watershed)) %>% st_drop_geometry() %>% as_tibble(), by = "UniqueID"))) }) ira_sthd_rch = ira_sthd_list %>% map_df(.id = "watershed", .f = function(x) { tibble(sthd_rch = list(x %>% inner_join(rch_cap %>% select(-c(GNIS_Name:Watershed)) %>% st_drop_geometry() %>% as_tibble(), by = "UniqueID"))) }) if(x == "juv_summer") { rch_pts = st_read("../output/gpkg/Sum_Juv_Capacity.gpkg", quiet = T) %>% st_transform(st_crs(rch_cap)) } else if(x == "juv_summer_dash") { rch_pts = st_read("../output/gpkg/Sum_Juv_Capacity_DASH.gpkg", quiet = T) %>% st_transform(st_crs(rch_cap)) } else if(x == "redds") { rch_pts = st_read("../output/gpkg/Redds_Capacity.gpkg", quiet = T) %>% st_transform(st_crs(rch_cap)) } else if(x == "juv_winter") { rch_pts = st_read("../output/gpkg/Win_Juv_Capacity.gpkg", quiet = T) %>% st_transform(st_crs(rch_cap)) } ira_chnk_pts = wtsd_chnk_poly %>% map_df(.id = "watershed", .f = function(x) { tibble(chnk_pts = list(st_intersection(rch_pts, x))) }) ira_sthd_pts = wtsd_sthd_poly %>% map_df(.id = "watershed", .f = function(x) { tibble(sthd_pts = list(st_intersection(rch_pts, x))) }) ira_chnk_rch %>% full_join(ira_sthd_rch, by = "watershed") %>% full_join(ira_chnk_pts, by = "watershed") %>% full_join(ira_sthd_pts, by = "watershed") %>% return() }) %>% left_join(wtsd_chnk_poly %>% map_df(.id = "watershed", .f = function(x) { tibble(wtsd_chnk_poly = list(x)) }), by = "watershed") %>% left_join(wtsd_sthd_poly %>% map_df(.id = "watershed", .f = function(x) { tibble(wtsd_sthd_poly = list(x)) }), by = "watershed") cap_df %<>% pivot_longer(cols = -(lifestage:watershed), names_to = "name", values_to = "value") %>% mutate(Species = if_else(grepl('chnk', name), "Chinook", if_else(grepl('sthd', name), "Steelhead", NA_character_))) %>% mutate(name = str_remove(name, "chnk_"), name = str_remove(name, "sthd_")) %>% pivot_wider(names_from = "name", values_from = "value") %>% left_join(tibble(Species = c("Chinook", "Steelhead"), spp_domain = list(chnk_domain, sthd_domain)), by = "Species")
cap_est = cap_df %>% nest(inputs = rch:spp_domain) %>% mutate(cap_rch = map2(inputs, Species, .f = function(y,z) { calc_watershed_cap(wtsd_polygon = y$wtsd_poly[[1]], capacity_sf = y$rch[[1]], capacity_name = if_else(z == "Chinook", "chnk_per_m", "sthd_per_m"), capacity_se_name = if_else(z == "Chinook", "chnk_per_m_se", "sthd_per_m_se"), # spp_range = y$spp_domain[[1]], by_stream = F) %>% rename(n = n_rchs) })) %>% mutate(cap_pts = map2(inputs, Species, .f = function(y,z) { calc_watershed_cap(wtsd_polygon = y$wtsd_poly[[1]], capacity_sf = y$pts[[1]], capacity_name = if_else(z == "Chinook", "chnk_per_m", "sthd_per_m"), capacity_se_name = if_else(z == "Chinook", "chnk_per_m_se", "sthd_per_m_se"), spp_range = y$spp_domain[[1]], by_stream = F) %>% select(-area) %>% rename(n = n_pts) })) %>% pivot_longer(cols = c(cap_rch, cap_pts), names_to = "source", values_to = "cap") %>% mutate(source = recode(source, "cap_rch" = "Reaches", "cap_pts" = "Points")) %>% unnest(cols = cap) %>% mutate(across(tot_length, ~ . / 1000)) cap_comp_df = cap_est %>% select(-inputs)
params = read_excel('../../MRA_QRF/data/IRA/USIRA Tables (all Upper Salmon).xlsx', "Escapement Data - Chinook", range = "B2:C7") %>% add_column(Species = 'Chinook', .before = 0) %>% bind_rows(read_excel('../../MRA_QRF/data/IRA/USIRA Tables (all Upper Salmon).xlsx', "Escapement Data - Steelhead", range = "B2:C7") %>% add_column(Species = 'Steelhead', .before = 0)) # set recovery goals rec_goals = tibble(Watershed = c("Upper Salmon", "Valley Creek", "Yankee Fork", "East Fork Salmon", "Pahsimeroi", "Lemhi", "North Fork Salmon", "Panther Creek"), Species = c("Chinook"), Mean = c(656, 506, 248, 283, 402, 347, NA, NA), Max = c(1419, 739, 343, 343, 822, 718, NA, NA), MAT = c(1000, 500, 500, 1000, 1000, 2000, 500, 750)) %>% bind_rows(tibble(Watershed = c("Upper Salmon", "Valley Creek", "Yankee Fork", "East Fork Salmon", "Pahsimeroi", "Lemhi", "North Fork Salmon", "Panther Creek"), Species = c("Steelhead"), Mean = c(92, 193, 95, 30, 1156, 337, 252, 449), Max = c(154, 278, 213, 54, 1614, 417, 349, 650), MAT = c(483, 247, 270, 1000, 1000, 1000, 500, 500))) %>% mutate(`MAT+25%` = MAT * 1.25) %>% pivot_longer(cols = -c(Watershed:Species), names_to = "Scenario", values_to = "Escapement") %>% inner_join(params %>% mutate(Parameter = recode(Parameter, "Female Ratio" = "prop_fem", "Redds/Female" = "redd_per_fem", "Fecundity" = 'fecud', 'Egg:Parr' = "egg_to_parr", "Parr:Presmolt" = "parr_to_presmolt")) %>% pivot_wider(id_cols = "Species", names_from = "Parameter", values_from = "Value"), by = "Species") %>% mutate(Redds = Escapement * prop_fem * redd_per_fem, Eggs = Redds * fecud, Parr = Eggs * egg_to_parr, Presmolts = Parr * parr_to_presmolt) %>% select(-c(prop_fem:parr_to_presmolt)) %>% pivot_longer(cols = c(Escapement, Redds:Presmolts), names_to = "Lifestage", values_to = "Abundance")
comp_df = cap_comp_df %>% rename(Lifestage = lifestage, Watershed = watershed) %>% mutate(Lifestage = recode(Lifestage, "juv_summer" = "Parr", "redds" = "Redds", "juv_winter" = "Presmolts")) %>% mutate(Scenario = paste0("QRF_cap_", source)) %>% rename(Abundance = tot_cap) %>% select(any_of(names(rec_goals))) %>% bind_rows(rec_goals) %>% pivot_wider(names_from = "Scenario", values_from = "Abundance") %>% select(Species, Watershed, Lifestage, Mean:`MAT+25%`, starts_with("QRF")) %>% # mutate(cap_def = `MAT+25%` - QRF_cap, # rel_def = cap_def / QRF_cap) %>% mutate(Lifestage = factor(Lifestage, levels = c("Redds", "Parr", "Presmolts"))) %>% arrange(Watershed, Lifestage) %>% filter(!is.na(Lifestage)) %>% pivot_longer(cols = Mean:QRF_cap_Points, names_to = "scenario", values_to = "value") %>% mutate(scenario = factor(scenario, levels = c("QRF_cap_Reaches", "QRF_cap_Points", "Mean", "Max", "MAT", "MAT+25%")), scenario = fct_rev(scenario)) bar_list_p = comp_df %>% split(list(.$Species, .$Lifestage)) %>% map(.f = function(df) { df %>% ggplot(aes(y = scenario, x = value, fill = scenario)) + geom_col(position = "dodge") + # scale_fill_brewer(palette = "Paired") + scale_fill_manual(values = c("QRF_cap_Reaches" = "black", "QRF_cap_Points" = "gray", "Mean" = "lightgreen", "Max" = "darkgreen", "MAT" = "lightblue", "MAT+25%" = "darkblue"), name = "Abundance") + facet_wrap(~ Watershed, scales = "free_x") + theme(legend.position = "bottom", axis.text.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_blank()) + labs(x = "Scenario", title = unique(df$Lifestage)) })
Introduction
When the IRA report came out, our QRF model included an extrapolation model based on GRTS master sample points [@See2021]. Since then, we have updated the extrapolation model to one that is based on a linear network of 200m reaches. This linear network is also based on an updated NHDPlus network, which often increases the length of various streams as the updated version has finer resolution.
To examine the impact of this change in extrapolation models, we re-compiled capacity estimates at the watershed scale for all IRA watersheds, across 3 lifestages and both species. We computed capacity using the master sample points, and also using the 200m reaches. We used the updated NHDPlus network for both, so the total stream kilometers for each watershed is the same. For some small tributaries, there may not be a master sample point on that tributary within the range of that species. In such cases, we have no estimate of the carrying capacity of that tributary, so we set it to zero.
Methods
We used the same life-history parameters used in the IRA addendum, shown in Table \@ref(tab:param-table).
params %>% kable(caption = "Life history parameters used to translate recovery goals into various life-stages.", digits = 2) %>% kable_styling()
Results
Figure \@ref(fig:all-comp-fig) shows the comparison between estimates from points and reaches, along with the 95% confidence intervals, broken down by species and lifestage. These results are shown in Table \@ref(tab:chnk-cap-tab) for Chinook, and Table \@ref(tab:sthd-cap-tab) for steelhead.
comp_p_list = cap_comp_df %>% rename(Lifestage = lifestage, Watershed = watershed, Length = tot_length, Capacity = tot_cap, SE = tot_cap_se) %>% mutate(Lifestage = recode(Lifestage, "juv_summer" = "Parr", "redds" = "Redds", "juv_winter" = "Presmolts")) %>% mutate(Lifestage = factor(Lifestage, levels = c("Parr", "Presmolts", "Redds"))) %>% select(-n) %>% pivot_wider(names_from = "source", values_from = c("Length", "Capacity", "SE"), names_glue = "{source}_{.value}") %>% split(list(.$Species, .$Lifestage)) %>% map(.f = function(z) { my_range = c(range(z$Reaches_Capacity), range(z$Points_Capacity)) z %>% ggplot(aes(x = Points_Capacity, y = Reaches_Capacity, color = Watershed)) + geom_abline(linetype = 2) + geom_errorbar(aes(ymin = Reaches_Capacity + Reaches_SE * qnorm(0.025), ymax = Reaches_Capacity + Reaches_SE * qnorm(0.975)), width = 0) + geom_errorbarh(aes(xmin = Points_Capacity + Points_SE * qnorm(0.025), xmax = Points_Capacity + Points_SE * qnorm(0.975)), height = 0) + geom_point(size = 3) + coord_cartesian(xlim = range(my_range), ylim = range(my_range)) + # scale_x_continuous(limits = range(my_range)) + # scale_y_continuous(limits = range(my_range)) + labs(x = "Points", y = "Reaches", title = paste(unique(z$Species), unique(z$Lifestage), "Capacity")) }) comp_p = ggarrange(plotlist = comp_p_list, nrow = 3, ncol = 2, common.legend = T, legend = "bottom") comp_p
The contemporary mean and max of recent escapements (translated to other lifestages), as well as the recovery goal of mean abundance threshold (MAT) and a 25% buffer on top of that (MAT+25%) are shown in Figures \@ref(fig:chnk-parr-fig) - \@ref(fig:sthd-redd-fig), together with estimates of current capacity from master sample points and reaches.
Chinook
bar_list_p$Chinook.Parr
bar_list_p$Chinook.Presmolts
bar_list_p$Chinook.Redds
cap_comp_df %>% filter(Species == "Chinook") %>% rename(Lifestage = lifestage, Watershed = watershed) %>% mutate(Lifestage = recode(Lifestage, "juv_summer" = "Parr", "redds" = "Redds", "juv_winter" = "Presmolts")) %>% mutate(Lifestage = factor(Lifestage, levels = c("Redds", "Parr", "Presmolts"))) %>% mutate(tot_cap_cv = tot_cap_se / tot_cap, Capacity = paste0(prettyNum(round(tot_cap), big.mark = ","), " (", round(tot_cap_cv, 2), ")")) %>% select(Lifestage:source, tot_length, Capacity) %>% rename(Length = tot_length) %>% pivot_wider(names_from = "source", values_from = c("Length", "Capacity"), names_glue = "{source}_{.value}") %>% kable(digits = 1, caption = "Estimated capacities (CVs) for Chinook in various lifestages using the master sample points and linear network reaches for extrapolation. Total watershed length for each method is shown as well.") %>% kable_styling(fixed_thead = T)
Steelhead
bar_list_p$Steelhead.Parr
bar_list_p$Steelhead.Presmolts
bar_list_p$Steelhead.Redds
cap_comp_df %>% filter(Species == "Steelhead") %>% rename(Lifestage = lifestage, Watershed = watershed) %>% mutate(Lifestage = recode(Lifestage, "juv_summer" = "Parr", "redds" = "Redds", "juv_winter" = "Presmolts")) %>% mutate(Lifestage = factor(Lifestage, levels = c("Redds", "Parr", "Presmolts"))) %>% mutate(tot_cap_cv = tot_cap_se / tot_cap, Capacity = paste0(prettyNum(round(tot_cap), big.mark = ","), " (", round(tot_cap_cv, 2), ")")) %>% select(Lifestage:source, tot_length, Capacity) %>% rename(Length = tot_length) %>% pivot_wider(names_from = "source", values_from = c("Length", "Capacity"), names_glue = "{source}_{.value}") %>% kable(digits = 1, caption = "Estimated capacities (CVs) for steelhead in various lifestages using the master sample points and linear network reaches for extrapolation. Total watershed length for each method is shown as well.") %>% kable_styling(fixed_thead = T)
Lemhi Specific
We broke down the Lemhi available capacity into the Upper and Lower Lemhi (defined by the confluence of Hayden Creek).
upper_lem_poly = upper_lem$basin lower_lem_poly = cap_est %>% filter(watershed == "Lemhi") %>% slice(1) %>% pull(inputs) %>% map("wtsd_poly") %>% extract2(1) %>% extract2(1) %>% st_difference(upper_lem_poly) lem_cap = cap_est %>% filter(watershed == "Lemhi") %>% select(lifestage:source) %>% mutate(upp_lem_cap = map(inputs, .f = function(y) { calc_watershed_cap(wtsd_polygon = upper_lem_poly, capacity_sf = y$rch[[1]], spp_range = y$spp_domain[[1]], by_stream = F) %>% rename(n = n_rchs) })) %>% mutate(low_lem_cap = map(inputs, .f = function(y) { calc_watershed_cap(wtsd_polygon = lower_lem_poly, capacity_sf = y$rch[[1]], spp_range = y$spp_domain[[1]], by_stream = F) %>% rename(n = n_rchs) })) %>% pivot_longer(cols = ends_with("lem_cap"), names_to = "area", values_to = "cap") %>% unnest(cols = cap) %>% mutate(area = recode(area, "upp_lem_cap" = "Upper Lemhi", "low_lem_cap" = "Lower Lemhi")) %>% select(-inputs) %>% mutate(across(tot_length, ~. / 1000)) lem_comp_df = lem_cap %>% rename(Watershed = watershed, Lifestage = lifestage) %>% mutate(Lifestage = recode(Lifestage, "juv_summer" = "Parr", "redds" = "Redds", "juv_winter" = "Presmolts")) %>% mutate(Scenario = paste0("QRF_cap_", source)) %>% select(Lifestage:area, Scenario, Abundance = tot_cap) %>% filter(source == "Points") %>% bind_rows(rec_goals %>% filter(Watershed == "Lemhi", Lifestage %in% c("Parr", "Redds", "Presmolts"))) %>% mutate(Lifestage = factor(Lifestage, levels = c("Redds", "Parr", "Presmolts")), Scenario = factor(Scenario, levels = c("QRF_cap_Points", "Mean", "Max", "MAT", "MAT+25%"))) %>% mutate(area = if_else(is.na(area), as.character(Scenario), area)) %>% filter(Scenario %in% c("QRF_cap_Points", "MAT+25%")) %>% mutate(Scenario = fct_expand(Scenario, "Available Capacity (QRF)", "Required Capacity\n(MAT +25%)"), Scenario = fct_recode(Scenario, "Available Capacity (QRF)" = "QRF_cap_Points", "Required Capacity\n(MAT +25%)" = "MAT+25%"), Scenario = fct_rev(Scenario)) %>% mutate(area = factor(area, levels = c("MAT+25%", "Upper Lemhi", "Lower Lemhi"))) %>% arrange(Species, Lifestage, area, Scenario) lem_bar_p = lem_comp_df %>% split(list(.$Species, .$Lifestage)) %>% map(.f = function(df) { df %>% ggplot(aes(y = Scenario, x = Abundance, fill = area)) + geom_col(position = "stack") + scale_fill_manual(values = c("Lower Lemhi" = "darkblue", "Upper Lemhi" = "lightblue", "MAT+25%" = "black"), breaks = c("MAT+25%", "Lower Lemhi", "Upper Lemhi"), name = NULL) + # scale_fill_brewer(palette = "Paired", # name = "Scenario") + theme(legend.position = "bottom", # axis.text.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_blank()) + scale_x_continuous(labels = function(x) format(x, scientific = F)) + labs(x = "Abundance", title = paste(unique(df$Species), unique(df$Lifestage))) })
Chinook
lem_bar_p$Chinook.Parr
lem_bar_p$Chinook.Presmolts
lem_bar_p$Chinook.Redds
Steelhead
lem_bar_p$Steelhead.Parr
lem_bar_p$Steelhead.Presmolts
lem_bar_p$Steelhead.Redds
walk(names(lem_bar_p), function(.x) { ggsave( path = "figures/", filename = paste0(.x, ".png"), plot = lem_bar_p[[.x]], width = 6.5, height = 2 ) })
Discussion
In general, we noted the following patterns in comparing QRF estimates of capacity using master sample points versus a linear network:
- Chinook parr usually had greater capacity estimates with linear network (with a couple exceptions)
- Chinook presmolts had similar capacity estimates (some smaller, some larger)
- Chinook redds had much smaller capacity estimates with linear network
- Steelhead parr had similar or slightly greater capacity estimates with linear network
- Steelhead presmolts had very similar capacity estimates
- Steelhead redds had much smaller capacity estimates with linear network
We do not believe there is anything inherently incorrect in either extrapolation method. Differences in estimates may arise due to differences in available covariates between master sample points and 200m reaches (and therefore differences in the extrapolation models), but also because the master sample point version must "smooth" across a linear network from a series of points (i.e. taking the average fish/m capacity from points on a stream, and multiplying ny the length of that stream), whereas the 200m reaches are already on a linear network, and we are making extrapolation estimates at every single reach. For this reason, we would recommend using the results from the 200m reach extrapolation going forward.
One of the goals of this exercise was to make as close to apple-to-apples comparisons between estimates based on master sample points and estimates based on 200m reaches on a linear network. It needs to be noted that there are a few streams in the upper Salmon area that are not included in the 200m reach linear network (e.g. Big Springs, Little Springs, Pratt Creek and Lee Creek in the Lemhi, Patterson Side Channels in the Pahsimeroi, Yellowbelly Creek in the Upper Salmon, etc.). Therefore, the estimates for the master sample points presented here may not match those in the IRA report exactly. But because of that, the lengths in Tables \@ref(tab:chnk-cap-tab) and \@ref(tab:sthd-cap-tab) are nearly identical, so any differences in capacity estimates here are not simply due to different species ranges being used; we tried to put them on a level playing field in that regard. This does mean that the linear network could be improved, especially in the Lemhi (and less so in the Pahsimeroi) to allow the 200m reach extrapolation to include those missing streams. That will involve generating not just shapefiles of the stream lines, but either calculating directly or finding a proxy for the various covariates that need to be attached to each reach for the extrapolation model.
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.