[projections] add adamont_v2

extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py

 import os
 import os.path as op
+import subprocess
 from collections import OrderedDict
 from datetime import datetime, timedelta
 from enum import Enum
@@ -9,21 +10,28 @@ import numpy as np
 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_gcm_rcm_couples import get_gcm_rcm_couple_adamont_to_full_name
 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, \
     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"
+from extreme_data.utils import DATA_PATH
+
+ADAMONT_PATH = op.join(DATA_PATH, 'ADAMONT')
+ADAMONT_v2_PATH = op.join(DATA_PATH, 'ADAMONT_v2')
+ADAMONT_v2_WEBPATH = """https://climatedata.umr-cnrm.fr/public/dcsc/projects/ADAMONT_MONTAGNE_2020/INDICATEURS_ANNUELS/alp_flat/"""
 
 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
     YEAR_MAX = 2100
@@ -32,7 +40,9 @@ class AbstractAdamontStudy(AbstractStudy):
                  year_min=None, year_max=None,
                  multiprocessing=True, season=Season.annual,
                  french_region=FrenchRegion.alps,
-                 scenario=AdamontScenario.histo, gcm_rcm_couple=('CNRM-CM5', 'ALADIN53')):
+                 scenario=AdamontScenario.histo,
+                 gcm_rcm_couple=('CNRM-CM5', 'ALADIN53'),
+                 adamont_version=2):
         # 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
@@ -40,24 +50,81 @@ class AbstractAdamontStudy(AbstractStudy):
             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))
+                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,
                          multiprocessing=multiprocessing, season=season, french_region=french_region)
+        self.adamont_version = adamont_version
         self.gcm_rcm_couple = gcm_rcm_couple
-        self.gcm_rcm_full_name = gcm_rcm_couple_to_full_name[gcm_rcm_couple]
+        self.gcm_rcm_full_name = get_gcm_rcm_couple_adamont_to_full_name(version=self.adamont_version)[gcm_rcm_couple]
         self.scenario = scenario
         assert issubclass(self.variable_class, AbstractAdamontVariable)
         # Assert the massif_name are in the same order
         for i, massif_name in enumerate(self.all_massif_names()):
             assert massif_name == massif_number_to_massif_name[i + 1]
+        assert self.adamont_version in [1, 2]
 
     @property
     def variable_name(self):
         return scenario_to_str(self.scenario) + ' ' + super().variable_name
 
+    @cached_property
+    def year_to_annual_maxima(self) -> OrderedDict:
+        if self.adamont_version == 1:
+            return super().year_to_annual_maxima
+        else:
+            return self.load_year_to_annual_maxima_version_2()
+
+    # Loading part for adamont v2
+
+    def load_year_to_annual_maxima_version_2(self):
+        year_to_annual_maxima = OrderedDict()
+        for dataset in self.datasets:
+            annual_maxima = np.array(dataset.variables[self.variable_class.indicator_name_for_maxima])
+            assert annual_maxima.shape[1] == len(self.column_mask)
+            annual_maxima = annual_maxima[:, self.column_mask]
+            year_to_annual_maxima = OrderedDict()
+            year_min, year_max = get_year_min_and_year_max_from_scenario(self.scenario, self.gcm_rcm_couple)
+            years = list(range(year_min, year_max + 1))
+            for year, maxima in zip(years, annual_maxima):
+                if self.year_min <= year <= self.year_max:
+                    year_to_annual_maxima[year] = maxima
+        return year_to_annual_maxima
+
+    def _load_dataset(self, scenario):
+        scenario_name = scenario_to_str(scenario)
+        nc_filename = self.nc_filename_adamont_v2(scenario)
+        nc_folder = op.join(ADAMONT_v2_PATH, self.variable_folder_name, scenario_name)
+        nc_filepath = op.join(nc_folder, nc_filename)
+        # Assert that the file is present, otherwise download it
+        if not op.exists(nc_filepath):
+            self._download_year_to_annual_maxima_version_2(scenario, nc_folder)
+        # Load the file
+        dataset = Dataset(filename=nc_filepath)
+        return dataset
+
+    def _download_year_to_annual_maxima_version_2(self, scenario, path_folder):
+        scenario_name = self._scenario_to_str_adamont_v2(scenario)
+        directory = self.gcm_rcm_full_name + '_' + scenario_name
+        filename = self.nc_filename_adamont_v2(scenario)
+        full_path = op.join(ADAMONT_v2_WEBPATH, directory, filename)
+        # Download file
+        request = 'wget {} -P {}'.format(full_path, path_folder)
+        subprocess.run(request, shell=True)
+
+    def nc_filename_adamont_v2(self, scenario):
+        scenario_name = self._scenario_to_str_adamont_v2(scenario)
+        return '_'.join([self.variable_class.indicator_name_for_maxima, self.gcm_rcm_full_name, scenario_name]) + '.nc'
+
+    def _scenario_to_str_adamont_v2(self, scenario):
+        scenario_name = scenario_to_str(scenario)
+        if scenario is AdamontScenario.histo:
+            scenario_name += 'RICAL'
+        return scenario_name
+
     # Loading part
 
     @cached_property
@@ -99,7 +166,8 @@ class AbstractAdamontStudy(AbstractStudy):
         # Load efficiently the variable object
         # Multiprocessing is set to False, because this is not a time consuming part
         data_list_list = [year_to_data_list[year] for year in self.ordered_years]
-        year_to_variable_object = self.efficient_variable_loading(self.ordered_years, data_list_list, multiprocessing=False)
+        year_to_variable_object = self.efficient_variable_loading(self.ordered_years, data_list_list,
+                                                                  multiprocessing=False)
         return year_to_variable_object
 
     def load_variable_object(self, data_list):
@@ -110,22 +178,28 @@ class AbstractAdamontStudy(AbstractStudy):
     @cached_property
     def flat_mask(self):
         zs_list = [int(e) for e in np.array(self.datasets[0].variables['ZS'])]
+        zs_list[-10:] = [np.nan for _ in range(10)]
         if len(self.datasets) > 1:
             zs_list_bis = [int(e) for e in np.array(self.datasets[1].variables['ZS'])]
-            assert all([a == b for a, b in zip(zs_list, zs_list_bis)])
+            zs_list_bis[-10:] = [np.nan for _ in range(10)]
+            assert all([(a == b) or (np.isnan(a) and np.isnan(b)) for a, b in zip(zs_list, zs_list_bis)])
         return np.array(zs_list) == self.altitude
 
     @cached_property
     def study_massif_names(self) -> List[str]:
-        massif_ids = np.array(self.datasets[0].variables['MASSIF_NUMBER'])[self.flat_mask]
+        massif_key = 'MASSIF_NUMBER' if self.adamont_version == 1 else 'massif_number'
+        massif_ids = np.array(self.datasets[0].variables[massif_key])[self.flat_mask]
         if len(self.datasets) > 1:
-            massif_ids_bis = np.array(self.datasets[1].variables['MASSIF_NUMBER'])[self.flat_mask]
+            massif_ids_bis = np.array(self.datasets[1].variables[massif_key])[self.flat_mask]
             assert all(massif_ids == massif_ids_bis)
         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]
+        if self.adamont_version == 1:
+            return [Dataset(file_path) for file_path in self.nc_file_paths]
+        else:
+            return [self._load_dataset(scenario) for scenario in self.adamont_real_scenarios]
 
     # PATHS
 
@@ -152,7 +226,8 @@ class AbstractAdamontStudy(AbstractStudy):
     @property
     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):
+        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,

extreme_data/meteo_france_data/adamont_data/adamont/adamont_snowfall.py

@@ -17,11 +17,13 @@ class AdamontSnowfall(AbstractAdamontStudy):
                  year_min=None, year_max=None,
                  multiprocessing=True, season=Season.annual,
                  french_region=FrenchRegion.alps,
-                 scenario=AdamontScenario.histo, gcm_rcm_couple=('CNRM-CM5', 'ALADIN53')):
+                 scenario=AdamontScenario.histo, gcm_rcm_couple=('CNRM-CM5', 'ALADIN53'),
+                 adamont_version=2):
         super().__init__(SafranSnowfallSimulationVariable, altitude,
                          year_min, year_max,
-                         multiprocessing, season, french_region, scenario, gcm_rcm_couple)
+                         multiprocessing, season, french_region, scenario, gcm_rcm_couple, adamont_version)
+
 
 if __name__ == '__main__':
-    study = AdamontSnowfall(altitude=1800)
+    study = AdamontSnowfall(altitude=1800, adamont_version=2)
     print(study.year_to_annual_maxima)

extreme_data/meteo_france_data/adamont_data/adamont/adamont_variables.py

@@ -2,6 +2,7 @@ import numpy as np
 
 from extreme_data.meteo_france_data.scm_models_data.abstract_variable import AbstractVariable
 from extreme_data.meteo_france_data.scm_models_data.safran.safran_variable import SafranSnowfallVariable
+from root_utils import classproperty
 
 
 class AbstractAdamontVariable(AbstractVariable):
@@ -10,6 +11,9 @@ class AbstractAdamontVariable(AbstractVariable):
     def variable_name_for_folder_and_nc_file(cls):
         return cls.keyword()
 
+    @classmethod
+    def indicator_name_for_maxima(cls):
+        raise NotImplementedError
 
 class SafranSnowfallSimulationVariable(AbstractAdamontVariable):
     UNIT = SafranSnowfallVariable.UNIT
@@ -26,3 +30,8 @@ class SafranSnowfallSimulationVariable(AbstractAdamontVariable):
     @classmethod
     def variable_name_for_folder_and_nc_file(cls):
         return 'Snow'
+
+    @classproperty
+    def indicator_name_for_maxima(cls):
+        return 'max-1day-snowf'
+

extreme_data/meteo_france_data/adamont_data/adamont_studies.py

@@ -31,7 +31,7 @@ class AdamontStudies(object):
         return list(self.gcm_rcm_couple_to_study.values())
 
     @cached_property
-    def study(self) -> AbstractStudy:
+    def study(self) -> AbstractAdamontStudy:
         return self.study_list[0]
 
     @property
@@ -88,7 +88,8 @@ class AdamontStudies(object):
         ax.set_ylim((ylim_min, ylim_max * 1.5))
         ax.tick_params(axis='both', which='major', labelsize=13)
         handles, labels = ax.get_legend_handles_labels()
-        ax.legend(handles[::-1], labels[::-1], ncol=2, prop={'size': 7})
+        ncol = 2 if self.study.adamont_version == 1 else 3
+        ax.legend(handles[::-1], labels[::-1], ncol=ncol, prop={'size': 7})
         plot_name = 'Annual maxima of {} in {} at {} m'.format(ADAMONT_STUDY_CLASS_TO_ABBREVIATION[self.study_class],
                                                        massif_name.replace('_', ' '),
                                                         self.study.altitude)

extreme_data/meteo_france_data/adamont_data/get_list_gcm_rcm_couples_adamont_v2.py → extreme_data/meteo_france_data/adamont_data/adamont_v2_utils.py

@@ -75,11 +75,11 @@ s = set()
 for l in links.split('\n'):
     print(l)
     if '_' in l and 'SAFRAN' not in l:
-        gcm, rcm, _ = l.split('_')
-        s.add((gcm, rcm))
+        rcm, gcm, _ = l.split('_')
+        s.add((gcm, rcm, rcm + '_' + gcm))
 print(len(s))
-for c in s:
-    print(c)
+for t in s:
+    print("('{}', '{}'): '{}',".format(*t))
 
 """
 ('SMHI-RCA4', 'ICHEC-EC-EARTH')
@@ -106,12 +106,4 @@ for c in s:
 
 """
 
-def get_year_min_and_year_max_used_to_compute_quantile(rcm):
-    if rcm == 'MOHC-HadGEM2-ES':
-        reanalysis_years = (1988, 2011)
-        model_year = (1982, 2005)
-    else:
-        reanalysis_years = (1981, 2011)
-        model_year = (1975, 2005)
-    return reanalysis_years, model_year
 

extreme_data/meteo_france_data/scm_models_data/safran/safran_max_snowf.py

@@ -6,21 +6,21 @@ from cached_property import cached_property
 from netCDF4._netCDF4 import Dataset
 
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, Safran
+from extreme_data.utils import DATA_PATH
 
 
 class AbstractSafranSnowfallMaxFiles(SafranSnowfall1Day):
-    path = """/home/erwan/Documents/projects/spatiotemporalextremes/local/spatio_temporal_datasets"""
 
     def __init__(self, safran_year, **kwargs):
         super().__init__(**kwargs)
         self.year_max = max(safran_year, self.year_max)
-        self.nc_filepath = op.join(self.path, 'SAFRAN_montagne-CROCUS_{}'.format(safran_year),
+        self.nc_filepath = op.join(DATA_PATH, 'SAFRAN_montagne-CROCUS_{}'.format(safran_year),
                                    'max-1day-snowf_SAFRAN.nc')
         self.safran_year = safran_year
 
     @property
     def ordered_years(self):
-        return [i for i in list(range(1959, self.safran_year+1)) if self.year_min <= i <= self.year_max]
+        return [i for i in list(range(1959, self.safran_year + 1)) if self.year_min <= i <= self.year_max]
 
     @cached_property
     def year_to_annual_maxima(self) -> OrderedDict:
@@ -44,9 +44,12 @@ class SafranSnowfall2020(AbstractSafranSnowfallMaxFiles):
         super().__init__(2020, **kwargs)
 
 
-
-
 class SafranSnowfall2019(AbstractSafranSnowfallMaxFiles):
 
     def __init__(self, **kwargs):
         super().__init__(2019, **kwargs)
+
+
+if __name__ == '__main__':
+    study = SafranSnowfall2020(altitude=1800)
+    print(len(study.column_mask))

extreme_data/utils.py

+DATA_PATH = """/home/erwan/Documents/projects/spatiotemporalextremes/local/spatio_temporal_datasets"""

projects/projected_snowfall/projected_data/main_projection.py

@@ -18,11 +18,11 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua
 
 
 def main():
-    fast = None
-    adamont_scenario = AdamontScenario.rcp85_extended
+    fast = True
+    adamont_scenario = [AdamontScenario.histo, AdamontScenario.rcp85_extended][0]
     year_min = 1982 if adamont_scenario is AdamontScenario.rcp85_extended else 2006
     # Set the year_min and year_max for the comparison
-    if fast is True:
+    if fast is None:
         year_max = [2030][0]
         massif_names = ['Vanoise']
         altitudes = [1800]
@@ -31,7 +31,7 @@ def main():
         # year_min = [1951][0]
         year_max = [2005][0]
         massif_names = ['Vercors']
-        altitudes =  [900, 1200, 1500, 1800, 2100, 2400]
+        altitudes = [900, 1200, 1500, 1800, 2100, 2400][3:4]
     else:
         year_max = [2100][0]
         massif_names = None
@@ -46,7 +46,8 @@ def main():
         adamont_studies = AdamontStudies(adamont_study_class, gcm_rcm_couples,
                                          altitude=altitude, year_min=year_min,
                                          year_max=year_max, season=season,
-                                         scenario=adamont_scenario)
+                                         scenario=adamont_scenario,
+                                         adamont_version=2)
         adamont_studies.plot_maxima_time_series_adamont(massif_names)
 
 

test/test_extreme_data/test_meteo_france_data/test_adamont_study.py

@@ -7,22 +7,23 @@ from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import
 class TestAdamontStudy(unittest.TestCase):
 
     def test_load_adamont_snowfall(self):
-        study_list = [
-            AdamontSnowfall(altitude=900),
-            AdamontSnowfall(altitude=1800)
-        ]
-        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(len(study.massif_name_to_annual_maxima['Vanoise']))
-        self.assertTrue(True)
+        for version in [1, 2]:
+            study_list = [
+                AdamontSnowfall(altitude=900, adamont_version=version),
+                AdamontSnowfall(altitude=1800, adamont_version=version)
+            ]
+            study_list.extend([
+                AdamontSnowfall(altitude=900, scenario=AdamontScenario.rcp45, adamont_version=version),
+                AdamontSnowfall(altitude=900, scenario=AdamontScenario.rcp85, adamont_version=version),
+                AdamontSnowfall(altitude=900, scenario=AdamontScenario.rcp85_extended, adamont_version=version)
+            ])
+            study_list.extend([AdamontSnowfall(altitude=900, gcm_rcm_couple=gcm_rcm_couple, adamont_version=version)
+                               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(len(study.massif_name_to_annual_maxima['Vanoise']))
+            self.assertTrue(True)
 
 
 if __name__ == '__main__':