Use remote data; revise text

environment.yml

@@ -18,6 +18,7 @@ dependencies:
   - nc-time-axis
   - gcsfs
   - cftime
+  - pydap
   - pip
   - pip:
       - sphinx-pythia-theme

notebooks/climate-modes-xeofs.ipynb

@@ -72,14 +72,14 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Preparing the data"
+    "## Accessing and preparing the data"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We will use the monthly [Extended Reconstructed Sea Surface Temperature](https://www.ncei.noaa.gov/products/extended-reconstructed-sst) (ERSST) version 3b dataset from NOAA, which is included locally in this cookbook."
+    "We will use the [NOAA Extended Reconstructed Sea Surface Temperature version 5 (ERSSTv5)](https://www.psl.noaa.gov/data/gridded/data.noaa.oisst.v2.html) monthly gridded dataset, which is accessible using [OPeNDAP](https://www.opendap.org/). More information on [using OPeNDAP to access NOAA data can be found here](https://psl.noaa.gov/data/help/using_opendap.html)."
    ]
   },
   {
@@ -88,15 +88,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sst = xr.open_dataset('./data/sst.mnmean.nc')\n",
-    "sst"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Define a helpful function to take global averages:"
+    "data_url = 'https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.ersst.v5/sst.mnmean.nc'"
    ]
   },
   {
@@ -105,17 +97,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def global_average(data):\n",
-    "    weights = np.cos(np.deg2rad(data.lat))\n",
-    "    data_weighted = data.weighted(weights)\n",
-    "    return data_weighted.mean(dim=['lat', 'lon'], skipna=True)"
+    "sst = xr.open_dataset(data_url, engine='pydap').sst\n",
+    "sst"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Here we compute the SST anomaly by subtracting out the average of each month using Xarray's `.groupby()` method. This deseasonalizes the data, which is important because climate modes will get buried under the magnitude of the seasonal changes in SST."
+    "Check that the data looks as expected:"
    ]
   },
   {
@@ -124,15 +114,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sst_clim = sst.groupby('time.month')\n",
-    "ssta = sst_clim - sst_clim.mean(dim='time')"
+    "sst.isel(time=0).plot()"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "To see what happens when we skip this step, let's do an EOF analysis on the unmodified monthly data:"
+    "Before we modify the data, let's do an EOF analysis on the whole dataset:"
    ]
   },
   {
@@ -141,7 +130,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "s_model = EOF(sst.sst, n_modes=4, dim=['time'], weights='coslat')\n",
+    "s_model = EOF(sst, n_modes=4, dim=['time'], weights='coslat')\n",
     "s_model.solve()\n",
     "s_eofs = s_model.eofs()\n",
     "s_pcs = s_model.pcs()\n",
@@ -179,7 +168,52 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "For EOF1, which explains about 86% of the variance, we can see there is interhemispheric asymmetry, and PC1 has a period of one year. This is the seasonal cycle. EOF3 and EOF4 are clearly picking up the global warming signal with some other variability mixed in, and EOF2 seems to be some combination of both. Here we have identified something else we will want to remove from the data if we are currently interested in looking at internal variability—the long-term warming. We can detrend the data by removing the global average SST anomaly:"
+    "EOF1 explains 83% of the variance, and the map shows interhemispheric asymmetry. The corresponding PC has a period of one year, which we can see more clearly by only plotting a few years:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "s_pcs.sel(mode=1, time=slice('1900', '1903')).plot(figsize=(8, 3))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This mode is showing the seasonal cycle. This is interesting, but it obfuscates other modes. If we want to study the other ways Earth's climate varies, we should remove the seasonal cycle from our data. Here we compute this (calling it the SST anomaly) by subtracting out the average of each month using Xarray's `.groupby()` method:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sst_clim = sst.groupby('time.month')\n",
+    "ssta = sst_clim - sst_clim.mean(dim='time')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The remaining 3 EOFs show a combination of the long-term warming trend, the seasonal cycle (EOF analyses do not cleanly separate physical modes), and other internal variability. The warming trend is also interesting (see the [CMIP6 Cookbook](https://projectpythia.org/cmip6-cookbook)), but here we want to pull out some modes of internal/natural variability. We can detrend the data by removing the global average SST anomaly."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def global_average(data):\n",
+    "    weights = np.cos(np.deg2rad(data.lat))\n",
+    "    data_weighted = data.weighted(weights)\n",
+    "    return data_weighted.mean(dim=['lat', 'lon'], skipna=True)"
    ]
   },
   {
@@ -204,7 +238,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds_model = EOF(ssta_dt.sst, n_modes=4, dim=['time'], weights='coslat')\n",
+    "ds_model = EOF(ssta_dt, n_modes=4, dim=['time'], weights='coslat')\n",
     "ds_model.solve()\n",
     "ds_eofs = ds_model.eofs()\n",
     "ds_pcs = ds_model.pcs()\n",
@@ -242,7 +276,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Now we can see some modes of variability! EOF1 looks like ENSO or IPO, and EOF2 looks like it is picking up a pattern of recent warming where the Southern Ocean and southeastern Pacific are slightly cooling. EOF3 shows PDO in the North Pacific. There is a lot going on in each of these maps, so to get a clearer index of some of these modes, we can restrict our domain. "
+    "Now we can see some modes of variability! EOF1 looks like ENSO or IPO, and EOF2 probably picking up a pattern of the recent temperature trend where the Southern Ocean and southeastern Pacific are slightly cooling. EOF3 and EOF4 appear to be showing some decadal modes of variability (PDO and maybe AMO), among other things. There is a lot going on in each of these maps, so to get a clearer index of some modes, we can restrict our domain. "
    ]
   },
   {
@@ -251,7 +285,7 @@
    "source": [
     "## El Niño Southern Oscillation (ENSO)\n",
     "\n",
-    "Here we restrict our domain to the equatorial Pacific. You can [read more about ENSO here](https://www.ncei.noaa.gov/access/monitoring/enso/)."
+    "Here we restrict our domain to the equatorial Pacific. Note that ENSO is commonly defined using an index of SST anomaly over a region of the equatorial Pacific (e.g., the [Oceanic Niño Index (ONI)](https://www.ncei.noaa.gov/access/monitoring/enso/sst)) instead of an EOF. You can [read more about ENSO here](https://www.ncei.noaa.gov/access/monitoring/enso/)."
    ]
   },
   {
@@ -269,7 +303,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ep_model = EOF(ep_ssta_dt.sst, n_modes=4, dim=['time'], norm=True, weights='coslat')\n",
+    "ep_model = EOF(ep_ssta_dt, n_modes=4, dim=['time'], norm=True, weights='coslat')\n",
     "ep_model.solve()\n",
     "ep_eofs = ep_model.eofs()\n",
     "ep_pcs = ep_model.pcs()\n",
@@ -308,14 +342,16 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "fig, ax = plt.subplots(1, 1, figsize=(10, 2), dpi=130)\n",
     "plt.fill_between(ep_pcs.time, ep_pcs.isel(mode=0).where(ep_pcs.isel(mode=0) > 0), color='r')\n",
     "plt.fill_between(ep_pcs.time, ep_pcs.isel(mode=0).where(ep_pcs.isel(mode=0) < 0), color='b')\n",
     "plt.ylabel('PC')\n",
-    "plt.xlabel('Time')\n",
+    "plt.xlabel('Year')\n",
     "plt.xlim(ep_pcs.time.min(), ep_pcs.time.max())\n",
     "plt.grid(linestyle=':')\n",
     "plt.title('ENSO Index (detrended equatorial Pacific SSTA EOF1)')"
@@ -325,7 +361,39 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Pacific Decadal Oscillation (PDO)"
+    "Compare to the ONI:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(1, 1, figsize=(10, 2), dpi=130)\n",
+    "plt.fill_between(ep_pcs.time, ep_pcs.isel(mode=0).where(ep_pcs.isel(mode=0) > 0), color='r')\n",
+    "plt.fill_between(ep_pcs.time, ep_pcs.isel(mode=0).where(ep_pcs.isel(mode=0) < 0), color='b')\n",
+    "plt.ylabel('PC')\n",
+    "plt.xlabel('Year')\n",
+    "plt.xlim(ep_pcs.time.sel(time='1950-01').squeeze(), ep_pcs.time.max())\n",
+    "plt.grid(linestyle=':')\n",
+    "plt.title('ENSO Index (detrended equatorial Pacific SSTA EOF1)')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<img src=\"images/oni.png\" alt=\"ONI\"></img>"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Pacific Decadal Oscillation (PDO)\n",
+    "\n",
+    "Here we restrict our domain to the North Pacific. You can [read more about PDO here](https://www.ncei.noaa.gov/access/monitoring/pdo/)."
    ]
   },
   {
@@ -343,7 +411,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "np_model = EOF(np_ssta_dt.sst, n_modes=4, dim=['time'], norm=True, weights='coslat')\n",
+    "np_model = EOF(np_ssta_dt, n_modes=4, dim=['time'], norm=True, weights='coslat')\n",
     "np_model.solve()\n",
     "np_eofs = np_model.eofs()\n",
     "np_pcs = np_model.pcs()\n",
@@ -388,11 +456,11 @@
    "outputs": [],
    "source": [
     "fig, ax = plt.subplots(1, 1, figsize=(10, 2), dpi=130)\n",
-    "plt.fill_between(np_pcs.time, -np_pcs.isel(mode=0).where(np_pcs.isel(mode=0) < 0), color='r')\n",
-    "plt.fill_between(np_pcs.time, -np_pcs.isel(mode=0).where(np_pcs.isel(mode=0) > 0), color='b')\n",
-    "plt.plot(np_pcs.time, -np_pcs.isel(mode=0).rolling(time=48, center=True).mean(), color='k', linewidth=2)\n",
+    "plt.fill_between(np_pcs.time, np_pcs.isel(mode=0).where(np_pcs.isel(mode=0) > 0), color='r')\n",
+    "plt.fill_between(np_pcs.time, np_pcs.isel(mode=0).where(np_pcs.isel(mode=0) < 0), color='b')\n",
+    "plt.plot(np_pcs.time, np_pcs.isel(mode=0).rolling(time=48, center=True).mean(), color='k', linewidth=2)\n",
     "plt.ylabel('PC')\n",
-    "plt.xlabel('Time')\n",
+    "plt.xlabel('Year')\n",
     "plt.xlim(np_pcs.time.min(), np_pcs.time.max())\n",
     "plt.grid(linestyle=':')\n",
     "plt.title('PDO Index (detrended North Pacific SSTA EOF1)')"
@@ -410,25 +478,18 @@
    "metadata": {},
    "source": [
     "## Summary\n",
-    "Add one final `---` marking the end of your body of content, and then conclude with a brief single paragraph summarizing at a high level the key pieces that were learned and how they tied to your objectives. Look to reiterate what the most important takeaways were.\n",
+    "In this notebook, we demonstrated a basic workflow for performing an EOF analysis on gridded SST data using the `xeofs` package. We plotted the PCs associated with ENSO and PDO using deseasonalized, detrended SSTs.\n",
     "\n",
     "### What's next?\n",
-    "Let Jupyter book tie this to the next (sequential) piece of content that people could move on to down below and in the sidebar. However, if this page uniquely enables your reader to tackle other nonsequential concepts throughout this book, or even external content, link to it here!"
+    "The next section will focus on applications of EOF analysis to answer scientific questions. (Coming soon!)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "## Resources and references\n",
-    "Finally, be rigorous in your citations and references as necessary. Give credit where credit is due. Also, feel free to link to relevant external material, further reading, documentation, etc. Then you're done! Give yourself a quick review, a high five, and send us a pull request. A few final notes:\n",
-    " - `Kernel > Restart Kernel and Run All Cells...` to confirm that your notebook will cleanly run from start to finish\n",
-    " - `Kernel > Restart Kernel and Clear All Outputs...` before committing your notebook, our machines will do the heavy lifting\n",
-    " - Take credit! Provide author contact information if you'd like; if so, consider adding information here at the bottom of your notebook\n",
-    " - Give credit! Attribute appropriate authorship for referenced code, information, images, etc.\n",
-    " - Only include what you're legally allowed: **no copyright infringement or plagiarism**\n",
-    " \n",
-    "Thank you for your contribution!"
+    "Huang, B., Thorne, P. S., Banzon, V., Boyer, T. P., Chepurin, G. A., Lawrimore, J. H., Menne, M. J., Smith, T. J., Vose, R. S., & Zhang, H. (2017). Extended Reconstructed Sea Surface Temperature, Version 5 (ERSSTv5): Upgrades, Validations, and Intercomparisons. *Journal of Climate*, *30*(20), 8179–8205. https://doi.org/10.1175/jcli-d-16-0836.1"
    ]
   }
  ],

notebooks/eof-intro.ipynb

@@ -161,7 +161,7 @@
     "In this notebook, we covered some background on EOFs, the steps required to carry out an EOF analysis, and some of the math behind the analysis.\n",
     "\n",
     "### What's next?\n",
-    "An example of carrying out an EOF analysis on idealized data using NumPy."
+    "An example of carrying out an EOF analysis on idealized data using NumPy. (Coming soon!)"
    ]
   },
   {