import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import seaborn as sns
set(style="ticks", font_scale=1.5)
sns.from ipywidgets import *
14 Exercises
Import relevant packages
Import streamflow data from USGS’s National Water Information System. We will be using data from Urbana, IL.
# Drainage area: 4.78 square miles
= "USGS 03337100 BONEYARD CREEK AT LINCOLN AVE AT URBANA, IL.dat"
data_file = pd.read_csv(data_file,
df_q_2020 =31, # no headers needed, we'll do that later
header=True, # blank spaces separate between columns
delim_whitespace=["Bkw"] # substitute these values for missing (NaN) values
na_values
)= ['agency_cd', 'site_no','datetime','tz_cd','EDT','discharge','code'] # rename df columns with headers columns
df_q_2020.columns 'date_and_time'] = df_q_2020['datetime'] + ' ' + df_q_2020['tz_cd'] # combine date+time into datetime
df_q_2020['date_and_time'] = pd.to_datetime(df_q_2020['date_and_time']) # interpret datetime
df_q_2020[= df_q_2020.set_index('date_and_time') # make datetime the index
df_q_2020 'discharge'] = df_q_2020['discharge'].astype(float)
df_q_2020['discharge'] = df_q_2020['discharge'] * 0.0283168 # convert cubic feet to m3
df_q_2020[
= plt.subplots(figsize=(10,7))
fig, ax 'discharge'], '-o')
ax.plot(df_q_2020[
plt.gcf().autofmt_xdate()set(xlabel="date",
ax.=r"discharge (m$^3$/5min)"); ylabel
Import sub-hourly (5-min) rainfall data from NOAA’s Climate Reference Network Data website
= "Champaign - IL.txt"
data_file = pd.read_csv(data_file,
df_p_2020 =None, # no headers needed, we'll do that later
header=True, # blank spaces separate between columns
delim_whitespace=["-99.000", "-9999.0"] # substitute these values for missing (NaN) values
na_values
)= pd.read_csv("HEADERS_sub_hourly.txt", # load headers file
headers =1, # skip the first [0] line
header=True
delim_whitespace
)= headers.columns # rename df columns with headers columns
df_p_2020.columns # LST = local standard time
"LST_TIME"] = [f"{x:04d}" for x in df_p_2020["LST_TIME"]] # time needs padding of zeros, then convert to string
df_p_2020['LST_DATE'] = df_p_2020['LST_DATE'].astype(str) # convert date into string
df_p_2020['datetime'] = df_p_2020['LST_DATE'] + ' ' + df_p_2020['LST_TIME'] # combine date+time into datetime
df_p_2020['datetime'] = pd.to_datetime(df_p_2020['datetime']) # interpret datetime
df_p_2020[= df_p_2020.set_index('datetime') # make datetime the index df_p_2020
Plot rainfall and streamflow. Does this makes sense?
= plt.subplots(2, 1, figsize=(10,7))
fig, (ax1, ax2) =0.05)
fig.subplots_adjust(hspace
= "2020-10-18"
start = "2020-10-25"
end 'PRECIPITATION'])
ax1.plot(df_p_2020[start:end]['discharge'], color="tab:blue", lw=2)
ax2.plot(df_q_2020[start:end][
set(xticks=[],
ax1.=r"precipitation (mm)")
ylabelset(xlabel="date",
ax2.=r"discharge (m$^3$/5min)")
ylabel
# makes slated dates plt.gcf().autofmt_xdate()
Define smaller dataframes for \(p(t)\) and \(q(t)\), between the dates:
= "2020-10-20 14:00:00"
start = "2020-10-21 04:00:00" end
Don’t forget to convert the units to SI!
Calculate total rainfall \(P^*\) and total discharge \(Q^*\), in m\(^3\).
# Drainage area: 4.78 square miles
= 4.78 / 0.00000038610 # squared miles to squared meters
area = "2020-10-20 14:00:00"
start = "2020-10-21 04:00:00"
end
= df_p_2020.loc[start:end]['PRECIPITATION'].to_frame()
df_p = df_p_2020.loc[start:end]['PRECIPITATION'].to_frame()
df_p_mm = df_q_2020.loc[start:end]['discharge'].to_frame()
df_q
'PRECIPITATION'] = df_p['PRECIPITATION'].values * area / 1000 # mm to m3 in the whole watershed
df_p['PRECIPITATION'] = df_p['PRECIPITATION'] / 60 / 5 # convert m3 per 5 min to m3/s
df_p[
= df_p['PRECIPITATION'].sum() * 60 * 5
P = df_q['discharge'].sum() * 60 * 5
Q
print("total precipitation during event: Pstar = {:.1e} m3".format(P.sum()))
print("total streamflow during event: Qstar = {:.1e} m3".format(Q.sum()))
total precipitation during event: Pstar = 2.6e+05 m3
total streamflow during event: Qstar = 5.2e+04 m3
Make another graph of \(p(t)\) and \(q(t)\), now with SI units.
= plt.subplots(2, 1, figsize=(10,7))
fig, (ax1, ax2) =0.05)
fig.subplots_adjust(hspace
= "2020-10-18"
start = "2020-10-25"
end 'PRECIPITATION'])
ax1.plot(df_p['discharge'], color="tab:blue", lw=2)
ax2.plot(df_q[
set(xticks=[],
ax1.=r"precipitation (m$^3$/s)",
ylabel="Precipitation and discharge, Boneyard Creek at Urbana, IL\n 20-21 October 2020, 5-minute data")
titleset(xlabel="date",
ax2.=r"discharge (m$^3$/s)")
ylabel
# makes slated dates plt.gcf().autofmt_xdate()
It’s time for base flow separation! Convert \(q(t)\) into \(q^*(t)\)
from matplotlib.dates import HourLocator, DateFormatter
import matplotlib.dates as mdates
import matplotlib.ticker as ticker
= plt.subplots(1, 2, figsize=(10,7))
fig, (ax1, ax2) =0.05)
fig.subplots_adjust(wspace
'discharge'], color="black", lw=2)
ax1.plot(df_q[= pd.to_datetime("2020-10-20 16:40:00")
point1 = pd.to_datetime("2020-10-21 00:00:00")
point2 = df_q.loc[[point1, point2]]['discharge']
two_points 'o', color="tab:red")
ax1.plot(two_points,
= pd.DataFrame(data=two_points, index=two_points.index)
new
= (new.resample("5min") #resample
df_linear ='time') #interpolate by time
.interpolate(method
)
="tab:blue")
ax1.plot(df_linear, color
= df_q.loc[df_linear.index]
df_between_2_points 'discharge'],
ax1.fill_between(df_between_2_points.index, df_between_2_points[=df_linear['discharge'],
y2="tab:blue", alpha=0.3)
color
= df_q.loc[df_linear.index]['discharge'] - df_linear['discharge']
qstar = qstar.sum() * 60 * 5
Qstar
="black", lw=2)
ax2.plot(qstar, color
ax2.fill_between(qstar.index, qstar,=0.0,
y2="tab:blue", alpha=0.3)
color
set(xlim=[df_q.index[0],
ax1.-1]],
df_q.index[=r"discharge (m$^3$/s)",
ylabel=[0, 5.5],
ylim=[0,1,2,3,4],
yticks="total discharge, q(t)")
titleset(yticks=[],
ax2.=[0, 5.5],
ylim=[df_q.index[0],
xlim-1]],
df_q.index[="effective discharge, q*(t)"
title
)
# makes slated dates plt.gcf().autofmt_xdate()
We can calculate \(p^*\) now, using
\[ P^* = Q^* \]
One of the simplest methods is to multiply \(p(t)\) by a fixed constant (<1) to obtain \(p^*\), so that the equation above holds true.
= Qstar/ P
ratio = df_p['PRECIPITATION'] * ratio
pstar = pstar.sum() * 5 * 60
Pstar print(f"Qstar / P = {ratio:.2f}")
Qstar / P = 0.16
Calculate now the centroid (\(t_pc\)) for effective precipitation \(p^*\) and centroid (\(t_{qc}\)) of effective discharge \(q^*\). Calculate also the time of peak discharge (\(t_{pk}\)). Then, calculate the centroid lag (\(T_{LC}\)), the centroid lag-to-peak (\(T_{LPC}\)), and the time of concentration (\(T_c\)). Use the equations below:
\(T_{LPC} \simeq 0.60 \cdot T_c\)
Time of precipitation centroid:
\[ t_{pc} = \frac{\displaystyle \sum_{i=1}^n p_i^* \cdot t_i}{P^*} \]
Time of streamflow centroid:
\[ t_{qc} = \frac{\displaystyle \sum_{i=1}^n q_i^* \cdot t_i}{Q^*} \]
Centroid lag:
\[ T_{LC} = t_{qc} - t_{pc} \]
Centroid lag-to-peak: \[ T_{LPC} = t_{pk} - t_{pc} \]
Time of concentration: \[ T_{LPC} \simeq 0.60 \cdot T_c \]
# pstar centroid
# time of the first (nonzero) rainfall data point
= pstar[pstar != 0.0].index[0]
t0 # time of the last (nonzero) rainfall data point
= pstar[pstar != 0.0].index[-1]
tf # duration of the rainfall event, in minutes
= (tf-t0) / pd.Timedelta('1 min')
td # make time array, add 2.5 minutes (half of dt)
= np.arange(0, td+1, 5) + 2.5
time # create pi array, only with relevant data (during rainfall duration)
= pstar.loc[(pstar.index >= t0) & (pstar.index <= tf)]
pi # convert from m3/5min to m3/s
= pi.values * 60 * 5
pi # time of precipitation centroid
= (pi * time).sum() / pi.sum()
t_pc # add initial time
= t0 + pd.Timedelta(minutes=t_pc)
t_pc
t_pc
# qstar centroid
# time of the first (nonzero) discharge data point
= qstar[qstar != 0.0].index[0]
t0 # time of the last (nonzero) discharge data point
= qstar[pstar != 0.0].index[-1]
tf # duration of the discharge event, in minutes
= (tf-t0) / pd.Timedelta('1 min')
td # make time array, add 2.5 minutes (half of dt)
= np.arange(0, td+1, 5) + 2.5
time # create qi array, only with relevant data (during discharge duration)
= qstar.loc[(qstar.index >= t0) & (qstar.index <= tf)]
qi # convert from m3/5min to m3/s
= qi.values * 60 * 5
qi # time of discharge centroid
= (qi * time).sum() / qi.sum()
t_qc # add initial time
= t0 + pd.Timedelta(minutes=t_qc)
t_qc
t_qc
# time of peak discharge
= qstar.max()
max_discharge = qstar[qstar == max_discharge].index[0]
t_pk
# centroid lag
= t_qc - t_pc
T_LC
# centroid lag-to-peak
= t_pk - t_pc
T_LPC
# time of concentration
= T_LPC / 0.60
T_c
print(f"T_LC = {T_LC}")
print(f"T_LPC = {T_LPC}")
print(f"T_c = {T_c}")
T_LC = 0 days 00:53:03.186594
T_LPC = 0 days 01:22:59.857820
T_c = 0 days 02:18:19.763033333