[projections] add plot global temp. also centered the global temp on the...

[projections] add plot global temp. also centered the global temp on the winter year. and average on the past 30 years

extreme_data/meteo_france_data/adamont_data/adamont_gcm_rcm_couples.py

@@ -34,6 +34,21 @@ gcm_rcm_couple_to_color = {
     ('ERAINT', 'ALADIN62'): 'deeppink'
 }
 
+gcm_to_color = {
+    'CNRM-CM5': 'red',
+
+    'MPI-ESM-LR': 'blue',
+
+    'HadGEM2-ES': 'green',
+
+    'EC-EARTH': 'violet',
+
+    'IPSL-CM5A-MR': 'orange',
+
+    'NorESM1-M': 'yellow',
+}
+
+
 
 def get_year_min_and_year_max_used_to_compute_quantile(gcm):
     if gcm == 'HadGEM2-ES':

extreme_data/meteo_france_data/adamont_data/adamont_scenario.py

@@ -15,6 +15,18 @@ class AdamontScenario(Enum):
 adamont_scenarios_real = [AdamontScenario.histo, AdamontScenario.rcp26, AdamontScenario.rcp45, AdamontScenario.rcp85]
 
 
+def get_linestyle_from_scenario(adamont_scenario):
+    assert isinstance(adamont_scenario, AdamontScenario)
+    if adamont_scenario is AdamontScenario.rcp26:
+        return 'dashed'
+    elif adamont_scenario is AdamontScenario.rcp45:
+        return 'dashdot'
+    elif adamont_scenario is AdamontScenario.rcp85:
+        return 'dotted'
+    else:
+        return 'solid'
+
+
 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)

extreme_data/meteo_france_data/adamont_data/cmip5/climate_explorer_cimp5.py

@@ -24,13 +24,22 @@ def get_scenario_name(scenario):
         return str(scenario).split('.')[-1]
 
 
-def year_to_global_mean_temp(gcm, scenario):
+def year_to_global_mean_temp(gcm, scenario, year_min=None, year_max=None, rolling=None):
+    d = OrderedDict()
+    years, global_mean_temps = years_and_global_mean_temps(gcm, scenario, year_min, year_max, rolling=rolling)
+    for year, global_mean_temp in zip(years, global_mean_temps):
+        d[year] = global_mean_temp
+    return d
+
+
+def years_and_global_mean_temps(gcm, scenario, year_min=None, year_max=None, rolling=None):
     # Compute everything
     ensemble_member = 'r{}i1p1'.format(gcm_to_rnumber[gcm])
     scenario_name = get_scenario_name(scenario)
 
     # Standards
     mean_annual_column_name = 'Annual mean'
+    rolling_mean_annual_column_name = 'Rolling annual mean for window={}'.format(rolling)
     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')
@@ -40,27 +49,48 @@ def year_to_global_mean_temp(gcm, scenario):
         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)
+        dat_to_csv(csv_filepath, txt_filepath, mean_annual_column_name, rolling_mean_annual_column_name, rolling=rolling)
 
     # 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
+    if year_min is None:
+        year_min = df.index[0]
+    if year_max is None:
+        year_max = df.index[-1]
+    df = df.loc[year_min:year_max]
+    years = list(df.index)
+    assert years[0] >= year_min
+    assert years[-1] <= year_max
+    if rolling:
+        global_mean_temp = list(df[rolling_mean_annual_column_name])
+    else:
+        global_mean_temp = list(df[mean_annual_column_name])
+    return years, global_mean_temp
 
 
-def dat_to_csv(csv_filepath, txt_filepath, mean_annual_column_name):
+def dat_to_csv(csv_filepath, txt_filepath, mean_annual_column_name, rolling_mean_annual_column_name, rolling=None):
     d = OrderedDict()
     with open(txt_filepath, 'r') as f:
         for i, l in enumerate(f):
-            year, l = l[:8], l[8:]
+            year, l = int(l[:5]), l[8:]
             month_temp = [float(f) for f in l.split()]
             assert len(month_temp) == 12
-            d[year] = list(month_temp)
+            d[int(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_temp_until_july = df.iloc[1:, :7]
+    assert len(df_temp_until_july.columns) == 7
+    df_temp_after_august = df.iloc[:-1, 7:]
+    assert len(df_temp_after_august.columns) == 5
+    l = df_temp_until_july.sum(axis=1).values + df_temp_after_august.sum(axis=1).values
+    l /= 12
+    l = [np.nan] + list(l)
+    assert len(l) == len(df.index)
+    df[mean_annual_column_name] = l
+    # Computing the rolling
+    if rolling is not None:
+        df[rolling_mean_annual_column_name] = df[mean_annual_column_name].rolling(window=rolling).mean()
     df.to_csv(csv_filepath)
 
 
@@ -85,10 +115,17 @@ 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)
+            years, temps = years_and_global_mean_temps(gcm, scenario)
+            print(years)
+            print(temps)
 
+def test_rolling():
+    df = pd.DataFrame([1, 2, 3, 4, 5])
+    print(df)
+    df2 = df.rolling(window=3).mean()
+    print(df2)
 
 if __name__ == '__main__':
     # main_example()
+    # test_rolling()
     main_test_cmip5_loader()

extreme_data/meteo_france_data/adamont_data/cmip5/cmip5_global_temp.py

-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 'historical'
-    else:
-        return str(scenario).split('.')[-1]
-
-
-def get_year_min_and_year_max_for_global_temp(scenario):
-    if scenario is AdamontScenario.histo:
-        return 1951, 2005
-    else:
-        return 2006, 2100
-
-
-def get_periods(gcm, scenario):
-    if scenario is AdamontScenario.histo:
-        if gcm == 'EC-EARTH':
-            return ['195001-201212']
-        else:
-            return ['185001-200512']
-    else:
-        if gcm == 'CNRM-CM5':
-            return []
-        else:
-            return ['200601-210012']
-    
-    
-
-
-def year_to_global_mean_temp(gcm, scenario):
-    # Compute everything
-    periods = get_periods(gcm, scenario)
-    ensemble_member = 'r{}i1p1'.format(gcm_to_rnumber[gcm])
-    scenario_name = get_scenario_name(scenario)
-    year_min, year_max = get_year_min_and_year_max_for_global_temp(scenario)
-
-    # Standards
-    mean_annual_column_name = 'Annual mean'
-    zip_filepath = op.join(GLOBALTEMP_DATA_PATH, 'download.zip')
-
-
-    # Create a csv file for each period
-    for period in periods:
-        filename = 'tas_Amon_{}_{}_{}_{}'.format(gcm, scenario_name, ensemble_member, period)
-        nc_filepath = op.join(GLOBALTEMP_DATA_PATH, filename + '.nc')
-        csv_filepath = op.join(GLOBALTEMP_DATA_PATH, filename + '.csv')
-        # Download if needed
-        if not op.exists(nc_filepath):
-            download_nc(ensemble_member, gcm, period, scenario_name, zip_filepath)
-        # Transform nc file into csv file
-        if not op.exists(csv_filepath):
-            nc_to_csv(csv_filepath, mean_annual_column_name, nc_filepath, year_max, year_min)
-
-    # Concatenate all csv together into a single summary_csv_filepath
-
-
-    # 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 nc_to_csv(csv_filepath, mean_annual_column_name, nc_filepath, year_max, year_min):
-    dataset = Dataset(nc_filepath)
-    tas_list = np.array(dataset.variables['tas'])
-    tas_list = np.mean(tas_list, axis=1)
-    tas_list = np.mean(tas_list, axis=1)
-    # 'days since 1850-1-1 00:00:00'
-    time_list = np.array(dataset.variables['time'])
-    assert len(time_list) == len(tas_list)
-    start = datetime(year=1850, month=1, day=1, hour=0, minute=0, second=0)
-    date_list = [start + timedelta(days=time) for time in time_list]
-    winter_year_list = [date.year if date.month < 8 else date.year + 1 for date in date_list]
-    winter_year_to_tas_list = {winter_year: [] for winter_year in range(year_min, year_max + 1)}
-    for winter_year, tas in zip(winter_year_list, tas_list):
-        if year_min <= winter_year <= year_max:
-            winter_year_to_tas_list[winter_year].append(tas)
-    # we have monthly values
-    for tas_list in winter_year_to_tas_list.values():
-        assert len(tas_list) == 12
-    winter_year_to_mean_tas = OrderedDict()
-    for winter_year, t in winter_year_to_tas_list.items():
-        winter_year_to_mean_tas[winter_year] = np.mean(t)
-    s = pd.Series(winter_year_to_mean_tas)
-    df = pd.DataFrame({mean_annual_column_name: s})
-    df.to_csv(csv_filepath)
-
-
-def download_nc(ensemble_member, gcm, period, scenario_name, zip_filepath):
-    gcm_lower = '_'.join(gcm.lower().split('-'))
-    c = cdsapi.Client()
-    c.retrieve(
-        'projections-cmip5-monthly-single-levels',
-        {
-            'ensemble_member': ensemble_member,
-            'format': 'zip',
-            'experiment': scenario_name,
-            'variable': '2m_temperature',
-            'model': gcm_lower,
-            'period': period,
-        },
-        zip_filepath)
-    # unzip and delete
-    request_list = [
-        'unzip {} -d {}'.format(zip_filepath, op.dirname(zip_filepath)),
-        'rm {}'.format(zip_filepath)
-    ]
-    for request in request_list:
-        print(request)
-        subprocess.run(request, shell=True)
-
-
-def  main_example():
-    scenario = AdamontScenario.histo
-    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)[:]:
-            if gcm != 'CNRM-CM5':
-                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/adamont_data/cmip5/plot_temperatures.py

+import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
+
+from extreme_data.meteo_france_data.adamont_data.adamont_gcm_rcm_couples import get_gcm_list, gcm_to_color
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import get_linestyle_from_scenario, \
+    adamont_scenarios_real, AdamontScenario, scenario_to_str
+from extreme_data.meteo_france_data.adamont_data.cmip5.climate_explorer_cimp5 import year_to_global_mean_temp, \
+    years_and_global_mean_temps
+
+
+def main_plot_temperature():
+    rolling = 30
+    ax = plt.gca()
+    for gcm in get_gcm_list(adamont_version=2)[:]:
+        # Plot the historical part in solid line (this part is the same between the different scenarios)
+        linestyle = get_linestyle_from_scenario(AdamontScenario.histo)
+        plot_temperature_for_rcp_gcm(ax, gcm, AdamontScenario.rcp45, year_min=1951, year_max=2005, linestyle=linestyle,
+                                     label=gcm, rolling=rolling)
+        for scenario in adamont_scenarios_real[1:]:
+            plot_temperature_for_rcp_gcm(ax, gcm, scenario, year_min=2005, year_max=2100, rolling=rolling)
+
+    ax2 = ax.twinx()
+    legend_elements = [
+        Line2D([0], [0], color='k', lw=1, label=scenario_to_str(s),
+               linestyle=get_linestyle_from_scenario(s)) for s in adamont_scenarios_real
+    ]
+    ax2.legend(handles=legend_elements, loc='upper center')
+    ax2.set_yticks([])
+
+
+    ax.legend(loc='upper left')
+    ax.set_xlabel('Years')
+    add_str = ' averaged on the last {} years'.format(rolling) if rolling is not None else ''
+    ax.set_ylabel('Global mean Temperature{} (K)\n'
+                  'mean is taken on the year centered on the winter'.format(add_str))
+    plt.show()
+
+
+def plot_temperature_for_rcp_gcm(ax, gcm, scenario, year_min, year_max, linestyle=None, label=None, rolling=None):
+    years, global_mean_temp = years_and_global_mean_temps(gcm, scenario, year_min=year_min, year_max=year_max, rolling=rolling)
+    color = gcm_to_color[gcm]
+    if linestyle is None:
+        linestyle = get_linestyle_from_scenario(scenario)
+    ax.plot(years, global_mean_temp, linestyle=linestyle, color=color, label=label)
+
+
+if __name__ == '__main__':
+    main_plot_temperature()