Add Holman 1986 to list of runup functions

Makefile

@@ -110,6 +110,16 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte
 	--profile-type "prestorm" \
 	--output-file "./data/interim/twl_mean_slope_sto06.csv"
 
+./data/interim/twl_mean_slope_hol86.csv: ./data/interim/waves.csv ./data/interim/tides.csv ./data/interim/profiles.csv ./data/interim/sites.csv ./data/interim/profile_features_crest_toes.csv
+	$(PYTHON_CLI) create-twl-forecast \
+	--waves-csv "./data/interim/waves.csv" \
+	--tides-csv "./data/interim/tides.csv" \
+	--profiles-csv "./data/interim/profiles.csv" \
+	--profile-features-csv "./data/interim/profile_features_crest_toes.csv" \
+	--runup-function "hol86" \
+	--slope "mean" \
+	--profile-type "prestorm" \
+	--output-file "./data/interim/twl_mean_slope_hol86.csv"
 
 ### IMPACTS
 
@@ -132,10 +142,15 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte
 	--forecasted-twl-csv "./data/interim/twl_foreshore_slope_sto06.csv" \
 	--output-file "./data/interim/impacts_forecasted_foreshore_slope_sto06.csv"
 
+./data/interim/impacts_forecasted_mean_slope_hol86.csv: ./data/interim/profile_features_crest_toes.csv ./data/interim/twl_mean_slope_hol86.csv
+	$(PYTHON_CLI) create-forecasted-impacts \
+	--profile-features-csv "./data/interim/profile_features_crest_toes.csv" \
+	--forecasted-twl-csv "./data/interim/twl_mean_slope_hol86.csv" \
+	--output-file "./data/interim/impacts_forecasted_mean_slope_hol86.csv"
 
 ### GEOJSONs
 
-geojsons: ./data/interim/impacts_forecasted_mean_slope_sto06.geojson ./data/interim/impacts_forecasted_mean_slope_sto06_R_high.geojson ./data/interim/profile_features_crest_toes.geojson ./data/interim/sites.geojson
+geojsons: ./data/interim/impacts_forecasted_mean_slope_hol86.geojson ./data/interim/impacts_forecasted_mean_slope_hol86_R_high.geojson ./data/interim/impacts_forecasted_mean_slope_sto06.geojson ./data/interim/impacts_forecasted_mean_slope_sto06_R_high.geojson ./data/interim/profile_features_crest_toes.geojson ./data/interim/sites.geojson
 
 ./data/interim/impacts_forecasted_mean_slope_sto06.geojson: ./data/interim/impacts_forecasted_mean_slope_sto06.csv ./data/interim/impacts_observed.csv
 	$(PYTHON_CLI) impacts-to-geojson \
@@ -152,6 +167,21 @@ geojsons: ./data/interim/impacts_forecasted_mean_slope_sto06.geojson ./data/inte
 	--impacts-csv "./data/interim/impacts_forecasted_mean_slope_sto06.csv" \
 	--output-geojson "./data/interim/impacts_forecasted_mean_slope_sto06_R_high.geojson"
 
+./data/interim/impacts_forecasted_mean_slope_hol86.geojson: ./data/interim/impacts_forecasted_mean_slope_hol86.csv ./data/interim/impacts_observed.csv
+	$(PYTHON_CLI) impacts-to-geojson \
+	--sites-csv "./data/interim/sites.csv" \
+	--observed-impacts-csv "./data/interim/impacts_observed.csv" \
+	--forecast-impacts-csv "./data/interim/impacts_forecasted_mean_slope_hol86.csv" \
+	--output-geojson "./data/interim/impacts_forecasted_mean_slope_hol86.geojson"
+
+./data/interim/impacts_forecasted_mean_slope_hol86_R_high.geojson: ./data/interim/impacts_forecasted_mean_slope_hol86.csv
+	$(PYTHON_CLI) r-high-to-geojson \
+	--sites-csv "./data/interim/sites.csv" \
+	--profiles-csv "./data/interim/profiles.csv" \
+	--crest-toes-csv "./data/interim/profile_features_crest_toes.csv" \
+	--impacts-csv "./data/interim/impacts_forecasted_mean_slope_hol86.csv" \
+	--output-geojson "./data/interim/impacts_forecasted_mean_slope_hol86_R_high.geojson"
+
 ./data/interim/profile_features_crest_toes.geojson: ./data/interim/profile_features_crest_toes.csv
 	$(PYTHON_CLI) profile-features-crest-toes-to-geojson \
 	--sites-csv "./data/interim/sites.csv" \

src/analysis/forecast_twl.py

@@ -41,13 +41,16 @@ def forecast_twl(
         # Process each site_id with a different process and combine results at the end
         with Pool(processes=n_processes) as pool:
             results = pool.starmap(
-                foreshore_slope_for_site_id, [(site_id, df_twl, df_profiles) for site_id in site_ids]
+                foreshore_slope_for_site_id,
+                [(site_id, df_twl, df_profiles) for site_id in site_ids],
             )
         df_twl["beta"] = pd.concat(results)
 
     elif slope == "mean":
         df_temp = df_twl.join(
-            df_profile_features.query("profile_type=='{}'".format(profile_type)).reset_index(level="profile_type"),
+            df_profile_features.query(
+                "profile_type=='{}'".format(profile_type)
+            ).reset_index(level="profile_type"),
             how="inner",
         )
         df_temp["mhw"] = 0.5
@@ -59,19 +62,26 @@ def forecast_twl(
             df_temp.dune_toe_z.isnull(), "dune_crest_z"
         ]
         df_temp["top_x"] = df_temp["dune_toe_x"]
-        df_temp.loc[df_temp.dune_toe_x.isnull(), "top_x"] = df_temp.loc[df_temp.dune_toe_x.isnull(), "dune_crest_x"]
+        df_temp.loc[df_temp.dune_toe_x.isnull(), "top_x"] = df_temp.loc[
+            df_temp.dune_toe_x.isnull(), "dune_crest_x"
+        ]
 
         with Pool(processes=n_processes) as pool:
             results = pool.starmap(
                 mean_slope_for_site_id,
-                [(site_id, df_temp, df_profiles, "top_elevation", "top_x", "mhw") for site_id in site_ids],
+                [
+                    (site_id, df_temp, df_profiles, "top_elevation", "top_x", "mhw")
+                    for site_id in site_ids
+                ],
             )
 
         df_twl["beta"] = pd.concat(results)
 
     # Estimate runup
     R2, setup, S_total, S_inc, S_ig = runup_function(
-        Hs0=df_twl["Hs0"].tolist(), Tp=df_twl["Tp"].tolist(), beta=df_twl["beta"].tolist()
+        Hs0=df_twl["Hs0"].tolist(),
+        Tp=df_twl["Tp"].tolist(),
+        beta=df_twl["beta"].tolist(),
     )
 
     df_twl["R2"] = R2
@@ -80,7 +90,9 @@ def forecast_twl(
 
     # Estimate TWL
     df_twl["R_high"] = df_twl["tide"] + df_twl["R2"]
-    df_twl["R_low"] = df_twl["tide"] + 1.1 * df_twl["setup"] - 1.1 / 2 * df_twl["S_total"]
+    df_twl["R_low"] = (
+        df_twl["tide"] + 1.1 * df_twl["setup"] - 1.1 / 2 * df_twl["S_total"]
+    )
 
     # Drop unneeded columns
     # df_twl.drop(columns=["E", "Exs", "P", "Pxs", "dir"], inplace=True, errors="ignore")
@@ -89,7 +101,13 @@ def forecast_twl(
 
 
 def mean_slope_for_site_id(
-    site_id, df_twl, df_profiles, top_elevation_col, top_x_col, btm_elevation_col, profile_type="prestorm"
+    site_id,
+    df_twl,
+    df_profiles,
+    top_elevation_col,
+    top_x_col,
+    btm_elevation_col,
+    profile_type="prestorm",
 ):
     """
     Calculates the foreshore slope values a given site_id. Returns a series (with same indicies as df_twl) of
@@ -135,7 +153,9 @@ def foreshore_slope_for_site_id(site_id, df_twl, df_profiles):
     """
 
     # Get the prestorm beach profile
-    profile = df_profiles.query("site_id =='{}' and profile_type == 'prestorm'".format(site_id))
+    profile = df_profiles.query(
+        "site_id =='{}' and profile_type == 'prestorm'".format(site_id)
+    )
     profile_x = profile.index.get_level_values("x").tolist()
     profile_z = profile.z.tolist()
 
@@ -175,8 +195,12 @@ def foreshore_slope_from_profile(profile_x, profile_z, tide, runup_function, **k
     # Initalize estimates
     max_number_iterations = 30
     iteration_count = 0
-    averaged_accuracy = 0.03  # if slopes within this amount, average after max number of iterations
+    averaged_accuracy = (
-    acceptable_accuracy = 0.01  # if slopes within this amount, accept after max number of iterations
+        0.03
+    )  # if slopes within this amount, average after max number of iterations
+    acceptable_accuracy = (
+        0.01
+    )  # if slopes within this amount, accept after max number of iterations
     preferred_accuracy = 0.001  # if slopes within this amount, accept
     beta = 0.05
 
@@ -212,7 +236,15 @@ def foreshore_slope_from_profile(profile_x, profile_z, tide, runup_function, **k
         iteration_count += 1
 
 
-def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, method="end_points", top_x=None, btm_x=None):
+def slope_from_profile(
+    profile_x,
+    profile_z,
+    top_elevation,
+    btm_elevation,
+    method="end_points",
+    top_x=None,
+    btm_x=None,
+):
     """
     Returns a slope (beta) from a bed profile, given the top and bottom elevations of where the slope should be taken.
     :param x: List of x bed profile coordinates
@@ -260,7 +292,9 @@ def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, metho
                 end_points[end_type]["x"] = intersection_x[-1]
             else:
                 # For bottom elevation, take most landward intersection that is seaward of top elevation
-                end_point_btm = [x for x in intersection_x if x > end_points["top"]["x"]]
+                end_point_btm = [
+                    x for x in intersection_x if x > end_points["top"]["x"]
+                ]
                 if len(end_point_btm) == 0:
                     # If there doesn't seem to be an intersection seaward of the top elevation, return none.
                     return None
@@ -275,7 +309,10 @@ def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, metho
         return -(z_top - z_btm) / (x_top - x_btm)
 
     elif method == "least_squares":
-        profile_mask = [True if end_points["top"]["x"] < pts < end_points["btm"]["x"] else False for pts in profile_x]
+        profile_mask = [
+            True if end_points["top"]["x"] < pts < end_points["btm"]["x"] else False
+            for pts in profile_x
+        ]
         slope_x = np.array(profile_x)[profile_mask].tolist()
         slope_z = np.array(profile_z)[profile_mask].tolist()
         slope, _, _, _, _ = stats.linregress(slope_x, slope_z)
@@ -287,12 +324,28 @@ def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, metho
 @click.option("--tides-csv", required=True, help="")
 @click.option("--profiles-csv", required=True, help="")
 @click.option("--profile-features-csv", required=True, help="")
-@click.option("--runup-function", required=True, help="", type=click.Choice(["sto06"]))
+@click.option(
-@click.option("--slope", required=True, help="", type=click.Choice(["foreshore", "mean"]))
+    "--runup-function", required=True, help="", type=click.Choice(["sto06", "hol86"])
-@click.option("--profile-type", required=True, help="", type=click.Choice(["prestorm", "poststorm"]))
+)
+@click.option(
+    "--slope", required=True, help="", type=click.Choice(["foreshore", "mean"])
+)
+@click.option(
+    "--profile-type",
+    required=True,
+    help="",
+    type=click.Choice(["prestorm", "poststorm"]),
+)
 @click.option("--output-file", required=True, help="")
 def create_twl_forecast(
-    waves_csv, tides_csv, profiles_csv, profile_features_csv, runup_function, slope, profile_type, output_file
+    waves_csv,
+    tides_csv,
+    profiles_csv,
+    profile_features_csv,
+    runup_function,
+    slope,
+    profile_type,
+    output_file,
 ):
     logger.info("Creating forecast of total water levels")
     logger.info("Importing data")

src/analysis/runup_models.py

@@ -10,14 +10,20 @@ def sto06(Hs0, Tp, beta):
     :return: Float or list of R2, setup, S_total, S_inc and S_ig values
     """
 
-    df = pd.DataFrame({"Hs0": Hs0, "Tp": Tp, "beta": beta}, index=[x for x in range(0, np.size(Hs0))])
+    df = pd.DataFrame(
+        {"Hs0": Hs0, "Tp": Tp, "beta": beta}, index=[x for x in range(0, np.size(Hs0))]
+    )
 
     df["Lp"] = 9.8 * df["Tp"] ** 2 / 2 / np.pi
 
     # General equation
     df["S_ig"] = pd.to_numeric(0.06 * np.sqrt(df["Hs0"] * df["Lp"]), errors="coerce")
-    df["S_inc"] = pd.to_numeric(0.75 * df["beta"] * np.sqrt(df["Hs0"] * df["Lp"]), errors="coerce")
+    df["S_inc"] = pd.to_numeric(
-    df["setup"] = pd.to_numeric(0.35 * df["beta"] * np.sqrt(df["Hs0"] * df["Lp"]), errors="coerce")
+        0.75 * df["beta"] * np.sqrt(df["Hs0"] * df["Lp"]), errors="coerce"
+    )
+    df["setup"] = pd.to_numeric(
+        0.35 * df["beta"] * np.sqrt(df["Hs0"] * df["Lp"]), errors="coerce"
+    )
     df["S_total"] = np.sqrt(df["S_inc"] ** 2 + df["S_ig"] ** 2)
     df["R2"] = 1.1 * (df["setup"] + df["S_total"] / 2)
 
@@ -37,6 +43,30 @@ def sto06(Hs0, Tp, beta):
     )
 
 
+def hol86(Hs0, Tp, beta):
+
+    df = pd.DataFrame(
+        {"Hs0": Hs0, "Tp": Tp, "beta": beta}, index=[x for x in range(0, np.size(Hs0))]
+    )
+
+    df["Lp"] = 9.8 * df["Tp"] ** 2 / 2 / np.pi
+
+    df["setup"] = 0.2 * df["Hs0"]
+    df["R2"] = 0.83 * df["beta"] * np.sqrt(df["Hs0"] * df["Lp"]) + df["setup"]
+
+    df["S_ig"] = np.nan
+    df["S_inc"] = np.nan
+    df["S_total"] = np.nan
+
+    return (
+        float_or_list(df["R2"].tolist()),
+        float_or_list(df["setup"].tolist()),
+        float_or_list(df["S_total"].tolist()),
+        float_or_list(df["S_inc"].tolist()),
+        float_or_list(df["S_ig"].tolist()),
+    )
+
+
 def float_or_list(a):
     """
     If only one value in the array, return the float, else return a list