try to implement the method for the paper 2 on projected snowfall.

extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py

@@ -21,6 +21,8 @@ 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.utils import Season, FrenchRegion, french_region_to_str
 
+class WrongYearMinOrYearMax(Exception):
+    pass
 
 class AbstractAdamontStudy(AbstractStudy):
     YEAR_MIN = 1950
@@ -33,8 +35,13 @@ class AbstractAdamontStudy(AbstractStudy):
                  scenario=AdamontScenario.histo, gcm_rcm_couple=('CNRM-CM5', 'ALADIN53')):
         # Load the default year_min & year_max for the scenario if not specified
         year_min_scenario, year_max_scenario = get_year_min_and_year_max_from_scenario(scenario, gcm_rcm_couple)
+        # Raise exception
         if year_min is None:
             year_min = year_min_scenario
+        else:
+            if year_min < year_min_scenario:
+                raise WrongYearMinOrYearMax('year min is {} and should be larger than {}'.format(year_min, year_min_scenario))
+
         if year_max is None:
             year_max = year_max_scenario
         super().__init__(variable_class=variable_class, altitude=altitude, year_min=year_min, year_max=year_max,

extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py

@@ -32,6 +32,9 @@ SLEurocode = 'SL from max HS with ' + eurocode_snow_density
 SCM_STUDIES = [SafranSnowfall, CrocusSweTotal, CrocusDepth, CrocusSwe3Days]
 SCM_STUDIES_NAMES = [get_display_name_from_object_type(k) for k in SCM_STUDIES]
 SCM_STUDY_NAME_TO_SCM_STUDY = dict(zip(SCM_STUDIES_NAMES, SCM_STUDIES))
+
+
+# I keep the scm study separated from the adamont study (for the tests)
 SCM_STUDY_CLASS_TO_ABBREVIATION = {
     SafranSnowfall: 'SF3',
     SafranSnowfall1Day: 'daily snowfall',
@@ -62,6 +65,7 @@ SCM_STUDY_CLASS_TO_ABBREVIATION = {
 ADAMONT_STUDY_CLASS_TO_ABBREVIATION = {
     AdamontSnowfall: 'daily snowfall',
 }
+STUDY_CLASS_TO_ABBREVIATION = {**ADAMONT_STUDY_CLASS_TO_ABBREVIATION, **SCM_STUDY_CLASS_TO_ABBREVIATION}
 
 altitude_massif_name_and_study_class_for_poster = [
     (900, 'Chartreuse', CrocusSweTotal),

projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py

@@ -6,7 +6,7 @@ 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.visualization.main_study_visualizer import \
-    SCM_STUDY_CLASS_TO_ABBREVIATION
+    SCM_STUDY_CLASS_TO_ABBREVIATION, STUDY_CLASS_TO_ABBREVIATION
 from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import plot_against_altitude
 from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import StudyVisualizer
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
@@ -134,7 +134,8 @@ class AltitudesStudies(object):
                 y = study.massif_name_to_annual_maxima[massif_name]
                 label = '{} m'.format(altitude)
                 ax.plot(x, y, linewidth=2, label=label)
-        ax.xaxis.set_ticks(x[11::20])
+        ax.xaxis.set_ticks([e for e in x if e % 10 == 0][::2])
+        ax.set_xlim((x[0], x[-1]))
 
         # Plot for the paper 2
         if massif_name == "Vanoise":
@@ -144,10 +145,10 @@ class AltitudesStudies(object):
         ax.tick_params(axis='both', which='major', labelsize=20)
         handles, labels = ax.get_legend_handles_labels()
         ax.legend(handles[::-1], labels[::-1], prop={'size': 20})
-        plot_name = 'Annual maxima of {} in {}'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[self.study_class],
+        plot_name = 'Annual maxima of {} in {}'.format(STUDY_CLASS_TO_ABBREVIATION[self.study_class],
                                                        massif_name.replace('_', ' '))
-        # ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
-        # ax.set_xlabel('years', fontsize=15)
+        ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
+        ax.set_xlabel('years', fontsize=15)
         plot_name = 'time series/' + plot_name
         self.show_or_save_to_file(plot_name=plot_name, show=show, no_title=True, tight_layout=True)
         ax.clear()

projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py

@@ -5,24 +5,17 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_s
     AltitudesStudiesVisualizerForNonStationaryModels
 
 
-def load_visualizer_list(season, study_class, altitudes_list, massif_names):
+def load_visualizer_list(season, study_class, altitudes_list, massif_names, **kwargs_study):
     model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
     visualizer_list = []
     # Load all studies
     for altitudes in altitudes_list:
-        # if issubclass(study_class, SimulationStudy):
-        #     for ensemble_idx in list(range(14))[:1]:
-        #         studies = AltitudesStudies(study_class, altitudes, season=season,
-        #                                    ensemble_idx=ensemble_idx)
-        #         plot_studies(massif_names, season, studies, study_class)
-        # else:
-        studies = AltitudesStudies(study_class, altitudes, season=season)
+        studies = AltitudesStudies(study_class, altitudes, season=season, **kwargs_study)
         visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
                                                                       model_classes=model_classes,
                                                                       massif_names=massif_names,
                                                                       show=False,
                                                                       temporal_covariate_for_fit=None,
-                                                                      # temporal_covariate_for_fit=MeanAlpsTemperatureCovariate,
                                                                       )
         visualizer_list.append(visualizer)
     compute_and_assign_max_abs(visualizer_list)

projects/projected_snowfall/trends/main_projected_trends_separately.py

+import time
+from typing import List
+
+import matplotlib as mpl
+import numpy as np
+
+from extreme_data.meteo_france_data.adamont_data.abstract_adamont_study import WrongYearMinOrYearMax
+from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoad3Days
+from projects.altitude_spatial_model.altitudes_fit.plots.plot_coherence_curves import plot_coherence_curves
+from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude_studies import \
+    plot_histogram_all_models_against_altitudes, plot_histogram_all_trends_against_altitudes, \
+    plot_shoe_plot_changes_against_altitude
+
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+
+from projects.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import load_visualizer_list
+
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import altitudes_for_groups
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_individual_aic
+
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \
+    SafranSnowfall5Days, SafranSnowfall7Days, SafranDateFirstSnowfall, SafranPrecipitation1Day, SafranPrecipitation3Days
+from extreme_data.meteo_france_data.scm_models_data.utils import Season
+
+
+def main():
+    study_class = AdamontSnowfall
+    scenario = AdamontScenario.rcp85_extended
+    gcm_rcm_couples = [('CNRM-CM5', 'CCLM4-8-17'), ('CNRM-CM5', 'RCA4')][1:]
+
+    fast = True
+    if fast is None:
+        massif_names = None
+        altitudes_list = altitudes_for_groups[1:2]
+    elif fast:
+        massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:1]
+        altitudes_list = altitudes_for_groups[1:3]
+    else:
+        massif_names = None
+        altitudes_list = altitudes_for_groups[:]
+
+    # One plot for each massif name
+    for massif_name in massif_names:
+        print(massif_name)
+        main_loop(altitudes_list, [massif_name], gcm_rcm_couples, study_class, scenario)
+
+
+def main_loop(altitudes_list, massif_names, gcm_rcm_couples, study_class, scenario):
+    assert isinstance(altitudes_list, List)
+    assert isinstance(altitudes_list[0], List)
+    altitudes_for_return_levels = [altitudes[-2] for altitudes in altitudes_list]
+    print(altitudes_for_return_levels)
+    season = Season.annual
+    first_year_min, first_year_max = 1951, 2010
+    nb_years = first_year_max - first_year_min + 1
+    temporal_windows = [(first_year_min + i, first_year_max + i) for i in [30 * j for j in range(4)]]
+    all_changes_in_return_levels = []
+    for gcm_rcm_couple in gcm_rcm_couples:
+        print('Inner', gcm_rcm_couple)
+        changes_in_return_levels = []
+        for temporal_window in temporal_windows:
+            year_min, year_max = temporal_window
+            try:
+                visualizer_list = load_visualizer_list(season, study_class, altitudes_list, massif_names,
+                                                       scenario=scenario, gcm_rcm_couple=gcm_rcm_couple,
+                                                       year_min=year_min, year_max=year_max)
+
+                for visualizer in visualizer_list:
+                    visualizer.studies.plot_maxima_time_series(massif_names)
+                massif_name = massif_names[0]
+                changes_for_temporal_window = [
+                    v.massif_name_to_one_fold_fit[massif_name].changes_of_moment(altitudes=[a],
+                                                                                 year=year_max,
+                                                                                 nb_years=nb_years,
+                                                                                 order=None
+                                                                                 )[0]
+                    for a, v in zip(altitudes_for_return_levels, visualizer_list)]
+
+            except WrongYearMinOrYearMax:
+                changes_for_temporal_window = [np.nan for _ in altitudes_for_return_levels]
+
+            print(changes_for_temporal_window)
+            changes_in_return_levels.append(changes_for_temporal_window)
+        changes_in_return_levels = np.array(changes_in_return_levels)
+        all_changes_in_return_levels.append(changes_in_return_levels)
+    all_changes_in_return_levels = np.array(all_changes_in_return_levels)
+
+    return all_changes_in_return_levels, temporal_windows, altitudes_for_return_levels, nb_years
+
+
+if __name__ == '__main__':
+    main()