[contrasting project] add comparative_curve_wrt_altitude.py to provide box plot wrt the altitude

projects/contrasting_trends_in_snow_loads/gorman_figures/figure1/comparative_curve_wrt_altitude.py

+from collections import OrderedDict
+import matplotlib.pyplot as plt
+import numpy as np
+
+from cached_property import cached_property
+
+from projects.contrasting_trends_in_snow_loads.gorman_figures.figure1.study_visualizer_for_double_stationary_fit import \
+    StudyVisualizerForReturnLevelChange
+
+
+class ComparativeCurveWrtAltitude(object):
+
+    def __init__(self, altitude_to_study_visualizer):
+        assert isinstance(altitude_to_study_visualizer, OrderedDict)
+        self.altitude_to_study_visualizer = altitude_to_study_visualizer  # type: OrderedDict[int, StudyVisualizerForReturnLevelChange]
+
+    @cached_property
+    def visualizer(self):
+        return list(self.altitude_to_study_visualizer.values())[0]
+
+    @property
+    def altitudes(self):
+        return list(self.altitude_to_study_visualizer.keys())
+
+    def all_plots(self):
+        self.shape_shoe_plot()
+        self.return_level_shot_plot()
+
+    def shape_shoe_plot(self):
+        altitude_to_list_couple = OrderedDict()
+        for a, v in self.altitude_to_study_visualizer.items():
+            altitude_to_list_couple[a] = v.result_from_double_stationary_fit.shape_list_couple
+        label = 'shape parameter'
+        self.abstract_shoe_plot(altitude_to_list_couple, plot_name=label, ylabel=label, relative_change=False)
+
+    def return_level_shot_plot(self):
+        altitude_to_list_couple = OrderedDict()
+        for a, v in self.altitude_to_study_visualizer.items():
+            altitude_to_list_couple[a] = v.result_from_double_stationary_fit.return_level_list_couple
+        label = '{}-year return level ({})'.format(self.visualizer.return_period, self.visualizer.study.variable_unit)
+        self.abstract_shoe_plot(altitude_to_list_couple, plot_name='return level', ylabel=label)
+
+    def abstract_shoe_plot(self, altitude_to_list_couple, plot_name, ylabel, relative_change=True):
+        # Prepare axis
+        ax = plt.gca()
+        ax2 = ax.twinx()
+
+        # Prepare data
+        altitude_to_all_relative_differences = OrderedDict()
+        altitude_to_all_differences = OrderedDict()
+        altitude_to_mean_value_before = OrderedDict()
+        altitude_to_mean_value_after = OrderedDict()
+        for altitude in self.altitudes:
+            list_couple = altitude_to_list_couple[altitude]
+            before, after = [np.array(a) for a in zip(*list_couple)]
+            altitude_to_mean_value_before[altitude] = before.mean()
+            altitude_to_mean_value_after[altitude] = after.mean()
+            altitude_to_all_relative_differences[altitude] = 100 * (after - before) / before
+            altitude_to_all_differences[altitude] = after - before
+
+        # ax2: Boxplot
+
+        width = 100
+        x = altitude_to_all_relative_differences if relative_change else altitude_to_all_differences
+        ax2.boxplot(x.values(), positions=self.altitudes, widths=width)
+        ax2.set_xlim([min(self.altitudes) - width, max(self.altitudes) + width])
+        if relative_change:
+            ax2.set_ylabel('Relative change in {} (%)'.format((ylabel.split('(')[0])))
+        else:
+            ax2.set_ylabel('Change in {}'.format(ylabel))
+
+        # ax: Mean plots on top
+        values = [list(v.values()) for v in [altitude_to_mean_value_before, altitude_to_mean_value_after]]
+        labels = [
+            '{}-{}'.format(self.visualizer.year_min_before, self.visualizer.year_max_before),
+            '{}-{}'.format(self.visualizer.year_min_after, self.visualizer.year_max_after),
+        ]
+        for label, value in zip(labels, values):
+            ax.plot(self.altitudes, value, label=label, marker='o')
+        ax.set_ylabel('Mean {}'.format(ylabel))
+        ax.legend(loc='upper left')  # , prop={'size': size})
+        ax.set_xlabel('Altitude (m)')
+
+        # Set same limits to align axis
+        a = np.array(list(x.values()))
+        all_values = list(a.flatten()) + values[0] + values[1]
+        print(type(all_values))
+        epsilon = 0.03 * (max(all_values) - min(all_values))
+        ylim = [min(all_values) - epsilon, max(all_values) + epsilon]
+        ax.set_ylim(ylim)
+        ax2.set_ylim(ylim)
+        ax.grid()
+
+        self.visualizer.plot_name = plot_name + ' summary'
+        self.visualizer.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True,
+                                             dpi=500)
+        plt.close()

projects/contrasting_trends_in_snow_loads/gorman_figures/figure1/figure1_mean_ratio_return_level_ratio.py

+from collections import OrderedDict
+
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad1Day, CrocusSnowLoad3Days, \
+    CrocusSnowLoadTotal
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \
+    SafranSnowfall7Days, SafranSnowfall5Days
+from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import ALL_ALTITUDES_WITHOUT_NAN
+from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
+    TemporalMarginFitMethod
+from projects.contrasting_trends_in_snow_loads.gorman_figures.figure1.comparative_curve_wrt_altitude import \
+    ComparativeCurveWrtAltitude
+
+from projects.contrasting_trends_in_snow_loads.gorman_figures.figure1.study_visualizer_for_double_stationary_fit import \
+    StudyVisualizerForReturnLevelChange
+
+
+def load_altitude_to_study_visualizer(study_class, save_to_file=True) -> OrderedDict:
+    altitude_to_study_visualizer = OrderedDict()
+    # for altitude in ALL_ALTITUDES_WITHOUT_NAN[2:10]:
+    for altitude in ALL_ALTITUDES_WITHOUT_NAN[2:5]:
+        return_period = 30
+        study_visualizer = StudyVisualizerForReturnLevelChange(study_class=study_class,
+                                                               altitude=altitude,
+                                                               return_period=return_period,
+                                                               save_to_file=save_to_file,
+                                                               fit_method=TemporalMarginFitMethod.extremes_fevd_l_moments)
+        altitude_to_study_visualizer[altitude] = study_visualizer
+    return altitude_to_study_visualizer
+
+
+def plots():
+    for study_class in [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:1]:
+        altitude_to_study_visualizer = load_altitude_to_study_visualizer(study_class, save_to_file=True)
+        # for v in altitude_to_study_visualizer.values():
+        #     v.all_plots()
+        comparative_curve = ComparativeCurveWrtAltitude(altitude_to_study_visualizer)
+        comparative_curve.all_plots()
+
+
+if __name__ == '__main__':
+    plots()

projects/contrasting_trends_in_snow_loads/gorman_figures/result_from_stationary_fit.py → projects/contrasting_trends_in_snow_loads/gorman_figures/figure1/result_from_stationary_fit.py

@@ -16,6 +16,10 @@ class ResultFromDoubleStationaryFit(object):
         self.result_before = ResultFromSingleStationaryFit(study_before, fit_method, return_period)
         self.result_after = ResultFromSingleStationaryFit(study_after, fit_method, return_period)
 
+    @property
+    def massif_names(self):
+        return self.re
+
     @cached_property
     def massif_name_to_difference_return_level(self):
         return {m: v - self.result_before.massif_name_to_return_level[m]
@@ -26,6 +30,16 @@ class ResultFromDoubleStationaryFit(object):
         return {m: 100 * v / self.result_before.massif_name_to_return_level[m]
                 for m, v in self.massif_name_to_difference_return_level.items()}
 
+    @cached_property
+    def return_level_list_couple(self):
+        return [(v, self.result_after.massif_name_to_return_level[m])
+                for m, v in self.result_before.massif_name_to_return_level.items()]
+
+    @cached_property
+    def shape_list_couple(self):
+        return [(v, self.result_after.massif_name_to_shape[m])
+                for m, v in self.result_before.massif_name_to_shape.items()]
+
 
 class ResultFromSingleStationaryFit(object):
 
@@ -54,4 +68,4 @@ class ResultFromSingleStationaryFit(object):
 
     @property
     def massif_name_to_difference_return_level_and_maxima(self):
-        return {m: r - np.max(self.massif_name_to_maxima[m]) for m, r in self.massif_name_to_return_level.items() }
+        return {m: r - np.max(self.massif_name_to_maxima[m]) for m, r in self.massif_name_to_return_level.items()}

projects/contrasting_trends_in_snow_loads/gorman_figures/figure1_mean_ratio_return_level_ratio.py → projects/contrasting_trends_in_snow_loads/gorman_figures/figure1/study_visualizer_for_double_stationary_fit.py

-from collections import OrderedDict
-
 import matplotlib
-import numpy as np
-from cached_property import cached_property
 
 from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad1Day, CrocusSnowLoad3Days, \
     CrocusSnowLoadTotal
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days
 from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import get_shifted_map, \
-    get_colors, shiftedColorMap, ticks_values_and_labels_for_percentages
+    get_colors, ticks_values_and_labels_for_percentages
 from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import ALL_ALTITUDES_WITHOUT_NAN
 from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import StudyVisualizer
-from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
 from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
     TemporalMarginFitMethod, fitmethod_to_str
 
 import matplotlib.pyplot as plt
 
-from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends
-from projects.contrasting_trends_in_snow_loads.gorman_figures.result_from_stationary_fit import \
+from projects.contrasting_trends_in_snow_loads.gorman_figures.figure1.result_from_stationary_fit import \
     ResultFromDoubleStationaryFit
 
 
@@ -161,17 +155,3 @@ class StudyVisualizerForReturnLevelChange(StudyVisualizer):
         ax.set_xticks([])
         ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=15)
         ax.set_title('Fit method is {}'.format(fitmethod_to_str(self.fit_method)))
-
-
-if __name__ == '__main__':
-    # for altitude in [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:1]:
-    for altitude in ALL_ALTITUDES_WITHOUT_NAN:
-        for study_class in [SafranSnowfall1Day, SafranSnowfall3Days, CrocusSnowLoadTotal, CrocusSnowLoad3Days,
-                            CrocusSnowLoad1Day][:1]:
-            return_period = 30
-            study_visualizer = StudyVisualizerForReturnLevelChange(study_class=study_class,
-                                                                   altitude=altitude,
-                                                                   return_period=return_period,
-                                                                   save_to_file=True,
-                                                                   fit_method=TemporalMarginFitMethod.extremes_fevd_l_moments)
-            study_visualizer.all_plots()