{
"cells": [
{
"cell_type": "markdown",
"id": "8c73bb8b",
"metadata": {},
"source": [
"# Utilities"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "cfe495d7",
"metadata": {
"ExecuteTime": {
"end_time": "2022-01-06T08:32:58.847592Z",
"start_time": "2022-01-06T08:32:57.118727Z"
}
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"from jupyterthemes import jtplot\n",
"import matplotlib as mpl\n",
"from scipy.integrate import odeint\n",
"from numpy import linalg as LA\n",
"from scipy.optimize import fsolve\n",
"from scipy.optimize import curve_fit\n",
"mpl.rcdefaults() \n",
"from scipy.stats import uniform\n",
"import scipy.stats as st\n",
"from matplotlib import colors\n",
"from matplotlib import cm\n",
"import pandas as pd\n",
"import datetime\n",
"from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
"import cartopy.crs as ccrs\n",
"import cartopy\n",
"import json\n",
"import xarray as xr\n",
"import pickle\n",
"from scipy.stats import linregress\n",
"from EWS_functions import *\n",
"from scipy.optimize import curve_fit\n",
"from scipy.ndimage import gaussian_filter1d\n",
"\n",
"jtplot.style(context='paper', fscale=1.4, spines=True, grid=False, ticks=True,gridlines='--')\n",
"\n",
"fontsize=16\n",
"mpl.rcParams['xtick.direction'] = 'in'\n",
"mpl.rcParams['ytick.direction'] = 'in'\n",
"mpl.rcParams['xtick.top'] = True\n",
"mpl.rcParams['ytick.right'] = True\n",
"\n",
"mpl.rcParams['font.size'] = 16\n",
"mpl.rcParams['legend.fontsize'] = 'large'\n",
"mpl.rcParams['figure.titlesize'] = 'medium'\n",
"mpl.rcParams['axes.labelsize']= 'x-large'\n",
"mpl.rcParams['figure.facecolor']='white'\n",
"\n",
"mpl.rcParams['font.family'] = 'sans-serif'\n",
"mpl.rcParams['font.sans-serif'] = ['Arial']\n",
"hfont = {'fontname':'Arial'}\n",
"\n",
"mpl.rcParams['text.latex.preamble']= r'\\usepackage{amsmath}'\n",
"mpl.rcParams['axes.prop_cycle'] = mpl.cycler(color=['#1b9e77','#d95f02','#7570b3','#e7298a','#66a61e','#e6ab02','#a6761d','#666666']) "
]
},
{
"cell_type": "markdown",
"id": "71826856",
"metadata": {},
"source": [
"## EWS plot functions"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "0e524815",
"metadata": {
"ExecuteTime": {
"end_time": "2022-01-06T08:32:58.856694Z",
"start_time": "2022-01-06T08:32:58.849471Z"
},
"init_cell": true
},
"outputs": [],
"source": [
"def get_EWS(time,data,trend,ws):\n",
" linfits = []\n",
" ps = []\n",
" bound = ws // 2\n",
" \n",
" std = runstd(data - trend, ws)[bound:-bound]\n",
" p0, p1 = np.polyfit(time[bound : -bound][:-2], std[:-2], 1)\n",
" linfits.append([p0,p1])\n",
" ps.append(kendall_tau_test(std[:-2], 1000, p0))\n",
" \n",
" ar1 = runac(data - trend, ws)[bound : -bound]\n",
" p0, p1 = np.polyfit(time[bound : -bound][:-3], ar1[:-3], 1)\n",
" linfits.append([p0,p1])\n",
" ps.append(kendall_tau_test(ar1[:-2], 1000, p0))\n",
" \n",
" lam = run_fit_a_ar1(data-trend,ws)[bound:-bound]\n",
" p0, p1 = np.polyfit(time[bound : -bound][:-2], lam[:-2], 1)\n",
" linfits.append([p0,p1])\n",
" ps.append(kendall_tau_test(lam[:-2], 1000, p0))\n",
" \n",
" return std, ar1, lam, linfits, ps"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "9115bb07",
"metadata": {
"ExecuteTime": {
"end_time": "2022-01-06T08:32:58.959523Z",
"start_time": "2022-01-06T08:32:58.857978Z"
},
"init_cell": true
},
"outputs": [],
"source": [
"def plot_EWS(data, timess, ws=70, col='k',lbl='',alph=1,lw=1):\n",
" bound = ws // 2\n",
" popt, cov = curve_fit(funcfit3, timess, data, p0 = [-8.33097773e-01, 1.05507897e-02, 2.02518923e+03], maxfev = 1000000000)\n",
" trend = funcfit3(timess, *popt)\n",
" std, ar1, lam, linfits, ps = get_EWS(timess,data,trend,ws)\n",
"\n",
" ax1.plot(timess[bound : -bound],std,color=col,label=lbl,alpha=alph,lw=lw)\n",
" pv = kendall_tau_test(std[:-2],1000,linfits[1][0])\n",
" ax1.plot(timess[bound : -bound][:-2],linfits[0][0] * timess[bound : -bound][:-2] + linfits[0][1],linestyle='--',color=col,alpha=alph,lw=lw,label=\"p = {:.3f}\".format(pv))\n",
"\n",
" ax2.plot(timess[bound : -bound],ar1,color=col,label=lbl,alpha=alph,lw=lw)\n",
" pv = kendall_tau_test(ar1[:-2],1000,linfits[1][0])\n",
" ax2.plot(timess[bound : -bound][:-2],linfits[1][0] * timess[bound : -bound][:-2] + linfits[1][1],linestyle='--',color=col,alpha=alph,lw=lw,label=\"p = {:.3f}\".format(pv))\n",
"\n",
" ax3.plot(timess[bound : -bound],lam,color=col,label=lbl,alpha=alph,lw=lw)\n",
" p0, p1, p2 = np.polyfit(timess[bound : -bound][:-2], lam[:-2], 2)\n",
" pl0, pl1 = np.polyfit(timess[bound : -bound][:-2], lam[:-2], 1)\n",
" pv = kendall_tau_test(lam[:-2], 1000, pl0) # precentile of 1000 fourier surrogates have a larger linear slope\n",
"# ax3.plot(timess[bound : -bound][:-2], p0 * timess[bound : -bound][:-2]**2+p1 * timess[bound : -bound][:-2] + p2, color=col,linestyle='--',alpha=alph,lw=lw,label=\"p = {:.3f}\".format(pv))\n",
" ax3.plot(timess[bound : -bound][:-2], pl0 * timess[bound : -bound][:-2] + pl1, color=col,linestyle='--',alpha=alph,lw=lw,label=\"p = {:.3f}\".format(pv))"
]
},
{
"cell_type": "markdown",
"id": "7aa8b6d1",
"metadata": {},
"source": [
"# files"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "c4e8b8ac",
"metadata": {
"ExecuteTime": {
"end_time": "2022-01-06T08:32:58.994164Z",
"start_time": "2022-01-06T08:32:58.961716Z"
}
},
"outputs": [],
"source": [
"ds = xr.open_dataset('DAMIP_amoc.nc')\n",
"damip_strn = ds.amoc_damip.sel(latitudes='26.5N')\n",
"damip_index = ds.index_damip"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "549fce18",
"metadata": {
"ExecuteTime": {
"end_time": "2022-01-06T08:32:59.489276Z",
"start_time": "2022-01-06T08:32:58.995726Z"
}
},
"outputs": [],
"source": [
"aer = xr.open_dataarray('hist_aer_amoc_index.nc').groupby('time.year').mean('time')\n",
"nat = xr.open_dataarray('hist_nat_amoc_index.nc').groupby('time.year').mean('time')\n",
"ghg = xr.open_dataarray('hist_ghg_amoc_index.nc').groupby('time.year').mean('time')"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "673e5811",
"metadata": {
"ExecuteTime": {
"end_time": "2022-01-06T08:32:59.504402Z",
"start_time": "2022-01-06T08:32:59.490830Z"
}
},
"outputs": [],
"source": [
"ds = xr.open_dataset('CMIP6_amoc.nc')\n",
"strn26 = ds.strength_265N\n",
"strn35 = ds.strength_35N\n",
"index = ds.index\n",
"hist_cesm2 = index.sel(models='CESM2')"
]
},
{
"cell_type": "markdown",
"id": "e1a1af19",
"metadata": {},
"source": [
"# SST Index"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "49537377",
"metadata": {
"ExecuteTime": {
"end_time": "2022-01-06T08:32:59.510292Z",
"start_time": "2022-01-06T08:32:59.505950Z"
}
},
"outputs": [],
"source": [
"# ws = 70\n",
"# bound = ws // 2\n",
"\n",
"# fig, (ax0, ax1,ax2,ax3) = plt.subplots(4,1,figsize=(15,25))\n",
"\n",
"# i = 0\n",
"# ens_colors = ['#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a']\n",
"\n",
"# amoc = np.nan_to_num(hist_cesm2.isel(ensemble_members=i).values)\n",
"# yrs = hist_cesm2.year.values\n",
"# plot_EWS(amoc,yrs,ws=ws,col='k',lbl='Total',lw=2)\n",
"\n",
"# amoc = np.nan_to_num(ghg.sel(models='CESM2').isel(ensemble_members=i).values)[0]\n",
"# yrs = ghg.year.values\n",
"# plot_EWS(amoc,yrs,ws=ws,col='C0',lbl='Hist-ghg',lw=2)\n",
"\n",
"# amoc = np.nan_to_num(aer.sel(models='CESM2').isel(ensemble_members=i).values)[0]\n",
"# yrs = aer.year.values\n",
"# plot_EWS(amoc,yrs,ws=ws,col='C1',lbl='Hist-aer',lw=2)\n",
"\n",
"# amoc = np.nan_to_num(nat.sel(models='CESM2').isel(ensemble_members=i).values)[0]\n",
"# yrs = nat.year.values\n",
"# plot_EWS(amoc,yrs,ws=ws,col='C2',lbl='Hist-nat',lw=2)\n",
"\n",
"# ax0.set_ylabel('SST Index')\n",
"# hist_cesm2.isel(ensemble_members=i).plot(label='Total',ax=ax0,color='k',lw=2)\n",
"# ghg.isel(ensemble_members=i).plot(label='Hist-ghg',ax=ax0,color='C0',lw=2)\n",
"# aer.isel(ensemble_members=i).plot(label='Hist-aer',ax=ax0,color='C1',lw=2)\n",
"# nat.isel(ensemble_members=i).plot(label='Hist-nat',ax=ax0,color='C2',lw=2)\n",
"# ax0.set_title('')\n",
"# ax0.set_xlabel('')\n",
"# ax0.legend(bbox_to_anchor=(1.05, 0.95), loc='upper left', borderaxespad=0., fontsize=20,fancybox=False,frameon=False)\n",
"\n",
"\n",
"# ax1.set_ylabel('Variance')\n",
"# ax1.legend(bbox_to_anchor=(1.05, 0.95), loc='upper left', borderaxespad=0., fontsize=20,fancybox=False,frameon=False)\n",
"# ax1.tick_params(axis='x', which='both',bottom=True,top=True,labelbottom=False, labelsize=fontsize)\n",
"# ax2.set_ylabel('AR1')\n",
"# ax2.legend(bbox_to_anchor=(1.05, 0.95), loc='upper left', borderaxespad=0., fontsize=20,fancybox=False,frameon=False)\n",
"# ax2.tick_params(axis='x', which='both',bottom=True,top=True,labelbottom=False, labelsize=fontsize)\n",
"# ax3.set_ylabel('$\\lambda$')\n",
"# ax3.legend(bbox_to_anchor=(1.05, 0.95), loc='upper left', borderaxespad=0., fontsize=20,fancybox=False,frameon=False)\n",
"# ax3.set_xlabel('Time [yrs]')\n",
"\n",
"# fig.suptitle('Damip forcing comparisons\\nCESM2 ensemble member 1',fontsize=20, y=0.9)\n",
"# plt.show()"
]
},
{
"cell_type": "markdown",
"id": "1726a778",
"metadata": {},
"source": [
"# AMOC strength"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "c66f679f",
"metadata": {
"ExecuteTime": {
"end_time": "2022-01-06T08:32:59.673347Z",
"start_time": "2022-01-06T08:32:59.511422Z"
}
},
"outputs": [],
"source": [
"with open('matthew/JSON_data/Figure_AR6_DAMIP_AMOC_26N_1000m.json', 'r') as handle:\n",
" json_load = json.load(handle)\n",
"\n",
"amoc_damip6_ts = np.ma.asarray(json_load[\"amoc_damip6_ts\"]) # Note the use of numpy masked arrays (np.ma)\n",
"damip6_models = json_load[\"damip6_models\"]\n",
"year = np.asarray(json_load[\"year\"])\n",
"ens_names = ['r{}i1p1f1'.format(i) for i in range(1,11)]"
]
},
{
"cell_type": "code",
"execution_count": 40,
"id": "8f767d27",
"metadata": {
"ExecuteTime": {
"end_time": "2022-01-06T08:41:15.297580Z",
"start_time": "2022-01-06T08:41:15.272216Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"INDEX calculated:\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" r1i1p1f1 | \n",
" r2i1p1f1 | \n",
" r3i1p1f1 | \n",
" r4i1p1f1 | \n",
" r5i1p1f1 | \n",
" r6i1p1f1 | \n",
" r7i1p1f1 | \n",
" r8i1p1f1 | \n",
" r9i1p1f1 | \n",
" r10i1p1f1 | \n",
"
\n",
" \n",
" \n",
" \n",
" | HadGEM3-GC31-LL | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | CNRM-CM6-1 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | GISS-E2-1-G | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | MIROC6 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
"
\n",
" \n",
" | CanESM5 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
"
\n",
" \n",
" | BCC-CSM2-MR | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" | IPSL-CM6A-LR | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
"
\n",
" \n",
"
\n",
"