From 5b9b4141b33d7687e863bbc698b782240d1ec152 Mon Sep 17 00:00:00 2001 From: Chris Leaman Date: Tue, 4 Dec 2018 09:29:51 +1100 Subject: [PATCH] Fix bugs for forecasting impacts --- src/analysis/forecasted_storm_impacts.py | 6 +- src/analysis/observed_storm_impacts.py | 72 ++++++++++++------------ 2 files changed, 39 insertions(+), 39 deletions(-) diff --git a/src/analysis/forecasted_storm_impacts.py b/src/analysis/forecasted_storm_impacts.py index 28db35a..0ea92ec 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) + 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,7 +29,7 @@ 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[["dune_toe_z", "dune_crest_z"]], how="left", left_index=True, right_index=True + df_profile_features.query("profile_type=='prestorm'")[["dune_toe_z", "dune_crest_z"]], how="left", left_index=True, right_index=True ) # Compare R_high and R_low wirth dune crest and toe elevations @@ -73,7 +73,7 @@ def storm_regime(df_forecasted_impacts): 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]) + 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) diff --git a/src/analysis/observed_storm_impacts.py b/src/analysis/observed_storm_impacts.py index 6900b31..6c76e36 100644 --- a/src/analysis/observed_storm_impacts.py +++ b/src/analysis/observed_storm_impacts.py @@ -30,14 +30,15 @@ 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) + 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)) - prestorm_dune_toe_x = df_profile_features.loc[df_profile_features.index == site_id].dune_toe_x.tolist() - prestorm_dune_crest_x = df_profile_features.loc[df_profile_features.index == site_id].dune_crest_x.tolist() + 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() # We may not have a dune toe or crest defined, or there may be multiple defined. prestorm_dune_crest_x = return_first_or_nan(prestorm_dune_crest_x) @@ -61,14 +62,20 @@ def volume_change(df_profiles, df_profile_features, zone): for profile_type in ["prestorm", "poststorm"] ] ) + x_first_obs = max( + [ + min(df_zone.query("profile_type == '{}'".format(profile_type)).index.get_level_values("x")) + for profile_type in ["prestorm", "poststorm"] + ] + ) # Where we want to measure pre and post storm volume is dependant on the zone selected if zone == "swash": - x_min = prestorm_dune_toe_x + x_min = max(prestorm_dune_toe_x,x_first_obs) x_max = x_last_obs elif zone == "dune_face": - x_min = prestorm_dune_crest_x - x_max = prestorm_dune_toe_x + x_min = max(prestorm_dune_crest_x, x_first_obs) + x_max = min(prestorm_dune_toe_x,x_last_obs) else: logger.warning("Zone argument not properly specified. Please check") x_min = None @@ -89,13 +96,23 @@ def volume_change(df_profiles, df_profile_features, zone): x_max=x_max, ) - # No point keeping so many decimal places, let's round them - prestorm_vol = np.round(prestorm_vol, 2) - poststorm_vol = np.round(poststorm_vol, 2) + # 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 + 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, + ) df_vol_changes.loc[site_id, "prestorm_{}_vol".format(zone)] = prestorm_vol df_vol_changes.loc[site_id, "poststorm_{}_vol".format(zone)] = poststorm_vol - df_vol_changes.loc[site_id, "{}_vol_change".format(zone)] = prestorm_vol - poststorm_vol + df_vol_changes.loc[site_id, "{}_vol_change".format(zone)] = diff_vol + df_vol_changes.loc[site_id, "{}_pct_change".format(zone)] = diff_vol / prestorm_vol * 100 return df_vol_changes @@ -127,31 +144,14 @@ def storm_regime(df_observed_impacts): :return: """ logger.info("Getting observed storm regimes") - df_observed_impacts.loc[df_observed_impacts.dune_face_vol_change <= 5, "storm_regime"] = "swash" - df_observed_impacts.loc[df_observed_impacts.dune_face_vol_change > 5, "storm_regime"] = "collision" - return df_observed_impacts + 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) + + df_observed_impacts.loc[swash, "storm_regime"] = "swash" + df_observed_impacts.loc[collision, "storm_regime"] = "collision" + return df_observed_impacts -# -# if __name__ == "__main__": -# logger.info("Importing existing data") -# data_folder = "./data/interim" -# df_profiles = pd.read_csv(os.path.join(data_folder, "profiles.csv"), index_col=[0, 1, 2]) -# df_profile_features = pd.read_csv(os.path.join(data_folder, "profile_features.csv"), index_col=[0]) -# -# logger.info("Creating new dataframe for observed impacts") -# df_observed_impacts = pd.DataFrame(index=df_profile_features.index) -# -# 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(df_profiles, df_profile_features, zone="dune_face") -# df_observed_impacts = df_observed_impacts.join([df_swash_vol_changes, df_dune_face_vol_changes]) -# -# # Classify regime based on volume changes -# df_observed_impacts = storm_regime(df_observed_impacts) -# -# # Save dataframe to csv -# df_observed_impacts.to_csv(os.path.join(data_folder, "impacts_observed.csv")) @click.command() @@ -163,10 +163,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]) + 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) + 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 +177,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) + df_observed_impacts.to_csv(output_file, float_format='%.4f') logger.info("Saved to %s", output_file) logger.info("Done!")