[projection snowfall] add comparison_window_for_the_max.py

extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py

 def get_year_min_and_year_max_used_to_compute_quantile(gcm):
+    """pour les périodes de calcul de quantiles/quantiles, voici le recap :
+    1980-2011 pour la réanalyse
+    1974-2005 pour le modèle
+
+    sauf GCM=MOHC-HadGEM2-ES
+    1987-2011 pour la réanalyse
+    1981-2005 pour le modèle
+    """
     if gcm == 'HadGEM2-ES':
         reanalysis_years = (1988, 2011)
         model_year = (1982, 2005)

extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py

@@ -96,9 +96,9 @@ def get_interval_limits(is_temperature_interval, is_shift_interval):
         left_limit, right_limit = temp_min, temp_max
     else:
         shift = 25
-        nb = 3
-        year_min = [1959 + shift * i for i in range(nb)]
-        year_max = [2020 + shift * i for i in range(nb)]
+        nb = 1
+        year_min = [1950 + shift * i for i in range(nb)]
+        year_max = [2100 + shift * i for i in range(nb)]
         left_limit, right_limit = year_min, year_max
     if not is_shift_interval:
         min_interval_left = min(left_limit)
@@ -122,8 +122,52 @@ def get_ticks_labels_for_interval(is_temperature_interval, is_shift_interval):
     return ticks_labels
 
 
-if __name__ == '__main__':
+def plot_nb_of_data():
     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)
+
+
+def plot_year_for_temperature_interval():
+    left_limit, right_limit = get_interval_limits(True, False)
+
+    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_year_for_temperature_interval_one_line(ax, gcm, scenario, left_limit, right_limit, first_scenario)
+
+    ax.legend(loc='upper left')
+    ticks_labels = get_ticks_labels_for_interval(True, False)
+    ax.set_yticks(right_limit)
+    ax.set_yticklabels(ticks_labels)
+    ax.set_xlabel('Last year considered')
+    ax2 = ax.twinx()
+    legend_elements = [
+        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='lower right')
+    ax2.set_yticks([])
+    plt.show()
+
+
+def plot_year_for_temperature_interval_one_line(ax, gcm, scenario, left_limits, right_limits, first_scenario):
+    year_max_list = [get_year_min_and_year_max(gcm, scenario, left_limit, right_limit,
+                              True)[1] for left_limit, right_limit in zip(left_limits, right_limits)]
+
+    color = gcm_to_color[gcm]
+    linestyle = get_linestyle_from_scenario(scenario)
+
+    # For the legend
+    if (len(year_max_list) > 0) and first_scenario:
+        ax.plot(year_max_list[0], right_limits[0], linestyle='solid', color=color, label=gcm)
+
+    ax.plot(year_max_list, right_limits, linestyle=linestyle, color=color, marker='o')
+
+
+if __name__ == '__main__':
+    # plot_nb_of_data()
+    plot_year_for_temperature_interval()
+

extreme_trend/ensemble_fit/visualizer_for_sensitivity.py

@@ -48,7 +48,7 @@ class VisualizerForSensivity(object):
 
         self.left_limit_to_right_limit = OrderedDict(zip(self.left_limits, self.right_limits))
         self.right_limit_to_visualizer = {}  # type: Dict[float, VisualizerForProjectionEnsemble]
-
+        print(self.right_limits)
         for left_limit, right_limit in zip(self.left_limits, self.right_limits):
             interval_str_prefix = "{}-{}".format(left_limit, right_limit)
             interval_str_prefix = interval_str_prefix.replace('.', ',')
@@ -88,18 +88,24 @@ 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 only_mean in [True, False]:
+                self.sensitivity_plot_return_levels(merge_visualizer_str, only_mean)
             for relative in [True, False]:
                 self.sensitivity_plot_changes(merge_visualizer_str, relative)
 
-    def sensitivity_plot_return_levels(self, merge_visualizer_str):
+    def sensitivity_plot_return_levels(self, merge_visualizer_str, only_mean):
         ax = plt.gca()
         for altitudes in self.altitudes_list:
             altitude_group = get_altitude_class_from_altitudes(altitudes)
-            self.interval_plot_return_levels(ax, altitude_group, merge_visualizer_str)
+            self.interval_plot_return_levels(ax, altitude_group, merge_visualizer_str, only_mean)
+        # todo: dans la fonction return level de sensitivity visualizer faire en sorte de ne pas tracer que le retour de niveau 100,
+        #  mais en tracer d autres pour voir comment ils évoluent.
 
         ticks_labels = get_ticks_labels_for_interval(self.is_temperature_interval, self.is_shift_interval)
-        name = 'Return levels at the end of the interval'
+        if only_mean:
+            name = 'Mean return levels at the end of the interval'
+        else:
+            name = 'Return levels at the end of the interval'
         if len(self.altitudes_list) == 1:
             altitude_group = get_altitude_group_from_altitudes(self.altitudes_list[0])
             name += ' at {} m'.format(altitude_group.reference_altitude)
@@ -109,7 +115,7 @@ class VisualizerForSensivity(object):
         ax.set_xticklabels(ticks_labels)
         self.save_plot(merge_visualizer_str, name)
 
-    def interval_plot_return_levels(self, ax, altitude_class, merge_visualizer_str):
+    def interval_plot_return_levels(self, ax, altitude_class, merge_visualizer_str, only_mean):
 
         mean_return_levels = []
         massif_names = AbstractStudy.all_massif_names() if self.massif_names is None else self.massif_names
@@ -121,18 +127,21 @@ class VisualizerForSensivity(object):
                 massif_name_to_tuple_list[massif_name].append((r, return_level))
             mean_return_levels.append(np.mean(list(massif_name_to_return_level.values())))
 
-        for massif_id, massif_name in enumerate(AbstractStudy.all_massif_names()):
-            if len(massif_name_to_tuple_list[massif_name]) > 0:
-                x, y = zip(*massif_name_to_tuple_list[massif_name])
-            else:
-                x, y = [], []
-            color, linestyle, label = get_color_and_linestyle_from_massif_id(massif_id, massif_name)
-            ax.plot(x, y, color=color, linestyle=linestyle, label=label)
+        if not only_mean:
+            for massif_id, massif_name in enumerate(AbstractStudy.all_massif_names()):
+                if (massif_name in massif_name_to_tuple_list) and (len(massif_name_to_tuple_list[massif_name]) > 0):
+                    x, y = zip(*massif_name_to_tuple_list[massif_name])
+                else:
+                    x, y = [], []
+                color, linestyle, label = get_color_and_linestyle_from_massif_id(massif_id, massif_name)
+                ax.plot(x, y, color=color, linestyle=linestyle, label=label)
 
         ax.plot(self.right_limits, mean_return_levels, label="Mean", color='k', linewidth=4)
-        down, up = ax.get_ylim()
-        ax.set_ylim((down, up * 1.3))
-        ax.legend(prop={'size': 7}, ncol=3, loc="upper center")
+
+        if not only_mean:
+            down, up = ax.get_ylim()
+            ax.set_ylim((down, up * 1.3))
+            ax.legend(prop={'size': 7}, ncol=3, loc="upper center")
 
     def interval_plot_changes(self, ax, altitude_class, merge_visualizer_str, relative):
         linestyle = get_linestyle_for_altitude_class(altitude_class)

extreme_trend/one_fold_fit/one_fold_fit.py

@@ -388,6 +388,7 @@ class OneFoldFit(object):
         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()
+        # todo: améliorer the last temporal coordinate. on recupère la liste des rights_limits, puis on prend la valeur juste au dessus ou égale."""
         print('last temporal coordinate', last_temporal_coordinate)
         altitude = self.altitude_group.reference_altitude
         coordinate = np.array([altitude, last_temporal_coordinate])

projects/past_extreme_snowfall/section_data_and_results/main_altitudes_studies.py

@@ -41,9 +41,10 @@ def main():
     study_classes = [SafranSnowfall1Day
                      , SafranSnowfall3Days,
                      SafranSnowfall5Days, SafranSnowfall7Days][:1]
-    study_classes = [SafranSnowfall2019, SafranSnowfall2020, SafranSnowfallCenterOnDay1day,
+    study_classes = [SafranSnowfall2020, SafranSnowfall2019, SafranSnowfallCenterOnDay1day,
                      SafranSnowfallNotCenterOnDay1day,
-                     SafranSnowfallCenterOnDay1dayMeanRate, SafranSnowfall1Day][1:2]
+                     SafranSnowfallCenterOnDay1dayMeanRate, SafranSnowfall1Day][1:3]
+    study_classes = [SafranSnowfallNotCenterOnDay1day, SafranSnowfall2019]
     seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
 
     set_seed_for_test()
@@ -88,10 +89,10 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes, model_must_p
 
 def plot_visualizers(massif_names, visualizer_list):
     # plot_histogram_all_models_against_altitudes(massif_names, visualizer_list)
-    plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list, with_significance=False)
+    plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list, with_significance=True)
     # plot_shoe_plot_ratio_interval_size_against_altitude(massif_names, visualizer_list)
     for relative in [True, False]:
-        plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative, with_significance=False)
+        plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative, with_significance=True)
     # plot_coherence_curves(['Vanoise'], visualizer_list)
     pass
 

projects/projected_extreme_snowfall/comparison_window_for_the_max/comparison_window_for_the_max.py

+import matplotlib.pyplot as plt
+from extreme_data.meteo_france_data.scm_models_data.altitudes_studies import AltitudesStudies
+from extreme_data.meteo_france_data.scm_models_data.safran.gap_between_study import \
+    GapBetweenSafranSnowfall2019AndMySafranSnowfall2019
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfallCenterOnDay, \
+    SafranSnowfallCenterOnDay1day
+from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
+    STUDY_CLASS_TO_ABBREVIATION
+
+
+def plot_comparison(altitudes, massif_name, stud_class_list):
+    altitude_studies_list = [AltitudesStudies(study_class, altitudes) for study_class in stud_class_list]
+
+    ax = plt.gca()
+
+    for altitude_studies in altitude_studies_list:
+        plot_studies(ax, altitude_studies, massif_name)
+
+    labelsize = 10
+
+    ax.tick_params(axis='both', which='major', labelsize=labelsize)
+    handles, labels = ax.get_legend_handles_labels()
+    ax.legend(handles[::-1], labels[::-1], prop={'size': labelsize})
+
+    plot_name = 'comparaison for {} at {}'.format(massif_name, altitudes[0])
+    altitude_studies_list[0].show_or_save_to_file(plot_name=plot_name, show=False, no_title=True, tight_layout=True)
+    ax.clear()
+    plt.close()
+
+def plot_studies(ax, altitude_studies, massif_name):
+    for altitude, study in list(altitude_studies.altitude_to_study.items()):
+        x = study.ordered_years
+        if massif_name in study.massif_name_to_annual_maxima:
+            y = study.massif_name_to_annual_maxima[massif_name]
+            label = '{} m'.format(altitude)
+            ax.plot(x, y, linewidth=2, label=label)
+        study_class = type(study)
+        plot_name = 'Annual maxima of {} in {}'.format(STUDY_CLASS_TO_ABBREVIATION[study_class],
+                                                       massif_name.replace('_', ' '))
+        fontsize = 10
+        ax.set_ylabel('{} ({})'.format(plot_name, study.variable_unit), fontsize=fontsize)
+
+
+if __name__ == '__main__':
+    for study_class in [GapBetweenSafranSnowfall2019AndMySafranSnowfall2019, SafranSnowfall1Day, SafranSnowfallCenterOnDay1day][:1]:
+        for altitudes, massif_name in [([2100, 2400, 2700], 'Chablais'), ([1200, 1500, 1800], 'Vercors')][:]:
+            plot_comparison(altitudes, massif_name, [study_class])
\ No newline at end of file

projects/projected_extreme_snowfall/discussion/main_sensitivity.py

@@ -45,7 +45,7 @@ def main():
     set_seed_for_test()
     AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
 
-    fast = None
+    fast = False
     scenarios = [AdamontScenario.rcp85]
     scenarios = rcp_scenarios[1:]
     scenarios = rcm_scenarios_extended[2:][::-1]
@@ -59,7 +59,7 @@ def main():
             massif_names = None
             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]
         elif fast:
             massif_names = ['Vanoise', 'Haute-Maurienne']
             gcm_rcm_couples = gcm_rcm_couples[4:6]
@@ -78,11 +78,11 @@ def main():
         model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
         assert scenario is not AdamontScenario.histo
         remove_physically_implausible_models = True
-        temp_cov = True
+        temp_cov = False
         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/discussion/main_sensitivity_one_altitude_range.py

@@ -48,7 +48,7 @@ def main():
     fast = False
     scenarios = [AdamontScenario.rcp85]
     scenarios = rcp_scenarios[1:]
-    scenarios = rcm_scenarios_extended[:][::-1]
+    scenarios = rcm_scenarios_extended[2:][::-1]
 
     if fast in [None, True]:
         scenarios = scenarios[:1]
@@ -64,10 +64,10 @@ def main():
             massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors', 'Ubaye']
             gcm_rcm_couples = gcm_rcm_couples[4:6]
             AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
-            altitudes_list = altitudes_for_groups[1:3]
+            altitudes_list = altitudes_for_groups[1:2]
         else:
             massif_names = None
-            altitudes_list = altitudes_for_groups[:]
+            altitudes_list = altitudes_for_groups[2:3]
 
         assert isinstance(gcm_rcm_couples, list)
 

projects/projected_extreme_snowfall/results/main_projections_ensemble.py

@@ -24,7 +24,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
 
@@ -47,6 +47,7 @@ def main():
 
     fast = True
     scenarios = rcp_scenarios[::-1] if fast is False else [AdamontScenario.rcp85]
+    scenarios = rcm_scenarios_extended[::-1]
 
     for scenario in scenarios:
         gcm_rcm_couples = get_gcm_rcm_couples(scenario)
@@ -72,7 +73,7 @@ def main():
         print('Covariate is {}'.format(temporal_covariate_for_fit))
 
         model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
-        assert scenario in rcp_scenarios
+        assert scenario is not AdamontScenario.histo
 
         visualizer = VisualizerForProjectionEnsemble(
             altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,