add scenario rcp85_extended and how to handle it.

extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py

@@ -10,7 +10,8 @@ from netCDF4._netCDF4 import Dataset
 
 from extreme_data.meteo_france_data.adamont_data.adamont.adamont_variables import AbstractAdamontVariable
 from extreme_data.meteo_france_data.adamont_data.adamont_scenario import scenario_to_str, AdamontScenario, \
-    get_year_min_and_year_max_from_scenario, gcm_rcm_couple_to_full_name, get_suffix_for_the_nc_file
+    get_year_min_and_year_max_from_scenario, gcm_rcm_couple_to_full_name, get_suffix_for_the_nc_file, \
+    scenario_to_real_scenarios, get_year_max
 from extreme_data.meteo_france_data.adamont_data.utils.utils import massif_number_to_massif_name
 
 ADAMONT_PATH = r"/home/erwan/Documents/projects/spatiotemporalextremes/local/spatio_temporal_datasets/ADAMONT"
@@ -56,31 +57,37 @@ class AbstractAdamontStudy(AbstractStudy):
     def ordered_years(self):
         return list(range(self.year_min, self.year_max + 1))
 
-    @cached_property
-    def winter_year_for_each_time_step(self):
-        year_min, _ = get_year_min_and_year_max_from_scenario(self.scenario, self.gcm_rcm_couple)
+    def winter_years_for_each_time_step(self, real_scenario, dataset):
+        year_min, year_max = get_year_min_and_year_max_from_scenario(real_scenario, self.gcm_rcm_couple)
         # It was written in the dataset  for the TIME variable that it represents
         # "'hours since 1950-08-01 06:00:00'" for the HISTO scenario
         # "'hours since 2005-08-01 06:00:00'" for the RCP scenario
         start = datetime(year=year_min - 1, month=8, day=1, hour=6, minute=0, second=0)
-        hours_after_start = np.array(self.dataset.variables['TIME'])
+        hours_after_start = np.array(dataset.variables['TIME'])
         dates = [start + timedelta(hours=h) for h in hours_after_start]
-        winter_year = [date.year if date.month < 8 else date.year + 1 for date in dates]
-        # Remark. The last winter year for the HISTO scenario correspond to 2006.
-        # Thus, the last value is not taken into account
-        return np.array(winter_year)
+        winter_years = [date.year if date.month < 8 else date.year + 1 for date in dates]
+        return winter_years
 
     @cached_property
     def year_to_variable_object(self) -> OrderedDict:
         year_to_data_list = {}
         for year in self.ordered_years:
             year_to_data_list[year] = []
-        data_list = self.dataset.variables[self.variable_class.keyword()]
-        data_year_list = self.winter_year_for_each_time_step
-        assert len(data_list) == len(data_year_list)
+        # Load data & year list
+        data_list, data_year_list = [], []
+        for dataset, real_scenario in zip(self.datasets, self.adamont_real_scenarios):
+            data_list.extend(dataset.variables[self.variable_class.keyword()])
+            data_year_list.extend(self.winter_years_for_each_time_step(real_scenario, dataset))
+            # Remark. The last winter year for the HISTO scenario correspond to 2006.
+            # Thus, the last value is not taken into account
+            if data_year_list[-1] > get_year_max(real_scenario, self.gcm_rcm_couple):
+                data_year_list = data_year_list[:-1]
+                data_list = data_list[:-1]
+            assert len(data_list) == len(data_year_list)
         for year_data, data in zip(data_year_list, data_list):
             if self.year_min <= year_data <= self.year_max:
                 year_to_data_list[year_data].append(data)
+        # Load the variable object
         year_to_variable_object = OrderedDict()
         for year in self.ordered_years:
             variable_array = np.array(year_to_data_list[year])
@@ -98,8 +105,13 @@ class AbstractAdamontStudy(AbstractStudy):
         return [massif_number_to_massif_name[massif_id] for massif_id in massif_ids]
 
     @cached_property
+    def datasets(self):
+        return [Dataset(file_path) for file_path in self.nc_file_paths]
+
+    @property
     def dataset(self):
-        return Dataset(self.nc_file_path)
+        # todo: improve that
+        return self.datasets[0]
 
     # PATHS
 
@@ -107,22 +119,29 @@ class AbstractAdamontStudy(AbstractStudy):
     def variable_folder_name(self):
         return self.variable_class.variable_name_for_folder_and_nc_file()
 
-    @property
-    def scenario_name(self):
-        return scenario_to_str(self.scenario)
-
     @property
     def region_name(self):
         return french_region_to_str(self.french_region)
 
     @property
-    def nc_files_path(self):
-        return op.join(ADAMONT_PATH, self.variable_folder_name, self.scenario_name)
+    def adamont_real_scenarios(self):
+        return scenario_to_real_scenarios(self.scenario)
+
+    @property
+    def scenario_names(self):
+        return [scenario_to_str(scenario) for scenario in self.adamont_real_scenarios]
+
+    @property
+    def nc_files_paths(self):
+        return [op.join(ADAMONT_PATH, self.variable_folder_name, name) for name in self.scenario_names]
 
     @property
-    def nc_file_path(self):
-        suffix_nc_file = get_suffix_for_the_nc_file(self.scenario, self.gcm_rcm_couple)
-        nc_file = '{}_FORCING_{}_{}_{}_{}.nc'.format(self.variable_folder_name, self.gcm_rcm_full_name,
-                                                     self.scenario_name,
-                                                     self.region_name, suffix_nc_file)
-        return op.join(self.nc_files_path, nc_file)
+    def nc_file_paths(self):
+        file_paths = []
+        for scenario, scenario_name, files_path in zip(self.adamont_real_scenarios, self.scenario_names, self.nc_files_paths):
+            suffix_nc_file = get_suffix_for_the_nc_file(scenario, self.gcm_rcm_couple)
+            nc_file = '{}_FORCING_{}_{}_{}_{}.nc'.format(self.variable_folder_name, self.gcm_rcm_full_name,
+                                                         scenario_name,
+                                                         self.region_name, suffix_nc_file)
+            file_paths.append(op.join(files_path, nc_file))
+        return file_paths

extreme_data/meteo_france_data/adamont_data/adamont_scenario.py

@@ -6,12 +6,36 @@ class AdamontScenario(Enum):
     rcp26 = 1
     rcp45 = 2
     rcp85 = 3
+    rcp85_extended = 4
+
+
+adamont_scenarios_real = [AdamontScenario.histo, AdamontScenario.rcp26, AdamontScenario.rcp45, AdamontScenario.rcp85]
 
 
 def get_year_min_and_year_max_from_scenario(adamont_scenario, gcm_rcm_couple):
     assert isinstance(adamont_scenario, AdamontScenario)
+    year_min = get_year_min(adamont_scenario, gcm_rcm_couple)
+    year_max = get_year_max(adamont_scenario, gcm_rcm_couple)
+    return year_min, year_max
+
+
+def get_year_max(adamont_scenario, gcm_rcm_couple):
+    real_adamont_scenarios = scenario_to_real_scenarios(adamont_scenario)
     gcm, rcm = gcm_rcm_couple
-    if adamont_scenario == AdamontScenario.histo:
+    if any([scenario is not AdamontScenario.histo for scenario in real_adamont_scenarios]):
+        if gcm == 'HadGEM2-ES':
+            year_max = 2099
+        else:
+            year_max = 2100
+    else:
+        year_max = 2005
+    return year_max
+
+
+def get_year_min(adamont_scenario, gcm_rcm_couple):
+    real_adamont_scenarios = scenario_to_real_scenarios(adamont_scenario)
+    gcm, rcm = gcm_rcm_couple
+    if AdamontScenario.histo in real_adamont_scenarios:
         if gcm == 'HadGEM2-ES':
             year_min = 1982
         elif rcm == 'RCA4':
@@ -20,14 +44,9 @@ def get_year_min_and_year_max_from_scenario(adamont_scenario, gcm_rcm_couple):
             year_min = 1952
         else:
             year_min = 1951
-        year_max = 2005
     else:
         year_min = 2006
-        if gcm == 'HadGEM2-ES':
-            year_max = 2099
-        else:
-            year_max = 2100
-    return year_min, year_max
+    return year_min
 
 
 def load_gcm_rcm_couples_for_year_min_and_year_max(year_min, year_max, adamont_scenario=AdamontScenario.histo):
@@ -42,6 +61,7 @@ def load_gcm_rcm_couples_for_year_min_and_year_max(year_min, year_max, adamont_s
 
 
 def get_suffix_for_the_nc_file(adamont_scenario, gcm_rcm_couple):
+    assert isinstance(adamont_scenario, AdamontScenario)
     year_min, year_max = get_year_min_and_year_max_from_scenario(adamont_scenario, gcm_rcm_couple)
     return '{}080106_{}080106_daysum'.format(year_min - 1, year_max)
 
@@ -50,6 +70,16 @@ def scenario_to_str(adamont_scenario):
     return str(adamont_scenario).split('.')[-1].upper()
 
 
+def scenario_to_real_scenarios(adamont_scenario):
+    if adamont_scenario in adamont_scenarios_real:
+        return [adamont_scenario]
+    else:
+        if adamont_scenario is AdamontScenario.rcp85_extended:
+            return [AdamontScenario.histo, AdamontScenario.rcp85]
+        else:
+            raise NotImplementedError
+
+
 def gcm_rcm_couple_to_str(gcm_rcm_couple):
     return ' / '.join(gcm_rcm_couple)
 

extreme_data/meteo_france_data/adamont_data/adamont_studies.py

@@ -2,8 +2,10 @@ import matplotlib.pyplot as plt
 
 from collections import OrderedDict
 
+import numpy as np
 from cached_property import cached_property
 
+from extreme_data.meteo_france_data.adamont_data.abstract_adamont_study import AbstractAdamontStudy
 from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_full_name, \
     gcm_rcm_couple_to_str, get_color_from_gcm_rcm_couple
 from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
@@ -19,7 +21,7 @@ class AdamontStudies(object):
         if gcm_rcm_couples is None:
             gcm_rcm_couples = list(gcm_rcm_couple_to_full_name.keys())
         self.gcm_rcm_couples = gcm_rcm_couples
-        self.gcm_rcm_couple_to_study = OrderedDict()  # type: OrderedDict[int, AbstractStudy]
+        self.gcm_rcm_couple_to_study = OrderedDict()  # type: OrderedDict[int, AbstractAdamontStudy]
         for gcm_rcm_couple in self.gcm_rcm_couples:
             study = study_class(gcm_rcm_couple=gcm_rcm_couple, **kwargs_study)
             self.gcm_rcm_couple_to_study[gcm_rcm_couple] = study
@@ -55,7 +57,16 @@ class AdamontStudies(object):
                 label = gcm_rcm_couple_to_str(gcm_rcm_couple)
                 color = get_color_from_gcm_rcm_couple(gcm_rcm_couple)
                 ax.plot(x, y, linewidth=linewidth, label=label, color=color)
-        if scm_study is not None:
+        if scm_study is None:
+            y = np.array([study.massif_name_to_annual_maxima[massif_name] for study in self.study_list
+                 if massif_name in study.massif_name_to_annual_maxima])
+            if len(y) > 0:
+                y = np.mean(y, axis=0)
+                label = 'Mean maxima'
+                color = 'black'
+                ax.plot(x, y, linewidth=linewidth * 2, label=label, color=color)
+        else:
+            # todo: otherwise display the mean in strong black
             try:
                 y = scm_study.massif_name_to_annual_maxima[massif_name]
                 label = 'Reanalysis'
@@ -64,10 +75,11 @@ class AdamontStudies(object):
             except KeyError:
                 pass
 
-        ax.xaxis.set_ticks(x[1::10])
+        ax.xaxis.set_ticks([year for year in y if year % 10 == 0])
+        ax.grid()
         ax.tick_params(axis='both', which='major', labelsize=13)
         handles, labels = ax.get_legend_handles_labels()
-        ax.legend(handles[::-1], labels[::-1])
+        ax.legend(handles[::-1], labels[::-1], ncol=2)
         plot_name = 'Annual maxima of {} in {}'.format(ADAMONT_STUDY_CLASS_TO_ABBREVIATION[self.study_class],
                                                        massif_name.replace('_', ' '))
         ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)

projects/projected_snowfall/comparison_with_scm/main_comparison_on_historical_period.py

+from collections import OrderedDict
+
 import matplotlib as mpl
 mpl.use('Agg')
 mpl.rcParams['text.usetex'] = True
 mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
 
-from projects.projected_snowfall.comparison_with_scm.comparison_plot import individual_plot
+from projects.projected_snowfall.comparison_with_scm.comparison_plot import individual_plot, collective_plot
 
 from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
 from extreme_data.meteo_france_data.adamont_data.adamont_scenario import load_gcm_rcm_couples_for_year_min_and_year_max, \
@@ -34,15 +36,15 @@ def main():
     # Load visualizers
 
     # Individual plot
-    for altitude in altitudes:
-        v = load_visualizer_histo(altitude, massif_names, year_min)
-        individual_plot(v)
+    # for altitude in altitudes:
+    #     v = load_visualizer_histo(altitude, massif_names, year_min)
+    #     individual_plot(v)
 
     # # Collective plot
-    # altitude_to_visualizer = OrderedDict()
-    # for altitude in altitudes:
-    #     altitude_to_visualizer[altitude] = load_visualizer_histo(altitude, massif_names, year_min)
-    # collective_plot(altitude_to_visualizer)
+    altitude_to_visualizer = OrderedDict()
+    for altitude in altitudes:
+        altitude_to_visualizer[altitude] = load_visualizer_histo(altitude, massif_names, year_min)
+    collective_plot(altitude_to_visualizer)
 
 
 def load_visualizer_histo(altitude, massif_names, year_min):

projects/projected_snowfall/projected_data/main_projection.py

@@ -18,7 +18,7 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua
 
 
 def main():
-    fast = True
+    fast = False
     # Set the year_min and year_max for the comparison
     if fast:
         year_min = [2006][0]
@@ -29,18 +29,19 @@ def main():
         year_min = [2006][0]
         year_max = [2100][0]
         massif_names = None
-        altitudes = [900, 1800, 2700]
+        altitudes = [900, 1800, 2700, 3600][2:]
+    adamont_scenario = AdamontScenario.rcp85
 
     # Load studies
     for altitude in altitudes:
         adamont_study_class = AdamontSnowfall
         season = Season.annual
         gcm_rcm_couples = load_gcm_rcm_couples_for_year_min_and_year_max(year_min, year_max,
-                                                                         adamont_scenario=AdamontScenario.rcp85)
+                                                                         adamont_scenario=adamont_scenario)
         adamont_studies = AdamontStudies(adamont_study_class, gcm_rcm_couples,
                                          altitude=altitude, year_min=year_min,
                                          year_max=year_max, season=season,
-                                         scenario=AdamontScenario.rcp85)
+                                         scenario=adamont_scenario)
         adamont_studies.plot_maxima_time_series(massif_names)
 
 

test/test_extreme_data/test_meteo_france_data/test_adamont_study.py

@@ -11,14 +11,17 @@ class TestAdamontStudy(unittest.TestCase):
             AdamontSnowfall(altitude=900),
             AdamontSnowfall(altitude=1800)
         ]
-        study_list.extend([AdamontSnowfall(altitude=900, scenario=AdamontScenario.rcp45),
-                           AdamontSnowfall(altitude=900, scenario=AdamontScenario.rcp85)])
+        study_list.extend([
+            AdamontSnowfall(altitude=900, scenario=AdamontScenario.rcp45),
+            AdamontSnowfall(altitude=900, scenario=AdamontScenario.rcp85),
+            AdamontSnowfall(altitude=900, scenario=AdamontScenario.rcp85_extended)
+        ])
         study_list.extend([AdamontSnowfall(altitude=900, gcm_rcm_couple=gcm_rcm_couple)
                            for gcm_rcm_couple in gcm_rcm_couple_to_full_name.keys()])
         for study in study_list:
             annual_maxima_for_year_min = study.year_to_annual_maxima[study.year_min]
             # print(study.altitude, study.scenario_name, study.gcm_rcm_couple)
-            # print(annual_maxima_for_year_min)
+            # print(len(study.massif_name_to_annual_maxima['Vanoise']))
         self.assertTrue(True)