diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py index cdc112819713da0768c95b215757382b789160de..d4b20888120afe71359a9ace1749de4908ac963c 100644 --- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py +++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py @@ -30,7 +30,7 @@ def main(): massif_names = None altitudes_list = altitudes_for_groups[:2] elif fast: - massif_names = ['Mercantour', 'Vercors', 'Ubaye'] + massif_names = ['Vanoise', 'Haute-Tarentaise', 'Vercors'] altitudes_list = altitudes_for_groups[:2] else: massif_names = None @@ -54,7 +54,7 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes): def plot_visualizers(massif_names, visualizer_list): - # plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list) + plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list) # plot_histogram_all_models_against_altitudes(massif_names, visualizer_list) plot_coherence_curves(massif_names, visualizer_list) pass diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py index cd7e8d556284b2e7ee8d41091020f7ec3cf0c1f0..5fdb9258e85c18f5c0641d398edde73e4ca29213 100644 --- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py +++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py @@ -395,11 +395,11 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): nb_valid_massif_names = len(valid_massif_names) nbs = np.zeros(4) - relative_changes = [] + changes = [] for one_fold in [one_fold for m, one_fold in self.massif_name_to_one_fold_fit.items() if m in valid_massif_names]: - # Compute relative changes - relative_changes.append(one_fold.relative_change_in_return_level_for_reference_altitude) + # Compute changes + changes.append(one_fold.change_in_return_level_for_reference_altitude) # Compute nb of non stationary models if one_fold.change_in_return_level_for_reference_altitude == 0: continue @@ -410,10 +410,9 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): nbs[idx + 1] += 1 percents = 100 * nbs / nb_valid_massif_names - mean_relative_change = np.mean(relative_changes) - median_relative_change = np.median(relative_changes) - print('mean', mean_relative_change, relative_changes) - return [nb_valid_massif_names, mean_relative_change, median_relative_change] + list(percents) + mean_change = np.mean(changes) + median_change = np.median(changes) + return [nb_valid_massif_names, mean_change, median_change] + list(percents) def get_valid_names(self, massif_names): valid_massif_names = set(self.massif_name_to_one_fold_fit.keys()) diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py index c9b590607ca81282674966fe8005be5a93a2cd8d..a0db58c0959f55ef88405629194c288948e34ee6 100644 --- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py +++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py @@ -312,8 +312,8 @@ class OneFoldFit(object): standard_gumbel_quantiles = [standard_gumbel_distribution.quantile(i / (n + 1)) for i in range(1, n + 1)] return standard_gumbel_quantiles - def best_confidence_interval(self, altitude) -> EurocodeConfidenceIntervalFromExtremes: + def best_confidence_interval(self, altitude, year) -> EurocodeConfidenceIntervalFromExtremes: EurocodeConfidenceIntervalFromExtremes.quantile_level = 1 - (1 / self.return_period) return EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator_extremes=self.best_estimator, ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle, - coordinate=np.array([altitude, 2019])) + coordinate=np.array([altitude, year])) diff --git a/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py b/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py index 5298b24043b2b87fce570dfc24eacc7cdb479dd7..71ad064db4f7fa3d5192f6a58f74a3bf31a5b88a 100644 --- a/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py +++ b/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py @@ -16,19 +16,20 @@ def plot_coherence_curves(massif_names, visualizer_list: List[ for massif_name in all_valid_names: _, axes = plt.subplots(2, 2) axes = axes.flatten() - colors = ['blue', 'green'] - labels = ['Altitudinal model', 'Pointwise model'] - altitudinal_model = [True, False] - for color, global_label, boolean in list(zip(colors, labels, altitudinal_model))[:]: - plot_coherence_curve(axes, massif_name, visualizer_list, boolean, color, global_label) + colors = ['blue', 'yellow', 'green'] + labels = ['Altitudinal-temporal model in 2019', 'Altitudinal-temporal model in 1969', 'Stationary model'] + altitudinal_model = [True, True, False] + years = [2019, 1969, None] + for color, global_label, boolean, year in list(zip(colors, labels, altitudinal_model, years))[::2]: + plot_coherence_curve(axes, massif_name, visualizer_list, boolean, color, global_label, year) visualizer.plot_name = '{}/{}'.format(folder, massif_name.replace('_', '-')) - visualizer.show_or_save_to_file(add_classic_title=False, no_title=True) + visualizer.show_or_save_to_file(add_classic_title=False, no_title=True, dpi=200) plt.close() def plot_coherence_curve(axes, massif_name, visualizer_list: List[AltitudesStudiesVisualizerForNonStationaryModels], - is_altitudinal, color, global_label): - x_all_list, values_all_list, labels, all_bound_list = load_all_list(massif_name, visualizer_list, is_altitudinal) + is_altitudinal, color, global_label, year): + x_all_list, values_all_list, labels, all_bound_list = load_all_list(massif_name, visualizer_list, is_altitudinal, year) for i, label in enumerate(labels): ax = axes[i] # Plot with complete line @@ -58,12 +59,12 @@ def plot_coherence_curve(axes, massif_name, visualizer_list: List[AltitudesStudi ax.set_ylabel(label) - ax.legend(prop={'size': 10}) + ax.legend(prop={'size': 5}) if i >= 2: ax.set_xlabel('Altitude') -def load_all_list(massif_name, visualizer_list, altitudinal_model=True): +def load_all_list(massif_name, visualizer_list, altitudinal_model=True, year=2019): all_altitudes_list = [] all_values_list = [] all_bound_list = [] @@ -74,8 +75,8 @@ def load_all_list(massif_name, visualizer_list, altitudinal_model=True): min_altitude, *_, max_altitude = visualizer.massif_name_to_massif_altitudes[massif_name] one_fold_fit = visualizer.massif_name_to_one_fold_fit[massif_name] altitudes_list = list(range(min_altitude, max_altitude, 10)) - gev_params_list = [one_fold_fit.get_gev_params(altitude, 2019) for altitude in altitudes_list] - confidence_interval_values = [one_fold_fit.best_confidence_interval(altitude) for altitude in altitudes_list] + gev_params_list = [one_fold_fit.get_gev_params(altitude, year) for altitude in altitudes_list] + confidence_interval_values = [one_fold_fit.best_confidence_interval(altitude, year) for altitude in altitudes_list] else: assert OneFoldFit.return_period == 100, 'change the call below' altitudes_list, study_list_valid = zip(*[(a, s) for a, s in visualizer.studies.altitude_to_study.items() @@ -98,7 +99,6 @@ def load_all_list(massif_name, visualizer_list, altitudinal_model=True): all_altitudes_list.append(altitudes_list) all_bound_list.append(list(zip(*bound_list))) labels = ['location parameter', 'scale parameter', 'shape pameter', '{}-year return level'.format(OneFoldFit.return_period)] - labels = [l + ' in 2019' for l in labels] return all_altitudes_list, all_values_list, labels, all_bound_list diff --git a/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py b/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py index ffd4860e4d05f3f08f2750b782998e32efb33fd4..b524e119e3599a9e24be614b773e966048e6cba4 100644 --- a/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py +++ b/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py @@ -69,7 +69,7 @@ def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list: L visualizer = visualizer_list[0] all_trends = [v.all_trends(massif_names) for v in visualizer_list] - nb_massifs, mean_relative_changes, median_relative_changes, *all_l = zip(*all_trends) + nb_massifs, mean_changes, median_changes, *all_l = zip(*all_trends) plt.close() ax = plt.gca() @@ -101,14 +101,19 @@ def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list: L # PLot mean relative change against altitude ax_twinx = ax.twinx() - ax_twinx.plot(x, mean_relative_changes, label='Mean relative change', color='black') - # ax_twinx.plot(x, median_relative_changes, label='Median relative change', color='grey') + ax_twinx.plot(x, mean_changes, label='Mean change', color='black') + # ax_twinx.plot(x, median_changes, label='Median change', color='grey') ax_twinx.legend(loc='upper right', prop={'size': size}) - ylabel = 'Relative change of {}-year return levels ({})'.format(OneFoldFit.return_period, - visualizer.study.variable_unit) + ylabel = 'Change of {}-year return levels ({})'.format(OneFoldFit.return_period, visualizer.study.variable_unit) ax_twinx.set_ylabel(ylabel, fontsize=legend_fontsize) + low, up = ax_twinx.get_ylim() + low2, up2 = ax.get_ylim() + new_lim = min(low, low2), max(up, up2) + ax.set_ylim(new_lim) + ax_twinx.set_ylim(new_lim) + visualizer.plot_name = 'All trends' visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)