TD2.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",
    "# librairies to do nicer plots...\n",
    "from matplotlib.patches import Rectangle    # only to draw a rectangle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Get data:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select the SST product you whant to use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sstproduct = 'erssst'\n",
    "sstproduct = 'oisst'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### read SST and get land-sea mask"
   ]
  },
  {
   "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",
    "else:\n",
    "\n",
    "    print(\"wrong name of SST product\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute Weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = np.cos( np.deg2rad(mask.lat) )\n",
    "weights = mask*weights\n",
    "weights.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. linear trend:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Long term mean trend \n",
    "[linear regression](http://xarray.pydata.org/en/stable/computation.html?highlight=polyval#fitting-polynomials) with [polyfit](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.polyfit.html?highlight=polyfit) and [polyval](http://xarray.pydata.org/en/stable/generated/xarray.polyval.html?highlight=polyval) methonds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "linfit = data.sst.polyfit('time', 1)               # compute linear regression coefficients\n",
    "linfit.polyfit_coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope = linfit.polyfit_coefficients.isel(degree=0)   # Polynomial coefficients, highest power first\n",
    "ns_century = 1.e9*3600.*24.*365.*100.                # scaling: nono-seconds to century\n",
    "(slope*ns_century).where(mask == 1.).plot(vmin=-5,vmax=5,cmap='RdBu_r')\n",
    "plt.title('Trend Slope in deg/century')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Detrend SST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,2, figsize=(6,3), constrained_layout=True)\n",
    "a = data.sst.groupby('time.year').mean(dim='time')      # yearly mean\n",
    "a.weighted(weights).mean(dim=('lon','lat')).plot(ax=axes[0], label=\"org\")\n",
    "a = trend.groupby('time.year').mean(dim='time')         # yearly mean\n",
    "a.weighted(weights).mean(dim=('lon','lat')).plot(ax=axes[0], label=\"trend\")\n",
    "axes[0].legend()\n",
    "axes[0].set_title('With trend')\n",
    "#\n",
    "a = sst_detrend.groupby('time.year').mean(dim='time')   # yearly mean\n",
    "a.weighted(weights).mean(dim=('lon','lat')).plot(ax=axes[1], label=\"detrend\")\n",
    "axes[1].legend()\n",
    "axes[1].set_title('Without trend')"
   ]
  },
  {
   "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": [
    "Interannual anomaly (with/without trend)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = data.sst.weighted(weights).mean(dim=('lon','lat'))\n",
    "clim = a.groupby(\"time.month\").mean(\"time\")   # climatological months (seasonal cycle)\n",
    "anom = a.groupby(\"time.month\") - clim         # interannual annomaly\n",
    "anom.plot(label=\"org\")\n",
    "#\n",
    "sstanom.weighted(weights).mean(dim=('lon','lat')).plot(label=\"detrended\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. SST Regression onto Nino3.4:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interannual standard deviation\n",
    "with [nino3.4](https://www.ncdc.noaa.gov/teleconnections/enso/indicators/sst/) box (5S-5N and 170-120W)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,2, figsize=(10,5), constrained_layout=True)\n",
    "sst_detrend.std(dim='time').where(mask == 1.).plot(ax=axes[0],vmin=0,vmax=5,cmap='YlGnBu')\n",
    "axes[0].set_title('SST detrend STD')\n",
    "sstanom.std(dim='time').where(mask == 1.).plot(ax=axes[1],vmin=0,vmax=2,cmap='YlGnBu')\n",
    "axes[1].set_title('SST detrend + anom STD')\n",
    "# overplot nino3.4 rectangle (170W-120W, 5S-5N)\n",
    "ax = plt.gca()                                                              # Get the current reference\n",
    "rect = Rectangle((360-170,-5),50,10,linewidth=1,edgecolor='r',fill=False)   # Create a Rectangle patch\n",
    "ax.add_patch(rect)                                                          # Add the patch to the Axes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Equatorial hovmoller of the SST interannual anomaly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,2,figsize=(10, 10))\n",
    "wgteq = weights.sel(lat=slice(2.5,-2.5), lon=slice(100.5, 285.5))\n",
    "#\n",
    "ssteq = sst_detrend.sel(lat=slice(2.5,-2.5), lon=slice(100.5, 285.5))\n",
    "ssteq.weighted(wgteq).mean(dim='lat').plot(ax=axes[0])\n",
    "#\n",
    "ssteq = sstanom.sel(lat=slice(2.5,-2.5), lon=slice(100.5, 285.5))\n",
    "ssteq.weighted(wgteq).mean(dim='lat').plot(ax=axes[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute nino3.4 (5S-5N and 170-120W) index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights_nino34 = weights.sel(lon=slice(360-169.5,360-120.5), lat=slice(4.5,-4.5))\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": [
    "### Monthly std of nino3.4 index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nino34_index.groupby('time.month').std(dim='time').plot.step(where='mid')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nino/Nina dissymmetry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(nino34_index.where(nino34_index > 0.).std().data)\n",
    "print(nino34_index.where(nino34_index < 0.).std().data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SST anomaly maps for some Nino Events (DJF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nino34_index.where(nino34_index > 1.5, drop=True).time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(7,1, figsize=(6, 16), constrained_layout=True)\n",
    "dataarray = sstanom.sel(time=slice('1982-12','1983-02')).mean(dim='time').where(mask == 1.)\n",
    "dataarray.plot(ax=axes[0], vmin=-2, vmax=2, cmap='RdBu_r')\n",
    "dataarray = sstanom.sel(time=slice('1987-12','1988-02')).mean(dim='time').where(mask == 1.)\n",
    "dataarray.plot(ax=axes[1], vmin=-2, vmax=2, cmap='RdBu_r')\n",
    "dataarray = sstanom.sel(time=slice('1991-12','1992-02')).mean(dim='time').where(mask == 1.)\n",
    "dataarray.plot(ax=axes[2], vmin=-2, vmax=2, cmap='RdBu_r')\n",
    "dataarray = sstanom.sel(time=slice('1997-12','1998-02')).mean(dim='time').where(mask == 1.)\n",
    "dataarray.plot(ax=axes[3], vmin=-2, vmax=2, cmap='RdBu_r')\n",
    "dataarray = sstanom.sel(time=slice('2002-12','2003-02')).mean(dim='time').where(mask == 1.)\n",
    "dataarray.plot(ax=axes[4], vmin=-2, vmax=2, cmap='RdBu_r')\n",
    "dataarray = sstanom.sel(time=slice('2009-12','2010-02')).mean(dim='time').where(mask == 1.)\n",
    "dataarray.plot(ax=axes[5], vmin=-2, vmax=2, cmap='RdBu_r')\n",
    "dataarray = sstanom.sel(time=slice('2015-12','2016-02')).mean(dim='time').where(mask == 1.)\n",
    "dataarray.plot(ax=axes[6], vmin=-2, vmax=2, cmap='RdBu_r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SST anomalies regressed onto Nino3.4 SST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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')"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}