From 217e633ffc1777cd06fe79c0695f3b659ef671f0 Mon Sep 17 00:00:00 2001 From: Chris Leaman Date: Mon, 11 Mar 2019 12:07:22 +1100 Subject: [PATCH] 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 | 135 ++++++++++++++++++++--------------- 1 file changed, 76 insertions(+), 59 deletions(-) diff --git a/src/analysis/forecast_twl.py b/src/analysis/forecast_twl.py index cd02071..70edcab 100644 --- a/src/analysis/forecast_twl.py +++ b/src/analysis/forecast_twl.py @@ -49,42 +49,7 @@ 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 - - # 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"]) - .groupby("site_id") - .apply( - 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, - method="least_squares", - ) - ) - .rename("beta") - .to_frame() - ) + df_slopes = get_mean_slope(df_profile_features, df_profiles, profile_type) # Merge calculated slopes onto each twl timestep df_twl = df_twl.merge(df_slopes, left_index=True, right_index=True) @@ -92,26 +57,7 @@ def forecast_twl( elif slope == "intertidal": 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 - - # Calculate slopes for each profile - df_slopes = ( - df_profiles.xs(profile_type, level="profile_type") - .dropna(subset=["z"]) - .groupby("site_id") - .apply( - 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), - method="least_squares", - ) - ) - .rename("beta") - .to_frame() - ) + 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) @@ -134,12 +80,83 @@ def forecast_twl( 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 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"]) + .groupby("site_id") + .apply( + 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), + method="least_squares", + ) + ) + .rename("beta") + .to_frame() + ) + return df_slopes + + +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 mean (dune toe to MHW) slopes") + + # 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"]) + .groupby("site_id") + .apply( + 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, + method="least_squares", + ) + ) + .rename("beta") + .to_frame() + ) + return df_slopes + + def mean_slope_for_site_id( site_id, df_twl,