Skip to content

Commit 0770529

Browse files
author
Lexi Van Blunk
committed
this code takes the average of the entire domain for each time step rather than filtering out the cells below sea level
1 parent 7e122db commit 0770529

File tree

1 file changed

+221
-0
lines changed

1 file changed

+221
-0
lines changed

‎scripts/outwash_ms/average_interior_elevation.py‎

Lines changed: 221 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,221 @@
1+
# loading in the npz files and saving variables that we want to plot together
2+
3+
# Lexi Van Blunk
4+
# 5/15/2025
5+
# calculating the average interior elevation using DomainTS instead of the h_b_TS variable because I think it is
6+
# messed up
7+
8+
import numpy as np
9+
from matplotlib import pyplot as plt
10+
11+
# ---------------------------------- set model parameters that change per run ------------------------------------------
12+
rname_array = ["r025", "r035"]
13+
# rname_array = ["r035"]
14+
for rname in rname_array:
15+
storm_interval = 20 # 20 or 10 years
16+
config = 4 # 1, 2, 3, or 4
17+
18+
# Display stats on console/show plots
19+
migration_stats = False
20+
plotters = False
21+
geomoetry_stats = True
22+
dune_crest_stats = False
23+
flux_stats = False
24+
25+
# location of the npz files
26+
datadir_b3d = "C:/Users/Lexi/PycharmProjects/CASCADE/data/outwash_data/storms/slope0pt03/rerun_output/{0}/overwash_only/".format(rname)
27+
datadir_100 = "C:/Users/Lexi/PycharmProjects/CASCADE/data/outwash_data/storms/slope0pt03/rerun_output/{0}/outwash100/".format(rname)
28+
datadir_50 = "C:/Users/Lexi/PycharmProjects/CASCADE/data/outwash_data/storms/slope0pt03/rerun_output/{0}/outwash50/".format(rname)
29+
datadir_0 = "C:/Users/Lexi/PycharmProjects/CASCADE/data/outwash_data/storms/slope0pt03/rerun_output/{0}/outwash0/".format(rname)
30+
31+
# initialize empty arrays
32+
tmax_array_b3d = []
33+
tmax_array_100 = []
34+
tmax_array_50 = []
35+
tmax_array_0 = []
36+
37+
tmax_first_drown_100 = []
38+
tmax_first_drown_50 = []
39+
tmax_first_drown_0 = []
40+
41+
avg_elev_TS_array_b3d = []
42+
avg_elev_TS_array_100 = []
43+
avg_elev_TS_array_50 = []
44+
avg_elev_TS_array_0 = []
45+
46+
avg_elev_b3d_array = np.zeros(100, dtype=object)
47+
avg_elev_100_array = np.zeros(100, dtype=object)
48+
avg_elev_50_array = np.zeros(100, dtype=object)
49+
avg_elev_0_array = np.zeros(100, dtype=object)
50+
51+
avg_last_elev_b3d_array = np.zeros(100, dtype=object)
52+
avg_last_elev_100_array = np.zeros(100, dtype=object)
53+
avg_last_elev_50_array = np.zeros(100, dtype=object)
54+
avg_last_elev_0_array = np.zeros(100, dtype=object)
55+
56+
for storm_num in range(1, 101):
57+
# print(storm_num)
58+
# b3d variables
59+
filename_b3d = "config{0}_b3d_startyr1_interval{1}yrs_Slope0pt03_{2}.npz".format(config, storm_interval, storm_num)
60+
file_b3d = datadir_b3d + filename_b3d
61+
b3d = np.load(file_b3d, allow_pickle=True)
62+
b3d_obj = b3d["cascade"][0]
63+
tmax_b3d = b3d_obj.barrier3d[0].TMAX
64+
tmax_array_b3d.append(tmax_b3d)
65+
66+
for t in range(0, tmax_b3d):
67+
# average elevation
68+
domain_TS = b3d_obj.barrier3d[0].DomainTS[t] * 10
69+
avg_elev_TS = np.average(domain_TS)
70+
avg_elev_TS_array_b3d.append(avg_elev_TS)
71+
72+
last_domain_TS_avg = np.average(domain_TS)
73+
avg_elev_b3d = np.average(avg_elev_TS_array_b3d)
74+
avg_elev_b3d_array[storm_num-1] = avg_elev_b3d
75+
avg_last_elev_b3d_array[storm_num-1] = last_domain_TS_avg
76+
77+
# 100% variables
78+
filename_100 = "config{0}_outwash100_startyr1_interval{1}yrs_Slope0pt03_{2}.npz".format(config, storm_interval, storm_num)
79+
file_100 = datadir_100 + filename_100
80+
outwash100 = np.load(file_100, allow_pickle=True)
81+
outwash100_obj = outwash100["cascade"][0]
82+
tmax_100 = outwash100_obj.barrier3d[0].TMAX
83+
tmax_array_100.append(tmax_100)
84+
85+
# average elevation
86+
if tmax_100 < 101:
87+
max_t = tmax_100 + 1
88+
else:
89+
max_t = tmax_100
90+
for t in range(0, max_t):
91+
# average elevation
92+
domain_TS = outwash100_obj.barrier3d[0].DomainTS[t] * 10
93+
avg_elev_TS = np.average(domain_TS)
94+
avg_elev_TS_array_100.append(avg_elev_TS)
95+
96+
last_domain_TS_avg = np.average(domain_TS) #domain TS will be the domain before drowning
97+
avg_elev_100 = np.average(avg_elev_TS_array_100)
98+
avg_elev_100_array[storm_num-1] = avg_elev_100
99+
avg_last_elev_100_array[storm_num-1] = last_domain_TS_avg
100+
101+
# 50% variables
102+
filename_50 = "config{0}_outwash50_startyr1_interval{1}yrs_Slope0pt03_{2}.npz".format(config, storm_interval, storm_num)
103+
file_50 = datadir_50 + filename_50
104+
outwash50 = np.load(file_50, allow_pickle=True)
105+
outwash50_obj = outwash50["cascade"][0]
106+
tmax_50 = outwash50_obj.barrier3d[0].TMAX
107+
tmax_array_50.append(tmax_50)
108+
109+
# average elevation
110+
if tmax_50 < 101:
111+
max_t = tmax_50 + 1
112+
else:
113+
max_t = tmax_50
114+
for t in range(0, max_t):
115+
# average elevation
116+
domain_TS = outwash50_obj.barrier3d[0].DomainTS[t] * 10
117+
avg_elev_TS = np.average(domain_TS)
118+
avg_elev_TS_array_50.append(avg_elev_TS)
119+
120+
last_domain_TS_avg = np.average(domain_TS)
121+
avg_elev_50 = np.average(avg_elev_TS_array_50)
122+
avg_elev_50_array[storm_num-1] = avg_elev_50
123+
avg_last_elev_50_array[storm_num-1] = last_domain_TS_avg
124+
125+
# 0% variables
126+
filename_0 = "config{0}_outwash0_startyr1_interval{1}yrs_Slope0pt03_{2}.npz".format(config, storm_interval, storm_num)
127+
file_0 = datadir_0 + filename_0
128+
outwash0 = np.load(file_0, allow_pickle=True)
129+
outwash0_obj = outwash0["cascade"][0]
130+
tmax_0 = outwash0_obj.barrier3d[0].TMAX
131+
tmax_array_0.append(tmax_0)
132+
133+
# average elevation
134+
if tmax_0 < 101:
135+
max_t = tmax_0 + 1
136+
else:
137+
max_t = tmax_0
138+
for t in range(0, max_t):
139+
# average elevation
140+
domain_TS = outwash0_obj.barrier3d[0].DomainTS[t] * 10
141+
avg_elev_TS = np.average(domain_TS)
142+
avg_elev_TS_array_0.append(avg_elev_TS)
143+
144+
# last_domain_TS_avg = np.average(outwash0_obj.barrier3d[0].DomainTS[-1]*10)
145+
last_domain_TS_avg = np.average(domain_TS)
146+
avg_elev_0 = np.average(avg_elev_TS_array_0)
147+
avg_elev_0_array[storm_num-1] = avg_elev_0
148+
avg_last_elev_0_array[storm_num-1] = last_domain_TS_avg
149+
150+
151+
# printing geometry stats
152+
if geomoetry_stats:
153+
# print avg height and width stats
154+
print("avg geometry stats, {0}".format(rname))
155+
print("baseline \n avg interior height: {0} \n avg last interior height: {1}"
156+
.format(np.round(np.average(avg_elev_b3d_array), 2),
157+
np.round(np.average(avg_last_elev_b3d_array), 2)))
158+
print("100% outwash \n avg interior height: {0} \n avg last interior height: {1}"
159+
.format(np.round(np.average(avg_elev_100_array), 2),
160+
np.round(np.average(avg_last_elev_100_array), 2)))
161+
print("50% outwash \n avg interior height: {0} \n avg last interior height: {1}"
162+
.format(np.round(np.average(avg_elev_50_array), 2),
163+
np.round(np.average(avg_last_elev_50_array), 2)))
164+
print("0% outwash \n avg interior height: {0} \n avg last interior height: {1}"
165+
.format(np.round(np.average(avg_elev_0_array), 2),
166+
np.round(np.average(avg_last_elev_0_array), 2)))
167+
168+
169+
# if plotters:
170+
# storm_num = 1
171+
# fig9 = plt.figure()
172+
# # fig9.suptitle('overwash storm {0} - {1}'.format(storm_num, rname), weight="bold")
173+
# # fig9.suptitle('{0}'.format(rname), weight="bold")
174+
# ax1 = fig9.add_subplot(111)
175+
# ls = "solid"
176+
# ax1.plot(shoreline_pos_array_b3d[storm_num-1])
177+
# ax1.plot(shoreline_pos_array_100[storm_num-1], linestyle=ls)
178+
# ax1.plot(shoreline_pos_array_50[storm_num-1], linestyle=ls)
179+
# ax1.plot(shoreline_pos_array_0[storm_num-1], linestyle=ls)
180+
# ax1.legend(["baseline", "100% outwash to shoreface", "50% outwash to shoreface", "0% outwash to shoreface"],
181+
# prop={'size': 9}, loc="upper left")
182+
# ax1.set_ylabel("Shoreline Position (m)")
183+
# ax1.set_xlabel("Simulation Years")
184+
# ax1.set_ylim(bottom=-60, top=160)
185+
# ax1.set_title('{0}'.format(rname), weight="bold")
186+
# fig9.subplots_adjust(hspace=0.3)
187+
#
188+
# # plotting one year of overwash and outwash flux
189+
# fig10 = plt.figure()
190+
# # fig9.suptitle('overwash storm {0} - {1}'.format(storm_num, rname), weight="bold")
191+
# # fig9.suptitle('{0}'.format(rname), weight="bold")
192+
# ax1 = fig10.add_subplot(211)
193+
# ls = "solid"
194+
# ax1.plot(overwash_array_b3d[storm_num-1])
195+
# ax1.plot(overwash_array_100[storm_num-1], linestyle=ls)
196+
# ax1.plot(overwash_array_50[storm_num-1], linestyle=ls)
197+
# ax1.plot(overwash_array_0[storm_num-1], linestyle=ls)
198+
# ax1.set_ylabel("Overwash Flux (m3/m)")
199+
# ax1.set_xlabel("Simulation Years")
200+
# ax1.set_ylim(top=150)
201+
# ax1.set_title(rname, weight="bold")
202+
# ax1.legend(["baseline", "100% outwash to shoreface", "50% outwash to shoreface", "0% outwash to shoreface"],
203+
# prop={'size': 9}, loc="upper left")
204+
#
205+
# max_t_b3d = tmax_array_b3d[storm_num-1]
206+
# max_t_100 = tmax_array_100[storm_num - 1]
207+
# max_t_50 = tmax_array_50[storm_num - 1]
208+
# max_t_0 = tmax_array_0[storm_num - 1]
209+
#
210+
# outwash_array_b3d_storm = outwash_array_b3d[storm_num-1]
211+
# outwash_array_b3d_storm[max_t_b3d+1:] = np.nan
212+
# outwash_array_100_storm = outwash_array_100[storm_num - 1]
213+
# outwash_array_100_storm[max_t_100+1:] = np.nan
214+
# outwash_array_50_storm = outwash_array_50[storm_num - 1]
215+
# outwash_array_50_storm[max_t_50+1:] = np.nan
216+
# outwash_array_0_storm = outwash_array_0[storm_num - 1]
217+
# outwash_array_0_storm[max_t_0+1:] = np.nan
218+
219+
220+
221+

0 commit comments

Comments
 (0)