From a111a56a7c73f7cf0b708d271b29d1d598fbe9f9 Mon Sep 17 00:00:00 2001
From: "Daniel.Frisinghelli" <daniel.frisinghelli@eurac.edu>
Date: Mon, 27 Sep 2021 15:43:49 +0000
Subject: [PATCH] Compute monthly means.
---
Notebooks/eval_precipitation.ipynb | 317 +++++++++++++++++++++++++----
1 file changed, 280 insertions(+), 37 deletions(-)
diff --git a/Notebooks/eval_precipitation.ipynb b/Notebooks/eval_precipitation.ipynb
index b1ea5e2..4bdac8d 100644
--- a/Notebooks/eval_precipitation.ipynb
+++ b/Notebooks/eval_precipitation.ipynb
@@ -59,7 +59,10 @@
"SPREDICTORS = 'p'\n",
"DEM = 'dem'\n",
"DEM_FEATURES = ''\n",
- "DOY = 'doy'"
+ "DOY = 'doy'\n",
+ "# WET_DAY_THRESHOLD = 1\n",
+ "# LOSS = 'MSELoss'\n",
+ "# LOSS = 'BernoulliGammaLoss'"
]
},
{
@@ -86,12 +89,13 @@
"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, roc_curve, auc, classification_report\n",
"\n",
"# locals\n",
- "from climax.main.io import ERA5_PATH, OBS_PATH, TARGET_PATH\n",
+ "from climax.main.io import ERA5_PATH, OBS_PATH, TARGET_PATH, DEM_PATH\n",
"from pysegcnn.core.utils import search_files\n",
"from pysegcnn.core.graphics import plot_classification_report"
]
@@ -193,7 +197,21 @@
"PATTERN = '_'.join([PATTERN, SPREDICTORS]) if SPREDICTORS else PATTERN\n",
"PATTERN = '_'.join([PATTERN, DEM]) if DEM else PATTERN\n",
"PATTERN = '_'.join([PATTERN, DEM_FEATURES]) if DEM_FEATURES else PATTERN\n",
- "PATTERN = '_'.join([PATTERN, DOY]) if DOY else PATTERN"
+ "PATTERN = '_'.join([PATTERN, DOY]) if DOY else PATTERN\n",
+ "# PATTERN = '_'.join([PATTERN, '{:0d}mm_{}'.format(WET_DAY_THRESHOLD, LOSS)])\n",
+ "PATTERN"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ecaba394-f802-481f-8274-c44e2f5fdf1a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# digital elevation model\n",
+ "dem = xr.open_dataset(search_files(DEM_PATH, 'eu_dem_v11_stt.nc').pop())\n",
+ "dem = dem.Band1.to_dataset().rename({'Band1': 'elevation'})"
]
},
{
@@ -207,7 +225,18 @@
"source": [
"# model predictions and observations NetCDF \n",
"y_pred = xr.open_dataset(search_files(TARGET_PATH.joinpath(PREDICTAND), '.'.join([PATTERN, 'nc$'])).pop())\n",
- "y_true = xr.open_dataset(search_files(OBS_PATH.joinpath(PREDICTAND), '.nc$').pop())"
+ "y_true = xr.open_dataset(search_files(OBS_PATH.joinpath(PREDICTAND), 'OBS_pr(.*).nc$').pop())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "325f3086-85f0-4c28-b37d-7370d2d92405",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# reference dataset: ERA-5 precipitation\n",
+ "y_refe = xr.open_dataset(search_files(ERA5_PATH.joinpath('ERA5', 'total_precipitation'), '.nc$').pop())"
]
},
{
@@ -218,7 +247,9 @@
"outputs": [],
"source": [
"# subset to time period covered by predictions\n",
- "y_true = y_true.sel(time=y_pred.time) "
+ "y_true = y_true.sel(time=y_pred.time)\n",
+ "y_refe = y_refe.sel(time=y_pred.time).drop_vars('lambert_azimuthal_equal_area')\n",
+ "y_refe = y_refe.rename({'tp': 'precipitation'})"
]
},
{
@@ -229,9 +260,22 @@
"outputs": [],
"source": [
"# align datasets and mask missing values in model predictions\n",
- "y_true, y_pred_pr, y_pred_prob = xr.align(y_true, y_pred.precipitation.to_dataset(), y_pred.prob.to_dataset(), join='override')\n",
+ "y_true, y_refe, y_pred_pr, y_pred_prob = xr.align(y_true, y_refe, y_pred.precipitation.to_dataset(), y_pred.prob.to_dataset(), join='override')\n",
"y_pred_pr = y_pred_pr.where(~np.isnan(y_true.precipitation), other=np.nan)\n",
- "y_pred_prob = y_pred_prob.where(~np.isnan(y_true.precipitation), other=np.nan)"
+ "y_pred_prob = y_pred_prob.where(~np.isnan(y_true.precipitation), other=np.nan)\n",
+ "y_refe = y_refe.where(~np.isnan(y_true.precipitation), other=np.nan)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4e71997f-7808-463a-8b45-dcac639ebe88",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# align digital elevation model\n",
+ "_, dem = xr.align(y_true.precipitation.isel(time=0), dem, join='override')\n",
+ "dem = dem.where(~np.isnan(y_true.precipitation.isel(time=0)), other=np.nan)"
]
},
{
@@ -257,9 +301,9 @@
"metadata": {},
"outputs": [],
"source": [
- "# get predicted and observed values over entire time series and grid points\n",
- "y_pred_values = y_pred_pr[NAMES[PREDICTAND]].groupby('time.month').mean(dim='time').values.flatten()\n",
- "y_true_values = y_true[NAMES[PREDICTAND]].groupby('time.month').mean(dim='time').values.flatten()"
+ "# calculate monthly mean precipitation (mm / month)\n",
+ "y_pred_values = y_pred_pr[NAMES[PREDICTAND]].resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim='time').values\n",
+ "y_true_values = y_true[NAMES[PREDICTAND]].resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim='time').values"
]
},
{
@@ -283,7 +327,20 @@
"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"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a1b33431-9f2b-4e42-bbb0-f4fa258edd98",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# group timeseries by month and calculate mean over time and space\n",
+ "y_pred_ac = y_pred_pr.resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim=('y', 'x', 'time'), skipna=True)\n",
+ "y_true_ac = y_true.resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim=('y', 'x', 'time'), skipna=True)"
]
},
{
@@ -304,22 +361,33 @@
"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(0, 11)\n",
+ "interval = np.arange(0, 300, 50)\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], interval[0], s='Coefficient of determination R$^2$ = {:.2f}'.format(r2), ha='right', fontsize=14)\n",
+ "ax.text(interval[-1] - 2, interval[0] + 2, 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('Monthly mean {} (mm / day): 1991 - 2010'.format(NAMES[PREDICTAND]), 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 {} (mm / month)'.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=0.25)\n",
+ "axins.plot(y_pred_ac[NAMES[PREDICTAND]].values, ls='--', color='k', label='Predicted')\n",
+ "axins.plot(y_true_ac[NAMES[PREDICTAND]].values, ls='-', color='k', label='Observed')\n",
+ "axins.legend(frameon=False, fontsize=12, loc='lower center');\n",
+ "axins.set_yticks(np.arange(0, 200, 50))\n",
+ "axins.set_yticklabels(np.arange(0, 200, 50), fontsize=12)\n",
+ "axins.yaxis.tick_right()\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",
"\n",
"# save figure\n",
"fig.savefig('../Notebooks/Figures/{}_r2.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')"
@@ -350,10 +418,13 @@
"source": [
"# yearly average bias over reference period\n",
"y_pred_yearly_avg = y_pred_pr.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) / y_true_yearly_avg) * 100\n",
+ "bias_yearly_avg_ref = ((y_refe_yearly_avg - y_true_yearly_avg) / y_true_yearly_avg) * 100\n",
"for var in bias_yearly_avg:\n",
- " print('Yearly average relative bias of {}: {:.2f}%'.format(var, bias_yearly_avg[var].mean().item()))"
+ " print('(Model) Yearly average relative bias of {}: {:.2f}%'.format(var, bias_yearly_avg[var].mean().item()))\n",
+ " print('(ERA-5) Yearly average relative bias of {}: {:.2f}%'.format(var, bias_yearly_avg_ref.mean().to_array().values.item()))"
]
},
{
@@ -364,9 +435,11 @@
"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 = np.abs(y_pred_yearly_avg - y_true_yearly_avg)\n",
+ "mae_avg_ref = np.abs(y_refe_yearly_avg - y_true_yearly_avg)\n",
"for var in mae_avg:\n",
- " print('Yearly average MAE of {}: {:.2f} '.format(var, mae_avg[var].item()) + 'mm / day')"
+ " print('(Model) Yearly average MAE of {}: {:.2f} mm'.format(var, mae_avg[var].mean().item()))\n",
+ " print('(ERA-5) Yearly average MAE of {}: {:.2f} mm'.format(var, mae_avg_ref.mean().to_array().values.item()))"
]
},
{
@@ -378,8 +451,10 @@
"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()) + 'mm / day')"
+ " print('(Model) Yearly average RMSE of {}: {:.2f} mm / day'.format(var, rmse_avg[var].item()))\n",
+ " print('(ERA-5) Yearly average RMSE of {}: {:.2f} mm / day'.format(var, rmse_avg_ref.mean().to_array().values.item()))"
]
},
{
@@ -400,6 +475,74 @@
"print('Yearly average Pearson correlation coefficient for {}: {:.2f}'.format(var, np.asarray(r).mean()))"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ca683581-6b3b-4abd-ad46-157ebded19c6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# plot yearly average MAE of reference vs. prediction\n",
+ "vmin, vmax = 0, 5\n",
+ "fig, axes = plt.subplots(1, 3, figsize=(24, 8), sharex=True, sharey=True)\n",
+ "\n",
+ "# plot bias of ERA-5 reference\n",
+ "reference = bias_yearly_avg_ref.mean(dim='year').to_array().squeeze()\n",
+ "im1 = axes[0].imshow(reference.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)\n",
+ "axes[0].text(x=reference.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(reference.mean().item()), fontsize=14, ha='right')\n",
+ "\n",
+ "# plot MAE of model\n",
+ "prediction = bias_yearly_avg.mean(dim='year').to_array().squeeze()\n",
+ "im2 = axes[1].imshow(prediction.values, origin='lower', cmap='RdBu_r', vmin=-40, vmax=40)\n",
+ "axes[1].text(x=reference.shape[0] - 2, y=2, s='Average: {:.1f}%'.format(prediction.mean().item()), fontsize=14, ha='right')\n",
+ "\n",
+ "# plot topography\n",
+ "im_dem = axes[2].imshow(dem['elevation'].values, origin='lower', cmap='terrain', vmin=0, vmax=4000)\n",
+ "\n",
+ "# set titles\n",
+ "axes[0].set_title('ERA-5', fontsize=14, pad=10);\n",
+ "axes[1].set_title('DCEDN', fontsize=14, pad=10);\n",
+ "axes[2].set_title('Copernicus EU-DEM v1.1', fontsize=14, pad=10)\n",
+ "\n",
+ "# adjust axes\n",
+ "for ax in axes.flat:\n",
+ " ax.axes.get_xaxis().set_ticklabels([])\n",
+ " ax.axes.get_xaxis().set_ticks([])\n",
+ " ax.axes.get_yaxis().set_ticklabels([])\n",
+ " ax.axes.get_yaxis().set_ticks([])\n",
+ " ax.axes.axis('tight')\n",
+ " ax.set_xlabel('')\n",
+ " ax.set_ylabel('')\n",
+ " ax.set_axis_off()\n",
+ "\n",
+ "# adjust figure\n",
+ "# fig.suptitle('Average yearly mean absolute error: 1991 - 2010', fontsize=20);\n",
+ "fig.subplots_adjust(hspace=0, wspace=0, top=0.85)\n",
+ "\n",
+ "# add colorbar for bias\n",
+ "axes = axes.flatten()\n",
+ "cbar_ax_bias = fig.add_axes([axes[-1].get_position().x1 + 0.01, axes[-1].get_position().y0,\n",
+ " 0.01, axes[-1].get_position().y1 - axes[-1].get_position().y0])\n",
+ "cbar_bias = fig.colorbar(im_dem, cax=cbar_ax_bias)\n",
+ "cbar_bias.set_label(label='Elevation (m)', fontsize=14)\n",
+ "cbar_bias.ax.tick_params(labelsize=14, pad=10)\n",
+ "\n",
+ "# add colorbar for predictand\n",
+ "cbar_ax_predictand = fig.add_axes([axes[0].get_position().x0, axes[0].get_position().y0 - 0.1,\n",
+ " axes[-1].get_position().x0 - axes[0].get_position().x0,\n",
+ " 0.03])\n",
+ "cbar_predictand = fig.colorbar(im1, cax=cbar_ax_predictand, orientation='horizontal')\n",
+ "cbar_predictand.set_label(label='Relative mean error (%)', fontsize=14)\n",
+ "cbar_predictand.ax.tick_params(labelsize=14, pad=10)\n",
+ "\n",
+ "# add metrics: MAE and RMSE\n",
+ "#axes[1].text(x=ds.shape[0] - 2, y=2, s='MAE = {:.1f}'.format(mae_avg[NAMES[PREDICTAND]].item()) + 'mm day$^{-1}$', fontsize=14, ha='right')\n",
+ "#axes[1].text(x=ds.shape[0] - 2, y=12, s='RMSE = {:.1f}'.format(rmse_avg[NAMES[PREDICTAND]].item()) + 'mm$^2$ day$^{-2}$', fontsize=14, ha='right')\n",
+ "\n",
+ "# save figure\n",
+ "fig.savefig('../Notebooks/Figures/{}_rbias_ERA_vs_model.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
@@ -409,7 +552,7 @@
"source": [
"# plot average of observation, prediction, and bias\n",
"vmin, vmax = 0, 5\n",
- "fig, axes = plt.subplots(len(y_pred_yearly_avg.data_vars), 3, figsize=(24, len(y_pred_yearly_avg.data_vars) * 6),\n",
+ "fig, axes = plt.subplots(len(y_pred_yearly_avg.data_vars), 3, figsize=(24, len(y_pred_yearly_avg.data_vars) * 8),\n",
" sharex=True, sharey=True)\n",
"axes = axes.reshape(len(y_pred_yearly_avg.data_vars), -1)\n",
"for i, var in enumerate(y_pred_yearly_avg):\n",
@@ -457,8 +600,8 @@
"cbar_predictand.ax.tick_params(labelsize=14)\n",
"\n",
"# add metrics: MAE and RMSE\n",
- "axes[1].text(x=ds.shape[0] - 2, y=2, s='MAE = {:.1f}'.format(mae_avg[NAMES[PREDICTAND]].item()) + 'mm day$^{-1}$', fontsize=14, ha='right')\n",
- "axes[1].text(x=ds.shape[0] - 2, y=12, s='RMSE = {:.1f}'.format(rmse_avg[NAMES[PREDICTAND]].item()) + 'mm$^2$ day$^{-2}$', fontsize=14, ha='right')\n",
+ "axes[1].text(x=ds.shape[0] - 2, y=2, s='MAE = {:.1f}'.format(mae_avg[NAMES[PREDICTAND]].mean().item()) + 'mm day$^{-1}$', fontsize=14, ha='right')\n",
+ "axes[1].text(x=ds.shape[0] - 2, y=12, s='RMSE = {:.1f}'.format(rmse_avg[NAMES[PREDICTAND]].mean().item()) + 'mm day$^{-1}$', fontsize=14, ha='right')\n",
"\n",
"# save figure\n",
"fig.savefig('../Notebooks/Figures/{}_average_bias.png'.format(PREDICTAND), dpi=300, bbox_inches='tight')"
@@ -490,7 +633,22 @@
"# 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_pr.groupby('time.season').mean(dim='time')\n",
- "bias_snl = ((y_pred_snl - y_true_snl) / y_true_snl) * 100"
+ "y_refe_snl = y_refe.groupby('time.season').mean(dim='time')\n",
+ "bias_snl = ((y_pred_snl - y_true_snl) / y_true_snl) * 100\n",
+ "bias_snl_ref = ((y_refe_snl - y_true_snl) / y_true_snl) * 100"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2a70e6ac-54fb-47b8-ae76-a1b51d9d1f17",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# print average bias per season: ERA-5\n",
+ "for var in bias_snl_ref.data_vars:\n",
+ " for season in bias_snl_ref[NAMES[PREDICTAND]].season:\n",
+ " print('(ERA-5) Average bias of mean {} for season {}: {:.1f}%'.format(var, season.values.item(), bias_snl_ref[var].sel(season=season).mean().item()))"
]
},
{
@@ -500,10 +658,10 @@
"metadata": {},
"outputs": [],
"source": [
- "# print average bias per season\n",
+ "# print average bias per season: model\n",
"for var in bias_snl.data_vars:\n",
" for season in bias_snl[NAMES[PREDICTAND]].season:\n",
- " print('Average bias of mean {} for season {}: {:.1f}%'.format(var, season.values.item(), bias_snl[var].sel(season=season).mean().item()))"
+ " print('(Model) Average bias of mean {} for season {}: {:.1f}%'.format(var, season.values.item(), bias_snl[var].sel(season=season).mean().item()))"
]
},
{
@@ -597,7 +755,8 @@
"with warnings.catch_warnings():\n",
" warnings.simplefilter('ignore', category=RuntimeWarning)\n",
" y_pred_ex = y_pred_pr.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')"
]
},
{
@@ -609,8 +768,31 @@
"source": [
"# calculate bias in extreme quantile for each year\n",
"bias_ex = ((y_pred_ex - y_true_ex) / y_true_ex) * 100\n",
+ "bias_ex_ref = ((y_refe_ex - y_true_ex) / y_true_ex) * 100"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9a610a70-bb94-47e6-b470-e01bf4a295c0",
+ "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}%'.format(quantile * 100, var, bias_ex_ref[var].mean().item()))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "423962d7-453d-4c86-b93e-a4571bfec1c4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# bias of extreme quantile: Model\n",
"for var in bias_ex:\n",
- " print('Yearly average bias for P{:.0f} of {}: {:.1f}%'.format(quantile * 100, var, bias_ex[var].mean().item()))"
+ " print('(Model) Yearly average bias for P{:.0f} of {}: {:.1f}%'.format(quantile * 100, var, bias_ex[var].mean().item()))"
]
},
{
@@ -622,8 +804,31 @@
"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": "7331acdd-60bb-414e-a0e7-e81956c3a1bf",
+ "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} mm / day'.format(quantile * 100, var, mae_ex_ref[var].item()))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d7bf1a6c-e111-4c72-8774-64a3939cbe50",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# mae of extreme quantile: Model\n",
"for var in mae_ex:\n",
- " print('Yearly average MAE for P{:.0f} of {}: {:.1f} mm / day'.format(quantile * 100, var, mae_ex[var].item()))"
+ " print('(Model) Yearly average MAE for P{:.0f} of {}: {:.1f} mm / day'.format(quantile * 100, var, mae_ex[var].item()))"
]
},
{
@@ -633,10 +838,33 @@
"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": "712a2109-14b9-45a3-9194-e2b917c5ba3f",
+ "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} mm / day'.format(quantile * 100, var, rmse_ex_ref[var].item()))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1dfad687-c16f-4767-8ddd-96c3e56bd96a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# rmse of extreme quantile: Model\n",
"for var in rmse_ex:\n",
- " print('Yearly average RMSE for P{:.0f} of {}: {:.1f} mm / day'.format(quantile * 100, var, rmse_ex[var].item()))"
+ " print('(Model) Yearly average RMSE for P{:.0f} of {}: {:.1f} mm / day'.format(quantile * 100, var, rmse_ex[var].item()))"
]
},
{
@@ -722,7 +950,8 @@
"with warnings.catch_warnings():\n",
" warnings.simplefilter('ignore', category=RuntimeWarning)\n",
" y_true_ex_snl = y_true.groupby('time.season').quantile(quantile, dim='time')\n",
- " y_pred_ex_snl = y_pred_pr.groupby('time.season').quantile(quantile, dim='time')"
+ " y_pred_ex_snl = y_pred_pr.groupby('time.season').quantile(quantile, dim='time')\n",
+ " y_refe_ex_snl = y_refe.groupby('time.season').quantile(quantile, dim='time')"
]
},
{
@@ -733,7 +962,8 @@
"outputs": [],
"source": [
"# compute relative bias in seasonal extremes\n",
- "bias_ex_snl = ((y_pred_ex_snl - y_true_ex_snl) / y_true_ex_snl) * 100"
+ "bias_ex_snl = ((y_pred_ex_snl - y_true_ex_snl) / y_true_ex_snl) * 100\n",
+ "bias_ex_snl_ref = ((y_refe_ex_snl - y_true_ex_snl) / y_true_ex_snl) * 100"
]
},
{
@@ -743,10 +973,23 @@
"metadata": {},
"outputs": [],
"source": [
- "# print average bias in extreme per season\n",
+ "# print average bias in extreme per season: ERA-5\n",
+ "for var in bias_ex_snl_ref.data_vars:\n",
+ " for season in bias_ex_snl_ref[NAMES[PREDICTAND]].season:\n",
+ " print('(ERA-5) Average bias of P{:.0f} {} for season {}: {:.1f}%'.format(quantile * 100, var, season.values.item(), bias_ex_snl_ref[var].sel(season=season).mean().item()))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0671b4f3-52c6-4167-b35b-03902dbe11a3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# print average bias in extreme per season: Model\n",
"for var in bias_ex_snl.data_vars:\n",
" for season in bias_ex_snl[NAMES[PREDICTAND]].season:\n",
- " print('Average bias of P{:.0f} {} for season {}: {:.1f}%'.format(quantile * 100, var, season.values.item(), bias_ex_snl[var].sel(season=season).mean().item()))"
+ " print('(Model) Average bias of P{:.0f} {} for season {}: {:.1f}%'.format(quantile * 100, var, season.values.item(), bias_ex_snl[var].sel(season=season).mean().item()))"
]
},
{
--
GitLab