[paper 2] add plot curves (to highlight the decreasing trend) v2

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

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:],

experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_trend_curves.py

+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()

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)