Outwasher pull request - tools folder update

.gitignore

@@ -138,3 +138,9 @@ brie_venv
 Output/
 Run_Output/
 /brie_venv/
+/data/outwash_data/storms/
+/data/outwash_data/OLD/
+/notebooks/for_thesis
+/notebooks/old_outwash_runs
+/notebooks/post-thesis
+scripts/outwash_ms/OLD

cascade/beach_dune_manager.py

@@ -662,7 +662,6 @@ def update(
         if self._overwash_removal:
             # barrier3d saves the pre-storm interior for each time step
             pre_storm_interior = barrier3d.PreStorm_InteriorDomain
-
             pre_storm_interior, post_storm_interior = resize_interior_domain(
                 pre_storm_interior=pre_storm_interior,
                 post_storm_interior=barrier3d.DomainTS[self._time_index - 1],

cascade/tools/creating_input_domains.py

@@ -0,0 +1,177 @@
+# Lexi Van Blunk
+# last updated 8/3/2023
+
+# this code breaks up the "observed" configurations into interior, dune, and beach domains for input into CASCADE
+
+# NOTE: configurations 1 and 2 vary slightly here because I was still using the domains from the thesis rather than
+# the ones from the new larger DEM. In the future, I would suggest using these newer domains since they are consistent
+# with configurations 3 and 4
+
+import numpy as np
+
+berm_el = 0.11
+save_dir = "C:/Users/Lexi/PycharmProjects/CASCADE/scripts/outwash_ms/configurations/"
+save_now = False
+
+# ------------------------------------------------ load data -----------------------------------------------------------
+
+section1 = np.load(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config1\config1_observed_pre.npy"
+)
+section2 = np.load(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config2\config2_observed_pre.npy"
+)
+section3 = np.load(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config3\config3_observed_pre.npy"
+)
+section4 = np.load(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config4\config4_observed_pre.npy"
+)
+
+
+# --------------------------------------------- section 1 --------------------------------------------------------------
+section1_int = section1[0:30]
+section1_dunes = section1[30:35]
+section1_beach = section1[35:]
+
+# take only two rows for the dunes because B3D only allows two dune rows
+section1_dunes = section1_dunes[1:3]
+
+full_1 = np.append(section1_int, section1_dunes, 0)
+full_1 = np.append(full_1, section1_beach, 0)
+
+if save_now:
+    # save the domains as needed for B3D
+    save_dir_1 = save_dir + "config1/"
+
+    # flip the interior for correct orientation in B3D
+    interior_b3d_input_1 = np.flip(section1_int)
+    np.save(
+        save_dir_1 + "NCB-default-elevation-config1-damMHW.npy", interior_b3d_input_1
+    )
+
+    # flip the dunes and put them into one long array
+    dunes_b3d_1 = np.flip(section1_dunes) - berm_el
+    dunes_input_1 = np.append(dunes_b3d_1[0], dunes_b3d_1[1], 0)
+    np.save(save_dir_1 + "NCB-default-dunes-config1-dam.npy", dunes_input_1)
+
+    # since the beach is not used in B3D, it can be saved normally (outwash orientation)
+    np.save(save_dir_1 + "NCB-default-beach-config1-damMHW.npy", section1_beach)
+
+    # save the full input just because
+    np.save(save_dir_1 + "NCB-default_full_config1-damMHW.npy", full_1)
+
+
+# --------------------------------------------- section 2 --------------------------------------------------------------
+section2_int = section2[0:27]
+section2_dunes = section2[27:31]
+section2_beach = section2[31:]
+
+# manipulate the dunes into 2 rows
+section2_dunes[0:2, 0:4] += 0.1
+section2_dunes[0:2, 0:2] = section2_dunes[3, 0:2]
+section2_dunes[0, 2:5] = section2_dunes[2, 2:5]
+section2_dunes = section2_dunes[0:2]
+section2_dunes[0, 0:20] += 0.05
+section2_dunes[1, 0:16] += 0.1
+
+full_2 = np.append(section2_int, section2_dunes, 0)
+full_2 = np.append(full_2, section2_beach, 0)
+
+if save_now:
+    # save the domains as needed for B3D
+    save_dir_2 = save_dir + "config2/"
+
+    # flip the interior for correct orientation in B3D
+    interior_b3d_input_2 = np.flip(section2_int)
+    np.save(
+        save_dir_2 + "NCB-default-elevation-config2-damMHW.npy", interior_b3d_input_2
+    )
+
+    # flip the dunes and put them into one long array
+    dunes_b3d_2 = np.flip(section2_dunes) - berm_el
+    dunes_input_2 = np.append(dunes_b3d_2[0], dunes_b3d_2[1], 0)
+    np.save(save_dir_2 + "NCB-default-dunes-config2-dam.npy", dunes_input_2)
+
+    # since the beach is not used in B3D, it can be saved normally (outwash orientation)
+    np.save(save_dir_2 + "NCB-default-beach-config2-damMHW.npy", section2_beach)
+
+    # save the full input just because
+    np.save(save_dir_2 + "NCB-default_full_config2-damMHW.npy", full_2)
+
+# --------------------------------------------- section 3 --------------------------------------------------------------
+
+section3_int = section3[:23, :]
+section3_dunes = section3[23:28, :]
+section3_beach = section3[28:, :]
+
+# manipulate the dunes into 2 rows that represent the true dune line
+section3_stitched_dunes = np.zeros([2, np.shape(section3_dunes)[1]])
+section3_stitched_dunes[:, 0:4] = section3[26:28, 0:4]
+section3_stitched_dunes[:, 4:9] = section3[25:27, 4:9]
+section3_stitched_dunes[:, 9:22] = section3[24:26, 9:22]
+section3_stitched_dunes[:, 22:] = section3[23:25, 22:]
+
+# # filling in the beach under the dunes
+extra_beach = np.zeros([3, np.shape(section3_dunes)[1]])
+extra_beach[:, 22:] = section3[25:28, 22:]
+extra_beach[:, 0:4] = section3[28:31, 0:4]
+extra_beach[:, 4:9] = section3[27:30, 4:9]
+extra_beach[:, 9:22] = section3[26:29, 9:22]
+section3_beach[0:2, 0:9] = section3[30:32, 0:9]
+section3_beach = np.append(extra_beach, section3_beach, 0)
+
+full_3 = np.append(section3_int, section3_stitched_dunes, 0)
+full_3 = np.append(full_3, section3_beach, 0)
+
+if save_now:
+    # save the domains as needed for B3D
+    save_dir_3 = save_dir + "config3/"
+
+    # flip the interior for correct orientation in B3D
+    interior_b3d_input_3 = np.flip(section3_int)
+    np.save(
+        save_dir_3 + "NCB-default-elevation-config3-damMHW.npy", interior_b3d_input_3
+    )
+
+    # flip the dunes and put them into one long array
+    dunes_b3d_3 = np.flip(section3_dunes) - berm_el
+    dunes_input_3 = np.append(dunes_b3d_3[0], dunes_b3d_3[1], 0)
+    np.save(save_dir_3 + "NCB-default-dunes-config3-dam.npy", dunes_input_3)
+
+    # since the beach is not used in B3D, it can be saved normally (outwash orientation)
+    np.save(save_dir_3 + "NCB-default-beach-config3-damMHW.npy", section3_beach)
+
+    # save the full input just because
+    np.save(save_dir_3 + "NCB-default_full_config3-damMHW.npy", full_3)
+
+# --------------------------------------------- section 4 --------------------------------------------------------------
+
+section4_int = section4[:26]
+section4_dunes = section4[28:30]
+section4_beach = section4[31:]
+full_4 = np.append(section4_int, section4_dunes, 0)
+full_4 = np.append(full_4, section4_beach, 0)
+
+beach4_slope = (np.mean(section4_beach[0, :]) - np.mean(section4_beach[8, :])) / 9
+
+if save_now:
+    # save the domains as needed for B3D
+    save_dir_4 = save_dir + "config4/"
+
+    # flip the interior for correct orientation in B3D
+    interior_b3d_input_4 = np.flip(section4_int)
+    np.save(
+        save_dir_4 + "NCB-default-elevation-config4-damMHW.npy", interior_b3d_input_4
+    )
+
+    # flip the dunes and put them into one long array
+    dunes_b3d_4 = np.flip(section4_dunes) - berm_el
+    dunes_input_4 = np.append(dunes_b3d_4[0], dunes_b3d_4[1], 0)
+    np.save(save_dir_4 + "NCB-default-dunes-config4-dam.npy", dunes_input_4)
+
+    # since the beach is not used in B3D, it can be saved normally (outwash orientation)
+    np.save(save_dir_4 + "NCB-default-beach-config4-damMHW.npy", section4_beach)
+
+    # save the full input just because
+    np.save(save_dir_4 + "NCB-default_full_config4-damMHW.npy", full_4)

cascade/tools/lidar_to_array.py

@@ -0,0 +1,199 @@
+# Lexi Van Blunk
+# last updated 8/3/2023
+
+# this code loads the csv file of elevation points from ArcGIS and organizes them into python arrays
+# we also break up the domains into the 4 "observed" configurations and save them to the
+# scripts/outwash_ms/configurations folder
+
+# SEE NOTEBOOK: plotting observed configurations, which does the same thing but visualizes each configuration into
+# subplots displaying the prestorm, poststorm, and elevation change plots (and shows boxes where we calculate erosion)
+# we also calculate erosion in the notebook
+
+import os
+
+import numpy as np
+import pandas as pd
+from matplotlib import pyplot as plt
+
+plt.rcParams.update({"font.size": 15})
+
+
+# --------------------------------------- load the data as a dataframe ------------------------------------------------
+
+os.chdir("C:/Users/Lexi/Documents/Research/Outwasher Paper/GIS data/")
+
+# PRESTORM DATA
+prestorm_points_df = pd.read_csv("raster_points.csv")
+
+# get the column from the csv that contains the elevations
+prestorm_elevations = prestorm_points_df["grid_code"]
+
+alongshore_length_feature_class = (
+    4000  # in meters, along-shore side of rectangle from GIS (barrier length)
+)
+num_cols_domain = int(
+    alongshore_length_feature_class / 10
+)  # we specified 10x10 m cells
+num_rows_domain = int(
+    len(prestorm_elevations) / num_cols_domain
+)  # this is the width of the island
+MHW = 0.36  # meters NAVD88? (subtract 0.36 m from NAVD88 elevations to get MHW elevations)
+
+prestorm_elevations = prestorm_elevations - MHW  # convert elevations from NAVD88 to MHW
+
+# POSTSTORM DATA
+post_points_df = pd.read_csv("post_storm_points.csv")
+post_elevations = post_points_df["grid_code"]
+post_elevations = post_elevations - MHW  # convert elevations from NAVD88 to MHW
+
+
+# ---------------------------- organize the list of elevations into an array--------------------------------------------
+
+elev_array = np.zeros([num_rows_domain, num_cols_domain])
+start = 0
+end = start + num_cols_domain
+for row in range(num_rows_domain):
+    elev_array[row, :] = prestorm_elevations[start:end]
+    start += num_cols_domain
+    end += num_cols_domain
+
+post_elev_array = np.zeros([num_rows_domain, num_cols_domain])
+start = 0
+end = start + num_cols_domain
+for row in range(num_rows_domain):
+    post_elev_array[row, :] = post_elevations[start:end]
+    start += num_cols_domain
+    end += num_cols_domain
+
+# set all values less than 0 to -3 m
+for r in range(150):
+    for c in range(num_cols_domain):
+        if post_elev_array[r, c] <= -0.5:
+            post_elev_array[r, c] = -3
+
+
+# ------------------------------------------ plotting the domains ------------------------------------------------------
+xlabel = "alongshore distance (m)"
+ylabel = "cross-shore distance (m)"
+
+# PRESTORM
+fig1 = plt.figure()
+ax1 = fig1.add_subplot(111)
+mat = ax1.matshow(
+    elev_array,
+    cmap="terrain",
+    vmin=-3.0,
+    vmax=6.0,
+)
+cbar = fig1.colorbar(mat)
+cbar.set_label("m MHW", rotation=270, labelpad=10)
+ax1.set_title("Pre-Storm Elevation")
+ax1.set_ylabel(ylabel)
+ax1.set_xlabel(xlabel)
+plt.gca().xaxis.tick_bottom()
+
+# add the rectangles around the 4 configurations
+# (x,y), length, width (up = neg)
+ax1.add_patch(plt.Rectangle((98, 187), 50, -40, lw=2, ec="k", fc="none"))  # black
+ax1.add_patch(plt.Rectangle((148, 187), 50, -40, lw=2, ec="m", fc="none"))  # purple
+ax1.add_patch(
+    plt.Rectangle((198, 187), 50, -40, lw=2, ec="tab:orange", fc="none")
+)  # blue
+ax1.add_patch(
+    plt.Rectangle((248, 187), 50, -40, lw=2, ec="tab:blue", fc="none")
+)  # orange
+
+# edit the tick marks to be in meters
+xtick_max = np.shape(elev_array)[1]  # n_cols = x
+x_ticks = np.array(range(0, xtick_max, 50))
+x_tick_labels = x_ticks * 10
+ytick_max = np.shape(elev_array)[0]  # n_rows = y
+y_ticks = np.array(range(0, ytick_max, 25))
+y_tick_labels = y_ticks * 10
+plt.xticks(x_ticks, x_tick_labels)
+plt.yticks(y_ticks, y_tick_labels)
+plt.tight_layout()
+
+# POSTSTORM
+fig2 = plt.figure()
+ax2 = fig2.add_subplot(111)
+mat = ax2.matshow(
+    post_elev_array,
+    cmap="terrain",
+    vmin=-3.0,
+    vmax=6.0,
+)
+cbar = fig2.colorbar(mat)
+cbar.set_label("m MHW", rotation=270, labelpad=10)
+ax2.set_title("Post-Storm Elevation")
+ax2.set_ylabel(ylabel)
+ax2.set_xlabel(xlabel)
+plt.gca().xaxis.tick_bottom()
+
+# add the rectangles around the 4 configurations
+ax2.add_patch(plt.Rectangle((98, 187), 50, -40, lw=2, ec="k", fc="none"))  # black
+ax2.add_patch(plt.Rectangle((148, 187), 50, -40, lw=2, ec="m", fc="none"))  # purple
+ax2.add_patch(
+    plt.Rectangle((198, 187), 50, -40, lw=2, ec="tab:orange", fc="none")
+)  # blue
+ax2.add_patch(
+    plt.Rectangle((248, 187), 50, -40, lw=2, ec="tab:blue", fc="none")
+)  # orange
+
+# edit the tick marks to be in meters
+xtick_max = np.shape(post_elev_array)[1]  # n_cols = x
+x_ticks = np.array(range(0, xtick_max, 50))
+x_tick_labels = x_ticks * 10
+ytick_max = np.shape(post_elev_array)[0]  # n_rows = y
+y_ticks = np.array(range(0, ytick_max, 25))
+y_tick_labels = y_ticks * 10
+plt.xticks(x_ticks, x_tick_labels)
+plt.yticks(y_ticks, y_tick_labels)
+plt.tight_layout()
+
+# --------------------------- convert to dam and make the 4 configurations ---------------------------------------------
+elev_array = elev_array / 10
+section1_pre = elev_array[147:188, 98:148]
+section2_pre = elev_array[147:188, 148:198]
+section3_pre = elev_array[147:188, 248:298]
+section4_pre = elev_array[147:188, 198:248]  # validation configuration
+
+np.save(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config1\config1_observed_pre.npy",
+    section1_pre,
+)
+np.save(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config2\config2_observed_pre.npy",
+    section2_pre,
+)
+np.save(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config3\config3_observed_pre.npy",
+    section3_pre,
+)
+np.save(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config4\config4_observed_pre.npy",
+    section4_pre,
+)
+
+post_elev_array = post_elev_array / 10
+section1_post = post_elev_array[147:188, 98:148]
+section2_post = post_elev_array[147:188, 148:198]
+section3_post = post_elev_array[147:188, 248:298]
+section4_post = post_elev_array[147:188, 198:248]
+
+np.save(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config1_observed_post.npy",
+    section1_post,
+)
+np.save(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config2_observed_post.npy",
+    section2_post,
+)
+np.save(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config3_observed_post.npy",
+    section3_post,
+)
+np.save(
+    r"C:\Users\Lexi\PycharmProjects\CASCADE\scripts\outwash_ms\configurations\config4_observed_post.npy",
+    section4_post,
+)