larrybradley/region_photometry.ipynb

## region_photometry.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "80ad0a0b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-10-14T20:19:39.297205Z",
     "start_time": "2021-10-14T20:19:38.178282Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x13d0d9550>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD6CAYAAABuxZF5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVgklEQVR4nO3dW4xdd3XH8e+aMzffguOQWFYcERARiIcSpFEAwQNNCkopInlACIoqV7Lkl1YKAgmSVqqE1IfQBy4PVSurQfiBklAuShTRgmuCKqQq4JAACSFNSIOI5cQQYnydy5mz+nC26ez1/3v2njPnzJzx//eRRjN7n9s6+8yaPeu//xdzd0SkPBObHYCIbA4lv0ihlPwihVLyixRKyS9SKCW/SKEm29zJzF4AzgLLQNfd58xsD/AAcCPwAvBhd391NGGKyLBZm+v8VfLPuftvV+z7B+B37n6vmd0NXO3un17teaZtxmfZsc6QReRy5jnPoi9Ym/u2OvNfxh3Ae6qfjwDfB1ZN/ll28Ha7bR0vKSKredSPtb5v25rfge+a2WNmdqjat9fdT1Y/vwTszT3QzA6Z2XEzO77EQuvARGS02p753+3uJ8zsOuComf1i5Y3u7maWrR/c/TBwGOAq26O+xCJjotWZ391PVN9PAd8CbgFeNrN9ANX3U6MKUkSGrzH5zWyHme269DPwPuBJ4CHgQHW3A8CDowpSRIavzb/9e4Fvmdml+/+ru/+Hmf0I+JqZHQR+BXx4dGGKyLA1Jr+7Pw+8NbP/FUBN9yJblHr4iRRKyS9SKCW/SKGU/CKFUvKLFErJL1IoJb9IoZT8IoVS8osUSskvUiglv0ihlPwihVLyixRKyS9SKCW/SKGU/CKFUvKLFErJL1IoJb9IoZT8IoVS8osUSskvUiglv0ihlPwihVLyixRKyS9SKCW/SKGU/CKFUvKLFErJL1KoxiW6ZQsy2+wI8tw3OwJZQWd+kUIp+UUKpeQXKVTr5Dezjpk9bmYPV9uvN7NHzew5M3vAzKZHF6aIDNtazvx3AU+v2P4s8Hl3fyPwKnBwmIFJxWyAr4kx/RrgvcjItEp+M9sP/BnwL9W2AbcCX6/ucgS4cwTxiciItD3zfwH4FNCrtq8BTrt7t9p+Ebg+90AzO2Rmx83s+BIL64lVRIaoMfnN7APAKXd/bJAXcPfD7j7n7nNTzAzyFCIyAm06+bwL+KCZvR+YBa4CvgjsNrPJ6uy/HzgxujCvUIPUtDYRNls8h7X9B2+F+Ly9ATroeC/s6ISbM88ZH9PmGKnz0EAafyvc/R533+/uNwIfAb7n7h8DHgE+VN3tAPDgyKIUkaFbz3X+TwOfMLPn6LcB3DeckERkI6ypb7+7fx/4fvXz88Atww9JRDaCBvZspDb1a1NNH+v3TM1v8XUmJuIdmuMYRKy9e71wc/12s1zNX48taRdI2hFI34/aAFpR916RQin5RQql5BcplJJfpFBq8Bulpoa1TOebpIGvU+8YkzTmhdsBLO6L27mOQWttBMw1qsXGueXltW0DHhr8jPp9vJc5XzV1DFIDYJbO/CKFUvKLFErJL1Io1fzDFGvNtXbYAawT9sWaf2py1dsBmKzfx8J20ukH8u0Aq8kNyomderrd2raF7VzNz1K4Tzymy5lOPk0dgwiPURsAoDO/SLGU/CKFUvKLFEo1/6BaXBdvqvGT+h5gaqp+n+n6dlLPT2cmTQ7tAj4Z2gXiNuCxHSCGFi+l9zK1d7dew1vYjvW8Ly4mTxH7KPjiUv126tsAHl7GJnqr3p797ApsB9CZX6RQSn6RQin5RQqlmr+tASbiaKrxs/V6qPFtZmbV2302fQ6fqd+nN1v/mHuT6d98j/vi2w0lsXXTmn8i7JuYD9f5F0L9Pp/poxBq/LTvROZzCG0HaY2vCUFydOYXKZSSX6RQSn6RQin5RQqlBr9BtZiII+nEEzrwxMY7SBv4fFvY3l7fXt6RLoG2vG0ybNfjWJ5JY+9NhsExsf0vtJFNdNMGsc5C/U6di1Orb5/PTERyITSKhtuzzXBxVuDYESg08LWaEKQAOvOLFErJL1IoJb9IoVTzX05Dp57s6rgNE2fGQTpJBx4yNf7ObbXt7s56p56lq9J2g6Wd9ddd3FGPdXk2eQjL03GVn3CHUBJ3FjM1/3z9QdPn6/eZOle/faqTHsPJ3EQjK2Q/lVDzexx0lLQJZCYObRr8cwV2+tGZX6RQSn6RQin5RQqlmr+thkE7kC6okUy2GSfSzFznj9fxY42/uLu+vbA7jWPhNfV9i1fVb+9uT+vX3nSom5uu86fzcDB5of7+u2fq2zNhDJK3mDR0MtTvlpvAM04UOrX64iCe+exaDf65wujML1IoJb9IoRqT38xmzeyHZvYTM3vKzD5T7X+9mT1qZs+Z2QNmlhmcLiLjqk3NvwDc6u7nzGwK+IGZ/TvwCeDz7n6/mf0zcBD4pxHGuqmS6/oDXeev/33MTcQR++rH6/ixxr94Tfr3e/HqsL27Xr/2doWFMYDOtvq+iTAJZi/0h1+8mP7qdM+GSUPCGILcJCKR9UJfiOVwjX4pvUYfJwqNi4Ukn0vmmr3FhT8y64lcaRo/De87V21OVV8O3Ap8vdp/BLhzFAGKyGi0qvnNrGNmTwCngKPAL4HT7n7pT+yLwPUjiVBERqJV8rv7srvfDOwHbgHe3PYFzOyQmR03s+NLLAwWpYgM3Zpa+939NPAI8E5gt5ldKvL2Aycu85jD7j7n7nNTpH3ZRWRzNDb4mdm1wJK7nzazbcB7gc/S/yPwIeB+4ADw4CgDHancIJ5cR5DaQ9LHxNVmkk49cSWdmbSTT5yIIw7SSTrwhMY9gIVr661Vk9dcrG1f95rzyWOu2Xahtj09UW80W+zV43rl4vbkOV75/Y56HFPbwj3q72Wimx7jTug81Jmvv26cEATSWYGTVY06zasDJ02AyeefWx14aw/2adPavw84YmYd+v8pfM3dHzaznwP3m9nfA48D940wThEZssbkd/efAm/L7H+efv0vIluQeviJFEoDe9qKNWBu0onYqWdy9dVy40o6kE62GSfiiIN0YgceSGv81137am37TVedSh5zw+zvatvbw8idC716h6Rfz+9JnuOZ6etq278Kty8u1dsJJhbSYxgHB02dDx2FMsfMLjQc99gWk/vskhr/yu/lozO/SKGU/CKFUvKLFEo1/6ByfQOSwT/hb2us+TMDXeKCGnGyzTgRR26QTryOH2v8uV3/mzzmDdP1+1xl9d6YZ7zeQWvvVL2+zzm3WG8neOl8fbt7NtPPIbzfeDxyx6wz2VDTx8+lzYrLBdCZX6RQSn6RQin5RQql5BcplBr8hik2JE3ElW/rf2s903gVV8uNK+nEWXbjDDyQDtKJHXhi4x7ATZPnatu7Juq/Gmd79dtzXp59TW37hW31jkC/2bartt3LzF6cvN+4enDmmMXjmsy6pAa+LJ35RQql5BcplJJfpFCq+dtqsbpM83OE7cxTxpVy4mPi7XGWXUgn4oiDdGIHHkhr/J0TcSnf+cbniK8T40hmBM6dehreb3aZ3mGcwuLne+WP69GZX6RUSn6RQin5RQqlmr+t3hAma4zleeYp42q48THx9riSDqSTbcaJOOIgHchdx58Pt9fr9zO+M3mO+Doxjhhr8l6h8f3mjllubs01G8bnu8XozC9SKCW/SKGU/CKFUs0/THERh1BHWq9enFo3LVYnuqHv/mJ9O1xKz66WGxfUiJNttpmII53Mo17jP7+YPkd8nRjHcoh1OrwXyLzfcDxyxywe16R+3+KLa4yKzvwihVLyixRKyS9SKCW/SKHU4DeoXCNSbGiKDVHd+miRiUzjVWehvq8zX//7HFe06Z7NNPiF1XLjSjo5cSKOgVbsOVN/nRjHRIg1vheAznxo8AzHI3fM4nFNjrsaALN05hcplJJfpFBKfpFCqeZvy5MRNel9luu1p3frg2Es1vzz6eSbnYv1SS2nz9fr0+6ZMMHlTPr3e2FqW207rpYbV9KBdLLNOBFHHKQTO/BAWuN3X6nHMXO6Huv0meQpkvfbuRhq/swxi8c1Hvf4uWQ/u/j5FkBnfpFCKflFCtWY/GZ2g5k9YmY/N7OnzOyuav8eMztqZs9W368efbgiMixtav4u8El3/7GZ7QIeM7OjwF8Cx9z9XjO7G7gb+PToQh2h7HXfWAPWV4L13GNCbWmx9lwKbQALS8lTxJp/6lz97/NMKNdzq9bGWBeX6vV5XC0X0gU1ksk2w0QccZAOpNfxkxr/1fr9Z36f1tlT5+rHsHOx+ZjF40o47h7bYtpc549tAFdg34DGM7+7n3T3H1c/nwWeBq4H7gCOVHc7Atw5ohhFZATW1NpvZjcCbwMeBfa6+8nqppeAvZd5zCHgEMAsaQuxiGyO1g1+ZrYT+AbwcXevXaTx/v9R2f+L3P2wu8+5+9wU6dxxIrI5WiW/mU3RT/yvuPs3q90vm9m+6vZ9QLr6o4iMrcZ/+83MgPuAp939cytuegg4ANxbfX9wJBGOCY+z8lhzg1/S6WexPljG5usNcwCd8/V9U5240m/zykET3frf9ImF+nb3bLo6blwxN04KHGfRzc3CEwfqxE48sYFv5nSmwe9MvUGvc74+o5DNpy8cj2vT55CbqTd+viVoU/O/C/gL4Gdm9kS172/oJ/3XzOwg/U5kHx5JhCIyEo3J7+4/IL9CGsBtww1HRDaKeviJFEoDe9pKOn2k/wx53Bc79XRCjb+YdlixC2HyjonV/z5bL63fO6EEjrX4clyAF1ieDrHHlw1vP86yC+lEHHGQTuzAE+t7gMlzoV3kQlgNOHPMkk494bgnnXpyg3g0sEdESqHkFymUkl+kUKr5LyfWiRaut2euCxvhenJ8TKxXLW03aLqKPxlX/VnO1d71j3XqfBiUk5kApDcZ+xOEuEJJHFfSgczkoxfjdr0Wj9fwIa3x7WJ92xfSxyTHtam/RZtr+lfgQJ5IZ36RQin5RQql5BcplGr+QWWuC3vsEL8c6nOaa/5YacZ7WHzOpVDfkk4I0putf8y5CUA87osvHJtAsisMrz7ZZpyII9dPP17HT2r83HX+pfCYcIySGr/Aa/o5OvOLFErJL1IoJb9IoZT8IoVSg19bDZ1++vdZffCPx7a5OAlF7nXidsMqQJBpWLtQH1DUmUwnEfE4gKhhYI/lVr3pxtmLQ2xxwE3u/cdBOrGBbylt8IsNfOnn0KKBr4BOPZHO/CKFUvKLFErJL1Io1fyDytWIDYN/LKyCk7QBkHYE8jiQZ6phRVqAyfrHamGbzAQh1mJi0Jrc4JgQaxJb0+q5pBNxJINyYn0PSU3fOHCnwPo+R2d+kUIp+UUKpeQXKZRq/mFKasmGGj+38EfT6rCxTo6TggLWCdfK431y9X2u38JqcnVzrLWTer1hMQ0yk222mYij6Tq+avwsnflFCqXkFymUkl+kUKr5R6lpPECLCUHipKBu4e91rp6NtXS8rr/W+r6tGEu87p/c3ly/DzQRh2r8VnTmFymUkl+kUEp+kUIp+UUKpQa/jTTAhCBNHYMst1pw3BEbCduIHYHarHKTBLJ649xAHXayT6QGvkHozC9SKCW/SKEak9/MvmRmp8zsyRX79pjZUTN7tvp+9WjDFJFha3Pm/zJwe9h3N3DM3W8CjlXbslbuzV/JY3q1L19eTr96Xv/K3afpa6lb/xrkORriiO8lW+8Pcoyklcbkd/f/An4Xdt8BHKl+PgLcOdywRGTUBm3t3+vuJ6ufXwL2Xu6OZnYIOAQwy/YBX05Ehm3dDX7e77B92f+93P2wu8+5+9wUM+t9OREZkkHP/C+b2T53P2lm+4BTwwxKVhioph3TVWhVn4+VQc/8DwEHqp8PAA8OJxwR2ShtLvV9Ffhv4E1m9qKZHQTuBd5rZs8Cf1Jti8gW0vhvv7t/9DI33TbkWERkA6mHn0ihNLDnSqSGNWlBZ36RQin5RQql5BcplJJfpFBKfpFCKflFCqXkFymUkl+kUEp+kUIp+UUKpeQXKZSSX6RQSn6RQin5RQql5BcplJJfpFBKfpFCKflFCqXkFymUkl+kUEp+kUIp+UUKpeQXKZSSX6RQSn6RQin5RQql5BcplJJfpFBKfpFCKflFCqXkFynUupLfzG43s2fM7Dkzu3tYQYnI6A2c/GbWAf4R+FPgLcBHzewtwwpMREZrPWf+W4Dn3P15d18E7gfuGE5YIjJq60n+64Ffr9h+sdonIlvA5KhfwMwOAYcAZtk+6pcTkZbWc+Y/AdywYnt/ta/G3Q+7+5y7z00xs46XE5FhMncf7IFmk8D/ALfRT/ofAX/u7k+t8pjfAL8CXgv8dqAX3nhbJdatEidsnVi3Spzw/7G+zt2vbfOAgf/td/eumf018B2gA3xptcSvHnMtgJkdd/e5QV97I22VWLdKnLB1Yt0qccJgsa6r5nf3bwPfXs9ziMjmUA8/kUJtVvIf3qTXHcRWiXWrxAlbJ9atEicMEOvADX4isrXp336RQin5RQq1ock/zqMAzexLZnbKzJ5csW+PmR01s2er71dvZoyXmNkNZvaImf3czJ4ys7uq/WMVr5nNmtkPzewnVZyfqfa/3swerX4PHjCz6c2McyUz65jZ42b2cLU9lrGa2Qtm9jMze8LMjlf71vT5b1jyb4FRgF8Gbg/77gaOuftNwLFqexx0gU+6+1uAdwB/VR3LcYt3AbjV3d8K3AzcbmbvAD4LfN7d3wi8ChzcvBATdwFPr9ge51j/2N1vXnF9f22fv7tvyBfwTuA7K7bvAe7ZqNdvGeONwJMrtp8B9lU/7wOe2ewYLxP3g8B7xzleYDvwY+Dt9HuiTeZ+LzY5xv1V0twKPAzYGMf6AvDasG9Nn/9G/tu/FUcB7nX3k9XPLwF7NzOYHDO7EXgb8ChjGG/1b/QTwCngKPBL4LS7d6u7jNPvwReATwG9avsaxjdWB75rZo9Vg+dgjZ//yEf1XSnc3c1srK6LmtlO4BvAx939jJn94bZxidfdl4GbzWw38C3gzZsbUZ6ZfQA45e6Pmdl7NjmcNt7t7ifM7DrgqJn9YuWNbT7/jTzztxoFOGZeNrN9ANX3U5sczx+Y2RT9xP+Ku3+z2j228br7aeAR+v86764GhsH4/B68C/igmb1Af2KaW4EvMp6x4u4nqu+n6P9RvYU1fv4bmfw/Am6qWk+ngY8AD23g6w/iIeBA9fMB+rX1prP+Kf4+4Gl3/9yKm8YqXjO7tjrjY2bb6LdLPE3/j8CHqrttepwA7n6Pu+939xvp/25+z90/xhjGamY7zGzXpZ+B9wFPstbPf4MbKd5PfxjwL4G/3exGkxDbV4GTwBL92u4g/ZrvGPAs8J/Ans2Os4r13fRrvp8CT1Rf7x+3eIE/Ah6v4nwS+Ltq/xuAHwLPAf8GzGz2MQ1xvwd4eFxjrWL6SfX11KVcWuvnr+69IoVSDz+RQin5RQql5BcplJJfpFBKfpFCKflFCqXkFynU/wFY9mnrhblKRQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from astropy.modeling.models import Gaussian2D\n",
    "\n",
    "gm = Gaussian2D(100, 25, 25, 5, 5)\n",
    "y, x = np.mgrid[0:51, 0:51]\n",
    "data = gm(x, y)\n",
    "error = np.sqrt(data)\n",
    "plt.imshow(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d809b6e4",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-10-14T20:19:39.353463Z",
     "start_time": "2021-10-14T20:19:39.299008Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6161.286706670988"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# region \"exact\" photometry\n",
    "from regions import PixCoord, CirclePixelRegion\n",
    "region = CirclePixelRegion(PixCoord(25, 25), radius=5)\n",
    "mask = region.to_mask(mode='exact')\n",
    "weighted_flux = mask.get_values(data, mask=None)\n",
    "region_phot = np.sum(weighted_flux)\n",
    "region_phot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8ab7ce08",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-10-14T20:19:39.359357Z",
     "start_time": "2021-10-14T20:19:39.355297Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "78.49386413389895"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# with 1-sigma flux error\n",
    "weighted_var = mask.get_values(error**2, mask=None)\n",
    "region_err = np.sqrt(np.sum(weighted_var))\n",
    "region_err"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f126c1a4",
   "metadata": {},
   "source": [
    "### Equivalent photometry in photutils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5df5fa02",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-10-14T20:19:39.436477Z",
     "start_time": "2021-10-14T20:19:39.361813Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><i>QTable length=1</i>\n",
       "<table id=\"table5320055920\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>id</th><th>xcenter</th><th>ycenter</th><th>aperture_sum</th><th>aperture_sum_err</th></tr></thead>\n",
       "<thead><tr><th></th><th>pix</th><th>pix</th><th></th><th></th></tr></thead>\n",
       "<thead><tr><th>int64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th></tr></thead>\n",
       "<tr><td>1</td><td>25.0</td><td>25.0</td><td>6161.286706670988</td><td>78.49386413389895</td></tr>\n",
       "</table></div>"
      ],
      "text/plain": [
       "<QTable length=1>\n",
       "  id  xcenter ycenter    aperture_sum    aperture_sum_err\n",
       "        pix     pix                                      \n",
       "int64 float64 float64      float64           float64     \n",
       "----- ------- ------- ----------------- -----------------\n",
       "    1    25.0    25.0 6161.286706670988 78.49386413389895"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from photutils import CircularAperture, aperture_photometry\n",
    "aper = CircularAperture((25, 25), r=5)\n",
    "aperture_photometry(data, aper, error=error, method='exact')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2f5371e3",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-10-14T20:20:37.798074Z",
     "start_time": "2021-10-14T20:20:37.792809Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([6161.28670667]), array([78.49386413]))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# or without a table\n",
    "aper.do_photometry(data, error=error, method='exact')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "047c5b5f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-10-14T20:19:39.563346Z",
     "start_time": "2021-10-14T20:19:39.444650Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<matplotlib.patches.Circle at 0x13d816fa0>,)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD6CAYAAABuxZF5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjE0lEQVR4nO2deZxU1bXvv6uqq3qgmboZREABQcUoovRF43AleB3jNUYxGs29xHAfNw435iafF4fcT6b3jOa+T+anMb6owQyGRDEaNQ4BcYgJyiQigyjKJDN0Nz13Va33RxXddfY5dFU3PRVnfT+f+nTtM65zTv9qn7X22nuLqmIYRviI9LUBhmH0DSZ+wwgpJn7DCCkmfsMIKSZ+wwgpJn7DCClF+WwkIh8CB4AkkFDVKhGpAOYD44APgc+o6v6eMdMwjO5G8mnnz4i/SlX3ZC37b2Cfqt4jIrcDQ1X1to6OE5diLWHAYZpsGMahaKKeFm2WfLbNq+Y/BJ8CZmS+zwMWAx2Kv4QBnCHnH8YpDcPoiCW6MO9t8/X5FXhBRJaJyNzMspGquj3zfQcwMmhHEZkrIktFZGkrzXkbZhhGz5JvzX+Oqm4TkRHAiyKyLnulqqqIBPoPqvoA8ADAIKmwXGLD6CfkVfOr6rbM313AE8B0YKeIjALI/N3VU0YahtH95BS/iAwQkYEHvwMXAquBp4DZmc1mA0/2lJGGYXQ/+bz2jwSeEJGD2/9WVZ8TkTeB34vIHGAT8JmeM9MwjO4mp/hVdSNwasDyvYCF7g2jQLEMP8MIKSZ+wwgpJn7DCCkmfsMIKSZ+wwgpJn7DCCkmfsMIKSZ+wwgpJn7DCCkmfsMIKSZ+wwgpJn7DCCkmfsMIKSZ+wwgpJn7DCCkmfsMIKSZ+wwgpJn7DCCkmfsMIKSZ+wwgpJn7DCCkmfsMIKSZ+wwgpJn7DCCkmfsMIKSZ+wwgpJn7DCCkmfsMIKSZ+wwgpJn7DCCk5p+g2ChCRvrYgGNW+tsDIwmp+wwgpJn7DCCkmfsMIKXmLX0SiIrJCRJ7OlMeLyBIReU9E5otIvOfMNAyju+lMzX8rsDar/D3gh6o6EdgPzOlOw4wMIl34RPrppwvXYvQYeYlfRMYAnwR+kSkLMBN4LLPJPOCKHrDPMIweIt+a/0fA14BUplwJVKtqIlPeCowO2lFE5orIUhFZ2krz4dhqGEY3klP8InIZsEtVl3XlBKr6gKpWqWpVjOKuHMIwjB4gnySfs4HLReRSoAQYBPwYGCIiRZnafwywrefMPELpik8rEaeYxzEk3xe8LNzjprqQoKMpZ0HUWR1wTHeffO6RJQ91iZz/Fap6h6qOUdVxwLXAIlW9HngJmJXZbDbwZI9ZaRhGt3M47fy3AV8RkfdIxwAe7B6TDMPoDTqV26+qi4HFme8bgendb5JhGL2BdezpTfLxX3P59K7/HuDzi3ueSMTdILcdXcH1vVMpZ7V3vUiQz++1zRcX8MUR8F+PxQDywtJ7DSOkmPgNI6SY+A0jpJj4DSOkWMCvJ8kVWAtIvvEF+KLexBhfMM9ZDyDuMrcclBjU2SBgUFDNDc4lk50rA+oE/ATvNpoKqK9yJQZZADAQq/kNI6SY+A0jpJj4DSOkmM/fnbi+ZmcTdgCJOstcnz9W1OF6AIq824hT9iX9QHAcoCOCOuW4ST2JhKcsTjnI56fV2ca9p8mAJJ9ciUE4+1gMALCa3zBCi4nfMEKKid8wQor5/F0lj3bxXD6+z78HiMW828S9ZZ8/Hw8YNNmJC2iRExdwy4C6cQDXNLcpPRXgeye8Prw4Zdef15YW3yHcHAVtafWux1sGUOc0Ekl1uD7w2YUwDmA1v2GEFBO/YYQUE79hhBTz+fOlCwNx5PLxA/11x8eX4uIO12uJ/xha7N0mVeJ9zKki/2++usvcy3VcYkn4ff6IsyzS5LTzNzv+e1NAjoLj4/tzJwKegxM78Pv4NiBIEFbzG0ZIMfEbRkgx8RtGSDHxG0ZIsYBfV8ljIA5fEo+TwOMG78Af4NNSp1zmLScH+KdAS5YWOWWvHcliv+2pIqdzjBv/c2JkkYQ/IBZt9m4UbYx1XK4PGIikwQmKOusDw3DuqMBuIpAT4MtrQJAQYDW/YYQUE79hhBQTv2GEFPP5D0WOpJ7A2XFzDJzpdtLxJfAQ4OOXl3rKiXJvUk/rIH/coLXce96WAV5bkyW+XUjG3Vl+nA0clzjaEuDzN3l3itd7t4nVedfHov57WBQ00EgWgU/F8fnV7XTkiwkEDByaq/PPEZj0YzW/YYQUE79hhBQTv2GEFPP58yVHpx3wT6jhG2zTHUgzoJ3fbcd3ffyWId5y8xC/Hc2DvctaBnnXJ8r8/msq7vjNudr5/eNwUNTgvf5Erbdc7PRB0jwGDS1y/HcJGsDTHSg01vHkIBrw7PLq/HOEYeI3PCTrSkkeGEAqWYQmohBNEokliJY1UjSovq/NM7oRE3+ISTXFadx8NM3bRtK0YyTN20aQrB50yO2jA+soHr2L+JidlB61k5JjPiJa1tSLFhvdSU7xi0gJ8ApQnNn+MVX9poiMB34HVALLgH9R1YCXQaO/0bxjGDVvnMqBlZPRloAxBQ5B8kA5DevKaVg3gWpAihKUn7KewWe8RfHoHT1mr9Ez5FPzNwMzVbVORGLAayLyZ+ArwA9V9Xcicj8wB/hZD9rap/ja9bvUzu8VWtBAHG6uvtuO7/r4jZV+/7VlqFMekkJTQv3K4znw6lSaN47x256hlAbG8wFlNFBMMy3EaaSUDxlHHQO99ieKOLDiYxxY8THix25n4FmrKK9agxSl/eWU04cgaBARF0k5uRBJp42+1d9G7w4U6k4W4nsuAW324k78ETCfyJFGTvGrqgJ1mWIs81FgJnBdZvk84FscweIvZFp3DWXPby+i+YPRvnWTWcNMFjGNZVSxlMmspSggCSaFsJ4TWMY0ljGNxcxgJae1rW/ZNIq9m0ZR+8ppDLvuOYrH7O7RazIOn7x8fhGJkn61nwjcC7wPVKvqwZ/YrYD/P8voUzQlVP99KnsWnY22tteoMVr4NE9wE/dxHq/kdawIymTWMZl1fI7fAPAG/8B93MR8rqGJdCZi67YRbP/+9Qy+cAkj/uHNtrcAo/+RVzu/qiZVdSowBpgOnJjvCURkrogsFZGlrTR3zUqj07RWD2TLw1ez+/kZbcKP0cLt3M1mjmE+1+Yt/EMxnTf5JTewjdHcxZ2U0pBekYpS89xZbP5/n6V5d8XhXorRQ3QqyUdVq4GXgI8DQ0Tk4JvDGGDbIfZ5QFWrVLUqhj+X3eh+WvYMZctDn6Fpc/vL2Kms5A2mczd3chQ7O9x/I+NZyjT+ylm8SRUbmEgqOKsegAr2cyd38xancjavtS1v3jEibce2kYd/UUa3k0+0fzjQqqrVIlIKXAB8j/SPwCzSEf/ZwJM9aWiPEtSJJygRxLOLfx93thlfUo87k06xP8nHHYjD7aTjS+BxgnvNOyvZ8sgsUnVl6VPSwte5izv5LjGcQFiGxZzHM3ySZUxjOadTwxDfNuUcYCorqWIpF/ICF/E8EWdojUm8xyv8Iz/hS9zJd2mkjFRDKVseuYoRNz/O4KHeH51Iwn+Po057UbTJez/cAUHAPyqwb1ajaO7ZgX0hQN/zD5oduLA7+4jmuAARmUI6oBcl/abwe1X9johMIC38CmAF8DlV7fC9fpBU6BlyfrcY3q3kIX7fKD2usAkYirvU231OBpR5yqlB3jJAS4W3F1/TMO8/e8Nwr11Nw7P23TeYLQ9eQ7JuAJAW7B+5gvNZ5DtPLQOZx2x+xo2s5STf+lyMZyP/zs+Zw4MMY69v/ZtUcQl/Zi/DAJDSJsbOfoySUe2BwJKAmGDZbq/ISvZ4hR3f1+jbJ1Lb4ClrvbdMozcXIWiaMLeFwDfNd1DGXz8U/xJdSK3uy2u+9Zzi705M/D0n/mRjMZt/fj2t+wcDMJhqnuNizmSJ7xy/5npu5cfso9K3DiA2uI6ykfuIxBOkWqM07R5C877BgdsOoI7vcRs3c59v3RomM5NF7OQoAKIDGjjmi78mlskUNPF3P50Rv2X4HSHsfu68NuGX0sCf+Gef8Hcwkn/n5zzFpzzLI2VNjLngTSqr1jP4hC2UjNhPsfOqXLurgpr1x7BvxSS2PT+d1tr020U95dzCvTzGLB7iC4znw7Z9TmItL3AhM1jMfipI1pex66l/4ujrn8xrDhSjZzHx54vrAwYNOuEm9ThvB+5sue5MOuAfbNMdiMPtpNMyJEXDO+OpXfmxg2dhHrM5NyvwBunX8Et5lj20+wkDR+1m6uynOf6Sv3Jc5Uee7cucnjsNQ+NwwnK4HFpvi7P46UtYOe8y9mdyBxbzCaawyudmTOFtFnAlM1mEEqF+wwT2bziR8ulriDT776HbOShW7yQKBdwzachx391YTNCz8/n4R36Wj3XpLXCSDcXsnX9BW/lafsfVPObZ5nU+zvks9Aj/5Ktf4Jr5t3PyrIXEB3QuPz9W0sKJl73K1b/5Oqff8CQSTQuljoF8kmd4lks828/gZf6Dn7aV9y74BInq8k6d0+h+TPwFzv6n/pFkTTrtdgQ7+b/c4lm/ilO4lGc5QPqVITaojsvvv4tzb3uEWNnh5V1E4wnOuPkPXPnwtygfuQeAZkqYxWO8wrmebe/mDiayAQBtLGHvH/ph7CdkmPgLmETdAOqWtL/u/4wbqWRf2/omirmG+W1Nd/GhtZz+o3sZXbW2W+0YcdIHXPHg/6JkVPoHoJEyPsujVNMeJCyjkYf4ApJpMmtcPdESgPoY8/m7SlDEytf5x/ltdX3+gI4u7oQa7mCb2QNx7FvyMUilj3kOr3ElT3i2/QbfYR2TAYiVNXLFfXdTOWkLVQM/8J13QnyXpzxIvG8FtepN0BoZG+E9wMAP0PvvYsEN36Jx32A+YjRf5kf8khvaNjmX17iSBTzOrLT9b09h+GUvdXi97v0IumfRohw+vftcLNoIWM1fsGhSqHlzSlv5JqepbQnT+QFfaSuf9eXfUjlpS4/aNGj0bs6786G28jxm8wyXera5MavvV+2Kk0g1+5N2jN7BxF+g1K8/jmRt2tcfyQ6u4nHP+jv5LsnMi92Y6W8z+dMv+Y7RE4yfsYyJF72eKQl3cLdn/fks4kTSboc2F3Pgrby7iRjdjIm/QKld/rG273N4kHjW/HTrOIFFzEwXIinOu/OhXn3TPed/PkJRaboF4W2m8CrneNZ/kfvbvtcuP7n3DDM8mPgLEFVo3jKqrfw5fu1Z/zNu5OD0FkPPXMOgXu5bXzqkjuMv+Wtb+T5u8qxP25uOXTRvH5EeK9DodSzg15241asTaFInEKUBwSt3tlx3Jp1UXElUl5OsT6cGD6KGyaxr354Ij/CvbeWJVy1mbMk+zzHc4B7ApKI6T3lgxPuvcSDlXR/EzpL26P751y1gzYJ0c97jXEUNgxhMLQCV7GMCG9nIcZCM0rSvkuIxu4Kv1509OOCeuffVN+qSBfgCsZq/AGne2t5F9jRWeNatZTLVpLv6FVfWUHn6hl617SCjTvyAgRO3AtBKnGVM86yfxrK279nXY/QeJv4CpDmrf3y2iACPyAZP3oRE+q7zyZATN7V9d8VfxdK2783W379PMPEXIC0723vjnc5yzzqP+I/v2aa9XAw6of38rviz7W7ZEdy70OhZzOfPlzxml8l9DKcccEh3phx3H41AqqW9bXwYezzrNzCp7fvQiVuIRxK+TjpuAg/4ffzyiDuVrzf/P+gY7nkqJm1u+/4ux3vWVWaNAaCtsfbrDrheD0GPoTuqMPf5Hvn9eqzmL0S0tV2opXj7tzfQPkZArNzf9703Kcqa0KMR7zgF2XZnX4/Re5j4C5GsMebFGYBKOxhrr7fJjjeknH+1SNawWKr9x+YwYeIvQCTWPtBGE97X85Ks1/NkU/6z8fQEicb2/gDuG0r2m0D29Ri9h71v5Ys7rFOXjuGUAw7pzobr7iMpiBS1i6XaGWxzXNZIOvs3jmboGetpSHl/BNxOOhDUjt/krHdG9lF/f3z3PPuzZgbKtsu1OxJLtF93wPV6CHoM3TE1QHc83wLDav4CJDZ8f9v3tzjVsy676a/23bG9ZlMQtevbz+82SWbbHRu2H6P3MfEXIMVj2ofAXkqVZ122yGrWH9NrNgWRff6O8hGyr8foPUz8BUjx6HaxuO3nU1hFSca/btg2nJp1fVP77982guo1xwIgpDxJPeD90TLx9w3m83cn7lDOjh8pKa9zKgm/sxpJePeJtnjLkRaID6wmUtxMqrmYPQxnE8dwLOk29TitXM0f+FUmv//dx2ewZdpbnmP4BuIIwD+Yh9fH39jiP8aWpvaRef7+6JVtjfQX8gIjaO9cVE8Z6w7O+CYpSobt4mCKgO96nfsRdM/c++rz3/vhENv9Aav5CxCJeGv/P3C1Z332wB57F0+lsZcHy0w0xVj7xxmB9kC6o0+KdE+++Ij0/ABG72PiL1AGnto+Dt/9fNGz7kyWcHrGx9aWGK//8Ppete2Nn8+iqTo9YOixfMhlPO1Zn93Fd+CU7h1P0MgfE3+BUn7KeiIl6ea495nIc1zkWf9Nvs3BdrF3nzmXD1+d2it27Vg1kVW/aR+6++vc5ZnTbzmnsYQz0oVogoHTVveKXYYfE3+BEoknGHj6O21ld8CMy/kT1/HbtvLLd82hfveQHrWp+UAZL317LppK/1tdwAv8D37h2SZtZzqjr/zkDRT1cQpymLGAX1cJCiK5gSY3EJXw9haJBASvos3eZdEm7+9z9ow2laesoub1dLT/GT7J63ycs/hb2/qf8h8sYiY7GEXDnqE8fuPXOe0H90FAC2D2QBwQMGOPk8CTHdwDaKkv4fFb7qBm09EADKSWX/Bvnm3WciK/5nNt5YpT3/LN0BNtcgKezv0IumfuffXddwsABmI1fwETr6ym9OT3AUgR5QYepon27L0K9jOP2cRIC7n+w1Es//It1GzNHe3vDPV7hvCnG++g5p3xQLpp736+yDG0d+lNEuEGHqY5k44cP/YjSkdv71Y7jM5h4i9wKq9aiBSnm+Xe5QT+i//tWX8hL/II/0qUdES9YfNIfn/td3l7/gVo6vA61KjC+mfO5ndXf49da447uJSf8CWu41HPtt/nqyzhzHQhmmTYtS/Y6Fp9jIm/wCmqOEDFFS+3lX/If/I6H/dscy3zmc81xEn/SCSaSnjt/8zmyS/eybalkzv9FqwKO1cfx5//86ss+uaNtBxIz9gbIckDzOUW7vVsv5YT+QbfaSsPufhvxI/ei9G3iPai/zNIKvQM6YdztAVVQc6sre6gkO5MsAASd3rRlXp73MmAMk85NchbBmip8PZ7bxrmndSiYbjXrqbhaTFu+9WnaXh/HJAex/9VzmUS73m2XcbpfJ5fsppTPMtLjtnJ+CteZdi09QwYuwuJKvGIt+29OVlEw7Zh7Fsxic1PnePrNzCB93mQOczgZc/y7RzFObyWHqwTKD56J8f826NIVCkJGFS4bLfXXy/Z0+opx/f5A4SR2gZPWeu9ZRq9nZS0xRvPANCE93rVFycIiDX0w9jBEl1Ire7L653KAn5HACIw8vK/sOlnnyPVVMJOjuICXuRlzmvL/AOYxnKWMY1v803+m6+RIP3D0rR5JGt/kp5CK1razKCJWxlw1F4isQSp1iKa9g6m+t2xJOr8P1ZCipu5l3u4nQF4RbebYVzIC23Cl6IER13xPBLtf6IJI/baf4QQG3KA0dc/icTSNeUmxnEOr7GOEzzbxWnlLv6LVUzhJu5lYGY47YMkG4vZ//ZxbH1xOpufPYutL05nz/ITfMIvpYHP8zBLqeKnfMkn/C2M4VxebX/LiCQZdc3TFI+01/3+Qk7xi8hYEXlJRNaIyDsicmtmeYWIvCgiGzJ/h/a8uUZHlB7zEUdf+xRk+vtvZSxn81fm8xnftpNZx73cwjZGcx83chl/YhQfdXj8YezmIp7j+3yFrYzhYb7A6c7Q4QDPcRFn8nfWZ+XvD/v8s5Qf758g1Og7cvr8IjIKGKWqy0VkILAMuAL4PLBPVe8RkduBoap6W0fH6rc+fxBOHECizqwybhl/HEBKvINmSJm39tRyr38PkBji3aZlqDeO0FThPW9jpf/3u3r/WLY9ejna0r7vVTzGfdzk6WATxEeMYgWnsZ+hNFNMnBYGU8NUVnqa7oKoYRBf4Qc8xBc4mMgj0QRHzfozA096j7jTbb90r9+PLtnnbbOP7/f650XVjj8PSJ03DqANTgygydtJyfXvAUh6z6tOuT/690F0q8+vqtuB7ZnvB0RkLTAa+BQwI7PZPGAx0KH4jd6hbMIWxsx+jO2/v4xETTrH/nFm8TLn8X2+ynX8lqJDDE97NNs5ms61v6cQ/sgV3MqP2Up7IDA6oIGjrvozA47b3MHeRl/RKZ9fRMYBpwFLgJGZHwaAHUDgzAsiMldElorI0lb8wz0bPUPpmJ0ce9OvGDxtVduyPQxnNo9wDJv5Jt9iG0cf1jn2UMk93MZxvM9VLPAIv+y0dYy7eZ4Jvx+Td1OfiJQDLwN3qeoCEalW1SFZ6/eraod+v732985rf4vzFGp2jGXPoxeSzPS0O0gRrfwzf2Imi6hiKafyFqXO2H2e4xJjFVNYxjRe5jwWcGVbxt5BIuUNVF79FwZM3UC82mubvfb3PN3e1CciMeBx4DequiCzeKeIjFLV7Zm4gH/2R6NfUHriJkbfMY+aRVXU/e0UkrXp/v0JYjzBlTzBlQBESXASaziedymlkWKaaSFOI6VsZAKrOZkW/IN/AkQGNFJ+xmoGn/8mUeusUxDkFL+ICPAgsFZVf5C16ilgNnBP5u+TPWJhP8FN+hAJqAnc2sKtTZzkEmnyvz1E673LYlF3pt/cP+qRhDNGfnMESFA27e/o1Deo2TSBmiVTafzAm6iTpIi3mcLbTMl5joMUj9nO4DNWUn7yBuKtSagH6tPnj9c629Z4a/rian/NH6v1JvVE6721tjQFJOi4STs5nkPQSL2+pJ4QkE/NfzbwL8DbIrIys+xO0qL/vYjMATZBQHuS0e+QaIrykzdQfvIGWnZV0LBhHE0fjaR5+whad1d4JgQJoqiymuLROykZtZPS4zZTcnTWC1+rJesXEvlE+18jeIY0gAJx4I0g4iP2ER+xD4BUHFJNMZo/GkGyroxUoghNFCHRJBJLUFTWSHz0LqKl6Zo44q+AjQLD0nuNNiIlrZRO2Ab4J8j0TZ5hFDwm/nxxO3YEvB775pxr9UaVfS0GLV7/FkAanME7Ih23xkoq5lsWdWpld8CMpDsBL5CMO7a7p3Uu3x1lF/wDccTrveVYndf3dv17gKI6Jy7S4DQPB9wz3E45zn33tWgFdtIJ36+b5fYbRkgx8RtGSDHxG0ZIMZ//ULh+opPxF9QuLG6+vLuP668GDCKSq7GsyJ31Jxnke3sfa6ze+xufLPb/5qeK3HwCxy7HJXZn0oGAwUcb3bLXF3fb8MHv40ujk53X7N/Hd19z5Vvk06ZfIBl9h4PV/IYRUkz8hhFSTPyGEVLM5+8qAe3CB2eqaSPp+Ofk9vldT9PdQtxjtvr75UcbvW3/qRLvY04V+X/z1V3mntgNgQTOMOxMsNHk5Dk0e68/KE/fbcf3+fhB7fytzj7OPcprMM4QYjW/YYQUE79hhBQTv2GEFBO/YYQUC/jlS46kn/Q2HXf+UTc2FzBzjO88btnpxCLuDLUEBNYavB2KokX+QUTU7UCUo2OPuDPhgm+2XJ9tboeboOt3O+m4Ab5Wf8DPDfD5n0MeAb4QJPW4WM1vGCHFxG8YIcXEbxghxXz+rhLkI+bo/CMRJ/kkYN4MNxFI3Y48MaeTStAw1O4Q4u6MwgEDhLizEOckqHOMY6vPNrfsdsDBPxCHf0jt3ANx5Oy4E0L/Pgir+Q0jpJj4DSOkmPgNI6SYz9+d+HzJHD5+0MQfvjZqZxvXTw6aNizqtJW72wT590F5Cx0R5De7vnauKbCCfP4c1xvoz+dqxzcfPxCr+Q0jpJj4DSOkmPgNI6SYz9+T5OoPkMeAIO6goCrO73WQP+v60m67fmf9+3xxbXHb/X3rc/vvXRqIw3z8vLCa3zBCionfMEKKid8wQoqJ3zBCigX8epMuDAiSKzFIgmYLdhe4QcJ8cBOB8pnlxmdIx8G5LiXsBB7IAnxdwWp+wwgpJn7DCCk5xS8iD4nILhFZnbWsQkReFJENmb9De9ZMwzC6m3xq/l8CFzvLbgcWquokYGGmbHQW1dwf3z4pz0eTSf8npd5P0Da5Pq0J76crx8hhh3stgf5+V+6RkRc5xa+qrwD7nMWfAuZlvs8DruheswzD6Gm6Gu0fqarbM993ACMPtaGIzAXmApRQ1sXTGYbR3Rx2wE/TCduHfPdS1QdUtUpVq2IUH+7pDMPoJrpa8+8UkVGqul1ERgG7utMoI4su+bT9dBZa88/7FV2t+Z8CZme+zwae7B5zDMPoLfJp6nsU+BtwgohsFZE5wD3ABSKyAfinTNkwjAIi52u/qn72EKvO72ZbDMPoRSzDzzBCinXsORKxwJqRB1bzG0ZIMfEbRkgx8RtGSDHxG0ZIMfEbRkgx8RtGSDHxG0ZIMfEbRkgx8RtGSDHxG0ZIMfEbRkgx8RtGSDHxG0ZIMfEbRkgx8RtGSDHxG0ZIMfEbRkgx8RtGSDHxG0ZIMfEbRkgx8RtGSDHxG0ZIMfEbRkgx8RtGSDHxG0ZIMfEbRkgx8RtGSDHxG0ZIMfEbRkgx8RtGSDHxG0ZIOSzxi8jFIrJeRN4Tkdu7yyjDMHqeLotfRKLAvcAlwEnAZ0XkpO4yzDCMnuVwav7pwHuqulFVW4DfAZ/qHrMMw+hpDkf8o4EtWeWtmWWGYRQART19AhGZC8wFKKGsp09nGEaeHE7Nvw0Ym1Uek1nmQVUfUNUqVa2KUXwYpzMMozsRVe3ajiJFwLvA+aRF/yZwnaq+08E+u4FNwDBgT5dO3PsUiq2FYicUjq2FYie023qsqg7PZ4cuv/arakJEbgGeB6LAQx0JP7PPcAARWaqqVV09d29SKLYWip1QOLYWip3QNVsPy+dX1WeBZw/nGIZh9A2W4WcYIaWvxP9AH523KxSKrYViJxSOrYViJ3TB1i4H/AzDKGzstd8wQoqJ3zBCSq+Kvz/3AhSRh0Rkl4iszlpWISIvisiGzN+hfWnjQURkrIi8JCJrROQdEbk1s7xf2SsiJSLyhoi8lbHz25nl40VkSeb/YL6IxPvSzmxEJCoiK0Tk6Uy5X9oqIh+KyNsislJElmaWder595r4C6AX4C+Bi51ltwMLVXUSsDBT7g8kgK+q6knAmcDNmXvZ3+xtBmaq6qnAVOBiETkT+B7wQ1WdCOwH5vSdiT5uBdZmlfuzrZ9Q1alZ7fude/6q2isf4OPA81nlO4A7euv8edo4DlidVV4PjMp8HwWs72sbD2H3k8AF/dleoAxYDpxBOhOtKOj/oo9tHJMRzUzgaUD6sa0fAsOcZZ16/r352l+IvQBHqur2zPcdwMi+NCYIERkHnAYsoR/am3mNXgnsAl4E3geqVTWR2aQ//R/8CPgakMqUK+m/tirwgogsy3Seg04+/x7v1XekoKoqIv2qXVREyoHHgS+raq2ItK3rL/aqahKYKiJDgCeAE/vWomBE5DJgl6ouE5EZfWxOPpyjqttEZATwooisy16Zz/PvzZo/r16A/YydIjIKIPN3Vx/b04aIxEgL/zequiCzuN/aq6rVwEukX52HZDqGQf/5PzgbuFxEPiQ9MM1M4Mf0T1tR1W2Zv7tI/6hOp5PPvzfF/yYwKRM9jQPXAk/14vm7wlPA7Mz32aR96z5H0lX8g8BaVf1B1qp+Za+IDM/U+IhIKem4xFrSPwKzMpv1uZ0AqnqHqo5R1XGk/zcXqer19ENbRWSAiAw8+B24EFhNZ59/LwcpLiXdDfh94Ot9HTRxbHsU2A60kvbt5pD2+RYCG4C/ABV9bWfG1nNI+3yrgJWZz6X9zV5gCrAiY+dq4BuZ5ROAN4D3gD8AxX19Tx27ZwBP91dbMza9lfm8c1BLnX3+lt5rGCHFMvwMI6SY+A0jpJj4DSOkmPgNI6SY+A0jpJj4DSOkmPgNI6T8f+FZD0HcguCfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(data)\n",
    "region.plot(color='red', lw=4)\n",
    "aper.plot(color='blue', lw=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b667b5d6",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}