[projection swe] modify interval (and data constraint) for temperature...

[projection swe] modify interval (and data constraint) for temperature sensitivity. add invidual plots in sensitivity. add changes plot in sensitivity.

extreme_data/meteo_france_data/adamont_data/cmip5/plot_temperatures.py

@@ -88,6 +88,6 @@ def plot_temperature_for_rcp_gcm(ax, gcm, scenario, year_min, year_max, linestyl
 
 
 if __name__ == '__main__':
-    for anomaly in [True, False][:]:
+    for anomaly in [True, False][:1]:
         main_plot_temperature_with_spline_on_top(anomaly=anomaly)
         # main_plot_temperature_all(anomaly=True, spline=True)

extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py

@@ -18,8 +18,8 @@ def get_year_min_and_year_max(gcm, scenario, left_limit, right_limit, is_tempera
     else:
         year_min, year_max = left_limit, right_limit
 
-    # A minimum of 30 years of data is needed to find a trend
-    if year_max - year_min + 1 >= 30:
+    # A minimum of 10 years of data is needed for each GCM/RCM
+    if year_max - year_min + 1 >= 10:
         return year_min, year_max
     else:
         return None, None
@@ -70,10 +70,11 @@ def plot_nb_data(is_temperature_interval, is_shift_interval):
     ax = plt.gca()
     for gcm in get_gcm_list(adamont_version=2)[:]:
         for i, scenario in enumerate(rcp_scenarios[2:]):
+            first_scenario = i == 0
             plot_nb_data_one_line(ax, gcm, scenario, left_limit, right_limit,
-                                  i == 0, is_temperature_interval)
+                                  first_scenario, is_temperature_interval)
 
-    ax.legend()
+    ax.legend(loc='upper left')
     ticks_labels = get_ticks_labels_for_interval(is_temperature_interval, is_shift_interval)
     ax.set_xticks(right_limit)
     ax.set_xticklabels(ticks_labels)
@@ -84,15 +85,15 @@ def plot_nb_data(is_temperature_interval, is_shift_interval):
         Line2D([0], [0], color='k', lw=1, label=scenario_to_str(s),
                linestyle=get_linestyle_from_scenario(s)) for s in adamont_scenarios_real
     ]
-    ax2.legend(handles=legend_elements, loc='upper center')
+    ax2.legend(handles=legend_elements, loc='lower right')
     ax2.set_yticks([])
     plt.show()
 
 
 def get_interval_limits(is_temperature_interval, is_shift_interval):
     if is_temperature_interval:
-        temp_min = np.arange(0, 3, 1)
-        temp_max = temp_min + 2
+        temp_min = np.arange(0, 2, 0.5)
+        temp_max = temp_min + 1.5
         left_limit, right_limit = temp_min, temp_max
     else:
         shift = 25
@@ -103,7 +104,7 @@ def get_interval_limits(is_temperature_interval, is_shift_interval):
     if not is_shift_interval:
         min_interval_left = min(left_limit)
         left_limit = [min_interval_left for _ in right_limit]
-    return left_limit, right_limit
+    return left_limit[:2], right_limit[:2]
 
 
 def get_ticks_labels_for_interval(is_temperature_interval, is_shift_interval):
@@ -117,7 +118,7 @@ def get_ticks_labels_for_interval(is_temperature_interval, is_shift_interval):
 
 
 if __name__ == '__main__':
-    for shift_interval in [False, True]:
+    for shift_interval in [False, True][:1]:
         for temp_interval in [False, True][1:]:
             print("shift = {}, temp_inteval = {}".format(shift_interval, temp_interval))
             plot_nb_data(is_temperature_interval=temp_interval, is_shift_interval=shift_interval)

extreme_trend/ensemble_fit/abstract_ensemble_fit.py

@@ -29,7 +29,7 @@ class AbstractEnsembleFit(object):
 
         # Set appropriate setting
         OneFoldFit.last_year = 2100
-        OneFoldFit.nb_years = 95
+        OneFoldFit.nb_years = OneFoldFit.last_year - 2005
 
     @property
     def altitudes(self):

extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py

@@ -31,7 +31,9 @@ class VisualizerForProjectionEnsemble(object):
                  confidence_interval_based_on_delta_method=False,
                  remove_physically_implausible_models=False,
                  gcm_to_year_min_and_year_max=None,
+                 interval_str_prefix='',
                  ):
+        self.interval_str_prefix = interval_str_prefix
         self.altitudes_list = altitudes_list
         self.temporal_covariate_for_fit = temporal_covariate_for_fit
         self.scenario = scenario
@@ -119,14 +121,14 @@ class VisualizerForProjectionEnsemble(object):
                                ]
             if key in merge_keys:
                 for v in visualizer_list:
-                    v.studies.study.gcm_rcm_couple = (key, "merge")
+                    v.studies.study.gcm_rcm_couple = ("{} {}".format(key, "merge"), self.interval_str_prefix)
             self.plot_for_visualizer_list(visualizer_list)
 
     def plot_together(self):
         visualizer_list = [together_ensemble_fit.visualizer
                            for together_ensemble_fit in self.ensemble_fits(TogetherEnsembleFit)]
         for v in visualizer_list:
-            v.studies.study.gcm_rcm_couple = ("together", "merge")
+            v.studies.study.gcm_rcm_couple = ("together merge", self.interval_str_prefix)
         self.plot_for_visualizer_list(visualizer_list)
 
     def ensemble_fits(self, ensemble_class):

extreme_trend/ensemble_fit/visualizer_for_sensitivity.py

@@ -2,6 +2,8 @@ from collections import OrderedDict
 import matplotlib.pyplot as plt
 from typing import List, Dict
 
+import numpy as np
+
 from extreme_data.meteo_france_data.adamont_data.cmip5.temperature_to_year import get_interval_limits, \
     get_year_min_and_year_max, get_ticks_labels_for_interval
 from extreme_data.meteo_france_data.scm_models_data.utils import Season
@@ -40,11 +42,14 @@ class VisualizerForSensivity(object):
         self.massif_names = massif_names
         self.left_limits, self.right_limits = get_interval_limits(self.is_temperature_interval,
                                                                   self.is_shift_interval)
+
         self.left_limit_to_right_limit = OrderedDict(zip(self.left_limits, self.right_limits))
         self.right_limit_to_visualizer = {}  # type: Dict[float, VisualizerForProjectionEnsemble]
 
         for left_limit, right_limit in zip(self.left_limits, self.right_limits):
-            print("Interval is", left_limit, right_limit)
+            interval_str_prefix = "{}-{}".format(left_limit, right_limit)
+            interval_str_prefix = interval_str_prefix.replace('.', ',')
+            print("Interval is", interval_str_prefix)
             # Build gcm_to_year_min_and_year_max
             gcm_to_year_min_and_year_max = {}
             gcm_list = list(set([g for g, r in gcm_rcm_couples]))
@@ -64,29 +69,61 @@ class VisualizerForSensivity(object):
                 massif_names=massif_names,
                 temporal_covariate_for_fit=temporal_covariate_for_fit,
                 remove_physically_implausible_models=remove_physically_implausible_models,
-                gcm_to_year_min_and_year_max=gcm_to_year_min_and_year_max
+                gcm_to_year_min_and_year_max=gcm_to_year_min_and_year_max,
+                interval_str_prefix=interval_str_prefix,
             )
             self.right_limit_to_visualizer[right_limit] = visualizer
 
     def plot(self):
-        # todo: before reactivating the subplot, i should ensure that we can modify the prefix
-        # todo: so that we can have all the subplot for the merge visualizer
-        # , and not just the plots for the last t_min
-        # for visualizer in self.temp_min_to_visualizer.values():
-        #     visualizer.plot()
+        for visualizer in self.right_limit_to_visualizer.values():
+            visualizer.plot()
         merge_visualizer_str_list = []
         if IndependentEnsembleFit in self.ensemble_fit_classes:
             merge_visualizer_str_list.extend([AbstractEnsembleFit.Median_merge, AbstractEnsembleFit.Mean_merge])
         if TogetherEnsembleFit in self.ensemble_fit_classes:
             merge_visualizer_str_list.append(AbstractEnsembleFit.Together_merge)
         for merge_visualizer_str in merge_visualizer_str_list:
-            self.sensitivity_plot(merge_visualizer_str)
+            self.sensitivity_plot_percentages(merge_visualizer_str)
+            for relative in [True, False]:
+                self.sensitivity_plot_changes(merge_visualizer_str, relative)
 
-    def sensitivity_plot(self, merge_visualizer_str):
+    def sensitivity_plot_changes(self, merge_visualizer_str, relative):
         ax = plt.gca()
         for altitudes in self.altitudes_list:
             altitude_class = get_altitude_class_from_altitudes(altitudes)
-            self.interval_plot(ax, altitude_class, merge_visualizer_str)
+            self.interval_plot_changes(ax, altitude_class, merge_visualizer_str, relative)
+
+        ticks_labels = get_ticks_labels_for_interval(self.is_temperature_interval, self.is_shift_interval)
+        name = 'relative changes' if relative else "changes"
+        ax.set_ylabel(name)
+        ax.set_xlabel('Interval used to compute the trends ')
+        ax.set_xticks(self.right_limits)
+        ax.set_xticklabels(ticks_labels)
+        lim_left, lim_right = ax.get_xlim()
+        ax.hlines(0, xmin=lim_left, xmax=lim_right, linestyles='dashed')
+        ax.legend(prop={'size': 7}, loc='upper center', ncol=2)
+        # ax.set_ylim((0, 122))
+        # ax.set_yticks([i * 10 for i in range(11)])
+        self.save_plot(merge_visualizer_str, name)
+
+    def interval_plot_changes(self, ax, altitude_class, merge_visualizer_str, relative):
+        linestyle = get_linestyle_for_altitude_class(altitude_class)
+
+        mean_changes = []
+        for v in self.right_limit_to_visualizer.values():
+            merge_visualizer = self.get_merge_visualizer(altitude_class, v, merge_visualizer_str)
+            changes, non_stationary_changes = merge_visualizer.all_changes(self.massif_names, relative=relative,
+                                                                           with_significance=False)
+            mean_changes.append(np.mean(changes))
+
+        label = altitude_class().formula
+        ax.plot(self.right_limits, mean_changes, linestyle=linestyle, label=label, color='darkgreen')
+
+    def sensitivity_plot_percentages(self, merge_visualizer_str):
+        ax = plt.gca()
+        for altitudes in self.altitudes_list:
+            altitude_class = get_altitude_class_from_altitudes(altitudes)
+            self.interval_plot_percentages(ax, altitude_class, merge_visualizer_str)
 
         ticks_labels = get_ticks_labels_for_interval(self.is_temperature_interval, self.is_shift_interval)
         ax.set_ylabel('Percentages of massifs (\%)')
@@ -96,16 +133,9 @@ class VisualizerForSensivity(object):
         ax.legend(prop={'size': 7}, loc='upper center', ncol=2)
         ax.set_ylim((0, 122))
         ax.set_yticks([i * 10 for i in range(11)])
-        merge_visualizer = self.first_merge_visualizer(merge_visualizer_str)
-        temp_cov = self.temporal_covariate_for_fit is AnomalyTemperatureWithSplineTemporalCovariate
-        merge_visualizer.plot_name = 'Sensitivity plot with ' \
-                                     'shift={} temp_interval={}, temp_cov={}'.format(self.is_shift_interval,
-                                                                                     self.is_temperature_interval,
-                                                                                     temp_cov)
-        merge_visualizer.show_or_save_to_file(no_title=True)
-        plt.close()
+        self.save_plot(merge_visualizer_str, "percentages")
 
-    def interval_plot(self, ax, altitude_class, merge_visualizer_str):
+    def interval_plot_percentages(self, ax, altitude_class, merge_visualizer_str):
         linestyle = get_linestyle_for_altitude_class(altitude_class)
         increasing_key = 'increasing'
         decreasing_key = 'decreasing'
@@ -130,16 +160,19 @@ class VisualizerForSensivity(object):
             color = label_to_color[label]
             ax.plot(self.right_limits, l, label=label_improved, color=color, linestyle=linestyle)
 
+    # Merge visualizer
+
     def get_merge_visualizer(self, altitude_class, visualizer_projection: VisualizerForProjectionEnsemble,
                              merge_visualizer_str):
         if merge_visualizer_str in [AbstractEnsembleFit.Median_merge, AbstractEnsembleFit.Mean_merge]:
             independent_ensemble_fit = \
-            visualizer_projection.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class][
-                IndependentEnsembleFit]
+                visualizer_projection.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class][
+                    IndependentEnsembleFit]
             merge_visualizer = independent_ensemble_fit.merge_function_name_to_visualizer[merge_visualizer_str]
         else:
             together_ensemble_fit = \
-            visualizer_projection.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class][TogetherEnsembleFit]
+                visualizer_projection.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class][
+                    TogetherEnsembleFit]
             merge_visualizer = together_ensemble_fit.visualizer
         merge_visualizer.studies.study.gcm_rcm_couple = (merge_visualizer_str, "merge")
         return merge_visualizer
@@ -148,3 +181,15 @@ class VisualizerForSensivity(object):
         altitude_class = get_altitude_class_from_altitudes(self.altitudes_list[0])
         visualizer_projection = list(self.right_limit_to_visualizer.values())[0]
         return self.get_merge_visualizer(altitude_class, visualizer_projection, merge_visualizer_str)
+
+    # Utils
+
+    def save_plot(self, merge_visualizer_str, name):
+        merge_visualizer = self.first_merge_visualizer(merge_visualizer_str)
+        temp_cov = self.temporal_covariate_for_fit is AnomalyTemperatureWithSplineTemporalCovariate
+        merge_visualizer.plot_name = 'Sensitivity plot for {} with  ' \
+                                     'shift={} temp_interval={}, temp_cov={}'.format(name, self.is_shift_interval,
+                                                                                     self.is_temperature_interval,
+                                                                                     temp_cov)
+        merge_visualizer.show_or_save_to_file(no_title=True)
+        plt.close()

projects/projected_extreme_snowfall/discussion/main_sensitivity.py

@@ -44,21 +44,24 @@ def main():
     set_seed_for_test()
     AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
 
-    fast = False
-    scenarios = rcp_scenarios[-1:] if fast is False else [AdamontScenario.rcp85]
+    fast = None
+    scenarios = [AdamontScenario.rcp85]
+    scenarios = rcp_scenarios[1:]
 
     for scenario in scenarios:
         gcm_rcm_couples = get_gcm_rcm_couples(scenario)
         if fast is None:
+            scenarios = scenarios[:1]
             massif_names = None
             gcm_rcm_couples = gcm_rcm_couples[4:6]
             AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
-            altitudes_list = altitudes_for_groups[3:]
+            altitudes_list = altitudes_for_groups[1:3]
         elif fast:
+            scenarios = scenarios[:1]
             massif_names = ['Vanoise', 'Haute-Maurienne']
             gcm_rcm_couples = gcm_rcm_couples[4:6]
             AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
-            altitudes_list = altitudes_for_groups[:1]
+            altitudes_list = altitudes_for_groups[1:3]
         else:
             massif_names = None
             altitudes_list = altitudes_for_groups[:]
@@ -72,11 +75,11 @@ def main():
         model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
         assert scenario in rcp_scenarios
         remove_physically_implausible_models = True
-        temp_cov = False
+        temp_cov = True
         temporal_covariate_for_fit = AnomalyTemperatureWithSplineTemporalCovariate if temp_cov else TimeTemporalCovariate
         print('Covariate is {}'.format(temporal_covariate_for_fit))
 
-        for is_temperature_interval in [True, False][1:]:
+        for is_temperature_interval in [True, False][:1]:
             for is_shift_interval in [True, False][1:]:
                 visualizer = VisualizerForSensivity(
                     altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,

projects/projected_extreme_snowfall/results/main_projections_ensemble.py

@@ -45,7 +45,7 @@ def main():
     set_seed_for_test()
     AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
 
-    fast = False
+    fast = True
     scenarios = rcp_scenarios[::-1] if fast is False else [AdamontScenario.rcp85]
 
     for scenario in scenarios: