From a3dd9ba5e13000fa0b3def0d8a5e28fc0db42385 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Fri, 12 Feb 2021 18:00:56 +0100
Subject: [PATCH] [projections] implement preliminary plots for the adamont v2
 data.

---
 .../adamont_data/adamont_gcm_rcm_couples.py   | 102 +++----
 .../adamont_data/adamont_scenario.py          |  21 +-
 .../adamont_data/adamont_studies.py           |   9 +-
 .../scm_models_data/abstract_study.py         |   1 -
 .../preliminary_analysis.py                   | 265 ++++++++++--------
 .../comparison_historical_visualizer.py       |   6 +-
 .../comparison_with_scm/comparison_plot.py    |   3 +-
 .../test_adamont_study.py                     |  21 +-
 8 files changed, 241 insertions(+), 187 deletions(-)

diff --git a/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py b/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py
index 2f709304..4da9dab4 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py
+++ b/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py
@@ -1,54 +1,3 @@
-gcm_rcm_couple_to_color = {
-    ('CNRM-CM5', 'CCLM4-8-17'): 'darkred',
-    ('CNRM-CM5', 'RCA4'): 'red',
-    ('CNRM-CM5', 'ALADIN53'): 'lightcoral',
-    # Adamont v2
-    ('CNRM-CM5', 'ALADIN63'): 'orangered',
-    ('CNRM-CM5', 'RACMO22E'): 'firebrick',
-
-    ('MPI-ESM-LR', 'CCLM4-8-17'): 'darkblue',
-    ('MPI-ESM-LR', 'RCA4'): 'blue',
-    ('MPI-ESM-LR', 'REMO2009'): 'lightblue',
-
-    ('HadGEM2-ES', 'CCLM4-8-17'): 'darkgreen',
-    ('HadGEM2-ES', 'RCA4'): 'green',
-    ('HadGEM2-ES', 'RACMO22E'): 'lightgreen',
-    # Adamont v2
-    ('HadGEM2-ES', 'RegCM4-6'): 'chartreuse',
-
-    ('EC-EARTH', 'CCLM4-8-17'): 'darkviolet',
-    ('EC-EARTH', 'RCA4'): 'violet',
-    # adamont v2
-    ('EC-EARTH', 'RACMO22E'): 'mediumorchid',
-
-    ('IPSL-CM5A-MR', 'WRF331F'): 'darkorange',
-    ('IPSL-CM5A-MR', 'RCA4'): 'orange',
-    # adamont v2
-    ('IPSL-CM5A-MR', 'WRF381P'): 'moccasin',
-
-    ('NorESM1-M', 'HIRHAM5'): 'yellow',
-    # adamont v2
-    ('NorESM1-M', 'REMO2015'): 'gold',
-
-    # adamont v2
-    ('ERAINT', 'ALADIN62'): 'deeppink'
-}
-
-gcm_to_color = {
-    'CNRM-CM5': 'red',
-
-    'MPI-ESM-LR': 'blue',
-
-    'HadGEM2-ES': 'green',
-
-    'EC-EARTH': 'violet',
-
-    'IPSL-CM5A-MR': 'orange',
-
-    'NorESM1-M': 'yellow',
-}
-
-
 
 def get_year_min_and_year_max_used_to_compute_quantile(gcm):
     if gcm == 'HadGEM2-ES':
@@ -136,5 +85,56 @@ gcm_rcm_couple_adamont_v2_to_full_name = {
     # ('ERAINT', 'ALADIN62'): 'CNRM-ALADIN62_ECMWF-ERAINT',
 }
 
+gcm_rcm_couple_to_color = {
+    ('CNRM-CM5', 'CCLM4-8-17'): 'darkred',
+    ('CNRM-CM5', 'RCA4'): 'red',
+    ('CNRM-CM5', 'ALADIN53'): 'lightcoral',
+    # Adamont v2
+    ('CNRM-CM5', 'ALADIN63'): 'orangered',
+    ('CNRM-CM5', 'RACMO22E'): 'firebrick',
+
+    ('MPI-ESM-LR', 'CCLM4-8-17'): 'darkblue',
+    ('MPI-ESM-LR', 'RCA4'): 'blue',
+    ('MPI-ESM-LR', 'REMO2009'): 'lightblue',
+
+    ('HadGEM2-ES', 'CCLM4-8-17'): 'darkgreen',
+    ('HadGEM2-ES', 'RCA4'): 'green',
+    ('HadGEM2-ES', 'RACMO22E'): 'lightgreen',
+    # Adamont v2
+    ('HadGEM2-ES', 'RegCM4-6'): 'chartreuse',
+
+    ('EC-EARTH', 'CCLM4-8-17'): 'darkviolet',
+    ('EC-EARTH', 'RCA4'): 'violet',
+    # adamont v2
+    ('EC-EARTH', 'RACMO22E'): 'mediumorchid',
+
+    ('IPSL-CM5A-MR', 'WRF331F'): 'darkorange',
+    ('IPSL-CM5A-MR', 'RCA4'): 'orange',
+    # adamont v2
+    ('IPSL-CM5A-MR', 'WRF381P'): 'moccasin',
+
+    ('NorESM1-M', 'HIRHAM5'): 'yellow',
+    # adamont v2
+    ('NorESM1-M', 'REMO2015'): 'gold',
+
+    # adamont v2
+    # ('ERAINT', 'ALADIN62'): 'deeppink'
+}
+
+gcm_to_color = {
+    'CNRM-CM5': 'red',
+
+    'MPI-ESM-LR': 'blue',
+
+    'HadGEM2-ES': 'green',
+
+    'EC-EARTH': 'violet',
+
+    'IPSL-CM5A-MR': 'orange',
+
+    'NorESM1-M': 'yellow',
+}
+
+
 if __name__ == '__main__':
     print(get_gcm_list(adamont_version=2))
diff --git a/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py b/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py
index 155e82eb..c53dc397 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py
+++ b/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py
@@ -1,7 +1,6 @@
 from enum import Enum
 
-from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import gcm_rcm_couple_to_color, \
-    get_gcm_rcm_couple_adamont_to_full_name
+from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import get_gcm_rcm_couple_adamont_to_full_name
 
 
 class AdamontScenario(Enum):
@@ -13,6 +12,7 @@ class AdamontScenario(Enum):
 
 
 adamont_scenarios_real = [AdamontScenario.histo, AdamontScenario.rcp26, AdamontScenario.rcp45, AdamontScenario.rcp85]
+rcp_scenarios = [AdamontScenario.rcp26, AdamontScenario.rcp45, AdamontScenario.rcp85]
 
 
 def get_linestyle_from_scenario(adamont_scenario):
@@ -65,6 +65,20 @@ def get_year_min(adamont_scenario, gcm_rcm_couple):
     return year_min
 
 
+def get_gcm_rcm_couple_adamont_version_2(scenario):
+    s = set(list(get_gcm_rcm_couple_adamont_to_full_name(version=2).keys()))
+    scenario_to_list_to_remove = {
+        AdamontScenario.rcp26: [('EC-EARTH', 'CCLM4-8-17'), ('CNRM-CM5', 'ALADIN53'), ('CNRM-CM5', 'RCA4'),
+                                ('MPI-ESM-LR', 'RCA4'), ('HadGEM2-ES', 'CCLM4-8-17'), ('IPSL-CM5A-MR', 'RCA4'),
+                                ('CNRM-CM5', 'CCLM4-8-17'), ('IPSL-CM5A-MR', 'WRF381P'), ('NorESM1-M', 'HIRHAM5'),
+                                ('IPSL-CM5A-MR', 'WRF331F'), ('HadGEM2-ES', 'RCA4')],
+        AdamontScenario.rcp45:  [('NorESM1-M', 'REMO2015'), ('HadGEM2-ES', 'RegCM4-6')],
+        AdamontScenario.rcp85: [],
+    }
+    for couple_to_remove in scenario_to_list_to_remove[scenario]:
+        s.remove(couple_to_remove)
+    return list(s)
+
 def load_gcm_rcm_couples(year_min=None, year_max=None,
                          adamont_scenario=AdamontScenario.histo,
                          adamont_version=2):
@@ -104,6 +118,3 @@ def scenario_to_real_scenarios(adamont_scenario):
 def gcm_rcm_couple_to_str(gcm_rcm_couple):
     return ' / '.join(gcm_rcm_couple)
 
-
-def get_color_from_gcm_rcm_couple(gcm_rcm_couple):
-    return gcm_rcm_couple_to_color[gcm_rcm_couple]
diff --git a/extreme_data/meteo_france_data/adamont_data/adamont_studies.py b/extreme_data/meteo_france_data/adamont_data/adamont_studies.py
index fab2b67b..6f494a3e 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont_studies.py
+++ b/extreme_data/meteo_france_data/adamont_data/adamont_studies.py
@@ -6,8 +6,9 @@ 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_gcm_rcm_couples import get_gcm_rcm_couple_adamont_to_full_name
-from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_str, get_color_from_gcm_rcm_couple
+from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import get_gcm_rcm_couple_adamont_to_full_name, \
+    gcm_rcm_couple_to_color
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_str
 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, ADAMONT_STUDY_CLASS_TO_ABBREVIATION
@@ -59,8 +60,8 @@ class AdamontStudies(object):
         for gcm_rcm_couple, study in list(self.gcm_rcm_couple_to_study.items())[::-1]:
             if massif_name in study.massif_name_to_annual_maxima:
                 y = study.massif_name_to_annual_maxima[massif_name]
-                label = gcm_rcm_couple_to_str(gcm_rcm_couple)
-                color = get_color_from_gcm_rcm_couple(gcm_rcm_couple)
+                label = gcm_rcm_couple_to_color[gcm_rcm_couple]
+                color = gcm_rcm_couple_to_color[gcm_rcm_couple]
                 ax.plot(x, y, linewidth=linewidth, label=label, color=color)
         if scm_study is None:
             y = np.array([study.massif_name_to_annual_maxima[massif_name] for study in self.study_list
diff --git a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
index 8ce54c7b..3cd4d53b 100644
--- a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
+++ b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
@@ -891,7 +891,6 @@ class AbstractStudy(object):
     def _massif_name_to_stationary_gev_params_and_confidence(self, quantile_level=None,
                                                              confidence_interval_based_on_delta_method=True):
         """ at least 90% of values must be above zero"""
-        print('study computation')
         massif_name_to_stationary_gev_params = {}
         massif_name_to_confidence = {}
         for massif_name, annual_maxima in self.massif_name_to_annual_maxima.items():
diff --git a/projects/altitude_spatial_model/preliminary_analysis.py b/projects/altitude_spatial_model/preliminary_analysis.py
index 6c7f5781..a31dcba5 100644
--- a/projects/altitude_spatial_model/preliminary_analysis.py
+++ b/projects/altitude_spatial_model/preliminary_analysis.py
@@ -1,9 +1,15 @@
 import matplotlib.pyplot as plt
+import os.path as op
 from itertools import chain
 
 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.adamont_snowfall import AdamontSnowfall
+from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import gcm_rcm_couple_adamont_v2_to_full_name
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import adamont_scenarios_real, \
+    get_gcm_rcm_couple_adamont_version_2, rcp_scenarios
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
 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
@@ -25,7 +31,7 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
         self.altitudes_for_temporal_hypothesis = [600, 1500, 2400, 3300]
 
     def plot_gev_params_against_altitude(self):
-        legend = True
+        legend = False
         elevation_as_xaxis = False
         param_names = GevParams.PARAM_NAMES + [100]
         if legend:
@@ -110,6 +116,8 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
                 tight_layout = False
 
             plot_name = '{} change with altitude'.format(param_name)
+            if isinstance(self.study, AbstractAdamontStudy):
+                plot_name = op.join(plot_name, gcm_rcm_couple_adamont_v2_to_full_name[self.study.gcm_rcm_couple])
 
             # # Display the legend
             if legend:
@@ -142,6 +150,10 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
             idx = 8 if param_name == GevParams.SHAPE else 1
             the = ' the' if param_name in GevParams.PARAM_NAMES else ''
             label = 'Elevation gradient for\n{} {}'.format(the, value_label[:-idx] + '/100m)')
+            plot_name = label.replace('/', ' every ')
+            if isinstance(self.study, AbstractAdamontStudy):
+                plot_name = op.join(plot_name, gcm_rcm_couple_adamont_v2_to_full_name[self.study.gcm_rcm_couple])
+
             gev_param_name_to_graduation = {
                 GevParams.LOC: 0.5,
                 GevParams.SCALE: 0.1,
@@ -153,7 +165,7 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
             visualizer.plot_map(cmap=plt.cm.coolwarm,
                                 graduation=gev_param_name_to_graduation[param_name],
                                 label=label, massif_name_to_value=massif_name_to_linear_coef,
-                                plot_name=label.replace('/', ' every '), add_x_label=False,
+                                plot_name=plot_name, add_x_label=False,
                                 negative_and_positive_values=param_name == GevParams.SHAPE,
                                 add_colorbar=True,
                                 massif_name_to_text=massif_name_to_r2_score,
@@ -185,130 +197,143 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
 
     # Plot against the time
 
-    @property
-    def year_min_and_max_list(self):
-        l = []
-        year_min, year_max = 1959, 1989
-        for shift in range(0, 7):
-            l.append((year_min + 5 * shift, year_max + 5 * shift))
-        return l
-
-    @property
-    def min_years_for_plot_x_axis(self):
-        return [c[0] for c in self.year_min_and_max_list]
-
-    def plot_gev_params_against_time_for_all_altitudes(self):
-        for altitude in self.altitudes_for_temporal_hypothesis:
-            self._plot_gev_params_against_time_for_one_altitude(altitude)
-
-    def _plot_gev_params_against_time_for_one_altitude(self, altitude):
-        for param_name in GevParams.PARAM_NAMES[:]:
-            ax = plt.gca()
-            for massif_name in self.study.all_massif_names()[:]:
-                self._plot_gev_params_against_time_for_one_altitude_and_one_massif(ax, massif_name, param_name,
-                                                                                   altitude,
-                                                                                   massif_name_as_labels=True)
-            ax.legend(prop={'size': 7}, ncol=3)
-            ax.set_xlabel('Year')
-            ax.set_ylabel(param_name + ' for altitude={}'.format(altitude))
-            xlabels = ['-'.join([str(e) for e in t]) for t in self.year_min_and_max_list]
-            ax.set_xticks(self.min_years_for_plot_x_axis)
-            ax.set_xticklabels(xlabels)
-            # ax.tick_params(labelsize=5)
-            plot_name = '{} change /all with years /for altitude={}'.format(param_name, altitude)
-            self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False)
-            ax.clear()
-            plt.close()
-
-    def _plot_gev_params_against_time_for_one_altitude_and_one_massif(self, ax, massif_name, param_name, altitude,
-                                                                      massif_name_as_labels):
-        study = self.altitude_to_study[altitude]
-        if massif_name in study.study_massif_names:
-            gev_params_list = study.massif_name_to_gev_param_list(self.year_min_and_max_list)[massif_name]
-            params = [gev_params.to_dict()[param_name] for gev_params in gev_params_list]
-            # params = np.array(params)
-            # param_normalized = params / np.sqrt(np.sum(np.power(params, 2)))
-            # confidence_intervals = [gev_params.param_name_to_confidence_interval[param_name] for gev_params in
-            #                         gev_params_list]
-            massif_id = self.study.all_massif_names().index(massif_name)
-            plot_against_altitude(self.min_years_for_plot_x_axis, ax, massif_id, massif_name, params,
-                                  altitude, False,
-                                  massif_name_as_labels)
-            # plot_against_altitude(self.years, ax, massif_id, massif_name, confidence_intervals, True)
-
-    # plot for each massif against the time
-
-    def plot_gev_params_against_time_for_all_massifs(self):
-        for massif_name in self.study.all_massif_names():
-            self._plot_gev_params_against_time_for_one_massif(massif_name)
-
-    def _plot_gev_params_against_time_for_one_massif(self, massif_name):
-        for param_name in GevParams.PARAM_NAMES[:]:
-            ax = plt.gca()
-            for altitude in self.altitudes_for_temporal_hypothesis:
-                self._plot_gev_params_against_time_for_one_altitude_and_one_massif(ax, massif_name, param_name,
-                                                                                   altitude,
-                                                                                   massif_name_as_labels=False)
-            ax.legend()
-            ax.set_xlabel('Year')
-            ax.set_ylabel(param_name + ' for {}'.format(massif_name))
-            xlabels = ['-'.join([str(e) for e in t]) for t in self.year_min_and_max_list]
-            ax.set_xticks(self.min_years_for_plot_x_axis)
-            ax.set_xticklabels(xlabels)
-            plot_name = '{} change /with years /for {}'.format(param_name, massif_name)
-            self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False)
-            ax.clear()
-            plt.close()
-
-    # PLot for each massif the derivative against the time for each altitude
-
-    def plot_time_derivative_against_time(self):
-        for param_name in GevParams.PARAM_NAMES[:]:
-            ax = plt.gca()
-            for massif_name in self.study.all_massif_names()[:]:
-                self._plot_gev_params_time_derivative_against_altitude_one_massif(ax, massif_name, param_name)
-            ax.legend(prop={'size': 7}, ncol=3)
-            ax.set_xlabel('Altitude')
-            ax.set_ylabel(param_name)
-            plot_name = '{} change /time derivative with altitude'.format(param_name)
-            self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False)
-            ax.clear()
-            plt.close()
-
-    def _plot_gev_params_time_derivative_against_altitude_one_massif(self, ax, massif_name, param_name):
-        altitudes = []
-        time_derivatives = []
-        for altitude, study in self.altitude_to_study.items():
-            if (massif_name in study.study_massif_names) and ("Mercan" not in massif_name):
-                gev_params_list = study.massif_name_to_gev_param_list(self.year_min_and_max_list)[massif_name]
-                params = [gev_params.to_dict()[param_name] for gev_params in gev_params_list]
-                x = list(range(len(params)))
-                y = params
-                a = self.get_robust_slope(x, y)
-                time_derivatives.append(a)
-                altitudes.append(altitude)
-        massif_id = self.study.all_massif_names().index(massif_name)
-        plot_against_altitude(altitudes, ax, massif_id, massif_name, time_derivatives, fill=False)
-
-    def get_robust_slope(self, x, y):
-        a, *_ = fit_linear_regression(x=x, y=y)
-        a_list = [a]
-        for i in range(len(x)):
-            x_copy, y_copy = x[:], y[:]
-            x_copy.pop(i)
-            y_copy.pop(i)
-            a, *_ = fit_linear_regression(x=x_copy, y=y_copy)
-            a_list.append(a)
-        return np.mean(np.array(a_list))
-
-
-if __name__ == '__main__':
+    # @property
+    # def year_min_and_max_list(self):
+    #     l = []
+    #     year_min, year_max = 1959, 1989
+    #     for shift in range(0, 7):
+    #         l.append((year_min + 5 * shift, year_max + 5 * shift))
+    #     return l
+
+    # @property
+    # def min_years_for_plot_x_axis(self):
+    #     return [c[0] for c in self.year_min_and_max_list]
+    #
+    # def plot_gev_params_against_time_for_all_altitudes(self):
+    #     for altitude in self.altitudes_for_temporal_hypothesis:
+    #         self._plot_gev_params_against_time_for_one_altitude(altitude)
+    #
+    # def _plot_gev_params_against_time_for_one_altitude(self, altitude):
+    #     for param_name in GevParams.PARAM_NAMES[:]:
+    #         ax = plt.gca()
+    #         for massif_name in self.study.all_massif_names()[:]:
+    #             self._plot_gev_params_against_time_for_one_altitude_and_one_massif(ax, massif_name, param_name,
+    #                                                                                altitude,
+    #                                                                                massif_name_as_labels=True)
+    #         ax.legend(prop={'size': 7}, ncol=3)
+    #         ax.set_xlabel('Year')
+    #         ax.set_ylabel(param_name + ' for altitude={}'.format(altitude))
+    #         xlabels = ['-'.join([str(e) for e in t]) for t in self.year_min_and_max_list]
+    #         ax.set_xticks(self.min_years_for_plot_x_axis)
+    #         ax.set_xticklabels(xlabels)
+    #         # ax.tick_params(labelsize=5)
+    #         plot_name = '{} change /all with years /for altitude={}'.format(param_name, altitude)
+    #         self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False)
+    #         ax.clear()
+    #         plt.close()
+    #
+    # def _plot_gev_params_against_time_for_one_altitude_and_one_massif(self, ax, massif_name, param_name, altitude,
+    #                                                                   massif_name_as_labels):
+    #     study = self.altitude_to_study[altitude]
+    #     if massif_name in study.study_massif_names:
+    #         gev_params_list = study.massif_name_to_gev_param_list(self.year_min_and_max_list)[massif_name]
+    #         params = [gev_params.to_dict()[param_name] for gev_params in gev_params_list]
+    #         # params = np.array(params)
+    #         # param_normalized = params / np.sqrt(np.sum(np.power(params, 2)))
+    #         # confidence_intervals = [gev_params.param_name_to_confidence_interval[param_name] for gev_params in
+    #         #                         gev_params_list]
+    #         massif_id = self.study.all_massif_names().index(massif_name)
+    #         plot_against_altitude(self.min_years_for_plot_x_axis, ax, massif_id, massif_name, params,
+    #                               altitude, False,
+    #                               massif_name_as_labels)
+    #         # plot_against_altitude(self.years, ax, massif_id, massif_name, confidence_intervals, True)
+    #
+    # # plot for each massif against the time
+    #
+    # def plot_gev_params_against_time_for_all_massifs(self):
+    #     for massif_name in self.study.all_massif_names():
+    #         self._plot_gev_params_against_time_for_one_massif(massif_name)
+    #
+    # def _plot_gev_params_against_time_for_one_massif(self, massif_name):
+    #     for param_name in GevParams.PARAM_NAMES[:]:
+    #         ax = plt.gca()
+    #         for altitude in self.altitudes_for_temporal_hypothesis:
+    #             self._plot_gev_params_against_time_for_one_altitude_and_one_massif(ax, massif_name, param_name,
+    #                                                                                altitude,
+    #                                                                                massif_name_as_labels=False)
+    #         ax.legend()
+    #         ax.set_xlabel('Year')
+    #         ax.set_ylabel(param_name + ' for {}'.format(massif_name))
+    #         xlabels = ['-'.join([str(e) for e in t]) for t in self.year_min_and_max_list]
+    #         ax.set_xticks(self.min_years_for_plot_x_axis)
+    #         ax.set_xticklabels(xlabels)
+    #         plot_name = '{} change /with years /for {}'.format(param_name, massif_name)
+    #         self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False)
+    #         ax.clear()
+    #         plt.close()
+    #
+    # # PLot for each massif the derivative against the time for each altitude
+    #
+    # def plot_time_derivative_against_time(self):
+    #     for param_name in GevParams.PARAM_NAMES[:]:
+    #         ax = plt.gca()
+    #         for massif_name in self.study.all_massif_names()[:]:
+    #             self._plot_gev_params_time_derivative_against_altitude_one_massif(ax, massif_name, param_name)
+    #         ax.legend(prop={'size': 7}, ncol=3)
+    #         ax.set_xlabel('Altitude')
+    #         ax.set_ylabel(param_name)
+    #         plot_name = '{} change /time derivative with altitude'.format(param_name)
+    #         self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False)
+    #         ax.clear()
+    #         plt.close()
+    #
+    # def _plot_gev_params_time_derivative_against_altitude_one_massif(self, ax, massif_name, param_name):
+    #     altitudes = []
+    #     time_derivatives = []
+    #     for altitude, study in self.altitude_to_study.items():
+    #         if (massif_name in study.study_massif_names) and ("Mercan" not in massif_name):
+    #             gev_params_list = study.massif_name_to_gev_param_list(self.year_min_and_max_list)[massif_name]
+    #             params = [gev_params.to_dict()[param_name] for gev_params in gev_params_list]
+    #             x = list(range(len(params)))
+    #             y = params
+    #             a = self.get_robust_slope(x, y)
+    #             time_derivatives.append(a)
+    #             altitudes.append(altitude)
+    #     massif_id = self.study.all_massif_names().index(massif_name)
+    #     plot_against_altitude(altitudes, ax, massif_id, massif_name, time_derivatives, fill=False)
+    #
+    # def get_robust_slope(self, x, y):
+    #     a, *_ = fit_linear_regression(x=x, y=y)
+    #     a_list = [a]
+    #     for i in range(len(x)):
+    #         x_copy, y_copy = x[:], y[:]
+    #         x_copy.pop(i)
+    #         y_copy.pop(i)
+    #         a, *_ = fit_linear_regression(x=x_copy, y=y_copy)
+    #         a_list.append(a)
+    #     return np.mean(np.array(a_list))
+
+def main_paper2():
     altitudes = list(chain.from_iterable(altitudes_for_groups))
 
     # altitudes = paper_altitudes
     altitudes = [1800, 2100]
     visualizer = PointwiseGevStudyVisualizer(SafranSnowfall1Day, altitudes=altitudes)
     visualizer.plot_gev_params_against_altitude()
+
     # visualizer.plot_gev_params_against_time_for_all_altitudes()
     # visualizer.plot_gev_params_against_time_for_all_massifs()
     # visualizer.plot_time_derivative_against_time()
+
+
+def main_paper3():
+    altitudes = list(chain.from_iterable(altitudes_for_groups))
+    altitudes = [1200, 1500, 1800]
+    for scenario in rcp_scenarios[:2]:
+        for gcm_rcm_couple in get_gcm_rcm_couple_adamont_version_2(scenario):
+            visualizer = PointwiseGevStudyVisualizer(AdamontSnowfall, altitudes=altitudes, scenario=scenario,
+                                                     gcm_rcm_couple=gcm_rcm_couple)
+            visualizer.plot_gev_params_against_altitude()
+
+if __name__ == '__main__':
+    main_paper3()
diff --git a/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py b/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py
index a4fe1766..e2c26b2e 100644
--- a/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py
+++ b/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py
@@ -5,8 +5,8 @@ from matplotlib.lines import Line2D
 import os.path as op
 
 from extreme_data.meteo_france_data.adamont_data.abstract_adamont_study import AbstractAdamontStudy
-from extreme_data.meteo_france_data.adamont_data.adamont_scenario import get_color_from_gcm_rcm_couple, \
-    gcm_rcm_couple_to_str, gcm_rcm_couple_to_color, scenario_to_str
+from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import gcm_rcm_couple_to_color
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_str, scenario_to_str
 from extreme_data.meteo_france_data.adamont_data.adamont_studies import AdamontStudies
 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 \
@@ -139,7 +139,7 @@ class ComparisonHistoricalVisualizer(StudyVisualizer):
         width = 10
         positions = [i * width * 2 for i in range(len(values))]
         labels = ['Reanalysis'] + [gcm_rcm_couple_to_str(couple) for couple in gcm_rcm_couples]
-        colors = ['black'] + [get_color_from_gcm_rcm_couple(couple) for couple in gcm_rcm_couples]
+        colors = ['black'] + [gcm_rcm_couple_to_color[couple] for couple in gcm_rcm_couples]
         # Permute values, labels & colors, based on the mean values
         mean_values = np.mean(values, axis=1)
         index_to_sort = np.argsort(mean_values)
diff --git a/projects/projected_snowfall/comparison_with_scm/comparison_plot.py b/projects/projected_snowfall/comparison_with_scm/comparison_plot.py
index 93735b48..398b9f3c 100644
--- a/projects/projected_snowfall/comparison_with_scm/comparison_plot.py
+++ b/projects/projected_snowfall/comparison_with_scm/comparison_plot.py
@@ -6,7 +6,8 @@ from typing import Dict
 import numpy as np
 from matplotlib.lines import Line2D
 
-from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_color, gcm_rcm_couple_to_str, \
+from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import gcm_rcm_couple_to_color
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_str, \
     scenario_to_str
 from projects.projected_snowfall.comparison_with_scm.comparison_historical_visualizer import \
     ComparisonHistoricalVisualizer
diff --git a/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py b/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py
index b287d74f..618c0bac 100644
--- a/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py
+++ b/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py
@@ -1,7 +1,8 @@
 import unittest
 
 from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import get_gcm_rcm_couple_adamont_to_full_name
-from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, adamont_scenarios_real, \
+    rcp_scenarios, get_gcm_rcm_couple_adamont_version_2
 from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
 
@@ -30,7 +31,7 @@ class TestAdamontStudy(unittest.TestCase):
 
     def test_massifs_names_adamont_v2(self):
         year_min = 2004
-        adamont_version = 2 # this test will not pass with adamont version 1
+        adamont_version = 2  # this test will not pass with adamont version 1
         for altitude in [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600]:
             reanalysis_study = SafranSnowfall1Day(altitude=altitude, year_min=year_min)
             for gcm_rcm_couple in get_gcm_rcm_couple_adamont_to_full_name(adamont_version).keys():
@@ -38,6 +39,22 @@ class TestAdamontStudy(unittest.TestCase):
                                                 year_min=year_min, gcm_rcm_couple=gcm_rcm_couple)
                 assert set(adamont_study.study_massif_names) == set(reanalysis_study.study_massif_names)
 
+    def test_existing_gcm_rcm_couple_and_rcp(self):
+        altitude = 1800
+        for scenario in rcp_scenarios[:2]:
+            l = []
+            for gcm_rcm_couple in get_gcm_rcm_couple_adamont_version_2(scenario):
+                adamont_study = AdamontSnowfall(altitude=altitude, adamont_version=2,
+                                                year_min=2098, gcm_rcm_couple=gcm_rcm_couple,
+                                                scenario=scenario)
+                try:
+                    _ = adamont_study.year_to_annual_maxima[2098]
+                except FileNotFoundError:
+                    l.append(gcm_rcm_couple)
+            print(scenario, l)
+
+
+        self.assertTrue(True)
 
 if __name__ == '__main__':
     unittest.main()
-- 
GitLab