Simple WRF analysis/plotting examples

WRF_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a simple example for plotting WRF output using several standard Python libraries:\n",
    "\n",
    "1. [xray](http://xray.readthedocs.org/en/stable/) - netCDF reading and data manipulation\n",
    "2. [matplotlib](http://matplotlib.org/) - Python basic plotting library\n",
    "3. [cartopy](http://scitools.org.uk/cartopy/docs/latest/index.html) - Geographic plotting library build on matplotlib\n",
    "\n",
    "If you're already using Anaconda, then the easiest way to install these will be to simply go to your command line and execute:\n",
    "\n",
    "```bash\n",
    "conda install xray matplotlib cartopy\n",
    "```\n",
    "\n",
    "Alternatively, you could use `pip` or install them from source, but that's not recommended. Other libraries to consider would be [iris](http://scitools.org.uk/iris/docs/latest/index.html) for reading in/handling netCDF data, or [matplotlib basemap](http://matplotlib.org/basemap/) for geographic plotting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading in data\n",
    "\n",
    "Using `xray`, it's trivial to read in a netCDF dataset. I've attached a timeslice of output from a simulation of a mid-latitude squall line in the file [squall_slice.nc]() for this purpose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import xray\n",
    "\n",
    "dataset = xray.open_dataset(\"squall_slice.nc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normally, you'd use `ncdump -h squall_slice.nc` from the command line to inspect this file, but we can do something equivalent using xray."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xray.Dataset>\n",
      "Dimensions:          (Time: 1, bottom_top: 80, bottom_top_stag: 81, south_north: 100, west_east: 500)\n",
      "Coordinates:\n",
      "    XLAT             (Time, south_north, west_east) float32 0.0 0.0 0.0 0.0 ...\n",
      "    XLONG            (Time, south_north, west_east) float32 0.0 0.0 0.0 0.0 ...\n",
      "  * Time             (Time) int64 0\n",
      "  * bottom_top       (bottom_top) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 ...\n",
      "  * bottom_top_stag  (bottom_top_stag) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 ...\n",
      "  * south_north      (south_north) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 ...\n",
      "  * west_east        (west_east) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ...\n",
      "Data variables:\n",
      "    P                (Time, bottom_top, south_north, west_east) float32 11.1484 ...\n",
      "    PB               (Time, bottom_top, south_north, west_east) float32 96225.7 ...\n",
      "    PH               (Time, bottom_top_stag, south_north, west_east) float32 0.0 ...\n",
      "    PHB              (Time, bottom_top_stag, south_north, west_east) float32 0.0 ...\n",
      "    PSFC             (Time, south_north, west_east) float32 0.0 0.0 0.0 0.0 ...\n",
      "    QCLOUD           (Time, bottom_top, south_north, west_east) float32 0.0 ...\n",
      "    QRAIN            (Time, bottom_top, south_north, west_east) float32 0.0 ...\n",
      "    QVAPOR           (Time, bottom_top, south_north, west_east) float32 0.0196612 ...\n",
      "    T                (Time, bottom_top, south_north, west_east) float32 1.60981 ...\n",
      "    T2               (Time, south_north, west_east) float32 0.0 0.0 0.0 0.0 ...\n",
      "Attributes:\n",
      "    TITLE:  OUTPUT FROM WRF V3.6 MODEL\n",
      "    START_DATE: 2011-04-25_18:00:00\n",
      "    SIMULATION_START_DATE: 2011-04-25_18:00:00\n",
      "    WEST-EAST_GRID_DIMENSION: 501\n",
      "    SOUTH-NORTH_GRID_DIMENSION: 101\n",
      "    BOTTOM-TOP_GRID_DIMENSION: 81\n",
      "    DX: 1000.0\n",
      "    DY: 1000.0\n",
      "    STOCH_FORCE_OPT: 0\n",
      "    GRIDTYPE: C\n",
      "    DIFF_OPT: 2\n",
      "    KM_OPT: 2\n",
      "    DAMP_OPT: 0\n",
      "    DAMPCOEF: 0.003\n",
      "    KHDIF: 0.0\n",
      "    KVDIF: 0.0\n",
      "    MP_PHYSICS: 9\n",
      "    RA_LW_PHYSICS: 0\n",
      "    RA_SW_PHYSICS: 0\n",
      "    SF_SFCLAY_PHYSICS: 0\n",
      "    SF_SURFACE_PHYSICS: 0\n",
      "    BL_PBL_PHYSICS: 0\n",
      "    CU_PHYSICS: 0\n",
      "    SF_LAKE_PHYSICS: 0\n",
      "    SURFACE_INPUT_SOURCE: 1\n",
      "    SST_UPDATE: 0\n",
      "    GRID_FDDA: 0\n",
      "    GFDDA_INTERVAL_M: 0\n",
      "    GFDDA_END_H: 0\n",
      "    GRID_SFDDA: 0\n",
      "    SGFDDA_INTERVAL_M: 0\n",
      "    SGFDDA_END_H: 0\n",
      "    HYPSOMETRIC_OPT: 1\n",
      "    SF_URBAN_PHYSICS: 0\n",
      "    SHCU_PHYSICS: 0\n",
      "    MFSHCONV: 0\n",
      "    FEEDBACK: 1\n",
      "    SMOOTH_OPTION: 2\n",
      "    SWRAD_SCAT: 1.0\n",
      "    W_DAMPING: 0\n",
      "    DT: 5.0\n",
      "    RADT: 0.0\n",
      "    BLDT: 0.0\n",
      "    CUDT: 0.0\n",
      "    AER_OPT: 0\n",
      "    SWINT_OPT: 0\n",
      "    AER_TYPE: 1\n",
      "    AER_AOD550_OPT: 1\n",
      "    AER_ANGEXP_OPT: 1\n",
      "    AER_SSA_OPT: 1\n",
      "    AER_ASY_OPT: 1\n",
      "    AER_AOD550_VAL: 0.12\n",
      "    AER_ANGEXP_VAL: 1.3\n",
      "    AER_SSA_VAL: 1.4013e-45\n",
      "    AER_ASY_VAL: 1.4013e-45\n",
      "    MOIST_ADV_OPT: 1\n",
      "    SCALAR_ADV_OPT: 1\n",
      "    TKE_ADV_OPT: 1\n",
      "    DIFF_6TH_OPT: 0\n",
      "    DIFF_6TH_FACTOR: 0.12\n",
      "    OBS_NUDGE_OPT: 0\n",
      "    BUCKET_MM: -1.0\n",
      "    BUCKET_J: -1.0\n",
      "    PREC_ACC_DT: 0.0\n",
      "    SF_OCEAN_PHYSICS: 0\n",
      "    ISFTCFLX: 0\n",
      "    ISHALLOW: 0\n",
      "    ISFFLX: 0\n",
      "    ICLOUD: 0\n",
      "    ICLOUD_CU: 0\n",
      "    TRACER_PBLMIX: 1\n",
      "    SCALAR_PBLMIX: 0\n",
      "    GRAV_SETTLING: 0\n",
      "    DFI_OPT: 0\n",
      "    WEST-EAST_PATCH_START_UNSTAG: 1\n",
      "    WEST-EAST_PATCH_END_UNSTAG: 500\n",
      "    WEST-EAST_PATCH_START_STAG: 1\n",
      "    WEST-EAST_PATCH_END_STAG: 501\n",
      "    SOUTH-NORTH_PATCH_START_UNSTAG: 1\n",
      "    SOUTH-NORTH_PATCH_END_UNSTAG: 100\n",
      "    SOUTH-NORTH_PATCH_START_STAG: 1\n",
      "    SOUTH-NORTH_PATCH_END_STAG: 101\n",
      "    BOTTOM-TOP_PATCH_START_UNSTAG: 1\n",
      "    BOTTOM-TOP_PATCH_END_UNSTAG: 80\n",
      "    BOTTOM-TOP_PATCH_START_STAG: 1\n",
      "    BOTTOM-TOP_PATCH_END_STAG: 81\n",
      "    GRID_ID: 1\n",
      "    PARENT_ID: 0\n",
      "    I_PARENT_START: 1\n",
      "    J_PARENT_START: 1\n",
      "    PARENT_GRID_RATIO: 1\n",
      "    CEN_LAT: 0.0\n",
      "    CEN_LON: 0.0\n",
      "    TRUELAT1: 0.0\n",
      "    TRUELAT2: 0.0\n",
      "    MOAD_CEN_LAT: 0.0\n",
      "    STAND_LON: 0.0\n",
      "    POLE_LAT: 0.0\n",
      "    POLE_LON: 0.0\n",
      "    GMT: 0.0\n",
      "    JULYR: 0\n",
      "    JULDAY: 1\n",
      "    MAP_PROJ: 0\n",
      "    MAP_PROJ_CHAR: Cartesian\n",
      "    MMINLU: \n",
      "    NUM_LAND_CAT: 24\n",
      "    ISWATER: 0\n",
      "    ISLAKE: 0\n",
      "    ISICE: 0\n",
      "    ISURBAN: 0\n",
      "    ISOILWATER: 0\n",
      "    history: Wed Oct 14 10:17:31 2015: ncks -d Time,0 -v PH,PHB,T,P,PB,PSFC,T2,QRAIN,QVAPOR,QCLOUD,XLAT,XLONG squall3d_03z.nc -o squall_slice.nc\n",
      "    NCO: 4.4.7\n"
     ]
    }
   ],
   "source": [
    "print(dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can look into more detail at a given field, too, to see what metadata is attached to it. Note that the default container that `xray` provides, the **Dataset**, implements the standard dictionary interface, so it's easy to grab data fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xray.DataArray 'QRAIN' (Time: 1, bottom_top: 80, south_north: 100, west_east: 500)>\n",
      "[4000000 values with dtype=float32]\n",
      "Coordinates:\n",
      "    XLONG        (Time, south_north, west_east) float32 0.0 0.0 0.0 0.0 0.0 ...\n",
      "  * bottom_top   (bottom_top) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ...\n",
      "  * Time         (Time) int64 0\n",
      "  * south_north  (south_north) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ...\n",
      "    XLAT         (Time, south_north, west_east) float32 0.0 0.0 0.0 0.0 0.0 ...\n",
      "  * west_east    (west_east) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ...\n",
      "Attributes:\n",
      "    FieldType: 104\n",
      "    MemoryOrder: XYZ\n",
      "    description: Rain water mixing ratio\n",
      "    units: kg kg-1\n",
      "    stagger: \n"
     ]
    }
   ],
   "source": [
    "rain = dataset['QRAIN']\n",
    "print(rain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To illustrate simple data manipulation, let's compute total rain amount for each model column in the `y` dimension, and plot against `x`. This will give us an idea of where the squall line is located in the model domain. \n",
    "\n",
    "Rather than use giant loops over all the model dimensions, we'll use the convenience tools exposed by `xray`. See [here](http://xray.readthedocs.org/en/stable/data-structures.html?highlight=pipe) for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = 'west_east'\n",
    "y = 'south_north'\n",
    "z = 'bottom_top'\n",
    "t = 'Time'\n",
    "\n",
    "total_rain = (rain.load()\n",
    "                 .sum(dim=[y, z], keep_attrs=True)\n",
    "                 .squeeze())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting simple data like this is very straightforward, but as a reminder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x11309b810>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAADICAYAAADRNycrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VeWZwPFfgCwkIAgRIojI+uAKiAujFVHc61hrW51q\nVazVQgE7dVwGWy3aaWvdOnUBmQ6gtmirjor7VhBatbgVtYAPIBBC2AMYSAKBkPnjPefm5pLlJNwl\n597n+/nkk3vfe+65L68mT553zaqtrcUYY4xJlHaproAxxpj0ZoHGGGNMQlmgMcYYk1AWaIwxxiSU\nBRpjjDEJZYHGGGNMQnVIxYeKyEnA86rau5HXvwv8EugBzAOuVdVNB/B53YFJwEOqWtba+6QDa4s6\n1hZ1rC3qWFvUiVdbJDWjEZEsEfk+8CaQ3cg1xwHTgMuAQmADMOsAP7o78HPve6aztqhjbVHH2qKO\ntUWduLRFsrvObgNuAP4LyGrkmiuAF1T1Q1XdBdwKnCcihySpjsYYY+Io2YFmhqoOAz5q4hoBlvhP\nVHUrsNUrN8YYEzJJDTSquiHAZQVAZUxZJZAf/xoZY4xJtJRMBmhGQ0ElH9gZ5M3e4FVsf+JA73sf\nEak+sOqFnj8Bo7e1hbVFFGuLOtYWdfp43weK7NepVBZ0gkBbDDRLieomE5FCoJtXHsQk3OBVQ94+\nsKqllbmprkAbYm1Rx9qijrVFnVcaKLsTmBLkzW0x0DwFzBeRmcDHwK+BV1V1W8D3PwQ8GVPWB3h7\n9uzZFBUVxa+mxsTJI488wiGHHMKll16a6qoYE7FhwwauuOIKgLOAkpiXA093TmWgiZxPICLTAFR1\nvKp+KiLXATOBImABcE3Qm3qpXL0G8NPfoqIiDjvssDhU3Zj42bx5M3fffTcAV199Nd2726xa0+Z8\nqaqrW/vmlAQaVX0HtxjTfz4+5vVngGeSXC1jUqKysm7uy3PPPcd1112XwtoYE3+2BY0xKVZdXTfe\n/P7776ewJsYkhgUaY1IsOtCUlpamsCbGJIYFGmNSLDrQrFu3LoU1MSYxLNAY04jJkyfzzW9+s14g\nSITdu3dHHlugMemoRYFGRHJF5DAR6SUi7RNVKWNS7S9/+Qt33303L7zwAgsWLIiUT506laKiIubO\njd8Si+hAtnXrVqqqquJ2b2PagmYDjYj0FZHfiMhnuFX7a4C1QLWIfCIivxCRvomuqDHJ9PLLL0ce\nv/22W+e7atUqJkyYwMaNG5kzZ07cPis2Y1q/fn3c7m1MW9BooBGRbiIyA1gE9AN+B5wGHAkcA5yB\n275fgM9EZJa3/Ysxobd169bIYz/QLFq0KFK2ZcuWuH1WdNcZ2IQAk36aWkczF3cuzI9UdXcj1ywA\nHhKRg4CxwDvAsfGsoDGpEB1oPvnkE8rKyli+fHmkLJ5jKbEZjY3TmHTTVKD5F1UN1FmsquXAgyLy\nv/GpljGptW2b2/EoOzubPXv28NZbb7Fs2bLI64kMNJbRmHTTaNdZ0CAT857Y7f2NCSU/oxk7diwA\nU6ZM4Z///Gfk9aCBZsOGDdx4443Mnz+/0Wtiu84sozHpxqY3G9MAP6OZPHkyQ4YMQVVZuHBh5PWd\nO3eyY8eOZu8za9Ysfvvb3zJ69GhWrlzZ4DV+RpOTkwNYRmPST6NdZyJyLlEbXzZFVd+MW42MSbHa\n2tpIRlNUVMT06dM5/fTTARg5ciRbtmxhxYoVrFu3jgbO6KinuLg48njt2rX0799/v2v8QNOvXz9U\n1TIak3aaymgeAF4P+GVM2qiqqqK6upq8vDw6duzIqFGjePjhhznooIO444476NWrFxAs8ygpqdtZ\nvaKiosFrogNN0PsaEyZNTQYYgdtBuTdwiqruSk6VjEktP5s5+OCDI2UTJkzgRz/6EVlZWTzxxBNA\nsLGUtWvXRh43Fmj8MZqBAweSm5vLypUrKSkpoU+fPg1eb0zYNDUZYBdwGZBLwFPUjEkH/vhMt27d\n6pVnZWUBRDKaeAUaP6Pp2rUrF110EbW1tcyePbvlFTemjWpyMoA3i+wqWnCSmjFh11BGE613b3ek\nfHOBprKyst56nOYCTU5ODueffz4AixcvblmljWnDmj34TFU/xh2pbExGaCyj8QXNaFavXl3veXNd\nZ7m5ufTo4c4DjOfOA8akWotP2BSR/wSmq+q2Vrx3ODAdOApYDoxT1YUNXDcZGA8cBPwTuEFVP2np\n5xnTGs1lNEEnAyxdurTe8yAZTWFhIeCOdzYmXbRmHc1PgYZ/ApsgInnAS8AMoAvwIPCiiBTEXHcm\ncBNwpqp29d5jxzqbpIlXRrNkyZJ6z4MEmkMOOQSwQGPSSzIXbJ4B1KjqdFWtUdVZwEbggpjrdnrf\ns72jCPbhdo02JimCZjTr1q2jtrbxpWZ+oDnhhBOAYF1nfkZjXWcmnSQz0AwBlsSUqVdeV6D6AfAI\nsBjYBUwGrkhGBY2B5jOa/Px8unbtSnV1db3B/lj+3mjNBZrojKZz587k5ORQWVlJZaX9fWXSQ2sC\nzThgUyveV8D+mUklkB9dICLfBq4HTvDe89/A817XmzEJ5weaxjIaCNZ95p8rM3jwYCBYoMnKyop0\nn1lWY9JFoMkAInJ41NO/At1EpBtui5pqYLOq7mvmNhVAx5iyfCB2w6jvAY9GDf7fJSLXAWcBL9MM\n70yc2HNxejf3PmN8fpbSWEYDLtAsWbKE0tJSjj12/5Mxampq2LTJ/T3mbzvTXKDJzc0FoLCwkNLS\nUjZv3szhhx/e4HuMSbJ+IpITU1amqoGWvgSddbYSyPK+otV6ZdUi8n/A9ara8E8TLAUmxpQJELsy\nrQqIzV5qgD0B6zoJ+HnAa43ZTzwymrKyMmpqaujevTtdu3YFmh+j8TfV9Kc4b9iwoRW1NyYhGjq7\n/E4CLuYPGmh+gJttNhH4u1d2IvAQ8BgwH7gHuA83LbmxiuaKyETcFOcrgR7AGzHX/QmYISJ/Bj4H\nbsB18f0tYF0fAp6MKetNww1lzH6amwwAzQcaP0gUFRWRn+96h4N0nQGRLGbNmjUtrboxiXImEDuf\nP/BC/qCB5i7ge6q6IKrsbRG5FnhKVX8jIj8BXqORQKOq1SJyPvAo8CvcOpqLVLVKRKZ514xX1Tki\nUgQ8jesC+wdwXhOZUuznlBHTACJS3cjlxtRTW1sbCRI9e/Zs9LqWBJqCAjeDv7FAs2uX20bQ7zrr\n27cvUH/nZ2NSbJWqrm7tm4MGmi7sP5YCblZYofd4GzED+7FU9XPg1AbKx8c8n47LeoxJqh07dlBR\nUUF+fj4HHXRQo9c1tw1NSwLN9u3bASJdbH6gsYzGpIugs85eBx4VkaP8AhE5EngYeN0bJBoPLIp/\nFY1JHn+1f+/evSObaDakud0B/BlnPXv2bHWgsYzGpIugGc044CngnyJSiZsA0BG3av964DzgauCi\nRFTSmGTxMxQ/Y2lMc11ny5cvB2DAgAGRQNPYupjYQOOP0VigMekiUKDx9jU7T0QGA8fhZoAtVtUV\nACLyFlAUYIqzMW2an6H4gaQxRUVFgOsiq6mpoX379vVe9xdrDh48mLy8PLKysti9e/d+1+7Zs4eK\nigratWtHp06dADjssMNo164d69ato7q6OjJJwJiwCrqO5nxVfU1VlwHLosr7AY+oauw2MsaEUtCM\nxt+XbPPmzWzatIlDDz203uuqCoCIkJWVRUFBATt37qSiooLOnTszceJE+vXrx9ixYwGXzfhdddnZ\n2fTq1Yu1a9dSWloaOXnTmLAKOkbznIh8y38iItki8lPcNjE9ElIzY1LADxBBFko21n22fft2Nm3a\nRH5+fiRgRY/TFBcXM3XqVG6++ebIe2OnUts4jUknQQPN1cATInKliJwOfAr8BLgRt57GmLTw7rvv\nAjBy5Mhmr21s5tnKlSsBNz7Trp37EYsONFVVVZFrX331VaBufMZn4zQmnQQdo3laRLYCz+EmAcwE\nJqtq4zsKGhMyGzduZPny5RQUFDBs2LBmr4+debZjxw7KysoiB575WQnUDzR79+6NlE+ePBnYP9BY\nRmPSSaOBxhv4j7YGmAD8D1AMFIpIIYA3dmNMqH366acAjBgxgg4dmv8bLLrrbM+ePYwaNYqlS5cy\nbtw4AI444ojItdGBpqamZr97xQ74++MyK1asaPk/xJg2pqmfpi+aeO2/vC9w+521b+JaY0LBn5Is\nIoGujw40jz32GIsWuWVkM2fOBBrPaPbtc5MzO3XqxM6d7vil2MWZw4cPB+CDDz5o1b/FmLakqTGa\n/gG/BiS4jsYkhZ89DBw4MND10YHGH2sB14UGjQcaf+HmWWedxfz58ykqKmLKlCn17j106FBycnJQ\nVaZPt00yTLg1mtEcyL42xoRRSwONPxlg9erVDe603Fig8U/lzM/PZ9SoUZFdBKLl5OQwYsQI3n//\nfcaNG8cPfvCD/dbqGBMWyTxh05g2raWBRkRo3749S5cuZdu2bfutpYk+pyY60Pg7BPhljZk6dWrk\n8apVqwLVyZi2yAKNMbhdm/3ZYtGD+E0pKCioF0zOOeeceq937Nix3rVQP9D4xwc0ZtiwYZx11lkA\nLF26NFCdjGmLLNAYgzuobNeuXXTp0qXJXZtjnXzyyZHHo0eP5jvf+Q4A5513Xr3rGhqjaS7QABx5\n5JGABRoTbkE31TQmrfmzvlp6dPKkSZN47rnn2LNnD+eeey6XXHIJw4YN44c//GG966IDjb+IM0ig\n8WfA2TRnE2ZB9zpbRd2xzdH24TbYXAv8WVV/H9/qGZMcJSUlQMsDzdFHH01JSQnV1dV07twZgNtu\nu22/6/wFmdu2bYt0qQUJNN26dQPgq6++alG9jGlLgnadPQQcgjtm+d+9r9m4Q89eBd4GpojIfySi\nksYkmp/R9OnTp8Xvzc3NjQSZxvTo4bYE3LRpU+AxGiByX3/KtDFhFLTr7EpgvKr+Mapsjoh8Btyi\nqieIyCe4Y5rvb+wmIjIcd3LmUbijnMep6sIGrjsN+B0wCFgF/FhV5wWsqzEt5m/10ppAE0R0oPG7\n0ZqbdQYWaEx6CJrRCPBhA+WfAcd4jxU4tIFr3A1E8nAHpc3AHQ39IPCiiBTEXNcLmAP8QlU7A7/C\n7R6dG7CuxrTYkiVLABgyZEhC7t+zZ0/A7afWkskAFmhMOggaaD4BbhGRSAbkPf4JbidngFNxe6A1\n5gygRlWnq2qNqs4CNgKxZ9lcBbypqs8DqOqfvPfWBqyrMS22ePFiAI455phmrmwd6zozmSxo19kE\n4E1gtYj8A7e32XHe96+LyCjgMdyxzo0ZAiyJKVOvPNpwoFREngNG4Q5a+7GqVgesqzEtUl5ezpo1\na8jNzWXAgMTsqNS1a1eys7MpLy9n61a36XlLAo2/J5oxYRT0mIBPvd2cL8MFmGrgWeApVa0SkSOA\nE1X10yZuUwDEHppeCcT+tHXHZTnfBL6DC16viMhgVd3eXF1FpLt3j2hNH5doMta+ffsie4kdddRR\nCdvmJSsrix49elBaWsqXX34JEDm6uSmW0Zg2op+IxJ4pXqaqZUHeHHR6c19VLcYdERBd3kFEblHV\newLcpgJ3lk20fCD2J2gX8Iqqvu09nyYiN+O65l4J8DmTgJ8HuM5kuNraWi677DKeffZZgMhiy0Tx\nA42f0XTp0qXZ9+Tn59OuXTuqqqrYu3dvoOMLjEmAuQ2U3QlMCfLmoGM0f4s9n0ZEvgb8A7g94D2W\n4iYV1LsNDXen5cWUteTPzIe8+0Z/ndmC95sM8fHHH0eCTK9evSLnyCSKvybGF3vYWUOysrIimY91\nn5kUOpP9f68+FPTNQf88mgPMF5FzgFLgXmAs8AxwXhPvizYXyBWRibgpzlcCPYA3Yq77A/C+iFwA\nvI4bH8oFAk1v9lK5eumciNj4jtmPv+Pyueeey8svv5zwbCE2sATJaMB1n5WXl7Njx45AwcmYBFh1\nIDv6B8poVHUibo3MO7jM5GTgLFX9N1UtDXiPauB84Lu4QDABuMgb45kmItO86xYBF+EOVtuOm4X2\nr6oaO75jzAEpLy8H4OCDD05Kl1R0kOjYsSPZ2dmB3mfjNCbsAv90qeqdIrIRt/7litYsoFTVz3Fj\nLbHl42OevwW81dL7G9MS/rYuQTOLAxUdaFrymRZoTNg1GmhEZP/TmOq8IiJbvce1qtorvtUyJvH8\njKYluzUfiNYGGn+MxgKNCaumMprJAe9hCylNKKUyo2nJWIsfaPyFnsaETVOB5glV3deSm4lIlqpa\n4DGhEJaMxt/tuaqqKu51MiYZmpoMsFBELg5yExFpJyKXAh/Hp1rGJF5Yxmj8HQQsozFh1VRGcynw\niIg8CLyA24JmMbAFF6AKgWG4bWIuxe15ltgVb8bEUdgyGgs0JqwaDTSqugq4QEROAiYCM3HBJdoG\n4DXcNOWGdnc2ps1KZUbTkgPW/IzGus5MWDU7vVlVPwCuEpEsoC9ukeU+XJAptTEZE1bJzmgOPbTu\nFI3x48c3cWV91nVmwq4l62hqgdXelzGhl+yMpnv37syfP58ePXpQWBjbOdA46zozYWc79JmMleyM\nBmDUqFEtfo91nZmwC7qppjFppba2NukZTWtZRmPCzgKNyUi7d+9mz5495OTkkJvbtk8JtzEaE3YW\naExGCks2A9Z1ZsIv6MFnA4G7gROBbCAr6mXb68yETirGZ1rLus5M2AWdDDATOAS4j/1PxLTpzSZ0\n/EATpozGAo0Jq6CB5kTgJG+bf2NCz+86C0NGY11nJuyCjtGsAdr+T6QxAYUpo7GuMxN2QTOa/8Tt\ne3YXsAyodzSyqi6Ld8WMSaQwZjQWaExYBQ00/+d9f7aB12qB9kFuIiLDgenAUcByYJyqLmzi+jG4\nzTw721HOJp7CmNFY15kJq6BdZ/2b+BoQ5AYikge8BMwAuuCOhH5RRAoauf5g3CQEY+LOMhpjkidQ\nRqOqq+PwWWcANao63Xs+S0R+AlwAPNPA9dOAp4Bb4vDZxtQTpowmLy8PgF27dqW4Jsa0TqOBRkTW\nA8eoapn3uDFB19EMAZbElKlXHvvZV+AmH0zDAo1JgLKyMqBlRyqnSk5ODgB79uxh3759tGtn66xN\nuDSV0UwGKqIeH6gCIDb3rwTyowtE5HDgLuBUIK+lHyIi3YHuMcW9W3ofk942bNgA1N+6v63Kysoi\nNzeX3bt3s3v37siYjTFJ1E9EcmLKylS1LMibmzr47LGGHh+ACiD2JySfqAWgItIOeBz4qapuEJF+\n3ktZBDcJ+PmBVNSkv/XrXZJeVFSU4poEk5eXZ4HGpNLcBsruBKYEeXPQLWgOAiYAR1M3gSALl3EM\nU9V+jb03ylLcSZ31bg3Mjnp+GHAyMExEpkV91loR+bqqvhfgcx4Cnowp603DDWUyVJgyGiCy8aeN\n05gUORMojSkLlM1A8OnNvwdGA28BlwJ/AgYDx7F/8GjMXCBXRCbipjhfiTut8w3/AlVdQ1RXmoj0\nBVYBvYNOb/ZSuXoNICLVjVxuMlBNTQ2bNm0CoEePHimuTTA2IcCk2KoDmRQWdFTxHOBKVf0ebkD/\nAVU9CRcwgmQzqGo1cD7wXVwgmABcpKpVIjLNy2BiZWF7qZk427JlCzU1NRQWFkYG2ts6P6PZvXt3\nimtiTMsFzWg64rq+ABYDI4BPgKnAPOCnQW7i7ZV2agPlDR6g7kXQQItBjQnK7zYLy/gMWEZjwi1o\nRrMCOMV7vBQ3jgJujKbtr3gzJsq6deuAcAYay2hMGAXNaO4FHheRDsCfgUUiUguMBP6aqMoZkwhr\n164FoE+fPimuSXA2GcCEWaCMRlUfB84ClqqqAhcBRcC7wDWJq54x8ecHmsMOOyzFNQnOMhoTZkGn\nNz8N3O4FGVT1LdwMNGNCJ4yBxjIaE2ZBx2jGAHsTWRFjkiWMgcYyGhNmQcdoHgBmish/AyuBevuV\n23k0JkxKS926szAGGstoTBgFDTS/8L6f1sBrgc+jMaYt8DOa3r3DswWeraMxYRY00PRPaC2MSZId\nO3bw1VdfkZeXR7du3VJdncAsozFhlszzaIxJuehus6ysluzVmlo2GcCEmR1sYTJKGCcCgE0GMOFm\ngcZklLAHGstoTBhZoDEZJayBxiYDmDBr6ijnwNvaejszG9PmrVixAgjX9jNgGY0Jt6YmAwT9P9qm\nN5vQWLBgAQAjR45McU1axjIaE2ZNBZozk1YLY5Jg9erVrFq1ii5dujB06NBUV6dF/IymqqqqmSuN\naXsaDTSq+k6QG4hIoIPPjEmlr776iu9///sAnH322bRvH64kPD/fHTxrgcaEUdBNNY8B7geOxk0g\nyPK+cnHn0QT6qRWR4bhTOY8ClgPjVHVhA9ddB9wM9AQUuFFV/xbkM4xpyO233868efPo2rUr9913\nX6qr02IdO3YELNCYcAo662wa0BmYAhyM25LmCe+1/U7MbIiI5AEvATOALsCDwIsiUhBz3RnAL4Fv\nq2oX4GHgJREJzzJu0+a88cYbAMyePZu+ffumuDYt5weaysrKFNfEmJYLGmhGABNV9X9xRzgvVdVb\ngMnApID3OAOoUdXpqlqjqrOAjcAFMdf1Bu5R1c8AVPUJoAaXBRnTYuvXr2fZsmUUFBRw9tlnp7o6\nrWIZjQmzoHud7QO2eo8VGArMA14HfhPwHkOAJTFl6pXXFaj+Mfq5iJyKy6Zi32tMIB999BHgZppl\nZ2enuDatY2M0JsyCZjQfAdeJSBbwKXCeVy4EP6emAIjN+yuB/MbeICJHAc/iDl3b2th1xjTFXzsz\nZMiQZq5su6zrzIRZ0IzmVuA1XFYzC7hVRFYAhwKPBbxHBdAxpiwf2NHQxSJyDvAn4D5VvSfgZyAi\n3YHuMcXh2Q/exJ0faAYMGJDimrSedZ2ZFOvXwCL+MlUtC/LmQBmNNzOsL/BHL7M4ATdB4FqCj9Es\nxWVA0YQGusRE5BrgGdystF8FvL9vEq5LLvprbgvvYdKIH2gGDhyY4pq0nnWdmRSby/6/V4P+7g88\nvXkucImqbgRQ1XXA/SJyCK5b7fiAFc0VkYm4Kc5XAj2AN2I+awzwCHC2qr4b9B8S5SHgyZiy3liw\nyVjpEGis68yk2JlAaUxZoGwGmt7r7ELgFNx6mdHAnSJSEXPZIFym0yxVrRaR84FHgV/h1tFcpKpV\nIjINqFXVHwG3ANnA6yL1EqBvqeqbAT6njJgGEBHbiy1D7dmzh+LiYrKysujXL7xri6O7zmpra0N1\nlo5JC6sO5FyypjKaxcB/4AINuO6y6F/YtcBO4KqgH6aqn9PAuhtVHR/1+Nyg9zOmOcXFxdTU1NCn\nT5/INi5h1L59e3Jycqiurqa6ujqy95kxYdDUFjSrcGtfEJHHgBtUtTxJ9TImLtKh28zXsWNHqqur\nqaystEBjQiXoUc5jRaSd1512FG7LGQVeVVXbt9y0WV9++SWQPoHmq6++oqqqioMPPjjV1TEmsKCT\nAfoALwMDcAGmAzAQWCcio1U1dpDImDYhHaY2+/yZZzYhwIRN0AWbDwEbgMNVdYSqDsVNAlgF/Hei\nKmfMgUq3rjOwKc4mfIIGmjHAzdGr81V1C26H5XMSUTFj4iHdus7AAo0Jn6CBppyGt4rJx+2DZkyb\nU1NTEwk01nVmTOoEDTTPA49459IAICLH4RZWzklExYw5UKWlpVRXV9OzZ086deqU6uocsIICd6LG\nddddx+eff57i2hgTXKOBRkTuiDor5jbcPmefiUiFt3BzEVAM/Hviq2lMy33xxRdAenSbAfTq1Qtw\n406nnHJKimtjTHBNzTqbglvFX+GtnzlbRI4FjgR2AV+o6rLEV9EkU21tLQ8//DC9evXikksuCfUK\n9AULFgBw8sknp7gm8RF9YNvOnTtTWBNjWibo7s1AZGW/5expbMmSJdxwww0APP7441x1VeCNH9qc\nv/zlLwCceeaZKa5JfBx++OGproIxrdJcoMltYGvo/aiq7SWWJpYvXx55/OKLL4Y20JSXl/Phhx/S\noUMHRo0alerqxEXsEdTbt2+na9euKaqNMcE1NxmgGNdN1tSXzbVMIytXrow8/uCDD1JYkwOzYMEC\nampqOOmkk+jcuXOqqxMXsRlN9H8rY9qy5jKabwHbklER0zZE//IqKSlhw4YNFBUVpbBGrZNu3Wbg\nAk3fvn0pLi4GYPXq1Rx/fJATOoxJreYCzbuquikpNTFtwqpVq+o9X7x4cagDzZgxY1Jck/jp0KED\nS5cuZdKkScyYMYN169alukrGBBJ0HY3JEH5GM2LECKD+mE1YbNq0ic8//5y8vDxGjhyZ6urEVceO\nHTniiCMAt07ImDBoKtA8gRuDMRli3759kYzm3HPdsUBhDDTz5s0D4Gtf+1qoz6BpjL+exjIaExZN\nnUczNon1MG3A+vXr2b17N4WFhQwfPhwIZ6BJx26zaL179wYsozHh0aJ1NAdKRIYD03Fn2iwHxqnq\nwgau+y7wS6AHMA+41saKEs/PZvr378/gwYMB+OSTT9i7dy8dOiT1f5VWq6qq4vXXXwfSayJANMto\nTNgkbYxGRPKAl4AZQBfgQeDFqG1u/OuOA6YBlwGFuOMJZiWrnpnMH5/p378/xxxzDIMGDaK0tJTf\n/e53Ka5Z80pKSrjvvvsYPnw4JSUlDB48OG1nZPkZTXFxcWSHgOLiYkpKSlJZLWMalczJAGcANao6\nXVVrVHUWsBG4IOa6K4AXVPVD7/TOW4HzROSQJNY1Y+zbV7f59uLFiwEXaNq1a8ftt98OwE033cQt\nt9yS8Lp88cUXXHPNNbz33nuBrt+8eTOPPvoop59+Oocffjg333wzqkrv3r15+umnQ5OFtVS3bt0Y\nOXIklZWV3H///VRUVHD88cczdOhQO0LAtEnJ/EkcAiyJKVOvPJoAkd80qrpVRLZ65ZsTWsMM8+ST\nTzJx4kQGDx7M2rVrI33+fpfTlVdeSU1NDT/84Q+59957GTRoENdeey3t2sX/75OKigpOPPFEdu7c\nyTvvvIOqkpNTf1OKyspKFi1axF//+ldeeukl3nvvPWprawHIy8vjwgsv5PLLL+eCCy4gNzc37nVs\nS375y1+bPLPWAAALOUlEQVQyZswYpkyZQrt27di61R0V9cYbb3DxxRenuHbG1JfMQFMAxB6kUcn+\n59wEva7F3n777Xpnrfu/pPzv6V62a9cuKioq6N69O9u3b+dnP/sZe/fuZeHCumGyrl271tuyZezY\nsWzbto0bb7yR66+/nrvvvpsTTjiBwsJCCgoK6NSpU+R7p06dyM/PjwSihuoTW7d9+/ZRU1PDnDlz\nIt1Aq1ev5sILL4x03a1bt47S0lI2bNhQLwPLzs5mzJgxXH755XzjG9/goIMOIlOMHj2aQw89lPXr\n13PHHXdEyh944AFqamoiz6M3RfUfBy0zJl69AskMNBVAx5iyfGBHTFlDQSUfCLRdrYh0B7rHFPcB\nGD9+PHv37g1U2UzRoUMHBg4cSLdu3RgwYACnnHIKGzdurHfNt7/9bWpra3nwwQdZs2YNa9asSVhd\nLr/8cp555hnmzZsXmabsy87OZtCgQQwdOpTRo0dz2mmnRc5oKS8vp7y8PCH1aqtuuukmbr311npl\n77//Pu+//36KamTSTYcOHfytjwY0sO9lmaqWBblPVvRfmokkIucBj6jqgKiyz4A7VPWFqLK7gUNU\n9VrveSFuLKdQVZvdDkdEpgA/j3P1jTHG1Henqk4JcmEyM5q5uN2gJ+KmOF+Jm778Rsx1TwHzRWQm\n8DHwa+DVIEHG8xDwZEzZQOAV4Czgy9ZVP230w/23OBNY1cy16c7aoo61RR1rizoDgLeBrwMrYl4L\nlM1AEgONqlaLyPm4w9R+hVtHc5GqVonINO+a8ar6qYhcB8wEioAFwDUt+JwyYhpARPyHJaq6+kD/\nLWEWlf6WWltYW/isLepYW9SJaosVB3LQZVLnf3oHp53aQPn4mOfPAM8kq17GGGMSxzbVNMYYk1AW\naIwxxiRUpgSaMuBOWjB4lcasLepYW9SxtqhjbVEnLm2RtOnNxhhjMlOmZDTGGGNSxAKNMcaYhLJA\nY4wxJqEs0BhjjEkoCzTGGGMSygKNMcaYhLJAY4wxJqHS86zbKCIyHLdb9FG4jTzHqerCpt8VbiJy\nEvC8qvb2nh+M26T0DOAr3PbeM6Ou/zVwLe7/hyeAG1V13343DhkR+RpwP+501i3APar6P5nYHiJy\nKW7h3WFAMfBTVZ2TiW3hE5GewOfANar6Sia2hYjchNvkeHdU8Xm405Dj1hZpndGISB7wEjAD6AI8\nCLwoIgUprViCiEiWiHwfeBPIjnrp90A57liGbwP3iMjJ3nsmAhcAxwJH4jY9/Y9k1jsRvF8aLwK/\nVdWuwHeAX4vIGDKsPURkMO6XxjWq2hn4MfBn75DAjGqLGDOAboC/aj0T22IY8J+q2jnq613i3BZp\nHWhw0bhGVaerao2qzsIdonZBiuuVKLcBNwD/BWQBiEgn4BvAz1W1WlU/xJ3Xc5X3nitxv4w3qupG\n3Pk/Y5Nd8QQ4HHhJVf8EoKr/AOYBp5Bh7eFt795DVf8uIh1wx2+UA9VkWFv4RGQc7tTeEu95pv6c\nDAc+jS5IRFuke6AZgksBo6lXno5mqOow4KOoskHAnphzNZZR1wZC/TZa5pWFmqp+qqpX+8+9DOc0\nXADOxPaoFJF+wC5cV8dPcQcCZlxbeBnejUD08SQZ93MiIvm4f8OPRWS9iCwRkWtIQFuke6ApACpj\nyiqB/BTUJeFUdUMDxQVAVUxZJdAx6vXKmNfaNXA+eGiJSBdcF+pHuKwmU9tjDZCLO2n2AeBCMqwt\nvIzuCWBizKm9mfhz0gP4KzAV6ANcj/v/4uvEuS3SfTJABXWN48sHdqSgLqlSCeTFlOXjug381zvG\nvLZXVauTULeE8/6Kfxk3EeQy4GgytD1UtcZ7OE9E/g84gcxri9uBRar6pohkeWVZZODPiZexnBFV\n9DcR+QMwiji3RbpnNEvZP6WLTfvS3XIgR0T6RJVFt8FS6nclpk37iMjxwN+B11T1YlXdTQa2h4hc\nICJvxRTnAl+SYW0BXAr8m4hsA7bixvL+hBu3zai2EJERIjI5prgjLvONa1uke0YzF8j1ZklMxw1i\n9QDeSGmtkkhVd4jIHNyMq+uAY4DvAud7l/wRuFlE5gJ7gcnAH1JS2Tjypq6+Dtyrqvf65RnaHh8D\nJ4jI93CDuufh/r0n4X7RZkxbqOqR0c9FZBUwQVVfFZFhZFBb4CaE3C4iy4DncdnNZbiMpitxbIu0\nzmi8VO58XCOVAROAi1Q1tv8xHUUfNHQdbrrzWuBZ4CZvJgm4/tk5wAfAYlyf7QNJrGeiXAsUAneI\nyI6or1+QYe3hzQz6V9y05m3AFOAb3my0jGqLZmRUW6jqctzU5TtwQech4GpVXUSc28IOPjPGGJNQ\naZ3RGGOMST0LNMYYYxLKAo0xxpiEskBjjDEmoSzQGGOMSSgLNMYYYxLKAo0xxpiEskBjTJKIyCEi\nclmq6wEgIqNF5NhU18NkBgs0xiTPPcDFqa6EZy5waKorYTKDBRpjkier+UuSqq3Vx6Qp24LGGEBE\nFgFPqepvvOf/C1yoqkXe86OBf+D2T7sdt0FrLvA34AZVXeVddwnuhNN+uH2i7lXV/xGRKbg9pQBW\nq2r/AHXqhTt+/BzcXlRzgFtUtcJ7/QLgTtzRBzVeXX6gqqXeuSsPAt8COuF2sZ6kqktEZDVuM02A\nKap6V4sbzJgWsIzGGOc16p/NcQZQKCKDvOfnAO8BPwNGA98ERgLrcee75IlID9yW878DBgN3AdO8\nsZB7gaeBF4ATm6uMd1bK87gDqE4CLsGd7z7Te72f9/rjuC3bzwf6UxfMJnplX8ftvluG23UX3Dk0\n4Dabva+5uhhzoNL9mABjgnodmCgi7YHewMG4wHIa7gybc3DjGpOB01T1I4icPV+MyxyW4H6m1qlq\nCfAHEVkDlKpqhYjsAmpVtSxAfc7ABatTVXWv91nXAF+IyI1Ae+DHqvqod/0aEXkGFwTBZVRVQLGq\nbhaRCXhnM6nqFhEB2KaqsSfQGhN3FmiMcd4D9gEn485Mfxe3BfppIvJHXMC5GbfF/jsiEt3nnAcM\nVtXZIvIiMMc75+Rl4DFV3dqK+hwFHARs84KCr9b7rHkiUiUit+IyliOB43DHVQM8ijvka52IvIfr\ndpvZinoYc8Cs68wYQFX34DKWM4DTgfnAAlyAORXX9dTeu3w0MDTqawhuPARVvRg4HteldQqwUES+\n3ooqdcCdgDk05muwd89jgS+8z/o7MAnXDZbl1WMpcATuvJEvgFuAD0SkoBV1MeaAWEZjTJ3XcVlA\nH2AaoLhf1j8AXgVW4E4U7BnVddYBeAqY6h0P/H1VvQFYBNwlIq/hxldeof5hdM1ZAhwGlKvqFu+z\nBPgNMA64Hvi7qkbW5YjIT6IejwN2qOpsXIZ1B2486XjcQVXGJI0FGmPqvI7LTHYBH6vqPhH5DPg3\n4GJvnGUq8LCI7AFWAbfhsqCJQAFwvYhsB2bhZnYNA17y7r8DFy96q2ppM3V5C3c2+1MicgsuU3kU\nqFbVDSJSCnxLRP4F2OTV8UJccAToBvxCRDYDy4DLgUrqznbfCRwjIn9X1fJWtZYxAVnXmTEeVS0G\nVuIyhX1e8TtANfC29/wW3GyvJ3DTnfsB56jqRlVdCXwH+AZufGc28HtVneq99wmgLy7baa4utd59\nynHdeG/hJiV807vkQVzX3mvAQtwEhkuBQSKSj1sc+gdcwPvCe9+/Rk1E+C1uGvadwVrHmNazdTTG\nGGMSyrrOjEkyb41Mz2Yu2+JPazYm7CzQGJN8PYF1TbxeC/wL8EFyqmNMYlnXmTHGmISyyQDGGGMS\nygKNMcaYhLJAY4wxJqEs0BhjjEkoCzTGGGMSygKNMcaYhPp/GJoDnd+VzAcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fcc6390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, aspect=200.)\n",
    "\n",
    "lines = ax.plot(total_rain[x], total_rain, color='k', lw=2)\n",
    "ax.set_xlabel(x)\n",
    "ax.set_ylabel(\"Total rain ({})\".format(total_rain.units))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "That's the basic example. Suppose, now, we wanted to visualize contours of liquid water content in a 2D slice of the storm? That's easy to do using the aggregation features from `xray` and the same basic plotting interface. However, we'll also use the \"selection\" feature to choose a slice through the model domain; see [here](http://xray.readthedocs.org/en/stable/indexing.html) for more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "liq_water = dataset['QRAIN'] + dataset['QCLOUD']\n",
    "liq_water.attrs['long_name'] = \"Total water mixing ratio\"\n",
    "liq_water.attrs['units'] = \"kg/kg\"\n",
    "\n",
    "liq_water_slice = (liq_water.load()\n",
    "                       .isel(south_north=slice(30, 70),\n",
    "                             bottom_top=slice(0, 40),\n",
    "                             west_east=slice(100, 320))\n",
    "                       .mean(y, keep_attrs=True)\n",
    "                       .squeeze())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVIAAAD1CAYAAAACu30NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX3wHWWV5z8XkLyAGqL8CBDezOjB0VGCEXdGYwG6FoiI\nNYUTKV8QXFxXZaxxmRpwViClgkMUd8FlQYtYK6MFM2uxvKnjKuxMtEolE40OxBM2wTjJAAn8CC9G\nf0C4+0d33/Tt2/fefu+nu8+n6lby69v99Hm6b3/7POd5Ob1+v49hGIaRnQPqNiCMiLwEuBC4VlUf\nq9uequly/btcd7D6N73++9VtQISXAJf5/3aRLte/y3UHq3+j6++akBqGYTQOE1LDMIycmJAahmHk\nxKnOJuBA/9+jROSZWi2phyODfztY/y7XHaz+Zdd/f+Bw4F5VnSu6cNeE9Gj/3+/XakX93F23ATXS\n5bqD1b/s+q8Eflh0oa4JqQJ8Y+11LDlspm5bjJbT372N3qJjku+/6U4Aeq98x+B4fnQjc1vvG+wz\n72WvgpedNNhn6Pi7Lh3sO++t58Hhf0Rv0TH0N93J3Pe/BsBjW3fwkpcdyfwLvzVy7t4r3zGwITjW\nSMbDj+zkved/FOChMsp3TUj3Aiw5bIalRxxety1Gi+nPbqH/3DZ68w+lt3hZsv1nN4KsZD//t9mf\nv4ffP3IfHLRvv3mL5sHsRvY74oLR4xfNY87fd97sRjj0Rex3xL/j+R9tZO4g2KnbWHIgzBwE8x79\nJ/Z7zap9xx76Inrz99A/9EWg6+gtP5XeYntGMrC3jEKts8noJP3t67Ptr+voz24Zu9+cbhh7fPi7\n4P9BWTt12+C7nboNdN3Iufvb1w+2p7XfKBfXPFLDyE1/dsuI0PSWrhh4nmEhDPZL4pWGj+ktXjZ0\njkAIZ+QY5nQD835xC72lK4YEONhvRo4ZKitOfMNlhI8P8/wvbhmqX9p6GMVhQmo0kv7slljRGBFR\nXQeycqxgDoQKRjzNSaIU3jfqTc7IMZ7nGpwfTxiD/cL7TCQoI8KcbmBeeMOE+hnVYEJqNJaomA6J\naFikQmI6TmjimuvBtpFjxghcwE7dxgwwT1YmqcaQEI+U4f891muVEwf1C2w2Ma0ei5EajSOu6T5O\nRAcCFMQWY44dKmf7+sEnCXEiOEDXMacbhrzRKONiqtMIjhuqn66z2GlNmJAajSWpaGQVq7yiNOm8\nO3VbYrum7ReI9ZxumB4uMErBhNRoHGGBm9SDPu3Y3PtPEK2JnuoENjy0b3ROGrENE+6EMqphaoxU\nRP4MWA0sBbYBf62qt4nICuAnwJ7Q7p9T1c+XYqnRSsJCmHQ858RtkzyyUCwxNyV4fmERTUpUsIO4\naj80SsEon4lCKiKvANYCb1XVH4vIW4C7RORIYDlwl6q+swI7jRYSFcUsHSX97euHet6bRHQolNFc\nJjbtVXUzMOOL6AHAEuBJ4Bk8Id1YvomGUR5DYYIxzfi04YO8RJv3aQjCAf3t63n+F7cktt2buVVt\nPdvE1Bipqu4RkeOA3wNfx2vaP4UnpG8Uka0isk1E1ojIgRMLM4wKydrJNKm8ostMyrghUmNJMQMq\nGMkQiK+RnqSdTb/BC728FbhaRE4BdgK3A68CTgZOwYulGkYiRgbH1xjTa8KwoaTe6WC/0JCvSTSh\n7q6TaEC+qgZtjXtE5FvAu1T1rNAuD4rIFcAVwCVJyvSTXUXzsxwZt6/RHMYOYp9AHgGtMz5apHea\npaNpEuHB+n3GX+MONuePi2k5P5Y34d5Ej1RE3i4i/yeyeR7QE5EvisjBoe0LgN+lOPeFeMvmhT9d\nX4ux0WR5KNOKaG/xssLFsyyPLOsQqLLKmYg/oH/aoiwt4G5GdefCvIVO80j/GVghIu8DvgmcBpwO\n/DFes74vIhcDxwKfAm5Ice5r/TLDHImJaSsoa6pi1od84KFlpakD3Ysa7tUeTgV2RLblTv88UUhV\n9RERORP4EvDf8dT7LFXdJCJnAF8GHsUbS3q9ql6T9MS+Kz1UgY6mWGgE08Z7ttmLiXqsaUV5kkeZ\npVe+yCFTwdJ8uV80zeFBVf110YVOjZGq6g+B18ds3wy8rWiDjPZgC2hU1CyfQtJQSIfEtHBsiqiR\nmrp7ga2XebpA5xFEu77pMSE1ppKm2T60mrvRSJo6U6xOTEiNxKRZXm5wTEtjp3UNzM+CvdTKx4TU\nSER0KmUSgazqAe6aULgQdzWGMSE1SmeS6KYKG0wqJ6mYNnUYk+E0JqRGI8gSIrAeaKMqOimkbY3b\nVUmRzfssnVlNZcNDewufDloU9uLJTueE1EQ0Gy725Ka1aaSDKEUzv0mdS0b1dE5IA0xQq6XR19vi\nqsYUOiWk4XUXjfyEm/dVzGBq030rs3k/mEJq8+wrI0/OpkPw0pCcAjwBrFbVtWUaWyTBQ9n1KYwu\nUsnU0sDLzJHHqU3DkAbxURPfTExbRi/I2XSeqr4Q+ARwi7+W6Ffx0o7MAGcDV4nIG0q2Nxdt8mhc\nYmhBk5i4ZZbr3uhQgNE58uRsOgu4TFWfUdV78ZbE+0DZBheJCWt+6hx0PxKmcTyW6WpvfRRrpaUn\nU84m4A+AZyPLUW0Gji/DSMMN6uq5L9w7dVxwjeaRKWcT8A5GV8PfAywszrQSCVYDNwqjKLGzJr3R\nRLLmbFoBzI/sthB4OumJq87ZNHhATUALJ8gt31u8bOp00LiEdyP57Rucq95wHqdyNm0BDhSRo8K7\nA/elOLflbGopRYigxa+NknAqZ9NJwNHAlSJyAfBq4Bz/u6RUmrMpSKkwst1WcU9Fb+mKwkWuyutf\n1irwGx7ay4mH7194uUbhOJWzabMvoNcD2/Ga9Bf5vfeJcCJnk64Da0IWRhHxzaHl+ixeWjq9pSvo\ndyvc5VzOpseBVUUbVDX97evNIy2AIuKaTgqnnxjOMCbRiSmiTj6gDSHrtbPOoppJM0PJZjPlphNC\nahRLq0TSRMQogETDnwwjK3lEd7AeQpuE22gl5pEaY0nbrI/25sfOu88QKrChUPE0ZcppFzAh7VaP\nZfMZd7/sPho10gkhtV759ASeY7AwSHSBkFqa20WKpazsRnx0zNhp64Atlk4IqWEY+wheiCamxWFC\niv2gooS90TbSmXGhY7xRaO+9rQsTUiMVRTyAuV9cFg9Nx5jrZWJaHJ0U0s54JDlI8pBVGictWjwz\nlOdaahGb2+8OSXI2vQn4It7qTo8CV6nqV0RkBfATvHVIAz6nqp8vxdKCCES0rMUrjIJJklOpQg91\np27bl1yuAQS/82BxnokvyC50vpXERCH1E9zdDnxUVW8WkeXA90VkC/Ay4C5VfWcFdhoOUdl6oeEE\ndUY67JpVyrSm/dHAHap6M4Cq/gy4B/gT4ARgY7nmGa6TVFCbOgQtaMEkadZXPUB+2vnCIawm5bZq\nItOS321U1XODv30PdSWegC4H3igiW0Vkm4isiVl52ugQjY6Z1nCOIoQ3cxkmpoWSuLNJRF4M3AGs\n9//dhdfsfxVwMl5++9XFm2i4Tm/pikQi2lSvtOlM6ly1e1IMiRYt8bOI3gk8AKxS1T5eOuaAB0Xk\nCuAK4JKEZVaas8lIx7SV8CfFSaflbhp7zmmdIS3FiTnz3eloqj5nE4CInAj8GPiOqr5LVedEZJGI\nXC0iB4d2XcBoZtFJWM4mRynTS+ktXjaxfFvpKT1xQmxD/MZSfc4mETkM+C6wRlXXhL56EjgTeF5E\nLgaOBT4F3JDi3JXmbAL7cVVJVq90CFk5FMsL7t/QsLWOxvqyeLH92S2dvV4hqs/ZBHwIeClwqYhc\nGtr+X4EzgC/jjS3dA1yvqtckPbETOZuMwinLm836EpwnJ+Z+gZY5EL+MZv24sa5dDJvEUH3OJlUN\n4p7jeFux5hiuEAhiohlORXifrtEwzy1NFlObjFI8nZwiahTD0LJ6MbHP0nuEx4hdYg+0YWKZm67V\nt0JMSI3GMUkow1OAs5ZRN0XMoR8XjnC53k3GhBQbS1cFrWv6O4gTw6g6SmeF1GJE0+ktXlbLcCRX\nhkC56r0lFcwknWT2giuGzgrpgO4MRK6FtN7+yP5l35+WxQ3NK60HE1IjF/3t6xN7NVV7P656lGUz\nTky7ej2qoBNC2ubmS9l1q6t5XwVZhGXcWqRNXWS5v31967zyOuiEkAKt/LG0+QXRxJWkwmLqyuLP\niXvvW/h8VEl3hLSlNGG2io2KcAfX0qW0hW4LacM7mqoS0SKFMIkXHdfhZKMsDJfJk7PpEGAt3jqk\nTwCrVXVtmcYa+xgIkq6jv3RF6V7ftGX1mkqV0yVPPHx/61VvKRM90lDOpi+p6iLg3cCVIvIW4Kt4\nq0DNAGcDV4nIG0q2t3Aa3ZHix7XaKHDQ8HvjGhYDLZU8OZvOAi5T1WdU9V68JfE+UKaxWWmr0FRJ\nXo83fHym5n0DydLh1NTe/66TNWdTD3g2shzVZuD4Mow03GCch5hmLOk0ouWYV1oRDe8vqJusOZvu\nYXQ1/D3AwuJMy0d/dkurhwd1ycvuLV0x/UEf832Rg9BdGNJUSozVRDQ3mXI24SW8mx/ZbSHwdNIT\n15mzyXqAs5FmjdJxxwcvt/7sltTN9yIWaY7D63Bql5hE1ye1WU0D3MnZhCeoB4rIUeFdgftSnLvU\nnE397esTL0ps1EPQagiLq2GUTCk5m6b12gc5m76gqhcF21X1KeA2vB78BSLyeuAc4Bspzn0tnviG\nP6emM386Yx/OlnkgVRE3ZbSIOOY0EQ3OYa2JeDbO9sZ+Z4PwhziVUd25Nm+heXI2XQBcD2zHa9Jf\n5PfeJ6KqnE2TvNI2dWRkaSq7Si110XUTX64uxEfHEYjoxtker13cr9ka53EyZ9OqYs0xMqHroOKX\nQnSAflLhS5PfaWRfP6to1V5pVSKaZehT1BM1Ma2Hbk8RNTJRlbcYiGhdnrZrYYQ8PfbWvC8XE1Ij\nM72lK4qJj8Z0DI6MJ50mpkXP3Kkphl7lFNKx44Kt0y81nRbSpsYU+7Nbap/yV9S1CwvotFEWicaT\nNphARJOK6aQOpqS0qZ+gTjolpK411brKJBFONUZ1gqh2uSmbRGDtWSiWTgkp+D8gWdnoN7Frs5qq\n8uwH55GVTnimVfbkR73UIrzRgCY/C67QOSFtE3O6odYZK73Fy5wMj7TBG7Xl9ppFa4W0a29ZFwWt\nDCbdV9emQUYFvW5xjNoT+5ux5fYy0VohbStBR9OQaDjQzE1L8BCnfeFV+sKoeWX+cZ1PRTbrhzAR\nzUxnhbSpHpxr8dEySNR7n4A0TXzXvNmqaOpz4BqdFVKjuTTh4S8qTluZN2rkIs16pCeJyI7Q3ytE\nZK+IPBX6XFyOmdkoasC4UQ55BHFwX0sKa1Txu9k42zNhbAlJkt/1gPOAq4HwoiLLgbtU9Z0l2ZaL\nNHO6jXopMrFelT32ec4VFtCy58dH1yY1iieJR/op4M+Bz+KlGAlYDmwsw6giMBFtBlm90rjl/OJI\nLXYJO1yKEtFx3xftqU4aMdCEUInrJBHSG1X1BLwUI2GWA28Uka0isk1E1sSsPO0k1tx3n6Gpo2Ne\nioV3vDnQa11FU78N42xdY6qQqurDY77aiZeq+VXAyXj57VcXZllBdKGXuw1Me7lVvZp+1b+bOAG1\n+GlzSJSzKQ5VPSv054MiEqxdekmS42vL2dTAMZdtJ1c8O4EXuVO3FTqdM4tHFzStTRxrp5ScTZmE\n1E/L/GngUlUNEt4tYDSz6CQuBC7Lcn6j5fir1Q95hSHB7PsLPEPy8Z/jxDS8PShrXla7czJOZIPO\nqKpEuL99fZvjpnF54VYDl+cpNKtH+gRwJvC8P+TpWLxOqRtSlHEt8M3ItiNJkQBvUkqK/vb1Fgut\nkf7slth7kPgBneRpRmd2kcxLHLdPsD0qqGnLD7zOInvIw2lEjEI4FdgR2ZbLG4X0QtoHUNXnReQM\n4MvAo3g57a9X1WuSFpQ3Z1OSpqDFR+shEFEYvge9pStiX37h/aOe5qQUzEV3mkzyWo3WUH3OpjCq\n+n+BmdDfm4G3FW1QGlI3QSw+WitxHmrcCzEsnFlFNKt3GC53Ro7JHA91YdzmJDvmdAMLJhzbpmSK\nVZC5s6luBh7PGA9nBBPR0hn0qEe8ywF+3DO2pTAl5plG0MJjJvOIWts8UVcEvo00VkgnzYaxJr0D\nxMU4g23BSy2yTyCiebzCuIHnXRaQtEv3zekG5pnTkZrGCuk4ksRO29RkcSVlxJA3GuNdDtkZEtBp\nve5Jm/DTKKMjqG1E74U175PTutWfRobMODBbpe2MNOkZfSjjVvPP2jsekGWh5GizP/wpi7oXdM6C\nterS0XghDb8xx3qjITG1IVHlkPTBC8QzLKI7ddvgk5Q84lSUcEZF2egujRbSsaIY08kRpnULmjTM\n646KqCvk9XBdIu+4U3M40tFYIY2L3YTjc0Cneuqb+HLIKqJNaIa7KrBpsPhochorpJDwRofEtM1v\n2bpiWkOD6WFo1lHa5vo0yo5lhs+T5vxtEE0jH40W0ok0rLmbBVd67LPgUpM+jqg4xnVSNZ3walpG\nPlojpF38QQxijR0IYdQpXG0QzTjGTo7AmvVpSTyOVEROAm5V1SP9vw8B1uKtQ/oEsFpV15ZipeEk\nI836lpFXQF0T4JExtEGr7TWrqjemZUz1SEWkJyLnA98DXhD66qvAk3jz788GrhKRN5RiZQLSpvDt\nogdbFGObhGPCKa4345vK5b96tpBynv/FLYWU02Uy5WwSkYOBs4DLVPUZVb0Xb0m8D5RlaCI6EBc1\nDMM9suZsejnwbGQ5qs3A8QXaZqSgzSMSXGsi56HMbKFRD3XcWNIZOSZRR6V1RiUna86mgxhdDX8P\nsLAIo8qgzbE8wygDe2aSk3XRkj3A/Mi2hcDTMfvGUnjOJmvWO0eReZKMeC7/1bNcfvwLpu+YkRam\nHSklZ1PW4U8PAAeKyFGhbQLcl6KMCwGNfBKnGZlKi4cEJc1TVDktvuYuc/mvni2k4ync3G9xk/5u\nRnXnwryFZvJIVfUpEbkNuFJELgBeDZwDnJ6imNw5m4x9uOQ1TEoPYlBpIru0NHmSR0LcydnkcwFw\nPbAdr0l/kd97n4i8OZuy0IqEeI6FMCYtsG3UT5CBdCqO/a5KxLmcTY8DNpLXKJ0TD9+/VT33ThBK\naW3kpzVTRA1HsDip0UFMSBtMXfEsl+KxTSXa3C5ifGmq3vuwNxrjmVq4Jh0mpIZRM4GIljlYPy8t\n7sUvBCeFtL/b5ma7zjSvtAO9v7nJK5xZxo/a2N5yaHwW0axvyiY2T13wCsIrPgUjIKznPjsue6FG\ncpz0SNNiD3E1RJfNs+tePK4J6+Bl2UDHo0qcFNLeomTNj/CDHZfut22MCJesrLSXfHD+aWmuE9hk\nTcxiKHN6qJEcJ4U0NR0fD1fpJIMpvb1GfrJ4pSao9dIOITVKZdJyata8r580ImqLbJdDo4W0cw9x\nyAOsuld8JNV1Aoq0cSRNRstxLVZqTCZ3r72IXARcAcyFNp+mqj/KW3ZSksZGmzzPvs4e+yQiaj33\nRpcpYvjTCcDFqnp1AWXFEohIbM9hhjhdnT2QYUHMakeV3mhUwOd0w/D5dR19mv2ScpWoV+rqilFG\nMU375cDGAsqZSvihzuKhNdljcsH2wPNv++iIPJQZgrDmvrvkElIRWYi3oPMnROQhEblfRM4rxjSP\nTE1ahxfOyOKFTrsGZYrsuLJNTEcJRNTEtHvk9UhngHXAdcBRwIeBq0XktLyGhelvXz/4FFJezTOE\neouXFR5eCHrWi6zboKwsY0ZbOpY0qUieePj+Q58iKUJM7UVYLLlipP4CqaeENv1QRG4C3gV8d9Kx\nSXM2RcWzP7tlrAglyozYwMWdXWjWG/k8zaLXVE2yyr7LK/HXiFM5mwAQkdeJyCWRzQsYzTAax9Sc\nTXHelYmKUTdxglrH8Kwknmki73VC66GFU0PdydkU4kng0yKyGbgVzztdBbw5wbGZczZNbL46HB81\n2kMW4axypf/Ezf/o89L+58eJnE1DqOoDInI2cCXwP4HfAOeq6s8THJs5Z1MTm+d5qGuMZuCN9GF8\nnFTXTXz4ikyEZylHRsnafE8SBmvpM1ZvzqZxqOq3gW8XYMsQaTtN5smJXXibOouFXJJhL4N20tgp\nouEHN80A9aa+ZZPa3bSYVhN77qugjN7+EczxKIzGCmkWmry2Yn/7+rE//Ka+HIxiKCP/k5EOZ4U0\nELvEImFv11JeEOEyw56/pRLJThpPM+m+ucUzZm3bJjocdeF0qpHe4mXTY6Whmz+pU6ap3mh0VfoB\nNf3ouyqgda4+FZx7Wmw1TkyDmGxQRtJQirVy0uGsRxomyU0N5w8au0/DRDQJvaUrSq3X4EVmHn/t\npJ0pNW7KqnXMFo/THins80ptmTafiAdexcthcO1Levhm5JjOLTict/feev/dohke6RSxqEpQ6qC3\neFnjm1m2wHM5FHUtor+vNj9PZdEIIYV8gtJ0IXIBu4bNZ0aOmdqst/ucDeeb9mHCb8lwbvUuvj3r\nqHf4nEHs1MIt9ZKkiR/XwRTnhRrZaYxHGhAeFtWlmx+ua5X1HlzvyNJ/wd9jbSkxw6g178uli45J\nXnILqYgsF5GfisjTIvIzEXlDEYYlpQs3ve6XRxXXuIsznFx6IXThOSqTvMvozQfuAG4EXgxcA9wu\nIgcVYNtYBt5Qipvf9B9K0+036qGLL6g6yOuRngLsVdUbVHWvqn4NeAR4e37TjChtF9M0D71L3lyj\nsPGjpZBXSI8H7o9sU3+7UQIuiWndKVsMj6QvlWAY2khHk0O/qaaSV0gPAvZEtu0BFuYs12gISXvt\nk44lddErbaL320Sbm0ze4U+/xUstEmYh8NS0A8fkbDoK4OFHdmYypr97l/efXU8Of3HALnq/b4+2\nB/Wsu0793btGr3XA7rmRTXO/TVjw0mN4bGt0EfN4ZveWn5Po4WfKE6UjXnIA/7JzdPhS2nNGr0Nw\n/LOhaz5v95x3vw7YNbRv8Dty5XdVBiFNWVZGzqa8QroJ+HhkmwDfSHDshcBlcV+89/yP5jTLcJM0\nore0NCtSs7uB5wyO3xzeeD+jkbjO8f2YbauBy/MUmldI7wbmicjHgRuA9+OlaP6HBMfG5Wz6A+Au\n4K1AFwNwx+Fd01OBB2u2pWq6XHew+pdd//3xnLzfANGURrXnbHpGRE4HrgeuAB4A3qmqU7OIjsnZ\nFPz3X8vIq+I6oSbHjq7Vv8t1B6t/RfUvzTkrImfTL4E3FmCLYRhGI2ncFFHDMAzXMCE1DMPIiWtC\n+hheD1ru4G9D6XL9u1x3sPo3uv69ft8yDhqGYeTBNY/UMAyjcZiQGoZh5MSE1DAMIycmpIZhGDkx\nITUMw8iJCalhGEZOTEgNwzByYkJqGIaRE2fy2ovIcryl+P4QbxWpj6jqT+q1qjxE5CK8FbPCKyCf\nhrdg5Fq8fFhPAKtVdW31FhaPiJwE3KqqR/p/H8KEuorIlcCH8H6nXwc+qarPV254QcTUfwXwE4az\nTHxOVT8vIj2830fj6y8ibwK+iLeM3aPAVar6lTbdfyc80rqykdbMCcDFqvrC0OdHwFeBJ/HWdT0b\nuKrqFNdFIyI9ETkf+B7wgtBXY+vqr3H7duCPgFfirTD2n6u0uygm1H85cFfkN/B5/7uP0YL6+2J5\nO/AlVV0EvBu4UkTeQovuvxNCSjezkS4HNoY3iMjBwFnAZar6jKrei7f49QdqsK9IPgX8OfBZ/GXy\nE9T1/XgP3yOq+ghwJfDBqg0viJH6+4z8BkK0pf5HA3eo6s0Aqvoz4B7gT2jR/XdFSDuVjVREFuI1\ncz4hIg+JyP0ich7wcuDZyMK2m2n+dbhRVU8AwpnyptVVGP5NbPa3NZG4+oMnpG8Uka0isk1E1oQW\nOG5F/VV1o6qeG/zte6gr8V4orbn/rghp17KRzgDrgOvwEv59GLgaOAOIZhdo/HVQ1YdjNh9EfF0X\nhL7fE/luv5jEZc4zpv4AO/Gava8CTsZrma32v2tN/QNE5MV4Ibz1eF5pa+6/K51NmbORNhH/LXxK\naNMPReQm4M3A/MjuC4GnKzKtSvYwua7hhyr47jlVjebbaSyqelbozwdF5Aq8DqZLaFn9ReQ44E68\njuRVeC+P1tx/VzzSTYy67VHXvjWIyOtE5JLI5gV4ibkOFJGjwrsD91VmXHU8QHxdg3u+ieGQRqt+\nDyJyiIhc7ceKAxawz0trTf1F5ETgx8B3VPVdqjpHy+6/Kx5pnmykTeRJ4NMishm4Fc87XYXnkS7C\n69W8AHg1cA5wel2GloWqPiUitzG+rn8L/KWI3A08h+el3VSLseXwBHAm8LyIXAwci9cpdYP/fSvq\nLyKHAd8F1qjqmmB72+6/Ex6p766fjnchH8Mb+pEoG2kTUdUH8IZ7XIonqtcC56rqz4EL8IbIbAf+\nF3CR36PZFsIriU+q63XAbcBP8TzydXhx5KbTB/DHQ54BvAZvbOU/Abeo6jX+fm2p/4eAlwKXishT\noc9naNH9txXyDcMwcuKER2oYhtFkTEgNwzByYkJqGIaRExNSwzCMnJiQGoZh5MSE1DAMIycmpIZh\nGDlxZWYTACJyALC0bjsMw2gt21X1uaILdUpI8UT0wbqNMAyjtRwH/LroQl0TUgC+sfY6lhw2U6sN\n/d3baj3/EA/9sm4L9rH1p3VbMGBuqztruTy2dUfdJgz4l5176zYBgE27e9N3qojfsT/f58jSyndS\nSJccNsPSIw6v1Yb+/OjyqPXRf+5FdZuwj9l5dVswYM6hRDQvOLDw1mJm/m1/N4T04J47QkrJM+Gt\ns8kwDCMnJqSGYRg5MSE1DMPIiQmpYRhGTkxIDcMwcmJCahiGkRMTUsMwjJyYkBqGYeTEhNQwDCMn\nJqSGYRg5cXKK6MOP7KzbBPq7d9Vtwj52PVm3BfvYPVe3BQPmflu3Bft47Bl3HqXZvW5MzXy674Yd\n4M21LxN37r7HIoD3nv/Ruu0wGoE7D6qt/tgYFpVRqGtCGrg7p1L/cnrHAXebLU7aYbaYLVltKaVJ\n5ZqQBssFgmdGAAAGj0lEQVTW7FDVX9dpiIgcaLa4aYfZYrbksKWUpbGss8kwDCMnJqSGYRg5MSE1\nDMPIiWtC+hiw2v+3bswWd+0As2UcZks8pdrS6/dLXoPfMAyj5bjmkRqGYTQOE1LDMIycmJAahmHk\nxITUMAwjJyakhmEYOTEhNQzDyIkJqWEYRk4KXbRERJYDNwB/CDwAfERVf5KzzGOBG4HXAw8Bn1TV\nu/zv5gHXAe8CngWuUdUrYso4CbhVVY/MaMM5wOeAGeAe4EOqutP/7iLgCoZXlTlNVX8UKeNNwBcB\nAR4FrlLVr9RkS6L9qrAlVFZp9yi0z2HAL4Hzgt9Q5PvS71EKW0q/R0lt8fcp8xn6M7zB8kuBbcBf\nq+ptY8qp4jlKbE9AYR6piMwH7sATvRcD1wC3i8hBCY8/WUTiltr6e+DHwCHAJ4BvishR/nefA44C\njgXeBPwHEXl3qMyeiJwPfA94QcZ6vQb4H8Aq4KXAw8DXQrucAFysqi8MfaLCdQhwO/AlVV0EvBu4\nUkTeUrUtKfcr3ZaK7lHAjcBiYGQWSoX3aKotPlXco6m2lH1/ROQVwFo8EX8h3jN+i4gsjimn9HuU\nxp4wRXqkpwB7VfUG/++vichfAG8XkR8A/w14G7AHuF5V/2ZagSLySuDVwJtUdS/wXRH5R+A9wBrg\nfcB7VPUp4CkR+TLwQTzxBfgU3sX+LPBXoXIXp7DnvcD/VtV7/WP/CtglIoeq6i5gOfE/zjBHA3eo\n6s0AqvozEbkH+GMR2YD30qnKFsbtV8N1gWruESLyEeBp4F/HlFHVPUpiC1Rzj5LYUvb92SwiM6q6\nR0QOAJYATwLPxJRTxT1KY8+AImOkxwP3R7apv/0mvHUAjwVOBt4nIh/0K3GniDyO580eLSKP+5/3\n+Mf+WlXnomWKyCI8tzx8zs3+MQE3quoJwPqIXWPtiUHC51DVWWDWM10W+t9/QkQeEpH7ReS8aAGq\nulFVzx0U6L1ZVwIbgb+t0pYp+1V6XXxKvUd+nV8BfBL4T2OOr+QeJbWlinuU1BYquD++aB0H/B74\nOl5T+uloIVXdo6T2hClSSA/CexOE2QMcAZyGF9v8napuA74AXOAb/Q5VPQR4B/AbVT3E/9w8pszf\nAQv874h8vwdYGPyhqg9HjRSRJZPsSVGvhXhCvg4vTnsU8GHgahE5bUxZiMiL8V4a6/1P1baM2++c\nGmwp/R75XsXXgY+r6uNjjo+ev5R7lMKW0u9RUlsqeIYCfgPMA96KV9dTxtnk21DWc5TJniKb9r/F\nE7gwC4H/h5dcZ4uIBNv3Y3QVlrgEPHvGlPk0+y7EAv/v8HeTOHqcPeLFXjexL1b0Hxm9wIPzqLfq\nd/gC/1BEbsLr/Ppu9MT+W+5OvI64VcBrq7Zlwn4X1HVdYijsHgGfBn6uqt8LfTc22VOZ9yipLVXc\no6S2jKHI+wOAH7oDuEdEvoX3W7kn7uQl36PU9kCxQroJ+Hhkm+C9NZ4DZlT1WRi8TQ6OKSMa7N4E\nHCsiB6pqEKMQ4Aeq+riI7MRryu8KfXffFDsfGmePqu6I2iVeYFpCf78ULzC/SUReB7xNVa8MHRIW\n9nA5JwLfAW5S1Yv8bZXbMmG/rXgddpVelzEUdV1+hdeRcriIrPK/fhFws4h8RlWvipRT5j1KbEsF\n9yjVdYmhyN/t24G/UNV/HzpkHhDrJVfwHKWyJ6CwZfTEy4myFfg83hCo9+MN33gZnpj+ArgET/n/\nDvg3Vf1AgnLXAz8ALsVLonUL8EpV3SEiX8ALyp+N1/v2D8Bfquq3ImWcDPy9qh7q//2DpPaIyGuB\nfwTOAP4ZuBZYoqpnisjL8eIz7wduxfMibgXerKo/D5URDC9Zo6prIuVXbcvY/fCGlVRmS6S8kynh\nHsXs+yDwMVX9dmR76fcohS2l36OktoS+P5lynqEleA7ThcA38Zro3wROUtXNkXKqeI4S2xOmsBip\n7zGeDpyD12z/GPBOVd3jbzsM+DVeh9AOIGnO5T/Fc90fAa7G66Xf4X/3X/zyfoUXU/pKVERDhN8Y\nie1R1Y14Taq1vg1LgPP87x7AE/FL8Xr2rgXOjRGLD+EJ/aUi8lTo85mqbZmyX9XXJUrh9ygFpd+j\npFRxjzJSxjP0MHAm3jCjx4HLgbPGiFYVz1EaewbYws6GYRg5sSmihmEYOTEhNQzDyIkJqWEYRk5M\nSA3DMHJiQmoYhpETE1LDMIycmJAahmHkxITUMAwjJyakhmEYOfn/LLqxcmqffYgAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1157ea5d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.colors import BoundaryNorm\n",
    "\n",
    "levels = [0, 1e-5, 2e-5, 5e-5, 1e-4, 2e-4, 5e-4, 1e-3, 2e-3, ]\n",
    "cmap = plt.get_cmap(\"Oranges\")\n",
    "norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n",
    "\n",
    "fig = plt.figure(figsize=(5., 4.))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "cf = ax.contourf(liq_water_slice, \n",
    "                 levels=levels, cmap=cmap, norm=norm)\n",
    "cb = plt.colorbar(cf, ax=ax, orientation='horizontal', format=\"%1.0e\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pres = (dataset['P'] + dataset['PB']) * 0.01\n",
    "pres.attrs['long_name'] = 'Total pressure'\n",
    "pres.attrs['units'] = 'mb'"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

Hi, how can you convert the z-coordinate using python. Because wrf uses sigma coordinate in z direction. So how do you interpolate the sigma level to pressure level?

Good