import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
set(style="ticks", font_scale=1.5)
sns.from scipy.optimize import curve_fit
import matplotlib.patches as patches
12 Exercises
12.1 Tasks
- Google the following: web plot digitizer
- Load image “nassif-16percent-slope.png” (see below)
- Create four csv files, one for each data set. Call them whatever you want. Legend: white circle = 312 mm/h, triangle = 234 mm/h, x = 156 mm/h, black circle = 78 mm/h.
The image is the second panel of Fig. 8, from
Nassif, S. H., and E. M. Wilson, 1975, “THE INFLUENCE OF SLOPE AND RAIN INTENSITY ON RUNOFF AND INFILTRATION”, Hydrological Sciences Journal. download here
Import relevant packages
Load all four files you created. Use numpy’s function loadtxt
. Make sure that the first point in each table corresponds to the appropriate rainfall rate. You can normalize the data if it is not.
= np.loadtxt("input_rate_078mm_per_h_16percent_slope.csv", delimiter=',')
d1 = np.loadtxt("input_rate_156mm_per_h_16percent_slope.csv", delimiter=',')
d2 = np.loadtxt("input_rate_234mm_per_h_16percent_slope.csv", delimiter=',')
d3 = np.loadtxt("input_rate_312mm_per_h_16percent_slope.csv", delimiter=',')
d4 1] = 78 * d1[:,1] / d1[:,1].max()
d1[:,1] = 156 * d2[:,1] / d2[:,1].max()
d2[:,1] = 234 * d3[:,1] / d3[:,1].max()
d3[:,1] = 312 * d4[:,1] / d4[:,1].max() d4[:,
Reproduce the original figure, make it look good, something like this:
= plt.subplots(figsize=(10,7))
fig, ax 0], d4[:,1], 'o', markerfacecolor="None", label=r"water input = 312 mm h$^{-1}$")
ax.plot(d4[:,0], d3[:,1], '^', label=r"water input = 234 mm h$^{-1}$")
ax.plot(d3[:,0], d2[:,1], 'x', label=r"water input = 156 mm h$^{-1}$")
ax.plot(d2[:,0], d1[:,1], 'o', label=r"water input = 78 mm h$^{-1}$")
ax.plot(d1[:,set(xlabel="Time (min)",
ax.=r"Infiltration rate (mm h$^{-1}$)")
ylabel="upper right"); ax.legend(loc
12.2 Horton’s equation
\[ f = f_c+(f_0-f_c)e^{-\beta t} \]
- \(f\): infiltration rate
- \(f_c\): infiltration capacity at large \(t\)
- \(f_0\): initial infiltration capacity
- \(\beta\): best fit empirical parameter
Write a function called horton
, that receives time t
and the three parameters, and returns the right-hand side of the equation above. Plot one of the data sets, together with a guess of the parameters that should roughly fit the data.
def horton(t, fc, f0, beta):
return fc + (f0 - fc)*np.exp(-beta*t)
= plt.subplots(figsize=(10,7))
fig, ax = d1[:,0]
t = t - t[0]
t = d1[:,1]
f 'o', label="data")
ax.plot(t, f, 35, 80, 0.5), '-', label="horton")
ax.plot(t, horton(t, set(xlabel="time (min)",
ax.="infiltration rate (mm/h)")
ylabel="upper right"); ax.legend(loc
Find the best fit for the parameters \(f_c, f_0, \beta\). Calculate the \(R^2\) for each data set.
For the best fit, use scipy’s curve_fit
. Write a function to compute the R-squared of your fit.
def horton(t, fc, f0, beta):
return fc + (f0 - fc)*np.exp(-beta*t)
def best_fit(data):
= data[:,0]
t = t[0]
t0 = t - t0
t = data[:,1]
f # best fit
= curve_fit(f=horton, # model function
popt, pcov =t, # x data
xdata=f, # y data
ydata=(130, 800, 0.5), # initial guess of the parameters
p0
)return [popt, pcov]
def calculate_r_squared(data, popt):
= data[:,0]
t = t - t[0]
t = data[:,1]
f # Calculate residuals
= f - horton(t, *popt)
residuals # You can get the residual sum of squares (ss_res) with
= np.sum(residuals**2)
ss_res # You can get the total sum of squares (ss_tot) with
= np.sum((f - np.mean(f))**2)
ss_tot # And finally, the r_squared-value with,
= 1 - (ss_res / ss_tot)
r_squared return r_squared
def plot_best_fit(data, axis, marker, markercolor):
# calculate best fit parameters
= best_fit(data)
popt, pcov = data[:,0]
t = data[:,1]
f # plot data points
=markercolor, markeredgecolor="black")
ax.plot(t, f, marker, markerfacecolor# plot best fit line
= calculate_r_squared(data, popt)
r_squared = r"$f_c=$ {:.2f}, $f_0=$ {:.2f}, $\beta=$ {:.2f}, $R^2=$ {:.2f}".format(popt[0],popt[1],popt[2], r_squared)
labeltext -t[0], *popt), color=markercolor, label=labeltext)
ax.plot(t, horton(t
= plt.subplots(figsize=(10,7))
fig, ax 'o', "tab:red")
plot_best_fit(d1, ax, 'x', "tab:blue")
plot_best_fit(d2, ax, '^', "tab:orange")
plot_best_fit(d3, ax, 'd', "tab:green")
plot_best_fit(d4, ax, set(xlabel="time (min)",
ax.="infiltration rate (mm/h)")
ylabel; ax.legend()
Make a graph of the infiltration rate and of the runoff, as a function of time. Use any of the four data sets you have.
= plt.subplots(figsize=(10,7))
fig, ax = d4
data = data[:, 0]
t = data[:, 1]
f = np.concatenate([ [0], t])
t = np.concatenate([ [f[0]], f])
f = f[0] - f
runoff *0 + f[0], ls="--", color="black", label="rainfall")
ax.plot(t, f="tab:blue", lw=3, label=r"infiltration")
ax.plot(t, f, color="tab:orange", lw=3, label=r"runoff")
ax.plot(t, runoff, colorset(xlabel="Time (min)",
ax.=r"Rate (mm h$^{-1}$)")
ylabel="lower right"); ax.legend(loc
12.3 Green & Ampt
\[f = \frac{A}{F} + B\]
where * \(A = K_\text{sat}\cdot|\psi_f|\cdot \left( \phi - \theta_0 \right)\) * \(B= K_\text{sat}\)
Write a function that calculates the cumulative of the infiltration rate.
\[ F(t) = \int_0^t f(t) \text{ d}t \]
Use numpy’s trapz
function, that implements the “trapezoidal rule”
def cumulative_F(t, f):
= np.array([0])
F = t/60 # convert minute to hour
t for i in np.arange(2,len(t)+1):
= np.trapz(f[:i], t[:i])
area = np.concatenate([F, [area]])
F return F
= plt.subplots(1, 2, figsize=(10,7))
fig, (ax1, ax2) = d1[:,0], d1[:,1]
t, f = cumulative_F(t, f)
F ="f, rate")
ax1.plot(t, f, label="F, cumulative")
ax2.plot(t, F, labelset(xlabel="t (min)",
ax1.="f (mm/h)")
ylabelset(xlabel="t (min)",
ax2.="F (mm)")
ylabel"right") ax2.yaxis.set_label_position(
Plot \(f\) as a function of \(F\). Try to guess \(A\) and \(B\) that give reasonable results.
= plt.subplots(figsize=(10,7))
fig, ax = d1[:,0], d1[:,1]
t, f = cumulative_F(t, f)
F
ax.plot(F, f)=50; B=30;
A/F + B, 'o')
ax.plot(F, Aset(xlabel="F",
ax.="f") ylabel
/Users/yairmau/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: divide by zero encountered in true_divide
[Text(0.5, 0, 'F'), Text(0, 0.5, 'f')]
Use the curve_fit
to find the optimal values for \(A\) and \(B\).
def G_and_A(F, A, B):
return A/F + B
= curve_fit(f=G_and_A, # model function
popt, pcov =F[1:], # x data
xdata=f[1:], # y data
ydata=(50, 30), # initial guess of the parameters
p0
)
# popt, pcov = curve_fit(G_and_A, F[1:], f[1:], p0=(50, 30)) # p0 = initial guess
print(popt)
= plt.subplots(figsize=(10,7))
fig, ax
ax.plot(F, f)1:], popt[0]/F[1:] + popt[1], 'o')
ax.plot(F[set(xlabel="F",
ax.="f") ylabel
[24.12368526 36.34242813]
[Text(0.5, 0, 'F'), Text(0, 0.5, 'f')]
12.4 Homework
Go to Soil Texture Calculator, estimate the texture of “standard soil” in Nassif & Wilson, 1975.