grid masking

gridmasking.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:8478dbaf749c2a991d706fa3e2c08923ebb0eda9c0478b20a905552c4e8f3214"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from matplotlib.path import Path\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def makeQuadCoords(xarr, yarr, zpnt=None):\n",
      "    coords = np.vstack([\n",
      "        np.hstack([xarr[0,:], xarr[1,::-1]]),\n",
      "        np.hstack([yarr[0,:], yarr[1,::-1]])\n",
      "    ]).T\n",
      "\n",
      "    return coords\n",
      "\n",
      "def drawCells(xx, yy, ax, mask=None):\n",
      "    rows, cols = xx.shape\n",
      "    for ii in range(rows - 1):\n",
      "        for jj in range(cols - 1):\n",
      "            if mask is not None and mask[jj, ii]:\n",
      "                coords = None\n",
      "                \n",
      "            else:\n",
      "                coords = makeQuadCoords(\n",
      "                    xarr=xx[jj:jj+2, ii:ii+2],\n",
      "                    yarr=yy[jj:jj+2, ii:ii+2],\n",
      "                )\n",
      "\n",
      "            if coords is not None:\n",
      "                rect = plt.Polygon(coords, edgecolor='w', linewidth=1.25,\n",
      "                                   zorder=0, facecolor='0.875')\n",
      "                ax.add_artist(rect)\n",
      "                \n",
      "def points_inside_poly(points, verts):\n",
      "    poly = Path(verts)\n",
      "    return [ind for ind, p in enumerate(points) if poly.contains_point(p)]\n",
      "\n",
      "def mask_polygon(xc, yc, polyverts, mask=None, mask_value=0.0):\n",
      "    \"\"\"\n",
      "    Mask Cartesian points contained within the polygon defined by polyverts\n",
      "    A cell is masked if the cell center (x_rho, y_rho) is within the\n",
      "    polygon. Other sub-masks (mask_u, mask_v, and mask_psi) are updated\n",
      "    automatically.\n",
      "    mask_value [=0.0] may be specified to alter the value of the mask set\n",
      "    within the polygon.  E.g., mask_value=1 for water points.\n",
      "    \"\"\"\n",
      "    \n",
      "    if mask is None:\n",
      "        mask = np.zeros_like(xc)\n",
      "\n",
      "    polyverts = np.asarray(polyverts)\n",
      "    if polyverts.ndim != 2:\n",
      "        raise ValueError('polyverts must be a 2D array, or a '\n",
      "                         'similar sequence')\n",
      "\n",
      "    if polyverts.shape[1] != 2:\n",
      "        raise ValueError('polyverts must be two columns of points')\n",
      "\n",
      "    if polyverts.shape[0] < 3:\n",
      "        raise ValueError('polyverts must contain at least 3 points')\n",
      "\n",
      "    inside = points_inside_poly(\n",
      "        np.vstack([xc.flatten(), yc.flatten()]).T,\n",
      "        polyverts\n",
      "    )\n",
      "    if np.any(inside):\n",
      "        mask.flat[inside] = mask_value\n",
      "        \n",
      "    return mask\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def drawDemo(usemask=False):\n",
      "\n",
      "    xn, yn = np.meshgrid(np.arange(1, 10), np.arange(1, 10))\n",
      "    xc, yc = np.meshgrid(np.arange(1.5, 9), np.arange(1.5, 9))\n",
      "    xp, yp = [3.25, 3.25, 6.25, 6.25, 3.25], [3.25, 7.25, 7.25, 3.25, 3.25]\n",
      "    mask = mask_polygon(xc, yc, np.vstack([xp, yp]).T, mask_value=1)\n",
      "    \n",
      "    fig, ax = plt.subplots()\n",
      "    ax.plot(xn.flatten(), yn.flatten(), 'ko')\n",
      "    ax.plot(xc.flatten(), yc.flatten(), 'b.')\n",
      "    ax.plot(xp, yp, 'g-')\n",
      "    drawCells(xn, yn, ax, mask=mask if usemask else None)\n",
      "    \n",
      "    ax.set_xlim(0, 10)\n",
      "    ax.set_ylim(0, 10)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "drawDemo(usemask=False)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEACAYAAABF+UbAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF+RJREFUeJzt3X9s3PV9x/HXezjcNUMBNZMgHYkSWWNjVbMOsQmxVTkx\nlkNpxfrHtA1lsYWq7Z/NdjetGkZacSJN6f6oSs9jf3QDZDuU/UgrBDFKw3q9ZBJTCuVHKQljs5ws\nZkkgbAuCDZcs7/1xd47t2F/7vr77fl6f770eUiTH2HdP3Zl3vr67773N3SEiInH6idABIiKSnoa4\niEjENMRFRCKmIS4iEjENcRGRiGmIi4hELHGIm9ljZnbezF6b97mPm9lzZvammR0xsxs6nykiIktZ\n6Uj8cQD3LPrcAwCec/dbAHy38XcREQnAVjrZx8y2AnjG3T/V+PsbAHa4+3kzuwlAzd1/rtOhIiJy\ntTSPid/o7ucbH58HcGMbe0REpAVremLT64fxOm9fRCSQnhTfc97MbnL3c2a2CcDbS32RmWm4i4ik\n4O622q9NcyT+NID+xsf9AJ5KCNEfdzz00EPBG1j+6LbQbaHbIvlPq1Z6ieGTAJ4H8LNmdsbM7gfw\nFQC/bmZvArir8XcREQkg8eEUd79vmf90dwdaRESkRTpjMwOlUil0Ag3dFlfotrhCt0V6K75OPPUF\nm3mnLltEJK/MDN7hJzZFRISEhriISMQ0xEVEIqYhLiISMQ1xEZGIaYiLiERMQ1xEJGIa4iIiEdMQ\nFxGJmIa4iEjENMRFRCKmIS4iEjENcRGRiOV2iE9OTqJcLqNUKqFcLmNycrIrG1g61MDVwdDA0sHQ\nsCYdXDHkoRw6dMh7e3ubS5wdgPf29vqhQ4e6qoGlQw1cHQwNLB0MDYs1ZufqZ20rX9zSBQcc4jt3\n7lxwpzT/lMvlrmpg6VADVwdDA0sHQ8NirQ7xXD6cMjs7u+TnP/zww65qYOlQA1cHQwNLB0PDWiXu\n2IxVoVBY9r/NzMxk0uAJW42yamDpUANXB0MDS8dyDcViMZPrb4dcDvHBwUFMTU1hampq7nNbtmxB\nX18fLl26lElDf38/Tp8+jdOnTwdrYOlQA1cHQwNLx1INvb29GBgYyOT62yGXQ/yzn/0sAGB0dHTu\n16K+vj7cddddmTU0r2tiYmLuX/usG1g61MDVwdDA0tG8rvHxcQD1I/CBgYG5GRKDrliUPDMzk+kR\nxnyFQgGbNm0K2sDSoQauDoYGlo6enh7cfPPNQa57MS1KFhHpIhriIiIR0xAXEYmYhriISMQ0xEVE\nIpbLlxi2y/DwRkxPr0OxeBmVygVs2HC5azsYGlg6GBpYOhgamDpC0JF4gunpdTh+vIijR9djeHhj\nV3cwNLB0MDSwdDA0MHWEoCGeoFis/2u+ffss9u9/t6s7GBpYOhgaWDoYGpg6QtAQT1CpXMCuXR9g\nYuJ80F/PGDoYGlg6GBpYOhgamDpC0BmbHcZwNhpLhxq4OhgaWDp0xqaIiAShIS4iEjENcRGRiKUe\n4mY2bGavm9lrZvZNM1t+E4OIiHREqic2zWwrgCqAW9191sz+DsCz7j4272uCPrE5OTmJSqWC2dlZ\nuDv6+/szf7/karU69z7FoRpYOtTA1cHQwNJRrVYxNjYGM0OhUMDg4GDQ9xNv9YnNtGdsvgfgIwDr\nzez/AKwH8FbKy2q7yclJDA0NLdjs09zckdUPSLVaxb59+xZsDMm6gaVDDVwdDA0sHUs1NOdGLIsh\nUr/E0Mx+H8BXAfwvgO+4+55F/z3YkXi5XMaRI0eu+vyOHTtw4MCBTBp2796NY8eOBW1g6VADVwdD\nA0vHcg3lchmHDx/OpGGxTI7EzawXwBcBbAVwEcA/mNlud39i/teNjIzMfVwqlVAqldJcXcuW22AN\nILPXgpotfx9k+XpUhg41cHUwNLB0LNeQ5bb7Wq2GWq2W+vvTPpxyO4Dn3f1dADCzbwO4E8CyQzxL\ny227NzOcPXs24xq+BpYONXB1MDSwdGS57X7xAe7evXtb+v60Q/wNAH9mZh8D8CGAuwF8P+Vltd1y\n2+737NmTeJTeTn19fTh16tRVm7yzbGDpUANXB0MDS8dSDV2x7d7dXzWzcQAvArgM4CUA32hn2Fpo\n2z1Xhxq4OhgaWDq07T7pgvXeKQA43heCpUMNXB0MDSwdeu8UEREJQkNcRCRiGuIiIhHTEBcRiZiG\nuIhIxLTtPgHLBm2GDoYGlg6GBpYOhgamjhB0JJ6AZYM2QwdDA0sHQwNLB0MDU0cIGuIJWDZoM3Qw\nNLB0MDSwdDA0MHWEoCGegGWDNkMHQwNLB0MDSwdDA1NHCDpjs8MYzkZj6VADVwdDA0uHztgUEZEg\nNMRFRCKmlxjKnE8/+Wlc/PHF0BnScP211+OV+14JnSHkNMRlzsUfX8R0/3THLp/hsU+GhtV2bBvb\nlnGVxEgPp4iIRCy3R+KTk5OoVCqYnZ2Fu6O/vz/zN72vVqtzbzYfqoGlQw1cHQwNLB3VahVjY2Mw\nMxQKBQwODka1FCKXQ3xychJDQ0ML1rM11y9l9QNSrVaxb9++BWufsm5g6VADVwdDA0vHUg3NuRHL\nIM/lEK9UKgsGOFD/4RgfH8fOnTszaRgbG1vwgxGiIU1HT0/7fyQYbguGhjQdeb0/WDqWapiamsLo\n6KiGeEhJS1azekG/2fKv1c/ypIJWOzrRxnBbMDSk6cjr/cHSsVxDczdvDHI5xAuFwpKfNzOcPXs2\n4xq+hqSOLNsYbguGhqSObrs/WDqKxWLQ629FLof44OAgpqamFjyksmXLFuzZsyfxKL2d+vr6cOrU\nqQW/qmXdkKajE20MtwVDQ5qOvN4fLB1LNfT29mJgYCCT62+HXA7x5mNZo6Ojc78W9fX1ZfqkTfO6\nJiYm0HwPmawbWDrUwNXB0MDS0byu5itkisUiBgYGonk8HNAbYHVcbCeX6GQfno5uuD9YOvQGWCIi\nEoSGuIhIxDTERUQipiEuIhIxDXERkYjl8iWG7TI8vBHT0+tQLF5GpXIh2O4+hg6GBpYOhgaWDoYG\npo4QdCSeYHp6HY4fL+Lo0fUYHt7Y1R0MDSwdDA0sHQwNTB0haIgnKBbr/5pv3z6L/fvf7eoOhgaW\nDoYGlg6GBqaOEDTEE1QqF7Br1weYmDgf9Nczhg6GBpYOhgaWDoYGpo4Q9Jh4gg0bLuORR94JnUHR\nwdDA0sHQwNLB0MDUEYKOxEVEIqYhLiISsdRD3MxuMLODZnbSzE6Y2R3tDBMRkZWt5THxrwN41t1/\n08x6APxkm5pERGSVUg1xM7sewGfcvR8A3P0SgIvtDFsrbbvn6lADVwdDA0tHt2673wbgHTN7HMAv\nAPgBgCF3/5+2la2Btt1zdaiBq4OhgaUjD9vuUy2FMLPbAfwzgDvd/QUzexjAe+7+5XlfE2wpRLlc\nxpEjR676/I4dO3DgwIFMGnbv3o1jx44FbWi1Y/Ojm3HmC2eCNnQKQ0OrHXm+P1g6lmsol8s4fPhw\nJg2LtboUIu2R+AyAGXd/ofH3gwAeWPxFIyMjcx+XSiWUSqWUV9cabbtP35HX7eoMDWk68np/sHQw\nbLuv1Wqo1Wqpvz/VEHf3c2Z2xsxucfc3AdwN4PXFXzd/iGdJ2+7Td3TbdnWGhqSObrs/WDqy3Ha/\n+AB37969LX3/Wl6dMgDgCTO7FsAUgPvXcFltpW336Tvyul2doSFNR17vD5aOrt527+6vAvilNra0\njbbdc3WogauDoYGlQ9vuky5Y2+4BcGzyXm1HN2xXZ2hYbUc33B8sHdp2LyIiQWiIi4hETENcRCRi\nGuIiIhHTEBcRiZg2+yRg2aDN0MHQwNLB0MDSwdDA1BGCjsQTsGzQZuhgaGDpYGhg6WBoYOoIQUM8\nAcsGbYYOhgaWDoYGlg6GBqaOEDTEE7Bs0GboYGhg6WBoYOlgaGDqCEGPiSdg2aDN0MHQwNLB0MDS\nwdDA1BGCjsRFRCKmIS4iEjENcRGRiGmIi4hELLdPbGrbPVeHGrg6GBpYOrp12z01bbvn6lADVwdD\nA0tHHrbd53KIVyqVBQMcqP9wjI+PY+fOnZk0jI2NLfjBCNGQpqOnp/0/Egy3BUNDmo683h8sHUs1\nTE1NYXR0VEM8JG27T9+R1+3qDA1pOvJ6f7B0MGy7X6tcDnFtu0/f0W3b1Rkakjq67f5g6chy2/1a\n5XKIa9t9+o68bldnaEjTkdf7g6Wjq7fdM9O2e64ONXB1MDSwdGjbfdIFa9s9AI5N3qvt6Ibt6gwN\nq+3ohvuDpUPb7kVEJAgNcRGRiGmIi4hETENcRCRiGuIiIhHL5UsM24VlgzZDB0MDSwdDA0sHQwNT\nRwg6Ek/AskGboYOhgaWDoYGlg6GBqSMEDfEELBu0GToYGlg6GBpYOhgamDpC0BBPwLJBm6GDoYGl\ng6GBpYOhgakjBD0mnoBlgzZDB0MDSwdDA0sHQwNTRwg6EhcRiZiGuIhIxDTERUQitqYhbmbXmNnL\nZvZMu4JERGT11nokPgTgBACO95wVEekyqV+dYmY3A9gF4M8B/HHbitpkcnISlUoFs7OzcHf09/dn\n/qb31Wp17s3mQzWwdKiBq4OhgaWjWq1ibGwMZoZCoYDBwcGolkKs5SWGXwPwJQAb2tTSNpOTkxga\nGlqwnq25fimrH5BqtYp9+/YtWPuUdQNLhxq4OhgaWDqWamjOjVgGeaohbmafA/C2u79sZqX2Jq1d\npVJZMMCB+g/H+Pg4du7cmUnD2NjYgh+MEA1pOnp62n/qAMNtwdCQpiOv9wdLx1INU1NTGB0dzfcQ\nB3AngHvNbBeAIoANZjbu7n3zv2hkZGTu41KphFKplPLqWpO0ZDWrFUxmy29XynINVKsdnWhjuC0Y\nGtJ05PX+YOlYrqG5mzcLtVoNtVot9fenGuLu/iCABwHAzHYA+JPFAxxYOMSzVCgUlvy8meHs2bMZ\n1/A1JHVk2cZwWzA0JHV02/3B0lEsFjO7rsUHuHv37m3p+9v1uxrVq1MGBwcxNTW14CGVLVu2YM+e\nPYlH6e3U19eHU6dOLfhVLeuGNB2daGO4LRga0nTk9f5g6Viqobe3FwMDA5lcfzuseYi7+1EAR9vQ\n0jbNx7JGR0fnfi3q6+vL9Emb5nVNTEzA3YM0sHSogauDoYGlo3ldzVfIFItFDAwMRPN4OABY88Zr\n+wWbeacuu1UzMzO4dOlSkOsuFArYtGlT0IbVdmwb24bp/umgDZ3G0LDajm64P1g6enp6Mn0+IImZ\nwd2Xf8JgEZ12LyISMQ1xEZGIaYiLiERMQ1xEJGIa4iIiEdN6tgTDwxsxPb0OxeJlVCoXgu3uY+hg\naGDpYGhg6WBoYOoIQUfiCaan1+H48SKOHl2P4eGNXd3B0MDSwdDA0sHQwNQRgoZ4gmKx/q/59u2z\n2L//3a7uYGhg6WBoYOlgaGDqCEFDPEGlcgG7dn2AiYnzQX89Y+hgaGDpYGhg6WBoYOoIQY+JJ9iw\n4TIeeeSd0BkUHQwNLB0MDSwdDA1MHSHoSFxEJGI6EpcFto1tC50gDddfe33oBImAhrjM6eSbLQEc\nb3TE0MDUIfHTwykiIhHL7ZG4tt1zdaiBq4OhgaWjm7fd09K2e64ONXB1MDSwdORh230ul0KUy2Uc\nOXLkqs/v2LEDBw4cyKRh9+7dOHbsWNAGlg41cHUwNLB0LNdQLpdx+PDhTBoWa3UpRC6PxLXtnqtD\nDVwdDA0sHQzb7tcql0Nc2+7j6FADVwdDA0tHltvu1yqXQ1zb7rk61MDVwdDA0qFt96S07Z6rQw1c\nHQwNLB3adp90wdp2D4DnpA6GDjVwdTA0sHRo272IiAShIS4iEjENcRGRiGmIi4hETENcRCRiuXyJ\nYbuwbNBm6GBoYOlgaGDpYGhg6ghBR+IJWDZoM3QwNLB0MDSwdDA0MHWEoCGegGWDNkMHQwNLB0MD\nSwdDA1NHCBriCVg2aDN0MDSwdDA0sHQwNDB1hKAzNjuM4Ww0lg41cHUwNLB06IxNEREJQkNcRCRi\nGuIiIhFLNcTNbLOZfc/MXjezH5nZYLvDRERkZame2DSzmwDc5O6vmNl1AH4A4PPufnLe1wR9YlPb\n7rk61MDVwdDA0sG27T6THZvufg7AucbH75vZSQCfAHAy8Rszom33XB1q4OpgaGDp0LZ7AGa2FcBR\nAJ909/fnfV7b7rVRXA2EHQwNLB1dv+2+8VDKQQBD8wd408jIyNzHpVIJpVJpLVe3atp2z9WhBq4O\nhgaWDoZt97VaDbVaLfX3px7iZrYOwLcAHHD3p5b6mvlDPEvadh9Hhxq4OhgaWDqy3Ha/+AB37969\nLX1/qiFu9X++HgVwwt0fTnMZnaRt91wdauDqYGhg6cjDtvu0r075VQDHAPwQQPMCht398LyvCf7q\nlJDb7oH6kyahN4qzdKiBq4OhgaVj/itkGLbdt/qYuN47pcMY3heCpUMNXB0MDSwdeu8UEREJQkNc\nRCRiGuIiIhHTEBcRiZiGuIhIxLTtPgHLBm2GDoYGlg6GBpYOhgamjhB0JJ6AZYM2QwdDA0sHQwNL\nB0MDU0cIGuIJWDZoM3QwNLB0MDSwdDA0MHWEoCGegGWDNkMHQwNLB0MDSwdDA1NHCDpjs8MYzkZj\n6VADVwdDA0uHztgUEZEgNMRFRCKmIS4iEjENcRGRiOX2ZB9tu+fqUANXB0MDSwfbtvtW5XKIa9s9\nV4cauDoYGlg6tO0+6YK17T74Jm+WDjVwdTA0sHR0/bZ7Vtp2z9WhBq4OhgaWDoZt92uVyyc2l9t2\nn+UGa4YGlg41cHUwNLB0MDSsmbt35E/9osM4dOiQ9/b2OupLnB2A9/b2+qFDh7qqgaVDDVwdDA0s\nHQwNizVm56pnbS4fTmk+IdHcdh9igzVDA0uHGrg6GBpYOhga1iqXT2yKiMRK750iItJFNMRFRCKm\nIS4iEjENcRGRiGmIi4hETENcRCRiGuIiIhHTEBcRiZiGuIhIxDTERUQipiEuIhIxDXERkYhpiIuI\nRCz1EDeze8zsDTP7VzP703ZGiYjI6qQa4mZ2DYC/BHAPgJ8HcJ+Z3drOsDyp1WqhE2jotrhCt8UV\nui3SS3sk/ssA/s3dT7n7RwD+FsBvtC8rX/QDeoVuiyt0W1yh2yK9tEP8pwGcmff3mcbnREQkQ2mH\nuFb2iIgQSLWezczuADDi7vc0/j4M4LK7/8W8r9GgFxFJoZX1bGmHeA+AfwHwawD+A8D3Adzn7idb\nvjAREUkt1bZ7d79kZn8I4DsArgHwqAa4iEj2OrbtXkREOq8jZ2zqRKA6M9tsZt8zs9fN7EdmNhi6\nKSQzu8bMXjazZ0K3hGRmN5jZQTM7aWYnGs8xdSUzG278//GamX3TzAqhm7JiZo+Z2Xkze23e5z5u\nZs+Z2ZtmdsTMbljpcto+xHUi0AIfAfgjd/8kgDsA/EEX3xYAMATgBPTqpq8DeNbdbwWwHUBXPhRp\nZlsB/B6A29z9U6g/NPs7IZsy9jjqc3K+BwA85+63APhu4++JOnEkrhOBGtz9nLu/0vj4fdT/Z/1E\n2KowzOxmALsA/A2AVT/znjdmdj2Az7j7Y0D9+SV3vxg4K5T3UD/QWd94scR6AG+FTcqOu/8TgP9a\n9Ol7AYw1Ph4D8PmVLqcTQ1wnAi2hcdTxiwCOhy0J5msAvgTgcuiQwLYBeMfMHjezl8zsr81sfeio\nENz9PwF8FcC/o/4qt/92938MWxXcje5+vvHxeQA3rvQNnRji3f6r8lXM7DoABwEMNY7Iu4qZfQ7A\n2+7+Mrr4KLyhB8BtAP7K3W8D8AFW8StzHplZL4AvAtiK+m+o15nZ7qBRRLz+qpMV52knhvhbADbP\n+/tm1I/Gu5KZrQPwLQAH3P2p0D2B3AngXjObBvAkgLvMbDxwUygzAGbc/YXG3w+iPtS70e0Annf3\nd939EoBvo/6z0s3Om9lNAGBmmwC8vdI3dGKIvwjgZ8xsq5ldC+C3ATzdgeuhZ2YG4FEAJ9z94dA9\nobj7g+6+2d23of7EVdXd+0J3heDu5wCcMbNbGp+6G8DrAZNCegPAHWb2scb/K3ej/sR3N3saQH/j\n434AKx74pTrZJ4lOBFrgVwD8LoAfmtnLjc8Nu/vhgE0Muv0htwEATzQOcqYA3B+4Jwh3f7XxG9mL\nqD9X8hKAb4Styo6ZPQlgB4CfMrMzAL4M4CsA/t7MvgDgFIDfWvFydLKPiEi8tJ5NRCRiGuIiIhHT\nEBcRiZiGuIhIxDTERUQipiEuIhIxDXERkYhpiIuIROz/ASv2axKzHIALAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x105a70b38>"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "drawDemo(usemask=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEACAYAAABF+UbAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGIlJREFUeJzt3X9s3PV9x/HXeyT46lWAyiSgI1GQNTZWNeuqbEJsVU6M\ncShUrH+MbSjLWaga/2y2u2nVMNKKHWmi+6MqPU/7g40ixyGsW1qhNoY0Wa/nTGJKoYWWkrBslpPF\nXQiFbUGw4cLy3h9359iO72v7e3ffz+tz39dDinCMfzzl+94757v73tvcHSIiEqefCh0gIiLpaYiL\niERMQ1xEJGIa4iIiEdMQFxGJmIa4iEjEEoe4mX3ZzM6b2ctL3vchMztqZqfM7IiZXdP9TBERWc1a\nt8SfAHDXivc9COCou98M4FuNv4uISAC21sk+ZrYNwDfc/aONv78KYKe7nzez6wHU3P0Xuh0qIiKX\nS3Of+HXufr7x9nkA13WwR0RENqCtBza9fjNe5+2LiASyKcXnnDez6939NTO7AcDrq32QmWm4i4ik\n4O623o9Nc0v86wAGG28PAng6IUR/3PHwww8Hb2D5o5+Ffhb6WST/2ai1nmL4FIDnAPy8mZ01s/sB\nfB7Ab5rZKQC3N/4uIiIBJN6d4u73tfhfd3ShRURENkhnbGagWCyGTqChn8Ul+llcop9Fems+Tzz1\nFzbzbn1tEZFeZWbwLj+wKSIiJDTERUQipiEuIhIxDXERkYhpiIuIRExDXEQkYhriIiIR0xAXEYmY\nhriISMQ0xEVEIqYhLiISMQ1xEZGIaYiLiESsZ4f49PQ0SqUSisUiSqUSpqenc9nA0qEGrg6GBpYO\nhoa2dHHFkIdy6NAhHxgYaC5xdgA+MDDghw4dylUDS4cauDoYGlg6GBpWaszO9c/ajXzwhr5wwCF+\n5513LrtQmn9KpVKuGlg61MDVwdDA0sHQsNJGh3hP3p2ysLCw6vvffffdXDWwdKiBq4OhgaWDoaFd\niTs2Y9XX19fy/83Pz2fS4AlbjbJqYOlQA1cHQwNLR6uGQqGQyffvhJ4c4sPDw5idncXs7Ozi+7Zu\n3YpyuYz3338/k4bBwUGcOXMGZ86cCdbA0qEGrg6GBpaO1RoGBgYwNDSUyffvhJ4c4nfffTcAYGJi\nYvHXonK5jNtvvz2zhub3mpqaWvzXPusGlg41cHUwNLB0NL/Xvn37ANRvgQ8NDS3OkBjkYlHy/Px8\nprcwlurr68MNN9wQtIGlQw1cHQwNLB2bNm3CjTfeGOR7r6RFySIiOaIhLiISMQ1xEZGIaYiLiERM\nQ1xEJGI9+RTDThkdvRZzc5tRKFxEpfIGrrrqYm47GBpYOhgaWDoYGpg6QtAt8QRzc5tx/HgBMzP9\nGB29NtcdDA0sHQwNLB0MDUwdIWiIJygU6v+ab9++gEceeTPXHQwNLB0MDSwdDA1MHSFoiCeoVN7A\nrl3vYGrqfNBfzxg6GBpYOhgaWDoYGpg6QtAZm13GcDYaS4cauDoYGlg6dMamiIgEoSEuIhIxDXER\nkYilHuJmNmpmr5jZy2Z2wMxab2IQEZGuSPXAppltA1AFcIu7L5jZVwA84+6TSz4m6AOb09PTqFQq\nWFhYgLtjcHAw89dLrlari69THKqBpUMNXB0MDSwd1WoVk5OTMDP09fVheHg46OuJb/SBzbRnbL4F\n4D0A/Wb2fwD6Afwo5dfquOnpaYyMjCzb7NPc3JHVAVKtVrF3795lG0OybmDpUANXB0MDS8dqDc25\nEctiiNRPMTSzBwB8AcD/Avimu+9Z8f+D3RIvlUo4cuTIZe/fuXMn9u/fn0nD7t27cezYsaANLB1q\n4OpgaGDpaNVQKpVw+PDhTBpWyuSWuJkNAPgMgG0ALgD4BzPb7e5PLv24sbGxxbeLxSKKxWKab7dh\nrTZYA8jsuaBmrS+DLJ+PytChBq4OhgaWjlYNWW67r9VqqNVqqT8/7d0pOwA85+5vAoCZfQ3AbQBa\nDvEstdp2b2Y4d+5cxjV8DSwdauDqYGhg6chy2/3KG7jj4+Mb+vy0Q/xVAH9uZh8A8C6AOwB8J+XX\n6rhW2+737NmTeCu9k8rlMk6fPn3ZJu8sG1g61MDVwdDA0rFaQy623bv7981sH4AXAFwE8D0Aj3Uy\nrB3ads/VoQauDoYGlg5tu0/6wnrtFAAcrwvB0qEGrg6GBpYOvXaKiIgEoSEuIhIxDXERkYhpiIuI\nRExDXEQkYtp2n4BlgzZDB0MDSwdDA0sHQwNTRwi6JZ6AZYM2QwdDA0sHQwNLB0MDU0cIGuIJWDZo\nM3QwNLB0MDSwdDA0MHWEoCGegGWDNkMHQwNLB0MDSwdDA1NHCDpjs8sYzkZj6VADVwdDA0uHztgU\nEZEgNMRFRCKmpxjKoo899TFc+MmF0BnScPWVV+Ol+14KnSHkNMRl0YWfXMDc4FzXvj7DfZ8MDevt\nuGnypoyrJEa6O0VEJGI9e0t8enoalUoFCwsLcHcMDg5m/qL31Wp18cXmQzWwdKiBq4OhgaWjWq1i\ncnISZoa+vj4MDw9HtRSiJ4f49PQ0RkZGlq1na65fyuoAqVar2Lt377K1T1k3sHSogauDoYGlY7WG\n5tyIZZD35BCvVCrLBjhQPzj27duHO++8M5OGycnJZQdGiIY0HZs2df6QYPhZMDSk6ejVy4OlY7WG\n2dlZTExMaIiHlLRkNasn9Ju1fq5+licVbLSjG20MPwuGhjQdvXp5sHS0amju5o1BTw7xvr6+Vd9v\nZjh37lzGNXwNQH0hbFZaXR55a0jqaHVcZHmssBybDB1ZHxft6MkhPjw8jNnZ2WV3qWzduhV79uxJ\nvJXeSeVyGadPn172q1rWDa06BgYGMDQ0lFnDapdHHhtadSQdF904VpiPTYbraYjjoh09OcSb92VN\nTEws/lpULpczfdCm+b2mpqbQfA2ZrBuWdjSfAVAoFDA0NJTp/X0rL4+8NqzWAejYDNnBcP1ol14A\nq8tYTi5Zzwv82LjBH+a4zPJirZN9ev3kK5YOvQCWiIgEoSEuIhIxDXERkYhpiIuIRExDXEQkYj35\nFMNOGR29FnNzm1EoXESl8kaw3X0MHQ88AJw6BfT3AwcOANdck3kCTQdDA8BxXDA0MHWEoFviCebm\nNuP48QJmZvoxOnptrjtOnQJmZoBnn60PsVAYOhgaAI7jgqGBqSMEDfEEhUL9X/Pt2xfwyCNv5rqj\nv7/+3x07gMceC5JA08HQAHAcFwwNTB0haIgnqFTewK5d72Bq6nzQX88YOg4cAO69Fzh6NNzdBywd\nDA0Ax3HB0MDUEYLO2OwyhrPRAJ2xyUpnbHJ06IxNEREJQkNcRCRiqYe4mV1jZgfN7KSZnTCzWzsZ\nJiIia2vneeJfAvCMu/+2mW0C8NMdahIRkXVKNcTN7GoAn3D3QQBw9/cBXOhkWLu07X55R+ht3ksv\njzw3rOzQsRm+g+H60Y60t8RvAvBjM3sCwC8B+C6AEXf/n46VtUHb7pM7st7mvdrlkceGVh06NrXt\nvh2pnmJoZjsA/DOA29z9eTN7FMBb7v65JR8T7CmGpVIJR44cuez9O3fuxP79+zNp2L17N44dOxa0\nIamjVCrh8OHDy97XracYtro8VmvoFoaGpI7Vjostj2/B2U+f7XgD+7HJcD3N+rhYaqNPMUx7S3we\nwLy7P9/4+0EAD678oLGxscW3i8UiisViym+3Mdp2v3ZHltu8W10eeWtI6gC07T5EB8P1o1aroVar\npf78VEPc3V8zs7NmdrO7nwJwB4BXVn7c0iGeJW27X1veNs0zNCR1aNs9V0eWx8XKG7jj4+Mb+vx2\nnp0yBOBJM7sSwCyA+9v4Wh2lbffJHXncNM/Q0KpD2+617b4dqYe4u38fwK90sKVjtFH88g5tuw/f\nsFoHoGNT2+7bo9dO6TKG14UA9NoprPTaKRwdeu0UEREJQkNcRCRiGuIiIhHTEBcRiZiGuIhIxLTt\nPgHLBm2GDpYN7wwdDA0Ax3HB0MDUEYJuiSdg2aDN0MGy4Z2hg6EB4DguGBqYOkLQEE/AskGboYNl\nwztDB0MDwHFcMDQwdYSgIZ6AZYM2QwfLhneGDoYGgOO4YGhg6ghBZ2x2GcPZaIDO2GSlMzY5OnTG\npoiIBKEhLiISMQ1xEZGIaYiLiESsZ0/20Ubx5R2ht3kzbJpnaFjZoWMzfAfD9aMdPTnEtVE8uSOP\nm+YZGlp16NjUtvt29ORTDLXtfu0ObbvPviGpQ9vuw3Tkeds9NW27X7sjb5vmGRqSOgBtuw/RwXD9\naFdPDnFtu19b3jbNMzQkdWjbPVdH1sdFO3pyiGvbfXJHHjfNMzS06tC2e227b0dPDnFtFL+8Q9vu\nwzes1gHo2NS2+/b05AObK+m1U/TaKaz02ikcHXrtFBERCUJDXEQkYhriIiIR0xAXEYmYhriISMR6\n8imGncKyQZuhg2XDO0MHQwPAcVwwNDB1hKBb4glYNmgzdLBseGfoYGgAOI4LhgamjhA0xBOwbNBm\n6GDZ8M7QwdAAcBwXDA1MHSFoiCdg2aDN0MGy4Z2hg6EB4DguGBqYOkLQGZtdxnA2GqAzNlnpjE2O\nDp2xKSIiQWiIi4hETENcRCRibQ1xM7vCzF40s290KkhERNav3VviIwBOANCjYSIiAaQ+Y9PMbgSw\nC8BfAPiTjhV1yPT0NCqVChYWFuDuGBwczPxF76vV6uKLzYdqaHZMTk7CzNDX14fh4eHMX/R+6eWR\n54aVHTo2w3cwXD/a0c5p918E8FkAV3WopWOmp6cxMjKybAVWc/1SVgdItVrF3r17l619yrqhVUfz\n55LVgbra5ZHHhlYdOjbDdTBcP9qV6u4UM/skgNfd/UUA634+Y1YqlcqyKwlQPzj27duHTZs2ZfJn\ncnJy2YERoqFVx+zsLCYmJrp+OTStdnnksaFVR6vjAkDujk2G62mI46IdaW+J3wbgHjPbBaAA4Coz\n2+fu5aUfNDY2tvh2sVhEsVhM+e02JmnJalZP6Ddr/W9blicVtOpo7nfMQqvLI28NSR3A6sdFN44V\n9mMzyw6G60etVkOtVkv9+amGuLs/BOAhADCznQD+dOUAB5YP8Sz19fWt+n4zw7lz5zKu4WsA6gth\ns9Lq8shbQ1JHq+Miy2OF5dhk6MjyuFh5A3d8fHxDn9+pl6KlenbK8PAwZmdnl/3aunXrVuzZsyfx\nllAnlctlnD59etmvalk3tOoYGBjA0NBQZg2rXR55bGjVkXRcdONYYT42Ga6nIY6LdrQ9xN19BsBM\nB1o6pvmAxMTExOKvReVyOdMHbZrfa2pqCs3XkMm6YWlH8xkAhUIBQ0NDmT5os/LyyGvDah2Ajs2Q\nHQzXj3bpBbC6jOHFfQC9ABYrvQAWR8d6rh9Z0QtgiYjkiIa4iEjENMRFRCKmIS4iEjENcRGRiHXq\neeI9aXT0WszNbUahcBGVyhvBdvcxdDzwQH3Le39/fcdkqN2SDB0MDQDHccHQwNQRgm6JJ5ib24zj\nxwuYmenH6Oi1ue44dQqYmQGefbY+xEJh6GBoADiOC4YGpo4QNMQTFAr1f823b1/AI4+8meuO/v76\nf3fsAB57LEgCTQdDA8BxXDA0MHWEoCGeoFJ5A7t2vYOpqfNBfz1j6DhwALj3XuDo0XB3H7B0MDQA\nHMcFQwNTRwg6Y7PLGM5GA3TGJiudscnRoTM2RUQkCD07RZaxcbodH7l19ZVXh06QCGiIyyLdlZK9\n0HdlSPx0d4qISMR69pa4Noov7wi9zZth0zxDw8oOHZvhOxiuH+3oySGujeLJHXncNM/Q0KpDx6a2\n3bejJ59iWCqVcOTIkcvev3PnTuzfvz+Tht27d+PYsWNBG5I6SqUSDh8+nElDq8sjbw1JHTo2w3Qw\nXD9W2uhTDHvylri23a/dkbdN8wwNSR2Ajs0QHQzXj3b15BDXtvv1mZ+fz+T7JP1GlqeGpA6G44Kh\ngaUjy2337erJIa5t92t3lMvlzJ7aNjg4iDNnzuS+IalDx2aYDm27J6WN4lwdauDqYGhg6dC2+6Qv\nrNdOAcDxuhAsHWrg6mBoYOnQa6eIiEgQGuIiIhHTEBcRiZiGuIhIxDTERUQi1pNPMewUlg3aDB0M\nDSwdDA0sHQwNTB0h6JZ4ApYN2gwdDA0sHQwNLB0MDUwdIWiIJ2DZoM3QwdDA0sHQwNLB0MDUEYKG\neAKWDdoMHQwNLB0MDSwdDA1MHSHojM0uYzgbjaVDDVwdDA0sHTpjU0REgtAQFxGJmIa4iEjEUg1x\nM9tiZt82s1fM7IdmNtzpMBERWVuqBzbN7HoA17v7S2b2QQDfBfApdz+55GOCPrCpjeJcHWrg6mBo\nYOlg23afyY5Nd38NwGuNt982s5MAPgzgZOInZkQbxbk61MDVwdDA0qFt9wDMbBuAGQAfcfe3l7xf\n2+61UVwNhB0MDSwdud9237gr5SCAkaUDvGlsbGzx7WKxiGKx2M63WzdtFOfqUANXB0MDSwfDtvta\nrYZarZb681MPcTPbDOCrAPa7+9OrfczSIZ4lbbuPo0MNXB0MDSwdWW67X3kDd3x8fEOfn2qIW/2f\nr8cBnHD3R9N8jW7StnuuDjVwdTA0sHT0wrb7tM9O+XUAxwD8AEDzC4y6++ElHxP82Skht90D9QdN\nQm8UZ+lQA1cHQwNLx9JnyDBsu9/ofeJ67ZQuY3hdCJYONXB1MDSwdOi1U0REJAgNcRGRiGmIi4hE\nTENcRCRiGuIiIhHTtvsELBu0GToYGlg6GBpYOhgamDpC0C3xBCwbtBk6GBpYOhgaWDoYGpg6QtAQ\nT8CyQZuhg6GBpYOhgaWDoYGpIwQN8QQsG7QZOhgaWDoYGlg6GBqYOkLQGZtdxnA2GkuHGrg6GBpY\nOnTGpoiIBKEhLiISMQ1xEZGIaYiLiESsZ0/20bZ7rg41cHUwNLB0sG2736ieHOLads/VoQauDoYG\nlg5tu0/6wtp2H3yTN0uHGrg6GBpYOnK/7Z6Vtt1zdaiBq4OhgaWDYdt9u3rygc1W2+6z3GDN0MDS\noQauDoYGlg6Ghra5e1f+1L90GIcOHfKBgQFHfYmzA/CBgQE/dOhQrhpYOtTA1cHQwNLB0LBSY3au\ne9b25N0pzQckmtvuQ2ywZmhg6VADVwdDA0sHQ0O7evKBTRGRWOm1U0REckRDXEQkYhriIiIR0xAX\nEYmYhriISMQ0xEVEIqYhLiISMQ1xEZGIaYiLiERMQ1xEJGIa4iIiEdMQFxGJmIa4iEjEUg9xM7vL\nzF41s381sz/rZJSIiKxPqiFuZlcA+CsAdwH4RQD3mdktnQzrJbVaLXQCDf0sLtHP4hL9LNJLe0v8\nVwH8m7ufdvf3APwdgN/qXFZv0QF6iX4Wl+hncYl+FumlHeI/C+Dskr/PN94nIiIZSjvEtbJHRIRA\nqvVsZnYrgDF3v6vx91EAF939L5d8jAa9iEgKG1nPlnaIbwLwLwB+A8B/APgOgPvc/eSGv5iIiKSW\natu9u79vZn8E4JsArgDwuAa4iEj2urbtXkREuq8rZ2zqRKA6M9tiZt82s1fM7IdmNhy6KSQzu8LM\nXjSzb4RuCcnMrjGzg2Z20sxONB5jyiUzG21cP142swNm1he6KStm9mUzO29mLy9534fM7KiZnTKz\nI2Z2zVpfp+NDXCcCLfMegD92948AuBXAH+b4ZwEAIwBOQM9u+hKAZ9z9FgDbAeTyrkgz2wbgDwB8\n3N0/ivpds78XsiljT6A+J5d6EMBRd78ZwLcaf0/UjVviOhGowd1fc/eXGm+/jfqV9cNhq8IwsxsB\n7ALwtwDW/ch7rzGzqwF8wt2/DNQfX3L3C4GzQnkL9Rs6/Y0nS/QD+FHYpOy4+z8B+K8V774HwGTj\n7UkAn1rr63RjiOtEoFU0bnX8MoDjYUuC+SKAzwK4GDoksJsA/NjMnjCz75nZ35hZf+ioENz9PwF8\nAcC/o/4st/92938MWxXcde5+vvH2eQDXrfUJ3Rjief9V+TJm9kEABwGMNG6R54qZfRLA6+7+InJ8\nK7xhE4CPA/hrd/84gHewjl+Ze5GZDQD4DIBtqP+G+kEz2x00iojXn3Wy5jztxhD/EYAtS/6+BfVb\n47lkZpsBfBXAfnd/OnRPILcBuMfM5gA8BeB2M9sXuCmUeQDz7v584+8HUR/qebQDwHPu/qa7vw/g\na6gfK3l23syuBwAzuwHA62t9QjeG+AsAfs7MtpnZlQB+F8DXu/B96JmZAXgcwAl3fzR0Tyju/pC7\nb3H3m1B/4Krq7uXQXSG4+2sAzprZzY133QHglYBJIb0K4FYz+0DjunIH6g9859nXAQw23h4EsOYN\nv1Qn+yTRiUDL/BqA3wfwAzN7sfG+UXc/HLCJQd7vchsC8GTjRs4sgPsD9wTh7t9v/Eb2AuqPlXwP\nwGNhq7JjZk8B2AngZ8zsLIDPAfg8gL83s08DOA3gd9b8OjrZR0QkXlrPJiISMQ1xEZGIaYiLiERM\nQ1xEJGIa4iIiEdMQFxGJmIa4iEjENMRFRCL2/8b3ybkvdnC3AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x105934208>"
       ]
      }
     ],
     "prompt_number": 11
    }
   ],
   "metadata": {}
  }
 ]
}