diff --git a/src/analysis/forecast_twl.py b/src/analysis/forecast_twl.py index 82efc65..0fab400 100644 --- a/src/analysis/forecast_twl.py +++ b/src/analysis/forecast_twl.py @@ -23,7 +23,7 @@ def forecast_twl( runup_function, n_processes=MULTIPROCESS_THREADS, slope="foreshore", - profile_type='prestorm' + profile_type="prestorm", ): # Use df_waves as a base df_twl = df_waves.copy() @@ -46,19 +46,20 @@ def forecast_twl( 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') - , how="inner") + df_temp = df_twl.join( + df_profile_features.query("profile_type=='{}'".format(profile_type)).reset_index(level="profile_type"), + how="inner", + ) df_temp["mhw"] = 0.5 with Pool(processes=n_processes) as pool: results = pool.starmap( - mean_slope_for_site_id, [(site_id, df_temp, df_profiles, "dune_toe_z", "dune_toe_x", "mhw") for - site_id in site_ids] + mean_slope_for_site_id, + [(site_id, df_temp, df_profiles, "dune_toe_z", "dune_toe_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(df_twl, Hs0_col="Hs0", Tp_col="Tp", beta_col="beta") + 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()) df_twl["R2"] = R2 df_twl["setup"] = setup @@ -69,13 +70,14 @@ def forecast_twl( 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") + # df_twl.drop(columns=["E", "Exs", "P", "Pxs", "dir"], inplace=True, errors="ignore") return df_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'): +def mean_slope_for_site_id( + 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 foreshore slopes. This function is used to parallelize getting foreshore slopes as it is computationally @@ -100,7 +102,7 @@ def mean_slope_for_site_id(site_id, df_twl, df_profiles, top_elevation_col, top_ top_elevation=row[top_elevation_col], btm_elevation=row[btm_elevation_col], method="end_points", - top_x= row[top_x_col] + top_x=row[top_x_col], ), axis=1, ) @@ -130,7 +132,7 @@ def foreshore_slope_for_site_id(site_id, df_twl, df_profiles): profile_x=profile_x, profile_z=profile_z, tide=row.tide, - runup_function=runup_models.sto06_individual, + runup_function=runup_models.sto06, Hs0=row.Hs0, Tp=row.Tp, ), @@ -216,16 +218,14 @@ def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, metho end_points = {"top": {"z": top_elevation}, "btm": {"z": btm_elevation}} - - for end_type in end_points.keys(): # Add x coordinates if they are specified - if top_x and end_type == 'top': - end_points['top']['x'] = top_x + if top_x and end_type == "top": + end_points["top"]["x"] = top_x continue - if btm_x and end_type == 'top': - end_points['btm']['x'] = btm_x + if btm_x and end_type == "top": + end_points["btm"]["x"] = btm_x continue elevation = end_points[end_type]["z"] @@ -306,14 +306,15 @@ def crossings(profile_x, profile_z, constant_z): @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): +def create_twl_forecast( + 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") df_waves = pd.read_csv(waves_csv, index_col=[0, 1]) df_tides = pd.read_csv(tides_csv, index_col=[0, 1]) df_profiles = pd.read_csv(profiles_csv, index_col=[0, 1, 2]) - df_profile_features = pd.read_csv(profile_features_csv, index_col=[0,1]) + df_profile_features = pd.read_csv(profile_features_csv, index_col=[0, 1]) logger.info("Forecasting TWL") df_twl = forecast_twl( @@ -323,7 +324,7 @@ def create_twl_forecast(waves_csv, tides_csv, profiles_csv, profile_features_csv df_profile_features, runup_function=getattr(runup_models, runup_function), slope=slope, - profile_type=profile_type + profile_type=profile_type, ) df_twl.to_csv(output_file) diff --git a/src/analysis/forecasted_storm_impacts.py b/src/analysis/forecasted_storm_impacts.py index 51fb920..0a04ea1 100644 --- a/src/analysis/forecasted_storm_impacts.py +++ b/src/analysis/forecasted_storm_impacts.py @@ -20,7 +20,7 @@ def forecasted_impacts(df_profile_features, df_forecasted_twl): """ logger.info("Getting forecasted storm impacts") - df_forecasted_impacts = pd.DataFrame(index=df_profile_features.index.get_level_values('site_id').unique()) + df_forecasted_impacts = pd.DataFrame(index=df_profile_features.index.get_level_values("site_id").unique()) # For each site, find the maximum R_high value and the corresponding R_low value. idx = df_forecasted_twl.groupby(level=["site_id"])["R_high"].idxmax().dropna() @@ -29,9 +29,12 @@ def forecasted_impacts(df_profile_features, df_forecasted_twl): # Join with df_profile features to find dune toe and crest elevations df_forecasted_impacts = df_forecasted_impacts.merge( - df_profile_features.query("profile_type=='prestorm'").reset_index('profile_type')[["dune_toe_z", - "dune_crest_z"]], how="left", - left_on=['site_id'], right_on=['site_id'] + df_profile_features.query("profile_type=='prestorm'").reset_index("profile_type")[ + ["dune_toe_z", "dune_crest_z"] + ], + how="left", + left_on=["site_id"], + right_on=["site_id"], ) # Compare R_high and R_low wirth dune crest and toe elevations @@ -68,16 +71,34 @@ def storm_regime(df_forecasted_impacts): return df_forecasted_impacts -def twl_exceedence_time(df_profile_features, df_forecasted_twl, z_twl_col='R_high', z_exceedence_col='dune_toe_z'): +def twl_exceedence_time(df_profile_features, df_forecasted_twl, z_twl_col="R_high", z_exceedence_col="dune_toe_z"): + """ + Returns a dataframe of number of hours the twl exceeded a certain z elevation. + May need to use this https://stackoverflow.com/a/53656968 if datetimes are not consistent. + :param df_profile_features: + :param df_forecasted_twl: + :param z_twl_col: + :param z_exceedence_col: + :return: + """ logger.info("Getting twl exceedence time") - df_dune_toes = df_profile_features.query('profile_type=="prestorm"').reset_index('profile_type')[ - 'dune_toe_z'].to_frame() + # Get a dataframe of prestorm dune toes organised by site_id + df_dune_toes = ( + df_profile_features.query('profile_type=="prestorm"').reset_index("profile_type")[z_exceedence_col].to_frame() + ) - df_merged = df_forecasted_twl.merge(df_dune_toes,left_on=['site_id'],right_on=['site_id']) + # Merge dune toes into site_id + df_merged = df_forecasted_twl.merge(df_dune_toes, left_on=["site_id"], right_on=["site_id"]) - return (df_merged[z_twl_col] >= df_merged[z_exceedence_col]).groupby('site_id').sum().rename( - 'twl_{}_exceedance_hrs'.format(z_exceedence_col)).to_frame() + # Return the sum of hours that twl exceeded the level + return ( + (df_merged[z_twl_col] >= df_merged[z_exceedence_col]) + .groupby("site_id") + .sum() + .rename("twl_{}_exceedance_hrs".format(z_exceedence_col)) + .to_frame() + ) @click.command() @@ -87,14 +108,14 @@ def twl_exceedence_time(df_profile_features, df_forecasted_twl, z_twl_col='R_hig def create_forecasted_impacts(profile_features_csv, forecasted_twl_csv, output_file): logger.info("Creating observed wave impacts") logger.info("Importing existing data") - df_profile_features = pd.read_csv(profile_features_csv, index_col=[0,1]) + df_profile_features = pd.read_csv(profile_features_csv, index_col=[0, 1]) df_forecasted_twl = pd.read_csv(forecasted_twl_csv, index_col=[0, 1]) df_forecasted_impacts = forecasted_impacts(df_profile_features, df_forecasted_twl) - df_forecasted_impacts = df_profile_features.merge(twl_exceedence_time(df_profile_features, df_forecasted_twl), - left_on=['site_id'], - right_on=['site_id']) + df_forecasted_impacts = df_forecasted_impacts.merge( + twl_exceedence_time(df_profile_features, df_forecasted_twl), left_on=["site_id"], right_on=["site_id"] + ) df_forecasted_impacts.to_csv(output_file) logger.info("Saved to %s", output_file) diff --git a/src/analysis/observed_storm_impacts.py b/src/analysis/observed_storm_impacts.py index 6c76e36..6c4c0e1 100644 --- a/src/analysis/observed_storm_impacts.py +++ b/src/analysis/observed_storm_impacts.py @@ -30,13 +30,13 @@ def volume_change(df_profiles, df_profile_features, zone): """ logger.info("Calculating change in beach volume in {} zone".format(zone)) - df_vol_changes = pd.DataFrame(index=df_profile_features.index.get_level_values('site_id').unique()) + df_vol_changes = pd.DataFrame(index=df_profile_features.index.get_level_values("site_id").unique()) df_profiles = df_profiles.sort_index() sites = df_profiles.groupby(level=["site_id"]) for site_id, df_site in sites: logger.debug("Calculating change in beach volume at {} in {} zone".format(site_id, zone)) - query ="site_id=='{}'&profile_type=='prestorm'".format(site_id) + query = "site_id=='{}'&profile_type=='prestorm'".format(site_id) prestorm_dune_toe_x = df_profile_features.query(query).dune_toe_x.tolist() prestorm_dune_crest_x = df_profile_features.query(query).dune_crest_x.tolist() @@ -71,11 +71,11 @@ def volume_change(df_profiles, df_profile_features, zone): # Where we want to measure pre and post storm volume is dependant on the zone selected if zone == "swash": - x_min = max(prestorm_dune_toe_x,x_first_obs) + x_min = max(prestorm_dune_toe_x, x_first_obs) x_max = x_last_obs elif zone == "dune_face": x_min = max(prestorm_dune_crest_x, x_first_obs) - x_max = min(prestorm_dune_toe_x,x_last_obs) + x_max = min(prestorm_dune_toe_x, x_last_obs) else: logger.warning("Zone argument not properly specified. Please check") x_min = None @@ -98,16 +98,11 @@ def volume_change(df_profiles, df_profile_features, zone): # Volume change needs to be calculated including a tolerance for LIDAR accuracy. If difference between # profiles is less than 20 cm, consider them as zero difference. - prestorm_z = df_zone.query("profile_type=='prestorm'").reset_index('profile_type').z - poststorm_z = df_zone.query("profile_type=='poststorm'").reset_index('profile_type').z + prestorm_z = df_zone.query("profile_type=='prestorm'").reset_index("profile_type").z + poststorm_z = df_zone.query("profile_type=='poststorm'").reset_index("profile_type").z diff_z = prestorm_z - poststorm_z diff_z[abs(diff_z) < 0.2] = 0 - diff_vol = beach_volume( - x=diff_z.index.get_level_values("x"), - z=diff_z, - x_min=x_min, - x_max=x_max, - ) + diff_vol = beach_volume(x=diff_z.index.get_level_values("x"), z=diff_z, x_min=x_min, x_max=x_max) df_vol_changes.loc[site_id, "prestorm_{}_vol".format(zone)] = prestorm_vol df_vol_changes.loc[site_id, "poststorm_{}_vol".format(zone)] = poststorm_vol @@ -146,14 +141,13 @@ def storm_regime(df_observed_impacts): logger.info("Getting observed storm regimes") swash = (df_observed_impacts.dune_face_pct_change <= 2) & (df_observed_impacts.dune_face_vol_change <= 3) - collision = (df_observed_impacts.dune_face_pct_change >= 2) |(df_observed_impacts.dune_face_vol_change > 3) + collision = (df_observed_impacts.dune_face_pct_change >= 2) | (df_observed_impacts.dune_face_vol_change > 3) df_observed_impacts.loc[swash, "storm_regime"] = "swash" df_observed_impacts.loc[collision, "storm_regime"] = "collision" return df_observed_impacts - @click.command() @click.option("--profiles-csv", required=True, help="") @click.option("--profile-features-csv", required=True, help="") @@ -163,10 +157,10 @@ def create_observed_impacts(profiles_csv, profile_features_csv, output_file): logger.info("Creating observed wave impacts") logger.info("Importing data") df_profiles = pd.read_csv(profiles_csv, index_col=[0, 1, 2]) - df_profile_features = pd.read_csv(profile_features_csv, index_col=[0,1]) + df_profile_features = pd.read_csv(profile_features_csv, index_col=[0, 1]) logger.info("Creating new dataframe for observed impacts") - df_observed_impacts = pd.DataFrame(index=df_profile_features.index.get_level_values('site_id').unique()) + df_observed_impacts = pd.DataFrame(index=df_profile_features.index.get_level_values("site_id").unique()) logger.info("Getting pre/post storm volumes") df_swash_vol_changes = volume_change(df_profiles, df_profile_features, zone="swash") @@ -177,7 +171,7 @@ def create_observed_impacts(profiles_csv, profile_features_csv, output_file): df_observed_impacts = storm_regime(df_observed_impacts) # Save dataframe to csv - df_observed_impacts.to_csv(output_file, float_format='%.4f') + df_observed_impacts.to_csv(output_file, float_format="%.4f") logger.info("Saved to %s", output_file) logger.info("Done!") diff --git a/src/cli.py b/src/cli.py index 10b79e4..9b1d3e5 100644 --- a/src/cli.py +++ b/src/cli.py @@ -13,7 +13,9 @@ import analysis.observed_storm_impacts as observed_storm_impacts # Disable numpy warnings import warnings -warnings.simplefilter(action='ignore', category=FutureWarning) + +warnings.simplefilter(action="ignore", category=FutureWarning) + @click.group() def cli(): diff --git a/src/data/parse_mat.py b/src/data/parse_mat.py index e7afd84..b2e73c6 100644 --- a/src/data/parse_mat.py +++ b/src/data/parse_mat.py @@ -58,35 +58,41 @@ def parse_dune_crest_toes(df_sites, crest_mat, toe_mat): crest_data = loadmat(crest_mat) toe_data = loadmat(toe_mat) - for n, _ in enumerate(crest_data['xc1']): - rows.extend([{ - 'dune_crest_x': crest_data['xc1'][n], - 'dune_crest_z': crest_data['zc1'][n], - 'dune_toe_x': toe_data['xt1'][n], - 'dune_toe_z': toe_data['zt1'][n], - 'profile_type': 'prestorm', - 'site_no': n+1 - },{ - 'dune_crest_x': crest_data['xc2'][n], - 'dune_crest_z': crest_data['zc2'][n], - 'dune_toe_x': toe_data['xt2'][n], - 'dune_toe_z': toe_data['zt2'][n], - 'profile_type': 'poststorm', - 'site_no': n + 1 - }]) + for n, _ in enumerate(crest_data["xc1"]): + rows.extend( + [ + { + "dune_crest_x": crest_data["xc1"][n], + "dune_crest_z": crest_data["zc1"][n], + "dune_toe_x": toe_data["xt1"][n], + "dune_toe_z": toe_data["zt1"][n], + "profile_type": "prestorm", + "site_no": n + 1, + }, + { + "dune_crest_x": crest_data["xc2"][n], + "dune_crest_z": crest_data["zc2"][n], + "dune_toe_x": toe_data["xt2"][n], + "dune_toe_z": toe_data["zt2"][n], + "profile_type": "poststorm", + "site_no": n + 1, + }, + ] + ) df_profile_features = pd.DataFrame(rows) # Want the site_id instead of the site_no, so merge in df_sites df_sites.reset_index(inplace=True) - df_profile_features = df_sites[['site_no','site_id']].merge(df_profile_features, how='outer', on=['site_no']) - df_profile_features.drop(columns=['site_no'],inplace=True) - df_profile_features.set_index(['site_id','profile_type'], inplace=True) + df_profile_features = df_sites[["site_no", "site_id"]].merge(df_profile_features, how="outer", on=["site_no"]) + df_profile_features.drop(columns=["site_no"], inplace=True) + df_profile_features.set_index(["site_id", "profile_type"], inplace=True) df_profile_features.sort_index(inplace=True) df_profile_features = df_profile_features.round(3) return df_profile_features + def combine_sites_and_orientaions(df_sites, df_orientations): """ Replaces beach/lat/lon columns with the unique site_id. @@ -193,7 +199,7 @@ def parse_profiles_and_sites(profiles_mat): site_counter = 0 for i, site in enumerate(mat_data["site"]): - logger.debug('Processing site {} of {}'.format(i+1, len(mat_data['site']))) + logger.debug("Processing site {} of {}".format(i + 1, len(mat_data["site"]))) # Give each site a unique id if len(site_rows) == 0 or site_rows[-1]["beach"] != site: @@ -248,7 +254,6 @@ def parse_profiles_and_sites(profiles_mat): profile_rows.append( { "site_id": site_id, - "lon": lon[0], "lat": lat[0], "profile_type": profile_type, @@ -387,6 +392,7 @@ def create_profile_features(crest_mat, toe_mat, sites_csv, output_file): df_profile_features.to_csv(output_file) logger.info("Created %s", output_file) + @click.command(short_help="create profiles.csv") @click.option("--profiles-mat", required=True, help=".mat file containing beach profiles") @click.option("--profiles-output-file", required=True, help="where to save profiles.csv")