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 b2924f9b3cd5fbfbac7e7f9f982809f822d17055..f9fe22353e5d294be1d123957eec3e88e4a0a199 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
@@ -13,7 +13,7 @@ from extreme_data.meteo_france_data.adamont_data.adamont.adamont_variables impor
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, get_suffix_for_the_nc_file, \
- scenario_to_real_scenarios, get_year_max
+ scenario_to_real_scenarios, get_year_max, adamont_scenarios_real
from extreme_data.meteo_france_data.adamont_data.utils.utils import massif_number_to_massif_name
from extreme_data.utils import DATA_PATH
@@ -89,6 +89,12 @@ class AbstractAdamontStudy(AbstractStudy):
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))
+ if self.scenario in adamont_scenarios_real:
+ time = dataset.variables['time']
+ msg = 'len_years={} while len_time={},' \
+ 'check year_min and year_max, ' \
+ 'check in debug mode the time field of the daatset to see the starting date'.format(years, time)
+ assert len(years) == len(time), msg
for year, maxima in zip(years, annual_maxima):
if self.year_min <= year <= self.year_max:
year_to_annual_maxima[year] = maxima
@@ -110,11 +116,10 @@ class AbstractAdamontStudy(AbstractStudy):
scenario_name = self._scenario_to_str_adamont_v2(scenario)
directory = self.gcm_rcm_full_name + '_' + scenario_name
filename = self.nc_filename_adamont_v2(scenario)
- print(directory)
- print(filename)
full_path = op.join(ADAMONT_v2_WEBPATH, directory, filename)
# Download file
request = 'wget {} -P {}'.format(full_path, path_folder)
+ print(request)
subprocess.run(request, shell=True)
def nc_filename_adamont_v2(self, scenario):
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 8f2b78714454c2b9a1d0418a1e6ce45a8f48cb8f..ea0e65fe34f2ac7c1e9cf9d1e038fb275a088723 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py
+++ b/extreme_data/meteo_france_data/adamont_data/adamont_scenario.py
@@ -43,7 +43,8 @@ def get_year_min(adamont_scenario, gcm_rcm_couple):
year_min = 1982
elif rcm == 'RCA4':
year_min = 1971
- elif gcm_rcm_couple in [('NorESM1-M', 'HIRHAM5'), ('IPSL-CM5A-MR', 'WRF331F')]:
+ elif gcm_rcm_couple in [('NorESM1-M', 'HIRHAM5'), ('IPSL-CM5A-MR', 'WRF331F'), ('CNRM-CM5', 'ALADIN63'),
+ ('IPSL-CM5A-MR', 'WRF381P')]:
year_min = 1952
else:
year_min = 1951
diff --git a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_quantile_period.py b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_quantile_period.py
index 33af02fa8ff941498591f87b1867dd7f2ae55361..a4289f85b2d2adc4ecea3954fe8f8fd18eb62979 100644
--- a/projects/projected_snowfall/comparison_with_scm/main_comparison_on_quantile_period.py
+++ b/projects/projected_snowfall/comparison_with_scm/main_comparison_on_quantile_period.py
@@ -20,29 +20,44 @@ from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import Alti
def compute_bias_and_display_it(ax,
altitude_studies_reanalysis: AltitudesStudies,
adamont_altitude_studies: AltitudesStudies,
- gcm_rcm_couple
+ gcm_rcm_couple,
+ massif_names=None,
):
bias_in_the_mean_maxima = []
altitudes = []
for altitude, study_reanalysis in altitude_studies_reanalysis.altitude_to_study.items():
- altitudes.append(altitude)
adamont_study = adamont_altitude_studies.altitude_to_study[altitude]
- mean_maxima_adamont = adamont_study.mean_annual_maxima
- mean_maxima_reanalysis = study_reanalysis.mean_annual_maxima
- bias = mean_maxima_adamont - mean_maxima_reanalysis
- bias_in_the_mean_maxima.append(bias)
+ can_compute_biais = (massif_names is None) or any([m in adamont_study.study_massif_names for m in massif_names])
+ if can_compute_biais:
+
+ altitudes.append(altitude)
+
+ if massif_names is None:
+ mean_maxima_adamont = adamont_study.mean_annual_maxima
+ mean_maxima_reanalysis = study_reanalysis.mean_annual_maxima
+ else:
+ mean_maxima_reanalysis = np.mean(np.concatenate([study_reanalysis.massif_name_to_annual_maxima[m]
+ for m in massif_names
+ if m in study_reanalysis.massif_name_to_annual_maxima]))
+ mean_maxima_adamont = np.mean(np.concatenate([adamont_study.massif_name_to_annual_maxima[m]
+ for m in massif_names
+ if m in adamont_study.massif_name_to_annual_maxima]))
+
+ bias = mean_maxima_adamont - mean_maxima_reanalysis
+ bias_in_the_mean_maxima.append(bias)
color = gcm_rcm_couple_to_color[gcm_rcm_couple]
label = gcm_rcm_couple_to_str(gcm_rcm_couple)
ax.plot(bias_in_the_mean_maxima, altitudes, label=label, color=color)
- return np.array(bias_in_the_mean_maxima)
+ return np.array(bias_in_the_mean_maxima), altitudes
def main_comparaison_plot():
altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
- for adamont_version in [1, 2]:
- ax = plt.gca()
+ for adamont_version in [1, 2][1:]:
+ print('version:', adamont_version)
+
gcm_rcm_couples = load_gcm_rcm_couples(adamont_scenario=AdamontScenario.histo, adamont_version=adamont_version)
study_class = SafranSnowfall2020 if adamont_version == 2 else SafranSnowfall1Day
@@ -57,47 +72,67 @@ def main_comparaison_plot():
altitudes=altitudes,
year_min=1988,
year_max=2011)
- bias_in_the_mean = []
- for gcm_rcm_couple in gcm_rcm_couples:
- print(gcm_rcm_couple)
- gcm, rcm = gcm_rcm_couple
- years_reanalysis, years_model = get_year_min_and_year_max_used_to_compute_quantile(gcm)
- assert years_reanalysis[0] in [1981, 1988]
- if years_reanalysis[0] == 1981:
- reanalysis_altitude_studies = reanalysis_altitude_studies_1981
- else:
- reanalysis_altitude_studies = reanalysis_altitude_studies_1988
- adamont_altitude_studies = AltitudesStudies(study_class=AdamontSnowfall,
- altitudes=altitudes,
- year_min=years_model[0],
- year_max=years_model[1],
- scenario=AdamontScenario.histo,
- gcm_rcm_couple=gcm_rcm_couple,
- adamont_version=adamont_version)
- bias_in_the_mean.append(compute_bias_and_display_it(ax, reanalysis_altitude_studies,
- adamont_altitude_studies, gcm_rcm_couple))
-
- bias_in_the_mean = np.array(bias_in_the_mean)
- min_bias, median_bias, max_bias = [f(bias_in_the_mean, axis=0) for f in [np.min, np.median, np.max]]
-
- # Plot the range for the bias, and the median
- ax.yaxis.set_ticks(altitudes)
- color = 'k'
- ax.plot(median_bias, altitudes, label='Median bias', color=color, linewidth=4)
- # ax.fill_betweenx(altitudes, min_bias, max_bias, label='Range for the bias', alpha=0.2, color='whitesmoke')
- ax.vlines(0, ymin=altitudes[0], ymax=altitudes[-1], color='k', linestyles='dashed')
- study_str = STUDY_CLASS_TO_ABBREVIATION[type(reanalysis_altitude_studies.study)]
- plot_name = 'Bias for annual maxima of {}'.format(study_str)
- ax.set_ylim(top=altitudes[-1] + 1300)
- study = adamont_altitude_studies.study
- ax.legend(ncol=3, prop={'size': 7})
- ax.set_ylabel('Altitude (m)', fontsize=10)
- ax.set_xlabel('Bias in the mean annual maxima of {} for ADAMONT v{} members\n'
- ' against {} on the quantile mapping period ({})'.format(study_str, adamont_version,
- comparaison_study_class,
- study.variable_unit), fontsize=10)
- reanalysis_altitude_studies.show_or_save_to_file(plot_name=plot_name, no_title=True)
- plt.close()
+
+ if adamont_version == 1:
+ list_of_massis_names = [None]
+ else:
+ list_of_massis_names = [None] + [[m] for m in reanalysis_altitude_studies_1981.study.all_massif_names()]
+
+ for massif_names in list_of_massis_names[:]:
+ ax = plt.gca()
+ bias_in_the_mean = []
+ list_altitudes_for_bias = []
+ for gcm_rcm_couple in gcm_rcm_couples:
+ print(massif_names, gcm_rcm_couple)
+ gcm, rcm = gcm_rcm_couple
+ years_reanalysis, years_model = get_year_min_and_year_max_used_to_compute_quantile(gcm)
+ assert years_reanalysis[0] in [1981, 1988]
+ if years_reanalysis[0] == 1981:
+ reanalysis_altitude_studies = reanalysis_altitude_studies_1981
+ else:
+ reanalysis_altitude_studies = reanalysis_altitude_studies_1988
+ adamont_altitude_studies = AltitudesStudies(study_class=AdamontSnowfall,
+ altitudes=altitudes,
+ year_min=years_model[0],
+ year_max=years_model[1],
+ scenario=AdamontScenario.histo,
+ gcm_rcm_couple=gcm_rcm_couple,
+ adamont_version=adamont_version)
+ bias, altitudes_for_bias = compute_bias_and_display_it(ax, reanalysis_altitude_studies,
+ adamont_altitude_studies, gcm_rcm_couple, massif_names)
+ bias_in_the_mean.append(bias)
+ list_altitudes_for_bias.append(altitudes_for_bias)
+
+ # Assert the all the bias have been computed for the same altitudes
+ altitudes_for_bias = list_altitudes_for_bias[0]
+ for alti in list_altitudes_for_bias:
+ assert alti == altitudes_for_bias
+
+ bias_in_the_mean = np.array(bias_in_the_mean)
+ min_bias, median_bias, max_bias = [f(bias_in_the_mean, axis=0) for f in [np.min, np.median, np.max]]
+
+ # Plot the range for the bias, and the median
+ ax.yaxis.set_ticks(altitudes)
+ color = 'k'
+ ax.plot(median_bias, altitudes_for_bias, label='Median bias', color=color, linewidth=4)
+ # ax.fill_betweenx(altitudes, min_bias, max_bias, label='Range for the bias', alpha=0.2, color='whitesmoke')
+ ax.vlines(0, ymin=altitudes[0], ymax=altitudes[-1], color='k', linestyles='dashed')
+ study_str = STUDY_CLASS_TO_ABBREVIATION[type(reanalysis_altitude_studies.study)]
+ ax.set_ylim(top=altitudes[-1] + 1300)
+ study = adamont_altitude_studies.study
+ ax.legend(ncol=3, prop={'size': 7})
+ ax.set_ylabel('Altitude (m)', fontsize=10)
+ massif_str = 'all massifs' if massif_names is None else 'the {} massif'.format(massif_names[0])
+ title = 'Bias in the mean annual maxima of {} of {}\n' \
+ 'for ADAMONT v{}' \
+ ' against {} on the quantile mapping period ({})'.format(study_str, massif_str,
+ adamont_version,
+ comparaison_study_class,
+ study.variable_unit)
+ plot_name = title
+ ax.set_xlabel(title, fontsize=10)
+ reanalysis_altitude_studies.show_or_save_to_file(plot_name=plot_name, no_title=True)
+ plt.close()
if __name__ == '__main__':
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 990f5dad7d497a9018c18a05ff905f7f9706ad20..b287d74f25ea918a29e0505156271414abc3f442 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
@@ -3,6 +3,7 @@ import unittest
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 AdamontScenario
from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
class TestAdamontStudy(unittest.TestCase):
@@ -27,6 +28,16 @@ class TestAdamontStudy(unittest.TestCase):
# print(len(study.massif_name_to_annual_maxima['Vanoise']))
self.assertTrue(True)
+ def test_massifs_names_adamont_v2(self):
+ year_min = 2004
+ adamont_version = 2 # this test will not pass with adamont version 1
+ for altitude in [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600]:
+ reanalysis_study = SafranSnowfall1Day(altitude=altitude, year_min=year_min)
+ for gcm_rcm_couple in get_gcm_rcm_couple_adamont_to_full_name(adamont_version).keys():
+ adamont_study = AdamontSnowfall(altitude=altitude, adamont_version=adamont_version,
+ year_min=year_min, gcm_rcm_couple=gcm_rcm_couple)
+ assert set(adamont_study.study_massif_names) == set(reanalysis_study.study_massif_names)
+
if __name__ == '__main__':
unittest.main()