Skip to content
Snippets Groups Projects
Commit b4011a1d authored by Gieße, Dr. Céline's avatar Gieße, Dr. Céline
Browse files

clean up code (sea-ice)

parent 84d5f9fb
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 24 additions and 32 deletions
Sea-ice.ipynb
+ 24
− 32
View file @ b4011a1d
%% Cell type:code id:fb7d6d3e-7726-49d2-a8c5-85bad5854f10 tags:
``` python
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import pandas as pd
import math
from sklearn.linear_model import LinearRegression
from matplotlib.patches import Circle
from matplotlib import cm, colors
# Set font sizes for plots
SMALL_SIZE = 12
MEDIUM_SIZE = 14
BIGGER_SIZE = 16
plt.rc('font', size=MEDIUM_SIZE) # controls default text sizes
plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=MEDIUM_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
plt.rc('font', size=MEDIUM_SIZE)
# Global warming levels
Tglob = ["T0","T15","T27"]
Tglob_values = [0,1.5,2.7]
# Data path
data_path = '/work/uo1075/u241297/data/large_ensembles/CMIP6/'
```
%% Cell type:markdown id:6528eeee-cbb0-44f2-9621-1fb397b2fb3d tags:
# Global warming levels
%% Cell type:code id:c60ff72f-9206-4cf6-84fc-d7e67828a052 tags:
``` python
path=data_path+'MPI-ESM1-2-LR/global-mean/'
# Load global-mean surface air temperature (GSAT) data
GSAT = np.zeros((50, 251))
for i in range(50):
file = path + "tas_hist_ssp585_r" + str(i + 1) + "_global-mean.nc"
GSAT[i, :] = np.squeeze(xr.open_dataset(file)['tas']) - 273.15
GSAT_rm = np.array(pd.DataFrame(GSAT).rolling(20, axis=1, center=True).mean())
Tzero_GSAT = np.mean(GSAT[:, 0:51]) # 1850-1900
GSATa = GSAT - Tzero_GSAT
GSATa_rm = GSAT_rm - Tzero_GSAT
```
%% Cell type:code id:6fa849af-5817-47d9-a15a-76d6f1d82b00 tags:
``` python
# Determine the indices where GSAT crosses the global warming levels
GWLcross_ind = {T: np.empty(50, dtype=int) for T in Tglob}
GWLcross_ind["T0"] = 50 * [20]
for Ti in range(1, 3):
for T in range(1, 3):
for i in range(50):
j = 20
while GSATa_rm[i, j] < Tglob_values[Ti]:
while GSATa_rm[i, j] < Tglob_values[T]:
j += 1
GWLcross_ind[Tglob[Ti]][i] = j
GWLcross_ind[Tglob[T]][i] = j
GWL_ind = {T: np.zeros((50, 20), dtype=int) for T in Tglob}
for Ti in range(3):
for T in range(3):
for i in range(50):
GWL_ind[Tglob[Ti]][i, :] = np.arange(GWLcross_ind[Tglob[Ti]][i] - 10,
GWLcross_ind[Tglob[Ti]][i] + 10)
GWL_ind[Tglob[T]][i, :] = np.arange(GWLcross_ind[Tglob[T]][i] - 10,
GWLcross_ind[Tglob[T]][i] + 10)
```
%% Cell type:markdown id:a9134a75-ce3e-4062-b4ef-7652ba0b5a34 tags:
# Sea ice area
%% Cell type:code id:e99f932a-601d-409f-8fb2-892de7fba8b9 tags:
``` python
# Load Sea Ice Area (SIA) data
path=data_path+'MPI-ESM1-2-LR/sia_nh/'
SIA = np.zeros((50,251,12))
for i in range(50):
file = path+"sia_nh_MPI-ESM1-2-LR_r"+str(i+1)+"i1p1f1.nc"
SIA[i,:] = np.array(xr.open_dataset(file)['sia_nh']).reshape((251,12))
# Create GWL samples of September and March SIA
SIA_sep_samples = {Tglob[Ti]: [SIA[i, GWL_ind[Tglob[Ti]][i, :], 8] for i in range(50)] for Ti in range(3)}
SIA_mar_samples = {Tglob[Ti]: [SIA[i, GWL_ind[Tglob[Ti]][i, :], 2] for i in range(50)] for Ti in range(3)}
SIA_sep_samples = {Tglob[T]: [SIA[i, GWL_ind[Tglob[T]][i, :], 8] for i in range(50)] for T in range(3)}
SIA_mar_samples = {Tglob[T]: [SIA[i, GWL_ind[Tglob[T]][i, :], 2] for i in range(50)] for T in range(3)}
```
%% Cell type:markdown id:a5d79ef0-78e6-49a5-bc79-b56988fafac8 tags:
## Bias correction (based on Niederdrenk and Notz, 2018)
%% Cell type:code id:b7d94506-010a-4378-a4c3-59af4914014e tags:
``` python
# Observational regression parameters (Niederdrenk and Notz, 2018)
slope_obs_sep, intercept_obs_sep = -4.1, 7.97
slope_obs_mar, intercept_obs_mar = -1.6, 15.08
# Model regression line
GSATa_ensmean = np.mean(GSATa, axis=0)
SIA_sep_ensmean = np.mean(SIA[:, :, 8], axis=0)
SIA_mar_ensmean = np.mean(SIA[:, :, 2], axis=0)
model_sep = LinearRegression().fit(GSATa_ensmean[103:167].reshape((-1, 1)), SIA_sep_ensmean[103:167])
model_mar = LinearRegression().fit(GSATa_ensmean[103:167].reshape((-1, 1)), SIA_mar_ensmean[103:167])
# Bias correction function
def bias_correction(SIA, model, slope_obs, intercept_obs, GSATa):
correction = intercept_obs - model.intercept_ + (slope_obs - model.coef_[0]) * GSATa
corrected_SIA = SIA + correction #[:, np.newaxis]
corrected_SIA[corrected_SIA < 0] = 0
return corrected_SIA
# Apply bias correction for September and March
SIA_sep_corr = bias_correction(SIA[:, :, 8], model_sep, slope_obs_sep, intercept_obs_sep, GSATa)
SIA_mar_corr = bias_correction(SIA[:, :, 2], model_mar, slope_obs_mar, intercept_obs_mar, GSATa)
# Bias corrected samples for September and March
SIA_sep_corr_samples = {
Tglob[Ti]: [SIA_sep_corr[i, GWL_ind[Tglob[Ti]][i, :]] for i in range(50)]
for Ti in range(3)
Tglob[T]: [SIA_sep_corr[i, GWL_ind[Tglob[T]][i, :]] for i in range(50)]
for T in range(3)
}
SIA_mar_corr_samples = {
Tglob[Ti]: [SIA_mar_corr[i, GWL_ind[Tglob[Ti]][i, :]] for i in range(50)]
for Ti in range(3)
Tglob[T]: [SIA_mar_corr[i, GWL_ind[Tglob[T]][i, :]] for i in range(50)]
for T in range(3)
}
```
%% Cell type:code id:ece4a709-832a-4316-ab0e-dfa9699f0b76 tags:
``` python
# Define a function for detrending by subtraction of ensemble mean
def detrend_samples(samples, Tglob):
ensmean_samples = {T: np.mean(samples[T], axis=0) for T in Tglob}
detrend_samples = {
T: samples[T] - ensmean_samples[T][np.newaxis, :] + np.mean(ensmean_samples[T])
for T in Tglob
}
for T in Tglob:
detrend_samples[T][detrend_samples[T] < 0] = 0
return detrend_samples
# Detrend by ensemble mean
SIA_sep_corr_detrend_samples = detrend_samples(SIA_sep_corr_samples, Tglob)
SIA_mar_corr_detrend_samples = detrend_samples(SIA_mar_corr_samples, Tglob)
```
%% Cell type:markdown id:8771f460-7cb6-475f-865a-e8c5ac2c2fbc tags:
## Circle plots
%% Cell type:code id:adc41e36-aefc-41be-a496-4a50f4fe47f3 tags:
``` python
# Constants
AO_area = 14.06
center = (0.5, 0.5)
radius_AO = 0.5
# Load the colormap
cmap = cm.get_cmap('Blues', 101)
# Create subplots
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
# Loop over the three GWLs
for Ti in range(3):
for T in range(3):
circle_AO = Circle(center, radius_AO, fc=cmap(0.99))
ax[Ti].add_patch(circle_AO)
ax[T].add_patch(circle_AO)
# Draw sea ice area circles
for perc in range(100, -1, -1):
radius_SIA = radius_AO * math.sqrt(np.percentile(SIA_sep_corr_detrend_samples[Tglob[Ti]], perc) / AO_area)
radius_SIA = radius_AO * math.sqrt(np.percentile(SIA_sep_corr_detrend_samples[Tglob[T]], perc) / AO_area)
circle_SIA = Circle(center, radius_SIA, fc=cmap(perc / 100))
ax[Ti].add_patch(circle_SIA)
ax[T].add_patch(circle_SIA)
# Add a white circle at the end
circle_SIA = Circle(center, radius_SIA, fc='white')
ax[Ti].add_patch(circle_SIA)
ax[T].add_patch(circle_SIA)
ax[Ti].set_aspect('equal')
ax[Ti].axis('off')
ax[T].set_aspect('equal')
ax[T].axis('off')
# Add a colorbar
cax = fig.add_axes([0.94, 0.18, 0.015, 0.64])
norm = colors.Normalize(vmin=0, vmax=100)
sm = plt.cm.ScalarMappable(cmap=cmap.reversed(), norm=norm)
sm.set_array([])
cbar = plt.colorbar(sm, cax=cax)
cbar.set_label('Probability of sea-ice cover [%]')
# Show plot
plt.show()
```
%% Output
%% Cell type:code id:0f9fe742-19db-47fb-9b60-c7b6e26bda4f tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment