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

comments

parent 39a392b9
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 2 additions and 2 deletions
Sea-ice.ipynb
+ 2
− 2
View file @ 08f9f9f5
%% 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 T in range(1, 3):
for i in range(50):
j = 20
while GSATa_rm[i, j] < Tglob_values[T]:
j += 1
GWLcross_ind[Tglob[T]][i] = j
GWL_ind = {T: np.zeros((50, 20), dtype=int) for T in Tglob}
for T in range(3):
for i in range(50):
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[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[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[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 T in range(3):
# Draw outer circle
circle_AO = Circle(center, radius_AO, fc=cmap(0.99))
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[T]], perc) / AO_area)
circle_SIA = Circle(center, radius_SIA, fc=cmap(perc / 100))
ax[T].add_patch(circle_SIA)
# Add a white circle at the end
# Draw inner white circle at minimum SIA
circle_SIA = Circle(center, radius_SIA, fc='white')
ax[T].add_patch(circle_SIA)
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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment