[PAPER 1] add ticks to the colorbar for the trend. add histogram & marker to the uncertainty plot

experiment/meteo_france_data/plot/create_shifted_cmap.py

@@ -24,10 +24,13 @@ def get_shifted_map(vmin, vmax):
     return shifted_cmap
 
 
-def create_colorbase_axis(ax, label, cmap, norm):
+def create_colorbase_axis(ax, label, cmap, norm, ticks_values_and_labels=None):
     divider = make_axes_locatable(ax)
     cax = divider.append_axes('right', size='5%', pad=0.0)
-    cb = cbar.ColorbarBase(cax, cmap=cmap, norm=norm)
+    ticks = ticks_values_and_labels[0] if ticks_values_and_labels is not None else None
+    cb = cbar.ColorbarBase(cax, cmap=cmap, norm=norm, ticks=ticks)
+    if ticks_values_and_labels is not None:
+        cb.ax.set_yticklabels([str(t) for t in ticks_values_and_labels[1]])
     if isinstance(label, str):
         cb.set_label(label)
     return norm

experiment/meteo_france_data/scm_models_data/abstract_study.py

@@ -306,6 +306,7 @@ class AbstractStudy(object):
                         massif_name_to_hatch_boolean_list=None,
                         norm=None,
                         massif_name_to_marker_style=None,
+                        ticks_values_and_labels=None,
                         ):
         if ax is None:
             ax = plt.gca()
@@ -391,7 +392,7 @@ class AbstractStudy(object):
         # create the colorbar only at the end
         if add_colorbar:
             if len(set(values)) > 1:
-                create_colorbase_axis(ax, label, cmap, norm)
+                create_colorbase_axis(ax, label, cmap, norm, ticks_values_and_labels=ticks_values_and_labels)
         if axis_off:
             plt.axis('off')
 

experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py

@@ -38,9 +38,10 @@ def plot_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVisualizerFo
     """
     visualizer = list(altitude_to_visualizer.values())[0]
     # Subdivide massif names in group of 3
-    m = 3
+    m = 1
     uncertainty_massif_names = visualizer.uncertainty_massif_names
     n = (len(uncertainty_massif_names) // m) + 1
+    print('total nb of massif', n)
     for i in list(range(n))[:]:
         massif_names = uncertainty_massif_names[m * i: m * (i + 1)]
         print(massif_names)
@@ -97,7 +98,7 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
     """ Generic function that might be used by many other more global functions"""
     altitudes = list(altitude_to_visualizer.keys())
     visualizer = list(altitude_to_visualizer.values())[0]
-    colors = ['tab:green', 'tab:olive']
+    colors = ['tab:green', 'tab:brown'][::-1]
     alpha = 0.2
     # Display the EUROCODE return level
     eurocode_region = massif_name_to_eurocode_region[massif_name]()
@@ -106,14 +107,27 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
     for j, (color, uncertainty_method) in enumerate(zip(colors, visualizer.uncertainty_methods)):
         # Plot eurocode standards only for the first loop
         if j == 0:
+            # Plot eurocode norm
             eurocode_region.plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(ax,
                                                                                                    altitudes=altitudes)
+            # Plot bars of TDRL only in the non stationary case
+            if non_stationary_context:
+                tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in altitude_to_visualizer.values()]
+                # Plot bars
+                colors = [v.massif_name_to_tdrl_color[massif_name] for v in altitude_to_visualizer.values()]
+                ax.bar(altitudes, tdrl_values, width=150, color=colors)
+                # Plot markers
+                markers_kwargs = [v.massif_name_to_marker_style(markersize=4)[massif_name]
+                                  for v, tdrl_value in zip(altitude_to_visualizer.values(), tdrl_values)]
+                for altitude, marker_kwargs, value in zip(altitudes, markers_kwargs, tdrl_values):
+                    ax.plot([altitude], [value / 2], **marker_kwargs)
+
         # Plot uncertainties
         plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,
                                               non_stationary_context, uncertainty_method)
 
     ax.legend(loc=2)
-    ax.set_ylim([0.0, 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,

experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

@@ -4,6 +4,7 @@ import os.path as op
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 
+from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map
 from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
     StudyVisualizer
 from experiment.paper_past_snow_loads.result_trends_and_return_levels.eurocode_visualizer import \
@@ -23,6 +24,7 @@ def minor_result(altitude):
     """Plot trends for a single altitude to be fast"""
     visualizer = StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude), multiprocessing=True)
     visualizer.plot_trends()
+    plt.show()
 
 
 def intermediate_result(altitudes, massif_names=None,
@@ -67,16 +69,17 @@ def major_result():
 
 if __name__ == '__main__':
     # minor_result(altitude=1800)
-    intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
-                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
-                        non_stationary_uncertainty=[False])
+    # intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
+    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
+    #                                          ConfidenceIntervalMethodFromExtremes.ci_bayes],
+    #                     non_stationary_uncertainty=[True])
     # intermediate_result(altitudes=[1500, 1800], massif_names=None,
     #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
     #                     non_stationary_uncertainty=[False])
     # intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
     #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
     #                     non_stationary_uncertainty=[False])
-    # intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=['Vercors', 'Oisans', 'Devoluy'],
-    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
-    #                                          ConfidenceIntervalMethodFromExtremes.ci_bayes],
-    #                     non_stationary_uncertainty=[False, True])
+    intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
+                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
+                                             ConfidenceIntervalMethodFromExtremes.ci_bayes],
+                        non_stationary_uncertainty=[False, True])

experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py

@@ -7,6 +7,7 @@ import numpy as np
 from cached_property import cached_property
 
 from experiment.eurocode_data.utils import EUROCODE_QUANTILE
+from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map, get_colors
 from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
     StudyVisualizer
@@ -55,6 +56,8 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         self.non_stationary_trend_test = [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest]
         self.non_stationary_trend_test_to_marker = dict(zip(self.non_stationary_trend_test, ["s", "^", "D"]))
 
+        self.global_max_abs_tdrl = None
+
     # Utils
 
     @cached_property
@@ -91,8 +94,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         massif_name_to_trend_test_that_minimized_aic = {}
         for massif_name, (x, y) in self.massif_name_to_non_null_years_and_maxima.items():
             quantile_level = self.massif_name_to_eurocode_quantile_level_in_practice[massif_name]
-            non_stationary_trend_test = [t(years=x, maxima=y, starting_year=starting_year, quantile_level=quantile_level)
-                                         for t in self.non_stationary_trend_test]
+            non_stationary_trend_test = [
+                t(years=x, maxima=y, starting_year=starting_year, quantile_level=quantile_level)
+                for t in self.non_stationary_trend_test]
             trend_test_that_minimized_aic = sorted(non_stationary_trend_test, key=lambda t: t.aic)[0]
             massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
         return massif_name_to_trend_test_that_minimized_aic
@@ -103,32 +107,62 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     def max_abs_tdrl(self):
         return max(abs(min(self.massif_name_to_tdrl_value.values())), max(self.massif_name_to_tdrl_value.values()))
 
+    @cached_property
+    def _max_abs_tdrl(self):
+        return self.global_max_abs_tdrl if self.global_max_abs_tdrl is not None else self.max_abs_tdrl
+
     def plot_trends(self, max_abs_tdrl=None):
-        if max_abs_tdrl is None:
-            max_abs_tdrl = self.max_abs_tdrl
-        assert max_abs_tdrl > 0
+        if max_abs_tdrl is not None:
+            self.global_max_abs_tdrl = max_abs_tdrl
         self.study.visualize_study(massif_name_to_value=self.massif_name_to_tdrl_value,
-                                   vmin=-max_abs_tdrl, vmax=max_abs_tdrl,
-                                   replace_blue_by_white=False, axis_off=True, show_label=False,
+                                   replace_blue_by_white=False, axis_off=False, show_label=False,
                                    add_colorbar=True,
-                                   massif_name_to_marker_style=self.massif_name_to_marker_style,
-                                   show=self.show)
+                                   massif_name_to_marker_style=self.massif_name_to_marker_style(),
+                                   massif_name_to_color=self.massif_name_to_tdrl_color,
+                                   cmap=self.cmap,
+                                   show=self.show,
+                                   ticks_values_and_labels=self.ticks_values_and_labels,
+                                   label=self.label_trend_bar)
         self.plot_name = 'tdlr_trends'
         self.show_or_save_to_file(add_classic_title=False, tight_layout=False, no_title=True)
         plt.close()
 
+    @property
+    def label_trend_bar(self):
+        return 'Change in {} years for Eurocode quantile of SL (kN $m^-2$)'.format(AbstractGevTrendTest.nb_years_for_quantile_evolution)
+
+    @property
+    def ticks_values_and_labels(self):
+        graduation = 0.2
+        positive_ticks = []
+        tick = graduation
+        while tick < self._max_abs_tdrl:
+            positive_ticks.append(round(tick, 1))
+            tick += graduation
+        all_ticks_labels = [-t for t in positive_ticks] + [0] + positive_ticks
+        ticks_values = [(t / 2 * self._max_abs_tdrl) + 0.5 for t in all_ticks_labels]
+        return ticks_values, all_ticks_labels
+
     @cached_property
     def massif_name_to_tdrl_value(self):
         return {m: t.time_derivative_of_return_level for m, t in
                 self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
+    
+    @cached_property
+    def cmap(self):
+        return get_shifted_map(-self._max_abs_tdrl, self._max_abs_tdrl)
 
-    @property
-    def massif_name_to_marker_style(self):
+    @cached_property
+    def massif_name_to_tdrl_color(self):
+        return {m: get_colors([v], self.cmap, -self._max_abs_tdrl, self._max_abs_tdrl)[0]
+                for m, v in self.massif_name_to_tdrl_value.items()}
+
+    def massif_name_to_marker_style(self, markersize=5):
         d = {}
         for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items():
             d[m] = {'marker': self.non_stationary_trend_test_to_marker[type(t)],
                     'color': 'k',
-                    'markersize': 5,
+                    'markersize': markersize,
                     'fillstyle': 'full' if t.is_significant else 'none'}
         return d
 

extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py

 import numpy as np
 import pandas as pd
+import rpy2
 from rpy2 import robjects
 
 from extreme_fit.distribution.gev.gev_params import GevParams
@@ -48,7 +49,10 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
             qcov = r("make.qcov")(self.result_from_fit,
                                   **kwargs)
             common_kwargs['qcov'] = qcov
-        mean_estimate, confidence_interval = self._confidence_interval_method(common_kwargs, ci_method, return_period)
+        try:
+            mean_estimate, confidence_interval = self._confidence_interval_method(common_kwargs, ci_method, return_period)
+        except rpy2.rinterface.RRuntimeError:
+            mean_estimate, confidence_interval = np.nan, (np.nan, np.nan)
         return mean_estimate, confidence_interval
 
     def _confidence_interval_method(self, common_kwargs, ci_method, return_period):

extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py

@@ -24,10 +24,7 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
                 'method': method_name,
             # xrange = NULL, nint = 20
         }
-        try:
-            res = r.ci(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
-        except rpy2.rinterface.RRuntimeError:
-            return np.nan, (np.nan, np.nan)
+        res = r.ci(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
         if self.is_non_stationary:
             a = np.array(res)[0]
             lower, mean_estimate, upper, _ = a