library(elfgen)
library(stringr)
basepath='/var/www/R';
source("/var/www/R/config.R")
options(scipen = 999)
watershed_hydrocode = '02070006'
source(paste0(github_location,"/om/R/summarize/rseg_elfgen.R"))
ichthy.localpath = "c:/usr/local/home/git/r-dh-ecohydro/IcthyMapsDATA/Frimpong_etal_2015_IcthymapsVA_region.csv"
# Retrieve dataset of interest
# You may enter either a 6, 8, 10, or 12-digit HUC code
# *Notes:
# By default the ichthy dataset is downloaded to a temp directory, however this may be overridden by
# supplying a local path of interest using the input parameter "ichthy.localpath"
# 6-digit HUCs like the following example may take a few minutes to process with elfdata() due to the
# large number of contained IchthyMaps historical stream fish distribution data
watershed.df <- elfdata(watershed.code = watershed_hydrocode, ichthy.localpath = ichthy.localpath)
flow_col = 'Q08'
watershed_formatted = watershed.df[,c(flow_col, 'NT.TOTAL.UNIQUE', 'watershed.code')]
breakpt <- bkpt_pwit(
"watershed.df" = watershed_formatted,
"quantile" = 0.92, "blo" = 100, "bhi" = 800
)
breakpt = 530
# other options:
# ymax - can be good, but if there are low x, hi y outliers it is bad
# breakpt <- bkpt_ymax("watershed.df" = watershed.df)
# breakpt = 530 # or manually defined based on state averages
elf <- elfgen(
"watershed.df" = watershed_formatted,
"quantile" = 0.9,
"breakpt" = breakpt,
"yaxis_thresh" = 53,
"xlabel" = "Mean Aug. Flow (ft3/s)",
"ylabel" = "Fish Species Richness",
plot_title = paste("Species Richness vs. Flow", watershed_feature$name)
)
elf
flow_col = "Q08"
confidence <- elfgen_confidence(
elf,
rseg.name=watershed_feature$name,
outlet_flow,
yaxis_thresh = 53,
cuf,
flow_col
)
#Elf$plot post - posted underneath elfgen_richness_change_huc_level scenario property------------
dR10 <- richness_change(elf$stats, "pctchg" = 10)
dR20 <- richness_change(elf$stats, "pctchg" = 20)
dR30 <- richness_change(elf$stats, "pctchg" = 30)
dR40 <- richness_change(elf$stats, "pctchg" = 40)
message(paste("Richness Change at -10%", dR10))
message(paste("Richness Change at -20%", dR20))
message(paste("Richness Change at -30%", dR30))
message(paste("Richness Change at -40%", dR40))
flow_col = 'Q09'
watershed_formatted = watershed.df[,c(flow_col, 'NT.TOTAL.UNIQUE', 'watershed.code')]
elf <- elfgen(
"watershed.df" = watershed_formatted,
"quantile" = 0.9,
"breakpt" = breakpt,
"yaxis_thresh" = 53,
"xlabel" = "Mean September Flow (ft3/s)",
"ylabel" = "Fish Species Richness"
)
elf
dR10 <- richness_change(elf$stats, "pctchg" = 10)
dR20 <- richness_change(elf$stats, "pctchg" = 20)
dR30 <- richness_change(elf$stats, "pctchg" = 30)
dR40 <- richness_change(elf$stats, "pctchg" = 40)
message(paste("Richness Change at -10%", dR10))
message(paste("Richness Change at -20%", dR20))
message(paste("Richness Change at -30%", dR30))
message(paste("Richness Change at -40%", dR40))
flow_col = 'MAF'
watershed_formatted = watershed.df[,c(flow_col, 'NT.TOTAL.UNIQUE', 'watershed.code')]
elf <- elfgen(
"watershed.df" = watershed_formatted,
"quantile" = 0.9,
"breakpt" = breakpt,
"yaxis_thresh" = 53,
"xlabel" = "Mean August Flow (ft3/s)",
"ylabel" = "Fish Species Richness"
)
elf
#Elf$plot post - posted underneath elfgen_richness_change_huc_level scenario property------------
dR10 <- richness_change(elf$stats, "pctchg" = 10)
dR20 <- richness_change(elf$stats, "pctchg" = 20)
dR30 <- richness_change(elf$stats, "pctchg" = 30)
dR40 <- richness_change(elf$stats, "pctchg" = 40)
message(paste("Richness Change at -10%", dR10))
message(paste("Richness Change at -20%", dR20))
message(paste("Richness Change at -30%", dR30))
message(paste("Richness Change at -40%", dR40))
# OR - Use dh/vahydro EDAS query, which selects Q and NT column already
watershed_feature <- RomFeature$new(
ds, list(bundle='watershed', ftype='nhd_huc8', hydrocode=watershed_hydrocode),
TRUE
)
erom_varkey = 'erom_q0001e_sept'
watershed_edas_elfdata <- dh_elfdata(
watershed_feature, ws_varkey = erom_varkey,
bio_varkey = 'aqbio_nt_total', ds
)
nhd_col = elfgen_varkey_nhd_col(erom_varkey)
watershed_formatted = watershed_edas_elfdata[,c(nhd_col, 'y_metric', 'hydrocode')]
breakpt <- bkpt_pwit(
"watershed.df" = watershed_formatted,
"quantile" = 0.9, "blo" = 50, "bhi" = 500
)
elf <- elfgen(
"watershed.df" = watershed_formatted,
"quantile" = 0.9,
"breakpt" = breakpt,
"yaxis_thresh" = 53,
"xlabel" = "Mean Sep Flow (ft3/s)",
"ylabel" = "Fish Species Richness"
)
elf
cuf=0.05
outlet_flow = 51.0 # just a guess need to query somehow for baseline mo flow
confidence <- elfgen_confidence(
elf,watershed_feature$name,outlet_flow,yaxis_thresh=53,cuf, nhd_col
)
confidence
#Elf$plot post - posted underneath elfgen_richness_change_huc_level scenario property------------
dR10 <- richness_change(elf$stats, "pctchg" = 10)
dR20 <- richness_change(elf$stats, "pctchg" = 20)
dR30 <- richness_change(elf$stats, "pctchg" = 30)
dR40 <- richness_change(elf$stats, "pctchg" = 40)
message(paste("Richness Change at -10%", dR10))
message(paste("Richness Change at -20%", dR20))
message(paste("Richness Change at -30%", dR30))
message(paste("Richness Change at -40%", dR40))
nhd_cols = elfgen_varkey_nhd_col(FALSE)
rsq_vals = c()
mo_vars = c()
da_col = nhd_cols[["erom_q0001e_mean"]]$elfgen
breakpt <- bkpt_pwit(
"watershed.df" = watershed.df[,c(da_col, 'y_metric', 'hydrocode')],
"quantile" = quant, "blo" = 50, "bhi" = 5000
)
breakpt = 500
for (i in names(nhd_cols) ) {
ncol = as.character(nhd_cols[[i]]$elfgen)
watershed_formatted = watershed.df[,c(ncol, 'y_metric', 'hydrocode')]
watershed_formatted$da <- watershed.df[,da_col]
mo_elf <- elfgen(
"watershed.df" = watershed_formatted,
"quantile" = quant,
"breakpt" = breakpt,
"yaxis_thresh" = 53,
"xlabel" = ncol,
"ylabel" = "Fish Species Richness",
break_var = "da"
)
#print(mo_elf)
nhd_cols[[i]]$rsq = mo_elf$stats$rsquared_adj
rsq_vals = cbind(rsq_vals, mo_elf$stats$rsquared_adj)
mo_vars = cbind(mo_vars, str_split(i,"_")[[1]][3])
}
barplot(
rsq_vals,
names = mo_vars,
ylim=c(0,1.0),
main=watershed_feature$name,
xlab = "NT = f(Q) R^2",
ylab = "Month",
cex.lab=1,
cex.axis=1,
las=2
)
