Block Average of 2D arrays.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Block average of 2D arrays in python\n",
    "\n",
    "Demonstrate a method to block average 2D arrays in python, handling masked values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pylab as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two routines needed to perform block average."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def blockshaped(arr, nrows, ncols):\n",
    "    \"\"\"\n",
    "    from: http://stackoverflow.com/questions/16873441/form-a-big-2d-array-from-multiple-smaller-2d-arrays/16873755#16873755\n",
    "    Return an array of shape (n, nrows, ncols) where\n",
    "    n * nrows * ncols = arr.size\n",
    "\n",
    "    If arr is a 2D array, the returned array looks like n subblocks with\n",
    "    each subblock preserving the \"physical\" layout of arr.\n",
    "    \"\"\"\n",
    "    h, w = arr.shape\n",
    "    return (arr.reshape(h//nrows, nrows, -1, ncols)\n",
    "               .swapaxes(1,2)\n",
    "               .reshape(-1, nrows, ncols))\n",
    "\n",
    "def coarsen_field(img, by=3,):\n",
    "    \"\"\" coarsen the spatial resolution of an image by block averaging \"\"\"\n",
    "    # coarsen by block averaging \n",
    "    tsize = (np.array(img.shape)//by).astype(int)\n",
    "    c = blockshaped(img[:tsize[0]*by,:tsize[1]*by],by,by).mean(axis=1).mean(axis=1).reshape(tsize)\n",
    "    \n",
    "    return c\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare dummy data for demonstration (from https://scitools.org.uk/cartopy/docs/latest/tutorials/understanding_transform.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon = np.linspace(-80, 80, 35)\n",
    "lat = np.linspace(10, 70, 55)\n",
    "lon2d, lat2d = np.meshgrid(lon, lat)\n",
    "\n",
    "data = np.cos(np.deg2rad(lat2d) * 4) + np.sin(np.deg2rad(lon2d) * 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original shape:  (55, 35)\n"
     ]
    }
   ],
   "source": [
    "orig_shape = data.shape\n",
    "print(\"Original shape: \", orig_shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKwAAAEICAYAAADCyx+UAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeT0lEQVR4nO2de7AkV13Hv7/pmbt3X8nmwWPJBjbC8ohoQGN4WiIkVojGICUUoAKKRZWKBRrUBATZsqRWRMoqUIsoKYLh/Q4UiFspIoUVIA+SsJCEEEhgzSZLYnb3PnbvnZn++Uf3Zs/5npk+3ffR9557f9+qrZ0z3X369O+e6fnM9/z6HFFVmEypqLPSDTCZmsg6rCkpWYc1JSXrsKakZB3WlJSsw5qSknXYiETkLSLy70u9b426VESeVHPfd4jI1Utx3tWu7ko3oE2JyGsBXArgiQCOAPgsgMtV9dC4Y1T1nXXrb7LvSklEPghgv6r+9Uq3ZSFaN3dYEbkUwN8D+AsAJwN4NoAnANgrIhNjjllXH+gkpKpr/h+AkwBMA3g5vb8FwEEAf1CW3wHgUwCuRnEH/sPyvaudY14N4F4ADwF4G4B7AJzvHH91+XonAAXwGgA/BvAggLc69ZwH4HoAhwAcAPA+ABPOdgXwpDHXcxaA/wYwBWBveazbxk8CuB/AYQBfA/Cz5fuvB9AHMF/G4wvl+5cBuLus73sAfmul/2bj/q2XO+xzAUwC+Iz7pqpOA/gygAucty9B0Wm3Afiwu7+InA3gXwD8DoDtKO7UZ0TO/XwATwHwIgBvF5Gnle8PAfwZgNMBPKfc/sc1r+cjAG4qj/1bFB8KV18GsAvAowHcfPw6VPWK8vW7VHWLql5c7n83gF8ur2c3gKtFZHvNtrSq9dJhTwfwoKoORmw7UG4/rutV9XOqmqvqUdr3t1Hclb6uqvMA3o7iTlil3ap6VFVvBXArgHMAQFVvUtVvqOpAVe8B8H4AvxK7EBF5PIBfAvA2VZ1T1a8B+IK7j6peqapTqjqH4q5/joicPK5OVf2kqt5XXvPHAdyF4htg1Wm9dNgHAZw+hkm3l9uP6ycV9TzO3a6qsyjQoEr3O69nUWAIROTJIvJFEblfRI4AeCf8D05VGx5W1RnnvXuPvxCRTET2iMjdZb33lJvG1i0irxaRW0TkkIgcAvD0mm1pXeulw14PYA7AS903RWQzgBcDuNZ5u+qOeQDADuf4jQBOW2Cb/hXAHQB2qepJAN4CQGocdwDAKWXbj+vxzutXocCa81F8xe883tzyf+/6ROQJAP4NwBsAnKaq2wDsq9mW1rUuOqyqHkbBZu8VkQtFpCciO1H8ONkP4D9qVvUpABeLyHNLZ2E3Fv6H3Yrih920iDwVwB/VOUhV7wVwI4DdIjIhIs8HcLGzy1YUH86HAGxCced29QCAn3HKm1F04p8CgIj8Poo77KrUuuiwAKCq70JxF3s3io7yTRRf7y8qWa9OHd8F8KcAPobiTjeFwmWodTzpzSjuhlMo7nAfb3DsqwA8C8D/AfgbAB9ytn0IBSL8L4pf/N+gYz8A4Ozy6/9zqvo9AP+I4lvoAQA/B+B/Gl9NSxK1BO4FS0S2oLCldqnqj1a4OetC6+YOu1QSkYtFZFPJkO8G8B2c+GFjWmZZh22uSwDcV/7bBeAVal9TrcmQwJSUFnWHLX9x3ykiPxCRy5aqUSbTOC34DisiGYDvoxjW3A/gBgCvLH91jtTJp2b62B29R8qTknvbu+QQCZWH8Pefo6Yfzf0clpl8g7990PPK/WHmlXVADlVO7aHzBZHjj3/m75F1h155snti4G1TNu9t2yS+8TAp/rE98U8Wi9U8x0r9WMwMKVZDf/v8wB9z0SHFarh0sRo8+DCGUzMj7cLFZCOdB+AHqvpDABCRj6Hgu7Ed9rE7enj/NWc+Un5qb8bbflpno1fO6I9yOPdHSn/U97ffNucP639r6oleed/D/vD4fQ/5o5X9w/4fLZv26+9Qh1b+m230OwlO8keCt5067ZWffNpPH3n9iyff6237hY33eOWn9A575cdkfqx64n/4pvNjXvnegd9lvjP3OK98w/RZXvm2Q34sf/LQNq889/CkV86m/fN3+hQr+Mon/XdyJ1b3734vxmkxSHAG/GHM/RiRCCIirxeRG0XkxsMPDXmzydRIi+mwo27ZwZ1fVa9Q1XNV9dyTT8tGHGIy1ddikGA/gDOd8g4UVs9YqQqO5SfYaCb3+/cW8b9CM/pMzKn/lXtMfWadJWady/3L6+f+5zOnstBXfliGL/q48/45lZmZjzlMPUsMeYwYk3l9Tvv+ucCxGVKZmDWIFZ2PmHU4oJsNM2vTWPGXrXt8xc+qxdxhbwCwS0TOKsfVXwHgmkXUZzJFteA7rKoOROQNAL4CIANwZTnWbjItmxb1zJKqfgnAl5aoLSZTVK0+ZJdDMKsn2GlWfS6aJS7LyDdi5p2JMGzMS2TGDDiMXKrA1ur47WEfNyfOGxAHHhueaM9sxEOeUb/ts8SoPeL7Wa2OFTNtECvibb6WWKyE96e/nWYVPu4yMazJ1Lqsw5qSknVYU1JqlWGH6HhsFnKZP37eIZhh5mXOYw50GREA5gecO8BDr357eXgx9BZp+5DgizhuSOebc9rHDMk8Pks+6Qw1pkeD9zPkMYf1+bEKeZ8YNohVdWw6/s8RKMeq67fXY2Ie83brHbvFZFqFsg5rSkrWYU1JqV0fVsXjTuaqKQIhZtipvNp3nR3S9oFf7vN4eL86l4CZlse/2XvMeX+qfzj0y67XeZTaPjX00weP5H4639aOnz/bo1yCqUjuwPTQry+IFfmwiDBswPscKyqzD8uxHye7w5qSknVYU1KyDmtKSi3nEnQ8dmIu29TxfdgM1ePhR3Kf86Ypp3SOfNgBcRiPdwfMGvEW+ePOHNehx1KGxLTH+ifaNzNgPqeyEu/T9h5BYvD7gGLNvM+edb/PDBvxXSNlTvfneSS92FkugWmtyDqsKSlZhzUlpXZzCVQw5TBsjLOYYacCZvX3nxnQc1GR55IkYM5qpg3GxwnzAgam+tn3dfNz2Qedoms7QmWet2CCGJZ/H7Cvy/MQHKM5G+KxQnWZeZ+nfKDYuc3nOQ28esdvMplWn6zDmpKSdVhTUmo5l6DjeaXso07mPvh0aHycOY45b4a9xT49w0UMmfH4N3EXcxhzGj1GBUrvDb1I8n1dr5Pn/WJPmfl9c065BGSM8v4cq2li5pD3I3kWzLSR2DHDdpj/3f2NYU1rRdZhTUnJOqwpKbX8TJd4/h/7qJPic1kmMR+W8mGJy3geAvZBm3qJnT49W08fd+1GckTn/QPceQqYYdlT5lhNUd4F5xIEzBqJ1RzxvjLvzzd7hotjFTJsRayMYU1rRdZhTUnJOqwpKbX+TJfr/zFnbSAQysiHDXIHeF5+XsMgktPJ3EXLDAQcxttDhuXjq8uuL8w+KOfHcqw28TNdgQ/bjGF5HoIY78ditRgf1nIJTGtG1mFNSSnaYUXkShE5KCL7nPdOFZG9InJX+f8py9tMk6lQHYb9IID3wV8x+jIA16rqnnJBucsA/FWsolw7HjsdGZAPS+DToQf/Y+PhR/s+w/J8rDxPQGw8PMZpPF9UZ545rjrf1n3Gi31QZkxm0E0dPxY9gsxYrvBsn+Yh6EfyX2O5AxSrbL46VlnGz7stkQ+rql9Dscy5q0sAXFW+vgrAS2L1mExLoYUy7GNU9QAAlP8/etyO7jpdc4eOjtvNZKqlZf/R5a7TtWHbxvgBJlOFFurDPiAi21X1gIhsB3CwzkG5isdmzGUbBjS3FhlyzLzBeDivYcDMOh/L4WQGJR82mPOUlp8Mcgn8/Tnf1p2viuf9Yh82yC1gH5ZM4MMDfoYrkv8aZVhQmWMTiZVUr3Gw3D7sNQBeU75+DYDPL7Aek6mR6thaHwVwPYCniMh+EXkdgD0ALhCRu1Cs5r1neZtpMhWKIoGqvnLMphctcVtMpqhazyVw/T/msg2daoYNOK5Pz9ZzTuc8LYkezeH0y8xhsXzYLMgliOXHnij35/2D2VOeisSq2/HzYWOxYt83J4btxuYhiD3/FsmHDeZ0sGe6TGtR1mFNSck6rCkptb7WrOv/MVdFGTbwEheX/xrmDkR82Dk/t4HHx0NvMZJb4JSDNbwG1bkFMwTY3dxvG/u4nCsc5r82zRWgcpA7TLESqp8Zdt6e6TKtQVmHNSWlVpFAVTw7JXi0uet/jQVLd0am1wnTCWMpctWPdTACdPhrrupRZYSPRodDtSe282PVc3QtQawyPxZdSsWMxYptNOFljZoOzXKsGJ8oVlmXp41y2mJIYForsg5rSkrWYU1JqWWG9e0U5rIe2VzMZTP96ul1gnTCplYNMy0xa9YnLqPpM2NTFYWPzDgFnk4+MlS7IcKwHFsetg6m02yYehnEKthOsRo0sACNYU1rRdZhTUnJOqwpKbXuw7rDp+wN9jKfu5jLeH8eXmQvM/aYR5A+OF/NYUJMK/yYd2+8tzjq/C7T8pTsvOwQD9Uyo8YYlqceXXw6YTXvB7HiqYqI9zPnb2c+rGnNyDqsKSlZhzUlpdZ9WJfN2EedyPzHPPoEM8xhYTph9RSR4ZLp7CVWc1hnniug8/XYB455lw7D9plhq3MLmOd5ev0gdyC2nLwf+hG8z8xaHbtYrPj3ge/DjodYu8OakpJ1WFNSsg5rSkqtMixUvPH+wbD6sRAhLmOOiy6R3jCnMyz7YCdzBHpd3p8eK4/k27qMzW0dkKfcJ94/Rryf0bRJ/PuAc4Vjjw/FlzViZqVYUexAPrF0x/O++bCmNSPrsKakZB3WlJRaZlj4U0wS881llH9KMMMcp5El0pnDwikh2VskDiMvUfpUYc4+LDMscV5Fe9kHZcZkX7Y/9M+V5/618O+DnM8dzR2m2AxieRYUqzmGYM67oFi5nrdftb/f+E0m0+qTdVhTUqozP+yZIvJVEbldRL4rIm8s37elj0ytqw7DDgBcqqo3i8hWADeJyF4Ar0XTpY8UHpvlxGGDgQ8vQ2JY5jh+ToiXFYrmEgScRlzGPiwxrObMhX44w/oxtsxt43kC+NrZV807zPvVucJZNM+Cy379wgwbixUxrARToZ54Ro1SZ/39KrYVJ1I9oKo3l6+nANwO4AzY0kemFVAjhhWRnQCeCeCbqLn0kbvs0XB6ZpHNNa131e6wIrIFwKcBvElVj9Q9zl32KNuyeSFtNJkeUS0fVkR6KDrrh1X1M+XbzZc+UvG8x5yYc5ARExLMDMlbDPJfOaczwmUScBlVwL4rlYUYFnR86MPS+dxcAooF+7Dsow742nkJpmEk/5V5Pxa7IC+CzNJYrIhhEfC+U99i8mFFRAB8AMDtqvoeZ5MtfWRqXXXusM8D8HsAviMit5TvvQXFUkefKJdB+jGAly1LC00mR3WWPfo6xjsNtvSRqVW1m0sAQIYuw9J4Nz0TxZ8SnstKuMyMylw2rB4fj3mJ6FNSaM7TrhO3DaoZ1i0HeRDUdr52ZlTNmWEpenl1/UGsKtoKLCBWGTN1Rf6s5cOa1oqsw5qSknVYU1JqPR/W40zmMs7Z5CXLeXvAYZH5oYIyjYdTmTlLB1SBsq9KPmxQP7fXeU2nZl824H2umxAxmivMsYvwP3vWaBqrIbW3R8/vDZfIhzWZVpOsw5qSknVYU1Jqn2FdziRvkL3GMJmgejw8ZFqtLgdcFuGwCMNiSA1gzmQfeDj6NRAyLT/npJxXEbSl2rPuBGVEynwtfK2RWAnngZAP68TK5iUwrRlZhzUlJeuwpqTUKsMKfD5hrgrWveJkAh4PJ64LGJa9zYBp2SetZlgN8mGrGZg5Lzz/6Ncjy8z7zLSceRFh4OXmfY4Vzw/LsfE8cPNhTWtF1mFNSck6rCkprUA+rFMIOIt9WN7OdTVk2gi3IacK2GukMpNWwGURDnSZlq8liA1v5zI3Jsh/ra4/iF0QS/Zhed3dSOyU80JoO8d+jOwOa0pK1mFNSck6rCkptT8/rIsqzKwBwxIDMpcRt0U5jH3TgFmruYzL4TNnVB/Vz+d32xe2lfmcGJBN6qDuCN83jFVT3g9iwbnDOf8+cH1YjJXdYU1JyTqsKSlZhzUlpfZ9WDeXoBpzEFAiW3XMsHx84E3S/o25jBcioPqquAyA8AXqmNcYkQ8buzaOVWT/IFYRDzvMt2W+r46VRJjbj73lEpjWiKzDmpKSdVhTUmqfYV1/sIrpRok5KMZpzGXB/hGG5XKQhBrhsqD+8e2J8XfMp9Ws/rlG1h+UI9cSLRPj8ly6zLgcuzGyO6wpKVmHNSWlOjNwT4rIt0Tk1nKdrt3l+7ZOl6l11WHYOQAvVNXpcq2Dr4vIlwG8FE3X6SKF3ERcRqAV7k8VRhm4GYdpbHvArBFjuaocu7aYGu4fjV3Ep43yf2CqR/JdK2Lhqs46Xaqq02WxV/5T2DpdphVQLYYVkaxc3+AggL2qurB1umZsnS7T4lSrw6rqUFWfAWAHgPNE5Ol1T+Ct07XZ1ukyLU6NfFhVPSQi1wG4EAtapwuNWStan6Mo466wKtvX2IOOeNgxz7qhAl92keLfB1Xry7qq4xI8SkS2la83AjgfwB2wdbpMK6A6d9jtAK4SkQxFB/+Eqn5RRK6HrdNlall11um6DcWCyPz+Q7B1ukwty0a6TEnJOqwpKVmHNSUl67CmpNRuPqygvuFWtz5HgfW4lOdaAlW2L9bWwMTl7bH9FxcMXeTxLFlgfXaHNSUl67CmpGQd1pSUWn+my1XIdNWcFmXUKAdG1gGjMnMWc5x0+Hj6/Dc532L5u+H+0dhxmW9tfC2d2LVH7o0Vsahqhsm0qmUd1pSUrMOaklLrDKsd92F82rhYb5E5kJdjDeqLcFjAZZFywLRc//j2xJgyuJZOJD81eu3V2wPftXGZ1uXi2GSZfz7ePkZ2hzUlJeuwpqRkHdaUlNpnWJfbAm+Py8RpASOOr3tkfcx1WYRhibMkKDPD+ttB20MuHN/WeGy43JDvo7Hh+plBiVEjsQr+dkHZrW88z9od1pSUrMOakpJ1WFNSaj8f1v2IsJfI5YDraG6rqFfJ5RizMpdVbw8Zl5i1w2XKTeiMfl2Udey+xck4ds2Obxqr5rwf8bSreN9yCUxrRdZhTUnJOqwpKbXvw7rowh+XCMOyd8fz+ofcVXFuNOeymC8bjI9T/dr1y7mzPfA9OTa8hgGXeX+erjXYfzxPjywznzfke46tdDmXoN690+6wpqRkHdaUlKzDmpJSqwyroFyCCJfx8LUOI94iMyptz5kpmcuIq9D1wyM9WoCVOY25jLaH5x/9emQ58Kz9ogS+ayQXN8LzMf7maxWKFXpU5txh5v2u68NaLoFpjah2hy3XOfi2iHyxLNuyR6bW1eQO+0YAtzvly1Ase7QLwLVl2WRaVtViWBHZAeDXAfwdgD8v374EwAvK11cBuA6xdbqEuJW4SwKGJS4LuIq3g8pLy2U6jDBswLTsXbIPO/o1MCJ3gJk1y6nMa4hVe9Z5UB7ftqLM1xJhWC4zwzLvO7GqmpOh7h32nwD8JXw72pY9MrWuOoty/AaAg6p600JOYMsemZZSdZDgeQB+U0QuAjAJ4CQRuRoLWfbIZFqk6izKcTmAywFARF4A4M2q+rsi8g8oljvag7rLHgmg7hmJozrdnHf320Lbw7F74jIeu6erzZkxu9Xj3TKkCgIvkXzXoH7/cM+HZc+Y+Bx07Z1uhPdpf762gPcr2lZsZ96vZtKAaTlW5NN6nvgy+bB7AFwgIncBuKAsm0zLqqYrIV6Hwg2wZY9MK6IVSC888dUl/DXXISSgb4Y8Z+uFvwb5XHQ8Wzf0Nae96q8t4eXmOaWOv+boazNEFCe9MNL2AJ8y32Lr8NAseUND3h4ZCq5qK1AjVkPCNx4qZlvLrc8ekTGtFVmHNSUl67CmpNTyY97qDadmZM102dZiq4YZlvbPOZ0vsLH8Mls1zGVBilzOHEhWER2f96oZ1uXW0Maqtvw4Vsywee4z7jAYCo0MzfKl92g7XVuHY0UMGzy+VBUrSy80rRVZhzUlJeuwpqTU/lRFDjuxl9jt+uVMmMv8z9eA0wuZyyqYERjBZV3mMjYnq33YkMtiQ8Pj28ZDq8yssVhpz69v0CPeD843vm1FmRm0+trBw9jMsBP+ds+TxnjZHdaUlKzDmpKSdVhTUloBhj3BUj16bHoDcRn7sEMaH5/nR1gCpvVPPySm5HI+QV7lBIOcX+QUOz6+iQ8b5A7QtWTEsD3i/4x9WIoV+7hDZtrAZ63OJRjStWUcK62ewj5n3p9w6qu4jdod1pSUrMOakpJ1WFNSaj2XoOOwU5ceVd7QHXjlDnuLxGXHiHn7PeK4oBzhtKDsc1aHuYymiAyYNajPP9z1OrmtQozZ6/mxmaQy+7CsOeL9IXvW1LbQs+ZyNb93cs6HJd92YnysluIxb5NpVcg6rCkpWYc1JaVWGVYEyBzu3EAcxgzbFfIOCW7Yx+1zfmyQAwraXp3jGXAWh6vDPm41wwa+r8OF2qv2XdmjnozEinWUYjUXmeaJmTX0rKtjFyQjEJeyj+vFyvJhTWtF1mFNSck6rCkptcyw6nEnc9jGbt8rM5cNaD6fCeK6o+RdMheGTOm3j7ksGC/3dw+eU4rWH5RPtI/zILIu51k0ixVrpuufvENMm/fG56cW20Hlan4X8qyVcwlo/6H5sKa1KOuwpqRkHdaUlFr3YV3uZA7b0pvzyh16uicnM2+mN+HvH+R4Vo+X5/7hgdeYTVR/nnkK+GHEh63KJQC1fWKCmLXnx2pTd94rxxg2yD0gn3cwEYlV1FPmidAi0+UHPq7rw2Ks7A5rSkrWYU1Jqe4qMvcAmAIwBDBQ1XNF5FQAHwewE8A9AF6uqg8vTzNNpkJNGPZXVfVBp3x8na49InJZWa5c9khEvfwBZtjNmc9lnA/LDBuMp5N3Odc4P5aeIdsQYdiYD0uMHOScOufj/FfOHQhi1ZRhuxu8co8YOcjDiOb20rXGYiXVDOx61svlw16CYn0ulP+/ZBF1mUy1VLfDKoD/EpGbROT15XuN1+kaHJ5dfItN61p1keB5qnqfiDwawF4RuaPuCVT1CgBXAMDmJ2+vfo7DZIqoVodV1fvK/w+KyGcBnIcFrNPVgXrcyb7r5i75sMywygzrcx3nx3LOZ+gt+uUheYOE1OAvJF6qKGTW2Hi8O1du9fNt7Lsy73c7tKwoiRmY8zBmm3rWEQ87CA7btEEuwvh9XdVZCXGziGw9/hrArwHYB+AaFOtzAXXX6TKZFqk6d9jHAPisFL/yugA+oqr/KSI3APiEiLwOwI8BvGz5mmkyFaqzEuIPAZwz4n1bp8vUulrNJeiIYlPvBHsxs27JmjEsMzCPl89MUM4nPZsfZVoqM1wFDNuAWQFAnfF79kU5d2Arxeqk7lGv3IswLMeKn6dbfH6sXxb6W7G3yr8XLJfAtCZlHdaUlKzDmpJS+wzr+InMrCcTl3Ui4+OHuhu9cjA3F+fHBjmfsbm1+IzEoEEuASrLvO6AN1cu+aKcKxDj/V7Hv3bWVDbplTkPIyOGHcR82cjzazw/bJNY2TNdpjUj67CmpGQd1pSUWmbY3GPYk7rHvO1bMr+c0aICOa2LtYU4j73LMD+22bz+PF4e82HD55RoO+U2uDmw7Ity7gAz61aKVU+qGZYZ2PXDgRHzlDXNJY7EKpxbqyJWxrCmtSLrsKakZB3WlJREed7+5TyZyE8B3AvgdAAPRnZfKa3mtgHro31PUNVHjdrQaod95KQiN6rqua2fuIZWc9sAa58hgSkpWYc1JaWV6rBXrNB562g1tw1Y5+1bEYY1mRYqQwJTUrIOa0pKrXZYEblQRO4UkR+U83GtqETkShE5KCL7nPdOFZG9InJX+f8pK9S2M0XkqyJyu4h8V0TeuMraNyki3xKRW8v27W6jfa11WBHJAPwzgBcDOBvAK0Xk7LbOP0YfBHAhvXd8krtdAK4tyyuhAYBLVfVpAJ4N4E/KeK2W9s0BeKGqngPgGQAuFJFnL3v7VLWVfwCeA+ArTvlyAJe3df6Kdu0EsM8p3wlge/l6O4A7V7qNZVs+D+CC1dg+AJsA3AzgWcvdvjaR4AwAP3HK+8v3VptqTXLXpkRkJ4BnAvgmVlH7RCQTkVtQTFO1V1WXvX1tdthRWY7mqUUkIlsAfBrAm1T1yEq3x5WqDlX1GQB2ADhPRJ6+3Odss8PuB3CmU94B4L4Wz19XD5ST26HuJHfLJRHpoeisH1bVz6y29h2Xqh4CcB2K3wPL2r42O+wNAHaJyFkiMgHgFSgmlFttWhWT3EkxmdkHANyuqu9xNq2W9j1KRLaVrzcCOB/AHcvevpbh/CIA3wdwN4C3roIfCx8FcABAH8U3wOsAnIbi1+1d5f+nrlDbno8CmW4DcEv576JV1L6fB/Dtsn37ALy9fH9Z22dDs6akZCNdpqRkHdaUlKzDmpKSdVhTUrIOa0pK1mFNSck6rCkp/T+RARNquhHi7wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(data)\n",
    "plt.gca().set_title(\"Original data\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute block average of the data by a factor of 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKwAAAEICAYAAADCyx+UAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPVElEQVR4nO3de4xc5X3G8e/jXa/xZWFtbMB3jIJoIY1C5EITKIpCUZyAwG1BhQQSaFCSPyCERI2gikqVNlVUVcSpQpJSh0AFAaUEtRHXRgSXUC4CjCUwJmDMxYsXvMbYaxt82fWvf5zjMF5mdtd73t2Z134+0sozc878zm9mH589c3nPq4jALBcTmt2A2YFwYC0rDqxlxYG1rDiwlhUH1rLiwGZK0gpJl49w3U9K6h7rnsbDQRNYSZ+T9JSk7ZJ6JN0n6fRm95UbSZdKeqTZfTRyUARW0jeAZcA/AUcDC4AfAecl3o4kHRTPWbYiIusf4AhgO3DBEOtMogj0hvJnGTCpXDYduBvoBd4pL8+rue8K4LvA/wHvAR8CLgXWAduAV4DP16z/18CastYDwMKaZQF8FXipXH4DoBHe9yzgBWAr8EPgf4HLGzzeycDNZZ3ngb8BumuWXwO8XPb/PPDn5e1/COwEBsrndEt5+9nAM0AfsB74+6b9vpsduASBXQL0A+1DrPMd4HHgKGAW8CjwD+WyI4G/BKYAncB/Av81KLCvAycB7eV/kD7ghHL5bOCk8vJSYG35i28Hvg08OiiwdwNdFH8FeoElw90XmFlu83xgInB1+ZgbBfZ7wG+BGcB84LlBgb0AmEPxF/avgB3A7HLZpcAjg+p9Evijcv2PAG8BSx3Y0QX288Cbw6zzMvDZmuufBl5tsO5HgXcGBfY7NdenAlvKkE8edN/7gC/VXJ8AvEu5pywDe3rN8l8A1wx3X+ALwOM1ywR0DxHYdfv+I5TXv1wb2DrrrwLOaxTYOusvA77fjN/3wXA89jYwU1L7EOvMAV6ruf5aeRuSpkj6N0mvSeoDHga6JLXVrL9+34WI2EGxV/oq0CPpHkl/UC5eCPxA0hZJW4DNFOGaW1PrzZrL7wLTRnDfOYN6iNrrDR5v7fLax46kL0haVbOtD1PsxeuSdKqkhyT1StpaPvaG64+lgyGwj1Ecdy0dYp0NFIHYZ0F5G8A3gROAUyPicOCM8nbVrL/fV9oi4oGIOIvicOAF4N/LReuBr0REV83P5Ih4dASPY6j79lD8aS8ak1R7vY6eQcsX1Nx3YdnvFcCREdFFcciw7/HW+/rez4FfAfMj4gjgJ+z//Iyb7AMbEVuBvwNukLS03GNOlPQZSf9crnY78G1JsyTNLNe/tVzWSfFiaoukGcB1Q21P0tGSzpU0FdhF8eJkoFz8E+BaSSeV6x4h6YIRPpSh7nsPcJKkvyj/knwNOGaIWr8oa02XNA+4smbZVIpQ9pbbuYxiD7vPW8A8SR01t3UCmyNip6RTgM+N8DGl14zjkLH4oTiWfYriBcSbFL/kT5TLDgP+lWLP01NePqxcNofiOHU78CLwFYpfaHu8fwx7ec12ZlO8Qt9KcSy7AjixZvklwLO8/4r6ppplAXyo5vrNwD+O8L5Lyv5G8i7BFOA/yv7qvUvwXYpDjk3A9bW1gI7yudsMbCpvO5/isGIbxYvGHwK3NuP3rLIhsyxkf0hghxYH1rLiwFpWHFjLylBvtic3c0ZbLJyfbpPr90wbfqUDsLVvStJ6Uzt3Jqt1XMe2ZLUAXtx5RNJ6u97tGH6lEerfvJmBHTvqvs87roFdOL+dx++fl6ze1T2nJqsFcN+Di5PW++PTX0hW6+eLHkpWC+DTa85JWu/llUN9jnFg3lj2/YbLfEhgWXFgLSsOrGXFgbWsOLCWlUqBlbRE0u8krZV0TaqmzBoZdWDLLzjfAHwGOBG4SNKJqRozq6fKHvYUYG1ErIuI3cAdJB6lajZYlcDOZf9hGN3sPxQEAElfLs8X8NSmt/dW2JxZtcDW++jsA1+ujYgbI2JxRCyeeaRf41k1VRLUzf7jhubx/jgpszFRJbBPAsdLWlSO/7mQYqCa2ZgZ9ZdfIqJf0hUUZyhpoxh/tDpZZ2Z1VPq2VkTcC9ybqBezYflVkGXFgbWsOLCWFQfWsjKuQ2S27RUrdk5MVm/V2+mG2wBM3pj2dFGre4c6m9CBeXh2slIAvNI7I2m9yW+le+4m7BliWbKtmI0DB9ay4sBaVhxYy4oDa1lxYC0rDqxlxYG1rDiwlhUH1rLiwFpWHFjLigNrWXFgLSsOrGXFgbWsOLCWFQfWsuLAWlbGdUxX38AU7t/6kWT11m9IOy5pzvq0Z1fc0NOZrNY9Cz+arBbA3jfSzknW2Z3uuWvzmC47WDiwlhUH1rLiwFpWHFjLigNrWaky7dF8SQ9JWiNptaSrUjZmVk+V92H7gW9GxEpJncDTkn4dEc8n6s3sA0a9h42InohYWV7eBqyhzrRHZiklOYaVdCxwMvBEnWW/n6fr3S27UmzODmGVAytpGvBL4OsR0Td4ee08XVO6JlXdnB3iqk6OPJEirLdFxF1pWjJrrMq7BAJ+CqyJiOvTtWTWWJU97GnAJcCnJK0qfz6bqC+zuqpMLPcI9eebNRsz/qTLsuLAWlYcWMvKuA6R2d7fwWMbFyWrN7GnI1ktgGmv70hab/Ib05LVeqL32GS1ACb3pN1XTXst3XM3YVfj4Tbew1pWHFjLigNrWXFgLSsOrGXFgbWsOLCWFQfWsuLAWlYcWMuKA2tZcWAtKw6sZcWBtaw4sJYVB9ay4sBaVhxYy4oDa1kZ1zFde/rb6Nl0RLJ60zamPS1C+/pNSetN3jg1Wa3uTV3JagEc+VYkrTcx4XOn3f0Nl3kPa1lxYC0rDqxlxYG1rDiwlhUH1rKS4pTxbZKekXR3iobMhpJiD3sVxQwyZmOu6hwH84CzgeVp2jEbWtU97DLgW0DD083VTns0sC3t2QHt0FNlUo5zgI0R8fRQ69VOe9TWme6jSjs0VZ2U41xJrwJ3UEzOcWuSrswaqDJ157URMS8ijgUuBH4TERcn68ysDr8Pa1lJ8vXCiFgBrEhRy2wo3sNaVhxYy4oDa1lxYC0r4zqmiwGxd9vEZOU6tqYdlzSwMe2Yrkl9C5LVGuhLOyfZYVsGktYb6E333EW/x3TZQcKBtaw4sJYVB9ay4sBaVhxYy4oDa1lxYC0rDqxlxYG1rDiwlhUH1rLiwFpWHFjLigNrWXFgLSsOrGXFgbWsOLCWlfEd0xVCe9LNrTWh8dCfUYk9u5PWm7An3Zgz9aedk0z9DU84OSqxO+FzF42fN+9hLSsOrGXFgbWsOLCWFQfWsuLAWlaqziLTJelOSS9IWiPp46kaM6un6vuwPwDuj4jzJXUAUxL0ZNbQqAMr6XDgDOBSgIjYDaR9591skCqHBMcBvcDPyqk7l0v6wLxG+83TtX17hc2ZVQtsO/Ax4McRcTKwA7hm8Er7zdM1bVqFzZlVC2w30B0RT5TX76QIsNmYqTJP15vAekknlDedCTyfpCuzBqq+S3AlcFv5DsE64LLqLZk1VimwEbEKWJymFbPh+ZMuy4oDa1lxYC0rDqxlZXzHdCmIjnTjnAbSTl3FhMMOS1pvYFK6cVgxMe28WnsT9gYwYfLkZLX0XuP9qPewlhUH1rLiwFpWHFjLigNrWXFgLSsOrGXFgbWsOLCWFQfWsuLAWlYcWMuKA2tZcWAtKw6sZcWBtaw4sJYVB9ay4sBaVsZ1TJfag4ldO5PV2zU97cnlJhxzVNJ6O7vS7Q86ut5NVgvgvekfONFkJZ1Hz0pXrLtxLL2Htaw4sJYVB9ay4sBaVhxYy0rVaY+ulrRa0nOSbpeU9tQpZoOMOrCS5gJfAxZHxIeBNuDCVI2Z1VP1kKAdmCypnWKOrg3VWzJrrMocB28A/wK8DvQAWyPifwavt9+0R307Rt+pGdUOCaYD5wGLgDnAVEkXD15vv2mPDk/76YodeqocEvwZ8EpE9EbEHuAu4BNp2jKrr0pgXwf+RNIUSaKY9mhNmrbM6qtyDPsExWRyK4Fny1o3JurLrK6q0x5dB1yXqBezYfmTLsuKA2tZcWAtKw6sZWVch8h0tPezaNbmZPXWHZX2g4jd849MWu+9o9PVWjTr7XTFgPXHpB1etHvBjGS1YmNbw2Xew1pWHFjLigNrWXFgLSsOrGXFgbWsOLCWFQfWsuLAWlYcWMuKA2tZcWAtKw6sZcWBtaw4sJYVB9ay4sBaVhxYy4oDa1kZ1zFdne27+NNZa5PVe3luwql2gG0LJyWtt2vu7mS1zpiZ7nkDWD5nXtJ6fQvSnct64NnG+1HvYS0rDqxlxYG1rDiwlhUH1rLiwFpWhg2spJskbZT0XM1tMyT9WtJL5b/Tx7ZNs8JI9rA3A0sG3XYN8GBEHA88WF43G3PDBjYiHgYGn8HtPOCW8vItwNK0bZnVN9pj2KMjogeg/PeoRivWztO14510n/zYoWnMX3TVztM1dXrHWG/ODnKjDexbkmYDlP9uTNeSWWOjDeyvgC+Wl78I/HeadsyGNpK3tW4HHgNOkNQt6UvA94CzJL0EnFVeNxtzw369MCIuarDozMS9mA3Ln3RZVhxYy4oDa1lxYC0riojx25jUC7w2glVnApvGuJ3RauXeoLX7G2lvCyOi7oC9cQ3sSEl6KiIWN7uPelq5N2jt/lL05kMCy4oDa1lp1cDe2OwGhtDKvUFr91e5t5Y8hjVrpFX3sGZ1ObCWlZYKrKQlkn4naa2klhonJmm+pIckrZG0WtJVze5pMEltkp6RdHezexlMUpekOyW9UD6HHx9VnVY5hpXUBrxI8XXFbuBJ4KKIeL6pjZXKL6rPjoiVkjqBp4GlrdIfgKRvAIuBwyPinGb3U0vSLcBvI2K5pA5gSkRsOdA6rbSHPQVYGxHrImI3cAfFYMeWEBE9EbGyvLwNWAPMbW5X75M0DzgbWN7sXgaTdDhwBvBTgIjYPZqwQmsFdi6wvuZ6Ny0UiFqSjgVOBp5ociu1lgHfAvY2uY96jgN6gZ+VhyzLJU0dTaFWCqzq3NYaxys1JE0Dfgl8PSL6mt0PgKRzgI0R8XSze2mgHfgY8OOIOBnYwSjPZdFKge0G5tdcnwdsaFIvdUmaSBHW2yLirmb3U+M04FxJr1IcSn1K0q3NbWk/3UB3ROz7i3QnRYAPWCsF9kngeEmLyoPyCykGO7YESaI4BlsTEdc3u59aEXFtRMyLiGMpnrffRMTFTW7r9yLiTWC9pBPKm84ERvVidVxPGT+UiOiXdAXwANAG3BQRq5vcVq3TgEuAZyWtKm/724i4t3ktZeVK4LZyZ7QOuGw0RVrmbS2zkWilQwKzYTmwlhUH1rLiwFpWHFjLigNrWXFgLSv/D23E85Muqn1aAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New shape:  (11, 7)\n"
     ]
    }
   ],
   "source": [
    "bavg = coarsen_field(data,by=5)\n",
    "\n",
    "plt.imshow(bavg)\n",
    "plt.gca().set_title(\"Coarsened data\")\n",
    "plt.show()\n",
    "print(\"New shape: \", bavg.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demonstrate block averaging with masked cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.ma.asarray(data)\n",
    "data[5:25,15:25] = np.ma.masked\n",
    "bavg = coarsen_field(data,by=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALYAAAEICAYAAADyRw61AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfu0lEQVR4nO2de7AlVXXGv+885sEMD3npCApGwWcpKgESLUUFRdRAGSWSaDRFYpUVLY0aBY2WpNRgYuG7ohgV4hsfRKSMika0iCgI8gwgPgYhjIwMDDP3zsyd81j5o/syvb99Tu9z7jnT90xn/apu3bP7sXv3Out0f3v16r1pZnCcutFY7gY4zu7AHdupJe7YTi1xx3ZqiTu2U0vcsZ1asuyOTfJtJP9t2tuOUJeRfNSEddxE8viS9ZeR/OtJjjFtSJ5P8t1TqOdVJC+fRpsKdR5P8s5CeT3JE5ZSV2t6zcpOFsCbADwSwBYAFwE4y8w2D9vHzN47av3jbFsFZvb4xc8k3wXgUWb28uVrkbPI1K7YJN8E4H0A/h7AvgCOA3AYgEtJrhiyz1R/WI6zyFQcm+Q+AM4G8Doz+7aZdcxsPYDTkDn3y/Pt3kXyqyQ/R3ILgFflyz5XqOsvSd5OchPJdxRvR8VtSR6ey4lXkvwtyXtIvr1QzzEkryC5meQGkh8d9gOTc3kWyRsK5e+RvLJQvpzkqfnn9SRPIHkSgLcB+DOScySvK1R5GMn/JrmV5HdJHjjkuMeTvJPkW0huzNt8KsmTSf6C5L0k3zbK+THjA3k995O8nuQTBhxzb5I/IPnhfJ/HkLw0P9atJE8rbHsAyYtJbsnt8ciEHZ9O8sd5++7I7+YguZLk+/Pv7G6SHye5uqyuwvn+LD/+3STPLd3BzCb+A3ASgC6A1oB1FwD4Yv75XQA6AE5F9qNanS/7XL7+cQDmADwdwAoA78+3P6Gw/+K2hwMwAJ/M63kSgAUAj83XPxXZXaOVb3szgDcU2mXIpIO2dxWA7QAOzPf9HYC7AOydH2c7gAPybdcPaluhrssA/ArAkfm+lwE4Z4gNj89t+E4AbQB/A+D3AL6QH/vxAHYA+IPU+QF4HoCrAewHgAAeC2Bdvu58AO8GcACAKwG8O1++BsAdAP4qr/MpAO4B8Ph8/ZcAXJhv9wQA/wvg8iHn8nAAWwGcnp/LAQCOytd9EMDFAPbPz+ubAP6pYIM7C/UU7XsFgFfkn9cCOK7MJ6clRQ4EcI+ZdQes25CvX+QKM/sPM+ub2XbZ9iUAvmlml5vZTmRfciqZ5Wwz225m1wG4DpmDw8yuNrOfmFk3v3t8AsAzUydiZjsA/AzAMwAcDeB6AJcDeBoyR7rNzDal6inwGTP7RX6uFwI4qmTbDoD3mFkHmSMdCOBDZrbVzG4CcBOAJ45wfh1kTvMYADSzm81sQ+E4DwXwQwBfMbN/yJe9EMB6M/tMXuc1AL4G4CUkmwD+FMA7zWzezG5EdsEaxl8A+J6ZfTG/e28ys2tJEtkP9u/M7F4z2wrgvQBelrDh4jk9iuSBZjZnZj8p23haGvceAAeSbA1w7nX5+kXuKKnnocX1ZraNZMqJflf4vA3ZrxkkjwRwLjLn3AvZuV6dqGuRHyK/euSf70PmNAt5eRwGtm8Im8ysl39e/NHfXVi/HSOcn5n9F8mPAvgYgIeTvAjAm81sS17PC5DdGT9eqPswAMeS3FxY1gLwWQAH5Z+L393tJefxMGR3KuWgvK1XZz4OILujNEvqWuQMAP8I4BaSv0F2Qbtk2MbTumJfgexLf3FxIck1AJ4P4PuFxWVX4A0ADi3svxrZbWwp/CuAWwAcYWb7INPALN/lARYd+xn55x8ic+xnYrhjV50mWXp+ZvZhM3sqMglzJLJO/SKfBPBtAN/KvyMgc9ofmtl+hb+1ZvYaZJKoi8xhF3l4SdvuwGANfg+yH+fjC8fY18zKfuyL53ObmZ0O4GBkQYqvFtoeMRXHNrP7kXUeP0LyJJJtkocD+Aqyq95nR6zqqwBeRPKP847Q2RjdGZW9kYUc50g+BsBrxtj3xwAeDeAYAFfmMuAwAMcC+NGQfe4GcDjJqp4NDD0/kn9I8liSbQDzyLR5T/Z/LYBbAVySX0AuAXAkyVfk3187r+ex+V3k6wDeRXIvko8D8MqStn0ewAkkTyPZyjueR5lZH9mP6gMkD87begjJ56VOluTLSR6U17E5X6zn9ABT+xLM7J+RXTXej8zgP0X2y32OmS2MWMdNAF6HTF9uQNYB2YjsbjAubwbw53kdnwTw5VF3NLN5ANcAuCnX+kB2V7rdzDYO2e0r+f9NJK9ZQnvHpez89smX3YdMMmxC9r08gGW9sFcj+46+gUzDPheZ3r0LmYR6H4CV+S6vRSaDfoesA/qZYQ0zs98COBnZM417AVyLvO8D4K0AfgngJ8wiY99DdhFJcRKAm0jOAfgQgJfl/aGB0Gb4RQOSa5H9Oo8ws98sc3OcPYhlf6SukHxRfrtbg+wqcwOysI/jjMzMOTaAU5DdCu8CcASyW87s3lacmWSmpYjjLJWJrth5BORWkr8keea0GuU4k7LkK3b+NOoXAE5EFtK7CsDpZvY/w/bZd/+mPeTQ9gPlVewH61sS2aOUewi3X5Cmb++HqSDz/ZXh+m47KHd64XMB60pksS/tkeNFltPLRDPcotkKo1OrWrueZe3V3Bms24thIGgVw33bElVM2Wqn2spCW8z3xFa9cP3Obvgsz3piq970bNW95z70ts4vNcwLYLInj8cA+KWZ/RoASH4JmT4e6tgPObSNT1y8K8b/mPZ8sP6ARpgL05Qv7/5++AT+N51w/fULhwTlK7eGzwhuvG9dUL5r075BuXN/+OU258L6G+L4pt/t6tCZsE/4EHa//eeC8pEH/P6Bz0/dN3yQ95TV64Pyo9v3B+UHN0NbtRn+SOf6YSTs9m7oWjcsPDQoXzX3iKB8/ebQlnds2i8oL9y3Kig358LjNzpiK4T0V4VL+gVb/e7sj2BSJpEihyB8xHpnviyA5KvzrKyf3b9paDzdcabKJI496FYR3XHM7DwzO9rMjt73gFFSAhxnciaRIncizB04FFmIbihmxI7+Lu023w9/B2sZ3rqb8ttZsPBWv8NCTb1NNPVCPzy9Tj/8HfelTJEacRkhclnQ7ftSVk2/o6D5t4nG3SEaWPsTC9YJjwW1TU/KoqkjW8nxRFP3unJRUk09rq305l3cfwqBukmu2FcBOILkI/K8jpchy7N1nGVnyVdsM+uSfC2A7yBLO/x0nuvhOMvORPnYZvYtAN+aUlscZ2pU+jJtH8Q226Xttlmo27aJbmxKPE01+XxCY6disaqBI50o0bso3NcI26Nx8L7o0K7o1B29Xe3ZlojBz1vY9m2iodvS/9hm5bZSzR3ZSvoDei4pW1G3l+/OmiVx8GXW2I4zs7hjO7XEHdupJZVq7B4agXaMdWOYH9EQsaWaXHWo6tSihgWAnV3NDdFH5mF79bFwHJuV9T0Rh6Ize3K8hUL7VONqf2GbxJnnpTFtSc6Ylxh9XF9oq7g/Iho7slW5bRphdwmmtmqF7Q00u+YqLAG/Yju1xB3bqSXu2E4tqTaObQx0seq+rSLUVGNv7ZfHrbf1ZH03LHc036FTniuimlvzGzR229ftpf5eLywXY8Xbpe1be2Fa6pZ+mCa6dyPM325LrsjWRG7IXC+sL7KVxLGR0NhRf0RtJWWNY6vtJ8Wv2E4tccd2aok7tlNLKs4VaQTaTnXjXo0wjt1Eeb7Dln6oQ+ckp3lB4thd0YmazxBp6kRsVi8LqjPXv+6NmBWOk/JPrj0tKGvMv9NRjZ2IWyfK+lqKDl0a2M5zRRxnMO7YTi1xx3ZqSbW5IkZsLWjsKI4tmls19tZIU4fbz3flvcHEe3vsqMYu19xR/oPI0EhXzjA6jsgOGXMlbSuUl7U/okO2iO2KcW4dk2Qp+BXbqSXu2E4tccd2aknFuSKNINascehV/VCYNST/YYto6q2qsSXfYkdH3nGU3I2m5jeILlSdqDpSXjOEVWrNyZiT3JC4P5LIo1HNnbCdauyG9k+K27vGdpzBuGM7tcQd26klFb/zyCB+qnHoVQxzjJtMxbElH1t0o44jovnX48ZiGx0ZG0MuC9aabk7x7kRttSD9EdP+yM7x3nFUW8Uau0Szu8Z2nMG4Yzu1xB3bqSWVv/NYjJ9qHHqlCLWmxLGj3BCdN0XnmEnkFKsulGlgIp2o62ONjT2GuD9S/j6o9j9Stpokju25Io4zBHdsp5YkHZvkp0luJHljYdn+JC8leVv+/0G7t5mOMx6jqMLzAXwUwL8Xlp0J4Ptmdk4+cemZAN6aqqhvjUDbbelKHFuEWUMG7lBNrvkO2zuhxtbxqHWcj1S+Q0pH6nh0DZ1McYbZ1pFxRDqJ/OtUbojYqrmz3FZNmQOzV3Uc28x+BOBeWXwKgAvyzxcAOHXypjjO9Fiqxn6wmW0AgPz/wcM2LM7zuLB5+7DNHGeq7PbOY3Gex5X7rU7v4DhTYKmR17tJrjOzDSTXAdg4yk59Y6CxNddjZVfG7pOApmryKN9B55hRTb0zlUOssViJY0djPsu0yXtQrkiUf53U2JCy2iZhK5bPQTMrceyLAbwy//xKAN+YvCmOMz1GCfd9EcAVAB5N8k6SZwA4B8CJJG8DcGJedpyZISlFzOz0IaueM+W2OM7UqDxXpBg/1XFAVjbKNbZuP9eRsTE0p3hnqBvT+dZhWXViKh+7uQflimj+dV80dis1jkjq/dBEPnY0Jou/8+g4adyxnVriju3UksrnUi/GT1UjJzW2jhsyYf51nBuSiGMvhLkrmv8QzQ8+w8T51+Pmgkg5yl0XW1HqV4290995dJwk7thOLalUipgxCDPpq1zzrVBqRFNOJ4blitNUU6mX5a8zqfRo6O21bAiBGaezM7Qddbq7cR+pq61Utomtmi0dbq7QFpcijjMYd2ynlrhjO7WkYo0dhplUY7cl/NeSV8PmO+XDckVpquOGsFRzi6ZudkQ3ynR6e9IQZ9EwwWOm9Ea2itaLrbrlodFgf9fYjjMYd2ynlrhjO7Wk8jh28bG3xqHbzVBzq8bW7fWxsA59m3q9KUpL3VmuEymamzr8QnvPuU5MnqZa3h+JbKVDnEl/pFn47jyO7ThDcMd2aok7tlNLKo9jF6cy1jj0imYvKHdEbOnUG3GaavnQt9H0El2NxZbrxMZOrUCOtwdpbE3pZWj6Af0R1dTltkvZSvsvYRx7cpG9B30TjjM67thOLXHHdmpJtQMGGIN8jm4v/F3pEGUUjb0gceve2PnXkHJ5vkOjEwpPLogQben2ovlnmcRrc+np7lRTi63EdpBnEmxpXs8uW3oc23GG4I7t1BJ3bKeWVKyxEcSaNQ690JT8ZxFbHR3CTHOKNTarU7ilcopVJ0oslh2psK9x7D1HY6enLRHbdFN5NGKrBRXpmlcjw88VnxmEVS8Jv2I7tcQd26klo4yP/TCSPyB5M8mbSL4+X+5T4jkzyygauwvgTWZ2Dcm9AVxN8lIAr8K4U+IZgvhpvxfqrG43FFc90dj6np6+R9eIyuHh49is6kjRjRrHFo1tfc3/3nPGEU7n0Wg5tBVVY6dsJRqb0RDPu95nncabo6NMh7fBzK7JP28FcDOAQ+BT4jkzzFgam+ThAJ4M4KcYcUq84nR4vbn5CZvrOKMxsmOTXAvgawDeYGZbRt2vOB1ec+2apbTRccZmJFFIso3MqT9vZl/PF48/JZ4xiDX3RRN3m6JZRWz1JLckyr/WnOKEbmSkG6UCjVtLmaKxofvPMBrzT9ouyqORYHPKVqKxIf2RoL4q8rFJEsCnANxsZucWVvmUeM7MMsoV+2kAXgHgBpLX5svehmwKvAvz6fF+C+Clu6WFjrMERpkO73IMj8D4lHjOTFJ54JW9osYOlVBf3hnUX5OOlUctq4ZW3dgrz39IxWLRkaTkvk53IdvPMGqbyFb6Puiktmrq+6kl+duej+04g3HHdmqJO7ZTSyrPxw50sGhk1dz6zmOUfx3pxMT4c1FZ8h2krDrQulKBxFvZ3YPi2Gq7RP9EY/4Y11Y9sXVb3m/tVRzHdpw9EXdsp5a4Yzu1pHqNXdTBfZmXRDRynCyicevy8ecoceuoHOnGhE5MaGz09hyN3dD+SBTz17JoZD3XlK2oeT4Sxy5odh9XxHGG4I7t1BJ3bKeWVKqxiVA/adw5mjdRk0VEk8twcLHGlvWaK0LVjQmNbVE+drlGf95R7wjKCw9ZG5TnH7xrzp3tB4fntrC/jKmyr8ybuEpFcFiEjhMyH+ZqrOzpM4Hp9kfUVjo+tmr04BmCx7EdZzDu2E4tccd2asky5GMXCjpGW1/j2Lpe6xpTcyfi3OhLBRqrlbIqwUg3JnRqUfPruUS20fVa1sZofyRRf2S7yJYax9Z55RO2M837kfVq+wnxK7ZTS9yxnVriju3UkurHxy5KKdXUkcYWjaq6UXRlUidq3DnS1OW6UcvxO5lSn9Svxy+2L26r9h9Eo2qQP6o70f8Y01bj9kciW2juel/7L8U4NibGr9hOLXHHdmqJO7ZTS6qPYxdzRcplGCIVq6FO1di6fxTble3H1o3l+daluhEA9QRtyGfEtonOJWWrxPaRrRLPAOLcc+1/lNuKiT5BaHvPFXGcgbhjO7XEHdupJdVr7GJ8tUxzDkJ1WkpHqm6Mtk9obC2blhO6Map/eHtS/YNUnNuaox9rYP1ROXEuybJocB1LXDW42m5C/Irt1BJ3bKeWjDKjwSqSV5K8Lp/n8ex8uc/z6Mwso2jsBQDPNrO5fC6ay0n+J4AXY9x5HoVY14luFCEYby8VJjX6eDrRUusjTZ0IzJeVU+eWYsztk7ZLxLmT/ZPooUQi37rEFkthlHkezczm8mI7/zP4PI/ODDOSxibZzOef2QjgUjNb2jyP8z7Po1MNIzm2mfXM7CgAhwI4huQTRj1AMM/jGp/n0amGseLYZraZ5GUATsKS5nnEVPRTUF+BpAZfZkrbN3YMP/EMIBXzH5Morj0h2n+ZxvzpRUaJihxEcr/882oAJwC4BT7PozPDjHLFXgfgApJNZD+EC83sEpJXwOd5dGaUUeZ5vB7Akwcs3wSf59GZUfzJo1NL3LGdWuKO7dQSd2ynllSbj01MN2CpodpEebkpbV+qrVEQXNentp/MGDbh/gqnXJ/iV2ynlrhjO7XEHdupJZW/81gk1pzlOjKpoZM6NTGPpJRVB6rOZEP3l+vEOMebtH8w5vZJ22lZL4F6Lo3UuSeuoSW2WAp+xXZqiTu2U0vcsZ1aUrnGtkZxMA1ZOWlsVnWq/GxjXZnQiZFuTJQjza31D29PSvNG59JI5Ecnz718fRS3Hrss8zqqbZrhvJOm6yfEr9hOLXHHdmqJO7ZTS6rX2EVdGcVGtSw6MtKww+seWJ/qzmZCY4sOZFRWjR2uh6yPdevwtqZto+Ux+x9J22j9qpFFQydsFX13UblY3+R626/YTi1xx3ZqiTu2U0uqz8cu/pQ0FqvlSHfK2HnJWK+WU5padWP5+liDi6ZuaFlyTxqDP2dlG7ptdjC13Xj7j2ur8fsjiWcCZf0RzxVxnMG4Yzu1xB3bqSXVx7GL0kp/VgmNrbFPnXcl1oUlx8b4ujEV147yH6R+a4XlfmF9FDdW2+gcM1rW7XW46mj74Xp/YFn7D2P2P9S2bGmuyHSvsX7FdmqJO7ZTS9yxnVpSqcY2SK5IQjdqeoL1ErFZ1dCyvq+aV3Wj6D60QvOwLROMq45U3Sjr4+MP/jywHMX8wyKjuHUiFzzR30j1D/RcKbZCW8qau679kVYxju25Io4zkJEdO5+H5uckL8nLPh2eM7OMc8V+PYCbC+UzkU2HdwSA7+dlx5kJRtLYJA8F8AIA7wHwxnzxKQCOzz9fAOAypOZ5pOhq0YWMNLboxkj36XpIebq60XoJjR1pbo39ahx78GdgQG6IaupmX8o6B2V5zL8flYe3LSvruSQ0tpZVY2t/pGCraYy5OOoV+4MA3oIw7O/T4TkzyyiTK70QwEYzu3opB/Dp8JzlYBQp8jQAf0LyZACrAOxD8nNYynR4jlMRo0yudBaAswCA5PEA3mxmLyf5L8imwTsHo06HR8CKRxSd12j1dfOwLbI+zs0Q3ai5GXK2fdXArfJ8BvakgigWK3HrqP5w9yCOrTF36T9Azr3RSvRHZHs9t6g/UtK2bL32R8o1c6S51VYS5w6eKSxzHPscACeSvA3AiXnZcWaCcWfmvQxZ9MOnw3NmmmVIW911y6TeXhsiReSO1O9rSEpvv3os2V9DWnJ7tXb57ZJ9TQGQG57eXuV2HUujQtpqou2RbGuGoceGPlKXmFlP1yce4Ze1FRjBVj2RjfqIX8N9xfr81TDHGYw7tlNL3LGdWlLx8AsWPAZvSsiqpeE+DWGpxpbt+5omGoX3wrKGsFQ3RqmXfdWpEkKT/fvtco1d1NVxeK88FKq2Uo3d74cavBc9wk48UtdTb8t6ObeG2ko0dvTaXpmtPG3VcQbjju3UEndsp5ZUP8RZQdtpLLbVCstNqm4Mf4ddTVtV3ViiaYEBurGlulGDu+Vx7Fg3ph7pD2+bPhJXTZ2ylbXD+rpt6Y9ExxvetqysGrn83KHpB6qxV4Trg5g+Jsev2E4tccd2aok7tlNLlkFj79J6bRnOYKXoRo1j9yT/Yae+uhVp7vDwPdG8Wu6vkFjvChWaYVFTN3X/ceLYUW6InEtTNHZb+idNjWOLrTQO3lPNHcWpy3NFenJuTbWVlU8d0tf+yIpCfVO43PoV26kl7thOLXHHdmpJ5bkijYK2a8kQAitb3aDc0Nis6MYdosk7bdGZUTmhI6NyqAMbqhtl6NtIU0f1hbsXY8XaVooGbrdD26ySssaxlQXpj/Q05i9ti2P+Wi7vXzT6mo8tce8Vw21V5fALjrNH4Y7t1BJ3bKeWVKqxSaBZ0MUrRSeqxm5RYq8ivjQO3tH87CgHGbK+PMc40oFqrobGwcs1dhQ3L+hWa5fHrTXGvyphK2W72GohMTycauo45l9uuyjZRHSzxsEDW3k+tuMMxh3bqSXu2E4tqVhjW6CLVSeubnWCsurGrowDtkJ053aJ/apujTVv2D7VjVE+RLh59B5fsv6ovKt9mufSbGkezXi2UuZb4cEborn77eH50dl6SLm8f0GJ+Zvmisj2PY9jO04ad2ynlrhjO7Wk8jh2URerTlzbXgjKDXn7rS/B0Pn2inD7KMe4PB+iH+4exWqbK8p/9zr1Ri8Rxy7LFYG0fcUK0dTt0FZ7tXYG5ZTGjnJLJE7eXZGwVTImrwMtJqYpieLgxTg2Jsav2E4tccd2asmos4atB7AVQA9A18yOJrk/gC8DOBzAegCnmdl9u6eZjjMe42jsZ5nZPYXy4jyP55A8My+XTodHWpAfohp7TTPUjZqPrRo7ypeQ2O/C2PnZ8o7lyoTGTsWxRcNHOc+F42n+teaGRLYaV2O3Vgbltmj4KM8mmVsu55qyFcs1ejHmv9xx7FOQze+I/P+pE7fGcabEqI5tAL5L8mqSr86XjT3PY/f+bZO32HFGYFQp8jQzu4vkwQAuJXnLqAcws/MAnAcAa45cN43RqxwnyUiObWZ35f83krwIwDFYwjyPDVigizVuvaYlcWzV2KYaO9Sdmp+tOcdxbDYs9yS2KpIfeoPTKexiTZ3KtyiOFV7+/qfGrbU/0mrIdNiCanTNs9k2bsw/8QwgMo6GuaNck+HbLoVRZuZdQ3Lvxc8AngvgRgAXI5vfERh1nkfHqYhRrtgPBnARs15tC8AXzOzbJK8CcCHJMwD8FsBLd18zHWc8RpmZ99cAnjRguc/z6MwsleaKNGjYq71LG6qmXtscT2OrRtd8iPkVknMsY2skNbeUVfxFGnsMTQ0AVsjP0Liy5obsLbbap7U9KLcTGlttpe+bTp6fHZYp35XGprU/47kijjMC7thOLXHHdmpJ9Rq7EI9VTb2v6MZGIv9hc2t1UI7G/tP87CjnODV2nx5RNHKUK4LSss4LE4wVLnFlzQVJ9UfajfDcla3NVUFZ82yaorG7qbh24v1OHR97HFstd66I48ws7thOLXHHdmpJxRq7H2jsfVo7gvVrm2G5KZO+9GVexbWiQzX2G+dnjzfviuZDpOLY8Xt8sl5yV4o52BpX1twQ1dR7i63aLNfYqtGLzxOAAeMgjpvLnrBVPHZfia1cYzvOYNyxnVriju3UEprOq7I7D0b+HsDtAA4EcE9i8+ViltsG/P9o32FmdtAkFVTq2A8clPyZmR1d+YFHYJbbBnj7RsWliFNL3LGdWrJcjn3eMh13FGa5bYC3bySWRWM7zu7GpYhTS9yxnVpSqWOTPInkrSR/mY/3t6yQ/DTJjSRvLCzbn+SlJG/L/z9omdr2MJI/IHkzyZtIvn7G2reK5JUkr8vbd/Ysta8yxybZBPAxAM8H8DgAp5N8XFXHH8L5AE6SZYuDbR4B4Pt5eTnoAniTmT0WwHEA/ja316y0bwHAs83sSQCOAnASyeNmpn1mVskfgD8C8J1C+SwAZ1V1/JJ2HQ7gxkL5VgDr8s/rANy63G3M2/INACfOYvsA7AXgGgDHzkr7qpQihwC4o1C+M182a4w02GaVkDwcwJMB/BQz1D6STZLXIhve7lIzm5n2VenYg7JsPdaYgORaAF8D8AYz27Lc7SliZj0zOwrAoQCOIfmEZW7SA1Tp2HcCeFihfCiAuyo8/qjcnQ+yiVEH29xdkGwjc+rPm9nXZ619i5jZZgCXIeuvzET7qnTsqwAcQfIRJFcAeBmygS1njZkYbJPZYImfAnCzmZ1bWDUr7TuI5H7559UATgBwy6y0r+pOxskAfgHgVwDePgOdni8C2ACgg+yOcgaAA5D15m/L/++/TG17OjKpdj2Aa/O/k2eofU8E8PO8fTcCeGe+fCba54/UnVriTx6dWuKO7dQSd2ynlrhjO7XEHdupJe7YTi1xx3Zqyf8B8oinGC1gc8UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAL8AAAEICAYAAAAOUiV/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAARaElEQVR4nO3deZBc1XnG4d+rkQa0IoHYtCCEoZRAKhUoxYaAgTIhZgvg2K6AjQPErkDZGAxUHHAR4wCJndhgORVsIIAhBYEQILaD2ReBMTsCh0WyAbFoYEASICSEkDTDlz/uAVrDTM8wfaa7qfM+VSr17Xv73K9Pv3OX7j59FRGYlWhUqwswaxWH34rl8FuxHH4rlsNvxXL4rVgO/wAkzZf0lQbbOE/S39eZ/x1JlzWyjtwkHSXp7gztbCspJI3OUVdNuyFp+3T7EklnDbetIYdf0hckPSTpTUndkm6QtMdwV1yCiDg2Is4EkLS3pK5W12TvG1L4JZ0EzAP+CdgS2Ab4MXBIzmJU8d7ImmLQoEnaBDgD+FpEXBsRqyNifUT8b0T8bVpmI0nzJL2U/s2TtFGaN0XSdZKWSXo93Z5R0/58Sf8o6dfAW8B2ade7WNIqSc9K+mLN8n8taWFq6yZJs2rmhaRjJT2V5p8rSUN87L6SFkl6Q9K/Ae89rk9/bCxpjaSpafo0ST2SJqXpsyTNS7cvSdPjgRuAaWnP+aakaanJTkn/kZ7rE5Lm1nktQtJX0/NbJelMSR+TdK+klZKuktQ5xH4fsI/7rPP7ku6WtEn6d1Ha87+YnltHWq5D0g8kLZe0GDhwoOeRlp8p6dpU36upzwd9neq0NzU9xxWSXpP0q0E3pBFR9x+wH9ADjK6zzBnAfcAWwObAPcCZad5mwGeBccBE4L+Bn9U8dj7wArATMBrYBFgJzEnztwZ2SrcPBZ4Gfj8texpwT01bAVwHTKbaOy0D9hvsscDUtM7PAWOAE9Nz/soAz/cu4LPp9s3AM8D+NfM+k25fApyVbu8NdPVp5zvA28ABQAfwXeC+Ov0cwC+ASam/1gK3AdulfnsSOHKwfgfG1+njo4C7qTaM/w7cBIxL834GnJ8evwXwAHBMmncssAiYCWwK3JHq/UBu0nP9DfDD1NbGwB4f4jXevp/+/S5wXnr9xgCfBFQ320MI/xeBlwdZ5hnggJrpTwPPDbDsHwGv9wn/GTXT44EV6YUb2+exNwBfrpkeRbW3mFXTMXvUzL8KOGWwxwJ/VRs6qq1+FwOH/0zgX9OL8zJwAvC99CKuAaZ+iPDfWjO9I7BmkPDvXjP9MPB3NdNnA/MG6/dB+vgo4H7gv4BrgM50/5ZUf2xja5Y9HLgj3b4dOLZm3p8xcPh3o9ow9TdvKK9xf+E/A/j5u/OG8m8ox9evAlNV/6x9GvB8zfTz6T4kjZN0vqTnJa2k2jJOfnd3mSx590ZErAb+kmpL0i3pl5J+L82eBfwo7dpWAK9RBXV6TVsv19x+C5gwhMdO61ND1E73406qMO8CPAbcAuwF7Ao8HRHL6zy2r771bjxIX79Sc3tNP9MToH6/D9LHANtTnc/9Q0SsS/fNotqidtf04flUewDo04dsmIe+ZgLPR0RPP/OG8hr35/tUe4yb0+HcKYMsP6Tw30u1az60zjIvpaLftU26D+BkYA7wiYiYBOyZ7q89pt7gq6URcVNE7Eu1O15EtfuFqnOPiYjJNf/GRsQ9Q3ge9R7bTfWCVIVV5wkzB2qI6rBuDvAZ4M6IeDI95wOp/jD60+yvz9bt9zp9DLAQOBq4QdKcdN8Sqi3/1Jr+mxQRO6X5G/QhVX8MZAmwzQB/5MN6jSNiVUScHBHbAX8OnCRpn3qPGTT8EfEG8G3gXEmHpi3KGEn7S/qXtNgVwGmSNk8ngt8G3n3/eiLVFmmFpE2B0+utT9KWkg5OJ4lrgTeB3jT7POBUSTulZTeR9PnBnsMQHvtLYCdJf5FekOOBrer0yVtUhxxf4/2w3wMcw8DhfwXYTNUbCM0wYL8P0scARMQVwLeAWyV9LCK6qc5vzpY0SdKodLK9V3rIVcDxkmZImgLU2/I+QPXH8j1J41W9ibB7mjes11jSQZK2Txuulen59NZ7zJDeVoyIc4CTqE4+llH9dR5HdQIEcBbwEPB/VIcBC9J9UL1FOhZYTnVSfOMgqxtFtdV6iWqXtxfw1VTH/wD/DFyZduWPA/sP8TkM+Nh0mPJ5quP2V4EdgF8P0uSdVIcBD9RMT6Q6vOhv/YuoNhKL0y59Wn/LZTSPgft9wD7uU/OlVMfSt0valurcqJPqxPp14GqqPQe8f3L8G6rX/9qBCouIXqqt8/ZUb3Z0UR2GNfIa7wDcSvWHfC/w44iYX+8BSicLZsXxB0pWLIffiuXwW7EcfitW1q+bDmbqph0xa2a+VS5ZP2HwhT6EN1aOy9re+IlvZ2tru85V2doC+N3bed9xXftWZ7a2el57jd7Vq/v9blVOTQ3/rJmjue/GGYMvOEQndn8iW1sAN9w24HfKhuWP91iUra3/nH1HtrYAPr3woKztPbOg3meCH86L836Yra16fNhjxXL4rVgOvxXL4bdiOfxWrIbCL2k/Sb+V9PRQvj9t1k6GHf40GOVcqm/c7QgcLmnHXIWZjbRGtvwfpxq1tDiN9rmSzL/mYDaSGgn/dDYcttZFP0PNJP2Nqt/7eWj5q+80sDqzvBoJf38fP39gcEBEXBARcyNi7tTNfH5t7aORNHax4ZjNGbw/btes7TUS/geBHSTNVvVDSYdR/aaM2UfCsL/YFhE9ko6jGrfZAVwcEU9kq8xshDX0rc6IuB64PlMtZk3lM1ArlsNvxXL4rVgOvxWrqcMYV70j5r89Jlt7j76ab0gkwNileYeNXrnbBVnby+nZu87M2t7YV/L13aj12Zqqv57mrMas/Tj8ViyH34rl8FuxHH4rlsNvxXL4rVgOvxXL4bdiOfxWLIffiuXwW7EcfiuWw2/FcvitWA6/Fcvht2I5/FYsh9+K1dQxvCt7x3HjG3+Yrb0lL22arS2AaUvK+RXpd17Me83hiV35+q7DY3jNRpbDb8Vy+K1YDr8Vy+G3Yjn8VqxGLkU6U9IdkhZKekLSCTkLMxtpjbzP3wOcHBELJE0EHpZ0S0Q8mak2sxE17C1/RHRHxIJ0exWwkH4uRWrWrrIc80vaFtgZuL+fee9dh/etFWtzrM4si4bDL2kCcA3wjYhY2Xd+7XV4x03eqNHVmWXTUPgljaEK/uURcW2eksyao5F3ewRcBCyMiHPylWTWHI1s+XcHvgR8StKj6d8BmeoyG3GNXIT6biDvdXzMmsif8FqxHH4rlsNvxWrqMMY3ezq5d+nsbO2N6e7M1hbAhBdWZ22vnY3tzrvdm/B8vr4btbY5w0m95bdiOfxWLIffiuXwW7EcfiuWw2/FcvitWA6/Fcvht2I5/FYsh9+K5fBbsRx+K5bDb8Vy+K1YDr8Vy+G3Yjn8ViyH34rV1DG863s66F6+Sbb2JizN+7NBo5csz9peOxv3SmRtb0zGvtO6nmxt1eMtvxXL4bdiOfxWLIffiuXwW7EcfitWjssSdUh6RNJ1OQoya5YcW/4TqK7EaPaR0ug1uWYABwIX5inHrHka3fLPA74JDPizurWXIu1dVc6vIFv7a+SCdAcBSyPi4XrL1V6KtGPi+OGuziy7Ri9Id7Ck54ArqS5Md1mWqsyaYNjhj4hTI2JGRGwLHAbcHhFHZKvMbIT5fX4rVpavNEfEfGB+jrbMmsVbfiuWw2/FcvitWA6/FaupY3jpFe+sGpOtuc438o5D7V2adwzvrl84O1tb3XvnvTbtNiuyNkfvsnx9Fz0ew2s2ohx+K5bDb8Vy+K1YDr8Vy+G3Yjn8ViyH34rl8FuxHH4rlsNvxXL4rVgOvxXL4bdiOfxWLIffiuXwW7EcfiuWw2/Fau4Y3hBan+/auaMyD/WM9euytjdqfb4xxurJe81h9eQdExzrMvZd5B2bPRBv+a1YDr8Vy+G3Yjn8ViyH34rl8FuxGr0a42RJV0taJGmhpN1yFWY20hp9n/9HwI0R8TlJncC4DDWZNcWwwy9pErAncBRARKwD8n5KZDaCGjns2Q5YBvxU0iOSLpT0gWuNbnAd3jffbGB1Znk1Ev7RwC7ATyJiZ2A1cErfhTa4Du+ECQ2sziyvRsLfBXRFxP1p+mqqPwazj4RGrsP7MrBE0px01z7Ak1mqMmuCRt/t+TpweXqnZzFwdOMlmTVHQ+GPiEeBuXlKMWsuf8JrxXL4rVgOvxXL4bdiNXcMr4LozDc+s7czW1MAjNp446zt9W6Ub9xtjOnN1hbAOxlrAxg1dmy2trSmOdtkb/mtWA6/Fcvht2I5/FYsh9+K5fBbsRx+K5bDb8Vy+K1YDr8Vy+G3Yjn8ViyH34rl8FuxHH4rlsNvxXL4rVgOvxXL4bdiNXUMr0YHYya/na29tVPy/vDtqK22yNre25PzbVs6J7+VrS2ANVM+8IPaDZm45eb5GutqTiy95bdiOfxWLIffiuXwW7EcfitWo5ciPVHSE5Iel3SFpLw/eWY2goYdfknTgeOBuRHxB0AHcFiuwsxGWqOHPaOBsZJGU12D96XGSzJrjkauyfUi8APgBaAbeCMibu673AaXIl25eviVmmXWyGHPFOAQYDYwDRgv6Yi+y21wKdJJeT9VNGtEI4c9fwo8GxHLImI9cC3wJ3nKMht5jYT/BWBXSeMkiepSpAvzlGU28ho55r+f6sLTC4DHUlsXZKrLbMQ1einS04HTM9Vi1lT+hNeK5fBbsRx+K5bDb8Vq6jDGztE9zN78tWztLd4i74dm62ZulrW9NVvma2v25q/mawxYslXeIaDrttk0W1uxtCNbW/V4y2/FcvitWA6/Fcvht2I5/FYsh9+K5fBbsRx+K5bDb8Vy+K1YDr8Vy+G3Yjn8ViyH34rl8FuxHH4rlsNvxXL4rVgOvxWrqWN4J45eyyc3fzpbe89Mz3j5S2DVrI2ytrd2+rpsbe05NV+/AVw4bUbW9lZuk++6JL2PNWeb7C2/Fcvht2I5/FYsh9+K5fBbsRx+K9ag4Zd0saSlkh6vuW9TSbdIeir9P2VkyzTLbyhb/kuA/frcdwpwW0TsANyWps0+UgYNf0TcBfT9ddlDgEvT7UuBQ/OWZTbyhnvMv2VEdAOk/7cYaMHa6/Cufj3fJ55mjRrxE97a6/COn9I50qszG7Lhhv8VSVsDpP+X5ivJrDmGG/5fAEem20cCP89TjlnzDOWtziuAe4E5krokfRn4HrCvpKeAfdO02UfKoF9pjojDB5i1T+ZazJrKn/BasRx+K5bDb8Vy+K1YiojmrUxaBjw/hEWnAstHuJzhaufaoL3rG2ptsyIi7wDtfjQ1/EMl6aGImNvqOvrTzrVBe9fXbrX5sMeK5fBbsdo1/Be0uoA62rk2aO/62qq2tjzmN2uGdt3ym404h9+K1Vbhl7SfpN9KelpSW40LljRT0h2SFkp6QtIJra6pL0kdkh6RdF2ra+lL0mRJV0talPpwt5bX1C7H/JI6gN9RfUW6C3gQODwinmxpYUkatLN1RCyQNBF4GDi0XeoDkHQSMBeYFBEHtbqeWpIuBX4VERdK6gTGRcSKVtbUTlv+jwNPR8TiiFgHXEk1UL4tRER3RCxIt1cBC4Hpra3qfZJmAAcCF7a6lr4kTQL2BC4CiIh1rQ4+tFf4pwNLaqa7aKNw1ZK0LbAzcH+LS6k1D/gm8E6L6+jPdsAy4KfpsOxCSeNbXVQ7hV/93Ncex2Q1JE0ArgG+ERErW10PgKSDgKUR8XCraxnAaGAX4CcRsTOwmjb4rad2Cn8XMLNmegbwUotq6ZekMVTBvzwirm11PTV2Bw6W9BzV4eKnJF3W2pI20AV0RcS7e8qrqf4YWqqdwv8gsIOk2emE6DCqgfJtQZKojlkXRsQ5ra6nVkScGhEzImJbqn67PSKOaHFZ74mIl4Elkuaku/YBWv5GQVMvS1RPRPRIOg64CegALo6IJ1pcVq3dgS8Bj0l6NN33rYi4vnUlfaR8Hbg8bdgWA0e3uJ72eavTrNna6bDHrKkcfiuWw2/FcvitWA6/Fcvht2I5/Fas/wfAmbX10Zf7cwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(data)\n",
    "plt.gca().set_title(\"Original with masked cells\")\n",
    "plt.show()\n",
    "plt.imshow(bavg)\n",
    "plt.gca().set_title(\"Coarsened with masked cells\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many ways to develop this further, including using non-square kernels, but this is my use case."
   ]
  }
 ],
 "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}