M04N05_using_the_eofs_package.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ATMS 523\n",
    "## using the `eofs` python package\n",
    "\n",
    "For performing EOF analysis in python, there are several examples.  We'll use the `eofs` package by Andrew Dawson.  It is available on GitHub here: https://github.com/ajdawson/eofs\n",
    "\n",
    "I install it in my python environment using the instructions in the GitHub repository readme.\n",
    "\n",
    "`eofs` works directly with `xarray` and allows you to preserve metadata, which makes it easy to work with.  Let's go through an example using Pacific Sea Surface Temperatures that will help us identify the pattern of the El Nino Southern Oscillation in monthly climate data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "\n",
    "from eofs.xarray import Eof\n",
    "from eofs.examples import example_data_path\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read SST anomalies using the xarray module. The file contains November-March\n",
    "# averages of SST anomaly in the central and northern Pacific.\n",
    "filename = example_data_path('sst_ndjfm_anom.nc')\n",
    "sst = xr.open_dataset(filename)['sst']\n",
    "sst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an EOF solver to do the EOF analysis. Square-root of cosine of\n",
    "# latitude weights are applied before the computation of EOFs.\n",
    "coslat = np.cos(np.deg2rad(sst.coords['latitude'].values))\n",
    "wgts = np.sqrt(coslat)[..., np.newaxis]\n",
    "solver = Eof(sst, weights=wgts)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve the leading EOF, expressed as the correlation between the leading\n",
    "# PC time series and the input SST anomalies at each grid point, and the\n",
    "# leading PC time series itself.\n",
    "eof1 = solver.eofsAsCorrelation(neofs=2)\n",
    "pc1 = solver.pcs(npcs=2, pcscaling=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the leading EOF expressed as correlation in the Pacific domain.\n",
    "clevs = np.linspace(-1, 1, 11)\n",
    "ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=190))\n",
    "fill = eof1[0].plot.contourf(ax=ax, levels=clevs, cmap=plt.cm.RdBu_r,\n",
    "                             add_colorbar=False, transform=ccrs.PlateCarree())\n",
    "ax.add_feature(cfeature.COASTLINE, color='k', edgecolor='k')\n",
    "cb = plt.colorbar(fill, orientation='horizontal')\n",
    "cb.set_label('correlation coefficient', fontsize=12)\n",
    "ax.set_title('EOF1 expressed as correlation', fontsize=16)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the leading PC time series.\n",
    "plt.figure()\n",
    "pc1[:, 0].plot(color='b', linewidth=2)\n",
    "ax = plt.gca()\n",
    "ax.axhline(0, color='k')\n",
    "ax.set_ylim(-3, 3)\n",
    "ax.set_xlabel('Year')\n",
    "ax.set_ylabel('Normalized Units')\n",
    "ax.set_title('PC1 Time Series', fontsize=16)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the leading EOF expressed as correlation in the Pacific domain.\n",
    "clevs = np.linspace(-1, 1, 11)\n",
    "ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=190))\n",
    "fill = eof1[1].plot.contourf(ax=ax, levels=clevs, cmap=plt.cm.RdBu_r,\n",
    "                             add_colorbar=False, transform=ccrs.PlateCarree())\n",
    "ax.add_feature(cfeature.COASTLINE, color='k', edgecolor='k')\n",
    "cb = plt.colorbar(fill, orientation='horizontal')\n",
    "cb.set_label('correlation coefficient', fontsize=12)\n",
    "ax.set_title('EOF1 expressed as correlation', fontsize=16)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the second PC time series.\n",
    "plt.figure()\n",
    "pc1[:, 1].plot(color='b', linewidth=2)\n",
    "ax = plt.gca()\n",
    "ax.axhline(0, color='k')\n",
    "ax.set_ylim(-3, 3)\n",
    "ax.set_xlabel('Year')\n",
    "ax.set_ylabel('Normalized Units')\n",
    "ax.set_title('PC2 Time Series', fontsize=16)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varfrac = solver.varianceFraction()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the fraction of variance explained by each EOF\n",
    "plt.figure(figsize=(11,6))\n",
    "eof_num = range(1, 16)\n",
    "plt.plot(eof_num, varfrac[0:15], linewidth=2)\n",
    "plt.plot(eof_num, varfrac[0:15], linestyle='None', marker=\"o\", color='r', markersize=8)\n",
    "plt.axhline(0, color='k')\n",
    "plt.xticks(range(1, 16))\n",
    "plt.title('Fraction of the total variance represented by each EOF')\n",
    "plt.xlabel('EOF #')\n",
    "plt.ylabel('Variance Fraction')\n",
    "plt.xlim(1, 15)\n",
    "plt.ylim(np.min(varfrac), np.max(varfrac)+0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reconstruction = solver.reconstructedField(4)\n",
    "reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the leading EOF expressed as correlation in the Pacific domain.\n",
    "clevs = np.linspace(-3, 3, 21)\n",
    "ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=190))\n",
    "fill = reconstruction.sel(time='1998-01-15', method='nearest').plot.contourf(ax=ax, levels=clevs, cmap=plt.cm.RdBu_r,\n",
    "                             add_colorbar=False, transform=ccrs.PlateCarree())\n",
    "ax.add_feature(cfeature.COASTLINE, color='k', edgecolor='k')\n",
    "cb = plt.colorbar(fill, orientation='horizontal')\n",
    "cb.set_label('SST Anomalies (reconstructed)', fontsize=12)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the leading EOF expressed as correlation in the Pacific domain.\n",
    "clevs = np.linspace(-3, 3, 21)\n",
    "ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=190))\n",
    "fill = sst.sel(time='1998-01-15', method='nearest').plot.contourf(ax=ax, levels=clevs, cmap=plt.cm.RdBu_r,\n",
    "                             add_colorbar=False, transform=ccrs.PlateCarree())\n",
    "ax.add_feature(cfeature.COASTLINE, color='k', edgecolor='k')\n",
    "cb = plt.colorbar(fill, orientation='horizontal')\n",
    "cb.set_label('SST Anomalies observed', fontsize=12)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.7 ('py3')",
   "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.9.7"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "7a6b41cc1ffbe7f6292ade58ee9ab0c89bd7fa770a77f101cc95d5710b4e5fa9"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}