Extract function for getting slopes

This is so it can be called from observed_storm_impacts.py. It might be better to put these in a slope_utils.py function?

src/analysis/forecast_twl.py

@@ -49,26 +49,49 @@ def forecast_twl(
         df_twl["beta"] = pd.concat(results)
 
     elif slope == "mean":
-        logger.info("Calculating mean (dune toe to MHW) slopes")
+        df_slopes = get_mean_slope(df_profile_features, df_profiles, profile_type)
-        btm_z = 0.5  # m AHD
 
-        # When calculating mean slope, we go from the dune toe to mhw. However, in some profiles, the dune toe is not
+        # Merge calculated slopes onto each twl timestep
-        # defined. In these cases, we should go to the dune crest. Let's make a temporary dataframe which has this
+        df_twl = df_twl.merge(df_slopes, left_index=True, right_index=True)
-        # already calculated.
-        df_top_ele = df_profile_features.xs(profile_type, level="profile_type").copy()
-        df_top_ele.loc[:, "top_ele"] = df_top_ele.dune_toe_z
-        df_top_ele.loc[
-            df_top_ele.top_ele.isnull().values, "top_ele"
-        ] = df_top_ele.dune_crest_z
 
-        n_no_top_ele = len(df_top_ele[df_top_ele.top_ele.isnull()].index)
+    elif slope == "intertidal":
-        if n_no_top_ele != 0:
+
-            logger.warning(
+        logger.info("Calculating intertidal slopes")
-                "{} sites do not have dune toes/crests to calculate mean slope".format(
+        df_slopes = get_intertidal_slope(df_profiles, profile_type)
-                    n_no_top_ele
+
+        # Merge calculated slopes onto each twl timestep
+        df_twl = df_twl.merge(df_slopes, left_index=True, right_index=True)
+
+    # 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(),
+        r=df_twl.merge(df_grain_size, on="site_id").r.tolist(),
     )
+
+    df_twl["R2"] = R2
+    df_twl["setup"] = setup
+    df_twl["S_total"] = S_total
+
+    # 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"]
     )
 
+    return df_twl
+
+
+def get_intertidal_slope(df_profiles, profile_type, top_z=1.15, btm_z=-0.9):
+    """
+    Gets intertidal slopes
+    :param df_profiles:
+    :param profile_type:
+    :return:
+    """
+
+    # Calculate slopes for each profile
     df_slopes = (
         df_profiles.xs(profile_type, level="profile_type")
         .dropna(subset=["z"])
@@ -77,25 +100,44 @@ def forecast_twl(
             lambda x: slope_from_profile(
                 profile_x=x.index.get_level_values("x").tolist(),
                 profile_z=x.z.tolist(),
-                    top_elevation=df_top_ele.loc[x.index[0][0], :].top_ele,
+                top_elevation=top_z,
-                    btm_elevation=btm_z,
+                btm_elevation=max(min(x.z), btm_z),
                 method="least_squares",
             )
         )
         .rename("beta")
         .to_frame()
     )
+    return df_slopes
 
-        # Merge calculated slopes onto each twl timestep
-        df_twl = df_twl.merge(df_slopes, left_index=True, right_index=True)
 
-    elif slope == "intertidal":
+def get_mean_slope(df_profile_features, df_profiles, profile_type, btm_z=0.5):
+    """
+    Calculates the mean slopes for all profiles
+    :param df_profile_features:
+    :param df_profiles:
+    :param profile_type:
+    :param btm_z: Typically mean high water
+    :return:
+    """
 
-        logger.info("Calculating intertidal slopes")
+    logger.info("Calculating mean (dune toe to MHW) slopes")
-        top_z = 1.15  # m AHD = HAT from MHL annual ocean tides summary report
-        btm_z = -0.9  # m AHD = HAT from MHL annual ocean tides summary report
 
-        # Calculate slopes for each profile
+    # When calculating mean slope, we go from the dune toe to mhw. However, in some profiles, the dune toe is not
+    # defined. In these cases, we should go to the dune crest. Let's make a temporary dataframe which has this
+    # already calculated.
+    df_top_ele = df_profile_features.xs(profile_type, level="profile_type").copy()
+    df_top_ele.loc[:, "top_ele"] = df_top_ele.dune_toe_z
+    df_top_ele.loc[
+        df_top_ele.top_ele.isnull().values, "top_ele"
+    ] = df_top_ele.dune_crest_z
+    n_no_top_ele = len(df_top_ele[df_top_ele.top_ele.isnull()].index)
+    if n_no_top_ele != 0:
+        logger.warning(
+            "{} sites do not have dune toes/crests to calculate mean slope".format(
+                n_no_top_ele
+            )
+        )
     df_slopes = (
         df_profiles.xs(profile_type, level="profile_type")
         .dropna(subset=["z"])
@@ -104,40 +146,15 @@ def forecast_twl(
             lambda x: slope_from_profile(
                 profile_x=x.index.get_level_values("x").tolist(),
                 profile_z=x.z.tolist(),
-                    top_elevation=top_z,
+                top_elevation=df_top_ele.loc[x.index[0][0], :].top_ele,
-                    btm_elevation=max(min(x.z), btm_z),
+                btm_elevation=btm_z,
                 method="least_squares",
             )
         )
         .rename("beta")
         .to_frame()
     )
-
+    return df_slopes
-        # Merge calculated slopes onto each twl timestep
-        df_twl = df_twl.merge(df_slopes, left_index=True, right_index=True)
-
-    # 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(),
-        r=df_twl.merge(df_grain_size, on="site_id").r.tolist(),
-    )
-
-    df_twl["R2"] = R2
-    df_twl["setup"] = setup
-    df_twl["S_total"] = S_total
-
-    # 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"]
-    )
-
-    # Drop unneeded columns
-    # df_twl.drop(columns=["E", "Exs", "P", "Pxs", "dir"], inplace=True, errors="ignore")
-
-    return df_twl
 
 
 def mean_slope_for_site_id(