Extract slope & width change for observed impacts

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()
 
@@ -255,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(
@@ -271,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