Nothing
inst/doc/publication-reproduction.R
In rivr: Steady and Unsteady Open-Channel Flow Computation
## ----include = FALSE----------------------------------------------------------
library(knitr)
opts_chunk$set(echo = TRUE, eval = FALSE)
## ----the-code-----------------------------------------------------------------
# ## ---- load-libs -------------------------------------------------------
# library(knitr)
# opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = "fig/",
# warning = FALSE, message = FALSE, error = FALSE, results = "asis",
# out.width = '\\textwidth', cache = TRUE, cache.path = "cache/")
# library(rivr)
# library(ggplot2)
# library(scales)
# library(RColorBrewer)
# library(dplyr)
# library(xtable)
#
# ## ---- startparams -----------------------------------------------------
# plotopts = list(theme_bw(), xlab(expression(Distance~from~control~section~~(ft))))
# g = 32.2
# Cm = 1.486
# slope = 0.001
# mannings = 0.045
# flow = 250
# width = 100
# sideslope = 0
# # calculate control depth as 1ft above the normal depth
# depth.m1 = round(1 + normal_depth(slope, mannings, flow, 2, Cm,
# width, sideslope), 2)
# depth.m2 = round(1.1*critical_depth(flow, 1, g, width, sideslope), 2)
# rivdist = 3000
#
# ## ---- restest ---------------------------------------------------------
# # Test sensitivity of step size
# stepsizes = c(500, 100, 50, 10, 1)
# resolution.test.m1 = list()
# resolution.test.m2 = list()
# for(r in stepsizes){
# resolution.test.m1[[paste0('dx=', r)]] = compute_profile(slope, mannings,
# flow, depth.m1, Cm, g, width, sideslope, stepdist = r,
# totaldist = rivdist)
# resolution.test.m2[[paste0('dx=', r)]] = compute_profile(slope, mannings,
# flow, depth.m2, Cm, g, width, sideslope, stepdist = r,
# totaldist = rivdist)
# }
# resolution.plot = NULL
# for(lbl in names(resolution.test.m1)){
# f1 = resolution.test.m1[[lbl]]
# f2 = resolution.test.m2[[lbl]]
# f1['type'] = 'M1'
# f2['type'] = 'M2'
# f = rbind(f1, f2)
# f['run'] = lbl
# f['res'] = as.numeric(substr(lbl, 4, 10L))
# resolution.plot = rbind(resolution.plot, f)
# }
# resolution.plot['res'] = factor(resolution.plot$res, levels = rev(sort(stepsizes)))
# resolution.plot['run'] = factor(resolution.plot$run,
# levels = names(resolution.test.m1))
# resolution.plot['type'] = factor(resolution.plot$type, levels = c('M1', 'M2'))
# stepsize = 10
#
# ## ----resfig ----------------------------------------------------------
# ggplot(resolution.plot, aes(x = x, y = y + z)) + plotopts +
# geom_line(aes(color = res, linetype = res), size = 1) +
# facet_wrap(~type) + ylab(expression(River~stage~~(ft))) +
# scale_linetype_manual("resolution (ft)", labels = levels(resolution.plot$res),
# values = c("solid", "longdash", "dashed", "dotdash", "dotted")) +
# scale_color_manual("resolution (ft)", values = brewer.pal(7, "YlGnBu")[3:7],
# labels = levels(resolution.plot$res))
#
# ## ----roughness -------------------------------------------------------
# roughness.test.m1 = list()
# roughness.test.m2 = list()
# roughnesses = 0.045*seq(0.5, 1.5, length = 5)
# for(n in roughnesses){
# roughness.test.m1[[paste0('n=',n)]] = compute_profile(slope, n, flow,
# depth.m1, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
# roughness.test.m2[[paste0('n=',n)]] = compute_profile(slope, n, flow,
# depth.m2, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
# }
# roughness.plot = NULL
# for(lbl in names(roughness.test.m1)){
# f1 = roughness.test.m1[[lbl]]
# f2 = roughness.test.m2[[lbl]]
# f1['type'] = 'M1'
# f2['type'] = 'M2'
# f = rbind(f1, f2)
# f['run'] = lbl
# f['mannings'] = as.numeric(substr(lbl, 4, 10L))
# roughness.plot = rbind(roughness.plot, f)
# }
# roughness.plot['run'] = factor(roughness.plot$run,
# levels = names(roughness.test.m1))
# roughness.plot['mannings'] = factor(roughness.plot$mannings,
# levels = sort(roughnesses))
# roughness.plot['type'] = factor(roughness.plot$type, levels=c('M1', 'M2'))
#
# ## ---- roughplot -------------------------------------------------------
# ggplot(roughness.plot, aes(x = x, y = y + z)) +
# geom_line(aes(linetype = mannings, color = mannings), size = 1) + plotopts +
# facet_wrap(~ type) + ylab(expression(River~stage~~(ft))) +
# scale_linetype_manual("Bed roughness", labels = levels(roughness.plot$mannings),
# values = c("solid", "longdash", "dashed", "dotdash", "dotted")) +
# scale_color_manual("Bed roughness", labels = levels(roughness.plot$mannings),
# values = brewer.pal(6, "Oranges")[2:6])
#
# ## ---- relroughplot ----------------------------------------------------
# roughness.rel.m1 = list()
# roughness.rel.m2 = list()
# for(n in roughnesses){
# thisyn = normal_depth(slope, n, flow, 2, Cm, width, sideslope)
# roughness.rel.m1[[paste0('n=',n)]] = compute_profile(slope, n, flow,
# 1.25*thisyn, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
# roughness.rel.m2[[paste0('n=',n)]] = compute_profile(slope, n, flow,
# 0.75*thisyn, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
# roughness.rel.m1[[paste0('n=',n)]]['pd.yn'] =
# (roughness.rel.m1[[paste0('n=',n)]]$y - thisyn)/thisyn
# roughness.rel.m2[[paste0('n=',n)]]['pd.yn'] =
# (roughness.rel.m2[[paste0('n=',n)]]$y - thisyn)/thisyn
# }
# roughness.rel.plot = NULL
# for(lbl in names(roughness.rel.m1)){
# f1 = roughness.rel.m1[[lbl]]
# f2 = roughness.rel.m2[[lbl]]
# f1['type'] = 'M1'
# f2['type'] = 'M2'
# f = rbind(f1, f2)
# f['run'] = lbl
# f['mannings'] = as.numeric(substr(lbl, 4, 10L))
# roughness.rel.plot = rbind(roughness.rel.plot, f)
# }
# roughness.rel.plot['mannings'] = factor(roughness.plot$mannings,
# levels = sort(roughnesses))
# roughness.rel.plot['run'] = factor(roughness.rel.plot$run,
# levels = names(roughness.rel.m1))
# roughness.rel.plot['type'] = factor(roughness.rel.plot$type,
# levels = c('M1', 'M2'))
# ggplot(roughness.rel.plot, aes(x = x, y = pd.yn)) +
# geom_line(aes(linetype = mannings, color = mannings), size = 1) +
# plotopts + facet_wrap(~type, scales = 'free_y') +
# scale_y_continuous(expression(Percent~difference~from~normal~depth), labels = percent) +
# scale_linetype_manual("Bed roughness", labels = levels(roughness.rel.plot$mannings),
# values = c("solid", "longdash", "dashed", "dotdash", "dotted")) +
# scale_color_manual("Bed roughness", labels = levels(roughness.rel.plot$mannings),
# values = brewer.pal(6, "Oranges")[2:6])
#
#
# ## ---- loadup-kwm ------------------------------------------------------
# oldscipen = options('scipen')
# options(scipen = 1000)
#
# slope = 0.001
# extent = 150000
# mannings = 0.045
# B = 100
# SS = 0
# Cm = 1.486
# g = 32.2
# iflow = 250
#
# # keep Courant number at 0.7 to balance temporal and spatial resolution
# idepth = normal_depth(slope, mannings, 250, 10, 1.485, 100, SS)
# iarea = channel_geom(idepth, B, SS)["A"]
# cn = 0.7*iarea/iflow
# # define upstream boundary condition assuming a timestep in seconds
# bcfunc = function(x)
# ifelse(x < 9000, 250 + (750/pi)*(1 - cos(pi*x/(60*75))), 250)
# bctime = 76000
# xnodes = extent/c(25, 50, 125, 250, 500, 1000, 2000, 5000) + 1
# myxres = extent/(xnodes - 1)
# mytres = cn*myxres
# plotopts = list(theme_bw())
#
# # run the model
# modtime = list()
# modelresults = list()
# for(i in seq(length(xnodes))){
# numnodes = xnodes[i]
# xstep = myxres[i]
# tstep = mytres[i]
# bc = bcfunc(seq(0, bctime, by=tstep))
# mp = c(1, as.integer(50000/xstep + 1), as.integer(100000/xstep + 1), numnodes)
# mt = as.integer(round(seq(1, length(bc), length.out=10)))[seq(7)]
# if(xstep == min(myxres)){
# mt = as.integer(round(seq(1, length(bc), length.out=125)))
# mtslice = mt[round(seq(1, length(mt), length.out=7))]
# }
# lbl = paste0('dx=',xstep)
# modtime[[lbl]] <- system.time(modelresults[[lbl]] <-
# route_wave(slope, mannings, Cm, g, B, SS, iflow, bc,
# timestep=tstep, spacestep=xstep, numnodes=numnodes,
# monitor.nodes=mp, monitor.times=mt, engine="Kinematic"))
# modelresults[[lbl]]['deltax'] = xstep
# modelresults[[lbl]]['deltat'] = tstep
# modelresults[[lbl]]['computationtime'] = modtime[[lbl]][[3]]
# modelresults[[lbl]]['label'] = lbl
# }
# allresults = do.call(rbind.data.frame, modelresults)
# row.names(allresults) = NULL
# names(allresults)[which(names(allresults)=='flow')] = "Q"
# names(allresults)[which(names(allresults)=='distance')] = "x"
# names(allresults)[which(names(allresults)=='time')] = "t"
#
# ## ---- floodhydrograph -------------------------------------------------
# floodhydro = data.frame(time = seq(15000), flow=bcfunc(seq(15000)))
# ggplot(floodhydro, aes(x=time/60, y=flow)) + geom_line(color="darkblue",
# size = 1) + plotopts +
# xlab(expression(time~~(minutes))) + ylab(expression(flow~~(ft^3~~s^-1)))
#
# ## ---- routeplot-kinematic ---------------------------------------------
# # still plot
# surfaces = group_by(filter(allresults, monitor.type=='timestep',
# deltax==min(deltax)))
# ggplot(filter(surfaces, step %in% mtslice[2:5]),
# aes(x=x, y=Q, linetype=factor(round(t/60)), color=factor(round(t/60)))) +
# geom_line(size = 1) + plotopts +
# scale_color_manual("time\n(minutes)", values = brewer.pal(7, "YlGnBu")[3:6]) +
# scale_linetype_manual("time\n(minutes)", values = c("solid",
# "dashed", "dotdash", "dotted")) +
# xlab(expression(distance~~downstream~~(ft))) +
# ylab(expression(flow~~(ft^3~~s^-1)))
# # animation
# routeylims = c(min(surfaces$Q), max(surfaces$Q))
# routexlims = c(min(surfaces$x), max(surfaces$x))
# thisdt = unique(surfaces$deltat)
# for(n in unique(surfaces$step))
# print(
# ggplot(filter(surfaces, step==n), aes(x=x, y=Q)) +
# geom_line(color = "darkblue", size = 1) + plotopts +
# scale_linetype_discrete("time\n(minutes)") +
# scale_x_continuous(expression(distance~~downstream~~(ft)),
# limits=routexlims) +
# scale_y_continuous(expression(flow~~(ft^3~~s^-1)), limits=routeylims) +
# ggtitle(paste(format(round((n - 1)*thisdt/60), width = 4), "minutes"))
# )
#
# ## ---- xsectionplot-kinematic ------------------------------------------
# xsections = filter(allresults, monitor.type == "node",
# x %in% c(0,50000))[c("x","t","Q", "deltax")]
# xsections['deltax'] = factor(paste(xsections$deltax, "feet"),
# levels=paste(sort(unique(xsections$deltax)), "feet"))
# ggplot(xsections, aes(x=t/60, y=Q, linetype=factor(format(x, big.mark=",",
# trim=TRUE), levels = format(sort(unique(x)), big.mark=",", trim=TRUE)),
# color = deltax)) +
# geom_path(size = 1) + facet_grid(deltax~.) + plotopts +
# scale_y_continuous(expression(flow~~(ft^3~~s^-1)), breaks=c(300, 500,700)) +
# scale_linetype_manual("downstream\ndistance (ft)", values = c("solid",
# "dashed", "dotted")) +
# scale_x_continuous(expression(time~~(min)), limits=c(0,600)) +
# scale_color_manual(values = c("#7570B3", "#8DA0CB", "#66C2A5", "#A6D854",
# "#E5C494", "#E6AB02", "#FC8D62", "#E78AC3"), guide = FALSE)
#
# ## ---- peakstable-kinematic --------------------------------------------
# peaks = as.data.frame(summarize(group_by(filter(allresults, x==50000,
# monitor.type=='node'), computationtime, deltax, deltat, monitor.type),
# peak.flow=max(Q), time.to.peak=t[which(Q==peak.flow)]))
# real.peak = max(floodhydro$flow)
# peaks['flow.percent.error'] = (peaks$peak.flow - real.peak)/real.peak
# tbld = as.data.frame(peaks[c("deltax", "deltat", "flow.percent.error",
# "computationtime")])
# tbld["deltax"] = round(tbld$deltax)
# tbld["deltat"] = round(tbld$deltat, 2)
# tbld["flow.percent.error"] = round(100*tbld$flow.percent.error, 2)
# names(tbld) = c("$\\Delta x$ (ft)", "$\\Delta t$ (s)", "\\% error (peak flow)",
# "cost (s)")
# ptable = xtable(tbld[order(tbld[,1]),], digits = c(0,0,2,2,2))
# print(ptable, include.rownames = FALSE, sanitize.colnames.function = identity,
# sanitize.text.function = function(x){x}, floating=FALSE, hline.after=NULL,
# add.to.row=list(pos=list(-1,0, nrow(ptable)),
# command=c('\\toprule ', '\\midrule ', '\\bottomrule ')),
# format.args = list(big.mark = ","))
#
# ## ---- loadup-dwm ------------------------------------------------------
# oldscipen = options('scipen')
# options(scipen = 1000)
#
# slope = 0.001
# extent = 150000
# mannings = 0.045
# B = 100
# SS = 0
# g = 32.2
# Cm = 1.486
# iflow = 250
#
# # define upstream boundary condition assuming a timestep in seconds
# idepth = normal_depth(slope, mannings, 250, 10, Cm, 100, SS)
# iarea = channel_geom(idepth, B, SS)[["A"]]
# # keep Courant number at 0.06 to balance temporal and spatial resolution
# cn = 0.06*iarea/250
# bcfunc = function(x)
# ifelse(x < 9000, 250 + (750/pi)*(1 - cos(pi*x/(60*75))), 250)
# bctime = 76000
# xnodes = extent/c(50, 125, 250, 500, 1000) + 1
# myxres = extent/(xnodes - 1)
# mytres = cn*myxres
# plotopts = list(theme_bw())
#
# # run the model
# modtime = list()
# modelresults = list()
# for(i in seq(length(xnodes))){
# numnodes = xnodes[i]
# xstep = myxres[i]
# tstep = mytres[i]
# bc = bcfunc(seq(0, bctime, by=tstep))
# dc = rep(-1, length(bc))
# mp = c(1, as.integer(c(1000, 50000, 100000, 149000)/xstep + 1), numnodes)
# mt = as.integer(round(seq(1, length(bc), length.out=10)))[seq(7)]
# if(xstep == min(myxres)){
# mt = as.integer(round(seq(1, length(bc), length.out=125)))
# mtslice = mt[round(seq(1, length(mt), length.out=7))]
# }
# lbl = paste0('dx=',xstep)
# modtime[[lbl]] <- system.time(modelresults[[lbl]] <-
# route_wave(slope, mannings, Cm, g, B, SS, iflow, bc, dc, timestep=tstep,
# spacestep=xstep, numnodes=numnodes, monitor.nodes=mp, monitor.times=mt,
# engine="Dynamic", boundary.type="QQ"))
# modelresults[[lbl]]['deltax'] = xstep
# modelresults[[lbl]]['deltat'] = tstep
# modelresults[[lbl]]['computationtime'] = modtime[[lbl]][[3]]
# modelresults[[lbl]]['label'] = lbl
# }
# # combine results
# allresults = do.call(rbind.data.frame, modelresults)
# row.names(allresults) = NULL
# names(allresults)[which(names(allresults)=='flow')] = "Q"
# names(allresults)[which(names(allresults)=='distance')] = "x"
# names(allresults)[which(names(allresults)=='time')] = "t"
# allresults['t'] = (allresults$step - 1)*allresults$deltat
# allresults['x'] = (allresults$node - 1)*allresults$deltax
#
# ## ---- routeplot-characteristic ----------------------------------------
# # still plot
# surfaces = filter(allresults, monitor.type=='timestep',
# deltax==min(deltax))
# ggplot(filter(surfaces, step %in% mtslice[2:5]),
# aes(x=x, y=Q, linetype=factor(round(t/60)), color=factor(round(t/60)))) +
# geom_line(size = 1) + plotopts +
# scale_color_manual("time\n(minutes)", values = brewer.pal(7, "YlGnBu")[3:6]) +
# scale_linetype_manual("time\n(minutes)", values = c("solid",
# "dashed", "dotdash", "dotted")) +
# xlab(expression(distance~~downstream~~(ft))) +
# ylab(expression(flow~~(ft^3~~s^-1)))
# # animation
# routeylims = c(min(surfaces$Q), max(surfaces$Q))
# routexlims = c(min(surfaces$x), max(surfaces$x))
# thisdt = unique(surfaces$deltat)
# for(n in unique(surfaces$step))
# print(
# ggplot(filter(surfaces, step==n), aes(x=x, y=Q)) +
# geom_line(color = "darkblue", size = 1) + plotopts +
# scale_linetype_discrete("time\n(minutes)") +
# scale_x_continuous(expression(distance~~downstream~~(ft)),
# limits=routexlims) +
# scale_y_continuous(expression(flow~~(ft^3~~s^-1)), limits = routeylims) +
# ggtitle(paste(format(round((n - 1)*thisdt/60), width = 4), "minutes"))
# )
#
# ## ---- xsectionplot-characteristic -------------------------------------
# xsections = filter(allresults, monitor.type == "node",
# x %in% c(0, 50000, 150000))[c("x","t","Q", "deltax")]
# xsections['deltax'] = factor(paste(xsections$deltax, "feet"),
# levels=paste(sort(unique(xsections$deltax)), "feet"))
# data(waterolympics)
# realdat = NULL
# for(x in unique(xsections$deltax)){
# thisdat = waterolympics[waterolympics$t > 200*60,]
# thisdat['deltax'] = x
# realdat = rbind(realdat, thisdat)
# }
# realdat['deltax'] = factor(realdat$deltax)
# ggplot(xsections, aes(x=t/60, y=Q, linetype=factor(format(x, big.mark=",",
# trim=TRUE), levels = format(sort(unique(x)), big.mark=",", trim=TRUE)),
# color = deltax)) +
# geom_path(size = 1) + geom_point(data = realdat[realdat$t > 200*60,], color = "black") +
# facet_grid(deltax~.) +scale_y_continuous(expression(flow~~(ft^3~~s^-1)),
# breaks=c(300, 500,700)) +
# scale_linetype_manual("downstream\ndistance (ft)", values = c("solid",
# "dashed", "dotted")) + plotopts +
# scale_x_continuous(expression(time~~(min)), limits=c(0,1250)) +
# scale_color_manual(values = c("#8DA0CB", "#66C2A5", "#A6D854", "#E5C494",
# "#E6AB02"), guide = FALSE)
#
# ## ---- peakstable-characteristic ---------------------------------------
# peaks = as.data.frame(summarize(group_by(filter(allresults, x==50000,
# monitor.type=='node'), computationtime, deltax, deltat, monitor.type),
# peak.flow=max(Q), time.to.peak=t[which(Q==peak.flow)]))
# real.peak.flow = max(waterolympics$Q)
# real.time.to.peak = waterolympics$t[which(waterolympics$Q==real.peak.flow)][1]
# peaks['flow.percent.error'] =
# (peaks$peak.flow - real.peak.flow)/real.peak.flow
# peaks['time.percent.error'] =
# (peaks$time.to.peak - real.time.to.peak)/real.time.to.peak
# tbld = as.data.frame(peaks[c("deltax", "deltat", "flow.percent.error",
# "time.percent.error", "computationtime")])
# tbld["flow.percent.error"] = round(100*tbld$flow.percent.error, 2)
# tbld["time.percent.error"] = round(100*tbld$time.percent.error, 2)
# names(tbld) = c("$\\Delta x$ (ft)", "$\\Delta t$ (s)", "\\% error (peak flow)",
# "\\% error (time to peak)", "cost (s)")
# ptable = xtable(tbld[order(tbld[,1]),], digits = c(0,0,2,2,2,2))
# print(ptable, include.rownames=FALSE, sanitize.colnames.function=identity,
# sanitize.text.function = function(x){x}, floating=FALSE, hline.after=NULL,
# add.to.row=list(pos=list(-1,0, nrow(ptable)),
# command=c('\\toprule ', '\\midrule ', '\\bottomrule ')),
# format.args = list(big.mark = ","))
#
# ## ---- loadup-boundaries -----------------------------------------------
# oldscipen = options('scipen')
# options(scipen = 1000)
#
# plotopts = list(theme_bw())
# slope = 0.00008
# extent = 5000
# mannings = 0.013
# B = 6.1
# SS = 1.5
# g = 9.81
# Cm = 1
#
# ic = 126
# id = 5.79
# ia = channel_geom(id, B, SS)[["A"]]
# CN = 0.9
#
# dx = 10
# dt = round(dx*CN/(ic/ia + sqrt(id*g)), 2)
# numnodes = extent/dx + 1
#
# bctime = 2000
# bc = rep(id, round(bctime/dt) + 1)
# dc = rep(0, length(bc))
# dt = round(bctime/(length(bc) - 1), 2)
#
# CN = dt*(ic/ia + sqrt(id*g))/dx
#
# mp = c(1, as.integer(c(1500, 2500, 3000, 5000)/dx + 1))
# mt = as.integer(round(seq(0, length(bc)-1, by=25))) + 1L
# mtslice = c(1, as.integer(c(500, 1000, 1500, 2000)/dt + 1))
#
# dclose.lax = route_wave(slope, mannings, Cm, g, B, SS, ic, bc, dc,
# timestep=dt, spacestep=dx, numnodes=numnodes, monitor.nodes=mp,
# monitor.times=mt, engine="Dynamic", scheme="Lax", boundary.type="yQ")
# dclose.mac = route_wave(slope, mannings, Cm, g, B, SS, ic, bc, dc,
# timestep=dt, spacestep=dx, numnodes=numnodes, monitor.nodes=mp,
# monitor.times=mt, engine="Dynamic", scheme="MacCormack", boundary.type="yQ")
# dclose.lax['scheme'] = "Lax diffusive"
# dclose.mac['scheme'] = "MacCormack predictor-corrector"
# dclose = rbind(dclose.lax, dclose.mac)
# row.names(dclose) = NULL
# dclose["CN"] = dt*(dclose$flow/dclose$area + sqrt(dclose$depth*g))/dx
#
# ## ---- through-time ----------------------------------------------------
# dclose.times = filter(dclose, monitor.type=="timestep")
# ggplot(filter(dclose.times, step %in% mtslice),
# aes(x=distance, y=depth, linetype = factor(round(time)),
# color = factor(round(time)))) + geom_path(size = 1) +
# scale_y_continuous(expression(depth~~(ft))) + facet_wrap(~scheme) +
# scale_x_continuous(expression(distance~~downstream~~(ft))) + plotopts +
# scale_color_manual("time (s)", values = brewer.pal(7, "YlGnBu")[3:7]) +
# scale_linetype_manual("time (s)", values = c("solid", "longdash", "dashed",
# "dotdash", "dotted"))
# # animation
# routeylims = c(min(dclose.times$depth), max(dclose.times$depth))
# routexlims = c(min(dclose.times$distance), max(dclose.times$distance))
# for(n in unique(dclose.times$step))
# print(
# ggplot(filter(dclose.times, step==n), aes(x=distance, y=depth)) +
# geom_line(color = "darkblue", size = 1) + plotopts +
# scale_linetype_discrete("time\n(minutes)") +
# scale_x_continuous(expression(distance~~downstream~~(ft)),
# limits=routexlims) +
# scale_y_continuous(expression(depth~~(ft)), limits=routeylims) +
# ggtitle(paste(format(round((n - 1)*dt), width = 4), "seconds")) +
# facet_wrap(~scheme)
# )
Try the rivr package in your browser
Any scripts or data that you put into this service are public.
rivr documentation built on Jan. 21, 2021, 5:06 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.
## ----include = FALSE----------------------------------------------------------
library(knitr)
opts_chunk$set(echo = TRUE, eval = FALSE)
## ----the-code-----------------------------------------------------------------
# ## ---- load-libs -------------------------------------------------------
# library(knitr)
# opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = "fig/",
# warning = FALSE, message = FALSE, error = FALSE, results = "asis",
# out.width = '\\textwidth', cache = TRUE, cache.path = "cache/")
# library(rivr)
# library(ggplot2)
# library(scales)
# library(RColorBrewer)
# library(dplyr)
# library(xtable)
#
# ## ---- startparams -----------------------------------------------------
# plotopts = list(theme_bw(), xlab(expression(Distance~from~control~section~~(ft))))
# g = 32.2
# Cm = 1.486
# slope = 0.001
# mannings = 0.045
# flow = 250
# width = 100
# sideslope = 0
# # calculate control depth as 1ft above the normal depth
# depth.m1 = round(1 + normal_depth(slope, mannings, flow, 2, Cm,
# width, sideslope), 2)
# depth.m2 = round(1.1*critical_depth(flow, 1, g, width, sideslope), 2)
# rivdist = 3000
#
# ## ---- restest ---------------------------------------------------------
# # Test sensitivity of step size
# stepsizes = c(500, 100, 50, 10, 1)
# resolution.test.m1 = list()
# resolution.test.m2 = list()
# for(r in stepsizes){
# resolution.test.m1[[paste0('dx=', r)]] = compute_profile(slope, mannings,
# flow, depth.m1, Cm, g, width, sideslope, stepdist = r,
# totaldist = rivdist)
# resolution.test.m2[[paste0('dx=', r)]] = compute_profile(slope, mannings,
# flow, depth.m2, Cm, g, width, sideslope, stepdist = r,
# totaldist = rivdist)
# }
# resolution.plot = NULL
# for(lbl in names(resolution.test.m1)){
# f1 = resolution.test.m1[[lbl]]
# f2 = resolution.test.m2[[lbl]]
# f1['type'] = 'M1'
# f2['type'] = 'M2'
# f = rbind(f1, f2)
# f['run'] = lbl
# f['res'] = as.numeric(substr(lbl, 4, 10L))
# resolution.plot = rbind(resolution.plot, f)
# }
# resolution.plot['res'] = factor(resolution.plot$res, levels = rev(sort(stepsizes)))
# resolution.plot['run'] = factor(resolution.plot$run,
# levels = names(resolution.test.m1))
# resolution.plot['type'] = factor(resolution.plot$type, levels = c('M1', 'M2'))
# stepsize = 10
#
# ## ----resfig ----------------------------------------------------------
# ggplot(resolution.plot, aes(x = x, y = y + z)) + plotopts +
# geom_line(aes(color = res, linetype = res), size = 1) +
# facet_wrap(~type) + ylab(expression(River~stage~~(ft))) +
# scale_linetype_manual("resolution (ft)", labels = levels(resolution.plot$res),
# values = c("solid", "longdash", "dashed", "dotdash", "dotted")) +
# scale_color_manual("resolution (ft)", values = brewer.pal(7, "YlGnBu")[3:7],
# labels = levels(resolution.plot$res))
#
# ## ----roughness -------------------------------------------------------
# roughness.test.m1 = list()
# roughness.test.m2 = list()
# roughnesses = 0.045*seq(0.5, 1.5, length = 5)
# for(n in roughnesses){
# roughness.test.m1[[paste0('n=',n)]] = compute_profile(slope, n, flow,
# depth.m1, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
# roughness.test.m2[[paste0('n=',n)]] = compute_profile(slope, n, flow,
# depth.m2, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
# }
# roughness.plot = NULL
# for(lbl in names(roughness.test.m1)){
# f1 = roughness.test.m1[[lbl]]
# f2 = roughness.test.m2[[lbl]]
# f1['type'] = 'M1'
# f2['type'] = 'M2'
# f = rbind(f1, f2)
# f['run'] = lbl
# f['mannings'] = as.numeric(substr(lbl, 4, 10L))
# roughness.plot = rbind(roughness.plot, f)
# }
# roughness.plot['run'] = factor(roughness.plot$run,
# levels = names(roughness.test.m1))
# roughness.plot['mannings'] = factor(roughness.plot$mannings,
# levels = sort(roughnesses))
# roughness.plot['type'] = factor(roughness.plot$type, levels=c('M1', 'M2'))
#
# ## ---- roughplot -------------------------------------------------------
# ggplot(roughness.plot, aes(x = x, y = y + z)) +
# geom_line(aes(linetype = mannings, color = mannings), size = 1) + plotopts +
# facet_wrap(~ type) + ylab(expression(River~stage~~(ft))) +
# scale_linetype_manual("Bed roughness", labels = levels(roughness.plot$mannings),
# values = c("solid", "longdash", "dashed", "dotdash", "dotted")) +
# scale_color_manual("Bed roughness", labels = levels(roughness.plot$mannings),
# values = brewer.pal(6, "Oranges")[2:6])
#
# ## ---- relroughplot ----------------------------------------------------
# roughness.rel.m1 = list()
# roughness.rel.m2 = list()
# for(n in roughnesses){
# thisyn = normal_depth(slope, n, flow, 2, Cm, width, sideslope)
# roughness.rel.m1[[paste0('n=',n)]] = compute_profile(slope, n, flow,
# 1.25*thisyn, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
# roughness.rel.m2[[paste0('n=',n)]] = compute_profile(slope, n, flow,
# 0.75*thisyn, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
# roughness.rel.m1[[paste0('n=',n)]]['pd.yn'] =
# (roughness.rel.m1[[paste0('n=',n)]]$y - thisyn)/thisyn
# roughness.rel.m2[[paste0('n=',n)]]['pd.yn'] =
# (roughness.rel.m2[[paste0('n=',n)]]$y - thisyn)/thisyn
# }
# roughness.rel.plot = NULL
# for(lbl in names(roughness.rel.m1)){
# f1 = roughness.rel.m1[[lbl]]
# f2 = roughness.rel.m2[[lbl]]
# f1['type'] = 'M1'
# f2['type'] = 'M2'
# f = rbind(f1, f2)
# f['run'] = lbl
# f['mannings'] = as.numeric(substr(lbl, 4, 10L))
# roughness.rel.plot = rbind(roughness.rel.plot, f)
# }
# roughness.rel.plot['mannings'] = factor(roughness.plot$mannings,
# levels = sort(roughnesses))
# roughness.rel.plot['run'] = factor(roughness.rel.plot$run,
# levels = names(roughness.rel.m1))
# roughness.rel.plot['type'] = factor(roughness.rel.plot$type,
# levels = c('M1', 'M2'))
# ggplot(roughness.rel.plot, aes(x = x, y = pd.yn)) +
# geom_line(aes(linetype = mannings, color = mannings), size = 1) +
# plotopts + facet_wrap(~type, scales = 'free_y') +
# scale_y_continuous(expression(Percent~difference~from~normal~depth), labels = percent) +
# scale_linetype_manual("Bed roughness", labels = levels(roughness.rel.plot$mannings),
# values = c("solid", "longdash", "dashed", "dotdash", "dotted")) +
# scale_color_manual("Bed roughness", labels = levels(roughness.rel.plot$mannings),
# values = brewer.pal(6, "Oranges")[2:6])
#
#
# ## ---- loadup-kwm ------------------------------------------------------
# oldscipen = options('scipen')
# options(scipen = 1000)
#
# slope = 0.001
# extent = 150000
# mannings = 0.045
# B = 100
# SS = 0
# Cm = 1.486
# g = 32.2
# iflow = 250
#
# # keep Courant number at 0.7 to balance temporal and spatial resolution
# idepth = normal_depth(slope, mannings, 250, 10, 1.485, 100, SS)
# iarea = channel_geom(idepth, B, SS)["A"]
# cn = 0.7*iarea/iflow
# # define upstream boundary condition assuming a timestep in seconds
# bcfunc = function(x)
# ifelse(x < 9000, 250 + (750/pi)*(1 - cos(pi*x/(60*75))), 250)
# bctime = 76000
# xnodes = extent/c(25, 50, 125, 250, 500, 1000, 2000, 5000) + 1
# myxres = extent/(xnodes - 1)
# mytres = cn*myxres
# plotopts = list(theme_bw())
#
# # run the model
# modtime = list()
# modelresults = list()
# for(i in seq(length(xnodes))){
# numnodes = xnodes[i]
# xstep = myxres[i]
# tstep = mytres[i]
# bc = bcfunc(seq(0, bctime, by=tstep))
# mp = c(1, as.integer(50000/xstep + 1), as.integer(100000/xstep + 1), numnodes)
# mt = as.integer(round(seq(1, length(bc), length.out=10)))[seq(7)]
# if(xstep == min(myxres)){
# mt = as.integer(round(seq(1, length(bc), length.out=125)))
# mtslice = mt[round(seq(1, length(mt), length.out=7))]
# }
# lbl = paste0('dx=',xstep)
# modtime[[lbl]] <- system.time(modelresults[[lbl]] <-
# route_wave(slope, mannings, Cm, g, B, SS, iflow, bc,
# timestep=tstep, spacestep=xstep, numnodes=numnodes,
# monitor.nodes=mp, monitor.times=mt, engine="Kinematic"))
# modelresults[[lbl]]['deltax'] = xstep
# modelresults[[lbl]]['deltat'] = tstep
# modelresults[[lbl]]['computationtime'] = modtime[[lbl]][[3]]
# modelresults[[lbl]]['label'] = lbl
# }
# allresults = do.call(rbind.data.frame, modelresults)
# row.names(allresults) = NULL
# names(allresults)[which(names(allresults)=='flow')] = "Q"
# names(allresults)[which(names(allresults)=='distance')] = "x"
# names(allresults)[which(names(allresults)=='time')] = "t"
#
# ## ---- floodhydrograph -------------------------------------------------
# floodhydro = data.frame(time = seq(15000), flow=bcfunc(seq(15000)))
# ggplot(floodhydro, aes(x=time/60, y=flow)) + geom_line(color="darkblue",
# size = 1) + plotopts +
# xlab(expression(time~~(minutes))) + ylab(expression(flow~~(ft^3~~s^-1)))
#
# ## ---- routeplot-kinematic ---------------------------------------------
# # still plot
# surfaces = group_by(filter(allresults, monitor.type=='timestep',
# deltax==min(deltax)))
# ggplot(filter(surfaces, step %in% mtslice[2:5]),
# aes(x=x, y=Q, linetype=factor(round(t/60)), color=factor(round(t/60)))) +
# geom_line(size = 1) + plotopts +
# scale_color_manual("time\n(minutes)", values = brewer.pal(7, "YlGnBu")[3:6]) +
# scale_linetype_manual("time\n(minutes)", values = c("solid",
# "dashed", "dotdash", "dotted")) +
# xlab(expression(distance~~downstream~~(ft))) +
# ylab(expression(flow~~(ft^3~~s^-1)))
# # animation
# routeylims = c(min(surfaces$Q), max(surfaces$Q))
# routexlims = c(min(surfaces$x), max(surfaces$x))
# thisdt = unique(surfaces$deltat)
# for(n in unique(surfaces$step))
# print(
# ggplot(filter(surfaces, step==n), aes(x=x, y=Q)) +
# geom_line(color = "darkblue", size = 1) + plotopts +
# scale_linetype_discrete("time\n(minutes)") +
# scale_x_continuous(expression(distance~~downstream~~(ft)),
# limits=routexlims) +
# scale_y_continuous(expression(flow~~(ft^3~~s^-1)), limits=routeylims) +
# ggtitle(paste(format(round((n - 1)*thisdt/60), width = 4), "minutes"))
# )
#
# ## ---- xsectionplot-kinematic ------------------------------------------
# xsections = filter(allresults, monitor.type == "node",
# x %in% c(0,50000))[c("x","t","Q", "deltax")]
# xsections['deltax'] = factor(paste(xsections$deltax, "feet"),
# levels=paste(sort(unique(xsections$deltax)), "feet"))
# ggplot(xsections, aes(x=t/60, y=Q, linetype=factor(format(x, big.mark=",",
# trim=TRUE), levels = format(sort(unique(x)), big.mark=",", trim=TRUE)),
# color = deltax)) +
# geom_path(size = 1) + facet_grid(deltax~.) + plotopts +
# scale_y_continuous(expression(flow~~(ft^3~~s^-1)), breaks=c(300, 500,700)) +
# scale_linetype_manual("downstream\ndistance (ft)", values = c("solid",
# "dashed", "dotted")) +
# scale_x_continuous(expression(time~~(min)), limits=c(0,600)) +
# scale_color_manual(values = c("#7570B3", "#8DA0CB", "#66C2A5", "#A6D854",
# "#E5C494", "#E6AB02", "#FC8D62", "#E78AC3"), guide = FALSE)
#
# ## ---- peakstable-kinematic --------------------------------------------
# peaks = as.data.frame(summarize(group_by(filter(allresults, x==50000,
# monitor.type=='node'), computationtime, deltax, deltat, monitor.type),
# peak.flow=max(Q), time.to.peak=t[which(Q==peak.flow)]))
# real.peak = max(floodhydro$flow)
# peaks['flow.percent.error'] = (peaks$peak.flow - real.peak)/real.peak
# tbld = as.data.frame(peaks[c("deltax", "deltat", "flow.percent.error",
# "computationtime")])
# tbld["deltax"] = round(tbld$deltax)
# tbld["deltat"] = round(tbld$deltat, 2)
# tbld["flow.percent.error"] = round(100*tbld$flow.percent.error, 2)
# names(tbld) = c("$\\Delta x$ (ft)", "$\\Delta t$ (s)", "\\% error (peak flow)",
# "cost (s)")
# ptable = xtable(tbld[order(tbld[,1]),], digits = c(0,0,2,2,2))
# print(ptable, include.rownames = FALSE, sanitize.colnames.function = identity,
# sanitize.text.function = function(x){x}, floating=FALSE, hline.after=NULL,
# add.to.row=list(pos=list(-1,0, nrow(ptable)),
# command=c('\\toprule ', '\\midrule ', '\\bottomrule ')),
# format.args = list(big.mark = ","))
#
# ## ---- loadup-dwm ------------------------------------------------------
# oldscipen = options('scipen')
# options(scipen = 1000)
#
# slope = 0.001
# extent = 150000
# mannings = 0.045
# B = 100
# SS = 0
# g = 32.2
# Cm = 1.486
# iflow = 250
#
# # define upstream boundary condition assuming a timestep in seconds
# idepth = normal_depth(slope, mannings, 250, 10, Cm, 100, SS)
# iarea = channel_geom(idepth, B, SS)[["A"]]
# # keep Courant number at 0.06 to balance temporal and spatial resolution
# cn = 0.06*iarea/250
# bcfunc = function(x)
# ifelse(x < 9000, 250 + (750/pi)*(1 - cos(pi*x/(60*75))), 250)
# bctime = 76000
# xnodes = extent/c(50, 125, 250, 500, 1000) + 1
# myxres = extent/(xnodes - 1)
# mytres = cn*myxres
# plotopts = list(theme_bw())
#
# # run the model
# modtime = list()
# modelresults = list()
# for(i in seq(length(xnodes))){
# numnodes = xnodes[i]
# xstep = myxres[i]
# tstep = mytres[i]
# bc = bcfunc(seq(0, bctime, by=tstep))
# dc = rep(-1, length(bc))
# mp = c(1, as.integer(c(1000, 50000, 100000, 149000)/xstep + 1), numnodes)
# mt = as.integer(round(seq(1, length(bc), length.out=10)))[seq(7)]
# if(xstep == min(myxres)){
# mt = as.integer(round(seq(1, length(bc), length.out=125)))
# mtslice = mt[round(seq(1, length(mt), length.out=7))]
# }
# lbl = paste0('dx=',xstep)
# modtime[[lbl]] <- system.time(modelresults[[lbl]] <-
# route_wave(slope, mannings, Cm, g, B, SS, iflow, bc, dc, timestep=tstep,
# spacestep=xstep, numnodes=numnodes, monitor.nodes=mp, monitor.times=mt,
# engine="Dynamic", boundary.type="QQ"))
# modelresults[[lbl]]['deltax'] = xstep
# modelresults[[lbl]]['deltat'] = tstep
# modelresults[[lbl]]['computationtime'] = modtime[[lbl]][[3]]
# modelresults[[lbl]]['label'] = lbl
# }
# # combine results
# allresults = do.call(rbind.data.frame, modelresults)
# row.names(allresults) = NULL
# names(allresults)[which(names(allresults)=='flow')] = "Q"
# names(allresults)[which(names(allresults)=='distance')] = "x"
# names(allresults)[which(names(allresults)=='time')] = "t"
# allresults['t'] = (allresults$step - 1)*allresults$deltat
# allresults['x'] = (allresults$node - 1)*allresults$deltax
#
# ## ---- routeplot-characteristic ----------------------------------------
# # still plot
# surfaces = filter(allresults, monitor.type=='timestep',
# deltax==min(deltax))
# ggplot(filter(surfaces, step %in% mtslice[2:5]),
# aes(x=x, y=Q, linetype=factor(round(t/60)), color=factor(round(t/60)))) +
# geom_line(size = 1) + plotopts +
# scale_color_manual("time\n(minutes)", values = brewer.pal(7, "YlGnBu")[3:6]) +
# scale_linetype_manual("time\n(minutes)", values = c("solid",
# "dashed", "dotdash", "dotted")) +
# xlab(expression(distance~~downstream~~(ft))) +
# ylab(expression(flow~~(ft^3~~s^-1)))
# # animation
# routeylims = c(min(surfaces$Q), max(surfaces$Q))
# routexlims = c(min(surfaces$x), max(surfaces$x))
# thisdt = unique(surfaces$deltat)
# for(n in unique(surfaces$step))
# print(
# ggplot(filter(surfaces, step==n), aes(x=x, y=Q)) +
# geom_line(color = "darkblue", size = 1) + plotopts +
# scale_linetype_discrete("time\n(minutes)") +
# scale_x_continuous(expression(distance~~downstream~~(ft)),
# limits=routexlims) +
# scale_y_continuous(expression(flow~~(ft^3~~s^-1)), limits = routeylims) +
# ggtitle(paste(format(round((n - 1)*thisdt/60), width = 4), "minutes"))
# )
#
# ## ---- xsectionplot-characteristic -------------------------------------
# xsections = filter(allresults, monitor.type == "node",
# x %in% c(0, 50000, 150000))[c("x","t","Q", "deltax")]
# xsections['deltax'] = factor(paste(xsections$deltax, "feet"),
# levels=paste(sort(unique(xsections$deltax)), "feet"))
# data(waterolympics)
# realdat = NULL
# for(x in unique(xsections$deltax)){
# thisdat = waterolympics[waterolympics$t > 200*60,]
# thisdat['deltax'] = x
# realdat = rbind(realdat, thisdat)
# }
# realdat['deltax'] = factor(realdat$deltax)
# ggplot(xsections, aes(x=t/60, y=Q, linetype=factor(format(x, big.mark=",",
# trim=TRUE), levels = format(sort(unique(x)), big.mark=",", trim=TRUE)),
# color = deltax)) +
# geom_path(size = 1) + geom_point(data = realdat[realdat$t > 200*60,], color = "black") +
# facet_grid(deltax~.) +scale_y_continuous(expression(flow~~(ft^3~~s^-1)),
# breaks=c(300, 500,700)) +
# scale_linetype_manual("downstream\ndistance (ft)", values = c("solid",
# "dashed", "dotted")) + plotopts +
# scale_x_continuous(expression(time~~(min)), limits=c(0,1250)) +
# scale_color_manual(values = c("#8DA0CB", "#66C2A5", "#A6D854", "#E5C494",
# "#E6AB02"), guide = FALSE)
#
# ## ---- peakstable-characteristic ---------------------------------------
# peaks = as.data.frame(summarize(group_by(filter(allresults, x==50000,
# monitor.type=='node'), computationtime, deltax, deltat, monitor.type),
# peak.flow=max(Q), time.to.peak=t[which(Q==peak.flow)]))
# real.peak.flow = max(waterolympics$Q)
# real.time.to.peak = waterolympics$t[which(waterolympics$Q==real.peak.flow)][1]
# peaks['flow.percent.error'] =
# (peaks$peak.flow - real.peak.flow)/real.peak.flow
# peaks['time.percent.error'] =
# (peaks$time.to.peak - real.time.to.peak)/real.time.to.peak
# tbld = as.data.frame(peaks[c("deltax", "deltat", "flow.percent.error",
# "time.percent.error", "computationtime")])
# tbld["flow.percent.error"] = round(100*tbld$flow.percent.error, 2)
# tbld["time.percent.error"] = round(100*tbld$time.percent.error, 2)
# names(tbld) = c("$\\Delta x$ (ft)", "$\\Delta t$ (s)", "\\% error (peak flow)",
# "\\% error (time to peak)", "cost (s)")
# ptable = xtable(tbld[order(tbld[,1]),], digits = c(0,0,2,2,2,2))
# print(ptable, include.rownames=FALSE, sanitize.colnames.function=identity,
# sanitize.text.function = function(x){x}, floating=FALSE, hline.after=NULL,
# add.to.row=list(pos=list(-1,0, nrow(ptable)),
# command=c('\\toprule ', '\\midrule ', '\\bottomrule ')),
# format.args = list(big.mark = ","))
#
# ## ---- loadup-boundaries -----------------------------------------------
# oldscipen = options('scipen')
# options(scipen = 1000)
#
# plotopts = list(theme_bw())
# slope = 0.00008
# extent = 5000
# mannings = 0.013
# B = 6.1
# SS = 1.5
# g = 9.81
# Cm = 1
#
# ic = 126
# id = 5.79
# ia = channel_geom(id, B, SS)[["A"]]
# CN = 0.9
#
# dx = 10
# dt = round(dx*CN/(ic/ia + sqrt(id*g)), 2)
# numnodes = extent/dx + 1
#
# bctime = 2000
# bc = rep(id, round(bctime/dt) + 1)
# dc = rep(0, length(bc))
# dt = round(bctime/(length(bc) - 1), 2)
#
# CN = dt*(ic/ia + sqrt(id*g))/dx
#
# mp = c(1, as.integer(c(1500, 2500, 3000, 5000)/dx + 1))
# mt = as.integer(round(seq(0, length(bc)-1, by=25))) + 1L
# mtslice = c(1, as.integer(c(500, 1000, 1500, 2000)/dt + 1))
#
# dclose.lax = route_wave(slope, mannings, Cm, g, B, SS, ic, bc, dc,
# timestep=dt, spacestep=dx, numnodes=numnodes, monitor.nodes=mp,
# monitor.times=mt, engine="Dynamic", scheme="Lax", boundary.type="yQ")
# dclose.mac = route_wave(slope, mannings, Cm, g, B, SS, ic, bc, dc,
# timestep=dt, spacestep=dx, numnodes=numnodes, monitor.nodes=mp,
# monitor.times=mt, engine="Dynamic", scheme="MacCormack", boundary.type="yQ")
# dclose.lax['scheme'] = "Lax diffusive"
# dclose.mac['scheme'] = "MacCormack predictor-corrector"
# dclose = rbind(dclose.lax, dclose.mac)
# row.names(dclose) = NULL
# dclose["CN"] = dt*(dclose$flow/dclose$area + sqrt(dclose$depth*g))/dx
#
# ## ---- through-time ----------------------------------------------------
# dclose.times = filter(dclose, monitor.type=="timestep")
# ggplot(filter(dclose.times, step %in% mtslice),
# aes(x=distance, y=depth, linetype = factor(round(time)),
# color = factor(round(time)))) + geom_path(size = 1) +
# scale_y_continuous(expression(depth~~(ft))) + facet_wrap(~scheme) +
# scale_x_continuous(expression(distance~~downstream~~(ft))) + plotopts +
# scale_color_manual("time (s)", values = brewer.pal(7, "YlGnBu")[3:7]) +
# scale_linetype_manual("time (s)", values = c("solid", "longdash", "dashed",
# "dotdash", "dotted"))
# # animation
# routeylims = c(min(dclose.times$depth), max(dclose.times$depth))
# routexlims = c(min(dclose.times$distance), max(dclose.times$distance))
# for(n in unique(dclose.times$step))
# print(
# ggplot(filter(dclose.times, step==n), aes(x=distance, y=depth)) +
# geom_line(color = "darkblue", size = 1) + plotopts +
# scale_linetype_discrete("time\n(minutes)") +
# scale_x_continuous(expression(distance~~downstream~~(ft)),
# limits=routexlims) +
# scale_y_continuous(expression(depth~~(ft)), limits=routeylims) +
# ggtitle(paste(format(round((n - 1)*dt), width = 4), "seconds")) +
# facet_wrap(~scheme)
# )
Try the rivr package in your browser
Any scripts or data that you put into this service are public.
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.