Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

steepness plots #4

Open
SophieRothman opened this issue May 15, 2023 · 0 comments
Open

steepness plots #4

SophieRothman opened this issue May 15, 2023 · 0 comments

Comments

@SophieRothman
Copy link
Contributor

FIRST BOX OF CODE

Code Block 5

Original K_sp value is 1e-5

K_sp = 1.0e-5 # units vary depending on m_sp and n_sp
m_sp = 0.5 # exponent on drainage area in stream power equation
n_sp = 1.0 # exponent on slope in stream power equation

frr = FlowAccumulator(mg3, flow_director="FlowDirectorD8") # intializing flow routing
spr = StreamPowerEroder(
mg3, K_sp=K_sp, m_sp=m_sp, n_sp=n_sp, threshold_sp=0.0
) # initializing stream power incision

theta = m_sp / n_sp

initialize the component that will calculate channel steepness

sf = SteepnessFinder(mg3, reference_concavity=theta, min_drainage_area=10.0)

initialize the component that will calculate the chi index

cf = ChiFinder(
mg3, min_drainage_area=10.0, reference_concavity=theta, use_true_dx=True
)

SECOND BOX OF CODE

Code Block 11

calculate the chi index

cf.calculate_chi()

chi-elevation plots in the profiled channels

plt.figure(4)

for i, outlet_id in enumerate(prf.data_structure):
for j, segment_id in enumerate(prf.data_structure[outlet_id]):
if j == 0:
label = f"channel {i + 1}"
else:
label = "nolegend"
segment = prf.data_structure[outlet_id][segment_id]
profile_ids = segment["ids"]
color = segment["color"]
plt.plot(
mg3.at_node["channel__chi_index"][profile_ids],
mg3.at_node["topographic__elevation"][profile_ids],
color=color,
label=label,
)

plt.xlabel("chi index (m)")
plt.ylabel("elevation (m)")
plt.legend(loc="lower right")
#title_text = (

f"$K_{{sp}}$={K_sp}; $time$={total_time} yr; $dx$={dxy} m; concavity={theta}"

#)
#plt.title(title_text)

chi map

plt.figure(5)
imshow_grid(
mg3,
"channel__chi_index",
grid_units=("m", "m"),
var_name="Chi index (m)",
cmap="jet",
)
#title_text = (

f"$K_{{sp}}$={K_sp}; $time$={total_time} yr; $dx$={dxy} m; concavity={theta}"

#)
#plt.title(title_text)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant