[projection swe] add sensitivity for return levels.

extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py

@@ -69,7 +69,7 @@ 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:]):
+        for i, scenario in enumerate(rcp_scenarios[1:]):
             first_scenario = i == 0
             plot_nb_data_one_line(ax, gcm, scenario, left_limit, right_limit,
                                   first_scenario, is_temperature_interval)
@@ -93,7 +93,7 @@ def plot_nb_data(is_temperature_interval, is_shift_interval):
 def get_interval_limits(is_temperature_interval, is_shift_interval):
     if is_temperature_interval:
         temp_min = np.arange(0, 2, 0.5)
-        temp_max = temp_min + 1.5
+        temp_max = temp_min + 2
         left_limit, right_limit = temp_min, temp_max
     else:
         shift = 25

extreme_trend/ensemble_fit/visualizer_for_sensitivity.py

@@ -59,6 +59,7 @@ class VisualizerForSensivity(object):
                 if year_min_and_year_max[0] is not None:
                     gcm_to_year_min_and_year_max[gcm] = year_min_and_year_max
 
+            print(gcm_to_year_min_and_year_max)
             visualizer = VisualizerForProjectionEnsemble(
                 altitudes_list, gcm_rcm_couples, study_class, season, scenario,
                 model_classes=model_classes,
@@ -84,9 +85,54 @@ class VisualizerForSensivity(object):
             merge_visualizer_str_list.append(AbstractEnsembleFit.Together_merge)
         for merge_visualizer_str in merge_visualizer_str_list:
             self.sensitivity_plot_percentages(merge_visualizer_str)
+            self.sensitivity_plot_return_levels(merge_visualizer_str)
             for relative in [True, False]:
                 self.sensitivity_plot_changes(merge_visualizer_str, relative)
 
+    def sensitivity_plot_return_levels(self, merge_visualizer_str):
+        ax = plt.gca()
+        for altitudes in self.altitudes_list:
+            altitude_class = get_altitude_class_from_altitudes(altitudes)
+            self.interval_plot_return_levels(ax, altitude_class, merge_visualizer_str)
+
+        ticks_labels = get_ticks_labels_for_interval(self.is_temperature_interval, self.is_shift_interval)
+        name = 'Return levels at the end of the interval'
+        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_return_levels(self, ax, altitude_class, merge_visualizer_str):
+        linestyle = get_linestyle_for_altitude_class(altitude_class)
+
+        mean_return_levels = []
+        for v in self.right_limit_to_visualizer.values():
+            merge_visualizer = self.get_merge_visualizer(altitude_class, v, merge_visualizer_str)
+            mean_return_level = merge_visualizer.mean_return_level(self.massif_names)
+            mean_return_levels.append(mean_return_level)
+
+        label = altitude_class().formula
+        ax.plot(self.right_limits, mean_return_levels, linestyle=linestyle, label=label, color='orange')
+
+    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_changes(self, merge_visualizer_str, relative):
         ax = plt.gca()
         for altitudes in self.altitudes_list:

extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py

@@ -266,7 +266,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         moment = moment.replace('moment',
                                 '{} in {}'.format(OneFoldFit.get_moment_str(order=order), OneFoldFit.last_year))
         plot_name = 'Model {} annual maxima of {}'.format(moment,
-                                                            SCM_STUDY_CLASS_TO_ABBREVIATION[self.studies.study_class])
+                                                          SCM_STUDY_CLASS_TO_ABBREVIATION[self.studies.study_class])
         ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
         ax.set_xlabel('altitudes', fontsize=15)
         ax.tick_params(axis='both', which='major', labelsize=13)
@@ -323,9 +323,9 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         raise NotImplementedError
         self.plot_map(massif_name_to_value=massif_name_to_best_coef,
                       label='Coef/{} plot for {} {}'.format(coef_name,
-                                                              SCM_STUDY_CLASS_TO_ABBREVIATION[
-                                                                  type(self.study)],
-                                                              self.study.variable_unit),
+                                                            SCM_STUDY_CLASS_TO_ABBREVIATION[
+                                                                type(self.study)],
+                                                            self.study.variable_unit),
                       add_x_label=False, graduation=graduation, massif_name_to_text=self.massif_name_to_best_name,
                       negative_and_positive_values=negative_and_positive_values)
 
@@ -388,7 +388,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
                                                        massif_name.replace('_', ' '), self.study.variable_unit))
             peak_year_folder = 'Peak year ' + ylabel
             plot_name = '{}/Peak year for {}'.format(peak_year_folder,
-                                                       massif_name.replace('_', ''))
+                                                     massif_name.replace('_', ''))
             self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show, no_title=True, tight_layout=True)
             plt.close()
 
@@ -462,6 +462,14 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
         plt.close()
 
+    def mean_return_level(self, massif_names):
+        valid_massif_names = self.get_valid_names(massif_names)
+        return_levels = []
+        for one_fold in [one_fold for m, one_fold in self.massif_name_to_one_fold_fit.items()
+                         if m in valid_massif_names]:
+            return_levels.append(one_fold.return_level_last_temporal_coordinate)
+        return np.mean(return_levels)
+
     def all_trends(self, massif_names, with_significance=True, with_relative_change=False):
         """return percents which contain decrease, significant decrease, increase, significant increase percentages"""
         valid_massif_names = self.get_valid_names(massif_names)

extreme_trend/one_fold_fit/one_fold_fit.py

@@ -29,6 +29,7 @@ from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_retur
     EurocodeConfidenceIntervalFromExtremes
 from extreme_trend.one_fold_fit.altitude_group import DefaultAltitudeGroup, altitudes_for_groups
 from root_utils import NB_CORES, batch
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import TimeTemporalCovariate
 from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
     AnomalyTemperatureWithSplineTemporalCovariate
@@ -382,6 +383,16 @@ class OneFoldFit(object):
 
         return is_significant, altitude_and_year_to_return_level_mean_estimate, altitude_and_year_to_return_level_confidence_interval
 
+    @property
+    def return_level_last_temporal_coordinate(self):
+        df_temporal_covariate = self.dataset.coordinates.df_temporal_coordinates_for_fit(temporal_covariate_for_fit=self.temporal_covariate_for_fit,
+                                                                                         drop_duplicates=False)
+        last_temporal_coordinate = df_temporal_covariate.loc[:, AbstractCoordinates.COORDINATE_T].max()
+        print('last temporal coordinate', last_temporal_coordinate)
+        altitude = self.altitude_group.reference_altitude
+        coordinate = np.array([altitude, last_temporal_coordinate])
+        return self.get_return_level(self.best_function_from_fit, coordinate)
+
     @property
     def bootstrap_fitted_functions_from_fit(self):
         print('nb of bootstrap for confidence interval=', AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP)

projects/projected_extreme_snowfall/discussion/main_sensitivity.py

@@ -25,7 +25,7 @@ from extreme_fit.model.utils import set_seed_for_test
 
 from extreme_data.meteo_france_data.adamont_data.adamont.adamont_safran import AdamontSnowfall
 from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \
-    rcp_scenarios
+    rcp_scenarios, rcm_scenarios_extended
 from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
     TimeTemporalCovariate
 
@@ -44,24 +44,26 @@ def main():
     set_seed_for_test()
     AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
 
-    fast = None
+    fast = True
     scenarios = [AdamontScenario.rcp85]
+    scenarios = rcm_scenarios_extended[1:]
     scenarios = rcp_scenarios[1:]
 
+    if fast in [None, True]:
+        scenarios = 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[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:3]
+            altitudes_list = altitudes_for_groups[2:3]
         else:
             massif_names = None
             altitudes_list = altitudes_for_groups[:]
@@ -73,7 +75,7 @@ def main():
         print('Scenario is', scenario)
 
         model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
-        assert scenario in rcp_scenarios
+        assert scenario is not AdamontScenario.histo
         remove_physically_implausible_models = True
         temp_cov = True
         temporal_covariate_for_fit = AnomalyTemperatureWithSplineTemporalCovariate if temp_cov else TimeTemporalCovariate