IOD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 0. Import needed Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the 3 main libraries\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "from eofs.xarray import Eof as eof\n",
    "# librairies to do nicer plots...\n",
    "from matplotlib.patches import Rectangle    # only to draw a rectangle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Get data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sstproduct = 'erssst'\n",
    "#sstproduct = 'oisst025'\n",
    "sstproduct = 'oisst'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if sstproduct == \"erssst\":\n",
    "    \n",
    "    data = xr.tutorial.open_dataset('ersstv5')\n",
    "    data = data.sel(lat=slice(60., -60.))          # get rid of data out of 60S-60N\n",
    "    data = data.sel(time=slice('1970','2021'))     # select years between 1970 and 2021\n",
    "    mask = data.sst.isel(time=0)\n",
    "    mask.data = np.where(np.isnan(mask.data), 0., 1.)\n",
    "    \n",
    "elif sstproduct == \"oisst\":\n",
    "    \n",
    "    data = xr.open_dataset('sst.mnmean.nc')\n",
    "    data = data.sel(time=slice('1982','2022'))     # select years between 1982 and 2022\n",
    "    data_mask = xr.open_dataset('lsmask.nc')\n",
    "    mask = data_mask.mask                          # get the DataArray from the Dataset\n",
    "    mask = mask.squeeze()                          # get rid of the time dimension in the variable mask\n",
    "    \n",
    "elif sstproduct == \"oisst025\":\n",
    "    \n",
    "    data = xr.open_dataset('sst.mon.mean.nc')\n",
    "    data = data.sel(time=slice('1981-11','2023-11'))     # select years between 1982 and 2022\n",
    "    mask = data.sst.isel(time=0)\n",
    "    mask.data = np.where(np.isnan(mask.data), 0., 1.)    \n",
    "else:\n",
    "\n",
    "    print(\"wrong name of SST product\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = np.cos( np.deg2rad(mask.lat) )\n",
    "weights = mask*weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Detrend SST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linfit = data.sst.polyfit('time', 1)\n",
    "trend = xr.polyval(coord=data.time, coeffs=linfit.polyfit_coefficients)   # SST trend\n",
    "sst_detrend = data.sst - trend.values + data.sst.mean(dim='time')         # detrended SST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Detrended interannual anomaly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sstbymth = sst_detrend.groupby(\"time.month\")\n",
    "mthclim = sstbymth.mean(\"time\")       # detrended climatological months (seasonal cycle)\n",
    "sstanom = sstbymth - mthclim          # detrended interannual annomaly\n",
    "sstanom = sstanom.rename('sstanom')   # change variable name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute nino3.4 (5S-5N and 170-120W) index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if sstproduct == \"oisst025\":\n",
    "    latslice = slice(-4.5,4.5)\n",
    "else:\n",
    "    latslice = slice(4.5,-4.5)\n",
    "\n",
    "weights_nino34 = weights.sel(lon=slice(360-169.5,360-120.5), lat=latslice)\n",
    "nino34_index = sstanom.weighted(weights_nino34).mean(dim=('lon','lat'))\n",
    "nino34_index = nino34_index.rename('nino34')   # change variable name\n",
    "nino34_index.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### cut over Indian Ocean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if sstproduct == \"oisst025\":\n",
    "    latslice = slice(-20, 20)\n",
    "else:\n",
    "    latslice = slice(20, -20)\n",
    "\n",
    "data = data.sel(lon=slice(30,120),lat=latslice)\n",
    "sstanom = sstanom.sel(lon=slice(30,120),lat=latslice)\n",
    "mask = mask.sel(lon=slice(30,120),lat=latslice)\n",
    "weights = np.cos( np.deg2rad(mask.lat) )\n",
    "weights = mask*weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SST anomalies regressed onto Nino3.4 SST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sstanom.coords[\"time\"] = ((\"time\"), nino34_index.data)  # redefine 'time' coordinates with nino34 time series\n",
    "linfit = sstanom.polyfit('time', 1)                     # compute the regression with this new \"time\" \n",
    "sstanom.coords[\"time\"] = ((\"time\"), data.time.data)     # put back original time\n",
    "\n",
    "sstreg = linfit.polyfit_coefficients.isel(degree=0)\n",
    "sstcor = xr.corr(sstanom, nino34_index, dim='time')                      \n",
    "\n",
    "#  create a 2 pannels figure\n",
    "fig, axes = plt.subplots(1,2,figsize=(14, 4))\n",
    "sstreg.where(mask == 1.).plot(ax=axes[0])                                # plot the regression coefficient\n",
    "rect = Rectangle((360-170,-5),50,10,linewidth=1,edgecolor='w',fill=False)   # nino3.4 rectangle (170W-120W, 5S-5N)\n",
    "axes[0].add_patch(rect)                                                     # overlay \n",
    "axes[0].set_title('Regression coefficient')\n",
    "sstcor.where(mask == 1.).plot(ax=axes[1])                                # plot the correlation\n",
    "rect = Rectangle((360-170,-5),50,10,linewidth=1,edgecolor='w',fill=False)   # nino3.4 rectangle (170W-120W, 5S-5N)\n",
    "axes[1].add_patch(rect)                                                     # overlay \n",
    "axes[1].set_title('Correlation')\n",
    "axes[1].set_xlabel('lon')\n",
    "axes[1].set_ylabel('lat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. EOF of SST anomalies:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create an [Eof](https://ajdawson.github.io/eofs/latest/api/eofs.xarray.html#eofs.xarray.Eof) object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = eof(sstanom, weights=weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explained variance by the different EOFs/PCs couple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver.varianceFraction(neigs=10).plot.step(where='mid')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Plot the first 2 EOFs, PCs and Variance Fraction:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eof = linear_regression_coef(PC,sstanom) = Cov(PC,sstanom)/Var(PC)\n",
    "# pcscaling=1 (default) --> Var(PC) = 1 --> eof = Cov(PC,sstanom)\n",
    "eofs = solver.eofsAsCovariance(neofs=2, pcscaling=1)      # beware of syntaxe: A C in eofsAsCovariance\n",
    "pcs  = solver.pcs(npcs=2, pcscaling=1)\n",
    "varfrac = solver.varianceFraction(neigs=2)                # beware of syntaxe: F in varianceFraction\n",
    "\n",
    "# create a 4 pannels figure\n",
    "fig, axes = plt.subplots(2,2,figsize=(10, 5),constrained_layout=True)\n",
    "# plot EOF1 and PC1\n",
    "eofs.sel(mode=0).plot(ax=axes[0,0], cbar_kwargs={'label': '°C for 1 std of PC'})\n",
    "axes[0,0].set_title('EOF 1: '+str(int(varfrac.values[0]*100))+'%') \n",
    "pcs.sel(mode=0).plot(ax=axes[1,0])\n",
    "axes[1,0].set_title('PC 1')\n",
    "# plot EOF2 and PC2\n",
    "eofs.sel(mode=1).plot(ax=axes[0,1], cbar_kwargs={'label': '°C for 1 std of PC'})\n",
    "axes[0,1].set_title('EOF 2: '+str(int(varfrac.values[1]*100))+'%')\n",
    "pcs.sel(mode=1).plot(ax=axes[1,1])\n",
    "axes[1,1].set_title('PC 2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(solver.pcs(npcs=10, pcscaling=1).var(dim='time'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcs  = solver.pcs(npcs=2, pcscaling=1)\n",
    "xr.corr(pcs, nino34_index, dim='time').data                    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Indian Ocean Dipole indices\n",
    "    IOD west: 50°E to 70°E and 10°S to 10°N\n",
    "    IOD east: 90°E to 110°E and 10°S to 0°S\n",
    "IOD index = IOD west - IOD east"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if sstproduct == \"oisst025\":\n",
    "    latslice = slice(-10,10)\n",
    "else:\n",
    "    latslice = slice(10,-10)\n",
    "weights_iodw = weights.sel(lon=slice(50,70), lat=latslice)\n",
    "iodw = sstanom.weighted(weights_iodw).mean(dim=('lon','lat'))\n",
    "\n",
    "if sstproduct == \"oisst025\":\n",
    "    latslice = slice(-10,0)\n",
    "else:\n",
    "    latslice = slice(0,-10)\n",
    "weights_iode = weights.sel(lon=slice(90,110), lat=latslice)\n",
    "iode = sstanom.weighted(weights_iode).mean(dim=('lon','lat'))\n",
    "iod_index = iodw - iode"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### IOD and Nino34 plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,figsize=(10, 3))\n",
    "iod_index.plot(ax=axes, label='IOD')\n",
    "nino34_index.plot(ax=axes, label='Nino34')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### correlation of IOD with EOF2 and Nino34"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Correlation IOD-EOF2  :', xr.corr(iod_index, pcs.sel(mode=1),dim='time').data)\n",
    "print('Correlation IOD-Nino34:', xr.corr(iod_index, nino34_index,dim='time').data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### StdDev of IOD and Nino34"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iod_std = iod_index.std(dim='time')\n",
    "nino34_std = nino34_index.std(dim='time')\n",
    "print('IOD StdDev   :', iod_std.data)\n",
    "print('Nino34 StdDev:', nino34_std.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monthly std of IOD index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iod_index.groupby('time.month').std(dim='time').plot.step(where='mid')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## relation with the precipitation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data2 = xr.open_dataset('precip.mon.mean.nc')\n",
    "data2 = data2.sel(lat=slice(60., -60.))          # get rid of data out of 60S-60N\n",
    "if sstproduct == \"oisst025\":\n",
    "    data2 = data2.sel(time=slice('1981-11','2023-11'))     \n",
    "else:\n",
    "    data2 = data2.sel(time=slice('1982','2022'))     \n",
    "data2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure( figsize=(10, 3) )\n",
    "crspc = ccrs.PlateCarree(central_longitude=200)\n",
    "ax = plt.axes(projection=crspc)\n",
    "crspc = ccrs.PlateCarree(central_longitude=0)\n",
    "data2.precip.mean(dim='time').plot(ax=ax, transform=crspc, cmap='YlGnBu')\n",
    "ax.gridlines(draw_labels = True)\n",
    "ax.coastlines()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linfit = data2.precip.polyfit('time', 1)\n",
    "trend = xr.polyval(coord=data2.time, coeffs=linfit.polyfit_coefficients)   # precip trend\n",
    "precip_detrend = data2.precip - trend.values + data2.precip.mean(dim='time') # detrended precip\n",
    "precipbymth = precip_detrend.groupby(\"time.month\")\n",
    "mthclim = precipbymth.mean(\"time\")             # detrended climatological months (seasonal cycle)\n",
    "precipanom = precipbymth - mthclim             # detrended interannual annomaly\n",
    "precipanom = precipanom.rename('precipanom')   # change variable name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure( figsize=(10, 3) )\n",
    "crspc = ccrs.PlateCarree(central_longitude=200)\n",
    "ax = plt.axes(projection=crspc)\n",
    "crspc = ccrs.PlateCarree(central_longitude=0)\n",
    "precipanom.isel(time=0).plot(ax=ax, transform=crspc, cmap='seismic')\n",
    "ax.gridlines(draw_labels = True)\n",
    "ax.coastlines()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precip vs. Nino3.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "precipanom.coords[\"time\"] = ((\"time\"), nino34_index.data)  # redefine 'time' coordinates with nino34 time series\n",
    "linfit = precipanom.polyfit('time', 1)                     # compute the regression with this new \"time\" \n",
    "precipanom.coords[\"time\"] = ((\"time\"), data2.time.data)     # put back original time\n",
    "\n",
    "precipreg_nino = linfit.polyfit_coefficients.isel(degree=0)\n",
    "precipcor_nino = xr.corr(precipanom, nino34_index, dim='time')                      \n",
    "\n",
    "#  create a 2 pannels figure\n",
    "crspc = ccrs.PlateCarree(central_longitude=200)\n",
    "fig, axes = plt.subplots(2,1,figsize=(10, 5), subplot_kw=dict(projection=crspc))\n",
    "\n",
    "crspc = ccrs.PlateCarree(central_longitude=0)\n",
    "# plot the regression coefficient\n",
    "precipreg_nino.plot(ax=axes[0], transform=crspc, cmap='seismic')\n",
    "axes[0].set_title('Regression coefficient')\n",
    "axes[0].coastlines()\n",
    "# plot the correlation\n",
    "precipcor_nino.plot(ax=axes[1], transform=crspc, cmap='seismic')\n",
    "axes[1].set_title('Correlation')\n",
    "axes[1].coastlines()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precip vs. IOD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "precipanom.coords[\"time\"] = ((\"time\"), iod_index.data)  # redefine 'time' coordinates with iod time series\n",
    "linfit = precipanom.polyfit('time', 1)                     # compute the regression with this new \"time\" \n",
    "precipanom.coords[\"time\"] = ((\"time\"), data2.time.data)     # put back original time\n",
    "\n",
    "precipreg_iod = linfit.polyfit_coefficients.isel(degree=0)\n",
    "precipcor_iod = xr.corr(precipanom, iod_index, dim='time')                      \n",
    "\n",
    "#  create a 2 pannels figure\n",
    "crspc = ccrs.PlateCarree(central_longitude=200)\n",
    "fig, axes = plt.subplots(2,1,figsize=(10, 5), subplot_kw=dict(projection=crspc))\n",
    "\n",
    "crspc = ccrs.PlateCarree(central_longitude=0)\n",
    "# plot the regression coefficient\n",
    "precipreg_iod.plot(ax=axes[0], transform=crspc, cmap='seismic')\n",
    "axes[0].set_title('Regression coefficient')\n",
    "axes[0].coastlines()\n",
    "# plot the correlation\n",
    "precipcor_iod.plot(ax=axes[1], transform=crspc, cmap='seismic')\n",
    "axes[1].set_title('Correlation')\n",
    "axes[1].coastlines()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  create a 2 pannels figure\n",
    "crspc = ccrs.PlateCarree(central_longitude=200)\n",
    "fig, axes = plt.subplots(2,1,figsize=(10, 5), subplot_kw=dict(projection=crspc))\n",
    "\n",
    "crspc = ccrs.PlateCarree(central_longitude=0)\n",
    "# plot the regression coefficient\n",
    "(precipreg_nino*nino34_std).plot(ax=axes[0], transform=crspc, cmap='seismic',vmin=-2.5,vmax=2.5)\n",
    "axes[0].set_title('Precip/Nino Regression * Nino StdDev')\n",
    "axes[0].coastlines()\n",
    "# plot the correlation\n",
    "(precipreg_iod*iod_std).plot(ax=axes[1], transform=crspc, cmap='seismic',vmin=-2.5,vmax=2.5)\n",
    "axes[1].set_title('Precip/IOD Regression * IOD StdDev')\n",
    "axes[1].coastlines()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precip anom for pure IOD, Nino and Mx IOD-Nino"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crspc = ccrs.PlateCarree(central_longitude=200)\n",
    "fig, axes = plt.subplots(3,1,figsize=(10, 7), subplot_kw=dict(projection=crspc))\n",
    "\n",
    "crspc = ccrs.PlateCarree(central_longitude=0)\n",
    "# precipanom of pure IOD events\n",
    "precip_iod = precipanom.where( (iod_index > 2*iod_std) & (nino34_index < 0.5*nino34_std) ).mean(dim='time')\n",
    "precip_iod.plot(ax=axes[0], transform=crspc, cmap='seismic',vmin=-10,vmax=10)\n",
    "axes[0].set_title('\"Pure IOD\" precip')\n",
    "axes[0].coastlines()\n",
    "# precipanom of pure Nino events\n",
    "precip_nino = precipanom.where( (iod_index < 0.5*iod_std) & (nino34_index > 2*nino34_std) ).mean(dim='time')\n",
    "precip_nino.plot(ax=axes[1], transform=crspc, cmap='seismic',vmin=-10,vmax=10)\n",
    "axes[1].set_title('\"Pure Nino\" precip')\n",
    "axes[1].coastlines()\n",
    "# precipanom of IOD-Nino events\n",
    "precip_iodnino = precipanom.where( (iod_index > 2*iod_std) & (nino34_index > 2*nino34_std) ).mean(dim='time')\n",
    "precip_iodnino.plot(ax=axes[2], transform=crspc, cmap='seismic',vmin=-10,vmax=10)\n",
    "axes[2].set_title('\"IOD+Nino\" precip')\n",
    "axes[2].coastlines()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precip November Mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure( figsize=(10, 3) )\n",
    "crspc = ccrs.PlateCarree(central_longitude=200)\n",
    "ax = plt.axes(projection=crspc)\n",
    "crspc = ccrs.PlateCarree(central_longitude=0)\n",
    "novmean = data2.precip.groupby('time.month').mean(dim='time').sel(month=11)\n",
    "novmean.plot(ax=ax, transform=crspc, cmap='YlGnBu',vmin=0,vmax=10)\n",
    "ax.gridlines(draw_labels = True)\n",
    "ax.coastlines()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ratio precip anom/mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure( figsize=(10, 3) )\n",
    "crspc = ccrs.PlateCarree(central_longitude=200)\n",
    "ax = plt.axes(projection=crspc)\n",
    "crspc = ccrs.PlateCarree(central_longitude=0)\n",
    "(precip_iodnino /novmean).where(novmean > 1).plot(ax=ax, transform=crspc, cmap='seismic',vmin=-2,vmax=2)\n",
    "ax.gridlines(draw_labels = True)\n",
    "ax.coastlines()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}