Skip to content
Snippets Groups Projects
Commit 1e032657 authored by Frisinghelli Daniel's avatar Frisinghelli Daniel
Browse files

Notebook for Figures of Capstone project.

parent 68635fab
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Notebooks/pr_capstone.ipynb 0 → 100644
+ 313
− 0
View file @ 1e032657
%% Cell type:markdown id:4d872421-be6b-43cd-ac61-158ef7170c0f tags:
### Imports
%% Cell type:code id:e3b75a90-4ff0-4b3e-ad02-b3e838c502aa tags:
``` python
# builtins
import datetime
import warnings
import calendar
# externals
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import scipy.stats as stats
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import scipy.stats as stats
from IPython.display import Image
from sklearn.metrics import r2_score, roc_curve, auc, classification_report
from sklearn.model_selection import train_test_split
# locals
from climax.main.io import ERA5_PATH, OBS_PATH, TARGET_PATH, DEM_PATH
from climax.main.config import CALIB_PERIOD, VALID_PERIOD
from pysegcnn.core.utils import search_files
from pysegcnn.core.graphics import plot_classification_report
```
%% Cell type:code id:741ccf7a-80e1-46bb-926e-a68c14f0c9b2 tags:
``` python
# entire reference period
REFERENCE_PERIOD = np.concatenate([CALIB_PERIOD, VALID_PERIOD], axis=0)
```
%% Cell type:markdown id:efaaf766-e1d3-48a7-8702-7fecde6c8ca3 tags:
### Load observations
%% Cell type:code id:c8049468-d848-4dd3-9082-97ad51e5f5ff tags:
``` python
# model predictions and observations NetCDF
y_true_pr = xr.open_dataset(search_files(OBS_PATH.joinpath('pr'), 'OBS_pr(.*).nc$').pop())
y_true_tmax = xr.open_dataset(search_files(OBS_PATH.joinpath('tasmax'), 'OBS_tasmax(.*).nc$').pop())
y_true_tmin = xr.open_dataset(search_files(OBS_PATH.joinpath('tasmin'), 'OBS_tasmin(.*).nc$').pop())
```
%% Cell type:markdown id:69390143-4c30-443f-9ef5-d1dcde4b2592 tags:
### Load ERA-5 reference dataset
%% Cell type:code id:f1481f20-2d91-429c-b38d-a6f74f97072d tags:
``` python
# search ERA-5 reference dataset
y_refe_pr = xr.open_dataset(search_files(ERA5_PATH.joinpath('ERA5', 'total_precipitation'), '.nc$').pop())
y_refe_tmax = xr.open_dataset(search_files(ERA5_PATH.joinpath('ERA5', '2m_{}_temperature'.format('max')), '.nc$').pop())
y_refe_tmin = xr.open_dataset(search_files(ERA5_PATH.joinpath('ERA5', '2m_{}_temperature'.format('min')), '.nc$').pop())
```
%% Cell type:code id:45e11517-197d-4b9c-bd9f-9f96a11ded98 tags:
``` python
# convert to °C
y_refe_tmax = y_refe_tmax - 273.15
y_refe_tmin = y_refe_tmin - 273.15
```
%% Cell type:markdown id:04315e91-2099-4ca7-a324-7173bdcf7750 tags:
### Select time period
%% Cell type:code id:437711c9-94db-4f02-82d1-8dabd109931b tags:
``` python
# time period
PERIOD = REFERENCE_PERIOD
```
%% Cell type:code id:e9b7df28-7490-4d2f-99e3-660f1dc6e99f tags:
``` python
# subset observations to time period
y_true_pr = y_true_pr.sel(time=PERIOD).precipitation
y_true_tmax = y_true_tmax.sel(time=PERIOD).tasmax
y_true_tmin = y_true_tmin.sel(time=PERIOD).tasmin
```
%% Cell type:code id:95b5b7d6-63ea-4cb9-955f-4df9e0b53789 tags:
``` python
# subset Era-5 to time period
y_refe_pr = y_refe_pr.sel(time=PERIOD).drop_vars('lambert_azimuthal_equal_area').rename({'tp': 'precipitation'}).precipitation
y_refe_tmax = y_refe_tmax.sel(time=PERIOD).drop_vars('lambert_azimuthal_equal_area').rename({'t2m': 'tasmax'}).tasmax
y_refe_tmin = y_refe_tmin.sel(time=PERIOD).drop_vars('lambert_azimuthal_equal_area').rename({'t2m': 'tasmin'}).tasmin
```
%% Cell type:markdown id:98d1b503-1f1a-4a8e-a8e4-eac1e045f320 tags:
## Align datasets
%% Cell type:code id:ef839050-2355-4da2-ad19-1e43c4b46414 tags:
``` python
# precipitation
y_true_pr, y_refe_pr = xr.align(y_true_pr, y_refe_pr, join='override')
y_refe_pr = y_refe_pr.where(~np.isnan(y_true_pr), other=np.nan)
```
%% Cell type:code id:f628cf0c-5465-48e2-9d83-89ec407783e7 tags:
``` python
# tasmax
y_true_tmax, y_refe_tmax = xr.align(y_true_tmax, y_refe_tmax, join='override')
y_refe_tmax = y_refe_tmax.where(~np.isnan(y_true_tmax), other=np.nan)
```
%% Cell type:code id:0eaa502f-5865-402f-810b-664c1275999f tags:
``` python
# tasmin
y_true_tmin, y_refe_tmin = xr.align(y_true_tmin, y_refe_tmin, join='override')
y_refe_tmin = y_refe_tmin.where(~np.isnan(y_true_tmin), other=np.nan)
```
%% Cell type:markdown id:9045be71-d0ef-4457-9655-68b1f85d5cfe tags:
### Plot ERA-5 vs. Observed
%% Cell type:code id:58fa1934-c9fc-4e42-9f2b-942a5ccf5617 tags:
``` python
y_refe_values = y_refe_pr.resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim='time')
y_true_values = y_true_pr.resample(time='1M').sum(skipna=False).groupby('time.month').mean(dim='time')
```
%% Cell type:code id:40be10b7-c868-4cce-b1fb-d201366c17c0 tags:
``` python
bias_pr = ((y_refe_values - y_true_values) / y_true_values) * 100
```
%% Cell type:code id:c52eadeb-1e46-41e6-a7b7-9178f3d65536 tags:
``` python
# plot average of observation and reference
vmin, vmax = 0, 150
fig, axes = plt.subplots(1, 3, figsize=(24, 8), sharex=True, sharey=True)
axes = axes.flatten()
# plot Era-5 reanalysis
im1 = axes[0].imshow(y_refe_values.mean(dim='month'), origin='lower', cmap='viridis_r', vmin=vmin, vmax=vmax)
im2 = axes[1].imshow(y_true_values.mean(dim='month'), origin='lower', cmap='viridis_r', vmin=vmin, vmax=vmax)
im3 = axes[2].imshow(bias_pr.mean(dim='month'), origin='lower', cmap='RdBu_r', vmin=-60, vmax=60)
axes[0].set_title('ERA-5 reanalysis', fontsize=16, pad=10);
axes[1].set_title('Observations', fontsize=16, pad=10);
axes[2].set_title('Bias: ERA-5 - Observations', fontsize=16, pad=10)
# adjust axes
for ax in axes.flat:
ax.axes.get_xaxis().set_ticklabels([])
ax.axes.get_xaxis().set_ticks([])
ax.axes.get_yaxis().set_ticklabels([])
ax.axes.get_yaxis().set_ticks([])
ax.axes.axis('tight')
ax.set_xlabel('')
ax.set_ylabel('')
ax.set_axis_off()
# adjust figure
fig.suptitle('Average monthly precipitation (mm): 1980 - 2010', fontsize=20);
fig.subplots_adjust(hspace=0, wspace=0, top=0.85)
# add colorbar for bias
axes = axes.flatten()
cbar_ax_bias = fig.add_axes([axes[-1].get_position().x1 + 0.01, axes[-1].get_position().y0,
0.01, axes[-1].get_position().y1 - axes[-1].get_position().y0])
cbar_bias = fig.colorbar(im3, cax=cbar_ax_bias)
cbar_bias.set_label(label='Relative bias / (%)', fontsize=16)
cbar_bias.ax.tick_params(labelsize=14)
# add colorbar for predictand
cbar_ax_predictand = fig.add_axes([axes[0].get_position().x0, axes[0].get_position().y0 - 0.1,
axes[-1].get_position().x0 - axes[0].get_position().x0,
0.05])
cbar_predictand = fig.colorbar(im1, cax=cbar_ax_predictand, orientation='horizontal')
cbar_predictand.set_label(label='Precipitation / (mm month$^{-1}$)', fontsize=16)
cbar_predictand.ax.tick_params(labelsize=14)
# save figure
fig.savefig('../Notebooks/Figures/capstone_pr.png', dpi=300, bbox_inches='tight')
```
%% Cell type:code id:dc1a78dd-4912-4fc2-ab76-1a878f5e1b4c tags:
``` python
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment