xarray-tutorial-egu2017.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [SC57 - Working with big, multi-dimensional geoscientific datasets in Python: a tutorial introduction to xarray](http://meetingorganizer.copernicus.org/EGU2017/session/25651)  \n",
    "  \n",
    "  \n",
    "Original notebook by [Stephan Hoyer](http://stephanhoyer.com), Rossbypalooza, 2016.  \n",
    "Modified by Edward Byers, Matthew Gidden and [Fabien Maussion](http://fabienmaussion.info/) for EGU General Assembly 2017, Vienna, Austria\n",
    "\n",
    "\n",
    "  \n",
    "  Thursday, 27th April, 15:30–17:00 / Room -2.91  \n",
    "  \n",
    "  \n",
    "**Convenors**\n",
    "* [Dr Edward Byers](mailto:byers@iiasa.ac.at)    - International Institute for Applied Systems Analysis, Laxenburg, Austria\n",
    "* [Dr Matthew Gidden](mailto:gidden@iiasa.ac.at)  - International Institute for Applied Systems Analysis, Laxenburg, Austria\n",
    "* [Dr Fabien Maussion](mailto:fabien.maussion@uibk.ac.at) - University of Innsbruck, Innsbruck, Austria\n",
    "-------------\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# With\n",
    "\n",
    "![](./figures/dataset-diagram-logo.png)\n",
    "\n",
    "# you can reach\n",
    "\n",
    "![](./figures/facet-plot.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Structure of this tutorial\n",
    "\n",
    "1. Introduction to key features of `xarray`\n",
    "2. Basic operations in xarray: opening, inspecting, selecting and indexing data\n",
    "3. Selecting data with named dimensions\n",
    "4. Operations and computation\n",
    "5. Groupby and \"split-apply-combine\"\n",
    "6. Graphics\n",
    "7. Out-of-core computation\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Key features of `xarray`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is `xarray`?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*  `xarray` is an open source project and Python package\n",
    "*  `xarray` has been designed to perform **labelled** data analysis on **multi-dimensional** arrays\n",
    "* the xarray approach adopts the Common Data Model for **self-describing scientific data** in widespread use in the Earth sciences\n",
    "*  `xarray.Dataset` is an in-memory representation of a netCDF file.\n",
    "* `xarray` is built on top of the dataprocessing library [Pandas](http://pandas.pydata.org) (the best way to work with tabular data (e.g., CSV files) in Python)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Our data\n",
    "\n",
    "<img src=\"./figures/dataset.png\" width=\"50%\" align=\"right\"> \n",
    "\n",
    "- numeric\n",
    "- multi-dimensional\n",
    "- labelled\n",
    "- (lots of) metadata\n",
    "- sometimes (very) large"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is `xarray` good for?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Gridded, multi-dimensional and large datasets, commonly used in earth sciences, but also increasingly finance, engineering (signal/image processing), and biological sciences\n",
    "* Integration with other data analysis packages such as Pandas \n",
    "* I/O operations (NetCDF)\n",
    "* Plotting\n",
    "* Out of core computation and parallel processing\n",
    "* Extensions based on xarray\n",
    "* ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Where can I find more info?\n",
    "\n",
    "### For more information about xarray\n",
    "\n",
    "- Read the [online documentation](http://xarray.pydata.org/)\n",
    "- Ask questions on [StackOverflow](http://stackoverflow.com/questions/tagged/python-xarray)\n",
    "- View the source code and file bug reports on [GitHub](http://github.com/pydata/xarray/)\n",
    "\n",
    "### For more doing data analysis with Python:\n",
    "\n",
    "- Thomas Wiecki, [A modern guide to getting started with Data Science and Python](http://twiecki.github.io/blog/2014/11/18/python-for-data-science/)\n",
    "- Wes McKinney, [Python for Data Analysis](http://shop.oreilly.com/product/0636920023784.do) (book)\n",
    "\n",
    "### Packages building on xarray for the geophysical sciences\n",
    "\n",
    "For analyzing GCM output:\n",
    "\n",
    "- [xgcm](https://github.com/xgcm/xgcm) by Ryan Abernathey\n",
    "- [oogcm](https://github.com/lesommer/oocgcm) by Julien Le Sommer\n",
    "- [MPAS xarray](https://github.com/pwolfram/mpas_xarray) by Phil Wolfram\n",
    "- [marc_analysis](https://github.com/darothen/marc_analysis) by Daniel Rothenberg\n",
    "\n",
    "Other tools:\n",
    "\n",
    "- [windspharm](https://github.com/ajdawson/windspharm): wind spherical harmonics by Andrew Dawson\n",
    "- [eofs](https://github.com/ajdawson/eofs): empirical orthogonal functions by Andrew Dawson\n",
    "- [infinite-diff](https://github.com/spencerahill/infinite-diff) by Spencer Hill \n",
    "- [aospy](https://github.com/spencerahill/aospy) by Spencer Hill and Spencer Clark\n",
    "- [regionmask](https://github.com/mathause/regionmask) by Mathias Hauser\n",
    "- [salem](https://github.com/fmaussion/salem) by Fabien Maussion\n",
    "\n",
    "Resources for teaching and learning xarray in geosciences:\n",
    "- [Fabien's teaching repo](https://github.com/fmaussion/teaching): courses that combine teaching climatology and xarray\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Basic operations in `xarray`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-------------------\n",
    "\n",
    "## Import python packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# standard imports\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "\n",
    "import warnings\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "np.set_printoptions(precision=3, linewidth=80, edgeitems=1)  # make numpy less verbose\n",
    "xr.set_options(display_width=70)\n",
    "warnings.simplefilter('ignore') # filter some warning messages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic data arrays in numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "a = np.array([[1, 3, 9], [2, 8, 4]])\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a[1, 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a.mean(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "numpy is a powerful but \"low-level\" array manipulation tool. Axis only have numbers and no names (it is easy to forget which axis is what, a common source of trivial bugs), arrays can't carry metadata (e.g. units), and the data is unstructured (i.e. the coordinates and/or other related arrays have to be handled separately: another source of bugs).\n",
    "\n",
    "**This is where ``xarray`` comes in!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Properties of `xarray.Dataset` and `xarray.DataArray` objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start with the \"air_temperature\" tutorial dataset. This tutorial comes with the xarray package. Other examples [here](https://github.com/pydata/xarray/tree/master/examples)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ds = xr.tutorial.load_dataset('air_temperature')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.dims"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.attrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "type(ds.air.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.dims"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.attrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ds.air.attrs['tutorial-date'] = 27042017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.attrs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's Do Some Math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kelvin = ds.air.mean(dim='time')\n",
    "kelvin.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "centigrade = kelvin - 273.16\n",
    "centigrade.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice xarray has changed the colormap according to the dataset (borrowing logic from Seaborn).\n",
    "* With degrees C, the data passes through 0, so a diverging colormap is used\n",
    "* With Kelvin, the default colormap is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# ufuncs work too\n",
    "np.sin(centigrade).plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding Data to `DataSet`s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's add those kelvin and centigrade dataArrays to the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds['centigrade'] = centigrade\n",
    "ds['kelvin'] = kelvin\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.kelvin.attrs  # attrs are empty! Let's add some"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ds.kelvin.attrs['Description'] = 'Mean air tempterature (through time) in kelvin.'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.kelvin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ds.to_netcdf('new file.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Selecting data with named dimensions\n",
    "\n",
    "In xarray there are many different ways for selecting and indexing data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Positional indexing (old way)\n",
    "\n",
    "This is the \"old way\", i.e. like ``numpy``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air[:, 1, 2]  # note that the attributes, coordinates are preserved"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air[:, 1, 2].plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This selection implies prior knowledge about the structure of the data, and is therefore much less readable than the \"xarray methods\" presented below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Selection by index\n",
    "\n",
    "Selection based on the **index** of a coordinate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.isel(time=0).plot();  # like above, but with a dimension name this time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Selection by value\n",
    "\n",
    "Selection based on the **value** of a coordinate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.sel(lat=72.5, lon=205).plot(); "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Selection by value works well for time, too"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.sel(time='2013-01-02').plot(); # Note that we will extract 4 time steps! 3d data is plotted as histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.sel(time='2013-01-02T06:00').plot();  # or look at a single timestep"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Selecting a range of values\n",
    "The syntax is similar, but you'll need to use a [slice](https://docs.python.org/3/library/functions.html#slice):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.sel(lat=slice(60, 50), lon=slice(200, 270), time='2013-01-02T06:00:00').plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nearest neighbor lookup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.sel(lat=41.8781, lon=360-87.6298, method='nearest', tolerance=5).plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Operations and computation\n",
    "\n",
    "* We can do arithmetic directly on `Dataset` and `DataArray` objects. \n",
    "* Labels are preserved and dataArray dimensions automatically aligned."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Broadcasting\n",
    "\n",
    "<img src=\"./figures/broadcast.png\" width=\"50%\" align=\"left\"> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = xr.DataArray(np.arange(3), dims='time', \n",
    "                 coords={'time':np.arange(3)})\n",
    "b = xr.DataArray(np.arange(4), dims='space', \n",
    "                 coords={'space':np.arange(4)})\n",
    "a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Alignment\n",
    "\n",
    "<img src=\"./figures/align.png\" width=\"50%\" align=\"left\"> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "atime = np.arange(3)\n",
    "btime = np.arange(5) + 1\n",
    "atime, btime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = xr.DataArray(np.arange(3), dims='time', \n",
    "                 coords={'time':atime})\n",
    "b = xr.DataArray(np.arange(5), dims='time', \n",
    "                 coords={'time':btime})\n",
    "a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aggregation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.median(dim=['lat', 'lon']).plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Masking with `.where()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "means = ds.air.mean(dim=['time'])\n",
    "means.where(means > 273.15).plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Groupby and \"split-apply-combine\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Xarray implements the \"split-apply-combine\" paradigm with `groupby`. This works really well for calculating climatologies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.groupby('time.season').mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figures/split_apply-combine.png\" alt=\"split\" style=\"width: 50%;\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds.air.groupby('time.month').mean('time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "clim = ds.air.groupby('time.month').mean('time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also do arithmetic with groupby objects, which repeats the arithmetic over each group:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "anomalies = ds.air.groupby('time.month') - clim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "anomalies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "anomalies.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "anomalies.sel(time= '2013-02').plot();  # Find all the anomolous values for February"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resample adjusts a time series to a new resolution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tmin = ds.air.resample('1D', dim='time', how='min')  # Resample to one day '1D\n",
    "tmax = ds.air.resample('1D', dim='time', how='max')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(tmin.sel(time='2013-02-15') - 273.15).plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ds_extremes = xr.Dataset({'tmin': tmin, 'tmax': tmax})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds_extremes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Graphics "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "``xarray`` plotting functions rely on matplotlib internally, but they make use of all available metadata to make the plotting operations more intuitive and interpretable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1D plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zonal_t_average = ds.air.mean(dim=['lon', 'time']) - 273.15\n",
    "zonal_t_average.plot();  # 1D arrays are plotted as line plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2D plots "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_average = ds.air.mean(dim='time') - 273.15\n",
    "t_average.plot();  # 2D arrays are plotted with pcolormesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_average.plot.contourf();  # but you can use contour(), contourf() or imshow() if you wish"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Customizing 2d plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_average.plot.contourf(cmap='BrBG_r', vmin=-15, vmax=15);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_average.plot.contourf(cmap='BrBG_r', levels=22, center=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dealing with Outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "air_outliers = ds.air.isel(time=0).copy()\n",
    "air_outliers[0, 0] = 100\n",
    "air_outliers[-1, -1] = 400\n",
    "air_outliers.plot();  # outliers mess with the datarange and colorscale!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Using `robust=True` uses the 2nd and 98th percentiles of the data to compute the color limits.\n",
    "air_outliers.plot(robust=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Facet plots "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_season = ds.air.groupby('time.season').mean(dim='time') - 273.15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# facet plot allows to do multiplot with the same color mappings\n",
    "t_season.plot.contourf(x='lon', y='lat', col='season', col_wrap=2, levels=22);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting on maps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For plotting on maps, we rely on the excellent [cartopy](http://scitools.org.uk/cartopy/docs/latest/index.html) library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = plt.figure(figsize=(8, 4))\n",
    "# Define the map projection *on which* you want to plot\n",
    "ax = plt.axes(projection=ccrs.Orthographic(-80, 35))\n",
    "# ax is an empty plot. We now plot the variable t_average onto ax\n",
    "# the keyword \"transform\" tells the function in which projection the air temp data is stored \n",
    "t_average.plot(ax=ax, transform=ccrs.PlateCarree())\n",
    "# Add gridlines and coastlines to the plot\n",
    "ax.coastlines(); ax.gridlines(); "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Facet plots on maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# this time we need to retrieve the plots to do things with the axes later on\n",
    "p = t_season.plot(x='lon', y='lat', col='season', transform=ccrs.PlateCarree(),\n",
    "                  subplot_kws={'projection': ccrs.Orthographic(-80, 35)})\n",
    "for ax in p.axes.flat:\n",
    "    ax.coastlines()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Seaborn is Cool"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Statistical visualization with [Seaborn](https://stanford.edu/~mwaskom/software/seaborn/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "data = (ds_extremes\n",
    "        .sel_points(lat=[41.8781, 37.7749], lon=[360-87.6298, 360-122.4194],\n",
    "                    method='nearest', tolerance=3,\n",
    "                    dim=xr.DataArray(['Chicago', 'San Francisco'],\n",
    "                                     name='location', dims='location'))\n",
    "        .to_dataframe()\n",
    "        .reset_index()\n",
    "        .assign(month=lambda x: x.time.dt.month))\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "sns.violinplot('month', 'tmax', 'location', data=data, split=True, inner=None);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Out-of-core computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a quick demo of [how xarray can leverage dask](http://xarray.pydata.org/en/stable/dask.html) to work with data that doesn't fit in memory. This lets xarray substitute for tools like `cdo` and `nco`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's open 10 years of runoff data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`xarray`can open multiple files at once using string pattern matching.  \n",
    "  \n",
    "  In this case we open all the files that match our `filestr`, i.e. all the files for the 2080s. \n",
    "  \n",
    "  Each of these files (compressed) is approximately 80 MB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from glob import glob\n",
    "files = glob('data/*dis*.nc')\n",
    "runoff = xr.open_mfdataset(files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "runoff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`xarray` even puts them in the right order for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "runoff.time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How big is all this data uncompressed? Will it fit into memory?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "runoff.nbytes / 1e9  # Convert to gigiabytes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with Big Data\n",
    "\n",
    "* This data is too big for our memory. \n",
    "* That means we need to process it in chunks. \n",
    "* We can do this chunking in `xarray` very easily.  \n",
    "  \n",
    "  `xarray` computes data 'lazily'. That means that data is only loaded into memory when it is actually required. This also allows us to inspect datasets without loading all the data into memory.  \n",
    "  \n",
    "  To do this `xarray` integrates with `dask` to support streaming computation on datasets that don’t fit into memory.\n",
    "\n",
    "<img src=\"./figures/dask-array.png\" width=\"40%\" align=\"center\"> \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "runoff = runoff.chunk({'lat': 60})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "runoff.chunks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%time ro_seasonal = runoff.groupby('time.season').mean('time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import dask\n",
    "from multiprocessing.pool import ThreadPool\n",
    "dask.set_options(pool=ThreadPool(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%time ro_seasonal.compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dask.set_options(pool=ThreadPool(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%time ro_seasonal = runoff.groupby('time.season').mean('time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%time result = ro_seasonal.compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "brazil = dict(lat=slice(10.75, -40.75), lon=slice(-100.25, -25.25))\n",
    "result.dis.sel(**brazil).plot(col='season', size=4, cmap='Spectral_r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# xarray can do more!\n",
    "\n",
    "* concatentaion\n",
    "* open network located files with openDAP\n",
    "* import and export Pandas DataFrames\n",
    "* .nc dump to \n",
    "* groupby_bins\n",
    "* resampling and reduction\n",
    "\n",
    "For more details, read this blog post: http://continuum.io/blog/xray-dask\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Raw Cell Format",
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}