xarray_fast_binning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A fast weighted-binning algorithm in xarray\n",
    "\n",
    "Application: taking the zonal mean of data on unstructured grids\n",
    "\n",
    "Additional package and sample data are available at https://github.com/JiaweiZhuang/cubedsphere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import cubedsphere as cs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read Cubed-Sphere data from the GFDL-FV3 model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'temp' (tile: 6, pfull: 20, y: 48, x: 48)>\n",
       "dask.array<getitem..., shape=(6, 20, 48, 48), dtype=float64, chunksize=(1, 20, 48, 48)>\n",
       "Coordinates:\n",
       "  * x        (x) float64 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 ...\n",
       "  * y        (y) float64 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 ...\n",
       "  * pfull    (pfull) float64 7.673 47.07 100.7 152.7 204.4 255.8 307.2 358.6 ...\n",
       "    time     float64 10.0\n",
       "    lon      (tile, y, x) float32 305.783 307.37 308.986 310.631 312.306 ...\n",
       "    lat      (tile, y, x) float32 -34.8911 -35.5988 -36.2858 -36.951 ...\n",
       "    area     (tile, y, x) float32 2.35456e+10 2.39763e+10 2.43917e+10 ...\n",
       "Dimensions without coordinates: tile\n",
       "Attributes:\n",
       "    long_name:     temperature\n",
       "    units:         K\n",
       "    valid_range:   [ 100.  350.]\n",
       "    cell_methods:  time: point"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prefix = 'atmos_daily'\n",
    "maindir = \"sample_data/FV3diag_C48/\"\n",
    "ds = cs.open_FV3data(maindir, prefix)\n",
    "dr = ds['temp'].isel(time=0)\n",
    "dr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A slow binning method \n",
    "Suggested by http://xarray.pydata.org/en/stable/examples/monthly-means.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define global variables for convenience\n",
    "dlat=4\n",
    "lat_bins = np.arange(-90, 91, dlat)\n",
    "lat_center = np.arange(-89, 90, dlat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 7.53 s, sys: 206 ms, total: 7.74 s\n",
      "Wall time: 7.93 s\n"
     ]
    }
   ],
   "source": [
    "def zonal_mean_slow(dr):\n",
    "    \n",
    "    area = dr['area']\n",
    "    area_grouped = area.groupby_bins('lat', lat_bins, labels=lat_center)\n",
    "    \n",
    "    # Calculate regridding weights.\n",
    "    # This step takes 90% of the computing time because the group object \n",
    "    # \"area_grouped\" are converted back to a normal xarray object, calling\n",
    "    # tons of xr.concat()\n",
    "    weights = area_grouped/area_grouped.sum()\n",
    "    \n",
    "    zonal_mean = (dr * weights).\\\n",
    "                 groupby_bins('lat', lat_bins, labels=lat_center).\\\n",
    "                 sum(dim='stacked_tile_y_x')\n",
    "            \n",
    "    return zonal_mean\n",
    "\n",
    "%time dr_zm1 = zonal_mean_slow(dr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The unstructured grid data is binned by latitudes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray (pfull: 20, lat_bins: 45)>\n",
       "dask.array<getitem..., shape=(20, 45), dtype=float64, chunksize=(20, 45)>\n",
       "Coordinates:\n",
       "  * lat_bins  (lat_bins) int64 -89 -85 -81 -77 -73 -69 -65 -61 -57 -53 -49 ...\n",
       "  * pfull     (pfull) float64 7.673 47.07 100.7 152.7 204.4 255.8 307.2 ...\n",
       "    time      float64 10.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dr_zm1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the zonal profile of the temperature field. Looks physically meaningful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x112dcac50>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEaCAYAAAD0YyfJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+8XFV57/HPlxzyixASQEJJgID8MNCWAIJa7PVYlEJr\nwWovpdIWSr39gSK3tL0StCVYb1Fu1RfcK/bWSi9QLEZEE20LIRcPtyIQDAmGJEBaScBoDvKbGEJy\nTp77x95DhpO95sycmTmzZ/J9v17zyj7PrNmz5mTOM3vWXvtZigjMzKw37dXpDpiZWfs4yZuZ9TAn\neTOzHuYkb2bWw5zkzcx6mJO8mVkPc5K3cSPpUEkvSVKn+2K2p3CSt7aR9ISkX6r8HBFPRcT0KPHF\nGZI+JOlBSdsk3VBw/+mS1knaIun/Sjqsxr5mSvp63vYJSb/V3t6b7c5J3uz1NgF/BXxp5B2SDgC+\nBnwM2B9YAXylxr6uB7YBbwB+G/iCpHmt7rBZLU7y1haSbgIOA76ZD9H8maTDJe2UtFfe5tuS/krS\nvZJelrRY0v6S/lHSi5IeqD5SlvQmSUslPZsfTf/nVvc7Ir4REUuA5wrufh/wSETcHhHbgYXACZKO\nKXj9U/P2H4+IVyLiXmAx8Dut7rNZLU7y1hYR8bvAk8B78iGav6ncNaLpbwLnA4cARwHfJTuKngk8\nClwJryXNpcA/AgcC5wGfl/SmoueX9HlJz0t6rurfyvaqMb6s44GHq17jVuDf8/hIxwA7IuI/qmIP\nJ9qatY2TvLXbaCdZ/yEiNkTEy8C/Av8REd+OiJ3AV4ET83bvAZ6IiJsi8zBwO1B4NB8RH4qImRGx\nf9W/le35Y3wt04AXR8ReAvZNtH2pzrZmbdPX6Q7YHm+wavuVgp+n5duHA2+VVBlGETABuLntPdxl\nCzB9RGw/4OUm25q1jZO8tVMrZ9E8BQxExC/X01jSF8hOdo7sg4ANEfFzY+jDGuCCqufYB3hjHh/p\ncaBP0hurhmxOSLQ1axsP11g7bQaOHBEb6xz5bwHHSPptSX2S9pb05tSYfET8cUTsm58PqL7tWyvB\nS5ogaTLZt4Q+SZMkTcjv/jpwvKRflzSJ7HzBqoh4vOD5t5INJ31C0lRJbwd+jfH95mHmJG9t9Sng\nL/ITnpflseoj67qP9CNiC3AG2QnXH+W3TwETW9TXio8DW4GPkp0Q3ko2ZZKIeAZ4P/DXZLNv3pz3\nBwBJCyT9c9W+PgRMBZ4mO2H8RxGxrsX9NatJJb4uxczMmuQjeTOzHuYkb2bWw7oqyUs6U9Kjkh6X\n9NFO98fMrOy6Zkw+vxT+ceB0spNuDwLnRcSjHe2YmVmJddM8+VOB9RGxEUDSrcA5ZJe+v0ZSd3xq\nmVkpRMSYS1/PkGLkJdA1bIyIuWN9rrHqpiQ/m+yCmIofkiX+3fxdvL4G1ASGW9KBCQy1ZD/DDf7a\nh5lQM/6thSt5z8ITR20/lNxP/f1J7btV7dtp2cL7eNfCt3XkuRt9Dzbevvi92ZfYT2r/lfiShQ9z\n9sITWtqXRjX6d5LeT/o9+Adq7rKFF4FP1tn249lV2+Oum5J83ZYsfK2GFMf2z+K4/jd0sDdmVhaP\nDWzmsYHB0Rs2YO+W7q31uinJbyIrXVsxJ4/tpnIEsktrjuTNrLsd238wx/Yf/NrP37rq+03vs+xJ\ntOz9q/YgcJSkw4Efk11p6JV2gGOq3rSWdmT/nE53oWsc2z+r013oGj6Sb5GIGJb0YbKa4nsBX/Il\n4plj+n+m013oCkf2H9rpLnSNY33gULcpne7AKLomyQNExB3AsZ3uh5lZRdmTaNn7Z2ZWamUfruma\ni6HqJSmYOeI1lWcWXyZ1Hjg18ywV35aIt/W/dEci3pppc3ueRo+zWpRSUjPDJxfEUl1Mxbvl7w3g\neTU1T15S3Fpn2/Nobk7+WPlI3sysCWU/kneSNzNrgpO8mVkPK3sSLXv/zMxKzVMozcx6mIdrzMx6\nWNmTaFctGmJmVjZ713krImmSpAckrZS0WtKVeXympKWSHpN0p6T9qh6zQNJ6SesknTFa/5zkzcya\n0FfnrUhEvAq8MyJOBOYDZ0k6FbgcWBYRxwJ3AwsAJB0HnAvMA84CrpdUc+69k7yZWROaOZIHiIit\n+eYkss+DIFsQ6cY8fiPw3nz7bODWiBiKiA3AehLralSUfThpbJ6/vU07Tl3tmZK6CrTRq0Y71b6R\nfVhrtepK2MR+ItH+laL2De674fYpZT+lmWk2ieZLm64A3gh8PiIelDQrIgYBImKzpIPy5rOB+6oe\nvimPta1/ZmZ7tKlNPj4idgInSpoOfF3S8exenGTMxUqc5M3MmpBKovcC321gPxHxkqQB4ExgsHI0\nL+lg4Om82SagumZ2cvGkCo/Jm5k1Ye++4lt/H1xRdSsi6cDKzBlJU4B3A+uAJcCFebMLgMX59hLg\nPEkTJR0BHAUsr9U/H8mbmTWhr94sWnw662eAG/Nx+b2Ar0TEv0i6H1gk6SJgI9mMGiJiraRFwFqy\nk2kXxyilhHuz1DBfa9PefeLVJ17HS5tPvDbUvpdPvL6/6VLDW/epr+3Un7rUsJlZ16n7SL5DSt49\nM7Ny23tSp3tQm5O8mVkzSp5FS949M7OSK3kWLXn3zMxKruRZtOTdMzMrubItXD5CbyZ5va++dl0z\ne7TdUygbed5Gp5Ha2KSmDzY6rbCdUyhLZiyTE1uRA0qeRUvePTOzkit5Fi1598zMSs5TKM3MeljJ\ns2jJu2dmVnI+8Wpm1sNKnkVL3j0zs5IreRYteffMzEqu5FnUi4aYmTVjQp23ApLmSLpb0hpJqyV9\nJI+fIOm7kh6WtFjStKrHLJC0XtI6SWeM1r2SfwaZmZXc5KYePQRcFhGr8kT+PUl3AV/M49+RdCHw\n34C/lHQc2QIi88iW/lsm6ehaC4f4SN7MrBlNHMlHxOaIWJVvbwEeBWYDR0fEd/Jmy4D359tnA7dG\nxFBEbADWA6fW6l5vHsm/f/QmYzLcov00Wl1gKHFZ+Y7U5eZTisON9t+LQJVPo3+xjU7va6SCQaN9\naXffx+K2FuyjRVlU0lxgPnA/sEbS2RGxhOzIfU7ebDZwX9XDNuWxdnevMZLmADcBs4CdwBcj4jpJ\nM4GvAIcDG4BzI+LF/DELgIvIUs+lEbG0E303M3udRBYd2AwDg/XtIh+quY0st22R9PvAdZL+gmzx\n7u0t7l7bjRyHWiFpKfB7wLKIuEbSR4EFwOVjGYcyMxsXiW8c/bOzW8VV3y9uJ6mPLMHfHBGLASLi\nMeCX8/uPBn41b74JOLTq4XPyWFJHxuQLxqHWkXX2HODGvNmNwHvz7YbHoczMxkVfnbe0G4C1EXFt\nJSDpDfm/ewEfB/42v2sJcJ6kiZKOAI4Clo/WvY4aMQ41KyIGIfsgkHRQ3qzhcSgzs3HRRBaVdBpw\nPrBa0kqy4sdXAMdI+lD+8+0R8X8AImKtpEXAWrK63xePNqLR0SRfMA41srNjG45Zs3DX9hv64aD+\nMe3GzHrM0wPwk4HW7rOJKpQRcS/FAz53ANclHnM1cHW9z9GxJF80DgUMSpoVEYOSDgaezuONjUMd\nv7D1HTaz7ndQ/+sP+tZd1fw+Oz4eUlsn58nvNg5FNt50Yb59AbC4Kt7QOJSZ2bhofky+rTo1hTI1\nDvVpYJGki4CNZDNqxjQOZWY2LlxqeHc1xqEA3pV4TEPjUGZm46LkwzUl756ZWcmVPIuWvHtmZiVX\n8ixa8u6N0ds63YEW683/Jeukbq5L1Mq/h1bUrvFC3mZmPazkWbTk3TMzKznPrjEz62Elz6Il756Z\nWcmVPIuWvHtmZiXn4Rozsx7W3Bqvbec1Xs3MmtFE7RpJcyTdLWmNpNWSPpLHT5B0n6SVkpZLenPV\nYxZIWi9pnaQz6umemZmNVXPDNSNXyfuepLuAa4ArI2KppLOA/wG8cyyr5DnJm5k1o4ksGhGbgc35\n9hZJjwKHkK19vV/ebAa7Squ/tkoesEFSZZW8B9rQPTMza1UWrVol7wHgT4A7JX0GEPALebOGV8nr\nzSR/VJ3tWvXqG91Pqn2j8dTXxL3b/LyN6HNFaACG1KL9tCi+IxEfbuNzdqqUQruftwWzawpWyfvj\nfPsbkn6DbP2Nd49l372Z5M3Mxksiiw6sgoGHR394YpW8CyLiUoCIuE3S3+fxxlbJS3fPzMzqkihQ\n1v+W7FZx1U3JPRStkrdJ0jsi4h5JpwPr8/gS4BZJnyMbphl1lTwneTOzZjSRRWuskvdfgOskTQC2\nAX8AY1slz0nezKwZzc2uqbVK3puLgo2ukuckb2bWjJJn0ZJ3z8ys5Fy7xsysh5U8i5a8e2ZmJVfy\nLFry7pmZlZzXeDUz62Elz6Il794YHdjk49tdLmByYlrr5O2F4SnTthbGJ05+tTA+aWLxfiZSHJ9E\n8X6K2k9o8BrxvsLr5Bs3oUX7adRwi86qDTW4n+HEm2c7EwvjryYOJ1Ptt26dWtx+2+7tt29LHKoW\ntAXSJRy2FYdbUnqhk0qeRUvePTOzcgvPrjEz613DJc+iJe+emVm5OcmbmfWwoQn1rqK6s639SHGS\nNzNrwvZJ9c6hfKWt/Uhxkjcza0KrZmC1i5O8mVkTGp0eO97qHUwyM7MCw/TVdSsiaY6kuyWtkbRa\n0iV5/FZJD+W3JyQ9VPWYBZLWS1on6YzR+ucjeTOzJjQ5XDMEXBYRq/J1XldIuisizqs0kPQ3wAv5\n9jzgXGAe2dJ/yyQdXWvhECd5M7MmNJPkI2IzsDnf3iJpHdmyfo9WNTsX6M+3zwFujYghYIOk9cCp\nwAOp5+jocI2kvfKvI0vyn2dKWirpMUl3Stqvqm1DX1HMzMbDMBPquo1G0lxgPlUJW9IvApsj4gd5\naDbwVNXDNuWxpE4fyV9Ktlbh9Pzny4FlEXGNpI8CC4DLJR1HI19Rmj0PkqpRMzkVb00tmn2nv1wc\np7H4lMRUrVT7qSRq4xTUtJmUqH+Tqi2TjjdWA6dbpMZeU3/kqfiriZoz2xM1arZSXIvmZfYtjk+d\nlojv3v6VxL5T9W+2vjylML6zr7g92xK1blLZqdFaN22Wqhu0fOAVHhyob9pkPlRzG3BpRGypuuu3\ngH9qpn8dS/KS5gC/Avx34LI8fA7wjnz7RmCALPGfTYNfUczMxkPqg/rk/mmc3L/rw/QLVz1f2E5S\nH1mCvzkiFlfFJwDvA06qar4JOLTq5zl5LKmTwzWfA/6cbHXyilkRMQivjVUdlMcb/opiZjYeWjBc\ncwOwNiKuHRF/N7AuIn5UFVsCnCdpoqQjgKOA5bV23pEjeUm/CgzmZ5T7azRNnjE2MyuDZubJSzoN\nOB9YLWklWc67IiLuAH6TEUM1EbFW0iKyYe4dwMW1ZtZA54ZrTgPOlvQrwBRgX0k3A5slzYqIQUkH\nA0/n7Rv7ivLFhbu2T+qHk/tb2HUz61orBuChgZbuMnUeph4RcS+Js4gR8XuJ+NXA1fU+h0b5EGg7\nSe8A/jQizpZ0DfBsRHw6P/E6MyIqJ15vAd5CNkxzF1B44lVScH+Tr8knXgGfeB2LrjnxSuLEa0H7\nlp14/WmDJ15Tb5FWnnh9q4iIRAdGJym+EyfX1fbtWtHUc41Vp2fXjPQpYJGki4CNZDNqxvQVxcxs\nPKRW3yqLjif5iLgHuCfffg54V6JdQ19RzMzGQ9lr13Q8yZuZdbNmxuTHQ7l7Z2ZWci41bGbWw5zk\nO2H3SSHFGn31DbafOLm4I1NTs2sSs18O4NlE/JnC+MysYN1uZiTiqeedVhCfmpi5MynxS2901k23\na3wWTWq2TPEMldRMlxeY0VD8GQ4sjD9bMHvq5cSMqkRXGB4qfq2vpGbXpGbLpCoCpN46HZqw5TF5\nM7Me5jF5M7Me5imUZmY9zMM1ZmY9rKuHayTtX+v+/OIlM7M9VrfPrllBVhWtqN5CAEe2vEdmZl2k\nq5N8RBwxXh0xM+tGzST5fPGkm4BZwE7gixFxXX7fJcDFZJND/zkiLs/jC4CL8vilEbG01nOMNlxz\nUq37I+Kh+l6KmVlvavLE6xBwWb62xjRghaSlwMHArwE/FxFDkg4EkDSPRpZCZfThms/UuC+AX6r/\ntZiZ9Z5U+ed65Cvgbc63t0haR1ZO/Q+AT+VLnhIRlSsfz6HBpVBHG65555h7b2a2B2jVmLykucB8\nsoT9N8B/kvTXZNf+/llErCD7ALiv6mGjLoVa19wfSb9bFI+Im+p5/LgbebV/ozOcUouDpC6n7ite\nB2BoR/F//oS9inc0JbF4R6rswIGJcgcHMVgYn/XaQluvlyqbMIPdFx5O9SW1mMjE5CIjxdeg93VJ\nuYPUV/TUdLrUBTOpxUFSi328wMzC+LMcUBgffG2Z5PoUJaxUEkstVJI0lFgvI1W+4KeJ+LbU/hvr\nTqu0Yp58PlRzG9kY+5Z8ce+ZEfFWSacAX2WME13qTX+nVG1PBk4HHiI7YWBmtsdKfbA/MfAkGwae\nHPXxeUK/Dbg5Ihbn4aeA2wEi4kFJw5IOIDtyP6zq4bWXQqXOJB8Rl4zo1Azg1noea2bWy1LfdA7r\nP4LD+ndNULznqntTu7gBWBsR11bFvkF2zvMeSccAEyPiWUlLgFskfZZsmOYoYHmt/o31Uq2fAp5e\naWZ7vCanUJ4GnA+slrSSbELLFcA/ADdIWk1WV/d3YWxLodY7Jv/N/MkB9gKOAxY1/IrMzHpMM0k+\nIu6F5A5+J/GYhpZCHW2e/KSIeJXsTG/FELAxIn5Y75OYmfWq1JoAZTHakfx9wEnAByOi8FPFzGxP\n1tVlDYCJkj4A/IKk9428MyJub0+3zMy6Q7cn+T8iOykwA3hPVVxkY/RO8ma2R+vqevIR8R3gO5LW\nAhOBt5Ml938DvtD+7pmZlVtX15Ov8p+AF4Hr8p8/QHYh1Lnt6JSZWbfo9uGaiuMj4riqn7+dH92b\nme3ReiXJPyTprRFxP4CktwDfa1+3mrR5xM+pV5mKp2rXpGpmJOyc1tjUqlTdllS9mFTNmUP4cWH8\nUJ4qjM9K1Lop2v8MXihsOzVRd2fiq68Wxie9urMwnjKhQ3VJhhv8Jv7qpL0K49snFb8XUvVfXmBG\nYTxVoyb1HpmQeE+lKicW9SdVRydVf2h4KJH0Un8/Lybiu5dOyqRq2nToPZKqP1QW9b6FTwa+K6lS\niOEw4LH8aqyIiJ9vS+/MzEquV8bkz2xrL8zMulRPDNdExMZ2d8TMrBv1RJI3M7NiXT1P3szMauuV\nMXkzMyuQWvWrLIrne5mZWV2GmFDXrYikOZLulrRG0mpJl+TxKyX9UNJD+e3MqscskLRe0jpJZ4zW\nPx/Jm5k1ocnhmiHgsohYla/zukLSXfl9n42Iz1Y3ljSPrNLAPLKl/5ZJOrrWwiEdO5KXtJ+kr+af\nRmskvUXSTElLJT0m6U5J+1W1b+jTy8xsPAwzoa5bkYjYHBGr8u0twDqyZf0gKwQ50jnArRExFBEb\ngPXAqbX618nhmmuBf4mIecAJwKPA5cCyiDgWuBtYACDpOHZ9ep0FXC8psfS7mdn4aSbJV5M0F5gP\nPJCHPixplaS/rzrgnQ2vu3R9E7s+FAp1ZLhG0nTgFyPiQoCIGAJelHQO8I682Y3AAFniP5v80wvY\nIKny6fUARUauWdVoWYNpiXjxlebp/cwovuPV7cUnaiZMLL4EPVUyIHUpe6pMwSH8qDCeKndwyEtP\n7xbrK66kAM8l4qlL2YurHXTs0vSU5B9I4o6+ScXlGvaZ/EphfOb+xfFZBxT/ogenF5eVmMj2wnjq\npGCqbEJRfGLiPys1TLF9W6Kcx5biMM8k4iPLk4y2n9R7p83vqVQCf3ngIbYMPFTXPvKhmtuASyNi\ni6TrgU9EREj6JPAZ4INj6V+nxuSPAJ6R9A9kR/HfA/4rMCsiBiH7GiPpoLz9bLJVqipG/fQyMxsP\nqZOqU/pPYUr/Ka/9PHjVlwrbSeojS/A3R8RigIj4SVWTLwLfzLc3AYdW3TcnjyV1arimj2xZwc9H\nxElkJYcuZ9di4RU1VyE3M+u07Uyq61bDDcDaiLi2EpB0cNX97wMeybeXAOdJmijpCOAoYHmtnXfq\nSP6HwFMRUalk+TWyJD8oaVZEDOYvsjJe0Nin178t3LV9WD+8sb9F3TazrrZxAJ4caOkumylrIOk0\nstX3VktaSXZgewXwAUnzgZ3ABuAPASJiraRFwFpgB3BxrZk10KEknyfxpyQdExGPA6cDa/LbhcCn\ngQuAxflDlgC3SPoc2TBN7U+vX1zYtr6bWRc7vD+7VXznqqZ32UxZg4i4Fwp3cEeNx1wNXF3vc3Ry\nnvxHyBL33sAPgN8je7GLJF0EbCRfeWosn15mZuPBZQ0SIuJh4JSCu96VaN/Qp5eZ2XhwFUozsx7m\nJG9m1sOGdzrJm5n1rFdTF3+VhJO8mVkTkguXl0RvJvl6yxqk/m+KF6dPX6KfLGtQHB4eSjwgUZY6\ndcl6qqzBARRfEp8qa3DYT3YvXwBkc55GSjRNPGV2mVuRRssaFFd8aL/UeyT1f546qNsnET8gsfuD\niuOzjyz+RQ+/obijqfIFMygujzCloIRGX6O//G2JN3LxU6bLF/w4EU/tJ9XNdpc1cJI3M+tdQzuc\n5M3MetbO4XKn0XL3zsys7DxcY2bWw5zkzcx62LZyr1/kJG9m1oySLXQzUieX/zMz635Ddd4KSJoj\n6e58nevVkj4y4v4/lbRT0v5VsYbWu/aRvJlZM3Y09egh4LKIWJUvAbhC0tKIeFTSHODdZBV5AZA0\nj13rXc8Blkk6ulZVXh/Jm5k1Y7jOW4GI2BwRq/LtLcA6di1t+jngz0c85Bzy9a4jYgNQWe86yUfy\nZmbNaNGYvKS5wHzgAUlnk62et1p63Yndhte7dpI3M2tGKsmvGoCHB+raRT5UcxtwKdlx/xVkQzVN\n680kP7J2TWoa696JeKLmTPI/c1oivqU4vD1V22NqcXhiotBLunbNM4XxQ7YmioEU1aiB7IvgSE8l\n2qZq17yYiDdauyal3TMbUn8hjdau2S8RT9SoSf7eEmbtU1xUaHBq8ROk3juTCuokTUiMNSTrqKem\nFKZqzhS/XdOrOCf+rnglEW933aNUTas39We3ipuKlxqU1EeW4G+OiMWSfhaYCzys7DB+DvCQpFPJ\nfiuHVT289nrX9GqSNzMbL80faNwArI2IawEi4hHg4Mqdkp4AToqI5yVV1rv+LPWsd42TvJlZc5pI\n8pJOA84HVktaCQRwRURUL+QdgGBs6107yZuZNaOJKZQRcS/pAeVKmyNH/NzQetdO8mZmzejUWgd1\ncpI3M2tGycsaOMmbmTUjNbumJJzkzcya4SN5M7Me5iRvZtbDnOTNzHpYc1Uo2643k/zmOttNTsRT\nJ1JSv63U5dqJ/aRWd09dPl50qTmkL02fmejQpB8V9ydZquCJgliqBELxVfXpy/N/mog3Oh2tVUdR\njf4lpGY275OIp8oapMpBpMo+JN6zkxLlEWYeUfxeSL13ikpopN6XScW7Tv+d/KTBeGo/rSqV0ShP\noTQz62EerjEz62GeQmlm1sM8Jm9m1sM8Jm9m1sNKPibvNV7NzJoxVOetgKQ5ku6WtEbSakmX5PFP\nSHpY0kpJd0iqri+/QNJ6SesknTFa95zkzcyasaPOW7Eh4LKIOB54G/BhSW8CromIEyLiROCfgSsB\nJB0HnAvMA84CrteIRWBH6liSl/Qnkh6R9H1Jt0iaKGmmpKWSHpN0p6T9qto39OllZjYuhuu8FYiI\nzRGxKt/eAqwDZufbFfsAO/Pts4FbI2IoIjaQLdJ5aq3udSTJSzoEuIRsSaufJzs38FvA5cCyiDgW\nuBtYkLdv+NPLzGxcbKvzNgpJc4H5wAP5z5+U9CTwAeAv82azef3li5vyWFInh2smAPvki9hOIevs\nOcCN+f03Au/Ntxv+9DIzGxfNDdcAIGka2WLel1aO4iPi4xFxGHAL2UHxmHRkdk1E/EjSZ4Anga3A\n0ohYJmlWRAzmbTZLqlysPRu4r2oXtT+9Rl4OvXeiXaqsQcq0RDx1GXdi9fidw8W/9tTl40WXmkP6\n0vQZqeu+B4vDPJmIF5UwKCp1AJAqmfBccfiVxJHNK6lL00tmyqREPPWe2j8RT5V3SEmVRziiODzt\niOL3yFS2FsaLSmhMaHT6SGpKYaPlDlLxVCmI1NFyu+exp17vTwbgmYFRH54f6N4G3BwRiwuafJls\nXH4hWe47tOq+OXksqSNJXtIMsqP2w8mqm3xV0vlkC9ZWq7lArZlZx6U+A2f2Z7eKR69K7eEGYG1E\nXFsJSDoqIv49//G9wKP59hLgFkmfIzvQPQpYXqt7nZon/y7gBxHxHICkrwO/AAxWjubzKUOVsleN\nfXr9dOGu7b37s5uZ2fYB2DHQ2n02MU9e0mnA+cBqSSvJDmyvAD4o6Viy7wkbgT8CiIi1khYBa8m+\no1wcETUPhjuV5J8E3ippMlntuNOBB4EtwIXAp4ELgMpXl8Y+vfZZ2KZum1lXm9if3SpeSR5d16+J\n4aCIuJfimqZ31HjM1cDV9T5Hp8bkl0u6DVhJ9itaCfwdsC+wSNJFZJ9e5+btG/70MjMbFy5rUCwi\nrgJGfow+RzaUU9S+oU8vM7Nx4SqUZmY9zFUozcx6mIdrzMx6WMmrUDrJm5k1w0nezKyHeUy+A0Ze\nGt/omFnqbHnqE7vB/e81oXhHqcvHpybqIyTLGry0pTDOjxMdSpUkeKog1kdWOajAcwX7GUz8zl5K\nPGXqV9ypv6NURYy+REenJ8oUzHqxOL5/qoxDqjzCYYl4Yv9FZQogXUKjyImsYiXzd4vPZyWrOHH3\nBzT6d5L6HaTiqb/PRv9uW8Vj8tZTGkjw1huKEjxQnOD3RJ5CaWbWwzxcY2bWwzxcY2bWwzy7xsys\nhznJm5n1sJKPyXdy+T8zs+43VOetgKQ5ku6WtEbSakmX5PFrJK2TtErS1yRNr3rMAknr8/vPGK17\nTvJmZp0zBFwWEccDbwM+LOlNwFLg+IiYTzZxeQGApOPISrDPA84CrpekWk/gJG9m1iERsTkiVuXb\nW4B1wOyg520WAAAHlUlEQVSIWBYRO/Nm95OthgdwNnBrRAxFxAayD4BTaz2Hk7yZWQlImgvMBx4Y\ncddFwL/k27N5/bXom/JYkk+8mpk1JXXm9Z78NjpJ04DbgEvzI/pK/GPAjoj4p7H2rjeTfL2vKtWu\n0XiqzsiU4vDUfRurRZOKT0vE+54uDMNgIp4qSVAQT5UvWJ84sZTqSqp2TerPpVOz1FL/5amaNtMT\n8ZcSL+DoxO9z//0SO0rUqEnWeUkYLlxWtFjD9W/anVUa3X/bs1zq3Xlafqv4ZGErSX1kCf7miFhc\nFb8Q+BXgl6qabwIOrfp5Th5L8nCNmVlTdtR5S7oBWBsR11YCks4E/hw4OyKqP8KXAOdJmijpCOAo\nYHmtnffmkbyZ2bgp/mZeD0mnAecDqyWtBAL4GHAdMBG4K588c39EXBwRayUtAtaSfXJcHBFR6zmc\n5M3MmjL2q6Ei4l4oHDs7usZjrgaurvc5nOTNzJpS7roGTvJmZk0pd10DJ3kzs6b4SN7MrIf5SN7M\nrIf5SN7MrIeNfQrleHCSNzNriodrxt+MET+nruBOlB3gwDr3W5G6BD2xnxlTXyiMz6Q4PiMRT7Xn\n2UR/nkvEU7UHCuIbEt9MU9dVpyopFBdkKPsx0S6pt06j5RqmJH6f+7e5fMFQIl5UqmBi4kmnsrX4\nSVN/J9MS8VT71N9hyqREvNYarM83+ByFPFxjZtbDfCRvZtbDfCRvZtbDfCRvZtbDnOTNzHpYuacL\nuJ68mVlThuq87U7SHEl3S1ojabWkj+Tx35D0iKRhSSeNeMwCSeslrZN0xmi985G8mVlTmhquGQIu\ni4hV+RKAKyQtBVYDvw787+rGkuYB5wLzyFaFWibp6Fo15X0k3wMeGNjW6S50hSc63YEusmFgY6e7\n0EXGfiQfEZsjYlW+vQVYB8yOiMciYj2gEQ85B7g1IoYiYgOwHji1Vu+c5HvAAwMNXiGzh9rQ6Q50\nkY1O8g1oevk/ACTNBeYDD9RoNht4qurnTXksycM1ZmZNaX6efD5UcxtwaX5E3zJO8mZmTUkdpf+A\negYJJfWRJfibI2LxKM03AYdW/TyHdFWRbP+jrAHbdST11gsys7aKiJHj3nWTtAE4vM7mGyNibsE+\nbgKeiYjLCu77NvBnEbEi//k44BbgLWTDNHcBNU+89lySNzPrFpJOA/4f2WyayG9XAJOB/0lWpu0F\nYFVEnJU/ZgHw+2RfIS6NiKU1n8NJ3sysd3l2jZlZD3OS71KSrpT0Q0kP5bczq+5r6Iq4PYGkMyU9\nKulxSR/tdH/KRNIGSQ9LWilpeR6bKWmppMck3SkptWqClZyHa7qUpCuBlyPisyPi84AvA6eQXxHH\nKCdmep2kvYDHgdOBHwEPAudFxKMd7VhJSPoBcHJEPF8V+zTwbERck38ozoyIyzvWSRszH8l3t6JZ\nAQ1fEbcHOBVYHxEbI2IHcCvZ78kyYvdccA5wY759I/Dece2RtYyTfHf7sKRVkv6+6ut0w1fE7QFG\n/k5+iH8n1QK4S9KDkj6Yx2ZFxCBkl94DB3Wsd9YUXwxVYpLuAmZVh8j+ID8GXA98IiJC0ieBzwAf\n3H0vZqM6LSJ+LOkNwFJJj5G9z6rtscN93c5JvsQi4t11Nv0i8M18u+Er4vYAm4DDqn7276RKRPw4\n//cnkr5BNrw1KGlWRAxKOpj0cu9Wch6u6VL5H17F+4BH8u0lwHmSJko6AjgKWD7e/SuZB4GjJB0u\naSJwHtnvaY8naWpeNwVJ+wBnkF2YswS4MG92ATDa5fZWUj6S717XSJoP7CQrsPiHABGxVtIiYC3Z\nFXEX78kzawAiYljSh4GlZAc2X4qIdR3uVlnMAr6elwPpA26JiKWSvgcsknQRsJGshrl1IU+hNDPr\nYR6uMTPrYU7yZmY9zEnezKyHOcmbmfUwJ3kzsx7mJG9m1sOc5M3MepiTvJWWpJdHuX8/SX88Spt3\nSPpm4r5vSZreTB/Nys5J3spstCv1ZgIXj3U/EfGeiHip4V6ZdREneSs9SftIWibpe/kKRr+W33U1\ncGS+Mtana+xiv/yo/VFJ11ft9wlJ++c1bdZK+jtJj0i6Q9KkvM1HJK3JSzp/uY0v06wtXNbASkvS\nSxExXdIEYEpEbJF0AHB/RBwt6XDgmxHx8zX28Q7gX4F5wJPAncDfRsTt+YpIbwb2JVtc5eSIWC3p\nK8DiiPiypE3A3IjYIWm6j/yt2/hI3rqBgKslPUy2nOEhkhpZxGJ5vipUAP8EvL1qvxVPRMTqfHsF\nMDfffhj4sqTzgeGxvgCzTnGSt25wPnAgcGJEnEhW23xyA4+vZwGMV6u2h9lVofVXgf8FnAQ8mK8X\na9Y1/Ia1Mqscae8HPB0ROyW9Ezg8j79MNtQymrfk4+57Ab8J/FuN5xrpsIi4B7gcmA5Mq7v3ZiXg\nJG9lVjnivgU4JR+u+W1gHUBEPAfcK+n7o5x4XU52NL4G+I+I+MaI/Y/cBkBSH/CP+fOuAK71mLx1\nG594NTPrYT6SNzPrYV7+z3qCpJ8FbmbXsIuAbRHxts71yqzzPFxjZtbDPFxjZtbDnOTNzHqYk7yZ\nWQ9zkjcz62FO8mZmPez/A6sEKLwKqdB6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x112695438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dr_zm1.plot(cmap='jet', yincrease=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A fast binning method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 464 ms, sys: 5.59 ms, total: 470 ms\n",
      "Wall time: 467 ms\n"
     ]
    }
   ],
   "source": [
    "def zonal_mean_fast(dr):\n",
    "\n",
    "    dr_grouped = dr.groupby_bins('lat', lat_bins, labels=lat_center)\n",
    "    \n",
    "    # Applying all operations within the group object, to avoid converting \n",
    "    # back and forth between normal DataArray and group view.\n",
    "    def mean_weighted_by_area(dr):\n",
    "        weights = dr['area']/dr['area'].sum()\n",
    "        mean_weighted = (dr*weights).sum(dim='stacked_tile_y_x')\n",
    "        return mean_weighted\n",
    "\n",
    "    zonal_mean = dr_grouped.apply(mean_weighted_by_area)\n",
    "    \n",
    "    return zonal_mean\n",
    "\n",
    "%time dr_zm2 = zonal_mean_fast(dr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is more than 10 times faster!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make sure results are identical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dr_zm1.identical(dr_zm2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}