Streamflow - code
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import seaborn as sns
sns.set(style="ticks", font_scale=1.5)
from ipywidgets import *
def linear_reservoir_response(t, p, q0, T):
return p + (q0 - p) * np.exp(-t/T)
def hydrograph_response(t, t0, p, T, Tp):
response = np.zeros_like(t)
# time for response starts at t0
index_t0 = np.where(t == t0)[0][0]
# we shift the time, because linear_reservoir_response assumes t starts at zero
t_response = t[index_t0:] - t0
# calculate rise part
rise = linear_reservoir_response(t_response, p, 0, T)
# calculate fall only if rainfall events ends before time array
if Tp < t[-1]:
# find index of Tp
index_Tp = np.where(t_response == Tp)[0][0]
# calculate initial condition of fall
q0 = rise[index_Tp]
# compute fall
fall = linear_reservoir_response(t_response, 0, q0, T)
# we save fall values in the indices of rise after Tp
rise[index_Tp:] = fall[:len(t_response)-index_Tp]
# save full response (stored on rise array) on response array, after t0=beginning of rain
response[index_t0:] = rise
return response
%matplotlib notebook
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,7))
t_bars = [0]
p_bars = [1]
ax1.bar(t_bars, p_bars, width=1, align='edge', linewidth=0, alpha=0.5)
t = np.arange(0, 10, 0.01)
T_watershed = 1.0
Tp = 1.0
response = np.zeros_like(t)
for i, prec in enumerate(p_bars):
h = hydrograph_response(t, t_bars[i], p_bars[i], T_watershed, 1)
response = response + h
# ax2.plot(t, h, color="black", alpha=0.3, ls="--", zorder=2)
ax2.plot(t, response, color="tab:red", lw=2, zorder=1)
ax1.set(ylabel=r"water input (L$^3$ T$^{-1}$)",
xlim=[-0.2, t[-1]],
ylim=[0, 1.1],
yticks=[0,0.5,1])
ax2.set(xlabel="time (h)",
ylabel="stream response (L$^3$ T$^{-1}$)",
xlim=[-0.2, t[-1]],
ylim=[0, 1.1],
yticks=[0,0.5,1])
ax1.text(0.99, 0.97, "hyetograph", transform=ax1.transAxes, ha='right', va='top', fontsize=20)
ax2.text(0.99, 0.97, "hydrograph", transform=ax2.transAxes, ha='right', va='top', fontsize=20)
height = response.max()
# ax2.annotate("",
# xy=(1, height*1.1), #xycoords='data',
# xytext=(2, height*1.1),# textcoords='data',
# size=6,
# arrowprops=dict(arrowstyle="|-|",
# shrinkA=0, shrinkB=0,
# connectionstyle="arc3",
# color='black'),
# )
# ax2.text(1.5, height*1.2, "T* = 1 h", ha="center", fontsize=16)
# ax2.annotate("",
# xy=(2, height), #xycoords='data',
# xytext=(2, height*(1.0-np.exp(-1))),# textcoords='data',
# size=6,
# arrowprops=dict(arrowstyle="|-|",
# shrinkA=0, shrinkB=0,
# connectionstyle="arc3",
# color='black'),
# )
# ax2.text(2.1, height*(1.0-np.exp(-1)/2), r"decay of 37% = exp($-1$)", va="center", fontsize=16)
ax1.set_title("response of linear watershed, T* = 1 h")
plt.savefig("hydrology_figures/hyetograph_hydrograph1.png")
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,7))
t_bars = [0]
p_bars = [1]
ax1.bar(t_bars, p_bars, width=1, align='edge', linewidth=0, alpha=0.5)
t = np.arange(0, 10, 0.01)
T_watershed = 2.0
Tp = 1.0
response = np.zeros_like(t)
for i, prec in enumerate(p_bars):
h = hydrograph_response(t, t_bars[i], p_bars[i], T_watershed, 1)
response = response + h
# ax2.plot(t, h, color="black", alpha=0.3, ls="--", zorder=2)
ax2.plot(t, response, color="tab:red", lw=2, zorder=1)
ax1.set(ylabel=r"water input (L$^3$ T$^{-1}$)",
xlim=[-0.2, t[-1]],
ylim=[0, 1.1],
yticks=[0,0.5,1])
ax2.set(xlabel="time (h)",
ylabel="stream response (L$^3$ T$^{-1}$)",
xlim=[-0.2, t[-1]],
ylim=[0, 1.1],
yticks=[0,0.5,1])
ax1.text(0.99, 0.97, "hyetograph", transform=ax1.transAxes, ha='right', va='top', fontsize=20)
ax2.text(0.99, 0.97, "hydrograph", transform=ax2.transAxes, ha='right', va='top', fontsize=20)
height = response.max()
ax2.annotate("",
xy=(1, height*1.15), #xycoords='data',
xytext=(3, height*1.15),# textcoords='data',
size=6,
arrowprops=dict(arrowstyle="|-|",
shrinkA=0, shrinkB=0,
connectionstyle="arc3",
color='black'),
)
ax2.text(2, height*1.2, "T* = 2 h", ha="center", fontsize=16)
ax2.annotate("",
xy=(3, height), #xycoords='data',
xytext=(3, height*(1.0-np.exp(-1))),# textcoords='data',
size=6,
arrowprops=dict(arrowstyle="|-|",
shrinkA=0, shrinkB=0,
connectionstyle="arc3",
color='black'),
)
ax2.text(3.1, height*(1.0-np.exp(-1)/2), r"decay of 37% = exp($-1$)", va="center", fontsize=16)
ax1.set_title("response of linear watershed, T* = 2 h")
plt.savefig("hydrology_figures/hyetograph_hydrograph2.png")
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,7))
t_bars = [0, 1, 2, 4]
p_bars = [3, 1, 2, 1]
ax1.bar(t_bars, p_bars, width=1, align='edge', linewidth=0, alpha=0.5)
t = np.arange(0, 10, 0.01)
T_watershed = 1.0
Tp = 1.0
response = np.zeros_like(t)
for i, prec in enumerate(p_bars):
h = hydrograph_response(t, t_bars[i], p_bars[i], T_watershed, 1)
response = response + h
ax2.plot(t, h, color="black", alpha=0.3, ls="-", zorder=2)
ax2.plot(t, response, color="tab:red", lw=2, zorder=1)
ax1.set(ylabel=r"water input (L$^3$ T$^{-1}$)",
xlim=[-0.2, t[-1]],
ylim=[0, 3.1],
yticks=[0,1,2,3])
ax2.set(xlabel="time (h)",
ylabel="stream response (L$^3$ T$^{-1}$)",
xlim=[-0.2, t[-1]],
ylim=[0, 3.1],
yticks=[0,1,2,3])
ax1.text(0.99, 0.97, "hyetograph", transform=ax1.transAxes, ha='right', va='top', fontsize=20)
ax2.text(0.99, 0.97, "hydrograph", transform=ax2.transAxes, ha='right', va='top', fontsize=20)
ax1.set_title("response of linear watershed, T* = 1 h")
plt.savefig("hydrology_figures/hyetograph_hydrograph3b.png")
%matplotlib widget
a = 1.0
b = 1.0
c = 1.0
x = np.linspace(0, 2 * np.pi)
def parabola(a):
return a*x**2 + b*x + c
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
line, = ax.plot(x, parabola(x))
def update_graph(a = 1.0):
line.set_ydata(parabola(a))
fig.canvas.draw_idle()
interact(update_graph);
%matplotlib widget
t_bars = [0, 1, 2, 4]
p_bars = [3, 1, 2, 1]
t = np.arange(0, 10, 0.01)
T_watershed = 1.0
Tp = 1.0
def resp(T_star):
response = np.zeros_like(t)
for i, prec in enumerate(p_bars):
h = hydrograph_response(t, t_bars[i], p_bars[i], T_star, 1)
response = response + h
return response
fig, ax = plt.subplots(figsize=(10,7))
ax.bar(t_bars, p_bars, width=1, align='edge', linewidth=0, alpha=0.5)
line, = ax.plot(t, resp(1.0), color="tab:red", lw=2, zorder=10)
ax.set(xlabel="time (h)",
ylabel="stream response (L$^3$ T$^{-1}$)",
xlim=[-0.2, t[-1]],
ylim=[0, 3.1],
yticks=[0,1,2,3],
title="response of linear watershed, T* = 1 h")
def update_graph(T_star = 1.0):
line.set_ydata(resp(T_star))
fig.canvas.draw_idle()
ax.set_title(f"response of linear watershed, T* = {T_star:.2f} h")
interact(update_graph, T_star=(0.01,5,0.01));
$$
t_{pc} = \frac{\sum_{i=1}^n P_i^* \cdot t_i}{P^*}
$$
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,7))
t_bars = np.array([0, 1, 2, 3, 4])
p_bars = np.array([1, 3, 2.5, 1.5, 1.8])
ax1.bar(t_bars, p_bars, width=1, align='edge', linewidth=0, alpha=0.5)
ti = t_bars + 0.5
t_pc = np.sum(p_bars*ti) / np.sum(p_bars)
print(t_pc)
ax1.plot([t_pc]*2, ax1.get_ylim(), color="black", ls="--")
ax1.text(t_pc, 0.1, r"$t_{pc}=$"+f"{t_pc:.2f}", fontsize=16)
t = np.arange(0, 16, 0.01)
T_watershed = 1.5
Tp = 1.0
response = np.zeros_like(t)
for i, prec in enumerate(p_bars):
h = hydrograph_response(t, t_bars[i], p_bars[i], T_watershed, 1)
response = response + h
ax2.plot(t, response, color="tab:red", lw=2, zorder=1)
ax1.set(ylabel=r"water input (L$^3$ T$^{-1}$)",
xlim=[-0.2, t[-1]],
ylim=[0, 3.1],
xticks=[],
yticks=[0,1,2,3])
ax2.set(xlabel="time (h)",
ylabel="stream response (L$^3$ T$^{-1}$)",
xlim=[-0.2, t[-1]],
ylim=[0, 3.1],
yticks=[0,1,2,3])
ti = t + 0.5 * (t[1] - t[0])
t_qc = np.sum(response*ti) / np.sum(response)
ax2.plot([t_qc]*2, ax2.get_ylim(), color="black", ls="--")
ax2.text(t_qc, 0.1, r"$t_{qc}=$"+f"{t_qc:.2f}", fontsize=16)
ax1.text(0.99, 0.97, "hyetograph", transform=ax1.transAxes, ha='right', va='top', fontsize=20)
ax2.text(0.99, 0.97, "hydrograph", transform=ax2.transAxes, ha='right', va='top', fontsize=20)
ax1.set_title(f"response of linear watershed, T* = {T_watershed:.1f} h")
# plt.savefig("hydrology_figures/hyetograph_hydrograph3b.png")
np.sum(p_bars)
np.sum(response)*0.01
import pandas as pd
data_file = "cincinnati1.dat"
df = pd.read_csv(data_file,
header=31, # no headers needed, we'll do that later
delim_whitespace=True, # blank spaces separate between columns
na_values=["Bkw"] # substitute these values for missing (NaN) values
)
# headers = pd.read_csv("HEADERS_sub_hourly.txt", # load headers file
# header=1, # skip the first [0] line
# delim_whitespace=True
# )
df.columns = ['agency_cd', 'stie_no','datetime','tz_cd','EDT','discharge','code'] # rename df columns with headers columns
# # LST = local standard time
# df["LST_TIME"] = [f"{x:04d}" for x in df["LST_TIME"]] # time needs padding of zeros, then convert to string
# df['LST_DATE'] = df['LST_DATE'].astype(str) # convert date into string
df['date_and_time'] = df['datetime'] + ' ' + df['tz_cd'] # combine date+time into datetime
df['date_and_time'] = pd.to_datetime(df['date_and_time']) # interpret datetime
df = df.set_index('date_and_time') # make datetime the index
df['discharge'] = df['discharge'].astype(float)
df['discharge'] = df['discharge'] * 0.0283168 # convert cubic feet to m3
# https://waterdata.usgs.gov/nwis/inventory?agency_code=USGS&site_no=03260015
%matplotlib notebook
plt.plot(df['discharge'])
two_points = df.loc[[pd.to_datetime("2021-03-18 02:00:00"),
pd.to_datetime("2021-03-18 22:00:00")]]['discharge']
plt.plot(two_points, 'o')
new = pd.DataFrame(data=two_points, index=two_points.index)
df_linear = (new.resample("15min") #resample hourly
.interpolate(method='time') #interpolate by time
)
# new.resample('H').interpolate(method='time')
plt.plot(df_linear)
qstar = df.loc[df_linear.index]['discharge'] - df_linear['discharge']
plt.plot(qstar)
plt.gcf().autofmt_xdate()
# plt.xlim(['2021-03-17 18:00:00', '2021-03-19 04:00:00'])
plt.xlim((pd.to_datetime("2021-03-17 18:00:00"),
pd.to_datetime("2021-03-19 04:00:00"))
)
plt.ylim([0, 130])
from matplotlib.dates import HourLocator, DateFormatter
import matplotlib.dates as mdates
import matplotlib.ticker as ticker
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,7))
fig.subplots_adjust(wspace=0.05)
ax1.plot(df['discharge'], color="black", lw=2)
two_points = df.loc[[pd.to_datetime("2021-03-18 02:00:00"),
pd.to_datetime("2021-03-18 22:00:00")]]['discharge']
ax1.plot(two_points, 'o', color="tab:red")
new = pd.DataFrame(data=two_points, index=two_points.index)
df_linear = (new.resample("15min") #resample hourly
.interpolate(method='time') #interpolate by time
)
ax1.plot(df_linear)
df_between_2_points = df.loc[df_linear.index]
ax1.fill_between(df_between_2_points.index, df_between_2_points['discharge'],
y2=df_linear['discharge'],
color="tab:blue", alpha=0.3)
ax1.annotate("event flow",
xy=(pd.to_datetime("2021-03-18 13:00:00"), 1.4), #xycoords='data',
xytext=(pd.to_datetime("2021-03-18 16:00:00"), 2.2),# textcoords='data',
arrowprops=dict(arrowstyle="->",
color='black'),
)
ax1.annotate("base flow",
xy=(pd.to_datetime("2021-03-18 17:00:00"), 0.1), #xycoords='data',
xytext=(pd.to_datetime("2021-03-18 17:00:00"), 1.1),# textcoords='data',
arrowprops=dict(arrowstyle="->",
color='black'),
)
qstar = df.loc[df_linear.index]['discharge'] - df_linear['discharge']
ax2.plot(qstar, color="black", lw=2)
ax2.fill_between(qstar.index, qstar,
y2=0.0,
color="tab:blue", alpha=0.3)
ax1.set(xlim=[pd.to_datetime("2021-03-17 18:00:00"),
pd.to_datetime("2021-03-19 04:00:00")],
ylabel=r"discharge (m$^3$/s)",
ylim=[0, 4],
yticks=[0,1,2,3,4])
ax2.set(yticks=[],
ylim=[0, 4],
xlim=[pd.to_datetime("2021-03-17 18:00:00"),
pd.to_datetime("2021-03-19 04:00:00")],
)
ax1.text(0.5, 0.97, "total flow rate = event + base\n" + r"$q(t) = q\!^*\!(t) + q_{bf}\,(t)$",
transform=ax1.transAxes, ha='center', va='top', fontsize=16)
ax2.text(0.5, 0.97, r"event flow rate, $q\!^*\!(t)$",
transform=ax2.transAxes, ha='center', va='top', fontsize=16)
ax1.xaxis.set_major_locator(HourLocator(byhour=[0, 8, 16], interval=1))
ax2.xaxis.set_major_locator(HourLocator(byhour=[0, 8, 16], interval=1))
plt.gcf().autofmt_xdate()
# dirty trick to get common title between the two panels:
# make a large invisible axes, give it a title
ax0 = fig.add_subplot(111, frame_on=False)
ax0.tick_params(labelcolor="none", bottom=False, left=False)
ax0.set_title("Pleasant Run Creek at Oak St. near Ludlow, KY\nMarch 2021, USGS 03260015 ")
data_file = "USGS 03337100 BONEYARD CREEK AT LINCOLN AVE AT URBANA, IL.dat"
df_q = pd.read_csv(data_file,
header=31, # no headers needed, we'll do that later
delim_whitespace=True, # blank spaces separate between columns
na_values=["Bkw"] # substitute these values for missing (NaN) values
)
df_q.columns = ['agency_cd', 'site_no','datetime','tz_cd','EDT','discharge','code'] # rename df columns with headers columns
df_q['date_and_time'] = df_q['datetime'] + ' ' + df_q['tz_cd'] # combine date+time into datetime
df_q['date_and_time'] = pd.to_datetime(df_q['date_and_time']) # interpret datetime
df_q = df_q.set_index('date_and_time') # make datetime the index
df_q['discharge'] = df_q['discharge'].astype(float)
df_q['discharge'] = df_q['discharge'] * 0.0283168 # convert cubic feet to m3
data_file = "Champaign - IL.txt"
df_p = pd.read_csv(data_file,
header=None, # no headers needed, we'll do that later
delim_whitespace=True, # blank spaces separate between columns
na_values=["-99.000", "-9999.0"] # substitute these values for missing (NaN) values
)
headers = pd.read_csv("HEADERS_sub_hourly.txt", # load headers file
header=1, # skip the first [0] line
delim_whitespace=True
)
df_p.columns = headers.columns # rename df columns with headers columns
# LST = local standard time
df_p["LST_TIME"] = [f"{x:04d}" for x in df_p["LST_TIME"]] # time needs padding of zeros, then convert to string
df_p['LST_DATE'] = df_p['LST_DATE'].astype(str) # convert date into string
df_p['datetime'] = df_p['LST_DATE'] + ' ' + df_p['LST_TIME'] # combine date+time into datetime
df_p['datetime'] = pd.to_datetime(df_p['datetime']) # interpret datetime
df_p = df_p.set_index('datetime') # make datetime the index
# plt.plot(df_p['PRECIPITATION'])
%matplotlib notebook
# Drainage area: 4.78 square miles
area = 4.78 / 0.00000038610 # squared miles to squared meters
start = "2020-10-20 14:00:00"
end = "2020-10-21 04:00:00"
df_p_small = df_p.loc[start:end].copy()
df_q_small = df_q.loc[start:end].copy()
df_p_small['PRECIPITATION'] = df_p_small['PRECIPITATION'].values * area / 1000 # mm to m3 in the whole watershed
total_p = df_p_small['PRECIPITATION'].sum()
df_p_small['PRECIPITATION'] = df_p_small['PRECIPITATION'] / 60 / 5 # convert m3 per 5 min to m3/s
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,7))
fig.subplots_adjust(hspace=0.05)
ax1.bar(df_p_small.index, df_p_small['PRECIPITATION'], width=1/24/12)
ax2.plot(df_q_small['discharge'], color="tab:blue", lw=2)
ax1.set_xlim([pd.to_datetime("2020-10-20 14:00:00"),
pd.to_datetime("2020-10-21 04:00:00")])
ax2.set_xlim([pd.to_datetime("2020-10-20 14:00:00"),
pd.to_datetime("2020-10-21 04:00:00")])
# total_p = df_p_small['PRECIPITATION'].sum()
total_q = df_q_small['discharge'].sum() * 60*5
# dirty trick to get common title between the two panels:
# make a large invisible axes, give it a title
ax0 = fig.add_subplot(111, frame_on=False)
ax0.tick_params(labelcolor="none", bottom=False, left=False)
ax0.set_ylabel(r"flow (m$^3$/s)")
ax1.text(0.99, 0.97, r"precipitation $p(t)$, hyetograph" + "\n" + f"P = {total_p:.2e}" + r" m$^3$",
transform=ax1.transAxes, ha='right', va='top', fontsize=20)
ax2.text(0.99, 0.97, r"discharge $q(t)$, hydrograph" + "\n" + f"Q = {total_q:.2e}" + r" m$^3$",
transform=ax2.transAxes, ha='right', va='top', fontsize=20)
ax1.set_xticks([])
locator = mdates.AutoDateLocator(minticks=3, maxticks=7)
formatter = mdates.ConciseDateFormatter(locator)
ax2.xaxis.set_major_locator(locator)
ax2.xaxis.set_major_formatter(formatter)
ax2.set_ylim([0,5.5])
ax1.set_title("Precipitation and discharge, Boneyard Creek at Urbana, IL\n 20-21 October 2020, 5-minute data")
fig.savefig("hydrology_figures/urbana_pq.png")
total_q/ total_p
from matplotlib.dates import HourLocator, DateFormatter
import matplotlib.dates as mdates
import matplotlib.ticker as ticker
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,7))
fig.subplots_adjust(wspace=0.05)
ax1.plot(df_q_small['discharge'], color="black", lw=2)
point1 = pd.to_datetime("2020-10-20 16:40:00")
point2 = pd.to_datetime("2020-10-21 00:00:00")
# two_points = df_q_small.loc[["2020-10-20 16:40:00", "2020-10-20 23:00:00"]]
two_points = df_q_small.loc[[point1, point2]]['discharge']
ax1.plot(two_points, 'o', color="tab:red")
new = pd.DataFrame(data=two_points, index=two_points.index)
df_linear = (new.resample("5min") #resample
.interpolate(method='time') #interpolate by time
)
ax1.plot(df_linear, color="tab:blue")
df_between_2_points = df_q_small.loc[df_linear.index]
ax1.fill_between(df_between_2_points.index, df_between_2_points['discharge'],
y2=df_linear['discharge'],
color="tab:blue", alpha=0.3)
ax1.annotate("event flow",
xy=(pd.to_datetime("2020-10-20 20:00:00"), 1.0), #xycoords='data',
xytext=(pd.to_datetime("2020-10-20 22:00:00"), 2.2),# textcoords='data',
size=16,
arrowprops=dict(arrowstyle="->",
color='black'),
)
ax1.annotate("base flow",
xy=(pd.to_datetime("2020-10-21 00:00:00"), 0.1), #xycoords='data',
xytext=(pd.to_datetime("2020-10-21 0:00:00"), 1.1),# textcoords='data',
size=16,
arrowprops=dict(arrowstyle="->",
color='black'),
)
qstar = df_q_small.loc[df_linear.index]['discharge'] - df_linear['discharge']
total_qstar = qstar.sum() * 60 * 5
ax2.plot(qstar, color="black", lw=2)
ax2.fill_between(qstar.index, qstar,
y2=0.0,
color="tab:blue", alpha=0.3)
ax1.set(xlim=[df_q_small.index[0],
df_q_small.index[-1]],
ylabel=r"discharge (m$^3$/s)",
ylim=[0, 5.5],
yticks=[0,1,2,3,4])
ax2.set(yticks=[],
ylim=[0, 5.5],
xlim=[df_q_small.index[0],
df_q_small.index[-1]],
)
ax1.text(0.5, 0.97, "total flow rate = event + base\n" + r"$q(t) = q\!^*\!(t) + q_{bf}\,(t)$",
transform=ax1.transAxes, ha='center', va='top', fontsize=16)
ax2.text(0.5, 0.97, r"event flow rate, $q\!^*\!(t)$",
transform=ax2.transAxes, ha='center', va='top', fontsize=16)
ax2.annotate(f"Q* = {total_qstar:.0f}" + r" m$^3$",
xy=(pd.to_datetime("2020-10-20 20:00:00"), 1.0), #xycoords='data',
xytext=(pd.to_datetime("2020-10-20 21:00:00"), 3),# textcoords='data',
size=16,
arrowprops=dict(arrowstyle="->",
color='black'),
)
locator = mdates.AutoDateLocator(minticks=3, maxticks=7)
formatter = mdates.ConciseDateFormatter(locator)
ax1.xaxis.set_major_locator(locator)
ax1.xaxis.set_major_formatter(formatter)
ax2.xaxis.set_major_locator(locator)
ax2.xaxis.set_major_formatter(formatter)
# dirty trick to get common title between the two panels:
# make a large invisible axes, give it a title
ax0 = fig.add_subplot(111, frame_on=False)
ax0.tick_params(labelcolor="none", bottom=False, left=False)
ax0.set_title("Boneyard Creek at Urbana, IL, USGS 03337100\n 20-21 October 2020, 5-minute data")
fig.savefig("hydrology_figures/urbana_q_qstar.png")
%matplotlib notebook
# Drainage area: 4.78 square miles
# area = 4.78 / 0.00000038610 # squared miles to squared meters
# start = "2020-10-20 14:00:00"
# end = "2020-10-21 04:00:00"
# df_p_small = df_p.loc[start:end].copy()
# df_q_small = df_q.loc[start:end].copy()
# df_p_small['PRECIPITATION'] = df_p_small['PRECIPITATION'].values * area / 1000 # mm to m3 in the whole watershed
# total_p = df_p_small['PRECIPITATION'].sum()
# df_p_small['PRECIPITATION'] = df_p_small['PRECIPITATION'] / 60 / 5 # convert m3 per 5 min to m3/s
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,7))
fig.subplots_adjust(hspace=0.05)
ratio = total_qstar/ total_p
pstar = df_p_small['PRECIPITATION'] * ratio
total_pstar = pstar.sum() * 5 * 60
ax1.bar(pstar.index, pstar, width=1/24/12)
ax2.plot(qstar, color="tab:blue", lw=2)
print(total_pstar, total_qstar, total_qstar/total_pstar)
ax1.set_xlim([pd.to_datetime("2020-10-20 14:00:00"),
pd.to_datetime("2020-10-21 04:00:00")])
ax2.set_xlim([pd.to_datetime("2020-10-20 14:00:00"),
pd.to_datetime("2020-10-21 04:00:00")])
# # total_p = df_p_small['PRECIPITATION'].sum()
# total_q = df_q_small['discharge'].sum() * 60*5
# dirty trick to get common title between the two panels:
# make a large invisible axes, give it a title
ax0 = fig.add_subplot(111, frame_on=False)
ax0.tick_params(labelcolor="none", bottom=False, left=False)
ax0.set_ylabel(r"flow (m$^3$/s)")
ax1.text(0.99, 0.97, r"effective precipitation $p\!*\!(t)$" + "\n" + f"P* = {total_pstar:.2e}" + r" m$^3$",
transform=ax1.transAxes, ha='right', va='top', fontsize=20)
ax2.text(0.99, 0.97, r"effective discharge $q\!*\!(t)$" + "\n" + f"Q* = {total_qstar:.2e}" + r" m$^3$",
transform=ax2.transAxes, ha='right', va='top', fontsize=20)
ax1.set_xticks([])
locator = mdates.AutoDateLocator(minticks=3, maxticks=7)
formatter = mdates.ConciseDateFormatter(locator)
ax2.xaxis.set_major_locator(locator)
ax2.xaxis.set_major_formatter(formatter)
ax2.set_ylim([0,5.5])
ax1.set_title("Effective precipitation and discharge, Boneyard Creek at Urbana, IL\n 20-21 October 2020, 5-minute data")
fig.savefig("hydrology_figures/urbana_pstar_qstar.png")
t0 = pstar[pstar != 0.0].index[0]
tf = pstar[pstar != 0.0].index[-1]
td = (tf-t0) / pd.Timedelta('1 min') # time difference in minutes
time = np.arange(0, td+1, 5) + 2.5
pi = pstar.loc[(pstar.index >= t0) & (pstar.index <= tf)]
pi = pi.values * 60 * 5
t_pc = (pi * time).sum() / pi.sum()
t_pc = t0 + pd.Timedelta(minutes=t_pc)
t_pc
# qstar centroid
t0 = qstar[qstar != 0.0].index[0]
tf = qstar[pstar != 0.0].index[-1]
td = (tf-t0) / pd.Timedelta('1 min') # time difference in minutes
time = np.arange(0, td+1, 5) + 2.5
qi = qstar.loc[(qstar.index >= t0) & (qstar.index <= tf)]
qi = qi.values * 60 * 5
t_qc = (qi * time).sum() / qi.sum()
t_qc = t0 + pd.Timedelta(minutes=t_qc)
t_qc
# time of peak discharge
max_discharge = qstar.max()
t_pk = qstar[qstar == max_discharge].index[0]
qstar
%matplotlib notebook
# Drainage area: 4.78 square miles
# area = 4.78 / 0.00000038610 # squared miles to squared meters
# start = "2020-10-20 14:00:00"
# end = "2020-10-21 04:00:00"
# df_p_small = df_p.loc[start:end].copy()
# df_q_small = df_q.loc[start:end].copy()
# df_p_small['PRECIPITATION'] = df_p_small['PRECIPITATION'].values * area / 1000 # mm to m3 in the whole watershed
# total_p = df_p_small['PRECIPITATION'].sum()
# df_p_small['PRECIPITATION'] = df_p_small['PRECIPITATION'] / 60 / 5 # convert m3 per 5 min to m3/s
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,7))
fig.subplots_adjust(hspace=0.05)
ratio = total_qstar/ total_p
pstar = df_p_small['PRECIPITATION'] * ratio
total_pstar = pstar.sum() * 5 * 60
ax1.bar(pstar.index, pstar, width=1/24/12)
ax2.plot(qstar, color="tab:blue", lw=2)
print(total_pstar, total_qstar, total_qstar/total_pstar)
ax1.set_xlim([pd.to_datetime("2020-10-20 14:00:00"),
pd.to_datetime("2020-10-21 04:00:00")])
ax2.set_xlim([pd.to_datetime("2020-10-20 14:00:00"),
pd.to_datetime("2020-10-21 04:00:00")])
# # total_p = df_p_small['PRECIPITATION'].sum()
# total_q = df_q_small['discharge'].sum() * 60*5
# dirty trick to get common title between the two panels:
# make a large invisible axes, give it a title
ax0 = fig.add_subplot(111, frame_on=False)
ax0.tick_params(labelcolor="none", bottom=False, left=False)
ax0.set_ylabel(r"flow (m$^3$/s)")
ax1.text(0.99, 0.97, r"effective precipitation $p\!*\!(t)$" + "\n" + f"P* = {total_pstar:.2e}" + r" m$^3$",
transform=ax1.transAxes, ha='right', va='top', fontsize=20)
ax2.text(0.99, 0.97, r"effective discharge $q\!*\!(t)$" + "\n" + f"Q* = {total_qstar:.2e}" + r" m$^3$",
transform=ax2.transAxes, ha='right', va='top', fontsize=20)
ax1.set_xticks([])
locator = mdates.AutoDateLocator(minticks=3, maxticks=7)
formatter = mdates.ConciseDateFormatter(locator)
ax2.xaxis.set_major_locator(locator)
ax2.xaxis.set_major_formatter(formatter)
ax2.set_ylim([0,5.5])
ax1.set_title("precipitation centroid and discharge centroid, Boneyard Creek at Urbana, IL\n 20-21 October 2020, 5-minute data")
ax1.plot([t_pc, t_pc], ax1.get_ylim(), color="tab:red", lw=2)
ax2.plot([t_qc, t_qc], ax2.get_ylim(), color="tab:red", lw=2)
ax2.plot([t_pk, t_pk], ax2.get_ylim(), color="tab:olive", lw=2)
ax1.text(t_pc, ax1.get_ylim()[1], r"$t_{pc}=$"+f"{t_pc.hour}:{t_pc.minute}", ha="right", va="top", fontsize=16)
ax2.text(t_qc, ax2.get_ylim()[0], r"$t_{qc}=$"+f"{t_qc.hour}:{t_qc.minute}", ha="right", va="bottom", fontsize=16)
ax2.text(t_pk, ax2.get_ylim()[0], r"$t_{pk}=$"+f"{t_pk.hour}:{t_pk.minute}", ha="left", va="bottom", fontsize=16)
ax2.plot([t_pc, t_pc], ax2.get_ylim(), color="black", lw=2, ls=':')
ax2.annotate("",
xy=(t_pc, ax2.get_ylim()[1]*0.9), #xycoords='data',
xytext=(t_qc, ax2.get_ylim()[1]*0.9),# textcoords='data',
size=16,
arrowprops=dict(arrowstyle="<->",
shrinkA=0, shrinkB=0,
connectionstyle="arc3",
color='black'),
)
ax2.annotate("",
xy=(t_pc, ax2.get_ylim()[1]*0.7), #xycoords='data',
xytext=(t_pk, ax2.get_ylim()[1]*0.7),# textcoords='data',
size=16,
arrowprops=dict(arrowstyle="<->",
shrinkA=0, shrinkB=0,
connectionstyle="arc3",
color='black'),
)
centroid_lag = t_qc - t_pc
centroid_lag_to_peak = t_pk - t_pc
centroid_lag_minutes = int(centroid_lag.total_seconds() / 60)
centroid_lag_to_peak_minutes = int(centroid_lag_to_peak.total_seconds() / 60)
ax2.text(t_pc, ax2.get_ylim()[1]*0.9, r"$T_{LC}=$"+f"{centroid_lag_minutes} minutes", ha="right", va="center", fontsize=16)
ax2.text(t_pc, ax2.get_ylim()[1]*0.7, r"$T_{LPC}=$"+f"{centroid_lag_to_peak_minutes} minutes", ha="right", va="center", fontsize=16)
fig.savefig("hydrology_figures/urbana_lags.png")
t_pk[0]
t_qc.hour