deblend_test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Deblending Test notebook.\n",
    "\n",
    "Jeff C. and Sasha B. for the SV bootcamp. First version on 05/31/18.\n",
    "\n",
    "In this notebook we use the HSC 100 sq deg. Deep dataset to test:\n",
    "\n",
    "    How to visualize the deblending of sources.\n",
    "    Comparing the \"truth\" as derived by COSMOS to the deblended sources.\n",
    "\n",
    "Borrowing some pieces from this notebook: https://github.com/RobertLuptonTheGood/notebooks/blob/master/Demos/Deblender.ipynb\n",
    "\n",
    "For info on the deblender functions available: https://github.com/lsst/analysis/blob/master/python/lsst/analysis/deblender.py"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "You'll need the lsst.analysis package to access the deblender class:\n",
    "\n",
    "1. Make a directory to contain the package (e.g., $HOME/lsst_repos/)\n",
    "2. $ cd $HOME/lsst_repos/\n",
    "3. $ git clone https://github.com/lsst/analysis\n",
    "4. $ cd analysis/\n",
    "5. $ setup -r .\n",
    "6. $ scons (build the package)\n",
    "7. $ cd ~/notebooks/\n",
    "8. If it exists, edit .user_setups. If not, create it and add the following line: “setup -r $HOME/lsst_repos/analysis/”\n",
    "9. $ source .user_setups (if no errors, then you’re good to go!)\n",
    "10. Now within your notebooks, you can import lsst.analysis.deblender (for example).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/lsst/software/stack/python/miniconda3-4.3.21/lib/python3.6/site-packages/matplotlib/__init__.py:1405: UserWarning: \n",
      "This call to matplotlib.use() has no effect because the backend has already\n",
      "been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n",
      "or matplotlib.backends is imported for the first time.\n",
      "\n",
      "  warnings.warn(_use_error_msg)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import glob\n",
    "from collections import OrderedDict\n",
    "import importlib\n",
    "\n",
    "# LSST stack imports\n",
    "import lsst.daf.persistence as dafPersist\n",
    "import lsst.analysis.deblender as deblender\n",
    "import lsst.afw.display as afwDisplay\n",
    "\n",
    "# Firefly client imports\n",
    "from firefly_client import FireflyClient\n",
    "\n",
    "# Standard libraries in support of Firefly display\n",
    "from urllib.parse import urlparse, urlunparse, ParseResult\n",
    "from IPython.display import display, Markdown, IFrame\n",
    "import os\n",
    "\n",
    "# Plotting setup:\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = (13, 8)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Grab the HSC DEEP data from the repo with the butler:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "depth = 'DEEP'\n",
    "dataPath = \"/datasets/hsc/repo/rerun/DM-13666/DEEP\"\n",
    "\n",
    "butler = dafPersist.Butler(dataPath)\n",
    "skymap = butler.get(\"deepCoadd_skyMap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function to RA and DEC center of a patch in a tract."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getPatchCenter(skymap, tract, patch):\n",
    "    '''\n",
    "    Find the center of the patch BBox in pixel coordinates and convert with tract WCS\n",
    "    An example of the format is:\n",
    "    tract = 1800\n",
    "    patch = '4,4'\n",
    "    '''\n",
    "    \n",
    "    # Convert strings to ints\n",
    "    patch = list(map(int, patch.split(','))) \n",
    "\n",
    "    # Get WCS for this tract\n",
    "    wcs = skymap.generateTract(tract).getWcs()\n",
    "\n",
    "    # Get four corners of patch\n",
    "    bbox = skymap.generateTract(tract).getPatchInfo(patch).getInnerBBox() \n",
    "\n",
    "    # Manually calculate the center of the patch and then return the values converted to degrees\n",
    "    sky = wcs.pixelToSky(np.mean((bbox.getBeginX(), bbox.getEndX())), \n",
    "                         np.mean((bbox.getBeginY(), bbox.getEndY())))\n",
    "    return sky[0].asDegrees(), sky[1].asDegrees()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure out which HSC tracts and patches overlap the COSMOS field:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 1303 patches\n"
     ]
    }
   ],
   "source": [
    "# (borrowed from https://github.com/lsst-com/notebooks/blob/master/WL_nulltest.ipynb)\n",
    "\n",
    "# In the generation 2 butler it isn't possibe to introspect the metadata for any step after the initial\n",
    "# data ingest.  This code parses the names of the files to build a dataframe for later inspection.\n",
    "\n",
    "infiles = glob.glob('/datasets/hsc/repo/rerun/DM-13666/{:s}/deepCoadd-results/HSC-I/*/*/calexp-*.fits'\n",
    "                    .format(depth))\n",
    "\n",
    "# Get the skymap for the repo so we can understand where the tracts and patches are\n",
    "skymap = butler.get('deepCoadd_skyMap')\n",
    "\n",
    "# Now get the names of the tracts and patches based on the file parsing\n",
    "tract_array = []\n",
    "patch_array = []\n",
    "for i in range(0, len(infiles)):\n",
    "    parts = infiles[i].split('/')\n",
    "    tract_array.append(int(parts[-3]))\n",
    "    patch_array.append(parts[-2])\n",
    "\n",
    "number_patches = len(patch_array)\n",
    "\n",
    "# Create arrays to hold the RA and DEC of the centers of the patches\n",
    "ra_center_array = np.empty(number_patches, dtype=np.float)\n",
    "dec_center_array = np.empty(number_patches, dtype=np.float) \n",
    "\n",
    "# Loop through the patches and figure out where their centers are in degrees\n",
    "for index, (tract, patch) in enumerate(zip(tract_array, patch_array)):\n",
    "    ra_center_array[index], dec_center_array[index] = getPatchCenter(skymap, tract, patch)\n",
    "    \n",
    "print('Found {:d} patches'.format(number_patches))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a pandas dataframe with information on each tract and patch\n",
    "data = OrderedDict({'tract': tract_array,\n",
    "        'patch': patch_array, \n",
    "        'ra_center': ra_center_array, \n",
    "        'dec_center': dec_center_array})\n",
    "\n",
    "df_tract_patch = pd.DataFrame(data=data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f607efe1278>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAFdpJREFUeJzt3X2wHXd93/H3x7JCeBIukcDXlhWbjAjImBhZNR6YZFx3mNrCwe2UoWZKTBk6ql2gZEKaEtqahz86ITNJGpfUrhI8sUsK9TQ8qMIMoZQE0tQQ2ZElS4ZBUBxLvmADtYRqD7Gkb/84K3F8uQ97dPfce672/ZrZufvw2/19vd7R9+w5+91fqgpJUn+dtdwBSJKWl4lAknrORCBJPWcikKSeMxFIUs+ZCCSp51ongiSrkvxVkp2zbEuSW5IcSLInyeZuw5QkjcsodwTvBB6cY9s1wMZm2gbcusi4JElLpFUiSLIeeC3wB3M0uQ64swbuAc5JMtVRjJKkMTq7Zbt/D/wa8Nw5tp8PPDy0fLBZNz3cKMk2BncMPPvZz77sJS95yUjBSlLf3Xvvvd+tqnVdHnPBRJDkWuDRqro3yZWL6ayqtgPbAbZs2VK7du1azOEkqXeSPNT1Mdt8NfRq4HVJvgV8DLgqyUdmtDkEXDC0vL5ZJ0macAsmgqr69apaX1UXAtcD/7Oq3jSj2Q7ghubpoSuAw1U1PfNYkqTJ0/Y3gh+T5EaAqroNuBvYChwAngDe0kl0kqSxGykRVNWfAn/azN82tL6At3UZmCRpaVhZLEk9ZyKQpJ4zEUhSz5kIJKnnTASS1HMmAknqOROBJPWciUCSes5EIEk9ZyLQilFV7HvkMFU15/yJEyfYcf8hTpw40ar98Hzb/kbtQ5p0JgKtGPunj3DTR+5j//SROed37p3mV/7r/ezcO92q/fB82/5G7UOadFmuTyyOR6BRVRX7p4+waWoNwKzzVcXOvdNce8kUSRZsPzyfpFV/o/Yx87jSYiS5t6q2dHpME4EkrRzjSAR+NSRJPWcikKSeMxFIUs+ZCCSp50wEktRzCyaCJD+Z5CtJ7k+yL8n7Z2lzZZLDSXY3083jCVeS1LU2dwQ/BK6qqp8DLgWuTnLFLO2+VFWXNtMHOo1SvTVqdfBiK4vH3YcVx5pECyaCGjjaLK5uJq9iLYlRq4MXW1k87j6sONZEOvkpZb4JWAXsBo4CH5xl+5XA94E9wGeAixc65mWXXVbSQk6cOFEPHHq8Tpw40Wr++PHj9andB+v48eMj79u2v8X0MbM/aVTArmrx7/Yo00iVxUnOAT4BvKOqHhhavwY4UVVHk2wFfreqNs6y/zZgG8CGDRsue+ihh04jdUlSfy17ZXFVPQ58Abh6xvoj1Xx9VFV3A6uTrJ1l/+1VtaWqtqxbt24RYUuSutLmqaF1zZ0ASZ4JvAb46ow256Z5s1aSy5vjfq/7cCVJXTu7RZsp4I4kqxj8A39XVe1MciNAVd0GvB64Kckx4Eng+hrlOydJ0rJZMBFU1R7gFbOsv21o/kPAh7oNTZK0FKwslqSeMxFIUs+ZCDTRyspiaexMBJpoVhZLS6DrCrW2k5XFasPKYunpWO7K4i45ZrEkjW7ZK4slSWceE4Ek9ZyJQJJ6zkQgST1nIpCknjMRSFLPmQg00crKYmnsTASaaFYWS0ug6wq1tpOVxWrDymLp6bCyWJL6zcpiSVLnTASS1HMmAknquQUTQZKfTPKVJPcn2Zfk/bO0SZJbkhxIsifJ5vGEK0nq2oKD1wM/BK6qqqNJVgN/nuQzVXXPUJtrgI3N9Erg1uavJGnCLXhH0DyxdLRZXN1MMx81ug64s2l7D3BOkqluQ5UkjUOr3wiSrEqyG3gU+FxVfXlGk/OBh4eWDzbrZh5nW5JdSXY99thjpxuzemS4irfN/LFjx3jvjr0cO3Zs5H3b9reYPmb2J02CNl8NUVXHgUuTnAN8IsnLquqBUTurqu3AdhjUEYy6v/rnZBXvSQvN3/vQ97njL/4agMt++vkj7fu6nzu/VX+L6WNmf9IkaJUITqqqx5N8AbgaGE4Eh4ALhpbXN+ukRbn2kqmn/V1ofuvFLwTg3259KWedddZI+7btb7F9zJyXltuClcVJ1gFPNUngmcCfAB+sqp1DbV4LvB3YyuBH4luq6vL5jmtlsSSNbhyVxW3uCKaAO5KsYvCbwl1VtTPJjQBVdRtwN4MkcAB4AnhLl0FKksZnwURQVXuAV8yy/rah+QLe1m1okqSlYGWxJPWciUCSes5EIEk9ZyKQpJ4zEWiilWMWS2NnItBEc8xiaQl0PfZl28kxi9WGYxZLT4djFktSvzlmsSSpcyYCSeo5E4Ek9ZyJQJJ6zkQgST1nIpCknjMRaKKVlcXS2JkINNGsLJaWQNcVam0nK4vVhpXF0tNhZbEk9duyVBYnuSDJF5LsT7IvyTtnaXNlksNJdjfTzV0GKUkanzaD1x8D3lVV9yV5LnBvks9V1f4Z7b5UVdd2H6IkaZwWvCOoqumquq+Z/wHwIHD+uAOTJC2NkZ4aSnIh8Argy7NsflWSPUk+k+TiOfbflmRXkl2PPfbYyMFKkrrXOhEkeQ7wx8AvV9XMZ9/uAzZU1cuB/wB8crZjVNX2qtpSVVvWrVt3ujFLkjrUKhEkWc0gCfxRVX185vaqOlJVR5v5u4HVSdZ2GqkkaSzaPDUU4MPAg1X123O0ObdpR5LLm+N+r8tApToDKoulSdTmjuDVwC8BVw09Hro1yY1JbmzavB54IMn9wC3A9eVVr46dCZXF0kTqukKt7WRlsUZ1JlQWS4uFlcWS1G+OWSxJ6pyJQJJ6zkQgST1nIpCknjMRqL/23AW/8zJ43zmDv3vuWu6IpGXR5u2j0plnz13w3/8FPPXkYPnww4NlgJe/YfnikpaBdwTqnarib/7kfT9KAic99ST1+Q90UgVcVhZrBTERaEWa6x/aNq9/2D99hLOPPjL7gQ8fPFUFvNg+rCzWSmEi0Iq0mNc/bJpaw7HnnDf7gZ+3nlvftJlNU2sW3cfJ4wzPS5PIymKtSCc/dZ/8x/XkfFWxc+80114yRZJZ2yT58d8IAFY/E37xllO/ESy6D2kMxlFZbCJQf+25Cz7/ATh8EJ63Hv7uzf5QrIk3jkTgU0Pqr5e/wX/4JfyNQJJ6z0QgST1nIpCknjMRSFLPmQgkqefaDF5/QZIvJNmfZF+Sd87SJkluSXIgyZ4km8cTrvpsrireUat+53v9w7j78HUTmkRt7giOAe+qqk3AFcDbkmya0eYaYGMzbQNu7TRKiTNj8HpfN6GJNOogx8CngNfMWPefgDcOLX8NmJrvOA5er1GdCYPXO5C9FovlHrw+yYXAF4GXVdWRofU7gd+oqj9vlj8P/Kuq2jVj/20M7hjYsGHDZQ899NDpZS9J6qllHbw+yXOAPwZ+eTgJjKKqtlfVlqrasm7dutM5hCSpY60SQZLVDJLAH1XVx2dpcgi4YGh5fbNOkjTh2jw1FODDwINV9dtzNNsB3NA8PXQFcLiqpjuMU5I0Jm1eOvdq4JeAvUl2N+veA2wAqKrbgLuBrcAB4AngLd2HKkkahwUTQfMD8LwvV29+yX5bV0FJkpaOlcWS1HMmAk20GrFa18piaXQmAk20Uat1rSyWTkPXFWptJyuL1cao1bpWFutMx3JXFnfJMYslaXTLWlksSTozmQgkqedMBJLUcyYCSeo5E4Ek9ZyJQJJ6zkSgiVZneGXxzGVpOZgINNHO9MpiK401EbquUGs7WVmsNs70ymIrjTUqrCyWpH6zsliS1DkTgST1nIlAknquzeD1tyd5NMkDc2y/MsnhJLub6ebuw5QkjUubwev/EPgQcOc8bb5UVdd2EpEkaUkteEdQVV8Evr8EsUiSlkGbO4I2XpVkD3AI+NWq2jdboyTbgG0AGzZs6Khrncmqiv2PHKEoNk2t4cHpH8w7/9Jzn8un936bi9Y+i4vPe9687U/UCUIGHQUuPu95AAv217aP/Y8c4ZvfPcprL5niwW//gBBeOvVc9k8fIYRN560Z9Dd9hE1Ta0iyXKdZPddFIrgP2FBVR5NsBT4JbJytYVVtB7bDoI6gg751hts/fYS33jGoN3nPa1/Cv/v0V+edv+FVG/itz36dc561mve+btO87Z86fuJUP6tXncWH/8ng0eyF+mvbx/t37OfxJ57i4ONPcudf/PWPtTnZ300fuY9b37T5VCKSltqiE0FVHRmavzvJf0yytqq+u9hjS5um1vDhN2859Sn7Z9Y+Z975l577XNaf86xTn9bnaz/zjmDT1OAT+kL9te3jRT/17FN3BL+wcd2pO4IXrX324I6g6e/WN20+NS8th1aVxUkuBHZW1ctm2XYu8J2qqiSXA/8N+Ola4MBWFkvS6MZRWbzgHUGSjwJXAmuTHATeC6wGqKrbgNcDNyU5BjwJXL9QEpAkTY4FE0FVvXGB7R9i8HipJGkFsrJYknrORCBJPWcikKSeMxFIUs+ZCDTxahHjAs83/8Chx9l36DAnTpwYedziNvNtxjaWJoGJQBNvMeMCzzf/T++4l7fesYude6dHHre4q7GNpYnQ9diXbSfHLFZbixkXeL75vQf/bz1w8PE6fvz4yOMWdzW2sTQqHLNYkvrNMYslSZ0zEUhSz5kIJKnnTASS1HMmAknqOROBJPWciUATr6wslsbKRKCJZ2WxNGZdV6i1nawsVltWFks/gpXFktRvy1JZnOT2JI8meWCO7UlyS5IDSfYk2dxlgJKk8WrzG8EfAlfPs/0aYGMzbQNuXXxYkqSlsmAiqKovAt+fp8l1wJ3N11f3AOckmeoqQEnSeHXx1ND5wMNDywebdT8mybYku5LseuyxxzroWpK0WEv6+GhVba+qLVW1Zd26dUvZtSRpDl0kgkPABUPL65t1kqQVoItEsAO4oXl66ArgcFVNd3BcCbCyWBq3No+PfhT438DPJjmY5K1JbkxyY9PkbuCbwAHg94F/PrZo1UtWFktj1nWFWtvJymK1ZWWx9CNYWSxJ/eaYxZKkzpkIJKnnTASS1HMmAknqOROBJPWciUCSes5EoBWlxlBZPPMR6q76mKuyWJo0JgKtKOOoLJ5Z4TvuymJp0lhQphWlqtg/fYRNU2sAFjW/75HDhLDpvDUk6byPqmLn3mmuvWSKJKfWD/cljWocBWUmAklaQawsliR1zkQgST1nIpCknjMRSFLPmQgkqedMBJLUcyYCSeq5VokgydVJvpbkQJJ3z7L9yiSHk+xuppu7D1VaXqO+YkJaKdoMXr8K+D3gGmAT8MYkm2Zp+qWqurSZPtBxnNKyG/UVE9JKcXaLNpcDB6rqmwBJPgZcB+wfZ2DSpNk0tYZb37T51Ksk2sxLK0Gbr4bOBx4eWj7YrJvpVUn2JPlMkos7iU6aIEm4+LznkaTVvLRStLkjaOM+YENVHU2yFfgksHFmoyTbgG0AGzZs6KhrSdJitLkjOARcMLS8vll3SlUdqaqjzfzdwOoka2ceqKq2V9WWqtqybt26RYQtSepKm0Twl8DGJBcl+QngemDHcIMk56a5F05yeXPc73UdrCSpewt+NVRVx5K8HfgssAq4var2Jbmx2X4b8HrgpiTHgCeB68vn5yRpRXA8AklaQRyPQJLUOROBJPWciUCSes5EIEk9ZyKQpJ4zEUhSz5kIJKnnTASS1HMmAknqOROBJPWciUCSes5EIEk9ZyKQpJ4zEUhSz5kIJKnnTASS1HMmAknqOROBJPWciUCSeq5VIkhydZKvJTmQ5N2zbE+SW5rte5Js7j5USdI4LJgIkqwCfg+4BtgEvDHJphnNrgE2NtM24NaO45QkjUmbO4LLgQNV9c2q+hvgY8B1M9pcB9xZA/cA5ySZ6jhWSdIYnN2izfnAw0PLB4FXtmhzPjA93CjJNgZ3DAA/TPLASNEuj7XAd5c7iBaMs1srIc6VECMYZ9d+tusDtkkEnamq7cB2gCS7qmrLUvZ/OoyzW8bZnZUQIxhn15Ls6vqYbb4aOgRcMLS8vlk3ahtJ0gRqkwj+EtiY5KIkPwFcD+yY0WYHcEPz9NAVwOGqmp55IEnS5Fnwq6GqOpbk7cBngVXA7VW1L8mNzfbbgLuBrcAB4AngLS363n7aUS8t4+yWcXZnJcQIxtm1zuNMVXV9TEnSCmJlsST1nIlAknquk0SQ5PYkj86sC0jyjiRfTbIvyW8Orf/15nUUX0vy9+Y45vOTfC7J15u/f2sccSZ5X5JDSXY309Zm/T8eWrc7yYkkl85yzFn3X8I4L0zy5ND62+Y4Zqfnc8QYX5Pk3iR7m79XzXHMZT2XzbaJuTaHtr0rSSVZ2yxP1LU5T5zLcm2eRpwTdX3OFWezbjzXZ1UtegJ+AdgMPDC07u8A/wN4RrP8gubvJuB+4BnARcA3gFWzHPM3gXc38+8GPjimON8H/OoC+10CfGOObQvuP844gQuH281zzE7P54gxvgI4r5l/GXBoQs/lRF2bzfoLGDyo8RCwdhKvzbniXK5r8zTinKjrc544x3Z9dnJHUFVfBL4/Y/VNwG9U1Q+bNo82668DPlZVP6yq/8PgSaPLZznsdcAdzfwdwN8fU5xtvJHBqzWWxCLinE+n53OUGKvqr6rqkWZxH/DMJM9YTP9tjXguJ/Ha/B3g14C5nuqYlGtzoTjns6znc0Kvz9nO59iuz3H+RvBi4OeTfDnJnyX52836uV5HMdML60e1CN8GXji+UHlHBm9NvX2O26h/BHx0Eft3Za5+LmpuVf8syc/Pse9Snc+FzsU/BO47+QHhNPYfZ5wTdW0muY7Bp9P752m27NfmAnFOzLXZ8nwu+/U5T5xjuz7HmQjOBp4PXAH8S+CuJDmdA9XgHmdcz7neCrwIuJTBu5F+a3hjklcCT1TVXO9Fmnf/JYhzGthQVZcCvwL8lyRr5jvQGM/nQufyYuCDwD87nf2XKs5RjOtcJnkW8B7g5nnaLPu1uUCcE3Nttjyfy359tolzFG3P5zgTwUHg4zXwFeAEg5c6tX0dxXfSvMG0+fvoLG0Wraq+U1XHq+oE8Pv8+K3W9czziavF/mONs7lN/F4zfy+D7w1fPMshxn4+5zsXSdYDnwBuqKpvjLr/EsU5SdfmzzD4Hvj+JN9qYrkvyblDbSbh2pwzzkm6NueLs+l3Uq7P+eIc2/U5zkTwSQY/GJPkxcBPMHiz3w7g+iTPSHIRgzEMvjLL/juANzfzbwY+NY4g8/TXZf8DYPjpkrOANzDPd7Dz7d+lufpJsi6DMSNI8iIG5/Obsxxi7OdznhjPAT7N4Aes/zXq/ksVJxN0bVbV3qp6QVVdWFUXMvhgtbmqvt38N0zEtTlfnJN0bS4Q58Rcnwv8fx/f9bnQr8ltJgafSqaBp5rA38rgH/6PMDhZ9wFXDbX/1ww+HXwNuGZo/R8AW5r5nwI+D3ydwdNHzx9TnP8Z2AvsaU7g1FD7K4F7ZjnOcJxz7r8UcTL4TnMfsLs5z7+4FOdzxBj/DfD/mhhPTi+YJcZlPZeTdm3O2P4thp4amqRrc644l+vaPI04J+r6XOD/+1iuT18xIUk9Z2WxJPWciUCSes5EIEk9ZyKQpJ4zEUhSz5kIJKnnTASS1HP/H+MYViA96nmJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5ec5078b00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make sure the readout of the tract patch metadata file looks reasonable\n",
    "plt.figure()\n",
    "plt.scatter(df_tract_patch.ra_center, df_tract_patch.dec_center, marker='.', s=1)\n",
    "# Centered on COSMOS field\n",
    "plt.xlim(160., 140.)\n",
    "plt.ylim(0., 4.)\n",
    "\n",
    "# Pick one patch\n",
    "\n",
    "#df_tract_patch.query('ra_center > 140 & ra_center < 160 & dec_center<4')\n",
    "selected_patch = df_tract_patch.query(\"tract==9813 and patch == '4,4'\")\n",
    "#selected_patch = df_tract_patch.query(\"tract==10056 and patch == '7,0'\")\n",
    "plt.scatter(selected_patch.ra_center, selected_patch.dec_center)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick patches that are roughly within the COSMOS footprint:\n",
    "cosmos = (df_tract_patch['ra_center'] > 149.5) & (df_tract_patch['ra_center'] < 151) & (df_tract_patch['dec_center'] > 1.5) & (df_tract_patch['dec_center'] < 2.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>tract</th>\n",
       "      <th>patch</th>\n",
       "      <th>ra_center</th>\n",
       "      <th>dec_center</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>900</th>\n",
       "      <td>9812</td>\n",
       "      <td>0,4</td>\n",
       "      <td>149.507475</td>\n",
       "      <td>2.231262</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>909</th>\n",
       "      <td>9812</td>\n",
       "      <td>0,5</td>\n",
       "      <td>149.507570</td>\n",
       "      <td>2.417912</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>930</th>\n",
       "      <td>9812</td>\n",
       "      <td>0,1</td>\n",
       "      <td>149.507190</td>\n",
       "      <td>1.671328</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>938</th>\n",
       "      <td>9812</td>\n",
       "      <td>0,2</td>\n",
       "      <td>149.507285</td>\n",
       "      <td>1.857966</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>946</th>\n",
       "      <td>9812</td>\n",
       "      <td>0,3</td>\n",
       "      <td>149.507380</td>\n",
       "      <td>2.044612</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1016</th>\n",
       "      <td>9813</td>\n",
       "      <td>0,4</td>\n",
       "      <td>150.995078</td>\n",
       "      <td>2.231262</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1018</th>\n",
       "      <td>9813</td>\n",
       "      <td>4,1</td>\n",
       "      <td>150.247887</td>\n",
       "      <td>1.671469</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1019</th>\n",
       "      <td>9813</td>\n",
       "      <td>1,5</td>\n",
       "      <td>150.808365</td>\n",
       "      <td>2.418002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1021</th>\n",
       "      <td>9813</td>\n",
       "      <td>5,2</td>\n",
       "      <td>150.061127</td>\n",
       "      <td>1.858114</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1024</th>\n",
       "      <td>9813</td>\n",
       "      <td>6,3</td>\n",
       "      <td>149.874323</td>\n",
       "      <td>2.044742</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1027</th>\n",
       "      <td>9813</td>\n",
       "      <td>3,1</td>\n",
       "      <td>150.434624</td>\n",
       "      <td>1.671461</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1028</th>\n",
       "      <td>9813</td>\n",
       "      <td>0,5</td>\n",
       "      <td>150.995173</td>\n",
       "      <td>2.417912</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1030</th>\n",
       "      <td>9813</td>\n",
       "      <td>4,2</td>\n",
       "      <td>150.247887</td>\n",
       "      <td>1.858124</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1033</th>\n",
       "      <td>9813</td>\n",
       "      <td>5,3</td>\n",
       "      <td>150.061103</td>\n",
       "      <td>2.044775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1035</th>\n",
       "      <td>9813</td>\n",
       "      <td>6,4</td>\n",
       "      <td>149.874276</td>\n",
       "      <td>2.231404</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1037</th>\n",
       "      <td>9813</td>\n",
       "      <td>2,1</td>\n",
       "      <td>150.621356</td>\n",
       "      <td>1.671434</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1039</th>\n",
       "      <td>9813</td>\n",
       "      <td>3,2</td>\n",
       "      <td>150.434647</td>\n",
       "      <td>1.858114</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1042</th>\n",
       "      <td>9813</td>\n",
       "      <td>4,3</td>\n",
       "      <td>150.247887</td>\n",
       "      <td>2.044786</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1043</th>\n",
       "      <td>9813</td>\n",
       "      <td>8,4</td>\n",
       "      <td>149.500696</td>\n",
       "      <td>2.231262</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1046</th>\n",
       "      <td>9813</td>\n",
       "      <td>5,4</td>\n",
       "      <td>150.061080</td>\n",
       "      <td>2.231440</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1047</th>\n",
       "      <td>9813</td>\n",
       "      <td>1,1</td>\n",
       "      <td>150.808081</td>\n",
       "      <td>1.671390</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1048</th>\n",
       "      <td>9813</td>\n",
       "      <td>2,2</td>\n",
       "      <td>150.621404</td>\n",
       "      <td>1.858084</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1050</th>\n",
       "      <td>9813</td>\n",
       "      <td>3,3</td>\n",
       "      <td>150.434671</td>\n",
       "      <td>2.044775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1052</th>\n",
       "      <td>9813</td>\n",
       "      <td>7,4</td>\n",
       "      <td>149.687480</td>\n",
       "      <td>2.231345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1055</th>\n",
       "      <td>9813</td>\n",
       "      <td>4,4</td>\n",
       "      <td>150.247887</td>\n",
       "      <td>2.231452</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1056</th>\n",
       "      <td>9813</td>\n",
       "      <td>8,5</td>\n",
       "      <td>149.500601</td>\n",
       "      <td>2.417912</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1057</th>\n",
       "      <td>9813</td>\n",
       "      <td>8,1</td>\n",
       "      <td>149.500981</td>\n",
       "      <td>1.671327</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1058</th>\n",
       "      <td>9813</td>\n",
       "      <td>0,1</td>\n",
       "      <td>150.994794</td>\n",
       "      <td>1.671328</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1059</th>\n",
       "      <td>9813</td>\n",
       "      <td>1,2</td>\n",
       "      <td>150.808152</td>\n",
       "      <td>1.858035</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1061</th>\n",
       "      <td>9813</td>\n",
       "      <td>2,3</td>\n",
       "      <td>150.621451</td>\n",
       "      <td>2.044742</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1064</th>\n",
       "      <td>9813</td>\n",
       "      <td>3,4</td>\n",
       "      <td>150.434695</td>\n",
       "      <td>2.231440</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1065</th>\n",
       "      <td>9813</td>\n",
       "      <td>7,5</td>\n",
       "      <td>149.687409</td>\n",
       "      <td>2.418002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1066</th>\n",
       "      <td>9813</td>\n",
       "      <td>7,1</td>\n",
       "      <td>149.687693</td>\n",
       "      <td>1.671390</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1067</th>\n",
       "      <td>9813</td>\n",
       "      <td>4,5</td>\n",
       "      <td>150.247887</td>\n",
       "      <td>2.418118</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1069</th>\n",
       "      <td>9813</td>\n",
       "      <td>8,2</td>\n",
       "      <td>149.500886</td>\n",
       "      <td>1.857966</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1070</th>\n",
       "      <td>9813</td>\n",
       "      <td>0,2</td>\n",
       "      <td>150.994889</td>\n",
       "      <td>1.857966</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1072</th>\n",
       "      <td>9813</td>\n",
       "      <td>1,3</td>\n",
       "      <td>150.808223</td>\n",
       "      <td>2.044688</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1075</th>\n",
       "      <td>9813</td>\n",
       "      <td>2,4</td>\n",
       "      <td>150.621499</td>\n",
       "      <td>2.231404</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1076</th>\n",
       "      <td>9813</td>\n",
       "      <td>6,5</td>\n",
       "      <td>149.874228</td>\n",
       "      <td>2.418066</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1077</th>\n",
       "      <td>9813</td>\n",
       "      <td>6,1</td>\n",
       "      <td>149.874418</td>\n",
       "      <td>1.671434</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1078</th>\n",
       "      <td>9813</td>\n",
       "      <td>3,5</td>\n",
       "      <td>150.434719</td>\n",
       "      <td>2.418105</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1080</th>\n",
       "      <td>9813</td>\n",
       "      <td>7,2</td>\n",
       "      <td>149.687622</td>\n",
       "      <td>1.858035</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1082</th>\n",
       "      <td>9813</td>\n",
       "      <td>8,3</td>\n",
       "      <td>149.500791</td>\n",
       "      <td>2.044612</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1083</th>\n",
       "      <td>9813</td>\n",
       "      <td>0,3</td>\n",
       "      <td>150.994983</td>\n",
       "      <td>2.044612</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1086</th>\n",
       "      <td>9813</td>\n",
       "      <td>1,4</td>\n",
       "      <td>150.808294</td>\n",
       "      <td>2.231345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1087</th>\n",
       "      <td>9813</td>\n",
       "      <td>5,5</td>\n",
       "      <td>150.061056</td>\n",
       "      <td>2.418105</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1088</th>\n",
       "      <td>9813</td>\n",
       "      <td>5,1</td>\n",
       "      <td>150.061151</td>\n",
       "      <td>1.671461</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1089</th>\n",
       "      <td>9813</td>\n",
       "      <td>2,5</td>\n",
       "      <td>150.621546</td>\n",
       "      <td>2.418066</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1091</th>\n",
       "      <td>9813</td>\n",
       "      <td>6,2</td>\n",
       "      <td>149.874371</td>\n",
       "      <td>1.858084</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1093</th>\n",
       "      <td>9813</td>\n",
       "      <td>7,3</td>\n",
       "      <td>149.687551</td>\n",
       "      <td>2.044688</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1146</th>\n",
       "      <td>9814</td>\n",
       "      <td>8,4</td>\n",
       "      <td>150.988300</td>\n",
       "      <td>2.231262</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1152</th>\n",
       "      <td>9814</td>\n",
       "      <td>8,5</td>\n",
       "      <td>150.988205</td>\n",
       "      <td>2.417912</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1153</th>\n",
       "      <td>9814</td>\n",
       "      <td>8,1</td>\n",
       "      <td>150.988584</td>\n",
       "      <td>1.671327</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1158</th>\n",
       "      <td>9814</td>\n",
       "      <td>8,2</td>\n",
       "      <td>150.988489</td>\n",
       "      <td>1.857966</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1165</th>\n",
       "      <td>9814</td>\n",
       "      <td>8,3</td>\n",
       "      <td>150.988394</td>\n",
       "      <td>2.044612</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      tract patch   ra_center  dec_center\n",
       "900    9812   0,4  149.507475    2.231262\n",
       "909    9812   0,5  149.507570    2.417912\n",
       "930    9812   0,1  149.507190    1.671328\n",
       "938    9812   0,2  149.507285    1.857966\n",
       "946    9812   0,3  149.507380    2.044612\n",
       "1016   9813   0,4  150.995078    2.231262\n",
       "1018   9813   4,1  150.247887    1.671469\n",
       "1019   9813   1,5  150.808365    2.418002\n",
       "1021   9813   5,2  150.061127    1.858114\n",
       "1024   9813   6,3  149.874323    2.044742\n",
       "1027   9813   3,1  150.434624    1.671461\n",
       "1028   9813   0,5  150.995173    2.417912\n",
       "1030   9813   4,2  150.247887    1.858124\n",
       "1033   9813   5,3  150.061103    2.044775\n",
       "1035   9813   6,4  149.874276    2.231404\n",
       "1037   9813   2,1  150.621356    1.671434\n",
       "1039   9813   3,2  150.434647    1.858114\n",
       "1042   9813   4,3  150.247887    2.044786\n",
       "1043   9813   8,4  149.500696    2.231262\n",
       "1046   9813   5,4  150.061080    2.231440\n",
       "1047   9813   1,1  150.808081    1.671390\n",
       "1048   9813   2,2  150.621404    1.858084\n",
       "1050   9813   3,3  150.434671    2.044775\n",
       "1052   9813   7,4  149.687480    2.231345\n",
       "1055   9813   4,4  150.247887    2.231452\n",
       "1056   9813   8,5  149.500601    2.417912\n",
       "1057   9813   8,1  149.500981    1.671327\n",
       "1058   9813   0,1  150.994794    1.671328\n",
       "1059   9813   1,2  150.808152    1.858035\n",
       "1061   9813   2,3  150.621451    2.044742\n",
       "1064   9813   3,4  150.434695    2.231440\n",
       "1065   9813   7,5  149.687409    2.418002\n",
       "1066   9813   7,1  149.687693    1.671390\n",
       "1067   9813   4,5  150.247887    2.418118\n",
       "1069   9813   8,2  149.500886    1.857966\n",
       "1070   9813   0,2  150.994889    1.857966\n",
       "1072   9813   1,3  150.808223    2.044688\n",
       "1075   9813   2,4  150.621499    2.231404\n",
       "1076   9813   6,5  149.874228    2.418066\n",
       "1077   9813   6,1  149.874418    1.671434\n",
       "1078   9813   3,5  150.434719    2.418105\n",
       "1080   9813   7,2  149.687622    1.858035\n",
       "1082   9813   8,3  149.500791    2.044612\n",
       "1083   9813   0,3  150.994983    2.044612\n",
       "1086   9813   1,4  150.808294    2.231345\n",
       "1087   9813   5,5  150.061056    2.418105\n",
       "1088   9813   5,1  150.061151    1.671461\n",
       "1089   9813   2,5  150.621546    2.418066\n",
       "1091   9813   6,2  149.874371    1.858084\n",
       "1093   9813   7,3  149.687551    2.044688\n",
       "1146   9814   8,4  150.988300    2.231262\n",
       "1152   9814   8,5  150.988205    2.417912\n",
       "1153   9814   8,1  150.988584    1.671327\n",
       "1158   9814   8,2  150.988489    1.857966\n",
       "1165   9814   8,3  150.988394    2.044612"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_tract_patch[cosmos]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick one patch within the COSMOS field (from the table above), and ingest its data:\n",
    "\n",
    "dataId = dict(tract=9813, patch='4,4', filter='HSC-I')\n",
    "coadd = butler.get(\"deepCoadd_calexp\", dataId, immediate=True)\n",
    "\n",
    "ss = butler.get('deepCoadd_meas', dataId)\n",
    "families = deblender.Families(ss, butler, nChildMin=2)\n",
    "### NOTE: nChildMin>2 seems to fail, for reasons we don't understand...\n",
    "\n",
    "nameNchild = \"deblend_nChild\"                \n",
    "if nameNchild not in ss.getSchema().extract(\"*\").keys():\n",
    "    nameNchild = \"deblend_nchild\"\n",
    "\n",
    "assert ss.getSchema().find(nameNchild) is not None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell is simply a way to identify indices of parents that have multiple children.\n",
    "# WARNING: Don't use print(families) -- this is too big and will cause problems...\n",
    "# Format of families: families[parent][children]\n",
    "\n",
    "#for i in range(len(families)):\n",
    "#    print(i,len(families[i][1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the firefly backend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_channel = '{}_test_channel'.format(os.environ['USER'])\n",
    "server = 'https://lsst-lspdev.ncsa.illinois.edu'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use this cell if you want an inline \"iframe\":\n",
    "\n",
    "#ff='{}/firefly/slate.html?__wsch={}'.format(server, my_channel)\n",
    "#IFrame(ff,800,600)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "Open your web browser to <a href=\"https://lsst-lspdev.ncsa.illinois.edu/firefly//slate.html?__wsch=jcarlin_test_channel\"\" target=\"_blank\">this link</a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "afwDisplay.setDefaultBackend('firefly')\n",
    "afw_display = afwDisplay.getDisplay(frame=1, name=my_channel)\n",
    "disp1  = afwDisplay.Display(1)\n",
    "disp2 = afwDisplay.Display(2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The showBlend utility should interactively overplot deblended parents+children on an image, but right now interactivity\n",
    "# doesn't seem to work, because callback from Firefly isn't currently working:\n",
    "\n",
    "#deblender.showBlend(coadd, families, display=disp2, imageDisplay=disp1, mtv=True, key='d')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's look at a deblended family!\n",
    "# The first index (\"2\" in this example) is the index number of the deblended object. The second index (\"0\" or \"1\" here)\n",
    "# corresponds to parent (0) and children (1). See the deblender.py code referenced above for more info.\n",
    "deblender.plotDeblendFamily(coadd, families[2][0], families[2][1], display=disp2,symbolSize=25,plotb=True,background=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "150.2725442471663 2.1339181894148607\n"
     ]
    }
   ],
   "source": [
    "# Check the coordinates of the ith deblended parent:\n",
    "# (where here i=0)\n",
    "print(families[0][0]['coord_ra'].asDegrees(),families[0][0]['coord_dec'].asDegrees())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the HDF5 file of COSMOS data created in this notebook: \n",
    "# https://github.com/lsst-com/notebooks/blob/master/hsc_cosmos_match.ipynb\n",
    "outfile_cosmos='/project/shared/data/COSMOS_catalogs/cosmos_acs_iphot_200709.h5'\n",
    "df_cosmos = pd.read_hdf(outfile_cosmos, 'df')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['number', 'mag_iso', 'magerr_iso', 'mag_isocor', 'magerr_isocor',\n",
       "       'mag_petro', 'magerr_petro', 'petro_radius', 'mag_aper', 'magerr_aper',\n",
       "       'mag_auto', 'magerr_auto', 'mag_best', 'magerr_best', 'flux_auto',\n",
       "       'fluxerr_auto', 'kron_radius', 'background', 'threshold', 'flux_max',\n",
       "       'flux_radius', 'isoarea_image', 'x_image', 'y_image', 'xmin_image',\n",
       "       'ymin_image', 'xmax_image', 'ymax_image', 'xpeak_image', 'ypeak_image',\n",
       "       'alphapeak_j2000', 'deltapeak_j2000', 'a_image', 'b_image', 'ra', 'dec',\n",
       "       'theta_image', 'mu_threshold', 'mu_max', 'isoarea_world', 'x_world',\n",
       "       'y_world', 'a_world', 'b_world', 'theta_world', 'flags', 'fwhm_image',\n",
       "       'fwhm_world', 'cxx_image', 'cyy_image', 'cxy_image', 'elongation',\n",
       "       'class_star', 'field', 'mu_class', 'x', 'y', 'z', 'spt_ind', 'cntr'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_cosmos.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick out the RA, Dec coordinates of all the parents:\n",
    "parents = [family[0] for family in families]\n",
    "parents_ra = [parent['coord_ra'].asDegrees() for parent in parents]\n",
    "parents_dec = [parent['coord_dec'].asDegrees() for parent in parents]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.62275 rad"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "parents[0]['coord_ra']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEJCAYAAABc/7oDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJztnX2UHNV5p5+fRhqWCMXYCmAsNAE74tPiw8xCesGmZRFjWG9IIuKQw5GMDRFwHBYSsgGRwLLWycqxY8fsZjEZA1m0kWNjS8ROYjDyHAaM1cKRiGCMhPiWLa0wtoyNVrtRS6N3/6hqq2jmo3umuutW1fucM2eqb92qrnv73vu773s/SmaG4ziO46TBtKwfwHEcxykOLiqO4zhOarioOI7jOKnhouI4juOkhouK4ziOkxouKo7jOE5qTCgqkuZKeljSZklPS7pulDiXSXpK0rCkdZJOi8NPkLQp8fe6pOvjc2+TtFbSc/H/tybut0zS85K2SrogzQQ7juM4nUMTrVORdDRwtJk9IWkWsBH4DTPbnIjz74AtZvaapAuB28zs7Kb79AA7gLPNbJukTwE/MbNPSroJeKuZ3SjpZODvgLOAdwDfAo43s5HUUu04juN0hAktFTPbaWZPxMe7gS3AnKY468zstfjjeuCYUW61EHjBzLbFny8G7o2P7wV+IxH+JTPba2YvAc8TCYzjOI4TONPbiSzpWOAM4PFxol0BPDBK+KVEFkiDo8xsZ3z8CnBUfDyHSJgabKdJxOJnWQosBZg5c+aZJ5544pu+cM8e2L0bZs2CmTPHeeICsmcPbN0KZiDB3LkwMlLOvHDKS3M9OOEEL/9JNm7c+GMzOyLNe7YsKpIOA1YD15vZ62PEWUAkKuc2hfcCvw4sG+06MzNJbe0XY2YDwABAf3+/bdiwoZ3LM6VWg6EhqFahUsn/9zhOqKxYAbfcEnWopk2DJUtg2SitUFnriqRtE8dqj5ZERdIMIkFZZWZrxohzKnAXcKGZ7Wo6fSHwhJn9MBH2Q0lHm9nOeNzm1Th8BzA3Ee+YOKwQ1GqwcCHU69DbC4ODnSvElUq5KojjNFOtRvWsUd+q1TfH6WadLAOtzP4ScDfRQPxnx4jTB6wBFpvZs6NE+V3e6PoC+Drwkfj4I8DXEuGXSjpE0nHAPOC74z3jnj0TpSIchoaiwjsyEv0fGsr6iRynuFQqkUgsXz62WHidTJdWLJVzgMXAsKRNcdjNQB+Amd0J3ArMBu6INIj9ZtYPIGkm8GvAVU33/SRwn6QrgG3Ah+P7PS3pPmAzsB/4+EQzv3bvjv4PD9TYtXqI2YuqzF8aZlejlZ7TZCmrCe844zGRxd7JOllGJpxSnAdOOqnf7vuD/867rlpIL3Xq9PLCXw8GKyydaPzdhHecyVPWDpmkjQ0DIC3amv0VKjNnwq7VQ5xEnemMYNTZtXoIAhWVTox1jGbCl6lyOM5U8PHH9CjMNi2zF1Wp08s+ethHL7MXVbN+pK7SMOF7etyEdxwnOwphqQDMX1phmMHgx1Q6RWNAsowmvOM44VCIMZW8rVNxHMcJgU6MqRTG/eU4TrjUatFCxFot6ydxOk1h3F9FpayzUpzikOeZiV7/2sdFJWDyXBkdp0FeZyZ6/Zsc7v4KGF/p6xSBvM5M9Po3OdxSCRhf6VteiuR2yevMRK9/k8NnfwVOkRoXpzXc7RIORa9/vqK+hDRW+jZmzxS1cDsHyesYRBHxlfbt46KSA7znWi7c7eLkGReVUKnV2LZyiEeoUqMSXM+16G6BLMnrGISTLaHUSReVEKnVGFmwkDl761xCL3dNH2T69KiUhNBzdcup8xTS7dLU6oXSCBaBkOqki0qIDA2hep2eeMfl944MccqVFfr6wqiA7vN32qap1Rv+3CALr68E0QgWgZDqpItKiFSrWG8v+/bW2Ucv35lRZcWScCpdFj5/79XmnKZWb9fqIer1ShCNYBEIaRzORSVEKhV6Hh5kezymsmJJJagK122ff0imfVlIXcSbWr3Zi6r0fjuMRrAIhDQO56ISKpUKv1ypsCTr5xiDbvr8QzLty0BHRLyp1ZtfqTA4P4xGsCiEMg43oahImgusBI4CDBgws9ub4lwG3AgI2A1cY2ZPxucOB+4C3h1f/zEzq0n6MnBCfIvDgZ+a2emSjgW2AFvjc+vN7OrxnnHPHl/DUWRCMu3LQMdEvKnVC6URdNKlFUtlP3CDmT0haRawUdJaM9uciPMScJ6ZvSbpQmAAODs+dzvwoJldIqkX+AUAM/udxsWSPgP8LHG/F8zs9FYTsXUr3HKLu0aKSkimfRlwEXemwoSiYmY7gZ3x8W5JW4A5wOZEnHWJS9YDxwBIegvwPuDyOF4dqCfvL0nAh4H3TzYRZu4aKToh9WqLPmmgHREvel447dPWmErsmjoDeHycaFcAD8THxwE/Av5G0mnARuA6M9uTiP9e4Idm9lwi7DhJm4islz81s2+P/1wwbZr3qpzOU5ZJA62IeFnyIs9kIfotb30v6TBgNXC9mb0+RpwFRKJyYxw0HXgP8HkzOwPYA9zUdNnvAn+X+LwT6IvdX38IfFHSL47yXUslbZC04cgjf8ry5V6onc7j26EfxPMibBqif8st0f9uvXWzJVGRNINIUFaZ2Zox4pxKNCB/sZntioO3A9vNrGHZfJVIZBrXTAd+C/hyI8zM9jauN7ONwAvA8c3fZ2YDZtZvZv3HHHM4y5a5oJSZbr2uNq/vBukEnhdhk5XotzL7S8DdwBYz++wYcfqANcBiM3u2EW5mr0j6gaQTzGwrsJDEWAxwPvCMmW1P3OsI4CdmNiLpncA84MVJpM0pCd10w2Q9aaCT7ox27511Xjjjk9WEi1bGVM4BFgPD8TgHwM1AH4CZ3QncCswG7og0iP2JPfqvBVbFM79eBD6auPelvNH1BdHA/ick7QMOAFeb2U/aTZgzOdJqtLrpy+32OpasJg10Ujwne++QJlA4byQr0W9l9tdjROtPxotzJXDlGOc2AaO+BMbMLh8lbDWRqy135H0mTFqNVrcHcMsyBbaT4ukLTItJFqLvK+pTIo8zYZpFMK2GJQvLoQxumE6KZ1mEOU3y3onsFC4qKZG3nt5oIphWw5JFA1UGN0wnxbMswpwWeexEdgsXlZTIW09vNBFctiydhsUbqM7RSfEsgzCnRd46kWMza2badyyEqLzyStRz6OaP2mz65q0hHUsE02pYvIFyikzeOpGjEU2/n/em5RpTRWaW9j27jtRvhx66oWsmaFFMX/cJh4P/FuEy1m+T999sxQq4+eZ+zDaMOxGrXQphqUB3TdCimL5uTYRBUTopRWS83ybv9SeyruxA2vdteZuW0OmmCeoriZ008e1OwmC0XRmK/NtEgvjcsxPFa5dCWCpz5sBXvtKlXkOtxjtWDnHbBVVeenuFJQG95tfJJ0Xwz+edsSyS4v82u/dMHKc9CiEqb3979wRlZMFC5uyt8/v0clHvICxxRXGmRt4meRSRsVza/tu0TyFEpWsMDaF6nR5GMOqcs2+IoaGw3h/v5JM8+OfzPjA9HuNZJHn4bULCRaUdqlWst5d9e+vso5fvzKiyopr1QzlO5yn6ZAK3SNLDRaUdKhV6Hh5k+8ohHqHKiiVupTjloCgzHsfDLZJ0cFFpl0qFX65UWJL1czhOFyn+gLWTFi4qjuNMiLuHnFYppagUecDRcTqFu4fCJLT2rHSiUvQBR8dxykOI7VlhVtS3SpFXyDqOUy6a27OVK9+8K0C3KZ2l4gOO+SQUEz+U53AceGN7Nn063HNPJDBZWi2lE5VxBxxrNbbF04Xn+XThYAjFxA/lORynQbI9+/734QtfyH7a94TuL0lzJT0sabOkpyVdN0qcyyQ9JWlY0jpJpyXOHS7pq5KekbRFUiUOv03SDkmb4r+LEtcsk/S8pK2SLkgrsQ0qleiFVM2CMrJgIXPuvIVL7lzIsmotUxPSOUgoLstQnmNCajW2XbOCldd0vgyPtgmj010a7dmSJWFsdNuKpbIfuMHMnpA0C9goaa2ZbU7EeQk4z8xek3QhMACcHZ+7HXjQzC6R1Av8QuK6vzSzv0h+maSTgUuBU4B3AN+SdLyZjUwqha3iW7AESyguy1CeY1wS+9NdQi8X3TPIig6VY7fcwiKUad8TioqZ7QR2xse7JW0B5gCbE3HWJS5ZDxwDIOktwPuAy+N4daA+wVdeDHzJzPYCL0l6HjgL6GxfyLdgCZZQKksozzEuXewclWGVfYiMN64XwrTvtsZUJB0LnAE8Pk60K4AH4uPjgB8BfxO7xDYC15lZY7vlayUtATYQWUOvEQnW+sT9tsdhzc+yFFgK0NfX104yRse3YAmaECpLSM8xJl3sHOXCcisYebAOWxYVSYcBq4Hrzez1MeIsIBKVcxP3fw9wrZk9Lul24CbgFuDzwHLA4v+fAT7W6vOY2QCRm43+/v503onsW7A4eaeLnaNcWG4FI33rcNbMdJ7sIC2JiqQZRIKyyszWjBHnVOAu4EIz2xUHbwe2m1nDsvkqkahgZj9MXPsF4B/jjzuAuYlbHxOHOU0UeXprx9NW5Jl+XewcBW+5ZUAny2761uEvzpr6UzVhZuP+AQJWAp8bJ04f8Dzw70Y5923ghPj4NuDT8fHRiTh/QDSOAtEA/ZPAIUTusxeBnvGe8cwzz7SysW6d2aGHmvX0RP/Xrcv6idKj42lbt872H3Ko7aPH9nConde7rlD552RHN+rlunVm//W/pnNvmLXFJtCAdv9asVTOARYDw5I2xWE3x0KCmd0J3ArMBu6QBLDfzPrjuNcCq+KZXy8CH43DPyXpdCL318vAVfH9npZ0H9FEgP3Ax63TM79SphsWRJEHSTuRtjf8Jj7Tz+kQ3aiX6VqHGbxO2MweI7JWxotzJXDlGOc2Af2jhC8e535/BvzZRM8WIt0aSCvyIGnaaWv+TR7/XJWTfaafMwHDAzV2rR5i9qIq85e2VomLXC9bpXQr6jtJrQa33QZ798KBA521IIo8SJp22pp7j/+4q8J8n+nnjMPwQI13XbWQk6hTf6iXYQZbEpYi18tWcVFJiUZvuCEo06Z1vqdS5EHSNNM2au/RZ/o547Br9RAnUWd67CLdtXoIWrRWilwvW6F0uxR3ikZvuCEo558f5hzy0Sj6VhuN3uPy5fn5TZxsmb2oSp1e9tHDPnqZvaia9SPlBrdUUqK5N3zbbflovPKwmCoNyt57dNpj/tIKwwy2PabiuKikRm58qU1T04o8i8xxpsL8pZWWXV7OQVxUUiT43vAoZkm1WintbJUiLx51nKxwUSkTo5gllWWVfFhYKVMWt5/jdBsXlTKRGPgZmd7Lqu9XmVfLgYXVAdzt5zidwWd/lYl44Gfb7y1noQ3ysS9UWLiwuLO+xqOhr1m/0MhxioZbKmWjUuGLQxUeGyl3Lz3UiRU+zuPkHReVEuJbSUSE5vbzcR6nCLj7q4T4YsDW6Pai0NHGeRwnb7ilUlJC66V3hCm8MyULq8EtSKcIuKg4xaRWY2TBQubsrXMJvVx0zyAr2tjePovZYaGO8zhOO7ioOMVkiu9MycpqKIUF6RQaF5UW8Vk5OaNaxabwzhS3GhxncsiiV/jmmv7+ftuwYUPH7u+zcnJKkd9D7zgpIGlj4i29qeCWSgv46uuc4u9MCQK38svFhFOKJc2V9LCkzZKelnTdKHEuk/SUpGFJ6ySdljh3uKSvSnpG0hZJlTj803HYU5Lul3R4HH6spP8naVP8d2eaCZ4Mvvr6IEV/94qTLg0r/5ZbKO3uDWWjFUtlP3CDmT0haRawUdJaM9uciPMScJ6ZvSbpQmAAODs+dzvwoJldIqkX+IU4fC2wzMz2S/pzYBlwY3zuBTM7fYppSw33r0e4G9BpF7fyy8eEomJmO4Gd8fFuSVuAOcDmRJx1iUvWA8cASHoL8D7g8jheHajHxw81XXPJFNLRcXxWjjcQTvt0axadu9jCoa0xFUnHAmcAj48T7Qrggfj4OOBHwN/ELrGNwHVmtqfpmo8BX058Pk7SJuBnwJ+a2bdHeZalwFKAvr6+dpLhTBJfnFdMOtkgd8PKdws6LFoWFUmHAauB683s9THiLCASlXMT938PcK2ZPS7pduAm4JbENX9C5GJbFQftBPrMbJekM4G/l3RK83ea2QCRm43+/v78T2HrEGk2GO4GLB7daJA7beW7BR0WLYmKpBlEgrLKzNaMEedU4C7gQjPbFQdvB7abWcOy+SqRqDSuuRz4ELDQ4rnNZrYX2Bsfb5T0AnA80Lk5wwWlEw2GuwGLRREaZLegw6KV2V8C7ga2mNlnx4jTB6wBFpvZs41wM3sF+IGkE+KghcRjMZI+CPwx8Otm9n8T9zpCUk98/E5gHvDiJNJWenyDwnzSzRl2uZnZWKux7ZoVrLym9qZ88Q1Sw2LCxY+SzgW+DQwDB+Lgm4E+ADO7U9JdwCJgW3x+f2NBjaTTiSyYXiJx+Gg8S+x54BCgYdWsN7OrJS0CPgHsi7/vP5vZP4z3jJ1e/JhX3NecP7L4zbo9yN3298X7uNneOnV6uai3vX3cnLHJZPGjmT0GaII4VwJXjnFuE/CmhzazXxkj/moiV5szRXwMJH9ktZFlt8rGpERzivu4Od3FV9QXHB8DyRdFHx+YlGhOcR83p7u4qDhOQIRkXQ4P1Ni1eojZi6rMX5rOg0xKNCsVeh4eZHu8j9sK38ctaHxDScdx3sTwQI13XbWQXqJxjBf+ejA1YfGFiuHgG0o6jtMVdq0e4iTqTI/HMXatHoKURMVdssXG31HvOB0g7xtvzl5UpU4v++hhH73MXlTt+jPkPQ/LilsqjpMyRZjKPX9phWEGUx1TacftVYQ8LCsuKo4zScZqJIuwSh0iYUnL5dWWSNRq7L1tiPfsrfKdA5Vc52EZKZyo+CCg0w3GaySLNi04jTrVstDGGXve3joPHejlfAZ5Ynol93lYJgolKm4yO91ivEYypGnBUyWtOtWy0MYZqwMjzKBOlSE2Wo4zsIQUSlSK4nZwwmeiRrIoM5zSqlMtC22csSP/Wmef9fIwVUZGvC7niUKJStHcDk645MEaScNtlWadaklo44zdvnKIj9xT5Z9HKl6Xc0bhFj+WZUylLOl0JkearuCsylrIZTzkZ2sHX/zYAkVxO4xHyGNHRalseSdNV3BWdSrUuhxy/QsBX/yYQ0J9T0qjst1yS/TfF61lRzffkzI8UGPoghUMD5TjBw+1/oVC4SyVMhDq2JFPlAiHbo35NPYIO4k69Yd6GSa9PcJCJdT6FwouKjkk1EFir2xh0Q33UZp7hOXFdRpq/QsFF5WcEqK/2Stb+Zi9qEr9oV6M+pT2CAtpnKIVcQux/oWCi4qTKl7ZykVae4SF4joNSdzyyoQD9ZLmSnpY0mZJT0u6bpQ4l0l6StKwpHWSTkucO1zSVyU9I2mLpEoc/jZJayU9F/9/a+KaZZKel7RV0gUTPeOePa0n2HGcdJm/tEL1m8umNJbSzYkF4+GD8FOnldlf+4EbzOxk4FeBj0s6uSnOS8B5ZjYfWA4MJM7dDjxoZicCpwFb4vCbgEEzmwcMxp+J730pcArwQeAOST3jPeDu3S2kwnGcYGm4Tpcvz9Y6CEXc8kzbix8lfQ34KzNbO8b5twLfM7M5kt4CbALeaU1fJGkrUDWznZKOBobM7ARJywDMbEUc75vAbWY25nzFk07qty1b/M2PjuNMnbxMGEiDzBc/SjoWOAN4fJxoVwAPxMfHAT8C/iZ2iW0ErjOzPcBRZrYzjvcKcFR8PAdYn7jf9jis+VmWAksB+vr62klGqQqN4zjt4eOCU6PlxY+SDgNWA9eb2etjxFlAJCo3xkHTgfcAnzezM4A9xG6uJLEV05bJZGYDZtZvZv1HHHFEy9f5Aj3HcZzO0ZKoSJpBJCirzGzNGHFOBe4CLjazXXHwdmC7mTUsm68SiQzAD2O3F/H/V+PwHcDcxK2PicMmTfK1pD4Q5ziO0zlamf0l4G5gi5l9dow4fcAaYLGZPdsIN7NXgB9IOiEOWghsjo+/DnwkPv4I8LVE+KWSDpF0HDAP+G5bqUrQbJnMnu0DcY7jOJ2ilTGVc4DFwLCkTXHYzUAfgJndCdwKzCaaqQWwPzH4cy2wSlIv8CLw0Tj8k8B9kq4AtgEfju/3tKT7iMRnP/BxMxuZbAKbLZNdu3yBnuM44ZPXsd/CbX3fjC9mejN5LayOUxa61W5lPvsrj/jWIQlqNbatHGLZPVUei19+5CLrONkxPFAbdTeCUHYYmAyFFxXwKYLAz7s+c/+1zjesl4UM8s/1Sq4KqzM6qViecYfjEarMW1LxMtEFxtvhOc+bs5ZCVBx+3vWZZiPMoM77NcSTvZVcFdZJU+AGMxU3Sa3GyIKFzNlb5xJ6ueieQVYMFSufQmS8HZ7z7GHxl3SVhcT+E9MO6eWEq6rlcH01Gsw7b+GSOxeyrFor1NqkyUyRT06xb9xE9ahxm0Gdc/YN+VT7LjB7UZU6veyjZ9QdnisVWLYsf3XULZWykOj69FSrLMlbSZ0scYPZE/cGowaz9V546JMa2nWTjGrZVKtYby/79kbb139nRpUVE9zHmTpp7fAcGi4qZWKCwaXQG9BJ0dRgPja9ysnfj9I6URrzMHOwXTfJqAPAyyr0PDzI9thFuKJgLsKQmb+0MumXmoWKi4oDhNmApiJylYMN5pdfqbL+GxW+8wW4996J05iXGTjtTEQZ07KpVPjlSoUlHXrGqVDIzk4b5C39LioBkWXhCa0BTVXk4gZzZAWM/EPraWzFtZS3Cp+3AeAQOzvdpNX0h1QOXVQCIevKE9oUxk6IXLtpnKgBzvo3myx5mmIfWmen27SS/tDKoYtKIGRdeULrwXZC5CaTxvEa4Kx/s1YJqRfbLqF1drpNK+kPrRy6qARCCJUnpB5sp0QuzTSG8JtNRGi92HYJrbPTbVpJf2jlsPB7f+WJPPcoy0rov9mKFdEO3SMj0c7cy5dHax+cYjHZctiJvb9cVBynwOTdUnE6i28oGTCh91g7SZnTHjpldx853cdFJQXK3Bssc9rzQkhjZU7x8b2/UqATryh+0/5MgeKvZ3YcJ4lbKimQ9uyLPPX+Q5t5Mhncfec46eGikgIt+a3baLlCm3c+Hnn32edJwB0nD0woKpLmAiuBowADBszs9qY4lwE3AgJ2A9eY2ZPxuZfjsBES766X9GXghPgWhwM/NbPTJR0LbAG2xufWm9nV4z3jnj0TpaLzjOu3jrdfV72O9fbS8/D4LVfeev959tnnScAdJw+0YqnsB24wsyckzQI2SlprZpsTcV4CzjOz1yRdCAwAZyfOLzCzHydvama/0ziW9BngZ4nTL5jZ6a0mYuvW1nadzYptK4eYszfafn3f3jrbVw7xy+M8bN57/3miWwLuLjanLEwoKma2E9gZH++WtAWYA2xOxFmXuGQ9cEyrDyBJwIeB97d6zZufMewe5iNUuYRejGj79UeoTrgbbJ57/3miGwLuLjanTLQ1+yt2TZ0BPD5OtCuABxKfDfiWpI2Slo4S/73AD83suUTYcZI2SXpE0nvHeJalkjZEfxa0i2jekgoX9Q5ym5ZzUe8g85Z4ixISnX7Dns+Qc8pEyyvqJR0GPAL8mZmtGSPOAuAO4Fwz2xWHzTGzHZKOBNYC15rZo4lrPg88b2afiT8fAhxmZrsknQn8PXCKmb0+1rOddFK/bdkS9op6d3+UF7dUnFDJbEW9pBnAamDVOIJyKnAXcGFDUADMbEf8/1VJ9wNnAY/G10wHfgs4MxF/L7A3Pt4o6QXgeGBM1Zg5s5VUZEsa7iwXpnziY2ROmWhl9peAu4EtZvbZMeL0AWuAxWb2bCJ8JjAtHouZCXwA+ETi0vOBZ8xse+KaI4CfmNmIpHcC84AX209asfDebr7xMTKnLLRiqZwDLAaGJW2Kw24G+gDM7E7gVmA2cEekQT+fOnwUcH8cNh34opk9mLj3pcDfNX3f+4BPSNoHHACuNrOfTCJthcKnvjpOuAwP1Ni1eojZi6rRe+dLTCuzvx4jWn8yXpwrgStHCX8ROG2c6y4fJWw1kavNSZC3tSuOUxaGB2q866qFnESd+kO9DDNYamHxvb9yQsMvv3y5u74cJyR2rR6ilzrTGWEGdXatHsr6kTLFt2nJEXn3y/tEA6cd8lJeZi+qUn/o4Dq02YuqWT9SprioBESR/bI+0cBphzyVl/lLKwwzWNi62y4uKoFQdL+sTzRw2iFv5WX+0goUqL5OBR9TCYSi+2UbEw16enyigTMxXl7yi1sqgVB0v6wvAHTaIevykpfxnBBpeZuWkOnv77cNG8LepqUV8jam4hXPKSJ5Gs+Zah3MbJsWpzt00y871cKYp4rnOO2Ql/GcUOugi0oOyUoQkt+bl4rnOO2Sl4XGodZBF5WckUbvZDKFsfl7P/e5N1a8D82use2aIR6hyrwllSAKt+NMhqzHc1olVPFzUckZafROJlMYm793166DFe9Ds2uc/B8XYnvrXEIvF90zyIohF5ZW8bGp8MjDQuNQxc9FJWek0TuZTGEc7Xt/XvFWDHGgHr0u2ahzzr4hhlxUWiJUv7iTD0IUPxeVnJFW76Tdwjju91arWG8v+/ZG06G/M6PKiurknqtshOoXd5zJ4lOKnXSo1di2MuwxlRDdTG6pOFnSiSnFLipOKQi58Q5R7DpNGdMcIr5OxXEmSchuphD94p0kZIF3pk4h9v565ZWooDqdpVaDFSvymde+l1Q4jCbwTnEohKWyY0fU8/EeT+fIe+8y1OmXZSTU9RVOOkxoqUiaK+lhSZslPS3pulHiXCbpKUnDktZJOi1x7uU4fJOkDYnw2yTtiMM3SboocW6ZpOclbZV0QSsJ8R5PZylC77JSgWXLXFCyxt9iWmxasVT2AzeY2ROSZgEbJa01s82JOC8B55nZa5IuBAaAsxPnF5jZj0e591+a2V8kAySdDFwKnAK8A/iWpOPNbGS8h/QeT2fpZu/SB3GLT9nGkcrEhKJiZjuBnfHxbklbgDnA5kScdYlL1gPHTOGZLga+ZGZ7gZckPQ+cBYzpyZ8zB77yFS+knaTHgytUAAARtElEQVRb7qO8u9kcp+y0NVAv6VjgDODxcaJdATyQ+GxE1sZGSUub4l4bu83ukfTWOGwO8INEnO1xWPOzLJW0QdKGnp4fjd/w1Gpsu2YFK6+p5XKQORS64T4qgpvNSY/hgRpDF6xgeMArbl5oeaBe0mHAauB6M3t9jDgLiETl3ETwuWa2Q9KRwFpJz5jZo8DngeVEorMc+AzwsVafx8wGiNxs9Pf3j73YplZjZMFC5vi+VEExlouriIO47s6bHEV/xXZRaUlUJM0gEpRVZrZmjDinAncBF5rZrka4me2I/78q6X4iV9ajZvbDxLVfAP4x/rgDmJu49TFx2OQYGkK+L1VQjOfiKtosLXfnTZ5dq4c4KX7FtjVese2iEjytzP4ScDewxcw+O0acPmANsNjMnk2Ez4wH95E0E/gA8L3489GJW/xmIxz4OnCppEMkHQfMA77bbsJ+TmNfKnp+vi9VEXq/eWYiF1eRZmm5O2/yzF5Upc7Bulu0V2wXlVYslXOAxcCwpE1x2M1AH4CZ3QncCswG7og0iP3x0v+jgPvjsOnAF83swfgen5J0OpH762Xgqvh+T0u6j2giwH7g4xPN/BqXSoWehwfZHu9LtSLQfanKRBFdXGNRprSmzfylFYYZzNUrttMmj65T3/vLyYQ8VpbJUqa0OunRDdep7/3lFIYyrVMoU1qd9Ah5v7rxKMTeX45TZPK855ozefK6X13hLBV3NThpk2WZ8tlj5SWvMyELJSpeAYtBSB2DrMtUXl0gTvuMVu7z6DotlKh4Bcw/WTfizWRdpnz22EFC6mykTWjlfioUSlRCrYBFrgxpk3Uj3kzWZSqvLpC0KVKjOxrN5X7lyvz+5oUSlRArYNErQ9pk3Yg3E0KZ6oQLJG8dndA6G2mTLPfTp8M990RpzWObUShRgfB8kEWvDGkTQiPeTGhlaqrksaMTWmcjbZLl/vvfhy98Ib9tRuFEJTSKXhk6wWiNeN561iGTx45OiJ2NtGmU+1oN7r03v22Gi0qHKUNl6DR57FmHTF47OkWzGMci722Gi0oX6GZlKGKPPo896zTo1G+Z90arDORZQF1UCkRRe/R57VlPhU7/lnlutCZDETtboeKiUhBqNbjtNti7Fw4cyFePfqIKX8aedVmts05Q1M5WqBRWVMrUM2lUmoagTJuWnx59qxW+bD3rVqyzMpXxqeAC3V0KKSpl65k0Kk1DUM4/P7Ja8pBmr/CjM651VquxbeUQy+6p8thIpRRlfCqU0X2aJYUUlbI1VM2VJi+CUqtFc/J7eqLPXuHfyKjWWdxjmvuvdb5hvSxkkH+uVwpfxqdCGd2nWVJIUSlbzySPlWZ4oMYDHx9i+EAVzajwe78HS5bk49kzJe4xTbMRZlDn/Rriyd5K4cv4VAnVfZqVC7PxvTBrZtr3LqSo5LGRnSqhVppRqdU48fcXcuv+OjfRywf2DdLXl91rnjtdsVO9f6LHNG16Lyd8tMpgF8TYx2/SJys3ffJ7Yd7xad9/QlGRNBdYSfS+eQMGzOz2pjiXATcCAnYD15jZk/G5l+OwEQ6+ux5Jnwb+A1AHXgA+amY/lXQssAXYGt9+vZld3W7CctXI5pApNTJDQ0wfqSNGMOq8f9oQ1Wo2P1anK3bq90/0mHqqVZZ0oZCXbYyyW2Tlpk9+L1GbnSqtWCr7gRvM7AlJs4CNktaa2eZEnJeA88zsNUkXAgPA2YnzC8zsx033XQssM7P9kv4cWEYkTAAvmNnpk0qR03Gm3MhUq+iQXmxvHXp6+e2/qjI/o0aq0xW7I/fvco+pbGOU3SIrN33ye0dGsLTvP6GomNlOYGd8vFvSFmAOsDkRZ13ikvXAMS3c96Gmay5p8ZmdjJlMI/NGyybqbWtoiBnVKvMzbKE6XbGLML5XhDSESFZu+uT33nzzc8+mfX+ZtS5UsWvqUeDdZvb6GHH+CDjRzK6MP78E/IzI/fXXZjYwyjX/AHzZzP42/o6ngefi6/7UzL493nP19/fbhg0bWk6HMzXatVRCd5/kakwlI4qQBufNSNrYGJJIi5YH6iUdBqwGrh9HUBYAVwDnJoLPNbMdko4E1kp6xsweTVzzJ0QutlVx0E6gz8x2SToT+HtJpzR/p6SlwFKAvr6+VpPhpEC7PazQ3Sed9iYVYXyvCGloFxfSydGSqEiaQSQoq8xszRhxTgXuAi40s12NcDPbEf9/VdL9wFlE1g6SLgc+BCy02GQys73A3vh4o6QXgOOBN5giscUzAJGl0mJ6nZRop5Fx94mTN0K3rkNm2kQRJAm4G9hiZp8dI04fsAZYbGbPJsJnxoP7SJoJfAD4Xvz5g8AfA79uZv83cc0Rknri43cC84AXJ5c8JwQals3y5V45nXwwmnXttEYrlso5wGJgWNKmOOxmoA/AzO4EbgVmA3dEGvTzqcNHAffHYdOBL5rZg/E9/go4hMglBgenDr8P+ISkfcAB4Goz+0mrCXKTNUzK6D4pI0Wpf2la10XJk1Zpa6A+VBoD9W6yOk52FK3+pSEGoedJJwbqJ3R/5Qk3WR0nO7pV/2o1WLEi+t9JKhVYtmxqIlDGNqlQ27T4gLDjZEc36l/oPf9mytgmFUpUyrjnV5aUzVfsjE836l/o09ObKWObVKgxFad75K3HWHTKIvBe7tIl08WPTrhk0aDkrcdYZMrU0Jax5583XFRyTlYNShl9xaFSNoH36elh46KSc7JqULzHGA4u8E5IuKjknCwbFO8xhoELvBMShRWVsgxceoPigAu8Ew6FFJUyDVyCNyiOk3eK1AkupKh0Y5yhSIWg6Phv5YRM0TrBhRSVTo8zFK0QFBn/rcpLXjoTRZu9V0hR6fQ4Q9EKQZHx36qc5KkzUbTZe4UUFUh/nCHZ6ylaISgy/lulR156/pCvzkTRJtsUVlTSZLReT5EKQZEpWoXNijz1/CF/nYkiTbYphKi88kpU6Dv1o4zW65nqlthO9yhShc2KPPX8wTsTWVIIUdmxI+pFdar3lLdej+OkTR7rgHcmsqEQogKd7T15r8cpO14HnFYpxNb3Ur8deuiG4P28juM4IZHJ64QlzZX0sKTNkp6WdN0ocS6T9JSkYUnrJJ2WOPdyHL5J0oZE+NskrZX0XPz/rYlzyyQ9L2mrpAsmesY5c8IfOAyNbr2S1Ukf/+2ckGnF/bUfuMHMnpA0C9goaa2ZbU7EeQk4z8xek3QhMACcnTi/wMx+3HTfm4BBM/ukpJvizzdKOhm4FDgFeAfwLUnHm9nIWA/49re7oLRD3mbyFJXJTNH1384JnQlFxcx2Ajvj492StgBzgM2JOOsSl6wHjmnhuy8GqvHxvcAQcGMc/iUz2wu8JOl54CzA+2UpkbeZPEVkNHGAiUXGfzsndNoaU5F0LPAo8G4ze32MOH8EnGhmV8afXwJ+BowAf21mA3H4T83s8PhYwGtmdrikvwLWm9nfxufuBh4ws682fc9SYGn88d3A91pOSLH5JaDZKmxi1kyYdzwgwOC5Z2H3ni48W7dpIS+yYs7b4e1zDn7+yY/grbNB08AOjP2bTPq3CzgvZs2EX5wFr+/uUjkMOC+6zglmNivNG7Y8+0vSYcBq4PpxBGUBcAVwbiL4XDPbIelIYK2kZ8zs0eR1ZmaS2poxEItTQ6A2pD3YlFc8Lw7ieXEQz4uDeF4cJDnOnRYTDtTHXzyDSFBWmdmaMeKcCtwFXGxmuxrhZrYj/v8qcD+RKwvgh5KOjq89Gng1Dt8BzE3c+pg4zHEcxwmcVmZ/Cbgb2GJmnx0jTh+wBlhsZs8mwmfGg/tImgl8gINuqq8DH4mPPwJ8LRF+qaRDJB0HzAO+227CHMdxnO7TivvrHGAxMCxpUxx2M9AHYGZ3ArcCs4E7Ig1if2xeHgXcH4dNB75oZg/G9/gkcJ+kK4BtwIfj+z0t6T6iiQD7gY+PN/MrZqCFdJQFz4uDeF4cxPPiIJ4XB0k9Lwqx+NFxHMcJg5bGVBzHcRynFVxUHMdxnNQIUlQkfTDeouX5eLV983lJ+m/x+ackvScOH3NLmfG2hQmZDuXFpyU9E8e/X9Lh3UzTZOhEPiSuvUGSSfqlbqRlqnQqLyRdG5eLpyV9qlvpmQodqh+nS1qveGspSWc13zdEppAX/0bSdyU9GefFf0lc0367aWZB/QE9wAvAO4Fe4Eng5KY4FwEPEC0A+1Xg8Tj8aOA98fEs4NnGtcCngJvi45uAP886rRnmxQeA6fHxn4eeF53KhzhsLvBNoskiv5R1WjMsEwuAbwGHxJ+PzDqtGebFQ8CFieuHsk5rh/NCwGHx8QzgceBX489tt5shWipnAc+b2YtmVge+RLR1S5KLgZUWsR44XNLRZrbTzJ6AaEsZoLGlTOOae+Pje4Hf6HRCUqAjeWFmD5nZ/vj6VrfVyZJOlQmAvwT+GMjLjJVO5cU1wCct2h4Ji9aVhU6n8sKAX4yP3wL8704nJAWmkhdmZv8njjMj/rPENW21myGKyhzgB4nP23ljI9BSHEVbypxBpLoAR1m0jxnAK0TTnUOnU3mR5GNEvZeQ6Ug+SLoY2GFmT6b7uB2lU2XieOC9kh6X9Iikf5viM3eKTuXF9cCnJf0A+AtgWWpP3DmmlBeSehQtGXkVWGtmk243QxSVKaMJtpSxyJbLS890SoyXF5L+hGgt0Kosnq2bNOeDpF8gWm91a7ZP1n3GKBPTgbcRuUX+E9EaMmX0iF1jjLy4BvgDM5sL/AHR4u9CY2YjZnY6kdfiLEnvHiVOS+1miKLSyjYtY8bR2FvKjLUtTMh0Ki+QdDnwIeCyuLCETCfy4V3AccCTkl6O4z8h6e2pP326dKpMbAfWxK6Q7wIHiDZeDJlO5cVHiHYIAfgKB7eWCpkp5UUDM/sp8DDwwTio/XYzrYGitP6IekwvElX4xoDTKU1x/j1vHHD6rh0ccFoJfG6U+36aNw44fSrrtGaYFx8k2rHgiKzTmGU+NF3/MvkYqO9Umbga+ER8fDyRm0RZpzejvNgCVOPjhcDGrNPa4bw4Ajg8Pj4U+Dbwofhz2+1m5pkxRgZdRDQb4wXgT+Kwq4GrEwXif8Tnh4H+OPxcIvPsKWBT/HdRfG42MAg8RzTL5W1ZpzPDvHg+bjQa4Xdmnc4s8qHp/i+TA1HpYJnoBf6WaG++J4D3Z53ODPPiXGAjUcP8OHBm1unscF6cCvxLnBffA25N3LPtdtO3aXEcx3FSI8QxFcdxHCenuKg4juM4qeGi4jiO46SGi4rjOI6TGi4qjuM4HUDSb8cbNB6Q1D9OvJclDTc2sEyEj7qZo6Rfk7QxvmajpPe38Cyr4s0mvyfpnniNTkdwUXEcx5kikqqS/mdT8PeA3wIebeEWC8zsdIvemNvgJmDQzOYRTett7Dz8Y+A/mNl8ooWa/6uF+68CTgTmE61FubKFayaFi4rjOE4HMLMtZrZ1CrcYdTNHM/sXM2tscvk0cKikQwAkfUBSTdITkr4Sb0ODmX3DYoDv0sFNZF1UHMdxssWAb8WurKWJ8FY2c1wEPGFmexW9D+hPgfPN7D3ABuAPk5Fjt9di4MG0E9Fgeqdu7DiOU3QkPQ4cAhwGvC3e6RfgRjP7Zou3OdfMdkg6Elgr6Rkze4PLzMxM0htWqks6heh9SB+Ig34VOBn4TrwXaC9Qa/quO4BHzezbLT5b27ioOI7jTBIzOxuiMRXgcjO7fBL32BH/f1XS/UQbWD5KvJmjme1s3sxR0jHA/cASM3uhEUy0bf3vjvY9kv4z0T5fV7X7jO3g7i/HcZyMkDRT0qzGMZHV8b349NeJBuKJ/38tjnc48E9EGz1+J3G79cA5kn4lce/j4+MrgQuA3zWzA51Mk4uK4zhOB5D0m5K2AxXgnyR9Mw5/h6RvxNGOAh6T9CTRAPo/mVljvOOTwK9Jeg44P/4M8PvArwC3xtOQN0k60sx+BFwO/J2kp4hcXyfG19wZf1ctjt+x9wj5hpKO4zhOaril4jiO46SGi4rjOI6TGi4qjuM4Tmq4qDiO4zip4aLiOI7jpIaLiuM4jpMaLiqO4zhOavx/1bSMLXF2X5gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5e948db5f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sanity check plot overplotting the HSC parents atop the COSMOS catalog:\n",
    "plt.figure()\n",
    "#plt.scatter(df_cosmos['ra'], df_cosmos['dec'], marker='.', s=1)\n",
    "plt.plot(df_cosmos['ra'], df_cosmos['dec'], 'b.')\n",
    "#plt.hist2d(df_cosmos['ra'], df_cosmos['dec'], bins=(300, 300), cmap='Blues')\n",
    "plt.plot(parents_ra,parents_dec,'r.')\n",
    "#plt.colorbar(label='Binned Counts')\n",
    "# NOTE: these coordinate ranges were chosen arbitrarily to zoom in on a small region:\n",
    "# Comment out the following two lines to see the whole COSMOS field (plus the HSC patch).\n",
    "plt.xlim(150.22,150.23)\n",
    "plt.ylim(2.25,2.27)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Picking a single parent (arbitrarily), so we can compare its family to COSMOS:\n",
    "pra = np.array(parents_ra)\n",
    "pdec = np.array(parents_dec)\n",
    "#pick_parent = (pdec > 2.255) & (pdec < 2.2575) & (pra > 150.222) & (pra < 150.224)\n",
    "#pick_parent_ind = [i for i, x in enumerate(pick_parent) if x][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Picking a single parent by specifying an RA, Dec and finding the closest parent to those coords.\n",
    "# ...so we can compare its family to COSMOS:\n",
    "desired_ra = 150.2229\n",
    "desired_dec = 2.25596\n",
    "total_sep = lambda target_ra, target_dec, candidate_ra, candidate_dec: ((target_ra - candidate_ra)*np.cos(np.radians(target_dec))) ** 2.0 + (target_dec - candidate_dec)**2.0\n",
    "pick_parent_ind = np.argmin([total_sep(desired_ra, desired_dec, pra[i], pdec[i]) for i in range(len(pra))])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "150.22291945562662 2.2559620121522577\n"
     ]
    }
   ],
   "source": [
    "print(pra[pick_parent_ind],pdec[pick_parent_ind])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#parents[pick_parent_ind]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What does the deblending look like for this source?\n",
    "deblender.plotDeblendFamily(coadd, families[pick_parent_ind][0], families[pick_parent_ind][1], display=disp2,symbolSize=25,plotb=True,background=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[150.22295034434663, 150.22279188143824, 150.2226893524708, 150.2231966152749] [2.2559250230721903, 2.256251958656348, 2.2553479670254704, 2.2566550472399265]\n"
     ]
    }
   ],
   "source": [
    "# Picking the children of our selected parent:\n",
    "cra = [child['coord_ra'].asDegrees() for child in families[pick_parent_ind][1]]\n",
    "cdec = [child['coord_dec'].asDegrees() for child in families[pick_parent_ind][1]]\n",
    "print(cra,cdec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEJCAYAAABc/7oDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X94VdWd7/H3V4IW8EeqILWEgFJLR0cFm0osQgu2IijVTm9tK49UIVJ7gwp2poj2On1up4zOtLZWab2Y6MhcaQcVp8wUFIvewU6ICogiv5Qo0FCo6EB1isWkfu8fax88CedXwvmZ83k9T55zztpr77POyUm+Z+/1XWuZuyMiIpINRxW6ASIi0nMoqIiISNYoqIiISNYoqIiISNYoqIiISNYoqIiISNYoqIiISNYoqIiISNYoqIiISNZUFLoBXdG/f38fOnRooZshIkXs4EHYsgU+9jHo1+/w7X/8I2zbBp/4BBxzTP7bVwhr1659090H5OO5SiqoDB06lDVr1hS6GSKSIy0tcM89sGgRvPkm9O8PV14JM2fCsGGZH2f5cpg6FS6/HOrqoLoadu6Ehobws3QpTJyYu9dRbMxsR76eS5e/RKQoLF8OtbXQpw80NYUzjqam8Li2NmzP1MSJ0NwcjjF6dDjG6NHhcXNzeQWUfLNSmlCypqbGdaYi0vO0tITAsXQpnH/+4dtXr4YvfCEEhK6csUhgZmvdvSYfz6UzFREpuHvugWuvTRxQIJTX1cH8+fltl3SdgoqIFNyiRTB9euo6dXWhnhQ3BRURKbg334QhQ1LXqa4O9aS4KaiISMH17w870uQn7dwZ6klxU1ARkYK78kpobExdp6Eh1JPipqBSQlpaYPZsGDgQevUKt7Nnh3KRUjZzJtx3X8jySmT16hBU6uvz2y7pOgWVEpHNHH6RYjNsGCxcGNKG584NX5Ta2sLt3LmhfOFCpROXgrRBxcwGm9nTZrbJzDaa2Y0J6kwxs5fMbIOZNZnZOXHbKs3sETPbYmabzez8Tvt+y8zczHS1NImWljA6eOlSmDcv/GFVVITbefNC+dSpOmOR0qYBiz1D2sGPZnYKcIq7rzOz44C1wOXuvimuzqeBze6+z8wmAt9191HRtgeBZ9y9wcyOBvq6+/5o22CgAfgE8El3T5nbUa6DH2fPDn9g8+YlrzN3bvjju/PO/LVLREpDUQ1+dPfd7r4uuv8OsBkY1KlOk7vvix42A1UAZnYCMBZojOq9FwsokR8B3wZKZ1h/ASiHX0RKRZf6VMxsKDASeDZFtelA7Ar/qcBe4AEze8HMGsysX3Ssy4Bd7v5iVxtdbpTDLyKlIuOgYmbHAo8Cs9z97SR1xhGCypyoqAI4F/iZu48E/gjcbGZ9gVuA2zJ43hlmtsbM1uzduzfT5ha9rmRyKYdfREpFRkHFzHoTAspD7r4kSZ2zCf0jl7n7W1FxK9Dq7rEzm0cIQWYY4SzmRTPbTrhcts7MPtL5uO6+wN1r3L1mwIC8LAeQc13N5FIOv4iUikw66g14EPgvd5+VpE418BQw1d2bOm17Bqhz961m9l2gn7v/Tac624Gacuio785srJrBVUSORFF11AOjgauA8Wa2PvqZZGbXmdl1UZ3bgJOAn0bb4//zXw88ZGYvASOAFDlMPV93ZmNVDr+IlAqtp5JnAweGS12pAkBLS8jP37Pn8PL58w9fFa++XgFFRJLL55mKgkqe9eoV+lAqUizk3NYW+lfa2/PXLhHpuYrt8pdkkTK5RKQnU1DJM2VyiUhPluIijOTCzJkhk2vy5OSZXA0NIZNLRKTU6Ewlz5TJJSI9mYJKAWg2VhHpqZT9JSLSwyn7S0RESpKCioiIZI2CioiIZI2CihSNriwHICLFSUFFikJXlwMQkeKk7C8pOE3tL5Jbyv6SstKd5QBEpDgpqEjBLVoE06enrlNXF+qJSHFTUJGCe/NNGDIkdZ3q6lBPRIpbWQcVZRsVBy0HINJzlG1QUbZR8dByACI9R1lmfynbqLjo9yGSW8r+yjFlGxUXLQcg0nOkDSpmNtjMnjazTWa20cxuTFBnipm9ZGYbzKzJzM6J21ZpZo+Y2RYz22xm50fl/xiVvWRmj5lZZXZfWnLKNio+Wg5ApGdIe/nLzE4BTnH3dWZ2HLAWuNzdN8XV+TSw2d33mdlE4LvuPira9iDwjLs3mNnRQF93329mFwFPuXu7md0B4O5zUrUlW5e/evUK/6wqUqx72dYW/rG1tx/x04mIFFRRXf5y993uvi66/w6wGRjUqU6Tu++LHjYDVQBmdgIwFmiM6r3n7vuj+yvcvb3zPvmgbCMRkdzoUp+KmQ0FRgLPpqg2HYjlTp0K7AUeMLMXzKzBzPol2Gda3D45p2wjEZHcyDiomNmxwKPALHd/O0mdcYSgEruMVQGcC/zM3UcCfwRu7rTPrUA78FCSY84wszVmtmbv3r2ZNjelmTPhvvtCVlEiq1eHoFJfn5WnExEpGxkFFTPrTQgoD7n7kiR1zgYagMvc/a2ouBVodffYmc0jhCAT2+dq4FJgiifp3HH3Be5e4+41AwYMyKS5aSnbSEQkNzLJ/jJCn8hmd78zSZ1qYAlwlbu/Eit39z3Ab81seFR0IbAp2udi4NvAF9z9wBG9im5QtpGISPZlkv11AfAMsAF4Pyq+BagGcPd7zawB+BIQ6/5uj2UamNkIwhnM0cBrwDVRltg24BggdlbT7O7XpWqLpr4XEem6fGZ/pUiqDdz9N4ClqVMH1CXZth447MW4+8cybKOIiJSIshxRLyIiuaGgIiIiWaOgIiIiWaOgIiIiWaOgIiIiWaOgIiIiWaOgIiIiWaOgIiIiWaOgIiIiWaOgIiIiWaOgItIFLS0wezYMHBhWEB04MDxuaSl0y0SKg4KKSIaWL4fa2jCjdVNTmNG6qSk8rq0N20XKXdpZiouJZimWQmlpCYFj6VI4//zDt69eHdbhaW7WOjxSfIpqjXoRgXvugWuvTRxQIJTX1cH8+fltl0ixUVARycCiRTB9euo6dXWhnkg5U1ARycCbb8KQIanrVFeHeiLlTEFFJAP9+8OOHanr7NwZ6omUMwUVkQxceSU0Nqau09AQ6omUs7TLCYsIzJwZsr8mT06e/dXQELK/RMqZzlREMjBsGCxcGNKG584NKcZtbeF27txQvnCh0olF0gYVMxtsZk+b2SYz22hmNyaoM8XMXjKzDWbWZGbnxG2rNLNHzGyLmW02s/Oj8hPN7EkzezW6/XB2X5pIdk2cGM5EDh6E0aPDoMfRo8Pj5uawXaTcpR38aGanAKe4+zozOw5YC1zu7pvi6nwa2Ozu+8xsIvBddx8VbXsQeMbdG8zsaKCvu+83s38A/svdbzezm4EPu/ucVG3R4EcRka4rqsGP7r7b3ddF998BNgODOtVpcvd90cNmoArAzE4AxgKNUb333H1/VO8y4MHo/oPA5Uf2UkREpNC61KdiZkOBkcCzKapNB2KzIJ0K7AUeMLMXzKzBzPpF2wa6++7o/h5gYFfaIiIixSfjoGJmxwKPArPc/e0kdcYRgkrsMlYFcC7wM3cfCfwRuLnzfh6uwSW8DmdmM8xsjZmt2bt3b6bNFRGRAsgoqJhZb0JAecjdlySpczbQAFzm7m9Fxa1Aq7vHzmweIQQZgN9H/TWxfps3Eh3X3Re4e4271wwYMCCT5oqISIFkkv1lhD6Rze5+Z5I61cAS4Cp3fyVW7u57gN+a2fCo6EIg1sG/FPh6dP/rwC+79QpERKRoZDL4cTRwFbDBzNZHZbcA1QDufi9wG3AS8NMQg2iPyzS4Hngoyvx6DbgmKr8dWGxm04EdwBVH/nJERKSQtJ6KiEgPV1QpxZIdWoZWRMqBgkoeaBlaESkXuvyVY1qGVkQKTZe/ehAtQysi5URBJce0DK2IlBMFlRzTMrQiUk4UVHJMy9CKSDlRUMmxQixDq/RlESkUBZUcmzkT7rsvZHklEluGtr4+O8+n9GXp6fSlqbgpqORYPpehbWmBqVND+vK8eeGYFRXhdt68UD51qv74pHTpS1Px0ziVPGlpCWnDixaFTvn+/cMlr/r67I1PmT07/HHNm5e8zty54Q/xzoRTg4oUL4356r58jlNRUOlBBg4M39pS/UG1tIR11ffsyV+7RLJBX5q6T0ElCQWV1Hr1Cn9QFSnmnm5rC3+Y7e35a5dINuhLU/dpRL10i9KXpSfTmK/SoKDSgxQifVkkX/SlqTQoqPQg+U5fFsknfWkqDQoqPUg+05dF8k1fmkqDgkoPM3FiSKk8eDB0WPbpE24PHgzlEycWuoWFpYFzpUtfmkqDgkoPNGxYSKncsydkee3ZEx6X+x+bBs6VPn1pKn5KKZayoIFzUs6KKqXYzAab2dNmtsnMNprZjQnqTDGzl8xsg5k1mdk5cdu2R+XrzWxNXPkIM2uOlZvZedl7WSIdabE0kfxIe6ZiZqcAp7j7OjM7DlgLXO7um+LqfBrY7O77zGwi8F13HxVt2w7UuPubnY67AviRuy83s0nAt939s6naojMV6S4NnJNyls8zlRRjrwN33w3sju6/Y2abgUHAprg6TXG7NANVGTy3A8dH908Afpdhm0W6TAPnRPIjbVCJZ2ZDgZHAsymqTQfiuzwd+LWZ/Rn4P+6+ICqfBTxhZj8gXIb7dFfaItIVsYFzqc5UNHBO5MhlnP1lZscCjwKz3P3tJHXGEYLKnLjiC9x9BDARqDezsVH5N4HZ7j4YmA0kHNZkZjOiPpc1e/fuzbS5Ih1o4JxIfmSU/WVmvYF/B55w94Tzf5rZ2cBjwER3fyVJne8C/+3uPzCzPwCV7u5mZsAf3P34RPvFqE9FukvZX1LOii37ywhnEZtTBJRqYAlwVXxAMbN+Uec+ZtYPuAh4Odr8O+Az0f3xwKvdfREi6WjgnEh+ZNKnMhq4CthgZuujsluAagB3vxe4DTgJ+GmIQbRHUXEg8FhUVgEscvfHo2NcC9xlZhXAn4AZWXlFIknEBs7Nnx+yvOIXS9MZikh2aPCjiEgPV1SXv0RERDKloCIiIlmjoCIiIlmjoCIiIlmjoCIiIlnTpWlaRKS47PzDTpa9uozWt1s50HaAvr37UnV8FZNOn0T1CdWFbp6UIQUVkRLU3NrM4o2L2b5/O5NOn8SEYRPo27svB9oOsPWtrcx6fBZDK4dyxZlXUFtVW+jmShlRUBEpIe7O/S/cz1Pbn2L6yOmMHTKWiqM6/hmPGTKGq0dczaodq7j7ubvZ+MZGpo2cRjQIWSSnFFRESsj9L9zP8797nvmT5lP5ocqk9SqOqmD8qeM595RzmfPkHMyMaSOn5bGlUq7UUS9SIppbm3lq+1Pc/rnbqfxQJS0tMHt2WICsV69wO3t2mM8spvJDldzx+TtY+fpKmlubC9d4KRsKKiIlYvHGxUwfOZ3KD1WyfHmYdblPn7Ci5cGD4bZPn1C+PG5Fo8oPVTJtxDQe3vhw4RovZUOXv6TklGPG0479O9i+fztjh4ylpQWmTj18Gv9hw2DePJg8+fBp/McOGcv85+ez8w87e+x7JMVBZypSMppbm7npiZuY9fgsKo6qYMKwCUw5awoThk2g4qgKZj0+i5ueuKlHXuZZvm05k06fRMVRFdxzD1x7beJ1YSCU19WF2ZhjevfqzaTTJ7Hs1WX5abCULZ2pSNFTxhO0vt3KhGETAFi0KFzqSqWuLkzvf2fcCkjDTxrOipYVOWyliIKKlABlPHHoMh+EdWCGDEldv7o61IvXp3cfDrQdyFELRQIFFSlqsYyndAEFQtbTPffAokWV7H3nDh78Yj0rTjmD73+zNv0CXDt3wrJl0NoKBw5A375QVQWTJoX/0AUWG9gIYWGxHTtSLyq2c2eoF+/dtncPBSaRXFGfihS1+IynVDpnQ733diX3fXMar/d5+LBsqA6am+Gmm2DWLKiogAkTYMoUmDCBnfyBe793Gd+ZO4qb/ulrfOep73DvmnvZ+Yed2X+haVQdX8XWt7YCYaXKxsbU9RsaQr14W9/ayqDjB+WohSKBgkoBZDK+QDpmPKUSnw01b174Bl9RAV89fyyD/vJ1FvzLTqZO7fT+uof/zHffDZdeCosXh46IMWNoHtjGTQceY9aHn6VixnVMOO9rTFm1nwlb2qiwXgVJCIh1sre/387MmXDffbB6deK6q1eHoFJf/0FZ25/bWPbqMi45/ZL8NFjKloJKnnVlfEG5SBZkF67+IOMplWTZULGMp98fv+ywbCjuvx+efz4Ujh8PFRW4O43rGrn7ubu59OOXsvjLi6n71DcY88VZfPLOnzPmxf3UvdiLxV9ezKUfv5S7n7ubxnWN5GNJ7uoTqhlaOZRVO1YxbBgsXBjShufODe9fW1u4nTs3lC9c2PHy2Kodqzi18lQGnzA4522VMufuJfPzyU9+0kvZtm3u/fu7NzUl3t7UFLZv25bfdhXSsmXhNc+dG153W1u4nTvXvc8lt/o//MuqtMc4+eTk79mq7av8Oyu/49u2uQ8cGBWuXu1+5ZXu+/Z1qNuwtsG/8W/f8H3v7jv8QO6h/pVXhv3dfd+7+3zG0hneuK4x05d7RFb/drVf+eiVh9q3bZv77NnhdfXqFW5nzz78vdj37j7/2iNf89W/XZ2XdkrxAdZ4nv5Pp68Ag4GngU3ARuDGBHWmAC8BG4Am4Jy4bduj8vWdXxhwPbAlOu4/pGtLqQeVWbPCP8tUbr45/GMoB+mC7FcemO2Vf7EmbZA96qgQjBJ5ftfzftPjN/l774V/vO4e3uCVKzvU6/wPO76Ns2aFwHXUUe5fqvy1P3XuTYfatO/dfX7lo1fm7R92w9oGn7F0RvLA10m+A58Up3wGlUwuf7UD33L3M4BaoN7MzuhU53XgM+5+FvA9YEGn7ePcfYS718QKzGwccFkUgM4EfpBBW0raokUwfXrqOnV1oV45SDeI72PVfbnk8gMdL1slEMuGSiSW8XQoG2rHDti+HcZ27KdJlBCQ6FLlHavHcvJ/v86XPrWT5cvzPwXKtJHTGFU1ivpl9ax8bSVtf25LWK/tz22sfG0l9cvqGVU1imtGXJOX9omkTSl2993A7uj+O2a2GRhEOHOJ1YkfitUMVGXw3N8Ebnf3g9Ex3uhCu0tSd8cX9FTpBvFVHV9F37Fb+cnVYzoM4usslg01b97h22IZT4eyoZYvD2nCFR989BMlBCSdCuUTveFvJrHklWWMmnodzc0wdmj+pkCJjb05Y8AZPLzxYeY/P59Jp09i+EnD6dO7D++2vcuWN7ewbNsyTqs8jevPu17rqUhedamj3syGAiOBZ1NUmw7Edzc78GszW2tmM+LKPw6MMbNnzew/zOxTXWlLKUr1jTom0fiCnipdkJ10+iSe27eMvW+1pzxOsmyoWMbTR9655INsqNZWGD68Q734KVBiUp5FDR/OacfsOtT5X4gpUGqravnhhB/y44t/TPv77axoWcHPN/ycFS0r+LP/mbsuvosfTvihAorkXcaDH83sWOBRYJa7v52kzjhCULkgrvgCd99lZicDT5rZFndfFT33iYRLap8CFpvZadH1v/hjzgBmAFQXwSC0I5HqG3VMovEFPUXniSA/NKkvf/9kFV//dOKJIKtPqObEo4ZywtmrgPFJjxufDVVXF36qq+EXz66i9eVTufbewR9kQ8UGNsaJnwIlJuVZVJ8+cOAAdf/zg6lQCjUFSvUJ1VxXc13en1ckmYyCipn1JgSUh9x9SZI6ZwMNwER3fytW7u67ots3zOwx4DxgFdAKLImCyHNm9j7QH9gbf1x3X0DUR1NTU5P73M0cmjkzXKOfPDnxN+DY+ILmHjYfYrKlb1uHH2D1f27lhbeTL33b/uIVDLz4bvb/6dyUAyAnTgzv2/z54R/93nf2c/QXG7n8lBv4edxsvfTtGwJLnPgpUGJSnkW9+y707dvhUqWmQBEJ0gYVCzPyNQKb3T3hlW0zqwaWAFe5+ytx5f2Ao6K+mH7ARcD/jjb/KzAOeNrMPg4cDfTo3oRk36h37gzBpKHh8PEFpczTTAT599+E2toxLPnXq2n76OETQa5eDcvvq+Wm/7uROU/O4Y7P35EysAwbFs4abpu3nzlPzmFU1eeYNrLT5Z+qKti6FcaMOVQUPwVKTMqpULZuhUGDOlyq1BQoIkEmfSqjgauA8Wa2PvqZZGbXmVnsvPs24CTgp9H2NVH5QOA3ZvYi8BzwK3d/PNp2P3Camb0M/AL4eudLXz1R7Bv1wYPhG3WfPuH24MFQPnFioVuYPfETQY4/dfxhgxhjQfavLq/gyQXj+eu/mE/zb5/jjhUPdBjEd/NFWcx4mjQpzPHV/kE/TfwUKDFJp0Jpawv7X3JJh0uVmgJFJLBS+j9eU1Pja9asSV9RCq65tZm7n7s744kg588P/Rh739nP0ZfXc/lHr+fvrus4EWRzazMPb3yY1/e/njLj6ctnfjl1B/VNN4WpWcaHfpqdf9jJrMdnsfjLiw8FvpaWcKmyc/YXK1fCsmWs/h8/PLQQVvXQNr7yyFe46+K7NGJdipKZrY0f0pFLmqVYciLTiSA/mFk49E/071/JBR+dRt9PPcywYR0DQ21VLbVVtYc6/Fe0rDjUHzLo+EHcdfFdmaX0XnFFmPPr3HOhsrLDFCjjTw2BJuGlyuP386e7Gvnnyhv427ipUFa+pilQRGIUVCTrMp0IcvnyMBbk2mtDptWQIaEfY0HjWH787/MZe/ROvn554qywI8p4qq2FjRthzhy44w6orOSKM6/g7ufu5txTPkgIiO/8n3j+fv76zTlsPu5z2PTaQ0v17v/TfhpfaOSGUTd0vz0iPYguf0nW3bvmXiqOqqDu3LqkdZJeXorMfbiBu+e382LjdblJXHCHBx4Il7OmTYOxY2ncsJDndj3XMSGgrQ1WrQoTUF54IVxzDUSrSe7/UywhYFSPWQxMeqZ8Xv7SLMWSda1vtzL8pOEp66SbomXSecM5Z8yutFO0dJtZCCbXXx863r/yFaatfZ9RBz5M/T9/lZW//DFtC34WLpUtWxbqTZsGZpoCJUu0BETPpMtfknWJxn10lm6Klj69+/CJsw6w6AZSTtFyxGprw8/OndiyZUx7uTdnHPwwD7f+gvl9DjKp/qsMP/3T9Oldwbs7ntEUKFmS7NJnY2P4dSxc2LMyIcuJgopkXaJxH52lm6Ll3bZ3GXhi3/zNg1ZdDdeFfpra6OeIEwIkoaTzqg0Ls01MnsyhzLqeMmarnCioSNbFxn2MGTImaZ1066xvfWsrfdoGFXQeNE2BkhvpLn2efz6H5lXL6Vmq5IT6VCTr4pe+TSbVOuuxiSD3rLqkx86DVs60BETPpqAiWRc/7iOZVOusr9qximPePZXFDYM7rLMu2VWojnItAdGzKahITlxx5hU0vtDI/j/tT7g92TrrL2zezw3/1MjyH3y5R82DVmwSLUDW1BQe19aG7bmiJSB6NgUVyYnaqlrGDx3PnCfnJA0snedB+1Dlfkb/7RwGHfwca39Zq+yfHInvKJ83LwTuiooPOsqXLg3bc3XGkurSZ0xPXgKip9PgR8kZd+eB9Q+w8vWVTBsxjbFDxtK7V+/D6rX9uY1VO1Zx//r7ufDUC7lmxDVYNMBQsm/27HBGkmpdn7lzQ7DPRUd5uoGvq1cr+yvb8jn4UUFFci5rE0FKVgwcGC51pfqH3dISzh737MlNG2LjVFItAaEz1exRUElCQaW0xcZ97Hp7V4dxH5NOT7zyo+RGr17hLKQixYCCtrZwNtOeeiXnIxI/O3WYTDRc8qqv1xlKtimoJKGgInLkiuFMRfJLc3+JSM6oo1xySSPqRcrMzJmho3zy5OQd5Q0NoaNcpKt0piJSZpKNEWppocMyzurXkO5QUBEpQ53HCPXpE24PHgzlyryS7lJQESlTw4aFcSh79oQsrz17wuNCnaFofZWeIW1QMbPBZva0mW0ys41mdmOCOlPM7CUz22BmTWZ2Tty27VH5ejM7LHXLzL5lZm5mmpRBpEwVctoYya5MOurbgW+5+zozOw5Ya2ZPuvumuDqvA59x931mNhFYAIyK2z7O3Q+bHs7MBgMXATu7/xJEpJRpfZWeJe2Zirvvdvd10f13gM3AoE51mtx9X/SwGajK8Pl/BHwbKJ3BMiKSVV1ZX0WKX5f6VMxsKDASeDZFtelA/MmqA782s7VmNiPuWJcBu9z9xa60QUR6Fq2v0rNkHFTM7FjgUWCWu7+dpM44QlCZE1d8gbuPACYC9WY21sz6ArcAt2XwvDPMbI2Zrdm7d2+mzRWREpGP9VWUBJA/GQUVM+tNCCgPufuSJHXOBhqAy9z9rVi5u++Kbt8AHgPOA4YBpwIvmtl2wuWydWb2kc7HdfcF7l7j7jUDBgzoymsTkRKQ6/VVlASQX5lkfxnQCGx294QTYZtZNbAEuMrdX4kr7xd17mNm/Qid8i+7+wZ3P9ndh7r7UKAVONfdNdOQSJnJ5bQxhV47phylnVDSzC4AngE2AO9HxbcA1QDufq+ZNQBfAmLfN9rdvcbMTiOcnUDINFvk7t9P8BzbgZpEGWLxeveu8fffX3NoNtOZM5UNIlLqcrm+SqHXjikWmqU4ibPOqvEXXljDjh3hm81992ndBZGeIFfrq2hG5kCzFCdxzDE6dRXpiXI1bUw+kgCko5I6U0m0nko5nLqKSPfoTCXQmUoXKH9dJP9KJUVXa8fkX8kHFZ26iuRXKaXozpwZ+l5Xr068PbZ2TH19ftvVk5V8UDmS/HUR6ZpSS9HV2jH5V/JBRaeuIvlTivN0ae2Y/CrpjvojyV8XkQ+0tISAsWhRuJycbCyYOr5Lkzrqkzh4UKeuItnWlT4SpehKOiV1ptK7d427fzCivr5eAUXkSHR1NLvOVEqTzlSSOOec4lj2VAqrVNJZS0FX+0iUoivplFRQESmldNZS0NW1TJSiK+mU1OWvRCPqpXzkcuLBctWrVwjMFSkWFm9rC0G7vT08ztU8XZI7uvwlkkApprMWu+6sZaIUXUlFZypSMtRJnHnqb6Y0NXx50JmKSALlns6ai/4k9ZFItimoSMnI9bKzxSxX06NoGhPJNgUVKRnlnM6ay/4k9ZFINqlS048WAAAMLklEQVRPRUpGOWd/qT9JjoT6VEQSKOdLNeXenySlQ0FFSkq5Xqop5/4kKS1pg4qZDTazp81sk5ltNLMbE9SZYmYvmdkGM2sys3Pitm2Pyteb2Zq48n80sy3Rfo+ZWWX2Xpb0ZMOGhfTWPXvKZ9qecu5PktKSyZlKO/Atdz8DqAXqzeyMTnVeBz7j7mcB3wMWdNo+zt1HdLqm9yTwl+5+NvAKMLdbr0CkDCj1V0pF2qDi7rvdfV10/x1gMzCoU50md98XPWwGqjI47gp3b+/KPiLlqpz7k6S0dKlPxcyGAiOBZ1NUmw7ED8Ny4NdmttbMZiTZZ1qnfUS6rKfPXlyu/UlSWjIOKmZ2LPAoMMvd305SZxwhqMyJK77A3UcAEwmXzsZ22udWwiW2h5Icc4aZrTGzNXv37s20uVJmymX24nLsT5LSktE4FTPrDfw78IS7J5wByMzOBh4DJrr7K0nqfBf4b3f/QfT4auAbwIXufiBdOzRORRIp5/ErIpkoqnEqZmZAI7A5RUCpBpYAV8UHFDPrZ2bHxe4DFwEvR48vBr4NfCGTgCKSjGYvFikeac9UzOwC4BlgA/B+VHwLUA3g7veaWQPwJSCWSd/u7jVmdhrh7AWgAljk7t+PjrsNOAZ4K9re7O7XpWqLzlQkEY02F0ktn2cqmqZFSl53FpoSKSdFdflLpNhptLlI8VBQkZKn0eYixSPFBQOR0jBzZsj+mjw5efZXQ0PI/hKR3NKZipQ8jTYXKR4KKtIjaLS5SHFQ9peISA+n7C8RESlJCioiIpI1CioiIpI1CioiIpI1CioiIpI1CioiIpI1CioiIpI1CioiIpI1CioiIpI1CioiIpI1CioiIpI1CioF0NICs2eHZXB79Qq3s2eHchGRUqagkmfLl4e1P/r0CeuqHzwYbvv0CeXLlxe6hSIi3adZivOopSUEjqVLky8m9YUvhKnatfaHiGRLUc1SbGaDzexpM9tkZhvN7MYEdaaY2UtmtsHMmszsnLht26Py9Wa2Jq78RDN70sxejW4/nL2XVZzuuQeuvTZxQIFQXlcH8+fnt10iItmSyeWvduBb7n4GUAvUm9kZneq8DnzG3c8Cvgcs6LR9nLuP6BQpbwZWuvvpwMrocY+2aBFMn566Tl1dqCciUorSBhV33+3u66L77wCbgUGd6jS5+77oYTNQlcFzXwY8GN1/ELg800aXqjffhCFDUteprg71RERKUZc66s1sKDASeDZFtelAfHezA782s7VmNiOufKC7747u7wEGdqUtpah/f9ixI3WdnTtDPRGRUlSRaUUzOxZ4FJjl7m8nqTOOEFQuiCu+wN13mdnJwJNmtsXdV8Xv5+5uZgkzBqJAFAtGB83s5UzbXED9gQTnG0MGf+xj/j7s3JV81+pBYGa2ozVnrftAknYWlVJoI6id2aZ2ZtfwfD1RRkHFzHoTAspD7r4kSZ2zgQZgoru/FSt3913R7Rtm9hhwHrAK+L2ZneLuu83sFOCNRMd19wVEfTRmtiZfGQxHQu3MnlJoI6id2aZ2Zld8klSuZZL9ZUAjsNnd70xSpxpYAlzl7q/Elfczs+Ni94GLgNiZxlLg69H9rwO/7O6LEBGR4pDJmcpo4Cpgg5mtj8puAaoB3P1e4DbgJOCnIQbRHkXvgcBjUVkFsMjdH4+OcTuw2MymAzuAK7LyikREpGDSBhV3/w1gaerUAXUJyl8Dzjl8D4gukV2YWTMP6ZyqXKzUzuwphTaC2pltamd25a2dJTWiXkREipvm/hIRkexx95z+ABcDW4FtwM0Jthvwk2j7S8C56fYF/hHYEtV/DKiM2zY3qr8VmBBX/klgQ7TtJ0RnaYVoJ/B5YG3UnrXA+Lh9/l90rPXRz8kFbOdQ4N24ttybyfuZ5zZOiWvfeuB9YEQB38vvRXXXAyuAjxbpZzNhOym+z2aydg6lG5/NArSzqD6fcdu/RRhD2P9IP5+HHTvVxiP9AXoBLcBpwNHAi8AZnepMIgyWNMI0MM+m25eQRVYR3b8DuCO6f0ZU7xjg1Gj/XtG256LjW/R8EwvYzpFxH7q/BHZ1+sOtKZL3cyjwcpK2JHw/893GTsc9C2gp8Ht5fNz+NxD9s6P4PpvJ2llsn81k7RxKFz+bhWhnsX0+o+2DgScICVL9j+Tzmegn15e/zgO2uftr7v4e8AvC9CzxLgMWetAMVEbjVpLu6+4r3L092j9+WpjLgF+4+0F3f50QWc+Ljne8uzd7eJcW0nFamLy2091fcPffReUbgT5mdkwRvp8JpXk/C9nGr0X7ZCJX7YwfGNyP8G0wdqxi+mwmbGcRfjaTvZ8JFdv72UnBP5+RHwHf7tTG7n4+D5ProDII+G3c41Y6zRuWok4m+wJM44NpYVIdqzVBeaHaGe9LwDp3PxhX9mA0q/P/isYJFbKdp0Zt+Q8zGxP3HMnez0K+l18Bft6pLO/vpZl938x+S7j0cVsGxyrIZzNJO+MVxWczRTu7+tksVDtjCv75NLPLCGefL3bhWKnez8OUdEe9md1KmEX5oUK3JZVk7TSzMwmXcr4RVzzF3c8ExkQ/VxWwnbuBancfAdwELDKz4/PVnkRSvJejgAPuHj+NT0HeS3e/1d0HR22cmY/n7I5U7Symz2aSdhbdZzPN+1nwz6eZ9SWMMUwU8LIm10FlF+H6XUxVVJZJnZT7mtnVwKWEX0zsNC7VsaoSlBeqnZhZFaHDeaq7H1pI2D+Y1uYdYBHhVLYg7YxOhd+K7q8lXGf9OKnfz7y/l5Gv0ulbYKHeyzgPEb7tpztWQT6bSdpZdJ/NRO3s5mcz7+2MUwyfz2GE/pIXzWx7VL7OzD6S5lip3s/DeYoOlyP9IQyufC16IbEOozM71bmEjp1Nz6Xbl5DZsAkY0OlYZ9Kxs+k1knc2TSpgOyujen+V4P2KdZz1Bh4BritgOwfEvX+nET5MJ6Z6P/PdxmjbUVHbTiuC9/L0uP2vBx4p0s9msnYW22czWTu7/NksRDuL7fPZaf/tcW3o1ucz4f/9TANEd38IGQqvEL5J3BqVXRd7A6OGzo+2byAuGyLRvlH5NsL1v0TphLdG9bfSMeujhjDvWAtwD4enGeatncB3gD/SMdXwZEIH31pCeuBG4K7YL7ZA7fxS1I71wDpgcibvZwF+558Fmju9T4V6Lx+N3peXgH8DBhXpZzNhOym+z2aydnbrs1mg3/tnKZLPZ6fjb6djSnG3Pp+dfzSiXkREsqakO+pFRKS4KKiIiEjWKKiIiEjWKKiIiEjWKKiIiKRhZl82s41m9r6ZJV0+2My2m9mGaJT8mrjyE83sSTN7Nbr9cFT+eTNbG+2z1szGR+V9zexXZrYlet7bM2jjFDN7KTpWk5klXMsq1xRURETimNlnzeyfOhW/DPwVsCqDQ4xz9xHece36m4GV7n46sDJ6DPAmIR36LMKy6v8ct88P3P0ThEk+R5vZxDTP+zrwmehY36NAC4gpqIiIpOHum9196xEc4jLgwej+g0STMnqSCTzd/YC7Px3VeY8wFqcKwMwGmNmjZvZ89DM6qtfk7vuiY6WdGDZXFFRERLLHgV9Hl7JmxJUPdPfd0f09wMAE+yaawBMzqwQmE85wIAyU/JG7fyrapyHBsaaTeNLVnEu7Rr2ISDkws2cJ05QcC5xoZuujTXPc/YkMD3OBu+8ys5OBJ81si7t3uGTm7m5mHUadx03geVGn8grCnGE/cffXouLPAWfETWp8vJkd6+7/He0zjhBULsiwzVmloCIiArj7KAh9KsDV7n51N44RmyTyDTN7jDBJ5Crg92Z2irvvjtYoeSO2T7IJPCMLgFfd/cdxZUcBte7+p87Pb2ZnE85cJno04Wa+6fKXiEgWmFk/Mzsudp9w1hGb6n4poSOe6PaXUb1K4FeEJX//s9Px/g44AZjV6alWECatjNUbEd1WA0uAq9z9ley9sq5RUBERScPMvmhmrcD5wK/M7Imo/KNmtiyqNhD4jZm9SJjZ91fu/ni07Xbg82b2KuHyVSxFeCbwMeC2KA15vZmdHJ293EpY5nddVF4X7XMDUBOlD28iTDIJYZ2Uk4Cfdk5pzidNKCkiIlmjMxUREckaBRUREckaBRUREckaBRUREckaBRUREckaBRUREckaBRUREckaBRUREcma/w+STDFY1iTTpgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5e94a954e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Zoom in, comparing detected sources in COSMOS (blue) to the parent (red) and children (green):\n",
    "plt.figure()\n",
    "#plt.scatter(df_cosmos['ra'], df_cosmos['dec'], marker='.', s=1)\n",
    "plt.plot(df_cosmos['ra'], df_cosmos['dec'], 'bo',markersize=10,fillstyle='none')\n",
    "#plt.hist2d(df_cosmos['ra'], df_cosmos['dec'], bins=(300, 300), cmap='Blues')\n",
    "plt.plot(pra[pick_parent_ind],pdec[pick_parent_ind],'r.',markersize=35,fillstyle='none',alpha=0.75)\n",
    "plt.plot(cra,cdec,'g.',markersize=35,fillstyle='none',alpha=0.75)\n",
    "#plt.colorbar(label='Binned Counts')\n",
    "plt.xlim(150.222,150.224)\n",
    "plt.ylim(2.25,2.265)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "cosmos_ra = np.array(df_cosmos['ra'])\n",
    "cosmos_dec = np.array(df_cosmos['dec'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now, we would like to mark the detected points from HSC and from COSMOS on top of the image.\n",
    "# This should work, in theory, but we ran out of time to finish this. \"afw_display.dot\" expects X,Y coordinates,\n",
    "#   so you need to convert from RA, Dec. \n",
    "# Note also that the \"Buffering\" portion here would cycle through all the objects before sending the XPA command\n",
    "#   to Firefly.\n",
    "\n",
    "disp1.mtv(coadd)\n",
    "# Note also that this displays the entire image -- ideally you would like a postage stamp around the object of interest.\n",
    "\n",
    "with afw_display.Buffering():\n",
    "#    for ind in range(0,2): #len(np.array(df_cosmos['ra']))):\n",
    "#        afw_display.dot('o', df_cosmos['ra'][ind], df_cosmos['dec'][ind], size=20, ctype='orange')\n",
    "        afw_display.dot('o', 1000.0, 1000.0, size=200, ctype='orange') # plot a circle at (X,Y)=(1000,1000)\n",
    "#        afw_display.dot('o', record.getX(), record.getY(), size=20, ctype='orange')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "150.22291945562662 2.2559620121522577\n"
     ]
    }
   ],
   "source": [
    "print(pra[pick_parent_ind],pdec[pick_parent_ind])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}