[projections] add global temp download from climate explorer. add std...

[projections] add global temp download from climate explorer. add std computation for the evaluation of ADAMONT v2

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,

extreme_data/meteo_france_data/adamont_data/cmip5/climate_explorer_cimp5.py

+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()

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)

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__':