diff --git a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
index 15ff9be2d7022cc6f0b13f1650bd7d24169ea7e5..d46c8ffd58f5cd873097f3108485d21c54899925 100644
--- a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
+++ b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
@@ -413,9 +413,15 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
b[2].append(diff)
b = [np.mean(np.array(e)) for e in b]
- # Return the concatenated results
- concatenated_result = list(a) + list(b)
- return concatenated_result
+ # Third array for curve of percentage of exceedance for return levels
+ c = [[] , [], []]
+ for _, eurocode, uncertainty in triplet:
+ metrics = [uncertainty.confidence_interval[0], uncertainty.mean_estimate, uncertainty.confidence_interval[1]]
+ for j, metric in enumerate(metrics):
+ relative_difference = 100 * (metric - eurocode) / eurocode
+ c[j].append(relative_difference)
+ c = np.array([np.mean(np.array(e)) if len(e) > 0 else 0 for e in c])
+ return a, b, c
# Part 3 - QQPLOT
diff --git a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
index 9053d8ce360c89bdfcdeeaa5db42239bc3e31514..6205582112fca1ac1d69981216b63729038bff0d 100644
--- a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
+++ b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
@@ -72,12 +72,12 @@ def intermediate_result(altitudes, massif_names=None,
_ = compute_minimized_aic(visualizer)
# Plots
- plot_trend_map(altitude_to_visualizer)
- plot_trend_curves(altitude_to_visualizer={a: v for a, v in altitude_to_visualizer.items() if a >= 900})
+ # plot_trend_map(altitude_to_visualizer)
+ # plot_trend_curves(altitude_to_visualizer={a: v for a, v in altitude_to_visualizer.items() if a >= 900})
plot_uncertainty_massifs(altitude_to_visualizer)
plot_uncertainty_histogram(altitude_to_visualizer)
- plot_selection_curves(altitude_to_visualizer)
- plot_intensity_against_gumbel_quantile_for_3_examples(altitude_to_visualizer)
+ # plot_selection_curves(altitude_to_visualizer)
+ # plot_intensity_against_gumbel_quantile_for_3_examples(altitude_to_visualizer)
# Additional plots
# uncertainty_interval_size(altitude_to_visualizer)
@@ -93,7 +93,7 @@ def major_result():
# study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1]
altitudes = [300, 600, 900, 1800, 2700][:2]
altitudes = [300, 600, 900, 1200, 1500, 1800]
- altitudes = paper_altitudes
+ # altitudes = paper_altitudes
# altitudes = [900, 1800, 270{{0][:1]
for study_class in study_classes:
print('new stuy class', study_class)
diff --git a/projects/exceeding_snow_loads/section_results/plot_diagnosis_risk.py b/projects/exceeding_snow_loads/section_results/plot_diagnosis_risk.py
index ea587bd2079e9a4d6efdf2ab561bcdc3e89bc3d8..f15a554d18d3e3673550e8d6f65b90d482fea040 100644
--- a/projects/exceeding_snow_loads/section_results/plot_diagnosis_risk.py
+++ b/projects/exceeding_snow_loads/section_results/plot_diagnosis_risk.py
@@ -23,7 +23,7 @@ def plot_diagnosis_risk(altitude_to_visualizer):
def plot_mean_exceedance(visualizers, altitudes, ax, model_subset_for_uncertainty, ci_method):
- l = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[3:] for v in visualizers]
+ l = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[1] for v in visualizers]
diff, diff_c, diff_e = zip(*l)
ax.plot(altitudes, diff, marker='o', color='red', label='diff')
ax.plot(altitudes, diff_c, marker='o', color='yellow', label='diff c')
diff --git a/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py b/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py
index cebdc81c014259ebbf1d1a551aad492fa14880a4..74767c3313784fd8e87a00b3b693e897f6fccbca 100644
--- a/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py
+++ b/projects/exceeding_snow_loads/section_results/plot_uncertainty_curves.py
@@ -122,10 +122,10 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, m
ax.set_xlim([200, 1900])
if massif_name in ['Maurienne', 'Chartreuse', 'Beaufortain']:
ax.set_ylim([-1.5, 13.5])
- ax.set_yticks([2 * i for i in range(7)])
+ ax.set_yticks([-1] + [2 * i for i in range(7)])
if massif_name in ['Vercors']:
- ax.set_ylim([-1, 10])
- ax.set_yticks([2 * i for i in range(6)])
+ ax.set_ylim([-1.5, 10.5])
+ ax.set_yticks([-1] + [2 * i for i in range(6)])
diff --git a/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py b/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py
index d6d31e60540a542c85467a039425b5f9931ac9de..923e9211fc30cf6b5cb44fef127d1dae093d9fbc 100644
--- a/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py
+++ b/projects/exceeding_snow_loads/section_results/plot_uncertainty_histogram.py
@@ -3,6 +3,8 @@ import matplotlib.pyplot as plt
import numpy as np
from extreme_data.eurocode_data.utils import EUROCODE_ALTITUDES
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
+ AbstractExtractEurocodeReturnLevel
from projects.exceeding_snow_loads.utils import dpi_paper1_figure
from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
@@ -44,7 +46,13 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
width = 100
else:
width = 200
- plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width=width)
+ legend_size = 10
+ # Plot histogram on the left axis
+ plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width, legend_size)
+
+ # Plot percentages of return level excess on the right axis
+ plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax.twinx(), fontsize_label, legend_size)
+
ax.set_xticks(altitudes)
ax.tick_params(labelsize=fontsize_label)
@@ -55,7 +63,7 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
# 'exceeds French standards (\%)', fontsize=fontsize_label)
ax.set_ylabel('Massifs exceeding French standards (\%)', fontsize=fontsize_label)
ax.set_xlabel('Altitude (m)', fontsize=fontsize_label)
- ax.set_ylim([0, 100])
+ ax.set_ylim([0, 110])
ax.yaxis.grid()
ax_twiny = ax.twiny()
@@ -66,20 +74,19 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
nb_massif_names = [len(v.intersection_of_massif_names_fitted) for v in altitude_to_visualizer.values()]
ax_twiny.set_xticklabels(nb_massif_names)
- ax_twiny.set_xlabel('Total number of massifs at each altitude (for the percentage)', fontsize=fontsize_label)
+ ax_twiny.set_xlabel('Number of massifs at each altitude (for the percentage and the mean)', fontsize=fontsize_label)
- ax.set_yticks([10 * i for i in range(11)])
+ ax.set_yticks([20 * i for i in range(6)])
visualizer.plot_name = 'Percentages of exceedance with {}'.format(
get_display_name_from_object_type(model_subset_for_uncertainty))
- # visualizer.show = True
- visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
+ visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure, tight_layout=True)
ax.clear()
ax_twiny.clear()
plt.close()
-def plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width):
- three_percentages_of_excess = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[:3] for v in
+def plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_method, ax, bincenters, width, legend_size=10):
+ three_percentages_of_excess = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[0] 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
@@ -90,4 +97,47 @@ def plot_histogram_ci_method(visualizers, model_subset_for_uncertainty, ci_metho
yerr = np.array([[d[1] - d[0], d[2] - d[1]] for d in three_percentages_of_excess]).transpose()
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)
+ ecolor = 'black'
+ label_name = 'Percentage of massifs exceeding'
+ # ax.bar(bincenters, y, width=width, color=color, yerr=yerr, label=label_name, ecolor=ecolor, capsize=5)
+ ax.bar(bincenters, y, width=width, color=color, label=label_name)
+ # Just to add something in the legend
+ label_confidence_interval = get_label_confidence_interval(label_name)
+ ax.errorbar(bincenters, y, yerr=yerr, label=label_confidence_interval,
+ fmt='none', color=ecolor, capsize=5)
+
+ ax.legend(loc='upper left', prop={'size': legend_size})
+
+
+def plot_percentage_of_excess(visualizers, model_subset_for_uncertainty, ci_method, ax, fontsize_label, legend_size=10):
+ l = [v.excess_metrics(ci_method, model_subset_for_uncertainty)[2] for v in visualizers]
+ lower_bound, mean, upper_bound = list(zip(*l))
+ other_mean = [e[1] for e in l]
+ altitudes = [v.altitude for v in visualizers]
+ # Display parameters
+ color = 'blue'
+ alpha = 0.2
+ full_label_name = 'Mean relative difference between\n' \
+ '50-year return levels and French standards (\%)'
+
+ label_name = 'Mean relative difference'
+ print('mean relative difference', mean)
+ ax.plot(altitudes, mean, linestyle='--', marker='o', color=color,
+ label=label_name)
+
+ label_confidence_interval = get_label_confidence_interval(label_name)
+ ax.fill_between(altitudes, lower_bound, upper_bound, color=color, alpha=alpha,
+ label=label_confidence_interval)
+
+ ax.tick_params(labelsize=fontsize_label)
+ ax.legend(loc='upper right', prop={'size': legend_size})
+ ax.set_ylabel(full_label_name, fontsize=fontsize_label)
+ ax.set_ylim([0, 110])
+
+
+def get_label_confidence_interval(label_name):
+ confidence_interval_str = '\n{}'.format(AbstractExtractEurocodeReturnLevel.percentage_confidence_interval)
+ confidence_interval_str += ' \% confidence interval'
+ label_confidence_interval = label_name + confidence_interval_str
+ return label_confidence_interval
+