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 "src/"
+	activate ./.venv && black --line-length 120 "src/"
 
 
 ###############################

src/analysis/forecast_twl.py

@@ -49,49 +49,26 @@ def forecast_twl(
         df_twl["beta"] = pd.concat(results)
 
     elif slope == "mean":
-        df_slopes = get_mean_slope(df_profile_features, df_profiles, profile_type)
+        logger.info("Calculating mean (dune toe to MHW) slopes")
-
+        btm_z = 0.5  # m AHD
-        # Merge calculated slopes onto each twl timestep
-        df_twl = df_twl.merge(df_slopes, left_index=True, right_index=True)
-
-    elif slope == "intertidal":
-
-        logger.info("Calculating intertidal slopes")
-        df_slopes = get_intertidal_slope(df_profiles, profile_type)
 
-        # Merge calculated slopes onto each twl timestep
+        # When calculating mean slope, we go from the dune toe to mhw. However, in some profiles, the dune toe is not
-        df_twl = df_twl.merge(df_slopes, left_index=True, right_index=True)
+        # 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
 
-    # Estimate runup
+        n_no_top_ele = len(df_top_ele[df_top_ele.top_ele.isnull()].index)
-    R2, setup, S_total, S_inc, S_ig = runup_function(
+        if n_no_top_ele != 0:
-        Hs0=df_twl["Hs0"].tolist(),
+            logger.warning(
-        Tp=df_twl["Tp"].tolist(),
+                "{} sites do not have dune toes/crests to calculate mean slope".format(
-        beta=df_twl["beta"].tolist(),
+                    n_no_top_ele
-        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"])
@@ -100,44 +77,25 @@ def get_intertidal_slope(df_profiles, profile_type, top_z=1.15, btm_z=-0.9):
                 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)
 
-def get_mean_slope(df_profile_features, df_profiles, profile_type, btm_z=0.5):
+    elif slope == "intertidal":
-    """
-    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 mean (dune toe to MHW) slopes")
+        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
 
-    # When calculating mean slope, we go from the dune toe to mhw. However, in some profiles, the dune toe is not
+        # Calculate slopes for each profile
-    # 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"])
@@ -146,15 +104,40 @@ def get_mean_slope(df_profile_features, df_profiles, profile_type, btm_z=0.5):
                 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)
+
+    # 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(

src/analysis/observed_storm_impacts.py

@@ -4,8 +4,7 @@ import pandas as pd
 from scipy.integrate import simps
 
 from logs import setup_logging
-from utils import crossings, get_i_or_default
+from utils import crossings
-from analysis.forecast_twl import get_mean_slope, get_intertidal_slope
 
 logger = setup_logging()
 
@@ -219,21 +218,9 @@ def overwrite_impacts(df_observed_impacts, df_raw_features):
     :param df_raw_profile_features:
     :return:
     """
-
+    df_observed_impacts.update(
-    # Get manually specified impacts from the profile features ./data/raw/ folder. Note that sites which need to be
+        df_raw_features.rename(columns={"observed_storm_regime": "storm_regime"})
-    # 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
 
 
@@ -256,7 +243,6 @@ 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(
@@ -273,118 +259,10 @@ 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)
 
-    # Calculate change in mean slope
+    # TODO Calculate change in slopes, shoreline and volume
-    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,10 +71,3 @@ 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