From 5ca08a78605e3e679b6ee37034faa08e03fab8ae Mon Sep 17 00:00:00 2001
From: "Daniel.Frisinghelli" <daniel.frisinghelli@eurac.edu>
Date: Fri, 1 Oct 2021 14:52:13 +0000
Subject: [PATCH] Refactor.
---
Notebooks/eval_precipitation.ipynb | 21 +-
.../eval_precipitation_sensitivity.ipynb | 111 +++++-
Notebooks/eval_temperature.ipynb | 344 ++++++++++++++++--
3 files changed, 412 insertions(+), 64 deletions(-)
diff --git a/Notebooks/eval_precipitation.ipynb b/Notebooks/eval_precipitation.ipynb
index 91b8a83..a0c652d 100644
--- a/Notebooks/eval_precipitation.ipynb
+++ b/Notebooks/eval_precipitation.ipynb
@@ -55,10 +55,12 @@
"PREDICTAND = 'pr'\n",
"MODEL = 'USegNet'\n",
"PPREDICTORS = 'ztuvq'\n",
+ "# PPREDICTORS = ''\n",
"PLEVELS = ['500', '850']\n",
+ "# PLEVELS = []\n",
"SPREDICTORS = 'p'\n",
- "DEM = 'dem'\n",
- "DEM_FEATURES = ''\n",
+ "DEM = ''\n",
+ "DEM_FEATURES = 'dem'\n",
"DOY = ''\n",
"WET_DAY_THRESHOLD = '2'\n",
"# LOSS = 'MSELoss'\n",
@@ -265,12 +267,13 @@
"outputs": [],
"source": [
"# align datasets and mask missing values in model predictions\n",
- "y_true, y_refe, y_pred = xr.align(y_true, y_refe, y_pred, join='override')\n",
- "y_pred_pr = y_pred.precipitation.where(~np.isnan(y_true.precipitation), other=np.nan)\n",
"if LOSS == 'BernoulliGammaLoss':\n",
- " y_pred_prob = y_pred.prob.where(~np.isnan(y_true.precipitation), other=np.nan)\n",
- "y_refe = y_refe.precipitation.where(~np.isnan(y_true.precipitation), other=np.nan)\n",
- "y_true = y_true.precipitation"
+ " y_true, y_refe, y_pred_pr, y_pred_prob = xr.align(y_true.precipitation, y_refe.precipitation, y_pred.precipitation, y_pred.prob, join='override')\n",
+ " y_pred_prob = y_pred_prob.where(~np.isnan(y_true), other=np.nan)\n",
+ "else:\n",
+ " y_true, y_refe, y_pred_pr = xr.align(y_true.precipitation, y_refe.precipitation, y_pred.precipitation, join='override')\n",
+ "y_pred_pr = y_pred_pr.where(~np.isnan(y_true), other=np.nan) \n",
+ "y_refe = y_refe.where(~np.isnan(y_true), other=np.nan)"
]
},
{
@@ -1046,7 +1049,7 @@
"outputs": [],
"source": [
"# minimum precipitation (mm / day) defining a wet day\n",
- "WET_DAY_THRESHOLD = 1"
+ "WET_DAY_THRESHOLD = float(WET_DAY_THRESHOLD)"
]
},
{
@@ -1274,7 +1277,7 @@
"outputs": [],
"source": [
"# true and predicted probability of precipitation\n",
- "p_true = (y_true > WET_DAY_THRESHOLD).values.flatten()\n",
+ "p_true = (y_true > float(WET_DAY_THRESHOLD)).values.flatten()\n",
"p_pred = y_pred_prob.values.flatten()"
]
},
diff --git a/Notebooks/eval_precipitation_sensitivity.ipynb b/Notebooks/eval_precipitation_sensitivity.ipynb
index eb55dd8..8bbfdb4 100644
--- a/Notebooks/eval_precipitation_sensitivity.ipynb
+++ b/Notebooks/eval_precipitation_sensitivity.ipynb
@@ -46,7 +46,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 1,
"id": "d75e3627-13c4-46b8-b749-d582e4b802ac",
"metadata": {},
"outputs": [],
@@ -56,7 +56,7 @@
"MODEL = 'USegNet'\n",
"PPREDICTORS = 'ztuvq'\n",
"PLEVELS = ['500', '850']\n",
- "SPREDICTORS = 'mslp'\n",
+ "SPREDICTORS = 'p'\n",
"DEM = 'dem'\n",
"DEM_FEATURES = ''\n",
"DOY = ''\n",
@@ -74,7 +74,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 2,
"id": "01ed9fbc-49c2-4028-a3f3-692a014dcf89",
"metadata": {},
"outputs": [],
@@ -101,7 +101,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 3,
"id": "073dfc6a-8b33-4aca-ab71-c06ef92bbb80",
"metadata": {},
"outputs": [],
@@ -120,7 +120,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 4,
"id": "030c221d-e0ad-4e67-b742-db34de549983",
"metadata": {},
"outputs": [],
@@ -135,10 +135,26 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 5,
"id": "7106fed2-9f63-453b-a9af-b5d490efbbb4",
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "['USegNet_pr_ztuvq_500_850_p_dem_0mm_BernoulliGammaLoss',\n",
+ " 'USegNet_pr_ztuvq_500_850_p_dem_05mm_BernoulliGammaLoss',\n",
+ " 'USegNet_pr_ztuvq_500_850_p_dem_1mm_BernoulliGammaLoss',\n",
+ " 'USegNet_pr_ztuvq_500_850_p_dem_2mm_BernoulliGammaLoss',\n",
+ " 'USegNet_pr_ztuvq_500_850_p_dem_3mm_BernoulliGammaLoss',\n",
+ " 'USegNet_pr_ztuvq_500_850_p_dem_5mm_BernoulliGammaLoss']"
+ ]
+ },
+ "execution_count": 5,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"# file pattern for different wet day thresholds\n",
"PATTERNS = ['_'.join([PATTERN, '{}mm_{}'.format(str(t).replace('.', ''), LOSS)]) for t in WET_DAY_THRESHOLD]\n",
@@ -147,7 +163,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 6,
"id": "0069c5d8-ee28-49a3-84f0-8af6b3bcb05d",
"metadata": {},
"outputs": [],
@@ -158,7 +174,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 7,
"id": "316a35e3-49e6-4d02-b865-25710a621daa",
"metadata": {},
"outputs": [],
@@ -170,7 +186,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 8,
"id": "58177982-7f86-4c9e-9980-563cf5d48852",
"metadata": {},
"outputs": [],
@@ -200,7 +216,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 9,
"id": "9758d58b-82e9-47b2-9cf7-1c93d5ed04a7",
"metadata": {},
"outputs": [],
@@ -211,10 +227,23 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 10,
"id": "57ca80c6-2c02-44d1-87eb-35c8a78ec71f",
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "(0.0 mm), r2: -1.11\n",
+ "(0.5 mm), r2: -1.43\n",
+ "(1.0 mm), r2: -0.72\n",
+ "(2.0 mm), r2: -1.22\n",
+ "(3.0 mm), r2: 0.30\n",
+ "(5.0 mm), r2: -0.32\n"
+ ]
+ }
+ ],
"source": [
"# calculate predicted monthly mean precipitation (mm / month)\n",
"for k, v in y_pred.items():\n",
@@ -248,7 +277,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 11,
"id": "02c99a93-530c-4570-a7ed-4cbf1940222f",
"metadata": {},
"outputs": [],
@@ -259,10 +288,35 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 12,
"id": "bc85224d-68d8-4791-8eec-ab39e16898b7",
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "(0.0 mm), relative bias = -50.56%\n",
+ "(0.0 mm), MAE = 1.50 mm\n",
+ "(0.0 mm), RMSE = 2.88 mm\n",
+ "(0.5 mm), relative bias = -53.85%\n",
+ "(0.5 mm), MAE = 1.60 mm\n",
+ "(0.5 mm), RMSE = 3.08 mm\n",
+ "(1.0 mm), relative bias = -44.11%\n",
+ "(1.0 mm), MAE = 1.32 mm\n",
+ "(1.0 mm), RMSE = 2.24 mm\n",
+ "(2.0 mm), relative bias = -48.03%\n",
+ "(2.0 mm), MAE = 1.43 mm\n",
+ "(2.0 mm), RMSE = 2.52 mm\n",
+ "(3.0 mm), relative bias = -14.26%\n",
+ "(3.0 mm), MAE = 0.62 mm\n",
+ "(3.0 mm), RMSE = 0.70 mm\n",
+ "(5.0 mm), relative bias = -30.53%\n",
+ "(5.0 mm), MAE = 0.97 mm\n",
+ "(5.0 mm), RMSE = 1.39 mm\n"
+ ]
+ }
+ ],
"source": [
"# yearly average bias over reference period\n",
"for k, v in y_pred.items():\n",
@@ -344,10 +398,23 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 13,
"id": "92733418-8840-46c2-bf33-adf93b78c646",
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "(0.0 mm), ROCSS = 0.52\n",
+ "(0.5 mm), ROCSS = 0.58\n",
+ "(1.0 mm), ROCSS = 0.60\n",
+ "(2.0 mm), ROCSS = 0.52\n",
+ "(3.0 mm), ROCSS = 0.62\n",
+ "(5.0 mm), ROCSS = 0.63\n"
+ ]
+ }
+ ],
"source": [
"# true and predicted probability of precipitation\n",
"for k, v in y_pred.items():\n",
@@ -365,6 +432,14 @@
" rocss = 2 * area - 1 # ROC skill score (cf. https://journals.ametsoc.org/view/journals/clim/16/24/1520-0442_2003_016_4145_otrsop_2.0.co_2.xml)\n",
" print('({:.1f} mm), ROCSS = {:.2f}'.format(k, rocss))"
]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d3fac03c-4215-465d-a0c5-358de4dd57f2",
+ "metadata": {},
+ "outputs": [],
+ "source": []
}
],
"metadata": {
diff --git a/Notebooks/eval_temperature.ipynb b/Notebooks/eval_temperature.ipynb
index 1cf8145..1193bc3 100644
--- a/Notebooks/eval_temperature.ipynb
+++ b/Notebooks/eval_temperature.ipynb
@@ -52,7 +52,7 @@
"outputs": [],
"source": [
"# define the model parameters\n",
- "PREDICTAND = 'tasmin'\n",
+ "PREDICTAND = 'tasmax'\n",
"MODEL = 'USegNet'\n",
"PPREDICTORS = 'ztuvq'\n",
"PLEVELS = ['500', '850']\n",
@@ -97,6 +97,7 @@
"import xarray as xr\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
+ "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
"import scipy.stats as stats\n",
"from IPython.display import Image\n",
"from sklearn.metrics import r2_score\n",
@@ -179,7 +180,7 @@
"outputs": [],
"source": [
"# model predictions and observations NetCDF\n",
- "y_pred = xr.open_dataset(search_files(TARGET_PATH.joinpath(PREDICTAND), '(.*)'.join([PATTERN, '.nc$'])).pop())\n",
+ "y_pred = xr.open_dataset(search_files(TARGET_PATH.joinpath(PREDICTAND), '.'.join([PATTERN, 'nc$'])).pop())\n",
"if PREDICTAND == 'tas':\n",
" # read both tasmax and tasmin\n",
" tasmax = xr.open_dataset(search_files(OBS_PATH.joinpath('tasmax'), '.nc$').pop())\n",
@@ -200,16 +201,59 @@
"y_true = y_true.sel(time=y_pred.time)"
]
},
+ {
+ "cell_type": "markdown",
+ "id": "06dcb26f-f8e0-4365-9279-871bc35492f9",
+ "metadata": {},
+ "source": [
+ "### Load ERA-5 reference dataset"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
- "id": "966d85fb-9185-408f-ac2b-1e4ca829ccd1",
+ "id": "8f20564f-3f1b-4da6-b625-f2d0c6392249",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# search ERA-5 reference dataset\n",
+ "y_refe = xr.open_dataset(search_files(ERA5_PATH.joinpath('ERA5', '2m_{}_temperature'.format(PREDICTAND.lstrip('tas'))), '.nc$').pop())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f7b2fe84-131d-485a-a6e7-3c47135447e1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# convert to °C\n",
+ "y_refe = y_refe - 273.15"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6be4c268-df97-4bd0-948a-cc1943a3c0fb",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# subset to time period covered by predictions\n",
+ "y_refe = y_refe.sel(time=y_pred.time).drop_vars('lambert_azimuthal_equal_area')\n",
+ "y_refe = y_refe.rename({'t2m': PREDICTAND})"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9c7efab4-0cd2-4f6d-afbd-5e284d2bda9d",
"metadata": {},
"outputs": [],
"source": [
"# align datasets and mask missing values in model predictions\n",
- "y_true, y_pred = xr.align(y_true, y_pred, join='override')\n",
- "y_pred = y_pred.where(~np.isnan(y_true), other=np.nan)"
+ "y_true, y_pred, y_refe = xr.align(y_true, y_pred, y_refe, join='override')\n",
+ "y_pred = y_pred.where(~np.isnan(y_true), other=np.nan)\n",
+ "y_refe = y_refe.where(~np.isnan(y_true), other=np.nan)"
]
},
{
@@ -228,6 +272,18 @@
"### Coefficient of determination"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bb8223e0-c5cb-477c-8c31-3dfb94b20ce1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# group timeseries by month and calculate mean over time and space\n",
+ "y_pred_ac = y_pred.groupby('time.month').mean(dim=('time', 'y', 'x'))\n",
+ "y_true_ac = y_true.groupby('time.month').mean(dim=('time', 'y', 'x'))"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
@@ -236,8 +292,8 @@
"outputs": [],
"source": [
"# get predicted and observed values over entire time series and grid points\n",
- "y_pred_values = y_pred[PREDICTAND].values.flatten()\n",
- "y_true_values = y_true[PREDICTAND].values.flatten()"
+ "y_pred_values = y_pred[PREDICTAND].groupby('time.month').mean(dim='time').values.flatten()\n",
+ "y_true_values = y_true[PREDICTAND].groupby('time.month').mean(dim='time').values.flatten()"
]
},
{
@@ -263,7 +319,8 @@
"outputs": [],
"source": [
"# calculate coefficient of determination\n",
- "r2 = r2_score(y_true_values, y_pred_values)"
+ "r2 = r2_score(y_true_values, y_pred_values)\n",
+ "r2"
]
},
{
@@ -277,29 +334,45 @@
"fig, ax = plt.subplots(1, 1, figsize=(10, 10))\n",
"\n",
"# plot only a subset of data: otherwise plot is overloaded ...\n",
- "subset = np.random.choice(np.arange(0, len(y_pred_values)), size=int(1e7), replace=False)\n",
- "ax.plot(y_true_values[subset], y_pred_values[subset], 'o', alpha=.5, markeredgecolor='grey', markerfacecolor='none', markersize=3);\n",
+ "# subset = np.random.choice(np.arange(0, len(y_pred_values)), size=int(1e7), replace=False)\n",
+ "# ax.plot(y_true_values[subset], y_pred_values[subset], 'o', alpha=.5, markeredgecolor='grey', markerfacecolor='none', markersize=3);\n",
"\n",
"# plot entire dataset\n",
- "# ax.plot(y_true_values, y_pred_values, 'o', alpha=.5, markeredgecolor='grey', markerfacecolor='none', markersize=3);\n",
+ "ax.plot(y_true_values, y_pred_values, 'o', alpha=.5, markeredgecolor='grey', markerfacecolor='none', markersize=3);\n",
"\n",
"# plot 1:1 mapping line\n",
- "interval = np.arange(-40, 45, 5)\n",
+ "interval = np.arange(-15, 45, 5)\n",
+ "#interval = np.arange(-30, 35, 5)\n",
"ax.plot(interval, interval, color='k', lw=2, ls='--')\n",
"\n",
"# add coefficient of determination: calculated on entire dataset!\n",
- "ax.text(interval[-1] - 1, interval[0] + 1, s='Coefficient of determination R$^2$ = {:.2f}'.format(r2), ha='right', fontsize=14)\n",
+ "ax.text(interval[-1] - 0.5, interval[0] + 0.5, s='Coefficient of determination R$^2$ = {:.2f}'.format(r2), ha='right', fontsize=18)\n",
"\n",
"# format axes\n",
"ax.set_ylim(interval[0], interval[-1])\n",
"ax.set_xlim(interval[0], interval[-1])\n",
"ax.set_xticks(interval)\n",
- "ax.set_xticklabels(interval, fontsize=14)\n",
+ "ax.set_xticklabels(interval, fontsize=16)\n",
"ax.set_yticks(interval)\n",
- "ax.set_yticklabels(interval, fontsize=14)\n",
- "ax.set_xlabel('Observed', fontsize=14)\n",
- "ax.set_ylabel('Predicted', fontsize=14)\n",
- "ax.set_title('{} (°C): 1991 - 2010'.format(NAMES[PREDICTAND].capitalize()), fontsize=16, pad=10);\n",
+ "ax.set_yticklabels(interval, fontsize=16)\n",
+ "ax.set_xlabel('Observed', fontsize=18)\n",
+ "ax.set_ylabel('Predicted', fontsize=18)\n",
+ "ax.set_title('Monthly mean {} (°C)'.format(NAMES[PREDICTAND]), fontsize=20, pad=10);\n",
+ "\n",
+ "# add axis for annual cycle\n",
+ "axins = inset_axes(ax, width=\"30%\", height=\"40%\", loc=2, borderpad=1)\n",
+ "axins.plot(y_pred_ac[PREDICTAND].values, ls='--', color='k', label='Predicted')\n",
+ "axins.plot(y_true_ac[PREDICTAND].values, ls='-', color='k', label='Observed')\n",
+ "axins.legend(frameon=False, fontsize=12, loc='lower center');\n",
+ "axins.yaxis.tick_right()\n",
+ "#axins.set_yticks(np.arange(-10, 11, 2))\n",
+ "#axins.set_yticklabels(np.arange(-10, 11, 2), fontsize=12)\n",
+ "axins.set_yticks(np.arange(-0, 22, 2))\n",
+ "axins.set_yticklabels(np.arange(0, 22, 2), fontsize=12)\n",
+ "axins.set_xticks(np.arange(0, 12))\n",
+ "axins.set_xticklabels([calendar.month_name[i + 1] for i in np.arange(0, 12)], rotation=90, fontsize=12)\n",
+ "# axins.set_ylabel('{} / (°C)'.format(NAMES[PREDICTAND].capitalize()), fontsize=14)\n",
+ "# axins.set_xlabel('Month', fontsize=14);\n",
"\n",
"# save figure\n",
"fig.savefig('../Notebooks/Figures/{}_r2.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')"
@@ -330,36 +403,43 @@
"source": [
"# yearly average bias over reference period\n",
"y_pred_yearly_avg = y_pred.groupby('time.year').mean(dim='time')\n",
+ "y_refe_yearly_avg = y_refe.groupby('time.year').mean(dim='time')\n",
"y_true_yearly_avg = y_true.groupby('time.year').mean(dim='time')\n",
"bias_yearly_avg = y_pred_yearly_avg - y_true_yearly_avg\n",
+ "bias_yearly_avg_ref = y_refe_yearly_avg - y_true_yearly_avg\n",
"for var in bias_yearly_avg:\n",
- " print('Yearly average bias of {}: {:.2f}'.format(var, bias_yearly_avg[var].mean().item()))"
+ " print('(Model) Yearly average bias of {}: {:.2f}°C'.format(var, bias_yearly_avg[var].mean().item()))\n",
+ " print('(ERA-5) Yearly average bias of {}: {:.2f}°C'.format(var, bias_yearly_avg_ref.mean().to_array().values.item()))"
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "2ef19dce-29ac-4f6f-9999-e67745a4afd1",
+ "id": "fa14fadb-0f76-4962-80cb-e69cd7f5fac5",
"metadata": {},
"outputs": [],
"source": [
"# mean absolute error over reference period\n",
"mae_avg = np.abs(y_pred_yearly_avg - y_true_yearly_avg).mean()\n",
+ "mae_avg_ref = np.abs(y_refe_yearly_avg - y_true_yearly_avg).mean()\n",
"for var in mae_avg:\n",
- " print('Yearly average MAE of {}: {:.2f}'.format(var, mae_avg[var].item()))"
+ " print('(Model) Yearly average MAE of {}: {:.2f}°C'.format(var, mae_avg[var].item()))\n",
+ " print('(ERA-5) Yearly average MAE of {}: {:.2f}°C'.format(var, mae_avg_ref.mean().to_array().values.item()))"
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "4c4dd156-f763-482c-84e6-c4329bfd3fe4",
+ "id": "d2c50119-3cee-43c5-838b-e84639af55db",
"metadata": {},
"outputs": [],
"source": [
"# root mean squared error over reference period\n",
"rmse_avg = ((y_pred_yearly_avg - y_true_yearly_avg) ** 2).mean()\n",
+ "rmse_avg_ref = ((y_refe_yearly_avg - y_true_yearly_avg) **2).mean()\n",
"for var in rmse_avg:\n",
- " print('Yearly average RMSE of {}: {:.2f}'.format(var, rmse_avg[var].item()))"
+ " print('(Model) Yearly average RMSE of {}: {:.2f}°C'.format(var, rmse_avg[var].item()))\n",
+ " print('(ERA-5) Yearly average RMSE of {}: {:.2f}°C'.format(var, rmse_avg_ref.mean().to_array().values.item()))"
]
},
{
@@ -463,27 +543,42 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "24aadff5-b19d-4f4b-a4c5-32ee656e64cd",
+ "id": "8b781300-2639-437d-9875-f57306715c80",
"metadata": {},
"outputs": [],
"source": [
"# group data by season: (DJF, MAM, JJA, SON)\n",
"y_true_snl = y_true.groupby('time.season').mean(dim='time')\n",
"y_pred_snl = y_pred.groupby('time.season').mean(dim='time')\n",
- "bias_snl = y_pred_snl - y_true_snl"
+ "y_refe_snl = y_refe.groupby('time.season').mean(dim='time')\n",
+ "bias_snl = y_pred_snl - y_true_snl\n",
+ "bias_snl_ref = y_refe_snl - y_true_snl"
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "2d0d5a46-0652-4289-9e1c-6d47aafcfef0",
+ "id": "b148ec7e-f832-40f2-b287-c0351e7e1aad",
"metadata": {},
"outputs": [],
"source": [
- "# print average bias per season\n",
+ "# print average bias per season: ERA-5\n",
+ "for var in bias_snl_ref.data_vars:\n",
+ " for season in bias_snl_ref[PREDICTAND].season:\n",
+ " print('(ERA-5) Average bias of mean {} for season {}: {:.1f}°C'.format(var, season.values.item(), bias_snl_ref[var].sel(season=season).mean().item()))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e54daae8-ea32-431f-8a34-f92d95e7627e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# print average bias per season: model\n",
"for var in bias_snl.data_vars:\n",
" for season in bias_snl[PREDICTAND].season:\n",
- " print('Average bias of {} for season {}: {:.2f}'.format(var, season.values.item(), bias_snl[var].sel(season=season).mean().item()))"
+ " print('(Model) Average bias of mean {} for season {}: {:.1f}°C'.format(var, season.values.item(), bias_snl[var].sel(season=season).mean().item()))"
]
},
{
@@ -612,7 +707,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "c3da76d3-7261-4084-b5ec-f65682fd6596",
+ "id": "8fd82b83-0824-4663-846c-53fc3afbc7d1",
"metadata": {},
"outputs": [],
"source": [
@@ -620,46 +715,116 @@
"with warnings.catch_warnings():\n",
" warnings.simplefilter('ignore', category=RuntimeWarning)\n",
" y_pred_ex = y_pred.groupby('time.year').quantile(quantile, dim='time')\n",
- " y_true_ex = y_true.groupby('time.year').quantile(quantile, dim='time')"
+ " y_true_ex = y_true.groupby('time.year').quantile(quantile, dim='time')\n",
+ " y_refe_ex = y_refe.groupby('time.year').quantile(quantile, dim='time')"
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "20db96c8-e04f-4acb-886b-abb740863fbb",
+ "id": "822984f8-cdaa-410b-96e8-e58da39ba40d",
"metadata": {},
"outputs": [],
"source": [
"# calculate bias in extreme quantile for each year\n",
"bias_ex = y_pred_ex - y_true_ex\n",
+ "bias_ex_ref = y_refe_ex - y_true_ex"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "efd902ee-1154-4b0c-bfbd-a64ace625e76",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# bias of extreme quantile: ERA-5\n",
+ "for var in bias_ex_ref:\n",
+ " print('(ERA-5) Yearly average bias for P{:.0f} of {}: {:.1f}°C'.format(quantile * 100, var, bias_ex_ref[var].mean().item()))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9343d8a4-07a6-4c9c-85f6-f670efe7c6cd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# bias of extreme quantile: Model\n",
"for var in bias_ex:\n",
- " print('Yearly average bias for P{:.0f} of {}: {:.2f}'.format(quantile * 100, var, bias_ex[var].mean().item()))"
+ " print('(Model) Yearly average bias for P{:.0f} of {}: {:.1f}°C'.format(quantile * 100, var, bias_ex[var].mean().item()))"
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "c7aeff55-deef-4d3c-a251-eac877c9afd9",
+ "id": "ba38c293-6a6c-4ce4-b5a2-8e468e69cb44",
"metadata": {},
"outputs": [],
"source": [
"# mean absolute error in extreme quantile\n",
"mae_ex = np.abs(y_pred_ex - y_true_ex).mean()\n",
+ "mae_ex_ref = np.abs(y_refe_ex - y_true_ex).mean()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d71dd376-f445-4c27-8ffc-7ec9f1583b52",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# mae of extreme quantile: ERA-5\n",
+ "for var in mae_ex_ref:\n",
+ " print('(ERA-5) Yearly average MAE for P{:.0f} of {}: {:.1f}°C'.format(quantile * 100, var, mae_ex_ref[var].item()))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "440cbb08-5a06-4708-baa5-3ca3bdb0675a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# mae of extreme quantile: Model\n",
"for var in mae_ex:\n",
- " print('Yearly average MAE for P{:.0f} of {}: {:.2f}'.format(quantile * 100, var, mae_ex[var].item()))"
+ " print('(Model) Yearly average MAE for P{:.0f} of {}: {:.1f}°C'.format(quantile * 100, var, mae_ex[var].item()))"
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "2dc6b7cf-7ae1-4b0e-8483-376eab59f5dd",
+ "id": "fd263446-b57f-4f4a-86fd-5b5596aa32ee",
"metadata": {},
"outputs": [],
"source": [
- "# root mean squared error over reference period\n",
+ "# root mean squared error in extreme quantile\n",
"rmse_ex = ((y_pred_ex - y_true_ex) ** 2).mean()\n",
+ "rmse_ex_ref = ((y_refe_ex - y_true_ex) ** 2).mean()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "90245a82-992f-4592-8749-d7a21c1f3a21",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# rmse of extreme quantile: ERA-5\n",
+ "for var in rmse_ex_ref:\n",
+ " print('(ERA-5) Yearly average RMSE for P{:.0f} of {}: {:.1f}°C'.format(quantile * 100, var, rmse_ex_ref[var].item()))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "92591caa-8ed4-4f91-a5ea-32f6b095cdfe",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# rmse of extreme quantile: Model\n",
"for var in rmse_ex:\n",
- " print('Yearly average RMSE for P{:.0f} of {}: {:.2f}'.format(quantile * 100, var, rmse_ex[var].item()))"
+ " print('(Model) Yearly average RMSE for P{:.0f} of {}: {:.1f}°C'.format(quantile * 100, var, rmse_ex[var].item()))"
]
},
{
@@ -728,6 +893,111 @@
"# save figure\n",
"fig.savefig('../Notebooks/Figures/{}_average_bias_p{:.0f}.png'.format(PREDICTAND, quantile * 100), dpi=300, bbox_inches='tight')"
]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b1550618-1314-4d45-a6b7-2f853dce3eb1",
+ "metadata": {},
+ "source": [
+ "### Compute ERA-5 daily minimum and maximum 2m temperature"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6a1a7bc8-2edf-46f3-9fb9-9f97ed262877",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# search ERA-5 hourly 2m temperature data\n",
+ "y_refe_list = sorted(search_files(ERA5_PATH.joinpath('Downloads', '2m_temperature'), '.nc$'))\n",
+ "y_refe_list"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b6d6e195-85e9-4ac0-9e50-79fd20042665",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "y_refe = xr.open_mfdataset(y_refe_list)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7cbdf926-9f01-49a7-9e61-d864af4a3c38",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "t_max = y_refe.resample(time='D').max(dim='time')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "71b38603-5ee3-4fd8-8360-67bdcb2b8f47",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "t_min = y_refe.resample(time='D').min(dim='time')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "33c09067-9060-4432-b33a-831b4ec59c66",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "t_max.to_netcdf('./tmax_tmp.nc', engine='h5netcdf')\n",
+ "t_min.to_netcdf('./tmin_tmp.nc', engine='h5netcdf')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a3a8d0fe-3d14-4a10-8314-59a5319551ee",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "grid = '/mnt/CEPH_PROJECTS/FACT_CLIMAX/OBS_DATA/grid_1km_laea'"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bb0de3d6-b4ef-4228-a2c5-a06b0e540af9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import pathlib\n",
+ "target = pathlib.Path('/mnt/CEPH_PROJECTS/FACT_CLIMAX/REANALYSIS/ERA5/2m_{}_temperature/ERA5_2m_{}_temperature_1981_2010.nc')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9083153e-bde5-4c06-88f8-a55ba693fdda",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from climax.core.utils import reproject_cdo"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cc973c07-0a24-4b3f-a5c8-5b9fc8347cdc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for name, ds in zip(['min', 'max'], ['./tmin_tmp.nc', './tmax_tmp.nc']):\n",
+ " out = pathlib.Path(str(target).format(name, name))\n",
+ " #out.parent.mkdir(parents=True, exist_ok=True)\n",
+ " reproject_cdo(grid, ds, out, mode='bilinear', overwrite=True)"
+ ]
}
],
"metadata": {
--
GitLab