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 abe1fbef39ba6c3d943374b6b7db22b18da4169b..1d20283ef685b1cee0e93802205c8480f5ee4e26 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -259,6 +259,10 @@ class AbstractStudy(object):
# Massif names that are present in the current study (i.e. for the current altitude)
return self.altitude_to_massif_names[self.altitude]
+ @property
+ def nb_study_massif_names(self) -> int:
+ return len(self.study_massif_names)
+
@property
def df_massifs_longitude_and_latitude(self) -> pd.DataFrame:
# DataFrame object that represents the massif coordinates in degrees extracted from the SCM data
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 d0d7be1a2b899d1d2e3e4cf91a288fc68e12a470..ca4f924bcae73705848abee13895936b144eb76a 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
@@ -4,8 +4,11 @@ import matplotlib as mpl
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoadEurocode, \
CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
+ ALL_ALTITUDES_WITHOUT_NAN
from experiment.paper_past_snow_loads.paper_main_utils import load_altitude_to_visualizer
from experiment.paper_past_snow_loads.paper_utils import paper_study_classes, paper_altitudes, ModelSubsetForUncertainty
+from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_trend_curves import plot_trend_curves
from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_curves import \
plot_uncertainty_massifs
from experiment.paper_past_snow_loads.result_trends_and_return_levels.plot_uncertainty_histogram import \
@@ -75,10 +78,13 @@ def intermediate_result(altitudes, massif_names=None,
else:
visualizer.plot_trends(None, add_colorbar=True)
+ # Plot trends
+ altitude_to_visualizer_for_plot_trend = {a: v for a, v in altitude_to_visualizer.items() if a >= 900}
+ plot_trend_curves(altitude_to_visualizer_for_plot_trend)
# Plot graph
- plot_uncertainty_massifs(altitude_to_visualizer)
- # Plot histogram
- plot_uncertainty_histogram(altitude_to_visualizer)
+ # plot_uncertainty_massifs(altitude_to_visualizer)
+ # # Plot histogram
+ # plot_uncertainty_histogram(altitude_to_visualizer)
def major_result():
@@ -98,11 +104,11 @@ def major_result():
if __name__ == '__main__':
- major_result()
- # intermediate_result(altitudes=[900, 1200], massif_names=['Chartreuse'],
- # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
- # ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
- # multiprocessing=True)
+ # major_result()
+ intermediate_result(altitudes=ALL_ALTITUDES_WITHOUT_NAN[2:], massif_names=None,
+ uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+ ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+ multiprocessing=True)
# intermediate_result(altitudes=[900, 1200], massif_names=['Maurienne'],
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
# ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_trend_curves.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_trend_curves.py
new file mode 100644
index 0000000000000000000000000000000000000000..4e3d7de039375a22f43179969dd8273ebcff59e2
--- /dev/null
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_trend_curves.py
@@ -0,0 +1,78 @@
+from typing import Dict
+import matplotlib.pyplot as plt
+
+from experiment.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
+from experiment.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
+from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure
+from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
+ StudyVisualizerForNonStationaryTrends
+
+
+def plot_trend_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+ """
+ Plot a single trend curves
+ :return:
+ """
+ visualizer = list(altitude_to_visualizer.values())[0]
+
+ ax = create_adjusted_axes(1, 1)
+ ax_twinx = ax.twinx()
+ ax_twiny = ax.twiny()
+
+ trend_summary_values = list(zip(*[v.trend_summary_values() for v in altitude_to_visualizer.values()]))
+ altitudes, percent_decrease, percent_decrease_signi, *mean_decreases = trend_summary_values
+
+ # parameters
+ width = 150
+ size = 20
+ legend_fontsize = 20
+ color = 'white'
+ labelsize = 15
+ linewidth = 3
+ ax.bar(altitudes, percent_decrease, width=width, color=color, edgecolor='blue', label='decreasing trend',
+ linewidth=linewidth)
+ ax.bar(altitudes, percent_decrease_signi, width=width, color=color, edgecolor='black',
+ label='significative decreasing trend',
+ linewidth=linewidth)
+ ax.legend(loc='upper left', prop={'size': size})
+
+ for ax_horizontal in [ax, ax_twiny]:
+ if ax_horizontal == ax_twiny:
+ ax_horizontal.plot(altitudes, [0 for _ in altitudes], linewidth=0)
+ else:
+ ax_horizontal.set_xlabel('Altitude', fontsize=legend_fontsize)
+ ax_horizontal.set_xticks(altitudes)
+ ax_horizontal.set_xlim([700, 5000])
+ ax_horizontal.tick_params(labelsize=labelsize)
+
+ # Set the number of massifs on the upper axis
+ ax_twiny.set_xticklabels([v.study.nb_study_massif_names for v in altitude_to_visualizer.values()])
+ ax_twiny.set_xlabel('Total number of massif at each altitude (for the percentage)', fontsize=legend_fontsize)
+
+ ax.set_ylabel('Percentage of massifs with decreasing trend (\%)', fontsize=legend_fontsize)
+ max_percent = int(max(percent_decrease))
+ n = 2 + (max_percent // 10)
+ ax_ticks = [10 * i for i in range(n)]
+ # upper_lim = max_percent + 3
+ upper_lim = n + 5
+ ax_lim = [0, upper_lim]
+ for axis in [ax, ax_twinx]:
+ axis.set_ylim(ax_lim)
+ axis.set_yticks(ax_ticks)
+ axis.tick_params(labelsize=labelsize)
+ ax.yaxis.grid()
+
+ label_curve = (visualizer.label).replace('change', 'decrease')
+ ax_twinx.set_ylabel(label_curve.replace('', ''), fontsize=legend_fontsize)
+ for region_name, mean_decrease in zip(AbstractExtendedStudy.region_names, mean_decreases):
+ if len(mean_decreases) > 1:
+ label = region_name
+ else:
+ label = 'Mean decrease'
+ ax_twinx.plot(altitudes, mean_decrease, label=label, linewidth=linewidth, marker='o')
+ ax_twinx.legend(loc='upper right', prop={'size': size})
+
+ # Save plot
+ visualizer.plot_name = 'Trend curves'
+ visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure)
+ plt.close()
diff --git a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
index efd0742ef5d8b823730ef6faaece0e56c499fe7f..0116fe33f7ea201cae951c24a2ac70aa908d00d0 100644
--- a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
@@ -66,7 +66,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
self.uncertainty_massif_names = uncertainty_massif_names
# Assign some default arguments
if self.model_subsets_for_uncertainty is None:
- self.model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
+ self.model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
ModelSubsetForUncertainty.non_stationary_gumbel_and_gev][:]
if self.uncertainty_methods is None:
self.uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
@@ -153,7 +153,8 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
# Extract the stationary model that minimized AIC
stationary_trend_test_that_minimized_aic = [t for t in sorted_trend_test if type(t) in
[GumbelVersusGumbel, GevStationaryVersusGumbel]][0]
- massif_name_to_stationary_trend_test_that_minimized_aic[massif_name] = stationary_trend_test_that_minimized_aic
+ massif_name_to_stationary_trend_test_that_minimized_aic[
+ massif_name] = stationary_trend_test_that_minimized_aic
# Extract the Gumbel model that minimized AIC
gumbel_trend_test_that_minimized_aic = [t for t in sorted_trend_test if type(t) in
[GumbelVersusGumbel, GumbelLocationTrendTest, GumbelScaleTrendTest,
@@ -300,7 +301,8 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
else:
raise ValueError(model_subset_for_uncertainty)
- def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, model_subset_for_uncertainty) \
+ def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method,
+ model_subset_for_uncertainty) \
-> Dict[str, EurocodeConfidenceIntervalFromExtremes]:
# Compute for the uncertainty massif names
arguments = [
@@ -390,7 +392,15 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
percentage_decrease = 100 * len(decreasing_trend_tests) / len(trend_tests)
significative_decrease_trend_tests = [t for t in decreasing_trend_tests if t.is_significant]
percentage_decrease_significative = 100 * len(significative_decrease_trend_tests) / len(trend_tests)
- massif_name_to_region_name = AbstractExtendedStudy.massif_name_to_region_name
- mean_relative = np.mean(np.array(list(self.massif_name_to_relative_change_value.values())))
-
- return self.altitude, percentage_decrease, percentage_decrease_significative
\ No newline at end of file
+ # For visualization at 2700m
+ if percentage_decrease_significative == percentage_decrease:
+ percentage_decrease += 0.4
+ compute_mean_decrease = lambda l: -np.mean(np.array(list(l)))
+ mean_decreases = [compute_mean_decrease(self.massif_name_to_relative_change_value.values())]
+ # Compute mean relatives per regions (for the moment i don't add the region means)
+ # massif_name_to_region_name = AbstractExtendedStudy.massif_name_to_region_name
+ # for region_name in AbstractExtendedStudy.real_region_names:
+ # change_values = [v for m, v in self.massif_name_to_relative_change_value.items()
+ # if massif_name_to_region_name[m] == region_name]
+ # mean_decreases.append(compute_mean_decrease(change_values))
+ return (self.altitude, percentage_decrease, percentage_decrease_significative, *mean_decreases)