add comparison with scm. add massif_name_to_daily_time_series to abstract_study.py

extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py

@@ -67,7 +67,8 @@ class AbstractAdamontStudy(AbstractStudy):
         data_year_list = self.winter_year_for_each_time_step
         assert len(data_list) == len(data_year_list)
         for year_data, data in zip(data_year_list, data_list):
-            year_to_data_list[year_data].append(data)
+            if self.year_min <= year_data <= self.year_max:
+                year_to_data_list[year_data].append(data)
         year_to_variable_object = OrderedDict()
         for year in self.ordered_years:
             variable_array = np.array(year_to_data_list[year])

extreme_data/meteo_france_data/adamont_data/adamont_scenario.py

@@ -27,31 +27,59 @@ def get_year_min_and_year_max_from_scenario(adamont_scenario, gcm_rcm_couple):
 
 def get_suffix_for_the_nc_file(adamont_scenario, gcm_rcm_couple):
     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)
+    return '{}080106_{}080106_daysum'.format(year_min - 1, year_max)
 
 
 def scenario_to_str(adamont_scenario):
     return str(adamont_scenario).split('.')[-1].upper()
 
 
+def get_color_from_gcm_rcm_couple(gcm_rcm_couple):
+    return gcm_rcm_couple_to_color[gcm_rcm_couple]
+
+
+gcm_rcm_couple_to_color = {
+    ('CNRM-CM5', 'CCLM4-8-17'): 'darkred',
+    ('CNRM-CM5', 'RCA4'): 'red',
+    ('CNRM-CM5', 'ALADIN53'): 'lightcoral',
+
+    ('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',
+
+    ('EC-EARTH', 'CCLM4-8-17'): 'darkviolet',
+    ('EC-EARTH', 'RCA4'): 'violet',
+
+    ('IPSL-CM5A-MR', 'WRF331F'): 'darkorange',
+    ('IPSL-CM5A-MR', 'RCA4'): 'orange',
+
+    ('NorESM1-M', 'DMI-HIRHAM5'): 'yellow',
+
+
+}
+
 gcm_rcm_couple_to_full_name = {
     ('CNRM-CM5', 'ALADIN53'): 'CNRM-ALADIN53_CNRM-CERFACS-CNRM-CM5',
     ('CNRM-CM5', 'RCA4'): 'SMHI-RCA4_CNRM-CERFACS-CNRM-CM5',
     ('CNRM-CM5', 'CCLM4-8-17'): 'CLMcom-CCLM4-8-17_CNRM-CERFACS-CNRM-CM5',
+
     ('EC-EARTH', 'CCLM4-8-17'): 'CLMcom-CCLM4-8-17_ICHEC-EC-EARTH',
+    ('EC-EARTH', 'RCA4'): 'SMHI-RCA4_ICHEC-EC-EARTH',
+
     ('MPI-ESM-LR', 'CCLM4-8-17'): 'CLMcom-CCLM4-8-17_MPI-M-MPI-ESM-LR',
+    ('MPI-ESM-LR', 'RCA4'): 'SMHI-RCA4_MPI-M-MPI-ESM-LR',
+    ('MPI-ESM-LR', 'REMO2009'): 'MPI-CSC-REMO2009_MPI-M-MPI-ESM-LR',
+
     ('HadGEM2-ES', 'CCLM4-8-17'): 'CLMcom-CCLM4-8-17_MOHC-HadGEM2-ES',
     ('HadGEM2-ES', 'RACMO22E'): 'KNMI-RACMO22E_MOHC-HadGEM2-ES',
+    ('HadGEM2-ES', 'RCA4'): 'SMHI-RCA4_MOHC-HadGEM2-ES',
 
     ('NorESM1-M', 'DMI-HIRHAM5'): 'DMI-HIRHAM5_NCC-NorESM1-M',
+
     ('IPSL-CM5A-MR', 'WRF331F'): 'IPSL-INERIS-WRF331F_IPSL-IPSL-CM5A-MR',
-    ('EC-EARTH', 'RCA4'): 'SMHI-RCA4_ICHEC-EC-EARTH',
     ('IPSL-CM5A-MR', 'RCA4'): 'SMHI-RCA4_IPSL-IPSL-CM5A-MR',
-    ('HadGEM2-ES', 'RCA4'): 'SMHI-RCA4_MOHC-HadGEM2-ES',
-    ('MPI-ESM-LR', 'RCA4'): 'SMHI-RCA4_MPI-M-MPI-ESM-LR',
-
-
-    ('MPI-ESM-LR', 'REMO2009'): 'MPI-CSC-REMO2009_MPI-M-MPI-ESM-LR',
-
-
 }

extreme_data/meteo_france_data/adamont_data/adamont_studies.py

+from collections import OrderedDict
+
+from cached_property import cached_property
+
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_full_name
+from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+
+
+class AdamontStudies(object):
+
+    def __init__(self, study_class, gcm_rcm_couples=None, **kwargs_study):
+        self.study_class = study_class
+        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_couples_to_study = OrderedDict()  # type: OrderedDict[int, AbstractStudy]
+        for gcm_rcm_couple in self.gcm_rcm_couples:
+            study = study_class(gcm_rcm_couple=gcm_rcm_couple, **kwargs_study)
+            self.gcm_rcm_couples_to_study[gcm_rcm_couple] = study
+
+    @property
+    def study_list(self):
+        return list(self.gcm_rcm_couples_to_study.values())
+
+    @cached_property
+    def study(self) -> AbstractStudy:
+        return self.study_list[0]
\ No newline at end of file

extreme_data/meteo_france_data/scm_models_data/abstract_study.py

@@ -105,7 +105,6 @@ class AbstractStudy(object):
 
     """ Time """
 
-
     @cached_property
     def year_to_first_index_and_last_index(self):
         year_to_first_index_and_last_index = OrderedDict()
@@ -292,6 +291,15 @@ class AbstractStudy(object):
             massif_name_to_annual_maxima[massif_name] = maxima
         return massif_name_to_annual_maxima
 
+    @cached_property
+    def massif_name_to_daily_time_series(self):
+        massif_name_to_daily_time_series = OrderedDict()
+        for i, massif_name in enumerate(self.study_massif_names):
+            a = [self.year_to_daily_time_serie_array[year][:, i] for year in self.ordered_years]
+            daily_time_series = np.array(list(chain.from_iterable(a)))
+            massif_name_to_daily_time_series[massif_name] = daily_time_series
+        return massif_name_to_daily_time_series
+
     @cached_property
     def massif_name_to_annual_maxima_ordered_years(self):
         massif_name_to_annual_maxima_ordered_years = OrderedDict()
@@ -813,11 +821,11 @@ class AbstractStudy(object):
     def massif_name_to_stationary_gev_params(self):
         d, _ = self._massif_name_to_stationary_gev_params_and_confidence(quantile_level=None)
         return d
-    
+
     @cached_property
     def massif_name_to_stationary_gev_params_and_confidence_for_return_level_100(self):
         return self._massif_name_to_stationary_gev_params_and_confidence(quantile_level=0.99)
-    
+
     def _massif_name_to_stationary_gev_params_and_confidence(self, quantile_level=None):
         """ at least 90% of values must be above zero"""
         massif_name_to_stationary_gev_params = {}
@@ -844,7 +852,7 @@ class AbstractStudy(object):
         for massif_name, annual_maxima in self.massif_name_to_annual_maxima.items():
             gev_params_list = []
             for index_min, index_max in self.index_min_and_max_list:
-                annual_maxima_sliced = annual_maxima[index_min: index_max+1]
+                annual_maxima_sliced = annual_maxima[index_min: index_max + 1]
                 gev_params_list.append(fitted_stationary_gev(annual_maxima_sliced))
             massif_name_to_stationary_gev_params[massif_name] = gev_params_list
         return massif_name_to_stationary_gev_params

extreme_data/meteo_france_data/scm_models_data/safran/safran.py

@@ -209,8 +209,11 @@ class SafranTemperature(Safran):
 if __name__ == '__main__':
     altitude = 600
     year_min = 1959
-    year_max = 2019
+    year_max = 1962
     study = SafranDateFirstSnowfall(altitude=altitude, year_min=year_min, year_max=year_max)
     print(study.study_massif_names)
+    print(study.massif_name_to_annual_maxima)
+    print(study.year_to_daily_time_serie_array[1959].shape)
+    print(study.massif_name_to_daily_time_series['Vanoise'].shape)
     # print(study.year_to_annual_maxima[1959])
     # print(study.ordered_years)

projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py

+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.lines import Line2D
+
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import get_color_from_gcm_rcm_couple
+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 \
+    SCM_STUDY_CLASS_TO_ABBREVIATION
+from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import StudyVisualizer
+from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
+
+
+class ComparisonHistoricalVisualizer(StudyVisualizer):
+
+    def __init__(self, scm_study: AbstractStudy,
+                 adamont_studies: AdamontStudies,
+                 show=False,
+                 massif_names=None,
+                 ):
+        super().__init__(scm_study, show=show, save_to_file=not show)
+        self.scm_study = scm_study
+        self.adamont_studies = adamont_studies
+        if massif_names is None:
+            massif_names = scm_study.study_massif_names
+        self.massif_names = massif_names
+
+    def get_values(self, study_method, massif_name):
+        """
+        Return an array of size (nb_ensembles + 1) x nb_observations
+        :param study_method:
+        :param massif_name:
+        :return:
+        """
+        values = [self.scm_study.__getattribute__(study_method)[massif_name]]
+        for study in self.adamont_studies.study_list:
+            values.append(study.__getattribute__(study_method)[massif_name])
+        return np.array(values)
+
+    def plot_comparison(self, plot_maxima=True):
+        if plot_maxima:
+            study_method = 'massif_name_to_annual_maxima'
+        else:
+            study_method = 'massif_name_to_daily_time_series'
+        value_name = study_method.split('to_')[1]
+        print('Nb massifs', len(self.massif_names))
+        for massif_name in self.massif_names:
+            values = self.get_values(study_method, massif_name)
+            plot_name = value_name + ' for {}'.format(massif_name.replace('_', '-'))
+            self.shoe_plot_comparison(values, plot_name)
+
+    def shoe_plot_comparison(self, values, plot_name):
+        ax = plt.gca()
+        width = 10
+        positions = [i * width * 2 for i in range(len(values))]
+        labels = ['Reanalysis'] + [' / '.join(couple) for couple in self.adamont_studies.gcm_rcm_couples]
+        colors = ['black'] + [get_color_from_gcm_rcm_couple(couple) for couple in self.adamont_studies.gcm_rcm_couples]
+        # Permute values, labels & colors, based on the mean values
+        print(values.shape)
+        mean_values = np.mean(values, axis=1)
+        print(mean_values.shape)
+        index_to_sort = np.argsort(mean_values)
+        colors = [colors[i] for i in index_to_sort]
+        labels = [labels[i] for i in index_to_sort]
+        values = [values[i] for i in index_to_sort]
+        # Add boxplot with legend
+        bplot = ax.boxplot(values, positions=positions, widths=width, patch_artist=True, showmeans=True)
+        for color, patch in zip(colors, bplot['boxes']):
+            patch.set_facecolor(color)
+        custom_lines = [Line2D([0], [0], color=color, lw=4) for color in colors]
+        ax.legend(custom_lines, labels, ncol=2)
+        ax.set_xticks([])
+        ax.set_xlim([min(positions) - width, max(positions) + width])
+        ylabel = 'Annual maxima' if 'maxima' in plot_name else 'daily values'
+        ax.set_ylabel('{} of {} ({})'.format(ylabel,
+                                             SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)],
+                                                        self.study.variable_unit))
+
+        self.plot_name = '{}'.format(plot_name)
+        self.show_or_save_to_file(add_classic_title=False, no_title=True, tight_layout=True)
+        ax.clear()
+        plt.close()
+
+

projects/projected_snowfall/comparison_with_scm/main_comparison_on_historical_period.py

+from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import gcm_rcm_couple_to_full_name, \
+    get_year_min_and_year_max_from_scenario, AdamontScenario
+from extreme_data.meteo_france_data.adamont_data.adamont_studies import AdamontStudies
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
+from projects.projected_snowfall.comparison_with_scm.comparison_historical_visualizer import \
+    ComparisonHistoricalVisualizer
+
+
+def load_historical_adamont_studies(study_class, year_min, year_max):
+    gcm_rcm_couples = []
+    for gcm_rcm_couple in gcm_rcm_couple_to_full_name.keys():
+        year_min_couple, year_max_couple = get_year_min_and_year_max_from_scenario(adamont_scenario=AdamontScenario.histo,
+                                                                                   gcm_rcm_couple=gcm_rcm_couple)
+        if year_min_couple <= year_min and year_max <= year_max_couple:
+            gcm_rcm_couples.append(gcm_rcm_couple)
+    return AdamontStudies(study_class, gcm_rcm_couples, year_min=year_min, year_max=year_max)
+
+
+def main():
+    fast = [True, False][0]
+    # Set the year_min and year_max for the comparison
+    if fast:
+        year_min = [1982, 1950][1]
+        massif_names = ['Vanoise']
+        altitudes = [1800]
+    else:
+        massif_names = None
+        year_min = [1982, 1950][1]
+        altitudes = [1800, 2100]
+
+    for altitude in altitudes:
+        plot(altitude, massif_names, year_min)
+
+
+def plot(altitude, massif_names, year_min):
+    year_min = max(1959, year_min)
+    year_max = 2005
+    study_class_couple = [(SafranSnowfall1Day, AdamontSnowfall)][0]
+    scm_study_class, adamont_study_class = study_class_couple
+    scm_study = scm_study_class(altitude=altitude, year_min=year_min, year_max=year_max)
+    adamont_studies = load_historical_adamont_studies(adamont_study_class, year_min, year_max)
+
+    visualizer = ComparisonHistoricalVisualizer(scm_study, adamont_studies, massif_names=massif_names)
+    for plot_maxima in [True, False][1:]:
+        visualizer.plot_comparison(plot_maxima)
+
+
+if __name__ == '__main__':
+    main()
\ No newline at end of file