Extract slope & width change for observed impacts

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?

Makefile

@@ -723,7 +723,7 @@ geojsons: ./data/interim/profile_features_crest_toes.geojson ./data/interim/site
 ###############################
 # Misc commands
 format: ./src/*.py ##@misc Check python file formatting
-	activate ./.venv && black --line-length 120 "src/"
+	activate ./.venv && black "src/"
 
 
 ###############################

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")
-        btm_z = 0.5  # m AHD
+        df_slopes = get_mean_slope(df_profile_features, df_profiles, profile_type)
 
-        # 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
+        # Merge calculated slopes onto each twl timestep
+        df_twl = df_twl.merge(df_slopes, left_index=True, right_index=True)
 
-        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
+    elif slope == "intertidal":
+
+        logger.info("Calculating intertidal slopes")
+        df_slopes = get_intertidal_slope(df_profiles, profile_type)
+
+        # 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,
-                    btm_elevation=btm_z,
+                top_elevation=top_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")
-        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
+    logger.info("Calculating mean (dune toe to MHW) slopes")
 
-        # 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,
-                    btm_elevation=max(min(x.z), btm_z),
+                top_elevation=df_top_ele.loc[x.index[0][0], :].top_ele,
+                btm_elevation=btm_z,
                 method="least_squares",
             )
         )
         .rename("beta")
         .to_frame()
     )
-
-        # 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
+    return df_slopes
 
 
 def mean_slope_for_site_id(

src/analysis/observed_storm_impacts.py

@@ -4,7 +4,8 @@ import pandas as pd
 from scipy.integrate import simps
 
 from logs import setup_logging
-from utils import crossings
+from utils import crossings, get_i_or_default
+from analysis.forecast_twl import get_mean_slope, get_intertidal_slope
 
 logger = setup_logging()
 
@@ -218,9 +219,21 @@ def overwrite_impacts(df_observed_impacts, df_raw_features):
     :param df_raw_profile_features:
     :return:
     """
-    df_observed_impacts.update(
-        df_raw_features.rename(columns={"observed_storm_regime": "storm_regime"})
-    )
+
+    # Get manually specified impacts from the profile features ./data/raw/ folder. Note that sites which need to be
+    # overwritten with a NaN, use the string 'none' in the csv. This is because when we use the df.update() command,
+    # it doesn't overwrite NaN values. So we'll put in the 'none' string, then overwrite that with the NaN.
+
+    df_overwritten_impacts = df_raw_features.rename(
+        columns={"observed_storm_regime": "storm_regime"}
+    ).storm_regime.to_frame()
+
+    df_observed_impacts.update(df_overwritten_impacts)
+
+    # Replace 'none' with nan
+    df_overwritten_impacts.loc[
+        df_overwritten_impacts.storm_regime == "unknown", "storm_regime"
+    ] = np.nan
     return df_observed_impacts
 
 
@@ -243,6 +256,7 @@ def create_observed_impacts(
         index=df_profile_features.index.get_level_values("site_id").unique()
     )
 
+    # TODO Review volume change with changing dune toe/crests
     logger.info("Getting pre/post storm volumes")
     df_swash_vol_changes = volume_change(df_profiles, df_profile_features, zone="swash")
     df_dune_face_vol_changes = volume_change(
@@ -259,10 +273,118 @@ def create_observed_impacts(
     df_raw_features = pd.read_csv(raw_profile_features_csv, index_col=[0])
     df_observed_impacts = overwrite_impacts(df_observed_impacts, df_raw_features)
 
-    # TODO Calculate change in slopes, shoreline and volume
+    # Calculate change in mean slope
+    df_prestorm_mean_slopes = get_mean_slope(
+        df_profile_features, df_profiles, profile_type="prestorm"
+    )
+    df_poststorm_mean_slopes = get_mean_slope(
+        df_profile_features, df_profiles, profile_type="poststorm"
+    )
+    df_diff_mean_slopes = df_poststorm_mean_slopes - df_prestorm_mean_slopes
+
+    # Calculate change in intertidal slope
+    df_prestorm_intertidal_slopes = get_intertidal_slope(
+        df_profiles, profile_type="prestorm"
+    )
+    df_poststorm_intertidal_slopes = get_intertidal_slope(
+        df_profiles, profile_type="poststorm"
+    )
+    df_diff_intertidal_slopes = (
+        df_poststorm_intertidal_slopes - df_prestorm_intertidal_slopes
+    )
+
+    # Rename slope columns and merge into observed impacts
+    renames = [
+        {"df": df_prestorm_mean_slopes, "new_col_name": "beta_prestorm_mean"},
+        {"df": df_poststorm_mean_slopes, "new_col_name": "beta_poststorm_mean"},
+        {"df": df_diff_mean_slopes, "new_col_name": "beta_diff_mean"},
+        {
+            "df": df_prestorm_intertidal_slopes,
+            "new_col_name": "beta_prestorm_intertidal",
+        },
+        {
+            "df": df_poststorm_intertidal_slopes,
+            "new_col_name": "beta_poststorm_intertidal",
+        },
+        {"df": df_diff_intertidal_slopes, "new_col_name": "beta_diff_intertidal"},
+    ]
+
+    for rename in renames:
+        rename["df"].rename(
+            {"beta": rename["new_col_name"]}, axis="columns", inplace=True
+        )
+
+    # Join all our slopes into the observed impacts
+    df_observed_impacts = pd.concat(
+        [
+            df_prestorm_mean_slopes,
+            df_poststorm_mean_slopes,
+            df_diff_mean_slopes,
+            df_prestorm_intertidal_slopes,
+            df_poststorm_intertidal_slopes,
+            df_diff_intertidal_slopes,
+            df_observed_impacts,
+        ],
+        axis=1,
+    )
+
+    # Calculate change in beach width
+    df_width_msl_prestorm = get_beach_width(
+        df_profile_features,
+        df_profiles,
+        profile_type="prestorm",
+        ele=0,
+        col_name="width_msl_prestorm",
+    )
+    df_width_msl_poststorm = get_beach_width(
+        df_profile_features,
+        df_profiles,
+        profile_type="poststorm",
+        ele=0,
+        col_name="width_msl_poststorm",
+    )
+    df_width_msl_change_m = (df_width_msl_poststorm - df_width_msl_prestorm).rename('df_width_msl_change_m')
+    df_width_msl_change_pct = (df_width_msl_change_m / df_width_msl_prestorm * 100).rename('df_width_msl_change_pct')
+
+    # Join beach width change onto observed impacts dataframe
+    df_observed_impacts = pd.concat(
+        [
+            df_observed_impacts,
+            df_width_msl_prestorm,
+            df_width_msl_poststorm,
+            df_width_msl_change_m,
+            df_width_msl_change_pct,
+        ],
+        axis=1,
+    )
 
     # Save dataframe to csv
     df_observed_impacts.to_csv(output_file, float_format="%.4f")
 
     logger.info("Saved to %s", output_file)
     logger.info("Done!")
+
+
+def get_beach_width(df_profile_features, df_profiles, profile_type, ele, col_name):
+    df_x_position = (
+        df_profiles.xs(profile_type, level="profile_type")
+        .dropna(subset=["z"])
+        .groupby("site_id")
+        .apply(
+            lambda x: get_i_or_default(
+                crossings(
+                    profile_x=x.index.get_level_values("x").tolist(),
+                    profile_z=x.z.tolist(),
+                    constant_z=ele,
+                ),
+                -1,
+                default=np.nan,
+            )
+        )
+        .rename("x_position")
+    )
+    df_x_prestorm_dune_toe = df_profile_features.xs(
+        "prestorm", level="profile_type"
+    ).dune_toe_x
+    df_width = (df_x_position - df_x_prestorm_dune_toe).rename(col_name)
+    return df_width

src/utils.py

@@ -71,3 +71,10 @@ def convert_coord_systems(
 
     g2 = transform(project, g1)  # apply projection
     return g2
+
+
+def get_i_or_default(l, i, default=None):
+    try:
+        return l[i]
+    except IndexError:
+        return default