From c6d9303f14d0843b1b480c2d89d90956353bdce4 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 2 Dec 2019 14:43:48 +0100
Subject: [PATCH] [PAPER 1] improve all uncertainty plots
---
experiment/eurocode_data/eurocode_region.py | 4 ++-
.../scm_models_data/abstract_study.py | 22 +++++++++-------
.../main_study_visualizer.py | 2 +-
.../main_result_trends_and_return_levels.py | 14 ++++++-----
.../plot_uncertainty_curves.py | 25 ++++++++++++++-----
.../plot_uncertainty_histogram.py | 21 ++++++++++------
.../confidence_interval_method.py | 15 +++++++++--
7 files changed, 70 insertions(+), 33 deletions(-)
diff --git a/experiment/eurocode_data/eurocode_region.py b/experiment/eurocode_data/eurocode_region.py
index 5fc80319..fb4ca7d9 100644
--- a/experiment/eurocode_data/eurocode_region.py
+++ b/experiment/eurocode_data/eurocode_region.py
@@ -44,7 +44,9 @@ class AbstractEurocodeRegion(object):
def lois_de_variation_1000_and_2000(self):
return 3.5, -2.45
- def plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(self, ax, altitudes, label='Eurocode standards', linestyle=None):
+ def plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(self, ax, altitudes,
+ label='French standards',
+ linestyle=None):
ax.plot(altitudes, [self.valeur_caracteristique(altitude) for altitude in altitudes],
label=label, color=self.eurocode_color, linewidth=5, linestyle=linestyle)
diff --git a/experiment/meteo_france_data/scm_models_data/abstract_study.py b/experiment/meteo_france_data/scm_models_data/abstract_study.py
index 173cfa44..23f39bc5 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -401,22 +401,26 @@ class AbstractStudy(object):
# Add legend for the marker
if massif_name_to_marker_style is not None:
- labels = ['\mathcal{M}_{\mu_1}', '\mathcal{M}_{\sigma_1}', '\mathcal{M}_{\mu_1, \sigma_1}']
- markers = ["s", "^", "D"]
- legend_elements = [
- Line2D([0], [0], marker=marker, color='w', label='${}$'.format(label),
- markerfacecolor='w', markeredgecolor='k', markersize=8)
- for label, marker in zip(labels, markers)
- ]
+ legend_elements = cls.get_legend_for_model_symbol(markersize=8)
ax.legend(handles=legend_elements, bbox_to_anchor=(0.01, 0.03), loc='lower left')
- ax.annotate("Filled symbol = significant trend ", xy=(0.05, 0.015), xycoords='axes fraction',
- fontsize=7)
+ ax.annotate("Filled symbol = significant trend ", xy=(0.05, 0.015), xycoords='axes fraction', fontsize=7)
if show:
plt.show()
return ax
+ @classmethod
+ def get_legend_for_model_symbol(cls, markersize):
+ labels = ['\mathcal{M}_{\mu_1}', '\mathcal{M}_{\sigma_1}', '\mathcal{M}_{\mu_1, \sigma_1}']
+ markers = ["s", "^", "D"]
+ legend_elements = [
+ Line2D([0], [0], marker=marker, color='w', label='${}$'.format(label),
+ markerfacecolor='w', markeredgecolor='k', markersize=markersize)
+ for label, marker in zip(labels, markers)
+ ]
+ return legend_elements
+
"""
CLASS ATTRIBUTES COMMON TO ALL OBJECTS
(written as object attributes/methods for simplicity)
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py
index 67ac17dd..0e32a1e6 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py
@@ -34,7 +34,7 @@ SCM_STUDY_CLASS_TO_ABBREVIATION = {
SafranSnowfall: 'SF3',
CrocusSweTotal: 'SWE',
CrocusSwe3Days: 'SWE3',
- CrocusSnowLoadEurocode: 'GSL_Eurocode',
+ CrocusSnowLoadEurocode: 'GSL from annual maximum of HS and {}'.format(eurocode_snow_density),
CrocusDepth: 'SD',
CrocusSnowLoadTotal: 'max GSL',
CrocusSnowLoad3Days: 'GSL3',
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
index 3c6df869..d472fa8e 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
@@ -7,7 +7,7 @@ from experiment.paper_past_snow_loads.paper_utils import paper_study_classes, pa
load_altitude_to_visualizer
from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_curves import \
plot_uncertainty_massifs
-from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoadEurocode
from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_histogram import \
plot_uncertainty_histogram
from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
@@ -58,15 +58,17 @@ def major_result():
ConfidenceIntervalMethodFromExtremes.ci_mle][:]
massif_names = None
non_stationary_uncertainty = [False, True][:]
- for study_class in paper_study_classes[2:]:
+ for study_class in paper_study_classes[:2]:
intermediate_result(paper_altitudes, massif_names, non_stationary_uncertainty, uncertainty_methods, study_class)
if __name__ == '__main__':
- # major_result()
- intermediate_result(altitudes=[900, 1200], massif_names=['Chartreuse', 'Vanoise', 'Mercantour'],
- uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
- ConfidenceIntervalMethodFromExtremes.ci_mle])
+ major_result()
+ # intermediate_result(altitudes=paper_altitudes, massif_names=['Chartreuse', 'Vanoise', 'Mercantour'],
+ # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+ # ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+ # non_stationary_uncertainty=[False, True][1:],
+ # study_class=CrocusSnowLoadEurocode)
# intermediate_result(altitudes=[900, 1200], massif_names=None)
# intermediate_result(ALL_ALTITUDES_WITHOUT_NAN)
# intermediate_result(paper_altitudes)
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
index bd9723a1..1192911f 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
@@ -4,6 +4,9 @@ import matplotlib.pyplot as plt
import numpy as np
from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR, EUROCODE_ALTITUDES
+from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
+ SCM_STUDY_CLASS_TO_ABBREVIATION
from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
@@ -102,20 +105,20 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes)
ax.legend(loc=2)
- ax.set_ylim([-1, 16])
+ # ax.set_ylim([-1, 16])
massif_name_str = massif_name.replace('_', ' ')
eurocode_region_str = get_display_name_from_object_type(type(eurocode_region))
is_non_stationary_model = non_stationary_context if isinstance(non_stationary_context,
bool) else 'Non' in non_stationary_context
if is_non_stationary_model:
- non_stationary_context = 'non-stationary'
+ non_stationary_context = 'selected non-stationary models'
else:
- non_stationary_context = 'stationary'
- title = '{} ({} Eurocodes area) with a {} model'.format(massif_name_str, eurocode_region_str,
- non_stationary_context)
+ non_stationary_context = 'the stationary model'
+ title = '{} massif with {}'.format(massif_name_str, non_stationary_context)
ax.set_title(title)
ax.set_xticks(altitudes)
- ax.set_ylabel(EUROCODE_RETURN_LEVEL_STR)
+ ylabel = EUROCODE_RETURN_LEVEL_STR.replace('GSL', SCM_STUDY_CLASS_TO_ABBREVIATION[type(visualizer.study)])
+ ax.set_ylabel(ylabel)
ax.set_xlabel('Altitude (m)')
ax.grid()
@@ -131,10 +134,20 @@ def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes):
edgecolor='black', hatch='//')
# Plot markers
markers_kwargs = [v.massif_name_to_marker_style[massif_name] for v in visualizers]
+ for k in markers_kwargs:
+ k['markersize'] = 7
for altitude, marker_kwargs, value in zip(valid_altitudes, markers_kwargs, tdrl_values):
# ax.plot([altitude], [value / 2], **marker_kwargs)
# Better to plot all the markers on the same line
ax.plot([altitude], 0, **marker_kwargs)
+ # Add a legend plot
+ legend_elements = AbstractStudy.get_legend_for_model_symbol(markersize=9)
+ ax2 = ax.twinx()
+ # ax2.legend(handles=legend_elements, bbox_to_anchor=(0.93, 0.7), loc='upper right')
+ # ax2.annotate("Filled symbol = significant trend ", xy=(0.85, 0.5), xycoords='axes fraction', fontsize=7)
+ ax2.legend(handles=legend_elements, loc='upper right')
+ ax2.annotate("Filled symbol = significant trend ", xy=(0.75, 0.97), xycoords='axes fraction', fontsize=10)
+ ax2.set_yticks([])
def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
index dbbfc221..468b538f 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
@@ -2,16 +2,18 @@ from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
-from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR
+from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR, EUROCODE_ALTITUDES
from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
-from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import ci_method_to_color
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import ci_method_to_color, \
+ ci_method_to_label
def plot_uncertainty_histogram(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
""" Plot one graph for each non-stationary context
:return:
"""
+ altitude_to_visualizer = {a: v for a, v in altitude_to_visualizer.items() if a in EUROCODE_ALTITUDES}
visualizer = list(altitude_to_visualizer.values())[0]
for non_stationary_context in visualizer.non_stationary_contexts:
plot_histogram(altitude_to_visualizer, non_stationary_context)
@@ -33,24 +35,27 @@ def plot_histogram(altitude_to_visualizer, non_stationary_context):
if len(visualizer.uncertainty_methods) == 2:
offset = -50 if j == 0 else 50
bincenters = altitudes + offset
-
width = 100
plot_histogram_ci_method(visualizers, non_stationary_context, ci_method, ax, bincenters, width=width)
ax.set_xticks(altitudes)
+ ax.legend()
ax.set_ylabel('Massifs exceeding French standards (\%)')
ax.set_xlabel('Altitude (m)')
ax.set_ylim([0, 100])
ax.set_yticks([10 * i for i in range(11)])
- # visualizer.plot_name = 'Percentage '
+ visualizer.plot_name = 'Percentages of exceedance with non_stationary={}'.format(non_stationary_context)
# visualizer.show = True
- # visualizer.show_or_save_to_file(no_title=True)
- plt.show()
+ visualizer.show_or_save_to_file(no_title=True)
def plot_histogram_ci_method(visualizers, non_stationary_context, ci_method, ax, bincenters, width):
three_percentages_of_excess = [v.three_percentages_of_excess(ci_method, non_stationary_context) for v in
visualizers]
+ epsilon = 0.5
+ three_percentages_of_excess = [(a, b, c) if a == b else (max(epsilon, a), b, c) for (a, b, c) in three_percentages_of_excess]
+ three_percentages_of_excess = [(a, b, c) if b == c else (a, b, min(100 - epsilon, c)) for (a, b, c) in three_percentages_of_excess]
y = [d[1] for d in three_percentages_of_excess]
yerr = np.array([[d[1] - d[0], d[2] - d[1]] for d in three_percentages_of_excess]).transpose()
- print(y, three_percentages_of_excess)
- ax.bar(bincenters, y, width=width, color=ci_method_to_color[ci_method], yerr=yerr)
+ label = ci_method_to_label[ci_method]
+ color = ci_method_to_color[ci_method]
+ ax.bar(bincenters, y, width=width, color=color, yerr=yerr, label=label, ecolor='black', capsize=5)
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
index e28d30fb..d46c9b5b 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
@@ -17,7 +17,18 @@ ci_method_to_method_name = {
ConfidenceIntervalMethodFromExtremes.ci_proflik: 'proflik',
}
+# ci_method_to_color = {
+# ConfidenceIntervalMethodFromExtremes.ci_mle: 'tab:brown',
+# ConfidenceIntervalMethodFromExtremes.my_bayes: 'tab:green'
+# }
+
ci_method_to_color = {
- ConfidenceIntervalMethodFromExtremes.ci_mle: 'tab:brown',
- ConfidenceIntervalMethodFromExtremes.my_bayes: 'tab:green'
+ ConfidenceIntervalMethodFromExtremes.ci_mle: 'yellowgreen',
+ ConfidenceIntervalMethodFromExtremes.my_bayes: 'darkgreen'
+}
+
+common_part_and_uncertainty = 'Return level and its uncertainty with the'
+ci_method_to_label = {
+ ConfidenceIntervalMethodFromExtremes.ci_mle: 'Bayesian procedure',
+ ConfidenceIntervalMethodFromExtremes.my_bayes: 'Delta method '
}
\ No newline at end of file
--
GitLab