From 17267d4250ba06779bdc91c6ee9aa71cfc920f2b Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 8 Feb 2021 15:52:09 +0100
Subject: [PATCH] [projections] add adamont_v2
---
.../adamont_data/abstract_adamont_study.py | 95 +++++++++++++++++--
.../adamont_data/adamont/adamont_snowfall.py | 8 +-
.../adamont_data/adamont/adamont_variables.py | 9 ++
.../adamont_data/adamont_studies.py | 5 +-
...ples_adamont_v2.py => adamont_v2_utils.py} | 16 +---
.../safran/safran_max_snowf.py | 13 ++-
extreme_data/utils.py | 1 +
.../projected_data/main_projection.py | 11 ++-
.../test_adamont_study.py | 33 +++----
9 files changed, 138 insertions(+), 53 deletions(-)
rename extreme_data/meteo_france_data/adamont_data/{get_list_gcm_rcm_couples_adamont_v2.py => adamont_v2_utils.py} (90%)
create mode 100644 extreme_data/utils.py
diff --git a/extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py b/extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py
index 63fc61f9..f059bc91 100644
--- a/extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py
+++ b/extreme_data/meteo_france_data/adamont_data/abstract_adamont_study.py
@@ -1,5 +1,6 @@
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,
diff --git a/extreme_data/meteo_france_data/adamont_data/adamont/adamont_snowfall.py b/extreme_data/meteo_france_data/adamont_data/adamont/adamont_snowfall.py
index 719e9293..e36570cc 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont/adamont_snowfall.py
+++ b/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)
diff --git a/extreme_data/meteo_france_data/adamont_data/adamont/adamont_variables.py b/extreme_data/meteo_france_data/adamont_data/adamont/adamont_variables.py
index e58b8489..81609475 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont/adamont_variables.py
+++ b/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'
+
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 548d8a3c..5d0e54f9 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont_studies.py
+++ b/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)
diff --git a/extreme_data/meteo_france_data/adamont_data/get_list_gcm_rcm_couples_adamont_v2.py b/extreme_data/meteo_france_data/adamont_data/adamont_v2_utils.py
similarity index 90%
rename from extreme_data/meteo_france_data/adamont_data/get_list_gcm_rcm_couples_adamont_v2.py
rename to extreme_data/meteo_france_data/adamont_data/adamont_v2_utils.py
index e46fafbc..17660905 100644
--- a/extreme_data/meteo_france_data/adamont_data/get_list_gcm_rcm_couples_adamont_v2.py
+++ b/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
diff --git a/extreme_data/meteo_france_data/scm_models_data/safran/safran_max_snowf.py b/extreme_data/meteo_france_data/scm_models_data/safran/safran_max_snowf.py
index 43881f41..05fd2d9e 100644
--- a/extreme_data/meteo_france_data/scm_models_data/safran/safran_max_snowf.py
+++ b/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))
diff --git a/extreme_data/utils.py b/extreme_data/utils.py
new file mode 100644
index 00000000..f2631eb5
--- /dev/null
+++ b/extreme_data/utils.py
@@ -0,0 +1 @@
+DATA_PATH = """/home/erwan/Documents/projects/spatiotemporalextremes/local/spatio_temporal_datasets"""
diff --git a/projects/projected_snowfall/projected_data/main_projection.py b/projects/projected_snowfall/projected_data/main_projection.py
index 71d33c7e..93ebb901 100644
--- a/projects/projected_snowfall/projected_data/main_projection.py
+++ b/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)
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 902639cc..f577dd69 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
@@ -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__':
--
GitLab