From ab3d7d0cf9f0a30081cd50283277fdee547064e5 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Fri, 2 Oct 2020 18:55:07 +0200
Subject: [PATCH] [contrasting] improve coherence plot & trend histogram
---
.../altitudes_fit/main_altitudes_studies.py | 4 +--
...es_visualizer_for_non_stationary_models.py | 13 +++++-----
.../one_fold_analysis/one_fold_fit.py | 4 +--
.../plots/plot_coherence_curves.py | 26 +++++++++----------
.../plots/plot_histogram_altitude_studies.py | 15 +++++++----
5 files changed, 33 insertions(+), 29 deletions(-)
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 cdc11281..d4b20888 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 cd7e8d55..5fdb9258 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 c9b59060..a0db58c0 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 5298b240..71ad064d 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 ffd4860e..b524e119 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)
--
GitLab