In mpdougherty/razviz: A package to visualize performance of HEC-RAS models
Purpose
The purpose of this project is to develop a method to streamline the production of longitudinal profile graphs for calibration runs of HEC-RAS hydrologic models. This project will focus on the Upper Mississippi River Flood Risk Management (UMR-FRM) Hydrologic Model. HEC-RAS supports the production of longitudinal profile graphs, but does not provide sufficient flexibility for producing these graphs for long river reaches. Therefore, exporting the model results to another application which permits more automation and more flexible visualization options is required. R has been selected for its flexibility and visualization capabilities.
Longitudinal profile graph required elements
The following graph elements are required to be displayed on the longitudinal profile graphs:
- HEC-RAS modelled water surface elevations by plan -- labels: plan DONE
- Modelled discharges by plan -- labels: plan, discharge, units, extent
- High water marks by event-- labels: none DONE
- Pool elevations -- labels: dam name, pool elevation
- Gages location -- labels: name, stages, river mile, gage zero elevation DONE
- Tributary confluences -- labels: name, river mile DONE
- Leveed area extents -- labels: name, leveed area Not needed, see below
- Levee authorized and surveyed elevations -- symbolize left/right descending bank, smoothing of heights; labels: name DONE
- Bridge elevations -- symbolize bridge underside min/max elevations DONE
- Y-axis Characteristics
- Minimize y-axis range -- Maximize the amount of y-axis resolution for interpreting differences between model results and high water marks.
- Common y-axis range -- Supports direct comparison of y-axis distances between plots. DONE
- Display minor grid lines every foot. DONE
Export HEC-RAS model results
The HEC-RAS Users Manual, Chapter 9 Viewing Results describes how to display tabular data (section 9-13). Tabular data displayed using the Detailed Output Tables, Profile Summary Tables, and User Defined Tables can be transferred to Excel using the "Windows Clipboard" method. Use the following steps to perform the export:
- In HEC-RAS, create the model results summary table.
- Use the menu xxxx|xxxx to copy the summary table to the clipboard.
- Open Excel and paste the clipboard contents into a new sheet.
- Edit the header row to create valid variable names (i.e., no special characters).
- Delete the row of units.
- The Excel file should be saved as
.csv format to facilitate import into R.
Alternatively, data can be exported from HEC-RAS to HEC-DSS format (section 9-34). The DSS-R interface project R package provides functions to read and write data from a HEC-DSS file by making calls to the Java methods exposed in HEC-DSSVue's scripting capability. This method has been tested, but not pursued for this project.
Import HEC-RAS model results into R
Import HEC-RAS model data into R.
library(tidyverse)
# Import csv file of RAS results
# 2014
hydroModel2014 <- read_csv(file = "Data/2014_results.csv")
hydroModel2014 <- add_column(hydroModel2014, Plan = "2014 Event", .after = 4)
hydroModel2014 <- add_column(hydroModel2014, Event = "2014")
# 2013
hydroModel2013 <- read_csv(file = "Data/2013_results.csv")
hydroModel2013 <- add_column(hydroModel2013, Plan = "2013 Event", .after = 4)
hydroModel2013 <- add_column(hydroModel2013, Event = "2013")
# 2012
hydroModel2012 <- read_csv(file = "Data/2012_results.csv")
hydroModel2012 <- add_column(hydroModel2012, Plan = "2012 Event", .after = 4)
hydroModel2012 <- add_column(hydroModel2012, Event = "2012")
#2008
hydroModel2008 <- read_csv(file = "Data/2008_results.csv")
hydroModel2008 <- add_column(hydroModel2008, Plan = "2008 Event", .after = 4)
hydroModel2008 <- add_column(hydroModel2008, Event = "2008")
# Append the events together
hydroModel <- bind_rows(hydroModel2014, hydroModel2013, hydroModel2012, hydroModel2008)
# Assign Plan as factor
hydroModel$Plan <- factor(hydroModel$Plan)
# Remove unneeded data frames
rm(hydroModel2014, hydroModel2013, hydroModel2012, hydroModel2008)
# Import test data
#hydroModel <- read_csv(file = "Data/UMRexample.csv")
# Print data frame
#hydroModel
Manufacture Plan data
Temporary measure to create plan data for testing until actual data is made available.
# 2013 Event
#hydroModel2013 <- add_column(hydroModel, Plan = "2013 Event", .after = 4)
#hydroModel2013 <- add_column(hydroModel2013, Event = "2013")
#hydroModel2013$WS_Elev <- hydroModel2013$WS_Elev + 2.5
# 2008 Event
#hydroModel2008 <- add_column(hydroModel, Plan = "2008 Event", .after = 4)
#hydroModel2008 <- add_column(hydroModel2008, Event = "2008")
#hydroModel2008$WS_Elev <- hydroModel2008$WS_Elev + 5
# 2001 Event
#hydroModel2001 <- add_column(hydroModel, Plan = "2001 Event", .after = 4)
#hydroModel2001 <- add_column(hydroModel2001, Event = "2001")
#hydroModel2001$WS_Elev <- hydroModel2001$WS_Elev - 4
# 1993 Event
#hydroModel1993 <- add_column(hydroModel, Plan = "1993 Event", .after = 4)
#hydroModel1993 <- add_column(hydroModel1993, Event = "1993")
#hydroModel1993$WS_Elev <- hydroModel1993$WS_Elev + 1.5
# Append the events together
#hydroModel <- bind_rows(hydroModel1993, hydroModel2001, hydroModel2008, hydroModel2013)
# Assign Plan as factor
#hydroModel$Plan <- factor(hydroModel$Plan)
# Remove unneeded data frames
#rm(hydroModel1993, hydroModel2001, hydroModel2008, hydroModel2013)
Import levee elevation data
Import levee elevation data into R.
library(dplyr)
levees <- read_csv(file = "Data/Levees_Authorized_Existing.csv")
# Create factors
levees$descending_bank <- factor(levees$descending_bank)
levees$levee <- factor(levees$levee)
# Sort data frame
levees <- arrange(levees, descending_bank, levee, elevation_type, river_mile)
Import high water marks
Import high water mark data into R.
library(lubridate)
high_water <- read_csv(file = "Data/UMR_high_water_marks.csv")
high_water$peak_date <- mdy(high_water$peak_date)
high_water$event <- as.character(year(high_water$peak_date))
high_water <- filter(high_water, event == "2014" | event == "2013" | event == "2012" | event == "2008")
Import bridge elevations
Import a data frame of bridge min/max elevations.
bridges <- read_csv(file = "Data/Bridge_elevations.csv")
bridges
Import gage locations
Import a data frame of gage locations.
gages <- read_csv(file = "Data/gage_locations.csv")
# readr::read_csv mangles newlines in strings. Given R's inscrutable handling of backslashes in strings, replace "\\n" with "\n".
gages$name <- gsub("\\n", "\n", gages$name, fixed=TRUE)
gages
Import river features
Import a data frame of river features.
features <- read_csv(file = "Data/River_features.csv")
# readr::read_csv mangles newlines in strings. Given R's inscrutable handling of backslashes in strings, replace "\\n" with "\n".
features$name <- gsub("\\n", "\n", features$name, fixed=TRUE)
features
Create gage labels
Create a new data frame to store gage label information based on the gage data frame.
# Set the interval for gage stage labels.
stage_interval_label <- 5
# Create the data frame to hold gage stage labels.
gage_labels <- data.frame(name = character(), river_mile = double(), elevation = double(), stage = double())
# Iterate through gages to create gage stage labels.
for (i in gages$name) {
name <- gages[gages$name == i,]$name
river_mile <- gages[gages$name == i,]$river_mile
# Identify the min and max stage label, rounded to the nearest stage_gage_label.
min_stage <- ceiling(gages[gages$name == i,]$min_stage / stage_interval_label) * stage_interval_label
max_stage <- floor(gages[gages$name == i,]$max_stage / stage_interval_label) * stage_interval_label
# Iterate through the min and max stages for each gage creating the labels.
for (j in seq(min_stage, max_stage, by = stage_interval_label)) {
elevation <- gages[gages$name == i,]$elevation + j
stage <- j
gage_labels <- rbind(gage_labels, data.frame(name = name, river_mile = river_mile, elevation = elevation, stage = stage))
}
}
#gage_labels
Create gage boxes
Use the gage data frame to create boxes that can be used to draw gages on plots.
# Set the interval for gage stage boxes.
stage_interval_boxes <- 1
# Set the width of the stage boxes (units: river miles)
box_width <- 0.1
# Create the data frame to hold gage stage boxes.
gage_boxes <- data.frame(gage_stage = character(), river_mile = double(), x = double(), y = double())
# Iterate through gages to create gage stage boxes.
for (i in gages$name) {
name <- gages[gages$name == i,]$name
river_mile <- gages[gages$name == i,]$river_mile
min_stage <- gages[gages$name == i,]$min_stage
max_stage <- gages[gages$name == i,]$max_stage
# Iterate through the min and max stages for each gage creating the boxes.
for (j in seq(min_stage, max_stage, by = stage_interval_boxes)) {
elevation <- gages[gages$name == i,]$elevation + j
stage <- j
gage_stage <- paste(name, stage)
# Calculate the box coordinates for the curent stage. Box coordinate 1 starts in the lower left corner.
gage_box_1 <- data.frame(gage_stage = gage_stage, river_mile = river_mile, x = river_mile + (box_width / 2), y = elevation)
gage_box_2 <- data.frame(gage_stage = gage_stage, river_mile = river_mile, x = river_mile + (box_width / 2), y = elevation + 1)
gage_box_3 <- data.frame(gage_stage = gage_stage, river_mile = river_mile, x = river_mile - (box_width / 2), y = elevation + 1)
gage_box_4 <- data.frame(gage_stage = gage_stage, river_mile = river_mile, x = river_mile - (box_width / 2), y = elevation)
gage_box_5 <- data.frame(gage_stage = gage_stage, river_mile = river_mile, x = river_mile + (box_width / 2), y = elevation)
# Append the box coordinates for the current stage to the main data frame, gage_boxes
gage_boxes <- rbind(gage_boxes, gage_box_1, gage_box_2, gage_box_3, gage_box_4, gage_box_5)
}
}
rm(stage_interval_boxes, box_width, gage_box_1, gage_box_2, gage_box_3, gage_box_4, gage_box_5)
#gage_boxes
Create flat pool elevations
Create a data frame to hold flat pool elevations.
pools <- tribble(
~name, ~start_river_mile, ~end_river_mile, ~elevation,
#----------------------------------------/------------------/----------------/----------
"Dam No. 19 Pool", 410.5, 364.3, 518.2,
"Dan No. 20 Pool", 364.2, 343.2, 480.0,
"Dam No. 21 Pool", 343.3, 324.9, 470.0,
"Dam No. 22 Pool", 325.0, 301.2, 459.5
)
Assign Plot Pages
The purpose of this step is to determine how many plot pages are needed and which river miles will appear on which plot pages.
# Max function that returns zero if vector contains all NAs
max0 <- function (x) ifelse( !all(is.na(x)), max(x, na.rm = TRUE), 0 )
# Min function that returns zero if vector contains all NAs
min0 <- function (x) ifelse( !all(is.na(x)), min(x, na.rm = TRUE), 0 )
# Specify number of river miles per plot
miles_per_plot <- 60
# Determine the number of plots needed
num_plots <- ceiling( (max(hydroModel$River_Sta) - min(hydroModel$River_Sta)) / miles_per_plot )
# Create an empty data frame of plots with fields for starting and ending river miles, min and max y-axis values.
plot_pages <- tibble(plot = 1:num_plots,
start_mile = rep(NA, num_plots),
end_mile = rep(NA, num_plots),
max_y = rep(0, num_plots),
mid_y = rep(0, num_plots),
min_y = rep(0, num_plots),
y_range = rep(0, num_plots),
plot_max_y = rep(0, num_plots),
plot_min_y = rep(0, num_plots),
plot_dif_y = rep(0, num_plots))
# Initialize beginning and ending river miles of the first plot
plot_pages[1,]$start_mile <- max(hydroModel$River_Sta)
plot_pages[1,]$end_mile <- max(hydroModel$River_Sta) - miles_per_plot
# Interate through the remaining plots specifiying their starting and ending river miles
for (j in 2:num_plots) {
plot_pages[j,]$start_mile <- plot_pages[j-1,]$end_mile
plot_pages[j,]$end_mile <- plot_pages[j,]$start_mile - miles_per_plot
}
# Iterage through plot pages determining the min and max y-axis values for each data series
for (k in 1:num_plots) {
# Set start and end mile for the current plot
start_mile <- plot_pages[plot_pages$plot == k,]$start_mile
end_mile <- plot_pages[plot_pages$plot == k,]$end_mile
# Subset relevant data frames for the current plot
hm <- hydroModel[hydroModel$River_Sta <= start_mile & hydroModel$River_Sta >= end_mile,]
hw <- high_water[high_water$river_mile <= start_mile & high_water$river_mile >= end_mile,]
#g <- gages[gages$river_mile <= start_mile & gages$river_mile >= end_mile,]
# Set plot_pages$max_y
aa <- c( max0(hm$WS_Elev), max0(hw$elevation_NAVD88) )
plot_pages[plot_pages$plot == k,]$max_y <- max0(aa[aa > 0])
# Set plot_pages$min_y
bb <- c( min0(hm$WS_Elev), min0(hw$elevation_NAVD88) )
plot_pages[plot_pages$plot == k,]$min_y <- min0(bb[bb > 0])
# Calculate y range
plot_pages[plot_pages$plot == k,]$y_range <- plot_pages[plot_pages$plot == k,]$max_y - plot_pages[plot_pages$plot == k,]$min_y
# Calculate mid y
plot_pages[plot_pages$plot == k,]$mid_y <- plot_pages[plot_pages$plot == k,]$max_y - ( plot_pages[plot_pages$plot == k,]$y_range / 2)
}
# Determine the max y range for all of the plots
plot_y_range <- ceiling(max(plot_pages$y_range))
# Iterate through plots specifying their plot_max_y, plot_min_y, and plot_dif_y
for (l in 1:num_plots) {
plot_pages[l,]$plot_max_y <- ceiling(plot_pages[l,]$mid_y + (plot_y_range / 2))
plot_pages[l,]$plot_min_y <- floor(plot_pages[l,]$mid_y - (plot_y_range / 2))
plot_pages[l,]$plot_dif_y <- plot_pages[l,]$plot_max_y - plot_pages[l,]$plot_min_y
}
# Display plot pages table
plot_pages
Profile Graph Function
library(ggplot2)
library(ggrepel)
library(scales)
library(dplyr)
# longitudinal_profile_plot function
longitudinal_profile_plot <- function(plot_number) {
# Get values from the plot_pages data frame for the current plot
start_mile <- plot_pages[plot_pages$plot == plot_number,]$start_mile
end_mile <- plot_pages[plot_pages$plot == plot_number,]$end_mile
plot_max_y <- plot_pages[plot_pages$plot == plot_number,]$plot_max_y
plot_min_y <- plot_pages[plot_pages$plot == plot_number,]$plot_min_y
# Subset data frames for the current plot
hm <- hydroModel[hydroModel$River_Sta <= start_mile & hydroModel$River_Sta >= end_mile,]
hw <- high_water[high_water$river_mile <= start_mile & high_water$river_mile >= end_mile,]
l <- levees[levees$river_mile <= start_mile & levees$river_mile >= end_mile,]
g <- gages[gages$river_mile <= start_mile & gages$river_mile >= end_mile,]
gl <- gage_labels[gage_labels$river_mile <= start_mile & gage_labels$river_mile >= end_mile,]
gb <- gage_boxes[gage_boxes$river_mile <= start_mile & gage_boxes$river_mile >= end_mile,]
f <- features[features$river_mile <= start_mile & features$river_mile >= end_mile,]
b <- bridges[bridges$river_mile <= start_mile & bridges$river_mile >= end_mile,]
# Create levee labels for the current plot
levee_labels <- summarize(group_by(l, levee, descending_bank),
river_mile = mean(river_mile),
elevation_NAVD88 = mean(elevation_NAVD88))
# Define colors. Inspired by palettes from https://wesandersonpalettes.tumblr.com/ using names from colors().
cols <- c("2008" = "darkslategray4",
"2012" = "cadetblue3",
"2013" = "coral3",
"2014" = "burlywood3",
"LEFT" = "palevioletred2",
"RIGHT" = "palevioletred4")
legend_labels <- c("2008", "2012", "2013", "2014", "Left Bank", "Right Bank")
# Create the plot
p <- ggplot(data = hm,
aes(x = River_Sta, y = WS_Elev, color = Event)) +
geom_line(size = 2) +
scale_color_manual(values = cols, labels = legend_labels) +
theme_bw() +
coord_cartesian(ylim = c(plot_min_y, plot_max_y)) +
scale_x_reverse(minor_breaks = seq(from = ceiling(start_mile), to = floor(end_mile), by = -1),
expand = c(0.01,0)) +
scale_y_continuous(minor_breaks = seq(from = plot_min_y, to = plot_max_y, by = 1),
expand = c(0.01,0),
oob = squish) +
theme(legend.position = c(.99, .99),
legend.justification = c("right", "top"),
legend.background = element_rect(fill = alpha('white', 0.6)),
legend.title = element_blank(),
panel.grid.major = element_line(colour = "grey10", size = 0.1)) +
labs(title = "Upper Mississippi River Hydraulic Model - Keokuk to Thebes",
x = "Miles Above Ohio River",
y = "Elevation (NAVD88 feet)") +
# Draw high water marks
geom_point(inherit.aes = FALSE,
data = hw,
aes(x = river_mile, y = elevation_NAVD88, color = event),
show.legend = FALSE, size = 4) +
# Draw raw levee lines, existing elevation
geom_line(inherit.aes = FALSE,
data = filter(l, elevation_type == "2016_NLD"),
aes(x = river_mile, y = elevation_NAVD88, group = levee, color = descending_bank),
show.legend = FALSE, size = 0.2, alpha = 0.5) +
# Draw smooth levee lines, existing elevation
geom_smooth(inherit.aes = FALSE,
data = filter(l, elevation_type == "2016_NLD"),
aes(x = river_mile, y = elevation_NAVD88, group = levee, color = descending_bank),
method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, size = 0.4) +
# Label levees
geom_label_repel(inherit.aes = FALSE,
data = levee_labels,
aes(x = river_mile, y = elevation_NAVD88, label = levee, color = descending_bank),
show.legend = FALSE, size = 2, force = 1, label.size = 0.1, segment.size = 0,
box.padding = unit(0.1, "lines"), fill = alpha("white", 0.6)) +
# Draw the gage boxes
geom_path(inherit.aes = FALSE,
data = gb,
aes(x = x, y = y, group = gage_stage),
size = 0.1) +
# Label gage stages
geom_text(inherit.aes = FALSE,
data = gl,
aes(x = river_mile, y = elevation, label = stage),
hjust = 0, nudge_x = 0.1, size = 2.5) +
# Label gages
geom_text_repel(inherit.aes = FALSE,
data = g,
aes(x = river_mile, y = (elevation + min_stage) + 2, label = name),
nudge_x = -0.4, angle = 90, size = 3, fontface = "bold", force = 0.1, segment.size = 0) +
# Label river features
geom_text_repel(inherit.aes = FALSE,
data = f,
aes(x = river_mile, y = rep(plot_min_y - 0, length(name)), label = name),
nudge_x = 0, angle = 90, size = 2.5, force = 0.01, segment.size = 0) +
# Draw bridge elevations
geom_errorbar(inherit.aes = FALSE,
data = b,
aes(x = river_mile, ymin = lowest_elevation, ymax = highest_elevation),
width = 0.5, size = 0.2, color = "red4")
return(p)
}
Create Profile Plots
# Iterate through plot pages to produce longitudinal profile plots
for (i in plot_pages$plot) {
# insert vertical white space so that the next figure is on a new page
cat('\\newpage')
# Set the current plot number
plot_number <- plot_pages[plot_pages$plot == i,]$plot
# Create the longitudinal profile plot
print(longitudinal_profile_plot(plot_number))
}
mpdougherty/razviz documentation built on April 1, 2021, 4:16 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.
Purpose
The purpose of this project is to develop a method to streamline the production of longitudinal profile graphs for calibration runs of HEC-RAS hydrologic models. This project will focus on the Upper Mississippi River Flood Risk Management (UMR-FRM) Hydrologic Model. HEC-RAS supports the production of longitudinal profile graphs, but does not provide sufficient flexibility for producing these graphs for long river reaches. Therefore, exporting the model results to another application which permits more automation and more flexible visualization options is required. R has been selected for its flexibility and visualization capabilities.
Longitudinal profile graph required elements
The following graph elements are required to be displayed on the longitudinal profile graphs:
- HEC-RAS modelled water surface elevations by plan -- labels: plan DONE
- Modelled discharges by plan -- labels: plan, discharge, units, extent
- High water marks by event-- labels: none DONE
- Pool elevations -- labels: dam name, pool elevation
- Gages location -- labels: name, stages, river mile, gage zero elevation DONE
- Tributary confluences -- labels: name, river mile DONE
- Leveed area extents -- labels: name, leveed area Not needed, see below
- Levee authorized and surveyed elevations -- symbolize left/right descending bank, smoothing of heights; labels: name DONE
- Bridge elevations -- symbolize bridge underside min/max elevations DONE
- Y-axis Characteristics
- Minimize y-axis range -- Maximize the amount of y-axis resolution for interpreting differences between model results and high water marks.
- Common y-axis range -- Supports direct comparison of y-axis distances between plots. DONE
- Display minor grid lines every foot. DONE
Export HEC-RAS model results
The HEC-RAS Users Manual, Chapter 9 Viewing Results describes how to display tabular data (section 9-13). Tabular data displayed using the Detailed Output Tables, Profile Summary Tables, and User Defined Tables can be transferred to Excel using the "Windows Clipboard" method. Use the following steps to perform the export:
- In HEC-RAS, create the model results summary table.
- Use the menu xxxx|xxxx to copy the summary table to the clipboard.
- Open Excel and paste the clipboard contents into a new sheet.
- Edit the header row to create valid variable names (i.e., no special characters).
- Delete the row of units.
- The Excel file should be saved as
.csvformat to facilitate import into R.
Alternatively, data can be exported from HEC-RAS to HEC-DSS format (section 9-34). The DSS-R interface project R package provides functions to read and write data from a HEC-DSS file by making calls to the Java methods exposed in HEC-DSSVue's scripting capability. This method has been tested, but not pursued for this project.
Import HEC-RAS model results into R
Import HEC-RAS model data into R.
library(tidyverse) # Import csv file of RAS results # 2014 hydroModel2014 <- read_csv(file = "Data/2014_results.csv") hydroModel2014 <- add_column(hydroModel2014, Plan = "2014 Event", .after = 4) hydroModel2014 <- add_column(hydroModel2014, Event = "2014") # 2013 hydroModel2013 <- read_csv(file = "Data/2013_results.csv") hydroModel2013 <- add_column(hydroModel2013, Plan = "2013 Event", .after = 4) hydroModel2013 <- add_column(hydroModel2013, Event = "2013") # 2012 hydroModel2012 <- read_csv(file = "Data/2012_results.csv") hydroModel2012 <- add_column(hydroModel2012, Plan = "2012 Event", .after = 4) hydroModel2012 <- add_column(hydroModel2012, Event = "2012") #2008 hydroModel2008 <- read_csv(file = "Data/2008_results.csv") hydroModel2008 <- add_column(hydroModel2008, Plan = "2008 Event", .after = 4) hydroModel2008 <- add_column(hydroModel2008, Event = "2008") # Append the events together hydroModel <- bind_rows(hydroModel2014, hydroModel2013, hydroModel2012, hydroModel2008) # Assign Plan as factor hydroModel$Plan <- factor(hydroModel$Plan) # Remove unneeded data frames rm(hydroModel2014, hydroModel2013, hydroModel2012, hydroModel2008) # Import test data #hydroModel <- read_csv(file = "Data/UMRexample.csv") # Print data frame #hydroModel
Manufacture Plan data
Temporary measure to create plan data for testing until actual data is made available.
# 2013 Event #hydroModel2013 <- add_column(hydroModel, Plan = "2013 Event", .after = 4) #hydroModel2013 <- add_column(hydroModel2013, Event = "2013") #hydroModel2013$WS_Elev <- hydroModel2013$WS_Elev + 2.5 # 2008 Event #hydroModel2008 <- add_column(hydroModel, Plan = "2008 Event", .after = 4) #hydroModel2008 <- add_column(hydroModel2008, Event = "2008") #hydroModel2008$WS_Elev <- hydroModel2008$WS_Elev + 5 # 2001 Event #hydroModel2001 <- add_column(hydroModel, Plan = "2001 Event", .after = 4) #hydroModel2001 <- add_column(hydroModel2001, Event = "2001") #hydroModel2001$WS_Elev <- hydroModel2001$WS_Elev - 4 # 1993 Event #hydroModel1993 <- add_column(hydroModel, Plan = "1993 Event", .after = 4) #hydroModel1993 <- add_column(hydroModel1993, Event = "1993") #hydroModel1993$WS_Elev <- hydroModel1993$WS_Elev + 1.5 # Append the events together #hydroModel <- bind_rows(hydroModel1993, hydroModel2001, hydroModel2008, hydroModel2013) # Assign Plan as factor #hydroModel$Plan <- factor(hydroModel$Plan) # Remove unneeded data frames #rm(hydroModel1993, hydroModel2001, hydroModel2008, hydroModel2013)
Import levee elevation data
Import levee elevation data into R.
library(dplyr) levees <- read_csv(file = "Data/Levees_Authorized_Existing.csv") # Create factors levees$descending_bank <- factor(levees$descending_bank) levees$levee <- factor(levees$levee) # Sort data frame levees <- arrange(levees, descending_bank, levee, elevation_type, river_mile)
Import high water marks
Import high water mark data into R.
library(lubridate) high_water <- read_csv(file = "Data/UMR_high_water_marks.csv") high_water$peak_date <- mdy(high_water$peak_date) high_water$event <- as.character(year(high_water$peak_date)) high_water <- filter(high_water, event == "2014" | event == "2013" | event == "2012" | event == "2008")
Import bridge elevations
Import a data frame of bridge min/max elevations.
bridges <- read_csv(file = "Data/Bridge_elevations.csv") bridges
Import gage locations
Import a data frame of gage locations.
gages <- read_csv(file = "Data/gage_locations.csv") # readr::read_csv mangles newlines in strings. Given R's inscrutable handling of backslashes in strings, replace "\\n" with "\n". gages$name <- gsub("\\n", "\n", gages$name, fixed=TRUE) gages
Import river features
Import a data frame of river features.
features <- read_csv(file = "Data/River_features.csv") # readr::read_csv mangles newlines in strings. Given R's inscrutable handling of backslashes in strings, replace "\\n" with "\n". features$name <- gsub("\\n", "\n", features$name, fixed=TRUE) features
Create gage labels
Create a new data frame to store gage label information based on the gage data frame.
# Set the interval for gage stage labels. stage_interval_label <- 5 # Create the data frame to hold gage stage labels. gage_labels <- data.frame(name = character(), river_mile = double(), elevation = double(), stage = double()) # Iterate through gages to create gage stage labels. for (i in gages$name) { name <- gages[gages$name == i,]$name river_mile <- gages[gages$name == i,]$river_mile # Identify the min and max stage label, rounded to the nearest stage_gage_label. min_stage <- ceiling(gages[gages$name == i,]$min_stage / stage_interval_label) * stage_interval_label max_stage <- floor(gages[gages$name == i,]$max_stage / stage_interval_label) * stage_interval_label # Iterate through the min and max stages for each gage creating the labels. for (j in seq(min_stage, max_stage, by = stage_interval_label)) { elevation <- gages[gages$name == i,]$elevation + j stage <- j gage_labels <- rbind(gage_labels, data.frame(name = name, river_mile = river_mile, elevation = elevation, stage = stage)) } } #gage_labels
Create gage boxes
Use the gage data frame to create boxes that can be used to draw gages on plots.
# Set the interval for gage stage boxes. stage_interval_boxes <- 1 # Set the width of the stage boxes (units: river miles) box_width <- 0.1 # Create the data frame to hold gage stage boxes. gage_boxes <- data.frame(gage_stage = character(), river_mile = double(), x = double(), y = double()) # Iterate through gages to create gage stage boxes. for (i in gages$name) { name <- gages[gages$name == i,]$name river_mile <- gages[gages$name == i,]$river_mile min_stage <- gages[gages$name == i,]$min_stage max_stage <- gages[gages$name == i,]$max_stage # Iterate through the min and max stages for each gage creating the boxes. for (j in seq(min_stage, max_stage, by = stage_interval_boxes)) { elevation <- gages[gages$name == i,]$elevation + j stage <- j gage_stage <- paste(name, stage) # Calculate the box coordinates for the curent stage. Box coordinate 1 starts in the lower left corner. gage_box_1 <- data.frame(gage_stage = gage_stage, river_mile = river_mile, x = river_mile + (box_width / 2), y = elevation) gage_box_2 <- data.frame(gage_stage = gage_stage, river_mile = river_mile, x = river_mile + (box_width / 2), y = elevation + 1) gage_box_3 <- data.frame(gage_stage = gage_stage, river_mile = river_mile, x = river_mile - (box_width / 2), y = elevation + 1) gage_box_4 <- data.frame(gage_stage = gage_stage, river_mile = river_mile, x = river_mile - (box_width / 2), y = elevation) gage_box_5 <- data.frame(gage_stage = gage_stage, river_mile = river_mile, x = river_mile + (box_width / 2), y = elevation) # Append the box coordinates for the current stage to the main data frame, gage_boxes gage_boxes <- rbind(gage_boxes, gage_box_1, gage_box_2, gage_box_3, gage_box_4, gage_box_5) } } rm(stage_interval_boxes, box_width, gage_box_1, gage_box_2, gage_box_3, gage_box_4, gage_box_5) #gage_boxes
Create flat pool elevations
Create a data frame to hold flat pool elevations.
pools <- tribble( ~name, ~start_river_mile, ~end_river_mile, ~elevation, #----------------------------------------/------------------/----------------/---------- "Dam No. 19 Pool", 410.5, 364.3, 518.2, "Dan No. 20 Pool", 364.2, 343.2, 480.0, "Dam No. 21 Pool", 343.3, 324.9, 470.0, "Dam No. 22 Pool", 325.0, 301.2, 459.5 )
Assign Plot Pages
The purpose of this step is to determine how many plot pages are needed and which river miles will appear on which plot pages.
# Max function that returns zero if vector contains all NAs max0 <- function (x) ifelse( !all(is.na(x)), max(x, na.rm = TRUE), 0 ) # Min function that returns zero if vector contains all NAs min0 <- function (x) ifelse( !all(is.na(x)), min(x, na.rm = TRUE), 0 ) # Specify number of river miles per plot miles_per_plot <- 60 # Determine the number of plots needed num_plots <- ceiling( (max(hydroModel$River_Sta) - min(hydroModel$River_Sta)) / miles_per_plot ) # Create an empty data frame of plots with fields for starting and ending river miles, min and max y-axis values. plot_pages <- tibble(plot = 1:num_plots, start_mile = rep(NA, num_plots), end_mile = rep(NA, num_plots), max_y = rep(0, num_plots), mid_y = rep(0, num_plots), min_y = rep(0, num_plots), y_range = rep(0, num_plots), plot_max_y = rep(0, num_plots), plot_min_y = rep(0, num_plots), plot_dif_y = rep(0, num_plots)) # Initialize beginning and ending river miles of the first plot plot_pages[1,]$start_mile <- max(hydroModel$River_Sta) plot_pages[1,]$end_mile <- max(hydroModel$River_Sta) - miles_per_plot # Interate through the remaining plots specifiying their starting and ending river miles for (j in 2:num_plots) { plot_pages[j,]$start_mile <- plot_pages[j-1,]$end_mile plot_pages[j,]$end_mile <- plot_pages[j,]$start_mile - miles_per_plot } # Iterage through plot pages determining the min and max y-axis values for each data series for (k in 1:num_plots) { # Set start and end mile for the current plot start_mile <- plot_pages[plot_pages$plot == k,]$start_mile end_mile <- plot_pages[plot_pages$plot == k,]$end_mile # Subset relevant data frames for the current plot hm <- hydroModel[hydroModel$River_Sta <= start_mile & hydroModel$River_Sta >= end_mile,] hw <- high_water[high_water$river_mile <= start_mile & high_water$river_mile >= end_mile,] #g <- gages[gages$river_mile <= start_mile & gages$river_mile >= end_mile,] # Set plot_pages$max_y aa <- c( max0(hm$WS_Elev), max0(hw$elevation_NAVD88) ) plot_pages[plot_pages$plot == k,]$max_y <- max0(aa[aa > 0]) # Set plot_pages$min_y bb <- c( min0(hm$WS_Elev), min0(hw$elevation_NAVD88) ) plot_pages[plot_pages$plot == k,]$min_y <- min0(bb[bb > 0]) # Calculate y range plot_pages[plot_pages$plot == k,]$y_range <- plot_pages[plot_pages$plot == k,]$max_y - plot_pages[plot_pages$plot == k,]$min_y # Calculate mid y plot_pages[plot_pages$plot == k,]$mid_y <- plot_pages[plot_pages$plot == k,]$max_y - ( plot_pages[plot_pages$plot == k,]$y_range / 2) } # Determine the max y range for all of the plots plot_y_range <- ceiling(max(plot_pages$y_range)) # Iterate through plots specifying their plot_max_y, plot_min_y, and plot_dif_y for (l in 1:num_plots) { plot_pages[l,]$plot_max_y <- ceiling(plot_pages[l,]$mid_y + (plot_y_range / 2)) plot_pages[l,]$plot_min_y <- floor(plot_pages[l,]$mid_y - (plot_y_range / 2)) plot_pages[l,]$plot_dif_y <- plot_pages[l,]$plot_max_y - plot_pages[l,]$plot_min_y } # Display plot pages table plot_pages
Profile Graph Function
library(ggplot2) library(ggrepel) library(scales) library(dplyr) # longitudinal_profile_plot function longitudinal_profile_plot <- function(plot_number) { # Get values from the plot_pages data frame for the current plot start_mile <- plot_pages[plot_pages$plot == plot_number,]$start_mile end_mile <- plot_pages[plot_pages$plot == plot_number,]$end_mile plot_max_y <- plot_pages[plot_pages$plot == plot_number,]$plot_max_y plot_min_y <- plot_pages[plot_pages$plot == plot_number,]$plot_min_y # Subset data frames for the current plot hm <- hydroModel[hydroModel$River_Sta <= start_mile & hydroModel$River_Sta >= end_mile,] hw <- high_water[high_water$river_mile <= start_mile & high_water$river_mile >= end_mile,] l <- levees[levees$river_mile <= start_mile & levees$river_mile >= end_mile,] g <- gages[gages$river_mile <= start_mile & gages$river_mile >= end_mile,] gl <- gage_labels[gage_labels$river_mile <= start_mile & gage_labels$river_mile >= end_mile,] gb <- gage_boxes[gage_boxes$river_mile <= start_mile & gage_boxes$river_mile >= end_mile,] f <- features[features$river_mile <= start_mile & features$river_mile >= end_mile,] b <- bridges[bridges$river_mile <= start_mile & bridges$river_mile >= end_mile,] # Create levee labels for the current plot levee_labels <- summarize(group_by(l, levee, descending_bank), river_mile = mean(river_mile), elevation_NAVD88 = mean(elevation_NAVD88)) # Define colors. Inspired by palettes from https://wesandersonpalettes.tumblr.com/ using names from colors(). cols <- c("2008" = "darkslategray4", "2012" = "cadetblue3", "2013" = "coral3", "2014" = "burlywood3", "LEFT" = "palevioletred2", "RIGHT" = "palevioletred4") legend_labels <- c("2008", "2012", "2013", "2014", "Left Bank", "Right Bank") # Create the plot p <- ggplot(data = hm, aes(x = River_Sta, y = WS_Elev, color = Event)) + geom_line(size = 2) + scale_color_manual(values = cols, labels = legend_labels) + theme_bw() + coord_cartesian(ylim = c(plot_min_y, plot_max_y)) + scale_x_reverse(minor_breaks = seq(from = ceiling(start_mile), to = floor(end_mile), by = -1), expand = c(0.01,0)) + scale_y_continuous(minor_breaks = seq(from = plot_min_y, to = plot_max_y, by = 1), expand = c(0.01,0), oob = squish) + theme(legend.position = c(.99, .99), legend.justification = c("right", "top"), legend.background = element_rect(fill = alpha('white', 0.6)), legend.title = element_blank(), panel.grid.major = element_line(colour = "grey10", size = 0.1)) + labs(title = "Upper Mississippi River Hydraulic Model - Keokuk to Thebes", x = "Miles Above Ohio River", y = "Elevation (NAVD88 feet)") + # Draw high water marks geom_point(inherit.aes = FALSE, data = hw, aes(x = river_mile, y = elevation_NAVD88, color = event), show.legend = FALSE, size = 4) + # Draw raw levee lines, existing elevation geom_line(inherit.aes = FALSE, data = filter(l, elevation_type == "2016_NLD"), aes(x = river_mile, y = elevation_NAVD88, group = levee, color = descending_bank), show.legend = FALSE, size = 0.2, alpha = 0.5) + # Draw smooth levee lines, existing elevation geom_smooth(inherit.aes = FALSE, data = filter(l, elevation_type == "2016_NLD"), aes(x = river_mile, y = elevation_NAVD88, group = levee, color = descending_bank), method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, size = 0.4) + # Label levees geom_label_repel(inherit.aes = FALSE, data = levee_labels, aes(x = river_mile, y = elevation_NAVD88, label = levee, color = descending_bank), show.legend = FALSE, size = 2, force = 1, label.size = 0.1, segment.size = 0, box.padding = unit(0.1, "lines"), fill = alpha("white", 0.6)) + # Draw the gage boxes geom_path(inherit.aes = FALSE, data = gb, aes(x = x, y = y, group = gage_stage), size = 0.1) + # Label gage stages geom_text(inherit.aes = FALSE, data = gl, aes(x = river_mile, y = elevation, label = stage), hjust = 0, nudge_x = 0.1, size = 2.5) + # Label gages geom_text_repel(inherit.aes = FALSE, data = g, aes(x = river_mile, y = (elevation + min_stage) + 2, label = name), nudge_x = -0.4, angle = 90, size = 3, fontface = "bold", force = 0.1, segment.size = 0) + # Label river features geom_text_repel(inherit.aes = FALSE, data = f, aes(x = river_mile, y = rep(plot_min_y - 0, length(name)), label = name), nudge_x = 0, angle = 90, size = 2.5, force = 0.01, segment.size = 0) + # Draw bridge elevations geom_errorbar(inherit.aes = FALSE, data = b, aes(x = river_mile, ymin = lowest_elevation, ymax = highest_elevation), width = 0.5, size = 0.2, color = "red4") return(p) }
Create Profile Plots
# Iterate through plot pages to produce longitudinal profile plots for (i in plot_pages$plot) { # insert vertical white space so that the next figure is on a new page cat('\\newpage') # Set the current plot number plot_number <- plot_pages[plot_pages$plot == i,]$plot # Create the longitudinal profile plot print(longitudinal_profile_plot(plot_number)) }
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.