Notebook3_EOFs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook #3: EOF analysis of LMR output\n",
    "\n",
    "This notebook walks you through simple steps to generate EOFs of LMR fields. It makes use of the eofs package:\n",
    " \n",
    "    Dawson, A., (2016). eofs: A Library for EOF Analysis of Meteorological, \n",
    "    Oceanographic, and Climate Data. Journal of Open Research Software. 4(1), \n",
    "    p.e14. DOI: http://doi.org/10.5334/jors.122\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import cartopy.util as cutil\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "from eofs.xarray import Eof\n",
    "import matplotlib.gridspec as gridspec\n",
    "# for graphics\n",
    "plt.style.use('ggplot')\n",
    "\n",
    "# Import a custom python function\n",
    "import sys\n",
    "sys.path.append('../python_functions')\n",
    "import filters as flt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LOAD DATA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the name of an experiment and open the data set using xarray.\n",
    "experiment_name = 'hackathon_PAGES2kv1'\n",
    "\n",
    "ds = xr.open_dataset('../data/'+experiment_name+'/tas_sfc_Amon_MCiters_ensemble_mean.nc',decode_times=False)\n",
    "print(ds)\n",
    "ds_ensMean = ds.mean(dim=('member')) # average over records\n",
    "\n",
    "tas     = ds_ensMean['tas']\n",
    "nt , ny, nx = tas.shape\n",
    "\n",
    "# glean coordinates\n",
    "lat   = tas['lat'].values\n",
    "lon   = tas['lon'].values\n",
    "years = tas['time'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EOF analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now let's perform EOF analysis per se. First, we need to define some weights proportional to the area of each grid cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#  Square-root of cosine of latitude weights are applied before the computation of EOFs.\n",
    "coslat = np.cos(np.deg2rad(lat))\n",
    "wgts = np.sqrt(np.abs(coslat[..., np.newaxis])) # ensure values near 0 are >0, and add extra axis for compatibility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an EOF solver to do the EOF analysis. the eofs package works natively with xarray, so there's very little we need to do. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tsolver = Eof(tas, weights=wgts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first let's inspect the eigenvalue spectrum\n",
    "pctvar = np.round(tsolver.varianceFraction(neigs=20).values*100) # eigenvalue spectrum\n",
    "f,ax = plt.subplots(2,1,sharex=True)\n",
    "ax[0].plot(pctvar,'o')\n",
    "ax[0].set_title('Eigenvalue spectrum',fontweight='bold')\n",
    "ax[0].set_ylabel('% variance')\n",
    "ax[1].plot(np.cumsum(pctvar),'o')\n",
    "ax[1].set_xlabel('Mode index')\n",
    "ax[1].set_ylabel('% total')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " You can see that the first mode accounts for ~50% of the variance, and the first 5, 90%. This is a standard \"red\" spectrum, as you'd expect for a geophysical field. \n",
    " \n",
    " Now let's retrieve the EOFs, expressed as the correlation between the leading PC time series and the input tas anomalies at each grid point, and the PC time series themselves. Let's keep 4 of those, since returns and interpretability diminish fast in EOF land."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eof = tsolver.eofsAsCovariance(neofs=4)\n",
    "pc  = tsolver.pcs(npcs=4, pcscaling=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Now let's plot the EOFs and PCs themselves. Because we're going to have to do this quite a bit, let's be smart and write a couple of functions to do that. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def contourLevels(DataArray,nc = 21):\n",
    "    '''\n",
    "    Computes evenly-spaced contour levels\n",
    "\n",
    "    '''\n",
    "    Li = max(abs(DataArray.max()),abs(DataArray.min()))\n",
    "    clevs = np.linspace(-Li/10, Li/10, nc)\n",
    "    return clevs, Li\n",
    "\n",
    "def eofPlotter_xr(eof,pc,pctvar,lon,proj,ttl,units,modes=(1,2),nc=31,cmap=plt.cm.RdBu_r,fnt='Open Sans'):\n",
    "    '''\n",
    "    Plot 2 EOF modes using the output of the xarray eof solver \n",
    "    '''\n",
    "    # declare figure\n",
    "    fig = plt.figure(1,figsize=(12,6))\n",
    "    gs = gridspec.GridSpec(2, 2, bottom=.05, top=.90, wspace=.1,width_ratios=[1, 1.3])\n",
    "    # plot EOF1\n",
    "    EOF1 = eof.sel(mode=modes[0]-1).values\n",
    "    EOF1c, lonc = cutil.add_cyclic_point(EOF1, lon) # make cyclic \n",
    "    ax1 = plt.subplot(gs[0, 0],projection=proj)\n",
    "    clevs, Li = contourLevels(EOF1)\n",
    "    #cf = ax1.contourf(lon,lat,-EOF1.values,levels=clevs, cmap=plt.cm.RdBu_r,transform=ccrs.PlateCarree())\n",
    "    cf = ax1.contourf(lonc,lat,EOF1c,nc, vmin=-Li,vmax=Li, cmap=cmap,transform=ccrs.PlateCarree())\n",
    "    ax1.set_global()\n",
    "    ax1.add_feature(cfeature.COASTLINE, edgecolor='k')\n",
    "    ax1.set_title(r'EOF '+str(modes[0]),fontname=fnt,fontsize=12)\n",
    "    #ax1.set_title(r'EOF '+str(modes[0])+' ('+units+', '+ str(pctvar[modes[0]-1])+ '% var)',fontname=fnt,fontsize=12)\n",
    "    cb = plt.colorbar(cf,ax=ax1, orientation='vertical')\n",
    "    # now plot PC1\n",
    "    pc1 = pc.sel(mode=modes[0]-1)\n",
    "    ax2 = plt.subplot(gs[0, 1])\n",
    "    ax2.plot(years,pc1); ax2.set_title('PC '+str(modes[0]))\n",
    "    plt.setp(ax2.get_xticklabels(), visible=False)\n",
    "\n",
    "    # plot EOF2\n",
    "    EOF2 = eof.sel(mode=modes[1]-1).values\n",
    "    EOF2c, lonc = cutil.add_cyclic_point(EOF2, lon) # make cyclic \n",
    "    clevs, Li = contourLevels(EOF2)\n",
    "    ax3 = plt.subplot(gs[1, 0],projection=proj)\n",
    "    cf = ax3.contourf(lonc,lat,EOF2c,nc, vmin=-Li,vmax=Li, cmap=cmap,transform=ccrs.PlateCarree())\n",
    "    ax3.set_global()\n",
    "    ax3.add_feature(cfeature.COASTLINE, edgecolor='k')\n",
    "    ax3.set_title(r'EOF '+str(modes[1]),fontname=fnt,fontsize=12)\n",
    "    #ax3.set_title(r'EOF '+str(modes[1])+' ('+units+', '+ str(pctvar[modes[1]-1])+ '% var)',fontname=fnt,fontsize=12)\n",
    "    cb = plt.colorbar(cf,ax=ax3, orientation='vertical')\n",
    "    # now plot PC2\n",
    "    pc2 = pc.sel(mode=modes[1]-1)\n",
    "    ax4 = plt.subplot(gs[1, 1])\n",
    "    ax4.plot(years,pc2); ax4.set_title('PC '+str(modes[1])); ax4.set_xlabel('Year CE')\n",
    "    \n",
    "    fig.suptitle(ttl,fontweight='bold',fontsize=16,fontname=fnt)    \n",
    "    ax = (ax1,ax2,ax3,ax3)\n",
    "    fn = ttl + '_modes'+str(modes)\n",
    "    fstr = fn.replace(' ','_').replace(',',\"\").replace('(',\"\").replace(')',\"\") # make filename from title\n",
    "\n",
    "    return fig, ax, fstr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot unfiltered EOFS\n",
    "proj = ccrs.Orthographic(central_longitude=0.0, central_latitude=90)  \n",
    "units = '\\xb0C'\n",
    "ttl = 'LMR unfiltered'\n",
    "  \n",
    "f1, ax1, fn = eofPlotter_xr(eof,pc,pctvar,lon,proj,ttl,units,modes=(1,2)) # invoke plotter\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ain't it pretty? \n",
    "The first mode clearly corresponds to a \"global warming\" mode with a strong 20th century expression and a spatial pattern suggestive of polar amplification. Mode 2 has a maximum over Alaska and a cold annulus throughout much of the northern midlatitudes.  \n",
    "\n",
    "Now let's look at modes 3 and 4:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f2, ax2, fn = eofPlotter_xr(eof,pc,pctvar,lon,proj,ttl,units,modes=(3,4))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By definition, all these modes are bi-orthogonal, and they get harder to interpret as we go down the list.\n",
    "\n",
    "Also, annual data are rather noisy. Let us zoom in on multidecadal variations (say, 20y and longer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Applying a low-pass filter\n",
    "\n",
    "Here we'll be using the filters package provided under \"python_functions\". We'll be using a Butterworth filter with default options. to see what they are, type: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "flt.butterworth?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# apply 20y low-pass filter\n",
    "Tc= 20; fc = 1./Tc\n",
    "ns = ny*nx  # number of spatial points\n",
    "tas_v = np.reshape(tas.values,(nt, ns)) # reshape into 2d\n",
    "tas_f = np.empty(tas_v.shape) # create recipient array\n",
    "\n",
    "for j in range(ns):  # loop over all grid points\n",
    "    tas_f[:,j] = flt.butterworth(tas_v[:,j],fc,pad='reflect') # apply Butterworth filter with default options\n",
    "\n",
    "tas_f = np.reshape(tas_f,(nt , ny, nx)) # reshape into 3d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the filer works on NumPy arrays, but doesn't \"know\" about xarray. We can apply EOF analysis on this new (NumPy) aray, but we need to invoke the eofs package in a slightly different way. We also need to modify our plotter routine so it handles NumPy arrays, not just xarrays. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run EOF solver\n",
    "from eofs.standard import Eof as SEof\n",
    "tsolver = SEof(tas_f, weights=wgts)\n",
    "pctvar10y = np.round(tsolver.varianceFraction(neigs=20)*100)\n",
    "# first let's  inspect the eigenvalue spectrum\n",
    "f,ax = plt.subplots(2,1,sharex=True)\n",
    "ax[0].plot(pctvar10y,'o')\n",
    "ax[0].set_title('Eigenvalue spectrum',fontweight='bold')\n",
    "ax[0].set_ylabel('% variance')\n",
    "ax[1].plot(np.cumsum(pctvar10y),'o')\n",
    "ax[1].set_xlabel('Mode index')\n",
    "ax[1].set_ylabel('% total')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is similar to before, except that the first mode now explains nearly 60% of the variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "eof10y = tsolver.eofsAsCovariance(neofs=2)\n",
    "pc10y  = tsolver.pcs(npcs=2, pcscaling=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def eofPlotter_np(eof,pc,pctvar,lon,proj,ttl,units,modes=(1,2),nc=31,cmap=plt.cm.RdBu_r,fnt='Open Sans'):\n",
    "    '''\n",
    "    Plot 2 EOF modes using the output of the 'standard' eof solver\n",
    "    '''\n",
    "    # declare figure\n",
    "    fig = plt.figure(1,figsize=(12,6))\n",
    "    gs = gridspec.GridSpec(2, 2, bottom=.05, top=.90, wspace=.1,width_ratios=[1, 1.3])\n",
    "    # plot EOF1\n",
    "    EOF1 = eof[modes[0]-1,...]\n",
    "    EOF1c, lonc = cutil.add_cyclic_point(EOF1, lon) # make cyclic\n",
    "    ax1 = plt.subplot(gs[0, 0],projection=proj)\n",
    "    clevs, Li = contourLevels(EOF1)\n",
    "    #cf = ax1.contourf(lon,lat,-EOF1.values,levels=clevs, cmap=plt.cm.RdBu_r,transform=ccrs.PlateCarree())\n",
    "    cf = ax1.contourf(lonc,lat,EOF1c,nc, vmin=-Li,vmax=Li, cmap=cmap,transform=ccrs.PlateCarree())\n",
    "    ax1.set_global()\n",
    "    ax1.add_feature(cfeature.COASTLINE, edgecolor='k')\n",
    "    ax1.set_title(r'EOF '+str(modes[0]),fontname=fnt,fontsize=12)\n",
    "    #ax1.set_title(r'EOF '+str(modes[0])+' ('+units+', '+ str(pctvar[modes[0]-1])+ '% var)',fontname=fnt,fontsize=12)\n",
    "    cb = plt.colorbar(cf,ax=ax1, orientation='vertical')\n",
    "    # now plot PC1\n",
    "    pc1 = pc[:,modes[0]-1]\n",
    "    ax2 = plt.subplot(gs[0, 1])\n",
    "    ax2.plot(years,pc1); ax2.set_title('PC '+str(modes[0]))\n",
    "    plt.setp(ax2.get_xticklabels(), visible=False)\n",
    "\n",
    "    # plot EOF2\n",
    "    EOF2 = eof[modes[1]-1,...]\n",
    "    EOF2c, lonc = cutil.add_cyclic_point(EOF2, lon) # make cyclic\n",
    "    clevs, Li = contourLevels(EOF2)\n",
    "    ax3 = plt.subplot(gs[1, 0],projection=proj)\n",
    "    cf = ax3.contourf(lonc,lat,EOF2c,nc, vmin=-Li,vmax=Li, cmap=cmap,transform=ccrs.PlateCarree())\n",
    "    ax3.set_global()\n",
    "    ax3.add_feature(cfeature.COASTLINE, edgecolor='k')\n",
    "    ax3.set_title(r'EOF '+str(modes[1]),fontname=fnt,fontsize=12)\n",
    "    #ax3.set_title(r'EOF '+str(modes[1])+' ('+units+', '+ str(pctvar[modes[1]-1])+ '% var)',fontname=fnt,fontsize=12)\n",
    "    cb = plt.colorbar(cf,ax=ax3, orientation='vertical')\n",
    "    # now plot PC2\n",
    "    pc2 = pc[:,modes[1]-1]\n",
    "    ax4 = plt.subplot(gs[1, 1])\n",
    "    ax4.plot(years,pc2); ax4.set_title('PC '+str(modes[1])); ax4.set_xlabel('Year CE')\n",
    "\n",
    "    fig.suptitle(ttl,fontweight='bold',fontsize=16,fontname=fnt)\n",
    "    ax = (ax1,ax2,ax3,ax3)\n",
    "    fn = ttl + '_modes'+str(modes)\n",
    "    fstr = fn.replace(' ','_').replace(',',\"\").replace('(',\"\").replace(')',\"\") # make filename from title\n",
    "\n",
    "    return fig, ax, fstr\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now plot\n",
    "ttl = 'LMR, 20y lowpass filter'\n",
    "f1, ax1, fn = eofPlotter_np(eof10y,pc10y,pctvar10y,lon,proj,ttl,units,modes=(1,2)) # invoke plotter\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The patterns are ostensibly the same as in the unfiltered case, but explain more variance overall, and the temporal signature is easier to see. Notice the practically flat curves over the first millennium BCE. Looking only at the ensemble mean, one would think that climate ground to a halt during this interval. This is, of course, be preposterous. The reality is that the ensemble mean is not a useful summary of what's going on. Let's therefore look at ensemble EOFs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ensemble EOFs\n",
    "\n",
    "Let's start by extracting the ensemble members and reshaping the matrix in a friendly form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_all = xr.open_dataset('../data/'+experiment_name+'/tas_sfc_Amon_MCiters_ensemble_full_all.nc',decode_times=False)\n",
    "print(ds_all)\n",
    "\n",
    "tas     = ds_all['tas']\n",
    "tas = np.reshape(tas.values, (11*20,2001, 42, 63))\n",
    "nens, nt, ny, nx = tas.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now loop over ensemble members and apply EOF analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npcs = 1  # keep only first PC\n",
    "nsample = 50\n",
    "ensPC = np.empty((nsample,nt,2)) # define receiver PC matrix\n",
    "select = np.random.randint(0,nens,nsample)\n",
    "for je in range(nsample):\n",
    "    print(\"Solving EOF problem for ensemble member \"+str(je+1)+\"/\"+ str(nsample))\n",
    "    EnsSolver = SEof(tas[select[je],...], weights=wgts)\n",
    "    eof = EnsSolver.eofsAsCovariance(neofs=1)\n",
    "    if (eof.max()>0.5):\n",
    "        sign = +1\n",
    "    else:\n",
    "        sign = -1 \n",
    "        \n",
    "    ensPC[je,...] = sign*EnsSolver.pcs(npcs=npcs, pcscaling=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot the result. Let's define a function to do that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def tsEnsPlot(t, X, ax, col='RoyalBlue',ns=7):\n",
    "    \"\"\"Plots results of an ensemble of timeseries X(t)\n",
    "       X: 2-dim Numpy Arrray\n",
    "       t: time axis (needs to have same dimesion as X # of rows or columns)\n",
    "       ax: axis to plot into\n",
    "       col: color of patch\n",
    "       ns: number of sample lines\n",
    "    \"\"\"\n",
    "    plt.style.use('ggplot')\n",
    "\n",
    "    nt = len(t)\n",
    "    tdim = X.shape.index(nt)\n",
    "    if tdim == 1:\n",
    "       X = X.T\n",
    "\n",
    "    nens = X.shape[1]\n",
    "    Xq = np.percentile(X,(2.5,25,50,75,97.5),axis=1).T # compute quantiles\n",
    "    # plot quantiles\n",
    "    ax.fill_between(t,Xq[:,0],Xq[:,4],color=col,alpha=0.3,label='95% HDR')\n",
    "    ax.fill_between(t,Xq[:,1],Xq[:,3],color=col,alpha=0.6,label='IQR')\n",
    "    if ns>0:\n",
    "        # plot individual members\n",
    "        ax.plot(t,X[:,np.random.randint(nens,size=ns)],color='Gray',linewidth=0.5,alpha=0.4)\n",
    "    # in any case, plot median\n",
    "    ax.plot(t,Xq[:,2],color='Black',linewidth=1,label='median')\n",
    "    ax.legend(loc='upper left',ncol=4)\n",
    "    return ax, Xq\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f,ax = plt.subplots(1,1)\n",
    "_ = tsEnsPlot(years, ensPC[:,:,0].T, ax=ax, col='RoyalBlue',ns=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we see that the first PC is quite consistent from iteration to iteration. Thus, using the ensemble mean for the EOF plots above wasn't as misguided as it might seem *a priori*. Breathe a huge sigh of relief, and please move on to the next activity. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}