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 0978312011050d01022d6ad9647eee93b53f9209..a292f9ccad1ae6887238b84d52aedb88b003a153 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
@@ -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
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 a7f1831a295436e342d36a60be65e0e4171bf36f..f95decdbe6676de9c288c7ede35e7051503301b3 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py
+++ b/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)
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 498ed6c8a0df7435636a0dbd3ac6d462b40d1178..e8a7d84e74fc058ad56f71e48b68cb1b4af15f42 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont_studies.py
+++ b/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)
diff --git a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_historical_period.py b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_historical_period.py
index 4231bdcf9edaaa4b270b4762a18bc5805fc9d6f0..f6ef4e35196f6a8baabc7b214a1bc59f01cd271e 100644
--- a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_historical_period.py
+++ b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_historical_period.py
@@ -1,9 +1,11 @@
+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):
diff --git a/projects/projected_snowfall/projected_data/main_projection.py b/projects/projected_snowfall/projected_data/main_projection.py
index b11886d727c4fced2576e32fc80d560cdcdbfc82..b540a1c43ebb4fe95f3197a6634b7e77a697eb3c 100644
--- a/projects/projected_snowfall/projected_data/main_projection.py
+++ b/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)
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 d72ae3c95a754890b0168b1bd0a699d92d93d16e..902639cc0e2d6fb41b58b66d050cceabdb6b21dc 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
@@ -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)