From 0adbf6863522ee4dbb460ce83eb935eb1e10b6a1 Mon Sep 17 00:00:00 2001 From: Chris Leaman Date: Fri, 8 Feb 2019 16:31:48 +1100 Subject: [PATCH] Add Holman 1986 to list of runup functions --- Makefile | 32 ++++++++++++- src/analysis/forecast_twl.py | 87 +++++++++++++++++++++++++++++------- src/analysis/runup_models.py | 36 +++++++++++++-- 3 files changed, 134 insertions(+), 21 deletions(-) diff --git a/Makefile b/Makefile index e9574a1..3bbede6 100644 --- a/Makefile +++ b/Makefile @@ -110,6 +110,16 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte --profile-type "prestorm" \ --output-file "./data/interim/twl_mean_slope_sto06.csv" +./data/interim/twl_mean_slope_hol86.csv: ./data/interim/waves.csv ./data/interim/tides.csv ./data/interim/profiles.csv ./data/interim/sites.csv ./data/interim/profile_features_crest_toes.csv + $(PYTHON_CLI) create-twl-forecast \ + --waves-csv "./data/interim/waves.csv" \ + --tides-csv "./data/interim/tides.csv" \ + --profiles-csv "./data/interim/profiles.csv" \ + --profile-features-csv "./data/interim/profile_features_crest_toes.csv" \ + --runup-function "hol86" \ + --slope "mean" \ + --profile-type "prestorm" \ + --output-file "./data/interim/twl_mean_slope_hol86.csv" ### IMPACTS @@ -132,10 +142,15 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte --forecasted-twl-csv "./data/interim/twl_foreshore_slope_sto06.csv" \ --output-file "./data/interim/impacts_forecasted_foreshore_slope_sto06.csv" +./data/interim/impacts_forecasted_mean_slope_hol86.csv: ./data/interim/profile_features_crest_toes.csv ./data/interim/twl_mean_slope_hol86.csv + $(PYTHON_CLI) create-forecasted-impacts \ + --profile-features-csv "./data/interim/profile_features_crest_toes.csv" \ + --forecasted-twl-csv "./data/interim/twl_mean_slope_hol86.csv" \ + --output-file "./data/interim/impacts_forecasted_mean_slope_hol86.csv" ### GEOJSONs -geojsons: ./data/interim/impacts_forecasted_mean_slope_sto06.geojson ./data/interim/impacts_forecasted_mean_slope_sto06_R_high.geojson ./data/interim/profile_features_crest_toes.geojson ./data/interim/sites.geojson +geojsons: ./data/interim/impacts_forecasted_mean_slope_hol86.geojson ./data/interim/impacts_forecasted_mean_slope_hol86_R_high.geojson ./data/interim/impacts_forecasted_mean_slope_sto06.geojson ./data/interim/impacts_forecasted_mean_slope_sto06_R_high.geojson ./data/interim/profile_features_crest_toes.geojson ./data/interim/sites.geojson ./data/interim/impacts_forecasted_mean_slope_sto06.geojson: ./data/interim/impacts_forecasted_mean_slope_sto06.csv ./data/interim/impacts_observed.csv $(PYTHON_CLI) impacts-to-geojson \ @@ -152,6 +167,21 @@ geojsons: ./data/interim/impacts_forecasted_mean_slope_sto06.geojson ./data/inte --impacts-csv "./data/interim/impacts_forecasted_mean_slope_sto06.csv" \ --output-geojson "./data/interim/impacts_forecasted_mean_slope_sto06_R_high.geojson" +./data/interim/impacts_forecasted_mean_slope_hol86.geojson: ./data/interim/impacts_forecasted_mean_slope_hol86.csv ./data/interim/impacts_observed.csv + $(PYTHON_CLI) impacts-to-geojson \ + --sites-csv "./data/interim/sites.csv" \ + --observed-impacts-csv "./data/interim/impacts_observed.csv" \ + --forecast-impacts-csv "./data/interim/impacts_forecasted_mean_slope_hol86.csv" \ + --output-geojson "./data/interim/impacts_forecasted_mean_slope_hol86.geojson" + +./data/interim/impacts_forecasted_mean_slope_hol86_R_high.geojson: ./data/interim/impacts_forecasted_mean_slope_hol86.csv + $(PYTHON_CLI) r-high-to-geojson \ + --sites-csv "./data/interim/sites.csv" \ + --profiles-csv "./data/interim/profiles.csv" \ + --crest-toes-csv "./data/interim/profile_features_crest_toes.csv" \ + --impacts-csv "./data/interim/impacts_forecasted_mean_slope_hol86.csv" \ + --output-geojson "./data/interim/impacts_forecasted_mean_slope_hol86_R_high.geojson" + ./data/interim/profile_features_crest_toes.geojson: ./data/interim/profile_features_crest_toes.csv $(PYTHON_CLI) profile-features-crest-toes-to-geojson \ --sites-csv "./data/interim/sites.csv" \ diff --git a/src/analysis/forecast_twl.py b/src/analysis/forecast_twl.py index 5298caa..9badbe2 100644 --- a/src/analysis/forecast_twl.py +++ b/src/analysis/forecast_twl.py @@ -41,13 +41,16 @@ def forecast_twl( # Process each site_id with a different process and combine results at the end with Pool(processes=n_processes) as pool: results = pool.starmap( - foreshore_slope_for_site_id, [(site_id, df_twl, df_profiles) for site_id in site_ids] + foreshore_slope_for_site_id, + [(site_id, df_twl, df_profiles) for site_id in site_ids], ) df_twl["beta"] = pd.concat(results) elif slope == "mean": df_temp = df_twl.join( - df_profile_features.query("profile_type=='{}'".format(profile_type)).reset_index(level="profile_type"), + df_profile_features.query( + "profile_type=='{}'".format(profile_type) + ).reset_index(level="profile_type"), how="inner", ) df_temp["mhw"] = 0.5 @@ -59,19 +62,26 @@ def forecast_twl( df_temp.dune_toe_z.isnull(), "dune_crest_z" ] df_temp["top_x"] = df_temp["dune_toe_x"] - df_temp.loc[df_temp.dune_toe_x.isnull(), "top_x"] = df_temp.loc[df_temp.dune_toe_x.isnull(), "dune_crest_x"] + df_temp.loc[df_temp.dune_toe_x.isnull(), "top_x"] = df_temp.loc[ + df_temp.dune_toe_x.isnull(), "dune_crest_x" + ] with Pool(processes=n_processes) as pool: results = pool.starmap( mean_slope_for_site_id, - [(site_id, df_temp, df_profiles, "top_elevation", "top_x", "mhw") for site_id in site_ids], + [ + (site_id, df_temp, df_profiles, "top_elevation", "top_x", "mhw") + for site_id in site_ids + ], ) df_twl["beta"] = pd.concat(results) # 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() + Hs0=df_twl["Hs0"].tolist(), + Tp=df_twl["Tp"].tolist(), + beta=df_twl["beta"].tolist(), ) df_twl["R2"] = R2 @@ -80,7 +90,9 @@ def forecast_twl( # 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"] + 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") @@ -89,7 +101,13 @@ def forecast_twl( def mean_slope_for_site_id( - site_id, df_twl, df_profiles, top_elevation_col, top_x_col, btm_elevation_col, profile_type="prestorm" + site_id, + df_twl, + df_profiles, + top_elevation_col, + top_x_col, + btm_elevation_col, + profile_type="prestorm", ): """ Calculates the foreshore slope values a given site_id. Returns a series (with same indicies as df_twl) of @@ -135,7 +153,9 @@ def foreshore_slope_for_site_id(site_id, df_twl, df_profiles): """ # Get the prestorm beach profile - profile = df_profiles.query("site_id =='{}' and profile_type == 'prestorm'".format(site_id)) + profile = df_profiles.query( + "site_id =='{}' and profile_type == 'prestorm'".format(site_id) + ) profile_x = profile.index.get_level_values("x").tolist() profile_z = profile.z.tolist() @@ -175,8 +195,12 @@ def foreshore_slope_from_profile(profile_x, profile_z, tide, runup_function, **k # Initalize estimates max_number_iterations = 30 iteration_count = 0 - averaged_accuracy = 0.03 # if slopes within this amount, average after max number of iterations - acceptable_accuracy = 0.01 # if slopes within this amount, accept after max number of iterations + averaged_accuracy = ( + 0.03 + ) # if slopes within this amount, average after max number of iterations + acceptable_accuracy = ( + 0.01 + ) # if slopes within this amount, accept after max number of iterations preferred_accuracy = 0.001 # if slopes within this amount, accept beta = 0.05 @@ -212,7 +236,15 @@ def foreshore_slope_from_profile(profile_x, profile_z, tide, runup_function, **k iteration_count += 1 -def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, method="end_points", top_x=None, btm_x=None): +def slope_from_profile( + profile_x, + profile_z, + top_elevation, + btm_elevation, + method="end_points", + top_x=None, + btm_x=None, +): """ Returns a slope (beta) from a bed profile, given the top and bottom elevations of where the slope should be taken. :param x: List of x bed profile coordinates @@ -260,7 +292,9 @@ def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, metho end_points[end_type]["x"] = intersection_x[-1] else: # For bottom elevation, take most landward intersection that is seaward of top elevation - end_point_btm = [x for x in intersection_x if x > end_points["top"]["x"]] + end_point_btm = [ + x for x in intersection_x if x > end_points["top"]["x"] + ] if len(end_point_btm) == 0: # If there doesn't seem to be an intersection seaward of the top elevation, return none. return None @@ -275,7 +309,10 @@ def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, metho return -(z_top - z_btm) / (x_top - x_btm) elif method == "least_squares": - profile_mask = [True if end_points["top"]["x"] < pts < end_points["btm"]["x"] else False for pts in profile_x] + profile_mask = [ + True if end_points["top"]["x"] < pts < end_points["btm"]["x"] else False + for pts in profile_x + ] slope_x = np.array(profile_x)[profile_mask].tolist() slope_z = np.array(profile_z)[profile_mask].tolist() slope, _, _, _, _ = stats.linregress(slope_x, slope_z) @@ -287,12 +324,28 @@ def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, metho @click.option("--tides-csv", required=True, help="") @click.option("--profiles-csv", required=True, help="") @click.option("--profile-features-csv", required=True, help="") -@click.option("--runup-function", required=True, help="", type=click.Choice(["sto06"])) -@click.option("--slope", required=True, help="", type=click.Choice(["foreshore", "mean"])) -@click.option("--profile-type", required=True, help="", type=click.Choice(["prestorm", "poststorm"])) +@click.option( + "--runup-function", required=True, help="", type=click.Choice(["sto06", "hol86"]) +) +@click.option( + "--slope", required=True, help="", type=click.Choice(["foreshore", "mean"]) +) +@click.option( + "--profile-type", + required=True, + help="", + type=click.Choice(["prestorm", "poststorm"]), +) @click.option("--output-file", required=True, help="") def create_twl_forecast( - waves_csv, tides_csv, profiles_csv, profile_features_csv, runup_function, slope, profile_type, output_file + waves_csv, + tides_csv, + profiles_csv, + profile_features_csv, + runup_function, + slope, + profile_type, + output_file, ): logger.info("Creating forecast of total water levels") logger.info("Importing data") diff --git a/src/analysis/runup_models.py b/src/analysis/runup_models.py index c508f65..374858f 100644 --- a/src/analysis/runup_models.py +++ b/src/analysis/runup_models.py @@ -10,14 +10,20 @@ def sto06(Hs0, Tp, beta): :return: Float or list of R2, setup, S_total, S_inc and S_ig values """ - df = pd.DataFrame({"Hs0": Hs0, "Tp": Tp, "beta": beta}, index=[x for x in range(0, np.size(Hs0))]) + df = pd.DataFrame( + {"Hs0": Hs0, "Tp": Tp, "beta": beta}, index=[x for x in range(0, np.size(Hs0))] + ) df["Lp"] = 9.8 * df["Tp"] ** 2 / 2 / np.pi # General equation df["S_ig"] = pd.to_numeric(0.06 * np.sqrt(df["Hs0"] * df["Lp"]), errors="coerce") - df["S_inc"] = pd.to_numeric(0.75 * df["beta"] * np.sqrt(df["Hs0"] * df["Lp"]), errors="coerce") - df["setup"] = pd.to_numeric(0.35 * df["beta"] * np.sqrt(df["Hs0"] * df["Lp"]), errors="coerce") + df["S_inc"] = pd.to_numeric( + 0.75 * df["beta"] * np.sqrt(df["Hs0"] * df["Lp"]), errors="coerce" + ) + df["setup"] = pd.to_numeric( + 0.35 * df["beta"] * np.sqrt(df["Hs0"] * df["Lp"]), errors="coerce" + ) df["S_total"] = np.sqrt(df["S_inc"] ** 2 + df["S_ig"] ** 2) df["R2"] = 1.1 * (df["setup"] + df["S_total"] / 2) @@ -37,6 +43,30 @@ def sto06(Hs0, Tp, beta): ) +def hol86(Hs0, Tp, beta): + + df = pd.DataFrame( + {"Hs0": Hs0, "Tp": Tp, "beta": beta}, index=[x for x in range(0, np.size(Hs0))] + ) + + df["Lp"] = 9.8 * df["Tp"] ** 2 / 2 / np.pi + + df["setup"] = 0.2 * df["Hs0"] + df["R2"] = 0.83 * df["beta"] * np.sqrt(df["Hs0"] * df["Lp"]) + df["setup"] + + df["S_ig"] = np.nan + df["S_inc"] = np.nan + df["S_total"] = np.nan + + return ( + float_or_list(df["R2"].tolist()), + float_or_list(df["setup"].tolist()), + float_or_list(df["S_total"].tolist()), + float_or_list(df["S_inc"].tolist()), + float_or_list(df["S_ig"].tolist()), + ) + + def float_or_list(a): """ If only one value in the array, return the float, else return a list