From 146e4da1cbe9e056e29b51a7f10c73922836a11a Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Thu, 22 Oct 2020 16:29:12 +0200
Subject: [PATCH] add comparison with scm. add massif_name_to_daily_time_series
to abstract_study.py
---
.../adamont_data/abstract_adamont_study.py | 3 +-
.../adamont_data/adamont_scenario.py | 46 ++++++++--
.../adamont_data/adamont_studies.py | 27 ++++++
.../scm_models_data/abstract_study.py | 16 +++-
.../scm_models_data/safran/safran.py | 5 +-
projects/projected_snowfall/__init__.py | 0
.../comparison_with_scm/__init__.py | 0
.../comparison_historical_visualizer.py | 84 +++++++++++++++++++
.../main_comparison_on_historical_period.py | 50 +++++++++++
9 files changed, 216 insertions(+), 15 deletions(-)
create mode 100644 extreme_data/meteo_france_data/adamont_data/adamont_studies.py
create mode 100644 projects/projected_snowfall/__init__.py
create mode 100644 projects/projected_snowfall/comparison_with_scm/__init__.py
create mode 100644 projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py
create mode 100644 projects/projected_snowfall/comparison_with_scm/main_comparison_on_historical_period.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 b03510e8..5657e7e7 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
@@ -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])
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 3e20e6b2..08701866 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py
+++ b/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',
-
-
}
diff --git a/extreme_data/meteo_france_data/adamont_data/adamont_studies.py b/extreme_data/meteo_france_data/adamont_data/adamont_studies.py
new file mode 100644
index 00000000..5a3ce2a1
--- /dev/null
+++ b/extreme_data/meteo_france_data/adamont_data/adamont_studies.py
@@ -0,0 +1,27 @@
+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
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 24432d06..e5ac7977 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
@@ -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
diff --git a/extreme_data/meteo_france_data/scm_models_data/safran/safran.py b/extreme_data/meteo_france_data/scm_models_data/safran/safran.py
index fab3a2d2..df4fdc01 100644
--- a/extreme_data/meteo_france_data/scm_models_data/safran/safran.py
+++ b/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)
diff --git a/projects/projected_snowfall/__init__.py b/projects/projected_snowfall/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/projects/projected_snowfall/comparison_with_scm/__init__.py b/projects/projected_snowfall/comparison_with_scm/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py b/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py
new file mode 100644
index 00000000..70cb8933
--- /dev/null
+++ b/projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py
@@ -0,0 +1,84 @@
+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()
+
+
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
new file mode 100644
index 00000000..79084a9f
--- /dev/null
+++ b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_historical_period.py
@@ -0,0 +1,50 @@
+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
--
GitLab