diff --git a/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py b/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py
index 13dc1f9ee96a8c6342e27105abfc768416f4d01a..7953d28f77c4adee834586378dc56838a02b3b3f 100644
--- a/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py
+++ b/extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py
@@ -59,8 +59,8 @@ def get_gcm_list(adamont_version):
gcm_to_rnumber = \
{
- 'MPI-ESM-LR': 2,
- 'CNRM-CM5': 11,
+ 'MPI-ESM-LR': 1,
+ 'CNRM-CM5': 1,
'IPSL-CM5A-MR': 1,
'EC-EARTH': 12,
'HadGEM2-ES': 1,
diff --git a/extreme_data/meteo_france_data/adamont_data/cmip5/climate_explorer_cimp5.py b/extreme_data/meteo_france_data/adamont_data/cmip5/climate_explorer_cimp5.py
new file mode 100644
index 0000000000000000000000000000000000000000..1ddb4aacd31acdec123b6d77c27ae50330f98337
--- /dev/null
+++ b/extreme_data/meteo_france_data/adamont_data/cmip5/climate_explorer_cimp5.py
@@ -0,0 +1,94 @@
+import calendar
+import os.path as op
+import pandas as pd
+import subprocess
+from datetime import datetime, timedelta
+
+import cdsapi
+import numpy as np
+from netCDF4._netCDF4 import Dataset, OrderedDict
+
+from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import gcm_to_rnumber, get_gcm_list
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import get_year_min_and_year_max_from_scenario, \
+ AdamontScenario, adamont_scenarios_real
+from extreme_data.utils import DATA_PATH
+
+GLOBALTEMP_WEB_PATH = "https://climexp.knmi.nl/CMIP5/Tglobal/"
+GLOBALTEMP_DATA_PATH = op.join(DATA_PATH, 'CMIP5_global_temp')
+
+
+def get_scenario_name(scenario):
+ if scenario is AdamontScenario.histo:
+ return 'historicalNat'
+ else:
+ return str(scenario).split('.')[-1]
+
+
+def year_to_global_mean_temp(gcm, scenario):
+ # Compute everything
+ ensemble_member = 'r{}i1p1'.format(gcm_to_rnumber[gcm])
+ scenario_name = get_scenario_name(scenario)
+
+ # Standards
+ mean_annual_column_name = 'Annual mean'
+ filename = 'global_tas_Amon_{}_{}_{}'.format(gcm, scenario_name, ensemble_member)
+ dat_filepath = op.join(GLOBALTEMP_DATA_PATH, filename + '.dat')
+ txt_filepath = op.join(GLOBALTEMP_DATA_PATH, filename + '.txt')
+ csv_filepath = op.join(GLOBALTEMP_DATA_PATH, filename + '.csv')
+ # Download if needed
+ if not op.exists(txt_filepath):
+ download_dat(dat_filepath, txt_filepath)
+ # Transform nc file into csv file
+ if not op.exists(csv_filepath):
+ dat_to_csv(csv_filepath, txt_filepath, mean_annual_column_name)
+
+ # Load csv file
+ df = pd.read_csv(csv_filepath, index_col=0)
+ d = OrderedDict(df[mean_annual_column_name])
+ print(gcm, scenario_name, np.mean(list(d.values())))
+ return d
+
+
+def dat_to_csv(csv_filepath, txt_filepath, mean_annual_column_name):
+ d = OrderedDict()
+ with open(txt_filepath, 'r') as f:
+ for i, l in enumerate(f):
+ year, l = l[:8], l[8:]
+ month_temp = [float(f) for f in l.split()]
+ assert len(month_temp) == 12
+ d[year] = list(month_temp)
+ df = pd.DataFrame.from_dict(d)
+ df = df.transpose()
+ df.columns = list(calendar.month_abbr)[1:]
+ df[mean_annual_column_name] = df.mean(axis=1)
+ df.to_csv(csv_filepath)
+
+
+def download_dat(dat_filepath, txt_filepath):
+ web_filepath = op.join(GLOBALTEMP_WEB_PATH, op.basename(dat_filepath))
+ dirname = op.dirname(dat_filepath)
+ requests = [
+ 'wget {} -P {}'.format(web_filepath, dirname),
+ 'tail -n +4 {} > {}'.format(dat_filepath, txt_filepath),
+ ]
+ for request in requests:
+ subprocess.run(request, shell=True)
+
+
+def main_example():
+ scenario = AdamontScenario.rcp45
+ gcm = 'EC-EARTH'
+ year_to_global_mean_temp(gcm, scenario)
+
+
+def main_test_cmip5_loader():
+ for scenario in adamont_scenarios_real[1:]:
+ for gcm in get_gcm_list(adamont_version=2)[:]:
+ print(gcm, scenario)
+ year_to_global_temp = year_to_global_mean_temp(gcm, scenario)
+ print(year_to_global_temp)
+
+
+if __name__ == '__main__':
+ # main_example()
+ main_test_cmip5_loader()
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 2dd13ecbe41fa077173baf882bc09f04804b6e33..8ce54c7b924be5ddb410c6f5b796d323d17505fe 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
@@ -927,9 +927,8 @@ class AbstractStudy(object):
massif_name_to_stationary_gev_params[massif_name] = gev_params_list
return massif_name_to_stationary_gev_params
- @property
- def mean_annual_maxima(self):
+ def aggregate_annual_maxima(self, f):
annual_maxima = []
for maxima in self.year_to_annual_maxima.values():
annual_maxima.extend(maxima)
- return np.mean(annual_maxima)
+ return f(annual_maxima)
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 8317e6e9f54391a8b8db6c5db309f8dca3f6076d..ad64be482fd83e50b11a59804be5d1a0f30a4c20 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
@@ -27,6 +27,7 @@ def compute_bias_and_display_it(ax,
gcm_rcm_couple,
massif_names=None,
relative_bias=False,
+ mean=True,
):
bias_in_the_mean_maxima = []
altitudes = []
@@ -37,16 +38,17 @@ def compute_bias_and_display_it(ax,
altitudes.append(altitude)
+ f = np.mean if mean else np.std
if massif_names is None:
- mean_maxima_adamont = adamont_study.mean_annual_maxima
- mean_maxima_reanalysis = study_reanalysis.mean_annual_maxima
+ mean_maxima_adamont = adamont_study.aggregate_annual_maxima(f)
+ mean_maxima_reanalysis = study_reanalysis.aggregate_annual_maxima(f)
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]))
+ mean_maxima_reanalysis = f(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 = f(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
if relative_bias:
@@ -63,7 +65,7 @@ def compute_bias_and_display_it(ax,
def main_comparaison_plot():
altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
- for adamont_version in [1, 2][:]:
+ 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)
@@ -87,67 +89,69 @@ def main_comparaison_plot():
else:
list_of_massis_names = [None] + [[m] for m in reanalysis_altitude_studies_1981.study.all_massif_names()]
- for relative_bias in [True, False][:1]:
- 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,
- relative_bias)
- 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])
- unit = '%' if relative_bias else study.variable_unit
- bias_name = 'Relative bias' if relative_bias else 'Bias'
- title = '{} in the mean annual maxima of {} of {}\n' \
- 'for ADAMONT v{}' \
- ' against {} on the quantile mapping period ({})'.format(bias_name,
- study_str, massif_str,
- adamont_version,
- comparaison_study_class,
- unit)
- folder = 'relative bias' if relative_bias else 'absolute bias'
- plot_name = op.join(folder, title)
- ax.set_xlabel(title, fontsize=10)
- reanalysis_altitude_studies.show_or_save_to_file(plot_name=plot_name, no_title=True)
- plt.close()
+ for relative_bias in [True, False][:]:
+ for mean in [True, False][1:]:
+ 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,
+ relative_bias, mean)
+ 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])
+ unit = '%' if relative_bias else study.variable_unit
+ bias_name = 'Relative bias' if relative_bias else 'Bias'
+ mean_str = 'mean' if mean else 'std'
+ title = '{} in the {} annual maxima of {} of {}\n' \
+ 'for ADAMONT v{}' \
+ ' against {} on the quantile mapping period ({})'.format(bias_name, mean_str,
+ study_str, massif_str,
+ adamont_version,
+ comparaison_study_class,
+ unit)
+ folder = 'relative bias' if relative_bias else 'absolute bias'
+ plot_name = op.join(folder, 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__':