import csv
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Read in data from 124 model runs that have already been compiled into a .csv file. Read into arrays for the disturbance-rate parameter, the weathering-rate parameter, and the resulting facet slope angle.
filename = 'data/slope_analysis20180826.csv'
# Count number of lines in file
num_lines = len(open(filename).readlines( ))
# Create data arrays
dist_param = np.zeros(num_lines - 1) # skip 1 header line
weath_param = np.zeros(num_lines - 1)
facet_angle = np.zeros(num_lines - 1)
# Read data
with open(filename, 'r') as csvfile:
myreader = csv.reader(csvfile)
i = 0
for row in myreader:
print(', '.join(row))
if i >= 1:
dist_param[i-1] = row[1]
weath_param[i-1] = row[2]
facet_angle[i-1] = row[3]
i += 1
Convert the data into a pandas DataFrame. This allows us to sort the data set according to $d'$ and $w'$, which in turn will make it possible to turn the data into 4 x 31 array in which each row represents one set of experiments at a fixed value of $d'$.
import pandas as pd
data = {'dprime' : dist_param, 'wprime': weath_param, 'angle': facet_angle}
df = pd.DataFrame(data)
df = df.sort_values(by=['dprime', 'wprime'])
Reshape into 4 x 31 arrays in which each row is a unique value of $d'$. Print the dist_param values to make sure it worked.
facet_angle = df['angle'].values.reshape((4, 31))
weath_param = df['wprime'].values.reshape((4, 31))
dist_param = df['dprime'].values.reshape((4, 31))
dist_param
Read facet angle from dissolution-only runs. This is another set of model runs that we want to include on the figure.
# Name of file
filename = 'data/dissolution_only_results20180717.csv'
# Count number of lines in file
num_lines = len(open(filename).readlines( ))
# Create data arrays
diss_param = np.zeros(num_lines - 2) # we'll skip 2 header lines
diss_facet_angle = np.zeros(num_lines - 2)
# Read data
with open(filename, 'r') as csvfile:
myreader = csv.reader(csvfile)
i = 0
for row in myreader:
print(', '.join(row))
if i >= 2:
diss_param[i-2] = row[0]
diss_facet_angle[i-2] = row[2]
i += 1
Plot the data. First, we'll define four plotting symbols for each of the four values of $d'$. Then we'll plot, for each of the four, the facet angle as a function of weathering parameter ($w'$).
Next, we add the dissolution runs (using the dissolution parameter, which corresponds to a weathering parameter in a world where $d'\rightarrow \infty$).
Then add the curve that represents the analytical solution for the dissolution-only case. Also plot the model's effective 30$^\circ$ angle of repose for reference.
Finally, add axis labels, legend, etc., and save the figure to a file.
psyms = ['k.', 'k+', 'k^', 'k*']
# Plot the weathering + disturbance runs
for d in range(4):
plt.semilogx(10.0**3 * weath_param[d,:], facet_angle[d,:], psyms[d])
# Plot the dissolution runs
plt.semilogx(10.0**3 * diss_param, diss_facet_angle, 'ko')
# Analytical solution for dissolution only
wprime = 10**3 * 10.0 ** np.arange(-4.0, -0.9, 0.1)
ang_pred = 60.0 - 2 * (180.0 / np.pi) * wprime
plt.plot(wprime, ang_pred, 'k')
plt.ylim([0, 60])
# Angle of repose
ww = np.array([0.1, 100.0])
angrep = np.array([30.0, 30.0])
plt.plot(ww, angrep, 'g--')
# Axis limits
plt.xlim([0.1, 100])
# Labels and legend
plt.xlabel(r"Dimensionless weathering rate parameter, $w'$", fontsize=14)
plt.ylabel('Facet angle (degrees)', fontsize=14)
plt.legend([r"$d'= 10^{-1}$", r"$d'= 10^0$", r"$d'= 10^1$", r"$d'= 10^2$", r"$d' \rightarrow \infty$", r"$\theta = 60^\circ - 360 w' / \pi$", r"$\theta = 30^\circ$"], ncol=2, fontsize=10)
plt.savefig('facet_angle_vs_wprime.pdf')
This next figure isn't meant for the paper, but is just here to check out what it looks like if you plot difference between fault dip and facet angle instead of the facet angle itself.
dang = 60.0 - facet_angle
plt.semilogx(weath_param, dang, 'k.')
plt.semilogx(diss_param, 60.0 - diss_facet_angle, 'ko')
w = 10.0 ** np.arange(-4.0, -0.9, 0.1)
ang_pred = w * 10**3 * 360.0 / np.pi
plt.plot(w, ang_pred, 'k')
plt.ylim([0, 60])
