{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# XMM-LSS master catalogue\n", "\n", "This notebook presents the merge of the various pristine catalogues to produce the HELP master catalogue on XMM-LSS." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This notebook was run with herschelhelp_internal version: \n", "017bb1e (Mon Jun 18 14:58:59 2018 +0100)\n", "This notebook was executed on: \n", "2019-02-07 12:43:54.870832\n" ] } ], "source": [ "from herschelhelp_internal import git_version\n", "print(\"This notebook was run with herschelhelp_internal version: \\n{}\".format(git_version()))\n", "import datetime\n", "print(\"This notebook was executed on: \\n{}\".format(datetime.datetime.now()))" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/pyenv/versions/3.7.2/lib/python3.7/site-packages/matplotlib/__init__.py:855: MatplotlibDeprecationWarning: \n", "examples.directory is deprecated; in the future, examples will be found relative to the 'datapath' directory.\n", " \"found relative to the 'datapath' directory.\".format(key))\n", "/opt/pyenv/versions/3.7.2/lib/python3.7/site-packages/matplotlib/__init__.py:846: MatplotlibDeprecationWarning: \n", "The text.latex.unicode rcparam was deprecated in Matplotlib 2.2 and will be removed in 3.1.\n", " \"2.2\", name=key, obj_type=\"rcparam\", addendum=addendum)\n", "/opt/pyenv/versions/3.7.2/lib/python3.7/site-packages/seaborn/apionly.py:9: UserWarning: As seaborn no longer sets a default style on import, the seaborn.apionly module is deprecated. It will be removed in a future version.\n", " warnings.warn(msg, UserWarning)\n" ] } ], "source": [ "%matplotlib inline\n", "#%config InlineBackend.figure_format = 'svg'\n", "\n", "import matplotlib.pyplot as plt\n", "plt.rc('figure', figsize=(10, 6))\n", "\n", "import os\n", "import time\n", "import gc\n", "from astropy import units as u\n", "from astropy.coordinates import SkyCoord\n", "from astropy.table import Column, Table, join\n", "import numpy as np\n", "from pymoc import MOC\n", "\n", "from herschelhelp_internal.masterlist import merge_catalogues, nb_merge_dist_plot, specz_merge\n", "from herschelhelp_internal.utils import coords_to_hpidx, ebv, gen_help_id, inMoc" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "TMP_DIR = os.environ.get('TMP_DIR', \"./data_tmp\")\n", "OUT_DIR = os.environ.get('OUT_DIR', \"./data\")\n", "SUFFIX = os.environ.get('SUFFIX', time.strftime(\"_%Y%m%d\"))\n", "\n", "try:\n", " os.makedirs(OUT_DIR)\n", "except FileExistsError:\n", " pass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## I - Reading the prepared pristine catalogues" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# SXDS now merged in with CANDELS in CFHT merge notebook\n", "#sxds = Table.read(\"{}/SXDS.fits\".format(TMP_DIR)\n", "# )[\"sxds_intid\", \n", "# \"sxds_ra\", \n", "# \"sxds_dec\",\n", "# # 'sxds_flag_merged'\n", "# ] # 1.11\n", "\n", "#WIRCAM:\n", "#cfht_wirds = Table.read(\"{}/CFHT-WIRDS.fits\".format(TMP_DIR)\n", "# )[\"wirds_intid\",\n", "# \"wirds_ra\",\n", "# \"wirds_dec\"] # 1.3\n", "\n", "#Megacam:\n", "#candels = Table.read(\"{}/CANDELS.fits\".format(TMP_DIR)\n", "# )[\"candels_id\",\n", "# \"candels_ra\",\n", "# \"candels_dec\"] # 1.1\n", "#cfhtls_wide = Table.read(\"{}/CFHTLS-WIDE.fits\".format(TMP_DIR)\n", "# )[\"cfhtls-wide_id\",\n", "# \"cfhtls-wide_ra\",\n", "# \"cfhtls-wide_dec\",\n", "# \"cfhtls-wide_stellarity\"] # 1.4a\n", "#cfhtls_deep = Table.read(\"{}/CFHTLS-DEEP.fits\".format(TMP_DIR)\n", "# )[\"cfhtls-deep_id\",\n", "# \"cfhtls-deep_ra\",\n", "# \"cfhtls-deep_dec\",\n", "# \"cfhtls-deep_stellarity\"] # 1.4b\n", "#We no longer use CFHTLenS as it is the same raw data set as CFHTLS-WIDE\n", "# cfhtlens = Table.read(\"{}/CFHTLENS.fits\".format(TMP_DIR)) # 1.5\n", "#sparcs = Table.read(\"{}/SpARCS.fits\".format(TMP_DIR)\n", "# )['sparcs_intid', \n", "# 'sparcs_ra', \n", "# 'sparcs_dec', \n", "# 'sparcs_stellarity'] # 1.12\n", "#\n", "#vipers = Table.read(\"{}/VIPERS.fits\".format(TMP_DIR)\n", "# )[\"vipers_id\",\n", "# \"vipers_ra\",\n", "# \"vipers_dec\"] # 1.15\n", "cfht = Table.read(\"{}/cfht_merged_catalogue_xmm-lss.fits\".format(TMP_DIR)\n", " )['cfht_intid', \n", " \"candels_id\",\n", " \"cfhtls-wide_id\",\n", " \"cfhtls-deep_id\",\n", " 'sparcs_intid',\n", " 'wirds_id',\n", " 'vipers_id',\n", " 'sxds_intid',\n", " 'sxds_b_id',\n", " 'sxds_v_id',\n", " 'sxds_r_id',\n", " 'sxds_i_id',\n", " 'sxds_z_id',\n", " 'cfht_ra', \n", " 'cfht_dec'] # 1.12\n", "\n", "\n", "\n", "#DECam:\n", "#decals = Table.read(\"{}/DECaLS.fits\".format(TMP_DIR)\n", "# )[\"decals_id\",\n", "# \"decals_ra\",\n", "# \"decals_dec\"] # 1.6\n", "decam = Table.read(\"{}/decam_merged_catalogue_xmm-lss.fits\".format(TMP_DIR)\n", " )[\"decam_intid\",\n", " \"decals_id\",\n", " \"des_id\",\n", " \"decam_ra\",\n", " \"decam_dec\"] # 1.6\n", "\n", "#Spitzer IRAC:\n", "#servs = Table.read(\"{}/SERVS.fits\".format(TMP_DIR)\n", "# )[\"servs_intid\",\n", "# \"servs_ra\",\n", "# \"servs_dec\"] # 1.8\n", "#swire = Table.read(\"{}/SWIRE.fits\".format(TMP_DIR)\n", "# )[\"swire_intid\",\n", "# \"swire_ra\",\n", "# \"swire_dec\"] # 1.7\n", "irac = Table.read(\"{}/irac_merged_catalogue_xmm-lss.fits\".format(TMP_DIR)\n", " )[\"irac_intid\",\n", " \"servs_intid\",\n", " \"swire_intid\",\n", " \"irac_ra\",\n", " \"irac_dec\"] # 1.7\n", "\n", "\n", "#Hyper Suprime Cam:\n", "#hsc_wide = Table.read(\"{}/HSC-WIDE.fits\".format(TMP_DIR)\n", "# )[\"hsc-wide_id\",\n", "# \"hsc-wide_ra\",\n", "# \"hsc-wide_dec\", \n", "# \"hsc-wide_stellarity\"] # 1.9a\n", "#hsc_deep = Table.read(\"{}/HSC-DEEP.fits\".format(TMP_DIR)\n", "# )[\"hsc-deep_id\",\n", "# \"hsc-deep_ra\",\n", "# \"hsc-deep_dec\", \n", "# \"hsc-deep_stellarity\"] # 1.9b\n", "#hsc_udeep = Table.read(\"{}/HSC-UDEEP.fits\".format(TMP_DIR)\n", "# )[\"hsc-udeep_id\",\n", "# \"hsc-udeep_ra\",\n", "# \"hsc-udeep_dec\", \n", "# \"hsc-udeep_stellarity\"] # 1.9c\n", "hsc = Table.read(\"{}/hsc_merged_catalogue_xmm-lss.fits\".format(TMP_DIR)\n", " )[\"hsc_intid\",\n", " \"hsc_ra\",\n", " \"hsc_dec\", \n", " \"hsc-wide_id\",\n", " \"hsc-deep_id\",\n", " \"hsc-udeep_id\"] # 1.9c\n", "\n", "\n", "#GPC1:\n", "ps1 = Table.read(\"{}/PS1.fits\".format(TMP_DIR)\n", " )[\"ps1_id\",\n", " \"ps1_ra\",\n", " \"ps1_dec\"] # 1.10\n", "\n", "# UKIDSS WFCAM:\n", "#dxs = Table.read(\"{}/UKIDSS-DXS.fits\".format(TMP_DIR)\n", "# )['dxs_id',\n", "# 'dxs_ra',\n", "# 'dxs_dec', \n", "# 'dxs_stellarity'] # 1.13\n", "#uds = Table.read(\"{}/UKIDSS-UDS.fits\".format(TMP_DIR)\n", "# )['uds_id',\n", "# 'uds_ra',\n", "# 'uds_dec',\n", "# 'uds_stellarity'] # 1.14\n", "ukidss = Table.read(\"{}/ukidss_merged_catalogue_xmm-lss.fits\".format(TMP_DIR)\n", " )['ukidss_intid',\n", " 'dxs_id',\n", " 'uds_id',\n", " 'ukidss_ra',\n", " 'ukidss_dec',\n", " # 'ukidss_stellarity'\n", " ] # 1.14\n", "\n", "\n", "\n", "#VIRCAM:\n", "#vhs = Table.read(\"{}/VISTA-VHS.fits\".format(TMP_DIR)\n", "# )[\"vhs_id\",\n", "# \"vhs_ra\",\n", "# \"vhs_dec\",\n", "# \"vhs_stellarity\"] # 1.16\n", "#video = Table.read(\"{}/VISTA-VIDEO.fits\".format(TMP_DIR)\n", "# )['video_id',\n", "# 'video_ra',\n", "# 'video_dec',\n", "# 'video_stellarity',\n", "# 'video_flag_gaia'] # 1.17\n", "#viking = Table.read(\"{}/VISTA-VIKING.fits\".format(TMP_DIR)\n", "# )[\"viking_id\",\n", "# \"viking_ra\",\n", "# \"viking_dec\",\n", "# \"viking_stellarity\",\n", "# \"viking_flag_gaia\"] # 1.18\n", "vircam = Table.read(\"{}/vista_merged_catalogue_xmm-lss.fits\".format(TMP_DIR)\n", " )[\"vircam_intid\",\n", " \"vhs_id\",\n", " 'video_id',\n", " \"viking_id\",\n", " \"vircam_ra\",\n", " \"vircam_dec\",\n", " #\"vircam_stellarity\",\n", " #\"vircam_flag_gaia\"\n", " ] " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## II - Merging tables\n", "\n", "We first merge the optical catalogues and then add the infrared ones. We start with PanSTARRS because it coevrs the whole field.\n", "\n", "At every step, we look at the distribution of the distances separating the sources from one catalogue to the other (within a maximum radius) to determine the best cross-matching radius." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add PanSTARRS" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "master_catalogue = ps1\n", "master_catalogue['ps1_ra'].name = 'ra'\n", "master_catalogue['ps1_dec'].name = 'dec'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### CANDELS\n", "\n", "We now use CANDELS-UDS which must be individually merged with the merged catalogues since it has measurements from different instruments" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#nb_merge_dist_plot(\n", "# SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n", "# SkyCoord(candels['candels_ra'], candels['candels_dec'])\n", "#)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Given the graph above, we use 0.8 arc-second radius\n", "#master_catalogue = merge_catalogues(master_catalogue, candels, \"candels_ra\", \"candels_dec\", radius=0.8*u.arcsec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add CFHT\n", "\n", "We independently merge all the CFHT Megacam and WIRCAM and some CANDELS in CFHT_Merge notebook " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEKCAYAAAAVaT4rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuYXFWd7vHvr2/V90u6q3O/kxBQIEAD4YgSR8HoMDCOjIByGUcn4zygo+PMPHp0xMF5zsHjc3RUdDCDOegMl2HAS1QUcUDjaAJJALnknoYk3emkO+lL+n79nT9qdyg6fanuru6qzn4/z1NPqtbeu2rtpPPW6rXXXsvcHRERCY+MVFdARESml4JfRCRkFPwiIiGj4BcRCRkFv4hIyCj4RURCRsEvIhIyCn4RkZBR8IuIhExWqiswnIqKCl+yZEmqqyEiMmPs2LHjuLtHE9k3LYN/yZIlbN++PdXVEBGZMczsYKL7qqtHRCRkFPwiIiGj4BcRCRkFv4hIyCj4RURCRsEvIhIyCn4RkZBR8IuIhIyCX0QkZNLyzt108+Azh04r+8Bli1JQExGRyVOLX0QkZBT8IiIho+AXEQkZBb+ISMiMGfxmttDMnjaznWb2ipn99TD7mJl93cz2m9mLZnZR3LbbzGxf8Lgt2ScgIiLjk8ionj7gU+7+nJkVATvM7El33xm3z7uBFcHjMuBfgMvMbBZwJ1AFeHDsJndvSupZiIhIwsZs8bt7nbs/FzxvBXYB84fsdh3wPY/ZCpSa2VzgXcCT7t4YhP2TwLqknoGIiIzLuPr4zWwJcCHwzJBN84HDca9rgrKRykVEJEUSDn4zKwQeAz7h7ieTXREzW29m281se0NDQ7LfXkREAgkFv5llEwv9B9z9+8PsUgssjHu9ICgbqfw07r7B3avcvSoaTWi9YBERmYBERvUY8B1gl7t/ZYTdNgG3BqN71gAt7l4HPAFcbWZlZlYGXB2UiYhIiiQyquctwC3AS2b2QlD2P4FFAO5+L/A48B5gP9ABfCjY1mhmXwS2Bcfd5e6Nyau+iIiM15jB7+7/DdgY+zhw+wjbNgIbJ1Q7ERFJOt25KyISMgp+EZGQUfCLiISMgl9EJGQU/CIiIaPgFxEJGQW/iEjIKPhFREJGwS8iEjIKfhGRkFHwi4iEjIJfRCRkFPwiIiGj4BcRCRkFv4hIyCj4RURCZsyFWMxsI3ANUO/ubx5m+98BH4x7v3OAaLD61mtAK9AP9Ll7VbIqLiIiE5NIi/9+YN1IG939y+6+2t1XA58Bfj1kecW3B9sV+iIiaWDM4Hf3zUCi6+TeBDw0qRqJiMiUSlofv5nlE/vN4LG4Ygd+YWY7zGx9sj5LREQmbsw+/nH4I+C3Q7p5rnD3WjOrBJ40s93BbxCnCb4Y1gMsWrQoidUSEZF4yRzVcyNDunncvTb4sx74AXDpSAe7+wZ3r3L3qmg0msRqiYhIvKQEv5mVAFcCP4orKzCzosHnwNXAy8n4PBERmbhEhnM+BKwFKsysBrgTyAZw93uD3d4L/MLd2+MOnQ38wMwGP+dBd/958qouIiITMWbwu/tNCexzP7Fhn/Fl1cAFE62YiIhMDd25KyISMgp+EZGQUfCLiISMgl9EJGQU/CIiIaPgFxEJGQW/iEjIKPhFREJGwS8iEjIKfhGRkFHwi4iEjIJfRCRkFPwiIiGj4BcRCRkFv4hIyCj4RURCZszgN7ONZlZvZsMum2hma82sxcxeCB6fj9u2zsz2mNl+M/t0MisuIiITk0iL/35g3Rj7/MbdVwePuwDMLBP4JvBu4FzgJjM7dzKVFRGRyRsz+N19M9A4gfe+FNjv7tXu3gM8DFw3gfcREZEkSlYf/+Vm9nsz+5mZvSkomw8cjtunJigblpmtN7PtZra9oaEhSdUSEZGhkhH8zwGL3f0C4BvADyfyJu6+wd2r3L0qGo0moVoiIjKcSQe/u59097bg+eNAtplVALXAwrhdFwRlIiKSQpMOfjObY2YWPL80eM8TwDZghZktNbMc4EZg02Q/T0REJidrrB3M7CFgLVBhZjXAnUA2gLvfC1wP/JWZ9QGdwI3u7kCfmd0BPAFkAhvd/ZUpOQsREUnYmMHv7jeNsf0e4J4Rtj0OPD6xqomIyFTQnbsiIiGj4BcRCRkFv4hIyCj4RURCRsEvIhIyCn4RkZBR8IuIhIyCX0QkZBT8IiIho+AXEQkZBb+ISMiMOVePDO/BZw6dVvaByxaloCYiIuOjFr+ISMgo+EVEQkbBLyISMmMGv5ltNLN6M3t5hO0fNLMXzewlM/udmV0Qt+21oPwFM9uezIqLiMjEJNLivx9YN8r2V4Er3f084IvAhiHb3+7uq929amJVFBGRZEpkBa7NZrZklO2/i3u5ldii6iIikqaS3cf/YeBnca8d+IWZ7TCz9Un+LBERmYCkjeM3s7cTC/4r4oqvcPdaM6sEnjSz3e6+eYTj1wPrARYt0nh4EZGpkpQWv5mdD9wHXOfuJwbL3b02+LMe+AFw6Ujv4e4b3L3K3aui0WgyqiUiIsOYdPCb2SLg+8At7r43rrzAzIoGnwNXA8OODBIRkekzZlePmT0ErAUqzKwGuBPIBnD3e4HPA+XAt8wMoC8YwTMb+EFQlgU86O4/n4JzEBGRcUhkVM9NY2z/CPCRYcqrgQtOP0JERFJJd+6KiISMgl9EJGQU/CIiIaPgFxEJGQW/iEjIKPhFREJGwS8iEjIKfhGRkFHwi4iEjIJfRCRkFPwiIiGj4BcRCRkF/wS5O509/amuhojIuCn4J6C1q5cHnjnEP/10J3UtnamujojIuCRt6cUwcHdeONzMT16so7d/AAd2H21lbkleqqsmIpIwtfjHYffRVv5zRw3Roggf+4MVzC3JZX99W6qrJSIyLgkFv5ltNLN6Mxt26USL+bqZ7TezF83sorhtt5nZvuBxW7IqngqvHm8nK8P4i7cuI1oU4axoIYcaO+jpG0h11UREEpZoi/9+YN0o298NrAge64F/ATCzWcSWaryM2ELrd5pZ2UQrm2pHWjqZXZxLZoYBsLyykP4B57UT7SmumYhI4hIKfnffDDSOsst1wPc8ZitQamZzgXcBT7p7o7s3AU8y+hdI2nJ36pq7mFeae6psSXkBmRnGgQZ194jIzJGsPv75wOG41zVB2UjlM05LZy+dvf1vuJCbk5XBoln5HFA/v4jMIGlzcdfM1pvZdjPb3tDQkOrqnOZIcxcA80py31B+VmUhR1q6aO/uS0W1RETGLVnBXwssjHu9ICgbqfw07r7B3avcvSoajSapWslzpKUTA+YMGbq5PFoIoO4eEZkxkhX8m4Bbg9E9a4AWd68DngCuNrOy4KLu1UHZjFPX3ElFYYScrDf+lc0vzSOSlaHgF5EZI6EbuMzsIWAtUGFmNcRG6mQDuPu9wOPAe4D9QAfwoWBbo5l9EdgWvNVd7j7aReK0VdfSxaLy/NPKMzOMZdFCjecXkRkjoeB395vG2O7A7SNs2whsHH/V0kdHdx/Nnb2sGeEO3bOiBeyqO8nBE+0sLi+Y5tqJiIxP2lzcTWdHWmIXdueW5g67fUlFLOxfONw8bXUSEZkoBX8CBidimzdCi7+iMIIBBxp0I5eIpD8FfwLqWrooycumIDJ8z1h2ZgZlBTm6wCsiM4KCPwFHmjuZWzJ8N8+gaGFEN3KJyIyg4B9DZ08/Da3dzCsdferlyqII1cfb6R/waaqZiMjEKPjHsPvoSRzGbvEXRejpG6C2SQuziEh6U/CPYVddKzDyhd1B0aIIAPsbWqe8TiIik6HgH8PBxtgc/CX52aPuNxj8B+o1skdE0puCfww1jZ2U5meTYTbqfvk5WVQU5ugOXhFJewr+MdQ0dVCWn5PQvsuihRrSKSJpT8E/hpqmTkoTDP7lCn4RmQEU/KPo6OnjRHsPZWP07w86q7KQpo5eTrR1T3HNREQmTsE/isGhmYl29SyPxubs0dQNIpLOFPyjqDkV/Im3+AFd4BWRtKbgH0VNUwcApQWJtfjnleSRl52pfn4RSWsK/lHUNHUSycqgaITJ2YbKyDCWRQsU/CKS1hIKfjNbZ2Z7zGy/mX16mO1fNbMXgsdeM2uO29Yft21TMis/1WqaOplfloeNMYY/3nKtxiUiaW7MpqyZZQLfBK4CaoBtZrbJ3XcO7uPun4zb/2PAhXFv0enuq5NX5elzuKmDBWWnL7c4mrMqC/nxi0fo7OknLydzimomIjJxibT4LwX2u3u1u/cADwPXjbL/TcBDyahcqtU0dbKgbPQ5eoZaHi3EHaqPq9UvIukpkeCfDxyOe10TlJ3GzBYDS4Gn4opzzWy7mW01sz+ecE2nWXt3H43tPeMP/koN6RSR9JbYVcvE3Qg86u79cWWL3b3WzJYBT5nZS+5+YOiBZrYeWA+waNGiJFdr/GqbY0M5F5Tl09bVl/BxS8oLMINqXeAVkTSVSIu/FlgY93pBUDacGxnSzePutcGf1cCveGP/f/x+G9y9yt2rotFoAtWaWoNDOcfb4s/NzmR+aR6vHleLX0TSUyLBvw1YYWZLzSyHWLifNjrHzFYBZcCWuLIyM4sEzyuAtwA7hx6bjgZv3hpv8AMsrShQ8ItI2hoz+N29D7gDeALYBTzi7q+Y2V1mdm3crjcCD7t7/NqD5wDbzez3wNPA3fGjgdLZ4Bj+aGFk3McujxZS3dDOG/8qRETSQ0J9/O7+OPD4kLLPD3n9hWGO+x1w3iTqlzI1TR3jHsM/aGlFAW3dfTS0dVNZNPqSjSIi00137o7gcGMnC8c5hn/Q0orYyJ5qjewRkTSk4B9BTVPHhPr3AZYFs3Sqn19E0pGCfxht3X00dfSO+67dQfNK8ohkZWhIp4ikJQX/MGonMaIHYpO1aWSPiKQrBf8wJjqGP97SigKqFfwikoYU/MN4fQz/xLp6IBb8h0500Ns/kKxqiYgkhYJ/GDVNHUSyMqgoTGwBluEsixbSN+CnvkRERNKFgn8Yg7NyTmQM/6DXh3TqAq+IpBcF/zBiwT/xbh54feF1XeAVkXSj4B/GZMbwDyrNz6EsP1vTM4tI2lHwDzHZMfzxlkULeVULsohImlHwD5GMoZyDllYUaNoGEUk7Cv4hahpjo3AWzpp8i39pRQH1rd20dSe+kIuIyFRT8A+RzBb/4AXe13SBV0TSiIJ/iJqmTnKzMygvmPgY/kFLKwoBOKAhnSKSRhT8QwwO5ZzMGP5Bi8vzycww9h1T8ItI+khoIRYzWwd8DcgE7nP3u4ds/zPgy7y+Fu897n5fsO024HNB+T+5+3eTUO8pU9M88aGcDz5z6LSyZRUF7D7aOtlqiYgkzZjBb2aZwDeBq4AaYJuZbRpmCcX/cPc7hhw7C7gTqAIc2BEc25SU2k+BmqZOVi8sTdr7rZpbzHMH0/Z0RSSEEunquRTY7+7V7t4DPAxcl+D7vwt40t0bg7B/Elg3sapOvdauXpqTNIZ/0Dlzi6ht7uRkV2/S3lNEZDISCf75wOG41zVB2VDvM7MXzexRM1s4zmMxs/Vmtt3Mtjc0NCRQreSrbZ7cPPzDOWdOMQB71N0jImkiWRd3fwwscffzibXqx92P7+4b3L3K3aui0WiSqjU+g2P4k9niXzW3CIBddSeT9p4iIpORSPDXAgvjXi/g9Yu4ALj7CXfvDl7eB1yc6LHpJJlj+AfNKc6lND+bXXVq8YtIekgk+LcBK8xsqZnlADcCm+J3MLO5cS+vBXYFz58ArjazMjMrA64OytJSMsfwDzIzVs0pYvdRtfhFJD2MOarH3fvM7A5igZ0JbHT3V8zsLmC7u28CPm5m1wJ9QCPwZ8GxjWb2RWJfHgB3uXvjFJxHUiRzDH+8VXOKeWT7YQYGnIyM5L63iMh4JTSO390fBx4fUvb5uOefAT4zwrEbgY2TqOO0mcwY/tGcO7eYjp5+DjV2sCRYoEVEJFV0526cw42dUxL8gxd41d0jIulAwR842dVLS2cvC5M4omfQytlFZBi6wCsiaUHBH6htSv5QzkG52ZksrSjQkE4RSQsK/kBNU/Jv3oq3am6x5uwRkbSg4A9MxRj+eOfMKeJQY4cWZRGRlFPwBw43dpKXncmsJI7hj3fO3MGpG9TdIyKppeAPVB9vY2lFQdLH8A9aFQT/Tl3gFZEUU/AH9h1rY8Xswil7/3klucwqyOGFQ81T9hkiIolQ8ANt3X3UNneyonLqgt/MWLNsFlurT+DuU/Y5IiJjUfADB+pjSyOumF00pZ+zZlk5tc2dHA5mARURSQUFP7D3WKzffSpb/ACXLysHYGv1iSn9HBGR0Sj4gf31beRkZrBoVvJv3op3VmUhFYU5bFHwi0gKKfiBffVtLIsWkJU5tX8dZsZly8rZckD9/CKSOgp+Yl09U92/P+jyZeUcPdnFayc6puXzRESGSmha5jNZR08fNU2dvL9q4dg7J8Hly1/v51+qKZpFZoQHnzl0WtkHLluUgpokR0LBb2brgK8RW4jlPne/e8j2vwE+Qmwhlgbgz939YLCtH3gp2PWQu1+bpLonxf5gRM/KKRzDH29ZRQHRoghbDpzgpktn7g+OyEzR2z9ATVMnrx1vp7a5k8b2Hhrbe+jo6aO2qZPcnEyKI9msnFNESV52qqs7LcYMfjPLBL4JXAXUANvMbJO774zb7Xmgyt07zOyvgP8D3BBs63T31Umud9LsOxYL/rMqp6erx8y4fFk5W4Lx/FN1p7BI2PT1D/DaiXb2HG1jz7FW9h1rZe+xVg6e6KBv4I3X1HKzM8jJzKCrb4CevoFT5QvL8jh/QSmXLZ015df8UimRFv+lwH53rwYws4eB64BTwe/uT8ftvxW4OZmVnEr76tvIzjQWl0/tiJ54a5aVs+n3R6g+3s7y6PT8piFypnB3aps72XcsFvB7j7ay+2gr+xvaToW4AbMKcphdnMsVZ1VQXhihojCHsvwc8iOZZGW8HuoD7hxv62bnkZO8fKSFn75Ux/OHmrjxkkVUFEVSdJZTK5Hgnw8cjntdA1w2yv4fBn4W9zrXzLYT6wa6291/OO5aTqF9x1pZVlFI9jR+uw/28285cELBLzKCgQHnSEss4PfVt7LvWBt769vYf6yV9p7+U/sV52YxpySXy5bMYnZJLrOLc6ksiiT8fzrDjMqiXCrPzmXt2ZXsPHKSx56r4Z6n93PtBfO4aHHZVJ1iyiT14q6Z3QxUAVfGFS9291ozWwY8ZWYvufuBYY5dD6wHWLRo+vq+99W3cd6Ckil7/+EuCrk780vz+MXOY9y8ZvGUfbbITDAw4NQ0dbL3WCv76l8P+f31bXT2vh7wRZEsosURzl9QSmVxhNlFsZDPy8lMan3OnVfM/LIVPLL9MI8+V0NXXz//Y3lFUj8j1RIJ/logfsjLgqDsDczsncBngSvdvXuw3N1rgz+rzexXwIXAacHv7huADQBVVVXTMsi9s6efw00d/MlF86fj404xM9530Xy+8fR+jjR3Mq90atYAEEkn7k5dS9ep7plYP/zpAV+cm0VlUS4XLiqlsiiXaFGE2cUR8nOmbxBiSV42H75iKQ88c4ifvlhHRWGEldM05Hs6JPI3uQ1YYWZLiQX+jcAH4ncwswuBbwPr3L0+rrwM6HD3bjOrAN5C7MJvWjjQ0IY7KfkHvf7ihXz9qf08tqOGj71jxbR/vshUau3qZc/RVnYdbWV33Un2BEHf2vX6QkSVRRGKc7O5aFEplUH3TGVR8lvwE5VhxvurFrBhczUPPXuIj165nNnFuamuVlKMGfzu3mdmdwBPEBvOudHdXzGzu4Dt7r4J+DJQCPxnMEplcNjmOcC3zWyA2M1idw8ZDZRS++qnZ46e4Swqz+fyZeU8suMwt7/9LDIyNLpHZp6BAedQYwe76k6y62gru+pOsvvoyTdMRJibncGc4lzOnVvM7OLc4DG9LfiJimRlcsuaxXzrVwf43pbXuH3tWeRH0r/eY0noDNz9ceDxIWWfj3v+zhGO+x1w3mQqOJX2HWsjK8NYXJ6aG6luuGQhn/iPF9j66okzrg9RzjxN7T3sPtrKnqMn2XOslV11seGSHcGF1gyDJRUFnL+glHPmFDOnOJc5JbmU5GXP6GHLpfk53LJmMRs2V/OTl+qm7WbPqTTzv7om4aXaFpZHC8nJSs143XVvnkPRj7L4z+01Cn5JG509/eyrjw2RHOyL33O0lfrWU5fuyMvOZG5JLqsXlp4K+NnFudM6Om46LZyVz5VnR3lqdz0XLCjl7Dkzu78/tMHf1dvPs682pvS269zsTK69YB6P7qjhC9e+KTR3DUp6GBxNs+voSXbXtbL76El2H23ltRPtDM4hmJVhzC7OZUFZHhcvLmNOcS6zS3IpimTN6Fb8RKxdGeWl2hZ++EItn5jh1+VCG/zPvtpId98Ab1sZTWk9brhkIQ88c4hNvz/CLRraKVOksb0ndoE16KYZbM0PjocfvOFpTkkubz+7MtaKL85lVmEOGSEL+JFkZWbwvgvn8+3N1Tyx8xgfumJpqqs0YaEN/s17G8jJymDN0vKU1uO8+SW8aV4x9/7qANdftCBtRjTIzOPu1Ld2c6C+jf0NsWGS++vb2HusjeNtb+ymmVOSy/kLSplTEgv42cW5KevynEkWlRewZnk5Ww+cYMfBRi5ePCvVVZqQ8Ab/vgYuXTIr5UFrZvzDNedy44at3PP0Pv7uXatSWh9Jf129/bx6vJ0DDW1UN7RT3dBG9fF2qhvaaet+fbhkJCuDyqIIi8vzuXRJGZVBK74oN3zdNMl09bmz2XnkJJ/9wcv85GNXzMg5fUIZ/HUtnew91sb1Fy9IdVWA2Nw9771wPhs2V/O+ixawTNM4CNDd18+B+nb2BpON7T3WxnOHmmhq7yH+DseSvGyihRHePL+EaFGEaGGEaFGEYgX8lIhkZXLN+XN54JlDfHfLQT48A7t8Qhn8v9l3HIC3rkhd//7QqRxWzSnilzszuXPTK3zvzy/Vf9iQaWjtZmfdydh4+LqTbDlwguNt3QxOKplhUFEYYV5pHqsXlp4K+IrCiLpoUuDcucWsPTvKV36xhz88by5zSmbWjV2hDP7NexuoLIqwKo2GZBXlZvOpq1fyhR/v5Kcv1XHN+fNSXSWZAv0DzqvH20+F/M4jJ3nuYBOtcV00JXnZr9/wFPTBlxfmvGFGSUktM+Mfr30TV311M//0053c84GLUl2lcQld8PcPOP+9/zjvWDU77VrVN69ZzKPP1fDpx15i0ax8zl9QmuoqySScaOtmTzBl8OBQyT1HW+kOpg7OzjTOqixixewi5pbExsLPLcmdEXe0CiwuL+D2tWfx1V/u5f1VDSkfITgeofsJe6m2heaOXt62Mv1umMrKzOBfb63iT+/dwq0bn+Xh9WtYNac41dWSMTR39MRmlTzWdqo//sWaljdcaC3IyWRuSR6XLJl1KuCjRRG14me4v7xyGZt+X8vfP/oiT3zibZTkz4x7cUIX/Jv3NmCW2v790cwtyePBj6zhT7/9O26+71ke+cs1utibJk60dQfTBsfmhB983hB3R2tOZgaVxRHOnl0UzA0fYU5xLoUhvOEpDHKzM/nnGy7kvd/6LZ/94Ut846YLZ8S/c6iCv7d/gB+9UMv580uYVZCT6uqMaFF5Pg98ZA03fHsLf/Ivv+Ou697MH50/d0b8QM107s7Rk12nxsDHQj42Lr6xvefUfoNDJReW5VO1uOzUzJIl+dm64SlkzltQwievWsmXn9jDO8+ZzR9fOL3TvE9EqIL/37ce5EBDO/96a1WqqzKsoSN9brt8Cf+1p56PP/Q8P3+5ji9e92bKC8/MpeCmW0/fAIca29lfHxsPf6ChjQP1bRwYMhY+LzuTyqIIyyoKuHxZOZVFsaGSM33iMUmuj165nKd31/MPP3yZqiVlLCibvqVcJyI0wd/Y3sNXn9zLW1dU8M5zKlNdnYRUFEV47KOX8+3N1fzzL/fym73HufHShfzZW5YyX4u3jMndOXaym+rjbbx2vOPUjU6vHm/nUGMH/XELcBfnZhEtio2Fj7XeYwGvLhpJRGaG8dUbVvPur/2GD/2/bTy0fg0VadxIC03wf/XJvbT39PMP15w7o/4jZ2VmcPvbz+Kqc2dzz1P72fjb19j429e4+tzZrHvzHNaurJwxF5SmQm//AHXNXRxq7OBQYwcHG9s5eLyD1060c/BExxtWdopkZVCWn0NFUYS3raigIrjRqaIwQm62psqQyVk4K59/vbWKD93/LDff9wwP/sWatO1SDkXw7z56kgeeOcgtaxbPuOXT4rt/1iwrZ9WcIpo6evjB80f42ctHycwwLl5UxoWLS7lwYSkXBNPkzqQvt5G4Oyc7+zh6sou6lk7qWro40txJbXMntU2d1DR1UtfSSVzDncwMi4V7YQ4XLSqlPLjJqaIwh+I89b/L1Lp8eTnfue0S/vz+bXzwvmd48COXUZaG4W/uYy9va2brgK8RW4HrPne/e8j2CPA94GLgBHCDu78WbPsM8GGgH/i4uz8x1udVVVX59u3bx3cmI3j1eDu3P/AcR1o6+dXfrqU0f/z/CMMtmJ5qA+7UNnWyq+4k+xvaqGvuoj/4tyyKZLG8spDl0UIWzspjQVk+80vzmF2c2u4Ld6etu4/mjl5aOntp7uilsaOHpvYeTrT3cLytm+Ot3Rxv66a+tZuG1u5TY94HGVCcl01pXjZlBTmU5WdTlp/DrILYQ+Eu02W0Kd1/vbeBv/judsoKsvncH57LNdMwOMPMdrh7Qhcwx2zxm1km8E3gKqAG2GZmm4YsofhhoMndzzKzG4EvATeY2bnE1uh9EzAP+KWZrXT3fqbYwIDzb1sP8r9/touczAy+esPqCYV+usowY+GsfBbOyudqgi6Pli5qmztpaO2ivrWbJ3cepbW7j6Hf7XnZmZTlZ1OSn0NJXhaFkWwKI5nkR7LIzcokkp1BJCuD7MwMMszIyjAGf2bdY186/e709zu9A05P3wA9fQN09/XT2dtPV28/7d39tHf30dbdR3tPHyc7+2jt6n1D6zyeAXk5mRRGsiiMZFFRGGFpeQFFuVkU52VTkpdNcV42xbnZZGqZSklzV66M8sgJ+0PYAAAIOElEQVRHL+ezP3iJjz30PA89e4i/uWolFy4qS4uf30S6ei4F9rt7NYCZPQxcB8QH/3XAF4LnjwL3WOzr7TrgYXfvBl41s/3B+21JTvVf5x5bVOKFw808f6iZLdUn2FV3kitXRvnS+86fcXNpjFd2ZgaLZuWzaNYbRxP09Q/Q0tlLU0cvrV29tHX30drVR2dPPx29/dS1dNHT10F33wDdvf30DTh9/X7qt4dEZGbEvhyyMoyc4AsjJyv2iGRlUl4QYX5pPrnZGeRmZZKfE3vk5WRRkBP7wsnLzkyL/xAiybJ6YSmb7riCB589xJd/vpvr791CWX42V66McvHiMuaV5jG3JI+KohzysjPJzc6cthXMEgn++cDhuNc1wGUj7RMszt4ClAflW4ccOyWDXHv7nXd85df09A0Qycrg/AUlfOl95/H+qoVnRH/3RGVlZlBeGBn3MNAB99hjIPY8nlnsN47Yg1D//YqMJjPDuGXNYq5bPY/Next4anc9v97TwA9fODLs/tGiCNs+O+wS5kmVNhd3zWw9sD542WZmeybzfnuJ/epx42QrFlMBHE/OW6WFM+184Mw7pzPtfOAMO6cPTsF7HgTscxM+POEl/BIJ/logfln5BUHZcPvUmFkWUELsIm8ixwLg7huADYlVe3qZ2fZEL5rMBGfa+cCZd05n2vnAmXlOM1UiHUrbgBVmttTMcog1ojcN2WcTcFvw/HrgKY8NF9oE3GhmETNbCqwAnk1O1UVEZCLGbPEHffZ3AE8QG8650d1fMbO7gO3uvgn4DvBvwcXbRoIelmC/R4hdCO4Dbp+OET0iIjKyhMbxh52ZrQ+6os4IZ9r5wJl3Tmfa+cCZeU4zlYJfRCRktAqEiEjIKPhHYWbrzGyPme03s0+nuj6TZWYbzazezF5OdV2SwcwWmtnTZrbTzF4xs79OdZ0my8xyzexZM/t9cE7/mOo6JYOZZZrZ82b2k1TXRRT8I4qbquLdwLnATcEUFDPZ/cC6VFciifqAT7n7ucAa4PYz4N+oG/gDd78AWA2sM7M1Ka5TMvw1sCvVlZAYBf/ITk1V4e49wOBUFTOWu28mNurqjODude7+XPC8lViwpP/yR6PwmLbgZXbwmNEX4sxsAfCHwH2provEKPhHNtxUFTM6VM5kZrYEuBB4JrU1mbygW+QFoB540t1n+jn9M/D3wMBYO8r0UPDLjGdmhcBjwCfc/WSq6zNZ7t7v7quJ3el+qZm9OdV1migzuwaod/cdqa6LvE7BP7KEp5uQ1DGzbGKh/4C7fz/V9Ukmd28GnmZmX5d5C3Ctmb1GrLv0D8zs31NbJVHwjyyRqSokhYKpv78D7HL3r6S6PslgZlEzKw2e5xFbB2N3ams1ce7+GXdf4O5LiP0fesrdb05xtUJPwT8Cd+8DBqeq2AU84u6vpLZWk2NmDxFbC+FsM6sxsw+nuk6T9BbgFmKtyBeCx3tSXalJmgs8bWYvEmt8POnuGgIpSaU7d0VEQkYtfhGRkFHwi4iEjIJfRCRkFPwiIiGj4BcRCRkFv4hIyCj4JeXMrD8Yg/9KMB3xp8wsI9hWZWZfH+XYJWb2gemr7Wmf3RnMq5MWzOyGYBpxjf2XESn4JR10uvtqd38TsTtV3w3cCeDu293946McuwRISfAHDgTz6iQsmPJ7Srj7fwAfmar3lzODgl/SirvXA+uBOyxm7WDr1cyujLtD93kzKwLuBt4alH0yaIX/xsyeCx7/Izh2rZn9ysweNbPdZvZAMOUDZnaJmf0u+G3jWTMrCmbI/LKZbTOzF83sLxOpv5n90Mx2BL+9rI8rbzOz/2tmvwcuH+Ez3xQ8fyH4zBXBsTfHlX978IsjWCjoueA9/iuJ/wxypnN3PfRI6QNoG6asGZgNrAV+EpT9GHhL8LwQyIrfHpTnA7nB8xXA9uD5WqCF2GR7GcSmrrgCyAGqgUuC/YqD910PfC4oiwDbgaVD6rgEeHlI2azgzzzgZaA8eO3A+4PnI33mN4APxu2TB5wTnHd2UP4t4FYgSmza8KXxnxt3rj8Z7u9aDz3cnaxxfk+IpNJvga+Y2QPA9929Jmi0x8sG7jGz1UA/sDJu27PuXgMQ9MsvIfZlUOfu2wA8mNbZzK4Gzjez64NjS4h9kbw6Rh0/bmbvDZ4vDI45EdTlsaD87BE+cwvw2WDhku+7+z4zewdwMbAtONc8YvP0rwE2u/urwXucMQvsyNRT8EvaMbNlxIKynliLFwB3v9vMfgq8B/itmb1rmMM/CRwDLiDWsu+K29Yd97yf0X/+DfiYuz8xjnqvBd4JXO7uHWb2KyA32Nzl7v2jHe/uD5rZM8RWq3o86F4y4Lvu/pkhn/VHidZLZCj18UtaMbMocC9wj7v7kG3L3f0ld/8SsZkrVwGtQFHcbiXEWtMDxGbuHOtC6h5grpldEnxGkZllEZuV9a+C+f4xs5VmVjDGe5UATUHoryLWKk/4M4MvvGp3/zrwI+B84L+A682sMth3lpktBrYCbzOzpYPlY9RN5BS1+CUd5AVdL9nEFlD/N2C4+fU/YWZvJ7aE3yvAz4Ln/cFF0/uJ9YE/Zma3Aj8H2kf7YHfvMbMbgG8E8993Emu130esK+i54CJwA/DHY5zHz4GPmtkuYuG+dZyf+X7gFjPrBY4C/8vdG83sc8AvgiGuvcDt7r41uHj8/aC8ntiIKJExaVpmkQmy2Dq/P3H3tFoaMehy+lt3vybVdZH0pK4ekYnrB0rS7QYuYr/1NKW6LpK+1OIXEQkZtfhFREJGwS8iEjIKfhGRkFHwi4iEjIJfRCRk/j+Zel4o5Vf7MAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nb_merge_dist_plot(\n", " SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n", " SkyCoord(cfht['cfht_ra'], cfht['cfht_dec'])\n", ")" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Given the graph above, we use 0.8 arc-second radius\n", "master_catalogue = merge_catalogues(master_catalogue, \n", " cfht, \n", " \"cfht_ra\", \n", " \"cfht_dec\", \n", " radius=0.8*u.arcsec)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(np.sum(master_catalogue['flag_merged'])/len(master_catalogue)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add HSC-PSS" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl03Gd97/H3dxaNlhlJ1mp5t+UtJglZnB1KwtYEaEJPKVuhvS00twuUtvT2wC2X3nLPvQfacykt0NtyCmUpIYQEaICQQBOHOLHj2E5CYsd2bHmVvGixJGuf7Xv/+P1GUWzZGkkzmplnvq9zdDL6zU/zeya2P3rm+zy/5xFVxRhjjFsChW6AMcaY3LNwN8YYB1m4G2OMgyzcjTHGQRbuxhjjIAt3Y4xxkIW7McY4yMLdGGMcZOFujDEOChXqwk1NTbpq1apCXd4YY0rS7t27e1W1eabzChbuq1atYteuXYW6vDHGlCQROZbNeVaWMcYYB1m4G2OMgyzcjTHGQRbuxhjjIAt3Y4xxkIW7McY4yMLdGGMcZOFujDEOsnA3xhgHFewO1UK6Z8fxaY+//4YVC9wSY4zJD+u5G2OMgyzcjTHGQRbuxhjjIAt3Y4xxkIW7McY4qOzCff/pc6TSWuhmGGNMXpVVuA+MxnnbP2zlR788WeimGGNMXpVZuCdIKzxz9CwdPcOFbo4xxuRNWYX78EQSAAF+8FwX8WS6sA0yxpg8KatwH/HD/dYNLZwdifPzl04XuEXGGJMfZRXuo/EUABsXx7hhdQPbOvo4cXa0wK0yxpjcK6twz5RlKkIBbn/NYsCbPWOMMa4pq3AfjXvhHgkFiISDRMIBxhNWdzfGuCercBeR20XkgIgcEpFPTPP8ChHZIiLPicgLIvK23Dd1/oYnvLJMRch725XhIOOJVCGbZIwxeTHjkr8iEgS+DLwF6AR2isiDqvrSlNM+Bdynqv9PRDYBDwGr8tDeWZu6vO/2jj4AIqEgAJWhIOM2Y8YY46Bseu7XA4dU9bCqxoF7gbvOO0eBWv9xHVCUdwnFkylCASEYEAAqwwHruRtjnJTNZh1LgRNTvu8EbjjvnP8J/ExEPgrUAG/OSetybCKZnizJgNeDHxpPFLBFxhiTH7kaUH0f8HVVXQa8DfiWiFzw2iJyt4jsEpFdPT09Obp09uLJNJEp4V4ZDlhZxhjjpGzCvQtYPuX7Zf6xqT4E3AegqtuBSqDp/BdS1a+o6mZV3dzc3Dy3Fs/D+T13G1A1xrgqm3DfCawTkdUiUgG8F3jwvHOOA28CEJHL8MJ94bvmM/B67sHJ7zPhrmqrRBpj3DJjuKtqEvgI8AiwD29WzF4R+YyI3Omf9nHg90Xkl8B3gP+iRZiYE8nUq8syoQBphUSq6JpqjDHzks2AKqr6EN70xqnHPj3l8UvALbltWu5NJNPUVoUnv4+Eg/7x1KvKNcYYU+rKKtEuHFD1wt3uUjXGuKaswv3CAVXvsQ2qGmNcU1bhfsGAqv94PGnhboxxS9mEezKVJqV6wVRIsLKMMcY9ZRPumV2Xzr+JCWDCyjLGGMeUTbhPTBvumZ67hbsxxi1lF+4VU2ruFaEAArYEgTHGOWUT7nF/0HRqzz0gQkXIVoY0xrinbMJ9sucefPVb9pYgsJ67McYtZRfukfD54W49d2OMe8om3F+ZLRN81XFvNyYLd2OMW8om3CeSr94/NaMyHGTCyjLGGMeUTbhPN88dvDKNlWWMMa4pm3CfSKYRIOTvn5phG3YYY1xUPuGeShMJBxA5L9xDttWeMcY9ZRPu8UT6gmmQ4PXcU2klkbKAN8a4o2zCfSKVvmCmDLyyYYeVZowxLimbcI8nUxfMcQevLAPYjBljjFPKJtwnLlGWAVvT3RjjlrIJ93gqfcE0SLA13Y0xbiqbcD9/i70M22rPGOOisgr36QZUM1vtTVhZxhjjkLIJ93gydZGeu5VljDHuKYtw9+ax67Q194iVZYwxDiqLcM/coDRduNuGHcYYF5VFuE+3xd5UtgSBMcY1ZRLuF26xN5UtHmaMcU1ZhHt8sud+8XC3O1SNMS4pi3CfuMha7hmV4YDdoWqMcUpZhPvFttjLiISsLGOMcUtZhPvFttjL8GruVpYxxrijTMI9i7KM9dyNMQ4pi3DPZkA1mdbJ84wxptSVRbhPzBTu/vGh8cSCtckYY/KpLMI9nkwTDgqB8/ZPzcjsxjQ0nlzIZhljTN6URbhfbEXIjMzKkBbuxhhXlEm4py46mAqvrOluZRljjCvKItzjF9moIyOz7O8567kbYxyRVbiLyO0ickBEDonIJy5yzrtF5CUR2Ssi9+S2mfPjlWVmDnfruRtjXBGa6QQRCQJfBt4CdAI7ReRBVX1pyjnrgE8Ct6hqv4i05KvBcxFPpqmJXKrm7gX/8IT13I0xbsim5349cEhVD6tqHLgXuOu8c34f+LKq9gOoandumzk/Mw2o2mwZY4xrsgn3pcCJKd93+semWg+sF5GnRORpEbk9Vw3MhfgMA6rBgBAOipVljDHOmLEsM4vXWQfcCiwDnhCRK1R1YOpJInI3cDfAihUrcnTpmU3MMKAKXt3deu7GGFdk03PvApZP+X6Zf2yqTuBBVU2o6hHgZbywfxVV/YqqblbVzc3NzXNt86yoessKXKrnDt7KkENWczfGOCKbcN8JrBOR1SJSAbwXePC8c36I12tHRJrwyjSHc9jOOUukFOXiW+xlVIYD1nM3xjhjxnBX1STwEeARYB9wn6ruFZHPiMid/mmPAH0i8hKwBfhvqtqXr0bPxkxb7GVUhoIMW83dGOOIrGruqvoQ8NB5xz495bECf+5/FZWZVoTMiFjP3RjjEOfvUI2n/HAPZtFzt5q7McYRzod7Isuee2U4wLD13I0xjnA+3OMpBSA8Q889Eg4yHE+STutCNMsYY/LK+XBP+GWZcHD6tdwzKkMBVGEkbr13Y0zpK6Nwn7nnDrYEgTHGDWUQ7l6ZZcYBVT/cbVDVGOOCMgh3r+cemqEsE7F9VI0xDimbcJ95KmQm3K3nbowpfc6He3yy5241d2NM+XA+3BNJJRgQgoEZZstYzd0Y4xD3wz2VnnEaJEwty1jN3RhT+soi3GeqtwOEQwFEsLtUjTFOcD7c46n0jHPcAQIiRCMhW9PdGOME58M9kdKswh0gFgnZgKoxxgnOh3syy5o7QKwybGUZY4wTnA/3eCpNeIYVITOilSGGJmxA1RhT+pwP90QqTTiQZVmmMmQ9d2OME9wP96Rm33O3mrsxxhHuh3sqTUXWNXebLWOMcYPz4Z7tVEjwBlTtJiZjjAucD/fELMI9GgkxnkhPLjZmjDGlyulwV9XZzXOvDAF2l6oxpvQ5He4Tmc2xs6y5RyN+uFvd3RhT4pwO97F4CiDr2TKxyjAA56zubowpcU6H+3jSD3cryxhjyozT4T7Zc59tuFtZxhhT4twO90Qm3GdXc7cbmYwxpc7pcB9PzLbn7tXc7UYmY0ypczrcx+LebJnZlmXsRiZjTKlzO9z9nns2OzEBREIBQgGxAVVjTMkri3DPtuYuIt76MhbuxpgS53S4j89ynjt4a7rbbBljTKlzOtzHZjmgChCL2OJhxpjSVxbhnm3NHfzdmKwsY4wpcU6He2YqZCjLmjtArYW7McYBTof7WCJFKCAEJPtwj0as5m6MKX1Oh/t4PDWrejvYhh3GGDc4He5jiVTW0yAzMrNlVDVPrTLGmPzLKtxF5HYROSAih0TkE5c47zdEREVkc+6aOHdjiex3YcqIRkIkUjq5FrwxxpSiGZNPRILAl4E7gE3A+0Rk0zTnxYCPATty3ci5GounqJjFHHfwBlTBFg8zxpS2bJLveuCQqh5W1ThwL3DXNOf9L+BzwHgO2zcv44nZ19yjtuyvMcYB2STfUuDElO87/WOTROQaYLmq/iSHbZu3udTcYxF/ZUgbVDXGlLB5D6iKSAD4PPDxLM69W0R2iciunp6e+V56RmNzmC0Ttd2YjDEOyCb5uoDlU75f5h/LiAGXA4+LyFHgRuDB6QZVVfUrqrpZVTc3NzfPvdVZmktZJrPs7zkLd2NMCcsm+XYC60RktYhUAO8FHsw8qaqDqtqkqqtUdRXwNHCnqu7KS4tnYTyRmtXSA/BKWcZq7saYUjZj8qlqEvgI8AiwD7hPVfeKyGdE5M58N3A+xhIpwqFZ1txtww5jjANC2Zykqg8BD5137NMXOffW+TcrN8YSKcIBq7kbY8qPs3eoptPKeCI9q7XcwVseuDIcsH1UjTElzdlwz9xhOtsBVYBoJGw3MRljSpqz4T7bLfamqq0KMTgWz3WTjDFmwTgf7rOdLQPQVBOhb9jC3RhTutwN9/jst9jLaIxWcHbEwt0YU7qcDffxOeyfmtEYraDPwt0YU8KcDffJmvss57kDNNZE6B+Nk0zZsr/GmNLkbLiPz6fmHq1AFfpH7UYmY0xpcjbcMzX30JzKMhEA+kYmctomY4xZKO6G+zymQjbWVADYjBljTMlyNtznU5ZpjHrh3jtsPXdjTGlyNtznNRWyxi/LWM/dGFOislo4rBSNJea+/EBdVZhgQGyuuzFmWvfsOH7BsfffsKIALbk4h8M9M6A6+5p7ICA01FTYgKoxZUpV6ewfo2tgjIAIwQAsb6imJVZZ6KZlzdlwH0+kqAwHCMjswx28QdVeK8sYUzY6+0fZerCX7R197Dx6llOD4696viIU4KO3rZ2cTVfsnA33sXiKqnBwVj8z9aNWMqUcOD3EPTuOF93HLWPM/I0nUuw4cpYt+7v58QsnJztzsUiIVU01bF65iGa/p55Ipfne7hM88GwnH379mjl3GheSu+GemH24T1UTCXK233ruxrika2CMLfu72bK/m20dfYz5n/BXNFRzw+pG1rZEaYlFkGnC++1XLOGBZzvZ3tHHLWubCtD62XE23McTKSor5h7u0UiIEduww5iSlkil2X2sny0Hunl8fw8HzgwBsKg6zGuX17GhtZY1zTVZTby4ZkU9e7oG+dlLp9nQGqMpVtzlGbfDPTSfnnuIiWSahK0vY0xJOTU4xi8O9PD4gR6eOtTL0ESScFC4blUDf3XtZYzEkzRHp++dX4qI8OtXL+ULj77MA892cvevrJn1aywkZ8N9LJGiah4995qI97/Geu/GFLfReJIdR87y5MFeth7s4eUzw4A3pXljW4wNrTHam6NE/DJt5t/2XNRWhbltQws/3XOagbEEi6orcvIe8sHdcI+nqK6Y+9uLToZ7KldNMsbkQDKV5sWuQZ461MvWg73sOtpPSpVQQFjdVMMdly9mXWuM1ovUzudrZUM1AKcGxi3cC2EskaahZv4992HruRtTUKrK0b5Rth7sYevBXp7u6JvcwP41S2q5eW0ja5ujrGrKrnY+X611lQhe+WfTktq8X2+unA338XmWZaJWljGmYMbiKbZ19HoDoQd66OwfA7yB0I1tXpllTXN08t/pQoqEgjRGKy6YB19snA13b5773H+L10S8XwzWczdmYZw4O8qWA908tr+bJw/2kkwrFcEA7c013PnaJaxridJQU1EUg5htdVV09o8WuhmX5G64z3Oee0UwQDgo1nM3Jk+GJ5LsONzHVn8gtKNnBIBVjdVcv7qBDYtjrG6smdOeDPnWVlfJi12DXidyHhWCfHI63Oczz11EqImErOduTI6MJ1LsPHqWbR19bO/o48WuQVJpJRz0BkLfdkUbG0tg/jh4PXeA0+fGWd1UU+DWTM/JcE+llXgyPa+eO/g3MsUt3I2ZC1XlwJkhtuzvYevBHp45cpZkWgkILFtUzevXNdHeHGVlQ3VR9s4vpa3eW5bg1OCYhftCmkh60xfnG+41FSGGJmwfVWOyNZFMsa2jj5+/dIYt+7snBx03Lo5x45pG2pujrGqqJjKPGwyLQSwSoiYS4tRA8Q6qOhnumY06Kucb7pEQp88V7x+eMcVgYDTOlgPd/PylM/ziQA8j8RQVwQBrW6LctKaRda0x6qrChW5mTokIbXWVnBocK3RTLsrNcE+80nNPpnXOrxONBBmZSKKqRTFCb0wxUFVeOnWOX7zcwy8O9LDrWD+ptNIci3DnVUupCAprmqMLMue8kNrqKtnW0UcqrQQDxZcPToZ7Zv/Uyoogw+Nzr5nXREIk08rwRJJYpVs9D2OypaocPzs6ORC6raNvcn/htrpKXr+2icvaalm6qKoklsLNlba6KlJppWdogsV1xbeJh5PhPhb3FvuqCs8v3DM3SPQNxy3cTVnpG55g68FenjrUy7aOProGvPJDLBJidXMNb1jfzLqWKLWOlVtmo63ulUFVC/cFMurPcJn3gGom3EcmWFWkI+LG5EI6rew5Ocij+7p5/EA3L3QNouotvnXTmkauWbmI9uaaOa2m6KqmaIRQQDg1OM7VhW7MNJwM98Exb4ZLfXWY42fn/jqZnrttt2dcNDyR5KlDvWzZ790V2j00gQDLFlXxpo0trG+NsaS+vEotsxEMCIvrKjlZpIOqTob7gB/u8x2hr5lSljGm1KXTyr7T59h6sJf7dp7gWN8oKVUioQDrWqL8yvpmNrTG5rUkbrlpq6tkT9c5VOc+cSNfnPxTHBx9pec+HzX+Ha59/uCRMaWmfyTOEwe9jSu2HuyZ/BTaWhvhlrWNrG+NsaKxmlDA7Zkt+bK4roqdR/s5N4+xvXxxMtwHxuIEAzLvFeNCwQCV4QB9I9ZzN6VBVdl3aojH9p/hsf3dPHd8AAWqK4Ksa4ly64YW1jaX90BoLjX467kPjBZfRrgZ7qMJ6qvCORn4qakIWbibotY/EufJQ7088XIPTxzs4cw575PmlcvquG1jCxtaY2U3TXGhZKoDmVJwMXEz3McS1M2zJJNRWxXm5EBxDpiY8hRPptl19CxbD/Xy5MFe9pz0ZrZUhYO0t0S5pb2JDYtjNn13AWTG9TKl4GKSVbiLyO3APwBB4F9V9bPnPf/nwIeBJNAD/J6qHstxW7M26Pfcc6EpGuFQ91BOXsuYuTpzbpzH9nezZX83j7/cQzyZJiCwoqGaN21sYW1LjGXWO19wleEgleEAA2PF9+l+xnAXkSDwZeAtQCewU0QeVNWXppz2HLBZVUdF5A+BvwXek48GZ2NgLE5LLDc3FTTHIuw8epazI3Eaaop3v0TjlnRaebFrkMf8aYovdg0CsLS+iquW1bO+NUZ7c83kps+mcOqrKhgo0Z779cAhVT0MICL3AncBk+GuqlumnP808IFcNnK2BkYTrG+J5eS1mqPe2tKHe4ZpqGnIyWsaM52h8QRbD3rzzrcc6KF32Jt3vryhml/d1MqGttq8bfps5q6+Ojx5b00xySbclwInpnzfCdxwifM/BPx0uidE5G7gboAVK1Zk2cTZGxzNXc292d84oKNnmM2rLNxN7qgqL58Z5vED3Xx35wmO9o2QVqgMB1jXEuO2Dc2st3nnRa++OsyxvuLbci+nf2tE5APAZuAN0z2vql8BvgKwefPmvMz6T6TSDE0kqa/KTQmlvjpMRSjAYX8LMGPmY3A0wVMd/syWl3s46a933lob4XVrm9mwOMaKhuqiXGXQTK++qoKxRIrhiWRBNuy+mGxa0gUsn/L9Mv/Yq4jIm4G/At6gqgW76+fcWG5uYMoIiLC6sYaOnuGcvJ4pL5l551sOeGu27DrajwKRUID25ig3rG5kXWuU+mobzylVmSrByYEx1rfmphycC9mE+05gnYisxgv19wLvn3qCiFwN/Atwu6p257yVszCQ43AHaG+pYf8pmzFjsjOeSLG9o49H95/hsX3dk73zy5fWcqtfalm2yHrnrsjMzOsqtXBX1aSIfAR4BG8q5NdUda+IfAbYpaoPAn8HRIHv+YM9x1X1zjy2+6Iyo9a53PmlvTnKI3vPEE+mqQjZbdrmQp39o5MDods6ehlPpCd3I7pxTSPrF8eotXnnTsp86iq2+2GyKhCp6kPAQ+cd+/SUx2/OcbvmbNCfb5rLj7lrmmtIpZXjZ0dYm6NZOKa0pdLKc8f7eXR/N99/tnPyrtCGmgquXr6IDYtjrG6qcX43IgOxyhABKdFwLyWZnnuubmICr+cO0NFj4V7OMlMV/3PfGR4/0MPZkTihgLCioZo7LvcC3dY7Lz8BEeqqwpwsss2y3Q33HNbcV/sbddigavk53jfKo/vP8Oi+brZ39JFSpSocZMPiGG/d1Mq6lhhVFXYjUbmrq6qY3K2qWLgX7mMJRMjpuhqxyjCttRGbDlkGMuWW/9zXzaP7znCw2/uF3t5cw81rG9m4uNamKpoL1FcX3xpUzoX74Gic2spwzv/xtTdHrefuqOGJJE+83POqcktAYFVTDW+/oo2Ni2M0+ncqGzOd+qowe7oGSaW1aH7xOxfuA2OJnJZkMtY01/Dg8ydRVaupOuDkwBiP7jvDN7cf43DvCKm0V25Z3xq1couZtbrqMMm00j00TltdVaGbA7gY7jlcEXKq9uYo58aT9I3EabJeXMlRVfafHuJne8/w832n2dN1DoDGmgpuWtPIxrYYKxtqiqbXZUpL5o74kwNjFu754q3lnvu7/dZkZsx0D1u4l4h0Wnm+c4CH95zme7tO0D+aeGUhrtcs5rI2m91iciNTLegaGOfalQVujM+5cB8cjbOyoTrnr9ve7M2YOdw7wg1rGnP++iY3VJUXOgf50S9P8tCLpzg5OE44KKxuquHWDS1stE0sTB5kqgXFNKjqXLjnq+a+pK6KynCAjm4bVC02qt7a5z958RTf29XJ2ZE4QRHWtUa5ZW0TGxfXWv3c5FUkHPTnulu450U6rQyO5afmHggIq5tsxkyxSKeV50708/Ce0/x0z2k6+8cIBbwe+m0bmtnUVmeBbhbUkvoqC/d8GRpPokpeau4Aly+p5WcvnSmq6U7lJJlKs+PIWX665xSP7D1Dz9AEQRHWtkT5jWuWcllbLdUVTv2VNiVkaX0lnf0W7nmR2ccwHz13gFs3tPC93Z08f6Kfa1faxh0LIZFKs72jj4dePMUje0/TP5qgKhzkDeubqa0Ks3FxjErbas4UgSX1VTxz5GyhmzHJrXDPw9IDAPfsOA7AWDxFQOCLjx3i6797fU6vYV6hqjx7vJ8fPNfFT144Rf9ogopQgMsWx7jj8jrWt8ZsdU5TdJbUV3FuPMnQeKIoBu3dCvc8rOU+VVVFkBUNNbx82tZ2z4djfSM88GwXP3iukxNnx6gMB3jrpsXUVoZZ1xq1FRZNUVta781v7xoYY+NiC/ecGhj1yjJ1OdpibzobFsd4ZO9pzpwbp7W2Mm/XKRdD4wkeevEUD+zu4pmjZxGgvSXKb167jE1ttUSs5GJKRGaBwaO9I2xcXFvg1jgW7oN57rkDbGj1wv0XB3p493XLZ/4Bc4F0Wnn6cB/f293JT/ecYjyRZk1TDW/d1MrVKxbldKMVYxbKK6vHFscCg06Fez52YTpfa22EuqowWw50W7jP0qnBMe7f1cl9u09Mll2uXFrPNSsXsXxRld0pakpaTSTEkrrKorkXxrlwj0ZCea3NigjrW2NsPdhLIpW2OvAMEqk0j+3v5t5njvOLl3tIK9yytpGb1jTxmiW19v/POKW9pXjuhXEr3MfiC/KRfkNrjJ1Hz7LraD83tdtSBNM50jvCd3ee4N+fPsbwRJLayhBvWN/MtSsbaKjJ35iIMYXU3hzl/t2dRbF6rFPhPjian6UHztfeUkM4KDx+oNvCfYqJZIqH95zmnh3H2XHkLMGA9yln88pFrG+N2Y1fxnntzTUMTyTpHpoo+IQLp8J9YCyxID33SCjIDasbeXjvaf7y9o1lH1pHe0f41A/38OzxfkbjKRpqKnjrplauWbmI2iKY72vMQmmfsnqshXsODYzG2bB4YTawfv8NK/ijbz/Lg7/s4tevXrYg1ywmqbSyZX8333z6GE+83ENA4LK2Wq5f3UB7c5SADY6aMtTe4od7zzA3r20qaFucCvfBsURe57hPdftrFrOprZYv/OdB3nHlkrIZGBwcS/DJB15g++E++kcT1FaGeNNlLVy3soFam8JoylxLLEI0EiqK6ZDOhLuqerswLUDNHbxVIv/8Lev58Dd38f1nO3nPdSsW5LqF0tEzzNefOsr9uzsZS6RY1VjD7Ze3samttuzLUsZkiAjtzTVFMWPGmXAfiadIpjVvi4ZN502XtfDa5fX846OHeOfVS4mE3LqbUlV58lAvX3vyCFsO9FARDHDnVUtYXFvJkvri2ErMmGKzpjnKjsN9hW6GO+GeWXpgoXru4P2W/ou3rueDX32G+3ae4IM3rVqwa+fTeCLFD57r4u9//jLdQxNEI17p5fpVDUWxIJIxxay9uYYfPNfFyESSmkjhItahcM/cnbowNffMSpGqyqrGaj778AGGxpP80W1rF+T6+XBqcIxvbj/Gd545zsBogra6St517TKuXFpHqEzGFIyZr8yMmSO9I1y+tK5g7XAm3BdiXZnpiAh3XbWUf/5FB996+hi/c/Oqgv62ni1V5bM/3c+2jj72nhxEFTYtqeVd1y5jdWNNwW/EMKbUTJ0xY+GeA0d6vdHptrqFn1vaWlvJ+65fwTe2HeVj9z7Pv3zw2qIfZBxPpPjxC6f4+rYj7Ok6R2U4wC3tTdy4ppFFdgepMXO2srGagFDwNWacCfftHX201VWyoqG6INdf3xrjHVe28aMXTvG/f7KP//GOy4qy13t6cJxv7zjGPTuO0zcSZ11LlLuuWsLVyxfZBhjG5EAkFGRFQ3XBp0M6Ee7ptLKto5c3bmwtaKDe1N5EYzTC1546wqnBMf72XVcWxQCkqrK9o49vPX2MR/aeRhU2Lo5x11VLaW+20osxudbeXPgFxJwI9/2nh+gfTXBzEazz8te/tokl9ZV87uED7P/SU/zTb13DZW2FWbj/5MAY33+2k/t3d3K0b5T66jC3rG3ihtWNtniXMXnU3hJl66FeUmktWInWiXDf1tELwM1rCx/u33nmBNFImN+7ZTX37jzOO774JNeuXMTn3/1ali3Kf8moe2icR/ac5t+eOsqR3hEUbxOBd127jCuW1pXNnbTGFFJ7cw3xZJqu/jFWNBamVOxIuPexpqmGtrriubFmdVPRhcx+AAAJSElEQVQNH7ltLY/u62b30X5u/bvHeefVS3n7FW3cuKaRqorc3PCUSKV5/sQAWw/28uTBHp47MYAqNEUj3LaxhWtWLLJeujELLPNp/amOXlY0Fubu9ZIP90QqzY7Dfbzz6qWFbsoFYpVh3nn1Um7d0MyZc+Pct8srkURCAa5f3cDGxTHWNEdZ3VRDU7SCuqoK6qrCrxrYTKeVsUSKEX8Z0ZMDY5waHOfAmSH2dg2y9+Q5kmlFgGWLqnjjhhYuX1pHSyxitXRjCuSKpXVc1lbL1586ynuvW16Qf4slH+4vdg0yEk9xc3thV2C7lPrqCuqrK/jEHRs52jvCgTNDHDwzzI4jZ4kn09P+TCggBAJy0ecrwwGW1Fdx05pGljdU094czdmnAWPM/IgIv3vzKv7SX2SvEPlU8uG+7ZBXby+FTTPCwQDrWmOsa/WWJU6rMjiaoG8kzkg8yWg8xVg8RSqtqCppVULBAJFQgHAwQDQSor46TF1VmGgkZD1zY4rYnVct4bMP7+ffnjpq4T4X2zr6uKyttiTrygERFtVU2E1DxjioMhzk/dev4MuPH+LE2VGWL/A9OFlNnRCR20XkgIgcEpFPTPN8RES+6z+/Q0RW5bqh0xlPpNh1rJ9bSqDXbowpPx+4cSVBEb6x7eiCX3vGcBeRIPBl4A5gE/A+Edl03mkfAvpVdS3w98Dnct3Q86kq9+/uJJ5MF8UUSGOMOd/iukruuKKN7+46Mbn+1ULJpud+PXBIVQ+rahy4F7jrvHPuAr7hP74feJPkqSCsqmw92MM7/2kbn/rhHjYujnHjGgt3Y0xx+r1bVjE8keR1n3uMT//HHvafPrcg180m3JcCJ6Z83+kfm/YcVU0Cg0BeEvdLjx3ig199hp5z4/ztb1zJjz/6OqorSn7owBjjqKtXLOL+P7iZN1/Wyr07T3D7F7bytSeP5P26C5qKInI3cLf/7bCIHJjrax0Dtv/3rE5tAnrnep0S4Pr7A3uPrnD6Pf7WLM790Oe8WvYcrczmpGzCvQtYPuX7Zf6x6c7pFJEQUAdcsM+Uqn4F+Eo2DcsVEdmlqpsX8poLyfX3B/YeXVEO77GYZFOW2QmsE5HVIlIBvBd48LxzHgR+x3/8LuAxVdXcNdMYY8xszNhzV9WkiHwEeAQIAl9T1b0i8hlgl6o+CHwV+JaIHALO4v0CMMYYUyBZ1dxV9SHgofOOfXrK43HgN3PbtJxZ0DJQAbj+/sDeoyvK4T0WDbHqiTHGuMcW9zbGGAc5G+4zLZlQ6kTkayLSLSJ7Ct2WfBGR5SKyRUReEpG9IvKxQrcpl0SkUkSeEZFf+u/vbwrdpnwRkaCIPCciPy50W8qFk+Ge5ZIJpe7rwO2FbkSeJYGPq+om4Ebgjx37c5wA3qiqrwWuAm4XkRsL3KZ8+Riwr9CNKCdOhjvZLZlQ0lT1CbyZSc5S1VOq+qz/eAgvHIpvV5Y5Uk9mF+Ww/+XcIJiILAPeDvxrodtSTlwN92yWTDAlxF9p9GpgR2Fbklt+ueJ5oBv4uao69f58XwD+Eph+5xmTF66Gu3GIiESBB4A/VdWFWXVpgahqSlWvwrvz+3oRubzQbcolEXkH0K2quwvdlnLjarhns2SCKQEiEsYL9m+r6vcL3Z58UdUBYAvujaPcAtwpIkfxyqNvFJF/L2yTyoOr4Z7NkgmmyPnLRn8V2Keqny90e3JNRJpFpN5/XAW8Bdhf2Fbllqp+UlWXqeoqvH+Hj6nqBwrcrLLgZLj7yw5nlkzYB9ynqnsL26rcEpHvANuBDSLSKSLzWGSuaN0CfBCvt/e8//W2Qjcqh9qALSLyAl6H5OeqalMFTU7YHarGGOMgJ3vuxhhT7izcjTHGQRbuxhjjIAt3Y4xxkIW7McY4yMLdGGMcZOFuFoyIpPy56nv9ZW4/LiIB/7nNIvKPl/jZVSLy/oVr7QXXHvPXgCkKIvIefzlrmxdvpmXhbhbSmKpepaqvwbsb8w7grwFUdZeq/sklfnYVUJBw93X4a8BkzV96Oi9U9bvAh/P1+qb0WbibglDVbuBu4CPiuTXTCxWRN0y5I/U5EYkBnwVe7x/7M783vVVEnvW/bvZ/9lYReVxE7heR/SLybX8ZA0TkOhHZ5n9qeEZEYv6qjH8nIjtF5AUR+a/ZtF9Efigiu/1PIXdPOT4sIv9XRH4J3HSRa77Gf/y8f811/s9+YMrxf8n8cvA3nnnWf41Hc/jHYFymqvZlXwvyBQxPc2wAaAVuBX7sH/sRcIv/OIq3kfvk8/7xaqDSf7wO2OU/vhUYxFssLoC3RMPrgArgMHCdf16t/7p3A5/yj0WAXcDq89q4Cthz3rEG/79VwB6g0f9egXf7jy92zS8CvzXlnCrgMv99h/3j/wT8NtCMt3z16qnXnfJefzzd/2v7sq/QLH8XGLMQngI+LyLfBr6vqp1+53uqMPAlEbkKSAHrpzz3jKp2Avh18lV4gX9KVXcCqL90sIi8FbhSRN7l/2wd3i+LIzO08U9E5Nf9x8v9n+nz2/KAf3zDRa65HfgrfxOL76vqQRF5E3AtsNN/r1V4a7zfCDyhqkf813B6gxaTOxbupmBEZA1eGHbj9VwBUNXPishPgLcBT4nIr07z438GnAFei9dDH5/y3MSUxyku/fdcgI+q6iOzaPetwJuBm1R1VEQeByr9p8dVNXWpn1fVe0RkB97uRA/5pSABvqGqnzzvWr+WbbuMmcpq7qYgRKQZ+GfgS6qq5z3Xrqovqurn8FZL3AgMAbEpp9Xh9YrTeCtHzjR4eQBoE5Hr/GvERCSEt3LoH/rrxiMi60WkZobXqgP6/WDfiNe7zvqa/i+1w6r6j8B/AFcCjwLvEpEW/9wGEVkJPA38ioiszhyfoW3GANZzNwuryi+ThPE2v/4WMN067X8qIrfhbcu2F/ip/zjlD1R+Ha8m/YCI/DbwMDByqQuralxE3gN80V87fQyv9/2veGWbZ/2B1x7gnTO8j4eBPxCRfXgB/vQsr/lu4IMikgBOA/9HVc+KyKeAn/nTQxPAH6vq0/6A7ff94914M42MuSRb8teYGYi3f+uPVbWotsDzy0N/oarvKHRbTPGxsowxM0sBdcV2ExPep5f+QrfFFCfruRtjjIOs526MMQ6ycDfGGAdZuBtjjIMs3I0xxkEW7sYY46D/D6ieVSmimgE/AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nb_merge_dist_plot(\n", " SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n", " SkyCoord(hsc['hsc_ra'], hsc['hsc_dec'])\n", ")" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Given the graph above, we use 0.8 arc-second radius\n", "master_catalogue = merge_catalogues(master_catalogue, hsc, \"hsc_ra\", \"hsc_dec\", radius=0.8*u.arcsec)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(np.sum(master_catalogue['flag_merged'])/len(master_catalogue)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add DECam (DES and DECaLS)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmUXGd55/HvU1Vd1V3V+2ZJ3Wq1JMuLjBeMvGETm90QYmdOGMAEmElMTBwgGxNCBg6ZkJwZGGY4CVuIA4RAsIGwZBwQcRJiY4xlW7JkS5aErL3Vkmz1vm/V9c4fdatdlqXu6u5bXVW3fp9z+rjq1u26T1nSr99+7nvfa845REQkWEKFLkBERPyncBcRCSCFu4hIACncRUQCSOEuIhJACncRkQBSuIuIBJDCXUQkgBTuIiIBFCnUgZubm11nZ2ehDi8iUpKefPLJXudcy0L7FSzcOzs72bFjR6EOLyJSkszseC77qS0jIhJACncRkQBSuIuIBJDCXUQkgBTuIiIBpHAXEQkghbuISAAp3EVEAkjhLiISQAteoWpmXwXeApxxzr3sHK//OvDHgAEjwN3Ouaf9LtQP9z7e9ZJt77yuowCViIjkVy4j968Bt87z+lHgZufc5cCfA/f4UJeIiCzDgiN359zDZtY5z+uPZj19DGhfflkiIrIcfvfc7wR+7PN7iojIIvm2KqSZvZp0uN80zz53AXcBdHSo1y0iki++jNzN7Argy8Dtzrm+8+3nnLvHObfFObelpWXB5YhFRGSJlh3uZtYBfB94t3Pu2eWXJCIiy5XLVMj7gFuAZjPrBv4UqABwzn0J+DjQBHzRzACSzrkt+SpYREQWlstsmTsWeP29wHt9q0hERJZNV6iKiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgAKdxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCSOEuIhJACncRkQBSuIuIBJDCXUQkgBa8zV7Q3ft41zm3v/O6jhWuRETEPxq5i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCaMFwN7OvmtkZM3vmPK+bmX3WzA6Z2W4zu9r/MkVEZDFyGbl/Dbh1ntffBGzyvu4C/nr5ZYmIyHIsGO7OuYeB/nl2uR34ukt7DKg3s9V+FZhPT50Y5ET/eKHLEBHxnR899zbgRNbzbm9b0fvh7lP8aM/pQpchIuK7FT2hamZ3mdkOM9vR09Ozkod+iWQqxfj0LF394wyOTxe0FhERv/kR7ieBtVnP271tL+Gcu8c5t8U5t6WlpcWHQy/d2NTs3OM9J4cKWImIiP/8CPf7gfd4s2auB4acc0Xf6xidTAJgwO5uhbuIBMuCq0Ka2X3ALUCzmXUDfwpUADjnvgRsBd4MHALGgd/IV7F+GpmaAWDzmlr2nhqmb3SKpupYgasSEfHHguHunLtjgdcd8H7fKlohmZH7DRub2HtqmD0nh7jl4tYCVyUi4o+yvUJ1dCod7msb4qxtqFLfXUQCpWzDfWQqSSwSoiIc4or2ek4PTdIzMlXoskREfFG24T46maQ6lu5KvaytDgP2nBwsbFEiIj4p33CfSlJTmQ73uqoKVtVV0qWrVUUkIMo33LNG7gA1lZG5PryISKkr23AfmZqhuvKFcE9EIy+6sElEpJSVZbgnZ1NMzqSojlXMbUvEIoxNJUnP7BQRKW1lGe6Z9ktNVlumOhYhmXJMz6YKVZaIiG/KOtxf1JaJhQHUmhGRQCjLcB/xrk7NPqGaiKYfj+mkqogEQFmG+1xb5kUjdy/cpxXuIlL6yjLcMyP3ROwc4a6Ru4gEQFmG++hUksqK9NIDGYmoeu4iEhxlG+7Z0yABopEQkZBp5C4igVCe4T4586J+O4CZpee6q+cuIgFQnuE+9eKlBzISsbDaMiISCGUZ7iOTyRfNcc9IRDVyF5FgKLtwn5lNMZVMvejq1IxETIuHiUgwlF24j57jAqaMRDSsE6oiEgjlF+7nWHogozoWYWbWMZ3U+jIiUtrKLtwzFzDVnDUVEnSVqogER9mF+3wjd12lKiJBUXbhPjI1A7ywCmQ2XaUqIkFRduE+OpmkqiJMJPTSj66Ru4gERfmF+9S557iDeu4iEhw5hbuZ3WpmB8zskJl95Byvd5jZg2a2y8x2m9mb/S/VH6OTyXPOcQeIRUKEtb6MiATAguFuZmHgC8CbgM3AHWa2+azdPgZ8xzn3cuAdwBf9LtQv843czcyb666eu4iUtlxG7tcCh5xzR5xz08C3gNvP2scBtd7jOuCUfyX663zrymRo8TARCYLzp9wL2oATWc+7gevO2ud/AP9qZh8EEsDrfKnOZynnmEqmqKp46UyZjEQsoraMiJQ8v06o3gF8zTnXDrwZ+IaZveS9zewuM9thZjt6enp8OnTuMleexuYL92hY68uISMnLJdxPAmuznrd727LdCXwHwDm3DagEms9+I+fcPc65Lc65LS0tLUureBkmZ9K99Fjk/B873ZZRz11ESlsu4b4d2GRm680sSvqE6f1n7dMFvBbAzC4lHe4rPzRfwFRm5L5AuE8nU3M/CEREStGC4e6cSwIfAB4A9pOeFbPXzD5hZrd5u30I+C0zexq4D/ivzjmXr6KXKhPulfO0Zaqj6dMQ/WPTK1KTiEg+5HJCFefcVmDrWds+nvV4H3Cjv6X5byqntkw6+PvHpllTX7UidYmI+K2srlCdnGvLzD9bBqBPI3cRKWFlFe7TSW/kXjHPyN1ry/SNTq1ITSIi+VBW4T45k9sJVVDPXURKW1mF+1Rm5D5PW6ayIkTI1JYRkdJWXuE+k6IibIRDdt59zIxELEL/qMJdREpXeYV7MkV0nlF7RiIa0chdREpaWYX7ZHKWynn67RmJWJj+MZ1QFZHSVVbhPjWTmnemTEY8GmFwfGYFKhIRyY/yCvdkat6TqRnxaJj+cbVlRKR0lVm4z847DTIjHo0wNDHDbKroVlAQEclJmYV7at51ZTLi0TDOwfCEWjMiUprKK9xnchu5z60vo9aMiJSo8gr3ZCrntgzAoMJdREpU2YT7dDJFMuXmvQtTRjya3mdgTG0ZESlNZRPumfuiLmbkrraMiJSqsgn30blwz33krraMiJSqsgn3kcncR+6xSIiKsNGvtoyIlKiyCfex6XS45zIV0syoj0c1cheRklU24T66iJE7QEO8ggGFu4iUqLIJ95FFnFAFaIhHNVtGREpW2YT73GyZHNoy4IW7Ru4iUqLKJtwX3ZZJVDCglSFFpESVTbhn2jLRRbRlBsencU6Lh4lI6SmbcB+dTBKLhAjZ+W+xl60hHiWZcnM/FERESknZhPvYVDLnlgxAfbwCgAHdbk9ESlBOaWdmt5rZATM7ZGYfOc8+bzOzfWa218zu9bfM5RudSuZ0dWpGYyIKoL67iJSkyEI7mFkY+ALweqAb2G5m9zvn9mXtswn4E+BG59yAmbXmq+ClGplK5nSLvYz6eCbcNXIXkdKTS9pdCxxyzh1xzk0D3wJuP2uf3wK+4JwbAHDOnfG3zOUbm0pSuZSRu9oyIlKCcgn3NuBE1vNub1u2i4CLzOznZvaYmd16rjcys7vMbIeZ7ejp6VlaxUs0OpnMeaYMpK9QBbVlRKQ0+XVCNQJsAm4B7gD+1szqz97JOXePc26Lc25LS0uLT4fOzehUkspFtGVqKysImVaGFJHSlEvanQTWZj1v97Zl6wbud87NOOeOAs+SDvuisdgTqqFQevGwfrVlRKQE5RLu24FNZrbezKLAO4D7z9rnn0iP2jGzZtJtmiM+1rkszjkv3Bf3i0p9vIJBtWVEpAQtmHbOuSTwAeABYD/wHefcXjP7hJnd5u32ANBnZvuAB4E/cs715avoxZqcSTGb4y32sjVo5C4iJWrBqZAAzrmtwNaztn0867ED/tD7Kjqji1wRMqMhHqV7YDwfJYmI5FVZXKG69HBXW0ZESlN5hPtk7ndhytaQiNKvxcNEpASVR7gvoy0znUwxMTObj7JERPKmzMJ9sSdUdSGTiJSmMgn3dDgvZm0ZSLdlQEsQiEjpKZNwT7dVltKWAS0eJiKlpzzCfaknVNWWEZESVR7hPjVDOGREQrndhSlDbRkRKVXlEe6TSapjESzHW+xl1FdlRu4KdxEpLeUR7lOzVMdyuhj3RSLhEDWVEY3cRaTklEm4zywp3CF90w713EWk1JRJuCeprlxauNfHo2rLiEjJKZNwX1pbBtIzZhTuIlJqyiPcJ5felmlKxOgbVbiLSGkpj3CfSi453C+ojdEzMkUqpcXDRKR0lEW4j03NLrnn3loTI5ly9Ks1IyIlJPDhnkqlb7GXWOLIvbW2EoAzw1N+liUikleBD/fR6fTSAzXLaMsAPD8y6VtNIiL5FvhwH/LmqNd5V5suVmtNeuTeo5G7iJSQwIf78GQ63GuXGO4tNemR+xmN3EWkhAQ+3Icmljdyr6wIU1dVwfMauYtICVlaI7qEDC8x3O99vGvucSwS4snjA77WJSKST4EfuQ9PpE+o1lYt/edYbWUFI5NaX0ZESkfgw325bRmAmsoII94NP0RESkFZhHvIWPIVqgA1lRWMTCZxTlepikhpyCnczexWMztgZofM7CPz7PdrZubMbIt/JS7P8OQMtVUVi75RR7aaygizzmnpXxEpGQuGu5mFgS8AbwI2A3eY2eZz7FcD/B7wuN9FLsfQxMyyWjKQDnfQdEgRKR25jNyvBQ45544456aBbwG3n2O/Pwc+BRRVAvoR7rWV6e/XEgQiUipyCfc24ETW825v2xwzuxpY65z70XxvZGZ3mdkOM9vR09Oz6GKXYnhiZi6clyozcn9+uKh+bomInNeyT6iaWQj4DPChhfZ1zt3jnNvinNvS0tKy3EPnxJ+2jDdyH9HIXURKQy7hfhJYm/W83duWUQO8DHjIzI4B1wP3F8tJ1aGJ5JKXHsiIRkJUVoToUbiLSInIJdy3A5vMbL2ZRYF3APdnXnTODTnnmp1znc65TuAx4Dbn3I68VLxI6dkyy78QtyZWobaMiJSMBcPdOZcEPgA8AOwHvuOc22tmnzCz2/Jd4HJMzswynUwtuy0DUFMVUVtGREpGTkNa59xWYOtZ2z5+nn1vWX5Z/vDj6tSM2kqN3EWkdAT6CtVMuC93tgykb/ZxZmRKV6mKSEkIdLgvdUXIc6mpqmA6mZpbiExEpJgFOtz9bMvMzXXXVaoiUgLKItyXOxUSdJWqiJSWQIe7r20ZrS8jIiUk0OE+lLlRR6UP89znliDQyF1Eil/Aw32GRDRMJLz8jxmLhKmORTRyF5GSEOhwH55c/roy2VprYuq5i0hJCHS4D03M+HIyNaO1NqaRu4iUBIX7IrTWVKrnLiIlYflnGovY8MQMaxvjvr3fmvoqfvzMaZKzKV/6+CJSOu59vOsl2955XUcBKslN4MPdz577xpYEM7OOEwMTrG9O+Pa+IrLyhiZm2H96mH2nhjnaO8apwQlODk4wMD5NZ1OC9c0JNrZUk4iVZkyWZtU58uNGHdk2tFQDcKRnVOEuUkL6x6bZc3KIZ7yvPSeH6B6YmHu9qiJMfbyCuqoKaisr2HVikMeP9hOLhHj/LRfSXBMrYPVLE9hwT86mGJue9WXRsIyNLelAP9Izxmsv9e1tRcRHA16QZ8J8d/cQJwdfCPLGRJS2+iouW13L6voqVtdVzt1tLWM25egeGOfr245z3/Yu7r55Y8m1YgMb7sOT6QuY6ny4UUdGfTxKUyLK4Z5R395TRJbuzMgk+04Ns/fU8DlH5E2JKGvqq7i8rY62hirW1FVRFQ0v+L7hkLGuKcFbX9HONx47zo/3PsevXLEmnx/Fd4EN97lFw+L+jdwBNrQkONIz5ut7isj8ppMpDveMsv/0ML94boT9p4fZ1TXI6NQLq7Q2ekF+2Zo62uqraKvPLcjnc+nqWl65sYlHD/dxYUs1l66uXe5HWTGBDfdhH9dyz7ahuZqf/OJ5X99TRF4wPDnDfm80vs874XnguRFmvXsphEPGBTUxLrqgmtV1Vayur2R17fKD/HxuvWwVx3rH+O6T3Xz4jRcTq8jPcfwW2HD3c7nfbBtbE3x7xzRD4zO+/1YgUm7OPtG599QwXf3jc683V8fYvKaW5uoYq+sqWVVXSXN1jHDIVqzGSDjEm69YzZd/dpRnz4xyeVvdih17ORTui7ShOT1j5nDvKFd3NPj63iJBNjqVZE/3ELu7B/nn3afpHhhncHxm7vXGRJQ1dZW8YfMFrK6rYk39S090FkpnU4J4NMy+U0MK90IbnvRvLfdsG7JmzCjcRc5tcmaWfaeH2dM9xNPdg+zuHuJwzyiZu1Q2xCtY2xDn+vVVizrRWSghMy5dVcve00MkUykioeKfORPYcM/XyH1tY5yKsGnGjIhnZjbFgedG2N09xJ6Tgzx9Yohnnx8hmUoneXUsQntDFa+5pJX2+jhtDVVUl+CFQZvX1PJk1wBHe8fY1FpT6HIWVHr/h3M0NDFDNBKi0ueTHxXhEB2NcY4o3KUMTc7M8uzzI3NTD585Ncz+08NMJ1NA+mKg9oYqbrqwmbaGKtob4tRWRjBbuR55vlzYWk1F2Nh3aljhXkjDE0nfZ8pkbGip1nRICbzx6ST7Tg3PXRC079QwB8+MMuuNyCsrQqypq+LazkbaG9JTDxsT0UAE+blUhENsaq1h/+lhbrtyTdF/zgCH+4yvFzBl29hSzUMHzmgBMQmMTGvl6e5Bnj6R7pE/+/wIXo5THYuwpr6SV13YzOr6KtbUVQY6yM9n8+pa9p0e5uTgBO0N/i1KmA85pZ+Z3Qr8FRAGvuyc++RZr/8h8F4gCfQAv+mcO+5zrYvi97oy2TZ4C4h1D0zQqTVmpMQ45zjRP8FTXpA/fWKQZ04NMTmTbq3Eo+nWyi0Xt85dDOT3xIRSdcmqGgzYd3q49MPdzMLAF4DXA93AdjO73zm3L2u3XcAW59y4md0N/G/g7fkoOFdDEzM0VUfz8t5za8z0jircpej1jk7NhfhT3UPsONbP+PQsAJGQsaa+ild0NNDeEKe9IditleWKxyJ0NifYd2qYN2xeVehy5pXLyP1a4JBz7giAmX0LuB2YC3fn3INZ+z8GvMvPIpdieHJmbtqi3+bmup8Z4zWX5OUQIksyOTPL3lPD7Ooa4KkTgzxyqHduLnnI0jec2by6lraGKtY2xLmgtnJFLwgKgs2ra/nRntP0j00XupR55RLubcCJrOfdwHXz7H8n8OPlFOWHfLZlGhJRGhNRjvRqxowUTvbMlT0n0xcH/eL0C1MQ2+rTs1Vu2JD+b1t9FdGIzhEtV2bQeLyvuCdV+HrG0czeBWwBbj7P63cBdwF0dOTvDiaplGN4YiZvs2UANjQnOHymuP9wJTgmpmfZeyqzjO0we08NcejM6FyQV1aEaK+Pc+OFzbR7o3L1yfPjgtpKYpHQi5ZJKEa5hPtJYG3W83Zv24uY2euAjwI3O+fOeaNR59w9wD0AW7ZscYuuNkdj00lSzv8LmLJtbNECYpIfzjmO9o6xs2uQnV0DPNU1yIHnR+amICZiEdrqK7nxwmbWeCc8G+IV6pOvkJAZ7Q1VgQj37cAmM1tPOtTfAbwzewczeznwN8Ctzrkzvle5SJleWH0eF/bavKaWb+84wYn+cV/v0yrlZ2wqydPdg+zqGmTn8QF2nRic+zsci4RY2xDnVZuaWeu1VmoCclFQKetojPPQgR7GppJFexu+BatyziXN7APAA6SnQn7VObfXzD4B7HDO3Q98GqgG/tH7S9flnLstj3XP60R/erH+fIbuDRubANh2uE/hLjnLHpU/eXyAXV0DHHhuhMyvsS3VMdY3Jbh5UwsdTXFaamKEFORFp6MxjgOe7h7klRubC13OOeX0I8c5txXYeta2j2c9fp3PdS1L5telDp9DN/vu5845ErEI9z7RxduuWTvPd0k5m0rO8szJIbYfG2DHsQEePdw7Nw2xsiI9Kn/1Ja2sbahibWOceLQ4R4HyYpkB3a6uEg/3UnO8f4xoOMQFtZV5O4aZsbElwZGeUZxz+jVZABiZnGFn1yDbj/bzxLF+nj4xyJS37sr65gSXrKplXVOcjkaNyktZPBqhuTrGzuMDhS7lvAIZ7if6x2lvrMr7/N2NzdXs7h7iSO8YG1uq83osKU7PDU2y/Vg/O471s/3YAL94bpiUS88pX11XxZZ1DaxrStDZnCjJlRDl/Doa4+w6MVi0g7tA/m3r6h/3vSVzLpn5ro8e7lO4lwHnHEd6x3jiaP/cyDxzM+ZoOER7Y/qS/c6mBGsbq4hFind9clm+jsY4O7sGONY3zvoivFI9cOHunON43/iK3EijMRGlvqqCbYd7eff16/J+PFlZydkU+0+PsP1YP9uP9fPwwV7GvBsyJ6JhOpsTXNlez7qmOKvr8v+bohSXzABy5/EBhftKGJqYYWQyuSIjdzNjQ0uCbYf7SKUcIf3jLmlDEzPs6hpg5/EBnuwaYFfX4NzJz7b6Ki5qraazOUFnU4Lmaq2/Uu5aa2NUxyLs7Brg117RXuhyXiJw4Z6vmTLns6Glmp3eRSaXrq5dkWPK8mWmJO447oX58QEOnkkvJxEOGRfUxriivZ7OpjjrmhJ5vSBOSlPIjKvW1rOza7DQpZxTcMO9aYXCvfmFvrvCvXhNJ1M8cyq9IuITRwfY2TUwd6FQVUWYjsY4r7v0AtY1pVdGVL9ccnF1Rz2ff/AQo1PJojthXlzV+OB4Xzrc167QWsv18SjrmxNsO9zLnTetX5FjysLGp5Ps6hrkce/k547j/czMpi8VakpE6WxKcPNFLaxrjNOsKYmyRFevayDl4KmuQW7aVFzz3QMX7if6x2mujq7oJcE3bGzin586xcxsigrdmakghsZn5k58bt1zmpODE6QcGLC6rpJrOhvTUxKb4tTkcUE5KS9bOhsJh4xtR3oV7vm2UtMgs71+8wXc+3gXW/ec5var2lb02OWqb3SKJ4728/jRfh470seB50dwDirCxpq6Kl61qYXOpgTrmuK+3yRdJKM6FuGK9jq2He4rdCkvEbhwP943zjWd+Z8Gme3mTS1saEnwtz87UhI3zi1FQxMzPH6kj0cP97HtcDrMIR3m6xoTvPaSC+hsjrO2Ia7fnmRFvXJjE1/66ZGi67sXTyU+mE6mOD00QUfjyo6eQyHjzpvW89EfPMMTR/u5bkPTih4/iKaTKXZ2DfDIwV4eOdTL0ycGcXhh3pTgDZsvYENzgraGuOaXS0HdsKGZLzx4mO3H+nn1xa2FLmdOoML9lNdnLcQqjb92dTv/54EDfPmRowr3Jchc/fnwsz387GAvjx3pY3x6lnDIuLK9jlsubmVja4KOhjgRjcyliLxiXQPRcIhth/sU7vmSmQa5rmnlrxarrAjz7uvX8bkHD3G0d6wor1grNgNj0/z8cC+PHOzlZwd7OTmYvpS/MRHl8rY6NrVWs6GlWj1zKWpV0TBXddQXXd89UOF+fIUvYDrbu25Yx5d+eoS/+/lRPnH7ywpSQzEbn07yxNF+th3u4+eHe9l7chhHeunbDc3VbOlsYFNrDY2JaKFLFVmUV25s4q9+cpCh8Rnq8niToMUIVLif6B8nGgnRWhMryPFbayq5/ao1/OOObu6+ZSOr66oKUkexmJlN8dSJQR452Mu2w33sOjHAzKwjGg7x8o56XnNpK5taa2ir17osUtpeubGZv/z3gzx+tI83XLaq0OUAAQv3rr70NMiVXuMl+yYeHY1xZp3j7n/Yybffd31ZXenonONY37jXN+/hsSP9jE4lMWBNfRU3bGhmY0uCdU0JohH1zSU4rlxbR2VFiEcPK9zzohBz3M/WVB3jrVe3c+8TXfz5D/fxF796eUHrybexqSTbDvfx02d7+OmzPXPnPRoTUTavqeXClmo2tlRTFS2fH3JSfmKRMNd0NvLYkeLpuwcm3J1zdPWPc+36xkKXwsva6njfL23gbx4+wlVrG3hrEa4Yt1SplGP/c8M8/GwvPzvYw/Zj6cv649EwN2xo4qq19WxqraapujCtMZFCuX5DE59+4AB9o1NF8fc/MOE+MD7D6FSyaG5W/UdvvJjd3UN89Ad7iEVC/MqVawpd0pJk1sd/9HAfjx5O9877vAW3VtVWcv36Ji5aVcO6Rk1RlPJ204XNfPqBAzyw93neeV1HocsJTrg/cTT969CFrcVxR6RIOMTn3/lyfuvrO/jgfbt45GAvf3rb5qK/AXIq5Xj2zAg7jg3w+NF+fnrgDMOT6RtU1FZG2NBSzasvaeXC1mpqtUaLyJwr2uu4vK2Ov/3ZEd5+zdqCTxIo7qRZhK88cpS2+ipu3Fg8FxA1Vcf49vtu4C///Vm++NBhth/v5w9edxFvvGxV0ZxQ7BudYnf3EE+dGOSpE4Ps7BpgxAvz1poYnc0J1jcn2NhcTZNuUCFyXmbG79yykbu/uZMfP3Oat1xR2N/WAxHuu7sH2X5sgI/98qVF0xrInkHTVh/nN29czw92neSD9+2ipSbGHdd28PpLL+DS1TUrUvPE9CyHe0Y53DPKs8+PsP/0CPtODfPc8CSQvqFzS02MS1bVsq4xzrqmOI0JhbnIYrzhslVsaE7w1w8d5pcvX13Qfz+BCPevPHKU6liEt12zttClnNfGlmr+8PUX0VZfxd9vO8bn/uMgn/3JQRLRMFeva+DiC2roaIrT0RinpSZGfTx9f9Z4NDzvXxDnHKNTSQbHZxgcn6F3dIrnhid5bmiSU4MTdPWP09U/znPDk7j0cuaELD0nf1VdJS/vqKe9Ic6a+sqymrYpkg/hkPHbN2/kw9/bzcMHe7n5opaC1VLy4X56aIIf7T7Ne27oLPoecMiM00OTvGHzKq5f38TRvjGO9Y5x6Mwo24/1MzmTesn3mEEsEiIWCVMRToe8c+kLhCaTKaaTL/0eSK9jXl0ZoTERZVVtJZvX1NJaU0lrTYymRLRofsMRCZpffXkbn/m3Z/nrhw4Vf7ib2a3AXwFh4MvOuU+e9XoM+DrwCqAPeLtz7pi/pZ7b3z96nJRz/MaNnStxON/UVlVwZXs9V7bXA+kR+MhUkoGxaUankoweoV1RAAAIBklEQVRPzzI+Pct0MkVyNsVMKkXKvfD9YTMqwkYkHCIWCRGPholHIySiYWqrKqiprCj4CR2RchSNhHjvq9bzFz/az4/3nOZNl68uSB0LhruZhYEvAK8HuoHtZna/c25f1m53AgPOuQvN7B3Ap4C356PgbKeHJrjviS7eeNmqopkCuVRmRm1lRdH/9iEiC7vj2g6+t/Mkd39zJ3fetJ4P33rxirc9cxm5Xwsccs4dATCzbwG3A9nhfjvwP7zH3wU+b2bmnHP4zDnHzq4B/u7nx/iXZ54D4H03b/T7MCIiS5aIRfjB77yS/7V1P1955CiPHenjPTesY6N3xXbDCiyOl0u4twEnsp53A9edbx/nXNLMhoAmoNePIrP945PdfPi7u6mpjPAbN3by7us76Wgq7VG7iARPZUWYP7v9Zdy0qYWPfG83f/y9PXOv3fVLG/jvb740r8df0ROqZnYXcJf3dNTMDizn/T7mfS1CM3n4gVMkgvzZINifL8ifDQL8+X59id/30U/BR5d+2HW57JRLuJ8EsucYtnvbzrVPt5lFgDrSJ1ZfxDl3D3BPLoXlg5ntcM5tKdTx8ynInw2C/fmC/Nkg+J+vWOUyH247sMnM1ptZFHgHcP9Z+9wP/Bfv8VuB/8hHv11ERHKz4Mjd66F/AHiA9FTIrzrn9prZJ4Adzrn7ga8A3zCzQ0A/6R8AIiJSIDn13J1zW4GtZ237eNbjSeA/+1taXhSsJbQCgvzZINifL8ifDYL/+YqSqXsiIhI8ugZdRCSAyiLczexWMztgZofM7COFrsdPZvZVMztjZs8Uuha/mdlaM3vQzPaZ2V4z+71C1+QnM6s0syfM7Gnv8/1ZoWvym5mFzWyXmf2w0LWUm8CHe9byCW8CNgN3mNnmwlblq68Btxa6iDxJAh9yzm0GrgfeH7A/uyngNc65K4GrgFvN7PoC1+S33wP2F7qIchT4cCdr+QTn3DSQWT4hEJxzD5OeoRQ4zrnTzrmd3uMR0iHRVtiq/OPSRr2nFd5XYE6CmVk78MvAlwtdSzkqh3A/1/IJgQmIcmFmncDLgccLW4m/vLbFU8AZ4N+cc0H6fH8JfBg497rUklflEO5S4sysGvge8PvOueFC1+Mn59ysc+4q0ld+X2tmLyt0TX4ws7cAZ5xzTxa6lnJVDuGey/IJUqTMrIJ0sH/TOff9QteTL865QeBBgnP+5EbgNjM7RroV+hoz+4fCllReyiHcc1k+QYqQpe8v+BVgv3PuM4Wux29m1mJm9d7jKtL3TPhFYavyh3PuT5xz7c65TtL/5v7DOfeuApdVVgIf7s65JJBZPmE/8B3n3N7CVuUfM7sP2AZcbGbdZnZnoWvy0Y3Au0mP+p7yvt5c6KJ8tBp40Mx2kx6E/JtzTlMGxRe6QlVEJIACP3IXESlHCncRkQBSuIuIBJDCXUQkgBTuIiIBpHAXEQkghbusGDOb9eaq7/WWuf2QmYW817aY2Wfn+d5OM3vnylX7kmNPeGvAFAUze7u3hLXmxcs5KdxlJU04565yzl1G+mrMNwF/CuCc2+Gc+915vrcTKEi4ew57a8DkzFtuOi+cc98G3puv95fSp3CXgnDOnQHuAj5gabdkRqFmdnPWFam7zKwG+CTwKm/bH3ij6Z+Z2U7v65Xe995iZg+Z2XfN7Bdm9k1vGQPM7Boze9T7reEJM6vxVmX8tJltN7PdZva+XOo3s38ysye930Luyto+amb/18yeBm44zzEv8x4/5R1zk/e978ra/jeZHw7ezWZ2eu/xEx//GCTInHP60teKfAGj59g2CFwA3AL80Nv2z8CN3uNq0jdyn3vd2x4HKr3Hm4Ad3uNbgCHSC8SFSC/NcBMQBY4A13j71XrvexfwMW9bDNgBrD+rxk7gmbO2NXr/rQKeAZq85w54m/f4fMf8HPDrWftUAZd6n7vC2/5F4D1AC+klq9dnHzfrs/7wXP+v9aWvyCJ/FoishJ8DnzGzbwLfd851e4PvbBXA583sKmAWuCjrtSecc90AXp+8k3Tgn3bObQdw3tLBZvYG4Aoze6v3vXWkf1gcXaDG3zWz/+Q9Xut9T59Xy/e87Ref55jbgI96N7P4vnPuoJm9FngFsN37rFWk13i/HnjYOXfUe49A3phF/Kdwl4Ixsw2kw/AM6ZErAM65T5rZj4A3Az83szee49v/AHgeuJL0CH0y67WprMezzP/33IAPOuceWETdtwCvA25wzo2b2UNApffypHNudr7vd87da2aPk75L0VavFWTA3zvn/uSsY/1KrnWJZFPPXQrCzFqALwGfd865s17b6Jzb45z7FOnVEi8BRoCarN3qSI+KU6RXjlzo5OUBYLWZXeMdo8bMIqRXC73bWzceM7vIzBILvFcdMOAF+yWkR9c5H9P7oXbEOfdZ4P8BVwA/Ad5qZq3evo1mtg54DPglM1uf2b5AbSKARu6ysqq8NkkF6ZtffwM41zrtv29mryZ9e7a9wI+9x7Peicqvke5Jf8/M3gP8CzA234Gdc9Nm9nbgc97a6ROkR99fJt222emdeO0BfnWBz/EvwG+b2X7SAf7YIo/5NuDdZjYDPAf8T+dcv5l9DPhXb3roDPB+59xj3gnb73vbz5CeaSQyLy35K7IAS9+/9YfOuaK6BZ7XHvpvzrm3FLoWKT5qy4gsbBaoK7aLmEj/9jJQ6FqkOGnkLiISQBq5i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIAP1/W31UmlMlKkUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nb_merge_dist_plot(\n", " SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n", " SkyCoord(decam['decam_ra'], decam['decam_dec'])\n", ")" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Given the graph above, we use 0.8 arc-second radius\n", "master_catalogue = merge_catalogues(master_catalogue, decam, \"decam_ra\", \"decam_dec\", radius=0.8*u.arcsec)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(np.sum(master_catalogue['flag_merged'])/len(master_catalogue)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add SXDS\n", "\n", "This now happens in CFHT merging notebook because CANDELS contains Suprime fluxes which need merging" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# nb_merge_dist_plot(\n", "# SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n", "# SkyCoord(sxds['sxds_ra'], sxds['sxds_dec'])\n", "# )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is strange that this does not peak at zero. This is bservable in the original band cross match. It implies there is a persistent offset. Perhaps each band should be astrometrically corrected before the original merge." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#sxds['flag_merged'].name = 'flag_merged_sxds'" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# # Given the graph above, we use 0.8 arc-second radius\n", "# master_catalogue = merge_catalogues(master_catalogue, sxds, \"sxds_ra\", \"sxds_dec\", radius=1.0*u.arcsec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add UKIDSS" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8nNV97/HPT7NpG0m2FsuWLMvYsokNBoPCEvYtkIRA22yEJl2SljYNzdJ0y70t917u0qRNSUKgCWuTNhAgCU0cMDjsOAY7XsDGC97kTbItWbK1a0az/O4fM2MGIVsjaaSZeeb3fr30YuaZR/Ocsc1XR+c553dEVTHGGOMsBZlugDHGmPSzcDfGGAeycDfGGAeycDfGGAeycDfGGAeycDfGGAeycDfGGAeycDfGGAeycDfGGAdyp3KSiNwAfBdwAQ+q6jdGOeeTwP8EFNisqree7j2rqqq0sbFxvO01xpi8tnHjxk5VrR7rvDHDXURcwL3AdUArsF5EVqjq9qRzmoCvA5eo6gkRqRnrfRsbG9mwYcNYpxljjEkiIgdSOS+VYZkLgD2q2qKqw8BjwM0jzvlT4F5VPQGgqh3jaawxxpj0SiXc64BDSc9b48eSLQIWicgaEVkbH8YxxhiTISmNuaf4Pk3AlUA98KqInK2q3cknichtwG0ADQ0Nabq0McaYkVLpubcBc5Oe18ePJWsFVqhqSFX3AbuIhf27qOr9qtqsqs3V1WPeDzDGGDNBqYT7eqBJROaLiBe4BVgx4pxfEOu1IyJVxIZpWtLYTmOMMeMwZrirahi4HVgF7ACeUNVtInKniNwUP20V0CUi24GXgL9R1a6parQxxpjTk0ztxNTc3Kw2FdIYY8ZHRDaqavNY59kKVWOMcSALd2OMcaB0TYXMWY+uO/ieY7deaNM0jTG5zXruxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQBbuxhjjQO5MN2C6PLruYKabYIwx0yalnruI3CAiO0Vkj4j8/Siv/5GIHBORN+Nff5L+phpjjEnVmD13EXEB9wLXAa3AehFZoarbR5z6uKrePgVtNMYYM06p9NwvAPaoaouqDgOPATdPbbOMMcZMRirhXgccSnreGj820sdEZIuI/ExE5o72RiJym4hsEJENx44dm0BzjTHGpCJds2V+BTSq6jLgOeBHo52kqverarOqNldXV6fp0sYYY0ZKJdzbgOSeeH382Emq2qWqwfjTB4Hz09M8Y4wxE5FKuK8HmkRkvoh4gVuAFckniMjspKc3ATvS10RjjDHjNeZsGVUNi8jtwCrABTysqttE5E5gg6quAL4kIjcBYeA48EdT2Oa0UFVWbD5MKKLMqSikrqKI+hnFuAok000zxphJS2kRk6quBFaOOHZH0uOvA19Pb9OmVjAcZd2+47gKhE0HFYCrFtdw3ZJZGW6ZMcZMXt6WH+gPhAH4veV1/N0NZzKj2MOx/uAY32WMMbkhf8M9GAv3Up+b8iIPM4q99A6FMtwqY4xJDwv3wtjIVFmRh96AhbsxxhnyPtxLfLFwLy/y0DsUIqqayWYZY0xa5G24DyTC3RvvuRe6ieo7x40xJpflbbj3B8MUe10npz6WFXkA6A1YuBtjcl9eh3tiSAagrDAe7nZT1RjjAHkd7qVJ4V4e77n3WLgbYxwgb8N9YES4lxa6KRBsxowxxhHyNtxHDssUiFDqc9uwjDHGEfIy3MORKIFQ9F09d0hMh7QbqsaY3JeX4T4wHAF4T7iXFXnosWEZY4wD5GW4J5ceSFZW6LFhGWOMI+RnuAcS4e561/HyIg/BcPRk+BtjTK7Ky3AfGFF6IKGsKPb8aE9g2ttkjDHplJfhPrJoWEJiIZOFuzEm1+VtuHtcgtf17o+fKEFwtNfC3RiT2/Iy3BMLmETevaVeoufebuFujMlxeRnuIxcwJXjdBRR5XDYsY4zJeXkb7iOnQSaUFbk5YuFujMlxFu4jlBd5bFjGGJPz8i7co6oMnGJYBmLj7nZD1RiT6/Iu3AOhCFF97+rUhLIiD539QUKR6DS3zBhj0ifvwv2d1amnGJYp9KAKHX3B6WyWMcakVf6F+/Doq1MTbJWqMcYJ8i/cA6OvTk04uZDJwt0Yk8NSCncRuUFEdorIHhH5+9Oc9zERURFpTl8T02vgFBUhE8oLbZWqMSb3jRnuIuIC7gU+BCwBPi0iS0Y5zw98GViX7kamU38wggDFXteorxd5XXjdBTYd0hiT01LpuV8A7FHVFlUdBh4Dbh7lvP8NfBPI6lTsD4Yp9rkpGFF6IEFEqC0rtGEZY0xOSyXc64BDSc9b48dOEpHzgLmq+nQa2zYlYnVlRu+1J8ws8XJ8YHiaWmSMMek36RuqIlIA3AV8LYVzbxORDSKy4dixY5O99IScbnVqQmWJly4Ld2NMDksl3NuAuUnP6+PHEvzAWcDLIrIfuAhYMdpNVVW9X1WbVbW5urp64q2ehFMVDUsW67nbPHdjTO5KJdzXA00iMl9EvMAtwIrEi6rao6pVqtqoqo3AWuAmVd0wJS2epIFgGP9Y4V4aG5ZR1WlqlTHGpNeY4a6qYeB2YBWwA3hCVbeJyJ0ictNUNzCdQpEowXB0zJ57VYmPUETps71UjTE56vQpF6eqK4GVI47dcYpzr5x8s6bG4HAEgGLv2MMyAF39wyc38DDGmFySVytUA6FYuBd6Tv+xZ5bGwt3G3Y0xuSqvwj14MtxPPxWyMqnnbowxuSivwj0QjpXx9bnH6LmXJHruFu7GmNyUX+Gecs/dB2Bz3Y0xOSuvwj2YYs+9yOuiyOOynrsxJmflV7in2HMHK0FgjMlteRXuiTF37xg9d4DKUitBYIzJXXkV7sFQBJ+74JQVIZNVWgkCY0wOy6twD4SiKQ3JAMws8XHcpkIaY3JUfoV7ODLmzdSEylIvnVZfxhiTo/Iq3IPj6rl7GQ5HGYiXLDDGmFySV+E+np77yYVMNjRjjMlBKRUOc4pgKEpFsTelc0+WIBgI0lBZPJXNMsbkgEfXHRz1+K0XNkxzS1KTdz33wvH23G06pDEmB+VVuI9nzN1KEBhjclneDMtEospwJIpvjHK/Ce+U/bVwNybfnGoIJpfkTbgPx1enFrpT67mXeF343AUW7sbkge7BYVZtO8pTW47QMxSi2OuiYWYxC2v8lI6xc1u2ys1WT0CqG3UkiAiVJV6r6W6MQ6kqr+3t4j9e38+Lb3cQiiiNlcXMLi9i48ETrG05TkWxh69dtxhXwdir2rNN/oR7OBbuvhR77hAbmumyEgTGOMrQcISfbWrlh2v2sffYAMVeFxfOr+Sc+grmVBQiInxk2Ww2H+rmpxtb2X6kl7PryjPd7HHLm3APhuLDMineUIV4CQIbljHGEU4MDPMfrx/gR6/v5/jAMMvqy/n4+fWcXVeOx/Xu3+gLRDhnbgXP7WhnXUuXhXs2e6fnnvoEocoSL3s7+qeqScaYadDRG+D+V1v40ev7CUWUM2v9fOy8ehori5HTFBEsEOGCxpn8ens7Hb0BasoKp6/RaZA34Z7ouac6Wwasprsxuexg1yD3r97LExtaiUSVZXXlXL6omlnjCOnmxpm8sKODdfuP89Flc6awtemXN+Ge6LmPb1jGy1AowtBwhCJv6t9njMmcrW093PdqC09vOYyrQPjYefV84coFrNnTNe73KvW5WVpXxhsHT3D9ktqU9oLIFnkT7ifH3MdxQzW5BEG910oQGJOtQpEoq7Yd5Uev7Wf9/hP43AVcsrCKSxZUUVbkmVCwJ1w4v5ItrT1sae2muXFmGls9tfIm3AOhCAUCHlfqU5qSSxDUz7BwNybbdPQFeOy3h3hk3QHae4PMnVnEh86qpXnezLT9tt1YWUyN38e6fcct3LNRIBzF53ad9gbKSJWlVoLAmGz0VmsP//CLt9ja1ktElaaaUq5fUsuiWn9KO62Nh4hwdl05L7zdwXA4mjNDM3kT7sFQJOUFTAmVVvbXmKyhqry6u5P7XtnLa3u78LkLuPCMmVw0v5Iqv29Kr514/+MDw9SW58asmZTCXURuAL4LuIAHVfUbI17/c+CLQAToB25T1e1pbuukJHruqUjUlUisan1+RzvBcDRrS3sa42TRqPKPv9zKSzs7ONwdoKzQzQ1La7lg/sxxTZCYjKp4IcHO/qBzwl1EXMC9wHVAK7BeRFaMCO9HVfUH8fNvAu4CbpiC9k5YMBQZ1zRIiM2Jd4kwEAxPUauMMacSikR5asthvv/yXna19zOzxMvvLa/j3IYK3AXTOzRSWZqYXJE7v8Wn0nO/ANijqi0AIvIYcDNwMtxVtTfp/BIg6zYeDYQj+H2ecX2PiFDiczEQtK32jJku/cEwj68/xEOrWzjcE6CpppRPNtdzdl1Fxmq8FHpclHhddPXnTjmSVMK9DjiU9LwVuHDkSSLyReCvAC9wdVpal0bBUJSq0vH/tC/xuRkYtp67MVPtcPcQf/fzLazff5xAKEpjZQl/cNG8KblJOhGVpT7H9dxToqr3AveKyK3APwB/OPIcEbkNuA2goWF6x68DociExudKfW76AhbuxkyVrW09PLC6hae2HEFVWTqnnEsWVtEwM7umH1eVetmTQ+VIUgn3NmBu0vP6+LFTeQz4/mgvqOr9wP0Azc3N0zp0EwxHU95iL5m/0EN7b2AKWmRM/lJVXtl1jAdWt7BmTxelPjd//IFGZpR4mZHiPsfTrbLUx6aD3TkzHTKVcF8PNInIfGKhfgtwa/IJItKkqrvjTz8C7CaLBMMRwlGdUM/dX+imPxgmqll3G8GYnDMcjrJi82EeeLWFne19GZn5MlHJK9ZnlxdluDVjGzPcVTUsIrcDq4hNhXxYVbeJyJ3ABlVdAdwuItcCIeAEowzJZFJ/fFhlPBUhE/yFbqIKg8N2U9WYieoLhPjJbw/y8G/2c7Q3wOJZfj5+Xj3L5pZP+8yXiUosauzsH3ZGuAOo6kpg5YhjdyQ9/nKa25VW/fGpjL4J9dxjM2z6AqG0tsmYfNDVH+SvntjM2pYuguEoZ1SV8IcXN7JoVum4Votng6qTixpzY8ZMXqxQTdwQHU/RsISyQve73sMYM7ZjfUEeWN3Cf75+gEAowtK6ci5vqsrpGk0+j4tSn5vOHJkxk1fhPt5FTPBOz713yHruxoylLxDigVdbeGD1PoLhCDefW8e8ymJq/LmxqnMslaXenJnrnifhHgvmid5QBeizVarGnFIgFOGRdQf511/vZHA4wtl15Vz3vllTXvNlulWV+NjV3pfpZqQkL8I9MeY+kamQHlcBhZ4CG3M3ZhShSJSfbWzl7hd2c6QnwILqEq5fWpvTwy+nU1nqpe9gOF7OJLtn9+RVuE/0L8Nf6LExd2OSBEIRfrqxlQdebeHg8UGWN1Twr584h/1dg5lu2pRKLgM+pyK7Z8zkRbi/c0N1YlOu/IW2StUYiM1+eXTdQX7wyl4GhiPUzyjisxfN48xav+ODHZLnulu4Z4W+QBh3geB2TSzcywo9HOgaSHOrjMkdO4708u9r9vGLNw8zHI6yaFYply+qZn5lSc5NaZyMk9Uhc+Cmap6Ee2hCC5gS/PH6MqqaV/+QTX6LRpX/sWIba/Z00tI5gMclLG+YwQfOqKSmzBmzX8bL53bhL3TTmQMb+ORFuPcHw5O6+eEvdBOOKr1DYcqLx1c22JhcMxAM87ONrfz7mn3s7xqkvMjDDUtraW6cQbE3LyLjtCpLcmM6ZF78TfUFwuPeYi9ZYq57R1/Awt04VntvgIfX7OMn6w7SGwizvKGCixdUsWR2WcbqqGejylIfO49m/3TIvAj3/kA45S32RpOY697RF6Rplj9dzTImK+zrHOC+V/by042tRKPKWXXZWXI3W1SWeOkPhrO+OmRehHtfMDypinPJPXdjnGJf5wB3v7CbX77ZhttVQPO8GVzWVM3MkuwsuZstEnnQHwwz0529f1b5Ee6BENWlE18pd7Ln3pv942zGjOWeF/fw0tsdvHHoBK4C4ZIFVVzaVHUytMzplfoS9aZCWf2DMC/CvT8YntSKOZ+7AI9L6OizcDe5q617iHte3MPj6w9SIMLFZ1Ry+aJqC/VxKo139vqzvCSJ48NdVWM3VCcxNiYi+As9Fu4mJ+3vHOC+V1v4+cZWAC6YP5MrFtVQXmShPhGJnruFe4YFQlEiUZ10HQh/oZsO227P5JAtrd3c90oLz2w9gttVwCea6/mLqxbyys5jmW5aTjsZ7lm+at3x4f5ORcjJ3dX2F3o4Zj13k+UiUeX5He08tHofv91/HJ+7gMuaqvnAgkr8hR4L9jRwFQjFXlfWV4p1frgnioZNYiokxHru+zutBIHJTsFwhCc3tXHfK3vZ3zVIXUUR/3jjEoSJlbo2p1fqc1vPPdN64ptsFE2y517mi22UPTgctlV6JmsMDod5ZO1B7n5xN32BMHUVRdzy/rksnVNuC4+mUGk8D7KZ41Oq92S4T7LnHr/51NEbpLHK8X9sJssNDUf48doD3PfqXjr7h1lQXcInzp/Lgur8KuSVKaWFblpPDGW6Gafl+JRK9NwLvZMfloHYKtXGqpJJt8uYiQiGI/xk3UHueWkvnf1BLl1YxVeubWJXe3+mm5ZX/NZzz7y09dxtlarJoHAkys83tXL3C3to6x5iflUJv7u8jvlVJRbsGVDqczMcjjIcjma6Kafk/HCP3/SYbLiX+WyVqpl+kajy1JbDfPu5XezvGuSc+nI+uHQWC6tLbfglg3JhIZPjw71nKEShp2DCG3UkFHldeF0FtpDJTItIVPlv//UWL73dQUdfkNqywpM7HlmoZ16pL15fJov3VnZ+uA+G0rIST0So9vtsWMZMqVAkylNbDnPPi3vYe2yAGr+PW94/l7PqyimwUM8a1nPPAj1DIcrSVDuj2u+zhUxmSgwOh3l8/SEeXL2Ptu4hzqz18+kLGlg6p8xCPQudLB5m4Z45vYH09NwB5lQU8vaR7C/Sb3JHZ3+Qv/npFta2dDEUijCvspg/uGgei2r9FupZLBdKEKQU7iJyA/BdwAU8qKrfGPH6XwF/AoSBY8DnVPVAmts6IT1DIWrTtN/jvMoSntveTjgSnfQYvslvB7oGuP/VFn62sZXhcJQzZ5dxeVMV8yptmm0uSJQgyOlhGRFxAfcC1wGtwHoRWaGq25NOewNoVtVBEfkC8M/Ap6aiwePVMxRicZp2T2qsLCYUUY70BJhru9SYCdh+uJfvv7KXp7ccxl1QwMfOr6O2rIhq/8T3GzCZke2rVFPpuV8A7FHVFgAReQy4GTgZ7qr6UtL5a4HPpLORk9EzFKIsTcMyjfFe1b7OAQt3My6bD3XzvRd38/yODnzuAi5dWMUHFlal7X6QmX6lPjd9OT4sUwccSnreClx4mvM/Dzwz2gsichtwG0BDQ0OKTZy4aFTpD4bTF+7xlakHugaA6rS8p3G2zYe6+fbzu3h55zHKizxc+74aLj6jiqJJrpg2mZftJQjSekNVRD4DNANXjPa6qt4P3A/Q3Nys6bz2aPoCYVRJ2w3VGr+PIo+LfZ2DaXk/41x3/XoXz+1oZ8eRXoq9Lq5fMosLz6i0Co0Oku0lCFIJ9zZgbtLz+vixdxGRa4H/DlyhqlkxXzBRV6a8yJOWZcIiwrzK4njP3Zj3aj0xyF2/3sV/vdGG113Ate+r4ZIFVZPeLMZkn0QJgqHhSFb+JpZKuK8HmkRkPrFQvwW4NfkEEVkO3AfcoKodaW/lBPXGV4+VFbrp7B9Oy3vOryphZ7tNhzTv1jMY4t9e3sO/v7YfgMuaqrh8UbWVh3awxEKmzv5gVt6DG/NfnqqGReR2YBWxqZAPq+o2EbkT2KCqK4B/AUqBn8aXRh9U1ZumsN0pSe65pyvc51WW8PyOdiJRtXrZhuFwlK8+/iYvvt1BIBRheUMF175vFhXF3kw3zUyxxFz3jr4cDXcAVV0JrBxx7I6kx9emuV1pcTLci9M3I2F+VWw65OHuoaz8CzXTIxpVVm49wr+s2smBrkEW1pTyobNqmV1elOmmmWlSGp/p1NmfFaPQ7+Ho3xkT4Z7O6WaJRSb7u2w6ZD5SVV7d3ck/P/s22w73cmatnz/6QCOL0rSWwuSORM/dwj0DepOGZdJlfnw65P7OAS5rsumQ+eT/PL2d57d3sL9rgBnFHj5xfj3nzK2wMgF56mS496VnyDfdHB3uPUMh3PFlwumSmA65v8umQ+aLDfuP8+3nd7FmTxf+QjcfXTab9zfOtBIUec5VIBR5XBzrz85KsY4P9/IiT1rrXyemQ+7vtOmQTrfp4Am+/dwuVu/upKrUy4fPns2F82fisVA3cf5Ct/XcM6E3kL7VqckaK0vY1WHTIZ3qW6t28tz2dna291HsdXHD0louOqMSr9tC3bxbqc9tY+6ZkM66Mskaq0p44W2bDuk0ezr6uOu5Xax86yhFHhcfXDKLixdU4nNn3wIVkx1KCy3cMyIxLJNuieqQNh3SGVpPDPKd53fz5KZWijwurj6zhksXVlmpADOmUp+blmPZOUTr6HDvHQrRMAXhmyggZtMhc9vxgWH+8tFNrN13HAE+sKCKKxZVU+Jz9P8WJo0S9WWysQSBo/8V9w6FKCtM/0dMlP616ZC5aWg4wsNr9vGDl/fSHwxz3rwZXHNmja0qNeOWKEFwrC9IQ2V2dfQcG+6qOmXDMrPKfBR6Cmw6ZI6JRpUn32jjX1a9TXtvkGvfN4ulc8qYlaadukz+8ccXSB7rD1i4T5fB4QjhqKY13B9dd/Dk44oiL2v2dKbtvc3U+r9P72DlW0do6x6ifkYRf3rZGScXpBkzUSfry/Rm301Vx4b7yYqQU9BzB6gs9XK0JzsXL5h3HO4e4p+eeZtfbT5MeZGHTzbXs6zeVpWa9PAnhmWycMaMY8O9ZwpKDySbU1HEtsO99AyG0lqYzKRHIBThwdUt3PvSXqKqXH1mDZc3VdtcdZNWJT43rgKxnvt06hmc2nBPzMLZdOgEVy2umZJrmPFTVe745TaefusIxweGWTqnjA+fNZsZJXaz1KRfgQhVpV46+rLvt3jnhvsU99zrZxRRIPDGAQv3bLGno587n9rOq7uOUeP38blL5rOwpjTTzTIOV+33cazPeu7Tpje+K/lU7S7vc7uoLStk48ETU/L+JnW9gRB3P7+bH762nyKvi4+cPZuLzqi01cNmWtT4C2nvtZ77tJnqnjtAQ2Uxbx7sJhyJWoXADIhGlZ9vauWbz75N18Awn2qey19fv5hfb2vPdNNMHqnx+3irrSfTzXgPR4e7yDt3s6dCw8wS1rYcZ2d7H0vnlE/Zdcx7bW3r4S8e2cTB44PMnVHEF65YQP2MYgt2M+2q/T66+oNZV2vKseHeOxSi1OemYAr/sE/eVD3YbeE+TboHh7nruV38eO0BijwuPnZeHcsbZtjURpMxNX4fUYWugSA1/uxZEOfocJ/KIRmAGcUeqv0+Nh04wWcvmjel18p30ajyxIZD/POqnXQPDvPZi+bRMLMk6+p5mPxTHQ/0jl4L92kxVaUHkokI5zfMYOMBu6k6lX677zhfefwNDncHaKws5vcvbLCNqE3WqCnzAWTdjBkL90k6b14Fz247yrG+INV+35RfL58cOj7IN555m6ffOkJ5kYdPNc9lWX15WnfWMmayqkst3KdVbyDEGVVTP8f5/HkzgNiWbNcvrZ3y6+WDgWCY77+8l/tXt1Ag8JVrm6go8trqUpOVEp26bFvI5Nhwn66e+9I55XhdBWw6YOE+WarK3/38LZ7deoTeQJhz51Zw/dLaafl7NGaiCj0uyos8dFjPfXr0DE1PzZdCj4uldWVsssVMk7K7vY9//OVW1rYcp66iiFsvaKCh0qo2mtyQjatUHRnuwXCEQCg6bT2+Cxpn8vCafVZEbAKGhiN894XdPLi6hRKfm985t47mRpvaaHJLjd+XdT33lAYxReQGEdkpIntE5O9Hef1yEdkkImER+Xj6mzk+vUOJ0gPT87Prw2fPJhRRnt12ZFqu5xSrdx/j+u+8yg9e2cvvLq/jxa9dwQXzZ1qwm5wTC/fsGnMfM9xFxAXcC3wIWAJ8WkSWjDjtIPBHwKPpbuBEJEoPTFUt95GW1ZfTWFnML988PC3Xy3Xdg8P87r1r+OxDv2VwOMyfXDqf5Q0zWGWrS02Oqikr5FhfEFXNdFNOSqVrewGwR1VbAETkMeBmYHviBFXdH38tOgVtHLfpqCuTTES46dw6vvfibtp7A7Zt22msfOsId/xyK8cHhrlqcTVXLq7BY3V5TI6rLvURCEXpC4anrFjheKXyf1UdcCjpeWv8WNZK7JA0nSF70zlzUIWnttjQzGg6+4N84ccb+YtHNlFbXsgXr1rIdUtqLdiNIyQWMmXTph3TekNVRG4DbgNoaGiYsusc7h4CYrslTZeFNaWcVVfGijfb+Pyl86ftutlOVfn6k2+xYvNhguEo1y+ZxaVN1VlVYMmYyUqe654tewikEu5twNyk5/XxY+OmqvcD9wM0NzdP2eBUW/cQfp97WoZlkjfNnjujmGe2HuXuF3bzpWuapvza2a6jL8Adv9jGs9uOUj+jiI+dV29DVsaRavzZt0o1lXBfDzSJyHxioX4LcOuUtmqS2rqHprXXnrCsvoJntx5lc2v3tF87m6gqT25q486ntjMUinDD0louWVhlvXXjWIniYdkU7mMOeKpqGLgdWAXsAJ5Q1W0icqeI3AQgIu8XkVbgE8B9IrJtKhs9lrYTQ9TNmP5wLy/y0FhVwuZD3Vl113w6HTo+yB//cD1f++lmFtaU8syXL+PyRTYMY5ytrNCNz12QVXPdUxpzV9WVwMoRx+5Ierye2HBNVjjcM8R58yoycu3zG2bws02tvLCjg2uXzMpIGzIhHInylz95g+d3tCMIHzl7NhcvqGRdy/FMN82YKSci1JRl1ypVx01VGAiG6R4MUVdRnJHrnzO3gpndBZ2sAAANyklEQVQlXu56blfe9N7fOHiCm+5ZwzNbj7KgupSvXNvEJQurbDGSySvVpdm1kMlx4f7OTJnM3LhzFQhXL65h+5Fexy/K6R4c5utPvsXvff81jg8Mc+sFDXz2onlUFHsz3TRjpl2NvzB/p0JOh9Z4uNdnYMw94Zy5FWw8eILvPL+LDy6ZNaVb/WVCJKr89RObWbX9KIFQhEsWVHHNmTX4PLYrkslfNWU+1u7rynQzTnJwzz1z4e4qEL58TRNvH+3j2W1HM9aOqbDp4Al+5941/NebbdT4fXzxqoV8+OzZFuwm79X4fXQPhgiGI5luCuDAcG87MYS7QDK+l+FHz5nDguoSvvP8LsKRrKjKMClHeob46uNv8nv/9hodfQE+1TyXP73sDNvuzpi46iyb6+64cD/cPURteWHGp965CoSvfXAxu9r7+cErezPalskYGo7w3ed3c/W3XuHpt47wF1cu4MWvXck5cytsuztjkiRGCw4dH8pwS2IcN+be1j1EXQaHZJJ96Kxablw2m+88v5vLF1WzrD4z0zMnIhpVnnyjjTt/tY3eQJiz6sq5YWktM0u8Vv3SmFEsmuUHYHdHHxcvqMxwaxwY7oe7A1w4f2amm3GyLMHyuTNYvbuTz/1wPa/+7VUUe7P/j/w3uzv5p2d2sO1wL/Uzirjl/Q00VtmuSMacTo3fR3mRh51H+zLdFMBhwzLhSJSjvYGMrE49lSKvi4+fX09n/zD/b+WOTDfntLa29fDZh9bxmYfW0T0Y4ru3nMufX7HAgt2YFIgIi2aVsqs9O8I9+7uR49DeFyQS1YzOlBnNgupSLl1YxY/XHuSMqlI+l2VVI7cf7uXuF3bz7LajFHtdfOTs2Vw4fyYDwYgtRDJmHBbN8vPUliOoasbvSTkq3NtOxG5kZMuYe7IPLp1Fqc/NnU9tJxSJ8mdXLMh0k3jzUDfff3kPq7a14/e5uWpxDZc1VVFo0xqNmZDFtX4eWXeQjr5gxiugOircs2GO+6m4Cwr43q3L+erjb/JPz7xNKBLli1ctnPaf7uFIlFXb2nnoNy1sOtiNv9DNl69p4nOXzOfpt2yjEWMmo6kmdlN159E+C/d0auvO3p47gMdVwHc+dS4eVwHf+vUu1u07zj/euOTkXfaptKejj59uaOXJN9o41hdkZomXG5fN5vyGGfg8Lgt2Y9Jg0azYRh272vu4fFF1RtviqHBvPTHEzBIvRd7sHVZwuwr41ifO4ey6cr7z/C4+9N3V/P6FDXyyeS5LZpelrVSBqrLtcC/PbW/nue3tbD/Si6tAuGpxDbPLC1lc67fxdGPSrLLUR1WpLytuqjoq3A9n0Rz30STv2lTocfGXVzfx/I52frz2AP/x+gEqS7xcsrCKpppS5lQUMaeiiJklXsqLPJQXeSj0FLxnGEdV6Q2EOdw9xJGeIXYc6eONgyd442A3XQPDCNAws5gPn1XLOXMr8GfJ5r3GONWiWaXsbO/PdDOcFe5t3UMsqM6daXslPjc3n1vH1WfWsKejn90d/by8s4MVm0dfJCQCPncBPrcLVSUUUYYjUSLRd5cWXlBdwlVn1qCqLK4to9TnqL9mY7Laoll+nthwiGhUM1o00DH/16sqh7uHuLwps+NcE+Ev9LC8YQbLG2YAEIpE6RkK0T0YYnA4TCAUZSgUIRSJxr8UAdwFQkGBUOx1UVHspaLIQ2WpNycWShnjVItr/QwOR2jrHmLuzMzsKwEOCvdYEEYyVsc9nTyuAqriY3fGmNySfFM1k+HumBWqbVlQx90YY5ris992ZvimqmPCfU9H7AZG/YzM/aQ0xpiyQg+zywvZneGbqo4J92e2HqHG7+N9s8sy3RRjTJ5bNMuf8QJijgj3vkCIl3Ye48Nnz854HXdjjFlc62fPsf6MbtTjiHB/bns7w+EoHz1nTqabYowxNNWUMhyOcuD4YMba4Ihw/9Xmw9RVFHFeQ+5shmGMca7z5sWmNf/qFGtWpkPOh3v34DCrd3dy47LZGS+xaYwxECvz/cEls3joN/voGQplpA05H+7Pbj1KOKrcuMyGZIwx2eNL1zTRFwjz72v2ZeT6OR/uT205wrzKYs6qs1kyxpjscVZdOdctmcXDGeq9pxTuInKDiOwUkT0i8vejvO4Tkcfjr68TkcZ0N3Q0nf1BXtvbyUeXzbEhGWNM1vnyNU30BsL8cM3+ab/2mOEuIi7gXuBDwBLg0yKyZMRpnwdOqOpC4NvAN9Pd0JHePtrL53+0gahis2SMMVkp0Xt/6DctHOkZmtZrp9JzvwDYo6otqjoMPAbcPOKcm4EfxR//DLhGpqgrHQhF+Naqndx49284dHyQ7316OYtrp36zC2OMmYivXNvEUCjCZd98iS8+sonX93ahqmN/4ySlEu51wKGk563xY6Oeo6phoAeoTEcDR7rnxT3c89Iebjp3Ds//1RXWazfGZLWlc8p57qtX8MeXNPKbPZ18+oG1PLC6ZcqvO61VIUXkNuC2+NN+Edk50ff6dvxrHKqAzoleLwfY58tdTv5s4PDP9/sT+J4/+yb82cQvOS+Vk1IJ9zZgbtLz+vix0c5pFRE3UA50jXwjVb0fuD+VhqWbiGxQ1eZMXHs62OfLXU7+bOD8z5etUhmWWQ80ich8EfECtwArRpyzAvjD+OOPAy/qdAwqGWOMGdWYPXdVDYvI7cAqwAU8rKrbROROYIOqrgAeAv5TRPYAx4n9ADDGGJMhKY25q+pKYOWIY3ckPQ4An0hv09IuI8NB08g+X+5y8mcD53++rCQ2emKMMc6T8+UHjDHGvFdehPtY5RNymYg8LCIdIrI1021JNxGZKyIvich2EdkmIl/OdJvSSUQKReS3IrI5/vn+V6bblG4i4hKRN0TkqUy3Jd84PtxTLJ+Qy34I3JDpRkyRMPA1VV0CXAR80WF/d0HgalU9BzgXuEFELspwm9Lty8COTDciHzk+3EmtfELOUtVXic1QchxVPaKqm+KP+4iFxMjV0TlLYxK7KHviX465CSYi9cBHgAcz3ZZ8lA/hnkr5BJPl4pVGlwPrMtuS9IoPW7wJdADPqaqTPt93gL8FMreRaB7Lh3A3OU5ESoGfA19R1d5MtyedVDWiqucSW/l9gYiclek2pYOI3Ah0qOrGTLclX+VDuKdSPsFkKRHxEAv2R1T1yUy3Z6qoajfwEs65f3IJcJOI7Cc2FHq1iPw4s03KL/kQ7qmUTzBZKF42+iFgh6relen2pJuIVItIRfxxEXAd8HZmW5Ueqvp1Va1X1UZi/8+9qKqfyXCz8orjwz1egjhRPmEH8ISqbstsq9JHRH4CvA4sFpFWEfl8ptuURpcAnyXW63sz/vXhTDcqjWYDL4nIFmKdkOdU1aYMmrSwFarGGONAju+5G2NMPrJwN8YYB7JwN8YYB7JwN8YYB7JwN8YYB7JwN8YYB7JwN9NGRCLxuerb4mVuvyYiBfHXmkXk7tN8b6OI3Dp9rX3PtYfiNWCygoh8Kl7C2ubFm1FZuJvpNKSq56rqUmKrMT8E/A8AVd2gql86zfc2AhkJ97i98RowKYuXm54Sqvo48CdT9f4m91m4m4xQ1Q7gNuB2ibky0QsVkSuSVqS+ISJ+4BvAZfFjX433pleLyKb41wfi33uliLwsIj8TkbdF5JF4GQNE5P0i8lr8t4bfiog/XpXxX0RkvYhsEZE/S6X9IvILEdkY/y3ktqTj/SLyryKyGbj4FNdcGn/8ZvyaTfHv/UzS8fsSPxzim81sir/HC2n8azBOpqr2ZV/T8gX0j3KsG5gFXAk8FT/2K+CS+ONSYhu5n3w9frwYKIw/bgI2xB9fCfQQKxBXQKw0w6WAF2gB3h8/ryz+vrcB/xA/5gM2APNHtLER2Dri2Mz4f4uArUBl/LkCn4w/PtU1vwf8ftI5RcD74p/bEz/+b8AfANXESlbPT75u0md9arQ/a/uyL/c4fxYYMx3WAHeJyCPAk6raGu98J/MA94jIuUAEWJT02m9VtRUgPk7eSCzwj6jqegCNlw4WkQ8Cy0Tk4/HvLSf2w2LfGG38koj8bvzx3Pj3dMXb8vP48cWnuObrwH+Pb2bxpKruFpFrgPOB9fHPWkSsxvtFwKuqui/+Ho7cmMWkn4W7yRgROYNYGHYQ67kCoKrfEJGngQ8Da0Tk+lG+/atAO3AOsR56IOm1YNLjCKf/dy7AX6rqqnG0+0rgWuBiVR0UkZeBwvjLAVWNnO77VfVREVlHbJeilfGhIAF+pKpfH3Gtj6baLmOS2Zi7yQgRqQZ+ANyjqjritQWq+paqfpNYtcQzgT7An3RaObFecZRY5cixbl7uBGaLyPvj1/CLiJtYtdAvxOvGIyKLRKRkjPcqB07Eg/1MYr3rlK8Z/6HWoqp3A78ElgEvAB8XkZr4uTNFZB6wFrhcROYnjo/RNmMA67mb6VUUHybxENv8+j+B0eq0f0VEriK2Pds24Jn440j8RuUPiY1J/1xE/gB4Fhg43YVVdVhEPgV8L147fYhY7/tBYsM2m+I3Xo8BvzPG53gW+HMR2UEswNeO85qfBD4rIiHgKPD/VPW4iPwD8Ov49NAQ8EVVXRu/Yftk/HgHsZlGxpyWlfw1ZgwS27/1KVXNqi3w4sNDf62qN2a6LSb72LCMMWOLAOXZtoiJ2G8vJzLdFpOdrOdujDEOZD13Y4xxIAt3Y4xxIAt3Y4xxIAt3Y4xxIAt3Y4xxoP8PeniMfZ//NY4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nb_merge_dist_plot(\n", " SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n", " SkyCoord(ukidss['ukidss_ra'], ukidss['ukidss_dec'])\n", ")" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Given the graph above, we use 0.8 arc-second radius\n", "master_catalogue = merge_catalogues(master_catalogue, ukidss, \"ukidss_ra\", \"ukidss_dec\", radius=0.8*u.arcsec)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(np.sum(master_catalogue['flag_merged'])/len(master_catalogue)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add VIRCAM" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8nNV97/HPbzSaGUmjfbNsSbaxRcCAMWAMhCRkIyEbtDe5KaVpb26buu0NTdLt3uR1e2mbtvcm6avpFtqGEpouEJJCAAVMuNwAWQh4t/GGsfCi1ZZs7dtolt/9Y2aUQcjWyHpGM/PM7/166ZWZZx7Nc8bEX585zzm/I6qKMcYYd/FkuwHGGGOcZ+FujDEuZOFujDEuZOFujDEuZOFujDEuZOFujDEuZOFujDEuZOFujDEuZOFujDEu5M3Whevq6nTNmjXZurwxxuSl3bt3n1XV+oXOy1q4r1mzhl27dmXr8sYYk5dE5FQ659mwjDHGuJCFuzHGuJCFuzHGuJCFuzHGuJCFuzHGuJCFuzHGuJCFuzHGuJCFuzHGuJCFuzHGuFDWVqhm20PbO+c9ftcNrcvcEmOMcZ713I0xxoUs3I0xxoUs3I0xxoUs3I0xxoXSCncRuU1EjopIh4h8fp7X/0pE9iV+XhORYeebaowxJl0LzpYRkSLgXuBWoBvYKSLtqno4eY6q/k7K+b8NXJOBthpjjElTOj33LUCHqh5X1RngYeCOC5z/i8C3nGicMcaYi5NOuK8CulKedyeOvYmIrAbWAs8tvWnGGGMultM3VO8EHlHV6HwvishWEdklIrsGBgYcvrQxxpikdMK9B2hJed6cODafO7nAkIyq3qeqm1V1c339gvu7GmOMuUjphPtOoE1E1oqIj3iAt889SUQuA6qBl5xtojHGmMVaMNxVNQLcDTwDHAG+o6qHROSLInJ7yql3Ag+rqmamqcYYY9KVVuEwVd0GbJtz7J45z//YuWYZY4xZCluhaowxLmThbowxLmThbowxLmThbowxLmThbowxLmThbowxLmThbowxLmThbowxLmThbowxLmThbowxLmThbowxLmThbowxLmThbowxLmThbowxLlTQ4T4eihCz8vPGGBcq2HCPxpSvPnuU517tz3ZTjDHGcQUb7hOhCNPhGNuPnyMSjWW7OcYY46iCDffxUASAiZkoh/pGs9waY4xxVsGHu0dg+/HBLLfGGGOclVa4i8htInJURDpE5PPnOefjInJYRA6JyEPONtN5yXC/bnU1J89NcGZ0OsstMsYY5ywY7iJSBNwLfADYAPyiiGyYc04b8AXgZlW9AvhcBtrqqIlEuN9yaQNFHmH7Ceu9G2PcI52e+xagQ1WPq+oM8DBwx5xzfh24V1WHAFQ156egjE9H8HqE6tJirlpVyd7OIUKRaLabZYwxjkgn3FcBXSnPuxPHUl0KXCoiL4rIyyJym1MNzJTxUISg34uIcMPaGkKRGK90j2S7WcYY4wivg+/TBrwTaAZ+JCJXqepw6kkishXYCtDa2urQpS/OeChCMBD/+K01pZT5vXQNTma1TcYY45R0eu49QEvK8+bEsVTdQLuqhlX1BPAa8bB/A1W9T1U3q+rm+vr6i22zI5I9dwARoSLgnb3Jaowx+S6dcN8JtInIWhHxAXcC7XPOeZx4rx0RqSM+THPcwXY6LjXcAYJ+C3djjHssGO6qGgHuBp4BjgDfUdVDIvJFEbk9cdozwDkROQw8D/yBqp7LVKOXKhZTJkIRylLCvczC3RjjImmNuavqNmDbnGP3pDxW4HcTPzlvZCpMTHlTz30iFEFVEZEsts4YY5auIFeonh0PAczeUIV4uIejysSMTYc0xuS/Ag33GWBOzz0R9GfHQllpkzHGOKlAwz3Rc58zLJP6mjHG5DML94SfhftMVtpkjDFOKshwPzc+g0egxFc0e8x67sYYNynIcD87HqLM58WTMiumzMLdGOMiBRvuqTNlAIo8QklxkYW7McYVCjTcZ96wgCkpGPBydszG3I0x+a9Awz30hpupSUG/l3MT1nM3xuS/ggt3Vb1guNtsGWOMGxRcuE/ORJkOx84f7raIyRjjAgUX7vPNcU8KBryMhSJMh60EgTEmvxVuuAfm77kDnJuwoRljTH4ruHAfSMyGmXe2jN/qyxhj3KHgwj05G+Z8Y+5gC5mMMfmv4ML97GzPvehNr1m4G2PcovDCfTxEZUkxXs+bP/ps2V+bDmmMyXMFF+7nJkLUBX3zvlZc5EnMdbeeuzEmvxVcuJ8dm6Eu6D/v67VBn/XcjTF5r/DCfTxEXfn5w70u6LfZMsaYvJdWuIvIbSJyVEQ6ROTz87z+SREZEJF9iZ9POd9UZ5wdD1FXNv+wDEBd0Gf1ZYwxee/N8wHnEJEi4F7gVqAb2Cki7ap6eM6p31bVuzPQRseEIlFGpyMXHJapC/rZeXJoGVtljDHOS6fnvgXoUNXjqjoDPAzckdlmZcbIZBiA6gv03GuDfoYmZ4hEY8vVLGOMcVw64b4K6Ep53p04NtdHReQVEXlERFrmeyMR2Soiu0Rk18DAwEU0d2nGQhFg/gVMSfVBH6owaCUIjDF5zKkbqt8D1qjqRuBZ4F/mO0lV71PVzaq6ub6+3qFLp28iEe7zlR5ISg7ZDNh0SGNMHksn3HuA1J54c+LYLFU9p6rJNLwfuM6Z5jlrIhSv9jjf6tSk5EyaczYd0hiTx9IJ951Am4isFREfcCfQnnqCiDSlPL0dOOJcE50zkcawTG1iPN4WMhlj8tmCs2VUNSIidwPPAEXAA6p6SES+COxS1XbgMyJyOxABBoFPZrDNF21iJo1hmUTP3cLdGJPPFgx3AFXdBmybc+yelMdfAL7gbNOcN55Gz73c78Xn9dgqVWNMXksr3PPdQ9s7AfjJsbMAbHulD3/x/OPuIkK9rVI1xuS5gio/EIrE564Xey/8sWvKfAxNWs/dGJO/Civcw1F8Xg8ekQueV1VazGBiwZMxxuSjwgr3SAx/0cIfuabMx5AtYjLG5LGCCveZaAzfAkMyANWlFu7GmPxWUOEeCsfwF6fXcx8LRZiJWH0ZY0x+Kqxwj8Twe8+/OjUpWVhs2G6qGmPyVEGF+0wkij+NYZma0ni4D9lNVWNMniqocA9F0hxzLysGrDKkMSZ/FVy4pzMsU1OW7LlbuBtj8lOBhfvihmWs526MyVcFE+4xVcJRTWtYpio55m7hbozJUwUT7slpjen03H1eD0G/l0EbljHG5KmCCffQbLgvPOYO8Zuq1nM3xuSrAgr3+C5M6fTcIT7ubvVljDH5qmDCfTHDMhBfyGQ9d2NMviqYcE8Oy/jSKD8AiZ67hbsxJk8VxGYdEK8rAwuPuSc39jgzOs3AeGj2+V03tGa2gcYY46AC6rkvbsy9zO9lJhIjHLXiYcaY/JNW0onIbSJyVEQ6ROTzFzjvoyKiIrLZuSY6I7TIMfdSX/xLzeRMNGNtMsaYTFkw6USkCLgX+ACwAfhFEdkwz3nlwGeB7U430gnJG6rpLGICKPXFh28mZyIZa5MxxmRKOkm3BehQ1eOqOgM8DNwxz3l/CnwZmHawfY4JRWII4EtjJyaAUn883CdC1nM3xuSfdJJuFdCV8rw7cWyWiFwLtKjqUw62zVEzkfj+qbLA/qlJZbPDMtZzN8bknyXfUBURD/BV4PfSOHeriOwSkV0DAwNLvfSixCtCpv9xk8MyEzbmbozJQ+mkXQ/QkvK8OXEsqRy4EnhBRE4CNwLt891UVdX7VHWzqm6ur6+/+FZfhHgt9/RKD0DKDdWQ9dyNMfknnXDfCbSJyFoR8QF3Au3JF1V1RFXrVHWNqq4BXgZuV9VdGWnxRUq33G9SkUcIFHus526MyUsLpp2qRoC7gWeAI8B3VPWQiHxRRG7PdAOdsthhGYiPu9uYuzEmH6W1QlVVtwHb5hy75zznvnPpzXLeTCRGVUnxon6n1Fdk89yNMXmpgFaoprd/aqoyv9fG3I0xeamgwj3dWu5Jpb4iG3M3xuSlggn3mUXeUIX4jBkbczfG5KOCCPdoLLF/aprlfpPKfEWEozpbusAYY/JFQYT7zCK32Esq9dsqVWNMfiqIcF9sud+kMlulaozJUwUS7osr95tkq1SNMfmqIMJ9sfunJiUrQ9pcd2NMvimIcJ/dP3WRY+7JypATNuZujMkzBRHuMxc55l7iK0KwnrsxJv8URLhf7Ji7R4RAcRETNuZujMkzadWWyXehRW6xl6rMb/VljDE/89D2zjcdu+uG1iy05MIKrOe+uDF3iI+7j1vP3RiTZwqk5x5FgOKi9LbYS1Ue8HJmNOR8o4wxOaejf4yfvn6O/V0j/PjYAD6vhw9d1cTq2rJsN23RCiTcY/iL098/NVUwUEzHwHgGWmWMyQVdg5O07+/le/t7efX0GAB1QR+1ZX7OjE7z9R8dZ8vaGt6/YQUlvsV/+8+Wggj3mfDiK0ImVQS8TIdjTIejBIrz5z+sMeb8JmciPH3gNP+xu4uXjw8CsLqmlI9sbOLypgoqS4oREUKRKP/v8Bl++vo5jg9M8Ln3tuG5iE5iNhREuIci0Yu6mQoQTNSXGRgL0VJT6mSzjDHLSFXZdWqIR3Z189SBPsZDEVbXlnLrhkauaamiqtT3pt/xe4v40MaVNFQEeGxvD73DUzRX50cOFEi4L36LvaTyQHz3pn4Ld2Py0rEzY7Tv76V9fy+nzk3iK/Jw1apKrl1dzZra0rSGay9vquDxvT28dmbMwj2XzFzELkxJ5YFkz33aySYZYzJEVTnWP87TB07z9ME+Xj09hkfgpnW1bF5dw5WrKhY9TBv0e1lZVcJrZ8Z592WNGWq5s9IKdxG5DfgboAi4X1W/NOf13wQ+DUSBcWCrqh52uK0XLRSJUe1/81eudPws3G3GjDG5SlU52DPK0wf7+P6h0xwfmECA1ppSPryxiatWVc5+C79YlzYGeeHoAFMz0by4sbpguItIEXAvcCvQDewUkfY54f2Qqv5j4vzbga8Ct2WgvRcldBG7MCWV+b0I8WEZY0zuUFUO943Svr+XJ/f30TM8RZFHuGFtDVeurGTDygoqlhjoqS5tLOf5owMc6x9jY3OVY++bKen03LcAHap6HEBEHgbuAGbDXVVHU84vA9TJRi7VUsbcPSIE/V7ruRuTI/pHp3l0Tw+P7ummo38cj8D6hiAfvXYVl6+omN1kx2nN1aUEij0cOzPumnBfBXSlPO8Gbph7koh8GvhdwAe825HWOWRmCeEOEAx4reduTBaFIlF+cKSfR3Z388PXBojGlM2rq7lj00quXFlJWYYCPVWRR2hrKOe1/jFU9aLWzSwnx/5EVPVe4F4RuQv4Q+C/zD1HRLYCWwFaW5enFkM4GiMS00WX+01VHrCeuzHLLRZT9nYN8djeHr63v4+RqTAVAS9vW1/Hda3V1JX7l71NlzYGOdAzwunRaZoqS5b9+ouRTrj3AC0pz5sTx87nYeAf5ntBVe8D7gPYvHnzsgzdJCs6LqXnXu4vpnt40qkmGWPOQ1U51DvK91752Ti63+vh/VesoLbMx7qGYFYXEbU1lAPw2plxV4T7TqBNRNYSD/U7gbtSTxCRNlU9lnj6IeAYOWLciXAPeDk7PkMspng8uf1VzJh8o6q80j3CtoN9PH3gNJ2Dk7Pj6P/5umYub6rImdXhFSXFNFUGeO3MGLdcWp/t5lzQguGuqhERuRt4hvhUyAdU9ZCIfBHYpartwN0i8l4gDAwxz5BMtkyEEht1LOH/HMGAl2hMGZycoS64/F8FjXEbVWVP5zDbDvTx/YOn6RmewusRbl5fx+bV1VzeVLEs4+gXo60hyE86zhKOxiguyt3Cumn96anqNmDbnGP3pDz+rMPtcszYdBiAwJJ67vHpVANjIQt3Y5bgSN8oj+/recPUxbaGIB+7tpnLmspnN6XPZY0VAWIKI5PhrIz7pyv3/ySXaDQZ7kvouZcnehD9YyEub3KkWcYUjIGxEH/Ufoi9nUP0jUzjkfjYda4NuaQrWYNmaHLGwj2bxqbjY+4lSwl3W6VqzKLMRGI892o/j+zu4vmj8amLzdUlfOTqlWxctTxTFzOlujT+TX54MpzlllxY/v4Jp2l0Kv4fwF+89GGZfqsvY8x5JWe6PLK7m/b9vQxOzNBQ7udTb19LwFtEY0Ug2010RHmgGI/Ee+65zP3hnui5L+Wrn8/rsVWqxpzHwFiIJ/b18Mjubl49PUaRR7i8qYKPbGxifUM5RS6bYVbkESpLii3cs210OozXI0u+q11f7rdVqsYkTIejPPdqP4/u7uaFxIrRTS1V3H71SjY2V+bFjdGlqCr12bBMto1NRxy5YVNf7reeuyloqsruU8kVo72MTkeoCHi5eV0t17ZW0+CSYZd0VJf66Ogfy3YzLsj14T46FSawhPH2pPpyP4d7Rxc+0RiX6egfo31fL4/t66FrcIpAcXzFaE1p9leMZkt1aTFj0xEi0RjeHJ3r7vpwd6rn3lDu54fWczcFond4iidf6eWJfb0c6h1FgHUNQT52XTNXNFUsaVGgG1SX+lBgZCpMbY6ufXF9uI9Ohx0K9wDjoQiTMxHXjyeawjQ8OcNTB/p4Yl8vO07EN42+uqWKD13VxFXNlY7WRs93VYnpkEOTFu5Z4+SYO8RnBqyudf0fmykQ0+Eoz7/az2N7e/jBkX6iqtQH/bz38kaubq7M2eDKturEQqbhHJ4x4/qUGp0KU1Wy9B5HQyLc+8dCrK4tW/L7GZMtsZiy4+Qgj+/t4akDfYxNR6gv93PTulo2tVTRVBnI+Vrl2VZRkvtz3d0f7tPhJa1OTUrtuRuTj14fGOexPT08treHnuEpfEUerlhZwdUtVayrD7puPnomFXmEikAxQzk8HdLV4T4TiTEdjjly86fBwt3koamZKE8d6OPhHZ3sOjWER+DtbfXcvL6ODU0V+JZQUK/Qxee6W889K5IVIUscmApZXerD6xErQWDywtHTYzy0/RTf3tXFdDhGXdDHbVesYFNrld0YdUh1aTHHz05kuxnn5fJwX3rpgSSPR6gL2kImk7vGQxG2Hejj2zu72H1qCF+Rh8uaytmytoa1tWU2ju6w6jIfo13DRGKxbDdlXq4OdyfK/aayEgQm10SiMV46fo6//L+vcah3hHBUqQv6+MCVK7i2tTqvqy/muqqSYhQYnYpkuynzcvV/eSd77hAv0t89ZHupmuxK7mLUvi8+2+Xs+AyBYg+bWqq5trWK1ppS66Uvg+qyn9V1z0WuDvdkuV8nyg8ArK4t5cWOs6iq/eUxy66jf4zH9/by4PZTDE3GC+JdtqKc921YwVtWlOf0lm9ulOtz3V0d7k733NfUljIVjjIwFiqoIkkmewbGQrTv7+XxvT0c6BnBI3BJfZD3XNbIhpX5t4uRm1SUeBHI2emQrg730dnZMs78BWhNLF46eW7Swt1kzHQ4yrOHz/C15zo41j9GTGFlVcDKAOQYr8dDRUkxQxN53HMXkduAvwGKgPtV9UtzXv9d4FNABBgAflVVTznc1kUbnY4ggmNzedfUlgJw8twEW9bWOPKexkB81ej2E4M8trebpw+cZiwUobKkmLe31bOppco1uxi5TVVpMcNTedpzF5Ei4F7gVqAb2Cki7ap6OOW0vcBmVZ0Ukd8CvgL8QiYavBijU2GCfq9jJUlXVZXg9QinzuXu3FaTXzr6x/junh4e2t7J8FQYn9fDlSsr2dRSxSX1ZQVZTjefVJf6OJmjeZBOz30L0KGqxwFE5GHgDmA23FX1+ZTzXwY+4WQjL9bYdMTRr7DeIg/N1SWcPGczZszFG5kK076vh+/s6p4dR1/fEOT9V6zgcls1mleqSosZ7Q7nZF33dMJ9FdCV8rwbuOEC5/8a8PR8L4jIVmArQGtra5pNvHij02HKA87eVlhdW2Y9d7NoqsquU0P8+VNHONgzQiSmNFUG+OBVTVzdXDm7CbvJL9UlPmIKp0enaa4uzXZz3sDR5BORTwCbgVvme11V7wPuA9i8ebM6ee35jE2HqXCgImSqNbWl7Dk1ZNMhTVpGpsI8vreHB7ef4rUz4/i9Hq5bXc3m1TWsrLLqi/kumS9n8jTce4CWlOfNiWNvICLvBf4ncIuq5sQyztGpCCurnL0Rtbq2jLFQhMGJGat1bealqrzSPcKD20/Rvr+X6XCMq5sr+cpHNzI5E7VhFxepKIlH6JnRnIi8N0gn3HcCbSKylnio3wnclXqCiFwDfB24TVX7HW/lRRoLhakIlDv6nqsTM2ZODU5auJs3GA9FeGJf/Obood5RfEUerm6pZMuaWlZVlxCJqQW7yyTv6Z0eyb2CgguGu6pGRORu4BniUyEfUNVDIvJFYJeqtgN/AQSB/0h8zexU1dsz2O60jE5FHBtzf2h7J8BsVchvbe/k2tZqR97b5LeDPSM8tKOTJ/b2MDET5bIV5dx+9Uo2tVTZIiOXK/UVUeQRzozmYbgDqOo2YNucY/ekPH6vw+1aMlXNyJh7TakPAc7l6MIFszzGQxHa9/XyrR2dHOgZobhIuGpVFVvW1tBSXWJj6QVCRKgIePM33PPRxEyUmOL4bBlvkYfK0mIGLdwL0sGeER7c3kn7vngv/S2N5XxkYxObWqop8VkvvRBVBIo5beG+fJIbdVQEiok5PC+ntszHufHcu4FiMmMiFKF9fy9fe66DnuGpn/XS11TTYhUYC15FSTH9eXpDNS8layyXB4oZcXh5cG2Zn4O9I46+p8ktqsqBnhG+taOT9n29TMxEaazwWy/dvElFwMvrA+M5Nz3ateE+23Mv8Toe7jVlPiZnooxMhqkstcUnbjIyGeYPHz/ArlND9I1Mz/bSr19TbXXSzbwqSoqZnIkyFnJ2RfxSuTbckxUh4yv/phx979pgvI7zqcEJNpZWOfreZvnFYspLx8/xnV1dPH3wNDORGCurAjbjxaQlGej9o9MW7sshWcu9wuEbqhAfloF46d+NzRbu+ap3eIr/9cRB9pwaYmgyTKDYw7WtVYnVoyXZbp7JE+WJhUynR0Ksb3B2Xc1SuDbck7swZaJmR01ie61TObzzuZlfJBrjB6/2860dnfzwtQFUYV19Ge/bsIINKytsNyOzaJXJhUw5NmPGveE+nbyh6vxH9Hk9VAS8nBq06pD5YmAsxMM7OnloRyd9I9M0Vvi5+13r8XuLZv+xNuZiJDuQuTbX3cXhHq+Nnanx0poyv1WHzAP7u4b55k9P0r6vl6gq6xuCfOKG1bxlRTlFHrs5apYu2dmzcF8mTtdyn6u+3M/R06PEYorHQiKnhKMxnj54mn9+8QR7O4cJ+r3ccEkNN66tpa7c6gEZ562oDFi4L5fRqXBGbqYmtdaUsPPkIMfPjufUTZRCNjQxw0M7Ovn6D19ndDpCbZmPD29s4rrWavw248VkUGNFgNM5tpDJveE+HaHc4boyqVpq4tUhd58asnDPstcHxnngJyd4dE830+EY6xuC/Nw1tVzaWG7b1Jll0VgR4NiZs9luxhu4NtzHpjPbc68L+qkqLWb3qSF+4frM7ypl3kg1Pjf9Gz8+wQ9e7cfn9fDzm1bxq29by+5TQ9lunikwKyoCDIyHiMY0Z+7luDbcR6fCrKzM3FxljwjXtlazp3M4Y9cwbxaKRHlyfx/f+MkJDveNUuYr4j2XNXDDJbUE/V4LdpMVjZUBojHl3HiIhgpnNwi6WK4N97Fp52q5n8+1rVU892o/w5MzVJXadLpMOjse4sGXO/m3l09xdjzE+oYgP3/NKja1VNncdJN1jYkb9adHpy3cM200A7Xc57p2dXyzjr2dw7zrsoaMXqtQ7esa5o/bD3GgZ4RoTLm0MchHNjaxviFodV5MzlhRGQ/0XNpuz5XhPhOJMR2OUe7P7Me7urmKIo+wp3PIwt1BoUiUbQf6+OZPT7G/axi/18OWNTXccEkNDeW50SsyJtWKRG89l1apujLcf1YRMrM99zK/l8ubym2c1yEDYyH+x6OvsP3EIBOhCHXBeInda1qrrXiXyWm1QT9FHqHfwj2zxjJYemCua1ureWR3N5FoDK+N/V6Ujv4x/ulHJ3hsXw8zkRiXrSjnpnW1rK+3oReTH4o8Qn3Qn1MbZaeVfiJyG/A3xDfIvl9VvzTn9XcAfw1sBO5U1UecbuhijKbswpRp162u5l9fOsXRM2NcsbIy49dzkz2dQ/zDC6/z7OEzBIo9fHxzMw3BgK0iNXmpsTKQX8MyIlIE3AvcCnQDO0WkXVUPp5zWCXwS+P1MNHKxBsbiNzWSddcz6drW+E3VPaeGLNzToKr8pOMs9zxxiBNnJygpLuLdlzVw0yW1lGX4HokxmdRY7ufUudwpJpjO36YtQIeqHgcQkYeBO4DZcFfVk4nXYhlo46J1Jao1JleRZlJzdQn15X52nxril29ak/Hr5atYTHn2yBn+/vkO9nePUBHw8sGrmrh+TTV+r42nm/y3ojLA9hOD2W7GrHTCfRXQlfK8G7ghM81xRufgFCXFRdQuQylXEeG61mp2nhzKuT0Uc0E0pjz5Si/3Pt/Ba2fGaa0p5f/8p6sIR+wehXGXxooAI1NhpsPRnJgAsKzfg0VkK7AVoLU1c0v2u4Yml3W/y3ddVs/3D51mX9cw1ySGaQpdKBLlsT09/MUzRzk3MUNDuZ+Pb27mqlVVqGLBblynMTkdcmSaNXVlWW4NpPM3rAdoSXnenDi2aKp6n6puVtXN9fX1F/MWaekanKSlZvm2SfvgVU0Eij08srt72a6Zq8ZDEe7/8XHe8ZXn+fx3DxAoLuKXbmjlM+9pY1NLdc7U3TDGaa2JYeBc2cQnnZ77TqBNRNYSD/U7gbsy2qolUFW6Bie5aV3tsl2zPFDMB65son1/L//rwxty4ivZcjs3HuL3/2M/Lx8fZCocZW1dGf/15iabzmgKxvqGIADHzoxxy6WZ67yma8FwV9WIiNwNPEN8KuQDqnpIRL4I7FLVdhG5HngMqAY+IiJ/oqpXZLTl5zE4McPETJSW6szfTH1oe+fs45oyH2PTEf6o/RBf/ujGjF87V3QPTfJPPzrOt3d1EQrH2LCygne01S/LzWxjcklNmY/aMh+vD4xwm71IAAANL0lEQVRnuylAmmPuqroN2Dbn2D0pj3cSH67Juq6hKeBnX5GWy9q6MqpKi9lTIKtVO/rH+fsXOmjf1wvAz1+zilXVJVYewBS0dQ1Bjp3Jo3DPJ53LOA0yVbIE8POv9tM3MkVTBssNZ9Oh3hH+xyOvcKh3FG+RsGVtDW9bX2dVMY0B2hqCPPlKX07MnHNduCfnuDdXL3+4XttazXOv9vPdPT18+l3rl/36mbTjxCB//0IHLxwdwO/1cMul9bx1fR1BW3hkzKy2hiAjU2HOjs9Qn+WV1q77m9k1OEld0JeV1Y41ZT7W1pXxyO5ufvOWdXk/MyQWU557tZ+v/+h1dp4corbMxx+8/y0EvEWU+ArvprExC0luuXmsf8zC3WldQ5M0L8PN1PO58ZJavrWjk3996ST/9ea1WWvHUoQiUZ7Y18s//eg4x/rHqSop5sMbm9i8ugaf1+anG3M+bY3xGTOv94/z1nV1WW2L+8J9cIpNLVVZu/6VKyu45dJ6/uKZo7zvihWsqsqfsfeRqTAPbj/FN188Sf9YiMubKmYXHuX7txBjlkNDuZ9yv5dj/dm/qeqqblgkGqNneGpZFzDNJSL82c9diSrc8/hBVDVrbUlX1+Akd/3Ty1z/5/+Pr3z/KBUlxXzyrWv4xA2ttvDImEUQEdY35saMGVf13PtGponGdNmnQc7VUlPK773vUv7sqSNsO3CaD21symp75qOq7Okc4hs/OcH3D54GYGNzFW9bX8fKPPq2YUyuaWsI8vzRgWw3w13hPlsNMotj7kmffOsantjXyx+1H+KKlRU5UWsC4lsQPn2wjwd+cmK2OuNv3LKOikAxlRneucqYQrC+Ich3dnUzPDmT1SnCrhqW6RrKzhz3+XiLPPzlx68mpsrH/vElDveOZrU9/aPT/NWzr3Hdnz7LZx/eR8/wNLdfvZLfvfUttFSXWrAb45C2xIyZjiyPu7uq5945OEmRR2iqzO4qydSyBL9y02r++cWT/Kd/eJF//7Ub2LymZtnaEYspL75+lod3dPHModNEYspbGsu58ZJa2hqDeKzmizGOS9aY6egfX9a/73O5Kty7BqdYWRXIqXKyDeUBfuMdl/DAiyf4pfu38+tvv4Stt1ySsS0AVZXXzozz1IE+/u2lkwxNhin1FXHD2hpuvKSW2qBtYWdMJq2qKiFQ7Mn6jBlXhXvn4GTWb6bOp6rUx9Z3rONgzwhfe76DB7ef4tPvWs9Hr22m2oENRaZmouw6NciLHed49vBpXh+YwCPxejfvu2IFVzRV5NQ/eMa4mccjrG8I2rCMk7qHJnnv5Y3Zbsa8gn4vN15SS0t1Kd8/1MefPXWEP3/qCNevreHWyxu5YlUF6+qDNJT7z1uTIhyNcWZ0mt7haV47M8aRvlEO943yStcIUVU8Aqtry7j96pVcsbKC8mXYINwY82br64PsPJndIoKuCfeJUISz4zM5cTP1QlZVl/CrN6+ld3iaw30j9I1M8+fbjsy+XuYrorKkmBJffIl/JKpMzkSZnIkyOBEiljJtvjzg5fKmCt66vpZ19UFW15bafqTG5IC2xnIe39fLRCiStY3fXRPu3YlSv7ke7hBf6LCquoRVieJmo1Nh+sdCDIyHODceYjocYyYaIxyJUeQRasp8NFZ4uGpVJVWl8SmL9eV+qkqKs155zhjzZsmbqq+eHuO61dnZetM14X6gZwSA1XkQ7nNVlBRTUVI8+38IY0x+27KmhpLiIh7a3pm1cHfFXbZoTPn6D1/n0sYgV62qzHZzjDEFrrrMx51bWnhiXw89w1NZaYMrwn3bgT6O9Y/z2+9uw2N1UIwxOeBTb78EgPt/fDwr18/7cI/FlL977hjrG4J88Krcq+FijClMq6pKuGPTKh7e0cXgxMyyXz+tcBeR20TkqIh0iMjn53ndLyLfTry+XUTWON3Q83n64GleOzPOb797vVUvNMbklN+85RKmwlH+5acnl/3aC4a7iBQB9wIfADYAvygiG+ac9mvAkKquB/4K+LLTDZ1PLKb87Q+Osa6+jA9vXLkclzTGmLS1NZZz64ZGvvnTkwxPLm/vPZ2e+xagQ1WPq+oM8DBwx5xz7gD+JfH4EeA9ksE5el2Dk3wzsZz/6JkxPvOeNuu1G2Ny0n975zrGpsO89UvP8flHX2F/1/Cy7POQTrivArpSnncnjs17jqpGgBGg1okGznXv8x28/SvP88ffO0z/2DSfe2+b9dqNMTnrmtZqHv/0zXx4YxNP7Ovljntf5Bs/OZHx6y7rPHcR2QpsTTwdF5GjS3m/U8BzwO8sfGodcHYp18px9vnyn9s/o6s/3y8t8vxf/zL8+sVfbnU6J6UT7j1AS8rz5sSx+c7pFhEvUAmcm/tGqnofcF86DXOSiOxS1c3Lfd3lYp8v/7n9M7r98+WidIZldgJtIrJWRHzAnUD7nHPagf+SePwx4DnNh81DjTHGpRbsuatqRETuBp4BioAHVPWQiHwR2KWq7cA3gH8TkQ5gkPg/AMYYY7IkrTF3Vd0GbJtz7J6Ux9PAf3a2aY5a9qGgZWafL/+5/TO6/fPlHLHRE2OMcZ+8Lz9gjDHmzVwd7guVTch3IvKAiPSLyMFstyUTRKRFRJ4XkcMickhEPpvtNjlJRAIiskNE9ic+359ku02ZICJFIrJXRJ7MdlsKiWvDPc2yCfnum8Bt2W5EBkWA31PVDcCNwKdd9t8wBLxbVa8GNgG3iciNWW5TJnwWOLLgWcZRrg130iubkNdU9UfEZye5kqr2qeqexOMx4gExd3V03tK45C7KxYkfV90EE5Fm4EPA/dluS6Fxc7inUzbB5IlEpdFrgO3ZbYmzEkMW+4B+4FlVddXnA/4a+O9ALNsNKTRuDnfjEiISBB4FPqeqo9luj5NUNaqqm4iv/N4iIldmu01OEZEPA/2qujvbbSlEbg73dMommBwnIsXEg/1BVf1uttuTKao6DDyPu+6h3AzcLiIniQ+LvltE/j27TSocbg73dMommByWKBv9DeCIqn412+1xmojUi0hV4nEJcCvwanZb5RxV/YKqNqvqGuJ//55T1U9kuVkFw7Xhnig9nCybcAT4jqoeym6rnCUi3wJeAt4iIt0i8mvZbpPDbgZ+mXiPb1/i54PZbpSDmoDnReQV4p2RZ1XVpgsaR9gKVWOMcSHX9tyNMaaQWbgbY4wLWbgbY4wLWbgbY4wLWbgbY4wLWbgbY4wLWbibZSMi0cRc9UOJMre/JyKexGubReRvL/C7a0TkruVr7ZuuPZWoAZMTROQXEqWsbV68mZeFu1lOU6q6SVWvIL4a8wPAHwGo6i5V/cwFfncNkJVwT3g9UQMmbYmy0xmhqt8GPpWp9zf5z8LdZIWq9gNbgbsl7p3JXqiI3JKyInWviJQDXwLenjj2O4ne9I9FZE/i562J332niLwgIo+IyKsi8mCijAEicr2I/DTxrWGHiJQnqjL+hYjsFJFXROQ30mm/iDwuIrsT30K2phwfF5G/FJH9wE3nueYVicf7EtdsS/zuJ1KOfz35j0Ni05k9iff4gYP/GYybqar92M+y/ADj8xwbBhqBdwJPJo59D7g58ThIfCP32dcTx0uBQOJxG7Ar8fidwAjxQnEe4uUZ3gb4gOPA9YnzKhLvuxX4w8QxP7ALWDunjWuAg3OO1ST+twQ4CNQmnivw8cTj813z74BfSjmnBLg88bmLE8f/HvgVoJ546eq1qddN+axPzvdnbT/2413kvwXGLIcXga+KyIPAd1W1O9H5TlUMfE1ENgFR4NKU13aoajdAYpx8DfHA71PVnQCaKB0sIu8DNorIxxK/W0n8H4sTC7TxMyLy84nHLYnfOZdoy6OJ4285zzVfAv5nYiOL76rqMRF5D3AdsDPxWUuI13i/EfiRqp5IvIdrN2cxzrJwN1kjIpcQD8N+4j1XAFT1SyLyFPBB4EURef88v/47wBngauI99OmU10Ipj6Nc+P/nAvy2qj6ziHa/E3gvcJOqTorIC0Ag8fK0qkYv9Puq+pCIbCe+Q9G2xFCQAP+iql+Yc62PpNsuY1LZmLvJChGpB/4R+Jqq6pzX1qnqAVX9MvFqiZcBY0B5ymmVxHvFMeKVIxe6eXkUaBKR6xPXKBcRL/Gqob+VqBuPiFwqImULvFclMJQI9suI967TvmbiH7Xjqvq3wBPARuAHwMdEpCFxbo2IrAZeBt4hImuTxxdomzGA9dzN8ipJDJMUE9/8+t+A+eq0f05E3kV8a7ZDwNOJx9HEjcpvEh+TflREfgX4PjBxoQur6oyI/ALwd4na6VPEe9/3Ex+22ZO48ToA/NwCn+P7wG+KyBHiAf7yIq/5ceCXRSQMnAb+t6oOisgfAv83MT00DHxaVV9O3LD9buJ4P/GZRsZckJX8NWYBEt+/9UlVzakt8BLDQ7+vqh/OdltM7rFhGWMWFgUqc20RE/FvL0PZbovJTdZzN8YYF7KeuzHGuJCFuzHGuJCFuzHGuJCFuzHGuJCFuzHGuND/B1ifXAeXf2q/AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nb_merge_dist_plot(\n", " SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n", " SkyCoord(vircam['vircam_ra'], vircam['vircam_dec'])\n", ")" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Given the graph above, we use 0.8 arc-second radius\n", "master_catalogue = merge_catalogues(master_catalogue, vircam, \"vircam_ra\", \"vircam_dec\", radius=0.8*u.arcsec)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(np.sum(master_catalogue['flag_merged'])/len(master_catalogue)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add IRAC" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEKCAYAAAD+XoUoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8nOWV6PHf0WhGXbIkS3K3ZCNjm2ZA2BBqKAFCgBSSEJJAdkkcdsOSbNolu/mQu2TZS5K7JJcAm3gD2TRiahKHEkMwHVxk3HCXq2RsS1axuqad+8e8IwZFWCNrpHlHc76fjz6eecu8z1jSmUfned7ziKpijDEmPWQkuwHGGGPGjgV9Y4xJIxb0jTEmjVjQN8aYNGJB3xhj0ogFfWOMSSMW9I0xJo1Y0DfGmDRiQd8YY9JIZrIbMNDEiRO1srIy2c0wxpiUsnbt2iOqWjbUca4L+pWVldTW1ia7GcYYk1JEZF88x1l6xxhj0ogFfWOMSSMW9I0xJo1Y0DfGmDRiQd8YY9KIBX1jjEkjFvSNMSaNWNA3xpg0YkHfGGPSSFx35IrIFcD/AzzAL1T17vc57hPA48BZqlrrbPsOcDMQAm5T1eWJaPhoenjV/kG337Boxhi3xBhjEmvIoC8iHuB+4DKgAVgjIstUdcuA4wqArwKrYrbNB64HTgKmAH8VkTmqGkrcWzDGGBOveNI7C4E6Vd2tqn5gKXDtIMd9H/gB0Buz7Vpgqar2qeoeoM55PWOMMUkQT9CfCtTHPG9wtvUTkTOA6ar69HDPNcYYM3ZGPJArIhnAPcA3RvAai0WkVkRqm5qaRtokY4wx7yOeoH8AmB7zfJqzLaoAOBl4SUT2AmcDy0SkJo5zAVDVJapao6o1ZWVDloM2xhhznOIJ+muAahGpEhEfkYHZZdGdqnpUVSeqaqWqVgIrgWuc2TvLgOtFJEtEqoBqYHXC34Uxxpi4DDl7R1WDInIrsJzIlM2HVHWziNwJ1KrqsmOcu1lEHgW2AEHgKzZzxxhjkieuefqq+gzwzIBtd7zPsRcNeH4XcNdxts8YY0wC2R25xhiTRizoG2NMGrGgb4wxacSCvjHGpBEL+sYYk0Ys6BtjTBqxoG+MMWnEgr4xxqQRC/pDaO8JsKupM9nNMMaYhIjrjtx0tmJ7I6v3tPCZhbZqljEm9VlPfwiN7ZE1YR6rreet/a1Jbo0xxoyMBf0hNHX0MX9yIYU5Xhb/upb6lu5kN8kYY46bBf1j6O4L0uUPUVmay43nzMQfDPPFX9USDmuym2aMMcfFgv4xNHX2AVBWkEV5QTbf+NCJbD/cwYG2niS3zBhjjo8F/WNo7IgG/WwATppSCMDOxo6ktckYY0bCZu8cQ1NHH5kZwoRcLwAb6o8C8MiaBg4djXwg3LDIZvUYY1KH9fSPoamjj4n5WWSIAJDj81CYndk/o8cYY1JNXEFfRK4Qke0iUicitw+y/xYR2SQi60XkNRGZ72yvFJEeZ/t6EflZot/AaGrq7KOsIOs92yoKszncYUHfGJOahgz6IuIB7geuBOYDn4kG9RgPq+opqroA+CFwT8y+Xaq6wPm6JVENH22BUJjWLv+gQb+xvY+w2gweY0zqiaenvxCoU9XdquoHlgLXxh6gqu0xT/OAlI+IRzr7UKB8QNAvL8giGFZau/zJaZgxxoxAPEF/KlAf87zB2fYeIvIVEdlFpKd/W8yuKhFZJyIvi8j5g11ARBaLSK2I1DY1NQ2j+aOnqePd6ZqxKgojM3kOW17fGJOCEjZ7R1XvB+4XkRuA7wI3AQeBGaraLCJnAn8UkZMG/GWAqi4BlgDU1NS44q+Epo4+BJiYP6CnXxh5frijj4E5LmNMenp41f5Bt7txdl88Pf0DwPSY59Ocbe9nKfBRAFXtU9Vm5/FaYBcw5/iaOraaOvuYkOvF63nvf1FWpofiXK/19I0xKSmeoL8GqBaRKhHxAdcDy2IPEJHqmKdXATud7WXOQDAiMguoBnYnouGjranjb2fuRJUXRAZzjTEm1QyZ3lHVoIjcCiwHPMBDqrpZRO4EalV1GXCriFwKBIBWIqkdgAuAO0UkAISBW1S1ZTTeSCKFVWnq6GN2Wf6g+ysKs6lr6iRkNXiMMSkmrpy+qj4DPDNg2x0xj7/6Puc9ATwxkgYmQ1t3gGBYKcsfvKdfUZhFKKw0d1pv3xiTWqwMwyDeb+ZOVP8Mng4L+sakm/cbtE0VFvQHEVtdczBlBVkINm3TmHTR1u3nj+sOsOVgO6/XNdPeE+DGcyqZWpyT7KYNmwX9QTR19JHr85CXNfh/j9eTQUmez2rwGDPOBUJhfrdyHz95YSdt3QEm5mcxIddLS5efZzcf5IvnzUp2E4fNgv4gOnoDFOV4j3lMRWE2h20GjzHj1r8/vYU/rXuHps4+TijL5/Nnz2RyUaRn/3rdEZ7edJC6xk5OKB98wodbWZXNQXT7Q+T4PMc8pqIwi+auPvqCoTFqlTFmLPQGQnz/qS08+OoeQqrcePZM/u7cyv6AD7CoqoQJOV6e23IITbE6XBb0B9HtD5HrO/YfQWUFWYQV6ltsFS1jxou3Dxzlqntf5cHX9rBoVgm3XVzN3MmFiFNePSrTk8El88ppaO1h8zvt7/Nq7mRBfxDd/iC5Q/T0S3J9ALZQujHjgKry8Kr9fPy/3qCrL8Rvb17ENadNxZf5/iFywfRiyvKzeH7r4ZS6Z8eC/gDhsNLjDw0d9J05/Puau8aiWcaYUdLjD/HxB97gX/6wiZkludx8XhX74+jMeTKES+dX0NTRx87DqbOEqg3kDtDeG0BhyPROns+Dz5PBfkvvGJOy6lu6+dKva9l+qIOL55Zz8dzy/pXy4lHtDOI2dfYxd7QamWAW9Ado7Q4AkaB+LCJCSZ4vrh6BMcZ9Xt3ZxD/9fh3hsHLjOZWcOKlg2K+R7fWQ4/XQkkLra1h6Z4DW7sg3b6j0DkBxno/9LZbeMSaVhMPKAy/VcdNDq6koyGbZrecdV8CPKs7z9seNVGA9/QHa+oP+0P81pXk+ave1oKp/M7pvjHGfJa/s5vG19ew43MkpU4v4+BlTeWNX84heszjXl1JVdy3oD9DaFUnvxNvT7w2Eaeroo9ypx2OMcaeVu5u5b8VOuv0hrl0whYWVJQnprBXn+th+qCNlOn8W9AdoHUZPPzptc39LtwV9Y1yqqy/ID/+yjV+9uY/SPB+3XFjJlAmJq5lTnOcjGFY6+oIUZh/7Tn43sKA/QGu3HwGyvEMPd5TmvRv0aypLRrllxpjhem3nEW5/ciMH2nr4wgcqqSzNO+bc++NRkhsJ9K1dfgv6qai1O0CuzxPXtK0JuV5EYF+zzeAxxk2WvLKbZzcdZF19G6V5Pr503iwqJ+aNyrWKnb/4W7v9zCwdnWskUlwfeSJyhYhsF5E6Ebl9kP23iMgmEVkvIq+JyPyYfd9xztsuIpcnsvGjoa3bH1dqByK3Yk8uzLa7co1xCVXl8bUN/Pj5HWxsOMoHTyzjtkuqRy3gQyS9A9DijAe63ZDRzVnj9n7gMqABWCMiy1R1S8xhD6vqz5zjrwHuAa5wgv/1wEnAFOCvIjJHVV1bpay1KxDXIG7U9JJcm6tvjAscbu/lO09uYsW2RmaW5PLR06f2L3g0mryeDAqyMvtn/rldPF3ahUCdqu4GEJGlwLVAf9BX1diKQ3lAtBDFtcBSVe0D9ohInfN6byag7aOitds/rKA/szSXF7c3jWKLjDFD+dP6A9zxp830BUN87+r5eD0Zw7qzdqSK83y0jKOgPxWoj3neACwaeJCIfAX4OuADLo45d+WAc6ceV0vHSFt3gKnDGNmfUZJLU0cfPXGUYzbGJFZXX5DPP7iKt/a3MaMkl+vOqCIrc+x/D4tzvSnzF3/ChrFV9X5VnQ38L+C7wzlXRBaLSK2I1DY1Ja/XrKrD7ulPL8kFoL41Nb7hxowXW95p5+r7XmPd/jYunlvO4gtmMfF9ljgdbcW5Po72BFKi2mY8Qf8AMD3m+TRn2/tZCnx0OOeq6hJVrVHVmrKysjiaNDp6AiH6gmFy32eZxMFER+ttBo8xYyNaBvmjD7xOZ2+Qm8+r4tJ5FWOazhmoOM9HWKG9x/2DufEE/TVAtYhUiYiPyMDsstgDRKQ65ulVwE7n8TLgehHJEpEqoBpYPfJmj45osbXh9PRnOD39VPnTzphU1tUX5GuPrOdf/rCJs2eV8uxXz2dWWfKXK4xO20yFvP6QXVpVDYrIrcBywAM8pKqbReROoFZVlwG3isilQABoBW5yzt0sIo8SGfQNAl9x98yd+IutRRXnesnPyrRpm8aMsv98bjtLV9dzpLOPy+ZXcOGcMpZvPpzsZgFQ4kzbbO3yQ/KSFXGJK4+hqs8AzwzYdkfM468e49y7gLuOt4Fjqa2/px9/ekdEmFGSa4upGDNKVJVfv7mP/3ppFzk+D39/XhWzXdC7j1WU40UgJapt2h25MYZTVjnWjJJcdjamzso5xqSK1i4//+uJjTy35TAnVhTwiTOnkT+MMbex4skQinK9/SliN3Pf/14StR1v0C/NZcX2RsJhJSPD/VX2jEkFr+08wjceW09Ll5/vXjWPbG985VGSpTjX158idjML+jGin9LDnW8/vSQXfzDM4Y5eJhclrnqfMemoLxjii/9Ty6t1RyjLz+LLF8weVso1WUpyfSnxF7/7/yfHUEuXn4KsTDIzhnf7QmVpZAbPvuZuC/rGjMD2Qx187ZH1bD3YzqKqEq48eXLCq2KOluI8L+29QQKhMF6Pe9vs3pYlQVu3nwl5wy+NWtk/V98Gc405HuGw8tBre7j6vtdo6ujlxrNncu2CqSkT8OHdaZttLs/rW08/Rmt3oP8bNxyTi7LxeoS9doOWMcP2TlsP33p8A6/XNXPpvHLu/sSpPOeSqZjDEVtiuSxJdwbHw4J+jLZuPxOOI+hnejKYXmzTNo0ZDlXl249v5M8b3yEcho8umMpZlcUpGfDh3RLLbp+2aUE/Rmt34Ljrbs8szWXvEevpGxOP1i4/3/3j2zy96SAzS3L5ZM30/hucUlV01l+P37X3nwIW9Hl41f7+x40dvce9qv3M0jxW72lJmcWRjUmWV3Y08c3HNtDa7efy+RWcP6fM1VMx4+X1ZJCZIfQELOinhFBY6Q2Eyc06vrKslaW5dPlDHOl0dz7PmGTp8Ye4+9mt/OrNfVSX5/PQF85iY8PRZDcroXJ8Huvpp4pufxAYXgmGWDMnvjuDx4K+Me/1w79s49Haeo50+jl3dikfOmnSuAv4ANlej+t7+qkzH2qUdTufzsO9GzcqOm3TZvAY8y5/MMw9z23nZy/vIhBSbj6viqtOneLqeewjkZMCQd96+o7jDfrRMYFQWMkQeHrjO/iDYW5YNCPhbTQmlbx94CjffGwD2w51cPr0CXzk1CnjfnW5HK+Hjj6bp58SekaY3vFkCBNyfTSnQO0NY0ZTbyDEfSvq+K+XdzEx38eDN9Vw+DgnSKSaHJ+Hxo7eZDfjmCzoO0aa3gEozfPR3GlB36SvtftaueU3a2nq7OOMGcVcdcrktAn4YOmdlBIN+nkjKOxUmu+jvr4NVfevk2lMIvUGQvxo+XYeen0PRdlevvCBSuZUFCS7WWMux+ehLxAmrOraaahxRTgRuQL4f0RWzvqFqt49YP/XgS8SWR2rCfh7Vd3n7AsBm5xD96vqNQlqe0J1+4NkZghez/F/o0rzsugNhPs/QIxJB+vr2/j6o+vZ3dTFZxfN4ISyfLK84zt3/36yvR4U6AuEXTt+MWTQFxEPcD9wGdAArBGRZaq6JeawdUCNqnaLyD8APwQ+7ezrUdUFCW53wnX5Q+T6PCO6sarUuaPQ8vomHQRDYe59YSf3vVjHpMJsfnvzIs6rnvieGx7TTY7zYdcTCKVu0AcWAnWquhtARJYC1xJZ9xYAVX0x5viVwOcS2cix0O0Pjbhmd0m+E/Q70yeHadLTgbYebliykn0t3ZwxIzIzZ39Ld1oHfIgJ+v4QHF9Fl1EXT5SbCtTHPG8AFh3j+JuBZ2OeZ4tILZHUz92q+sdht3IM9Cbgk7kk14dgPX0zvi3ffIhvPbaBvmCYT9dM57TpE5LdJNeIxhA3D+YmdCBXRD4H1AAXxmyeqaoHRGQWsEJENqnqrgHnLQYWA8yYkZz57b2BEBNyhl9LP1amJ4OiXC8tFvTNOBQMhfnh8u0seWU3p04r4rJ5FZTm293nsWLTO24VT9A/AEyPeT7N2fYeInIp8K/Aharan99Q1QPOv7tF5CXgdOA9QV9VlwBLAGpqapIy9aUnEGJSYfaIXycybdPSO2Z8+dnLu1i6up69zV0sqirhqlMmkzlO76odiWxv5P+k18WTOeL5rq0BqkWkSkR8wPXAstgDROR04OfANaraGLO9WESynMcTgXOJGQtwk95AiOwEDLyU5mVZeseMK2/sOsL9K+o40NbNJ8+cxrULplrAfx/jIr2jqkERuRVYTmTK5kOqullE7gRqVXUZ8CMgH3jMmf0SnZo5D/i5iISJfMDcPWDWjyuEVSNTrBIwzaw030e3P8TR7gBFuSNLFxmTTOGw8sBLddzz/A5K8rL4u3OrmFQ08r+GxzOfJ4MMSfGgD6CqzwDPDNh2R8zjS9/nvDeAU0bSwLHgD4ZRInNsR6o0L5Lj3NvcxWm5NsBlUlNjey/ffHwjr+xo4toFU1gwbULazr0fDhGJ3JWb4umdcS/6qZydgEWYJxZEpm3uauoc8WsZkwzPbjrI5T95hdV7mrnrYyfzk08vsIA/DDk+d5disDIMRPL5kLievkeEHYct6JvUcrQ7wBd+uZp19W1MK87hCx+oQhB+v7p+6JNNvxyvpz+muJEFfd7t6SfiDjpPhjCxwEddY8eIX8uYsfLXLYf5lz9s4khnHxfPLeeDJ5bjyXBn7Ri3c/tCKhb0gV5/GEhMTx+gvCDbevomJbR1+/nfyzbzx/XvMHdSAZ+smc7UCTnJblZKy/F5XH2vjgV93k3vJGL2DkB5YRZvv3OUHr9762+Y9Pbwqv1sO9TOH9YdoKsvyMVzy7noxDIyM2yYb6TcXl7Zgj7QG4zm9BPzA19ekI1qZDD35KlFCXlNYxKlvTfAE2sbWLu/lUmF2dx0TiVTrHefMNGcvltLrFvQ592cflZmYnrlFc7C6DsbOyzoG1d5o+4I33xsAweP9nLRnDIunltuN1olWI7PQ1gjU8HdyII+kVumfZkZCRu4Ks3PwusRdlpe37hEjz/ED/6yjf95Yy+zJuZxy4WzmV6Sm+xmjUvZLq+/Y0Ef6E3Q3bhRngyhamKeDeYaV7j72W08VltPc5efc2aXcvn8SfgScE+KGZzbi65Z0CfyzUlUPj+quryAze8cTehrGjMc/mCYn/x1Bz9/eRdFOV5uPq+K2WX5yW7WuNdff8eld+Va0CcykJuo6ZpRJ5Tn8+zbByOF3OxuRjPG6ho7uO3369lysJ2amcV8+JTJ9nM4RqynnwJ6AyEKsxNbHG1ORQFhZwbPSVNsMNeMDVXlNyv3cdfTW8nLyuS/b6yhqcNKfY+l96ye5UIW9Ink9MsLEtsLqq6I/Bld12hB34yNB1/dw+NvNbD1YDtzKvL5xBnTLOAngQ3kpoAef+Jz+pWleWRmCDsOWzkGM/o21Ldx34s7OdoT4MOnTObc2aU4Zc7NGMvyZiDg2vo7aR/0VXVU8u6+zAwqJ+bZtE0zqlSVX7+5j39/egt5vky+fIFNxUy2DBFX199J+6AfraWfyCmbUdXl+Ww/ZD19Mzo6+4Lc/sRGntp4kIvnlvOB2aXk+tL+V9oVcnzuramf9j8hvcHEFluLVV1RwPLNh2wGj0moh1ft59DRXh5evZ/mzj4un1/B+XPKyLB0jmu4uf5OXIlsEblCRLaLSJ2I3D7I/q+LyBYR2SgiL4jIzJh9N4nITufrpkQ2PhF6ElhLf6Dq8nzCCnuOdCX8tU16UlXe3HWEB16qozcQ4ubzqrjwxHIL+C7j5tWzhgz6IuIB7geuBOYDnxGR+QMOWwfUqOqpwOPAD51zS4DvAYuAhcD3RKQ4cc0fueiq9YkeyAU4cVIBANsOtSf8tU36ae7s44u/quXPGw8yuyyf2y6pZpbdbOVK2d4MegLurL0TT6RbCNSp6m5V9QNLgWtjD1DVF1W123m6EpjmPL4ceF5VW1S1FXgeuCIxTU+MRJdVjjVrYh45Xg+bGizom5F5YethLv/Jq7y68wgfOXUyN54zk/ystM/OulaOz72rZ8XzUzMViF0vrYFIz/393Aw8e4xzpw48QUQWA4sBZsyYEUeTEmc00zuZngzmTylk04G2hL+2SQ+/fG0PT286SO2+SBnkWy6czaSi7GQ3ywwhmtNXVddNnU1oV0FEPgfUABcO5zxVXQIsAaipqRnTItSjOZALcMrUIh6trScUVlt+zgzL2n2t/PTFOlq7/Fw4p4xLrAxyysjxegiFNVLM0WULKcUT9A8A02OeT3O2vYeIXAr8K3ChqvbFnHvRgHNfOp6GjpZ3F0VP7C/Tw6v2A5Fpdd3+EPe+sJOKwmxuWDS2f8mY1BMMhfnpijrue7GOwuxMvnT+LCon5iW7WWYYsp1Af7QnkJJBfw1QLSJVRIL49cANsQeIyOnAz4ErVLUxZtdy4D9iBm8/BHxnxK1OoF5/CK9HRm2ZuOh6owfaeqgotD/LzbHtburk649uYH19Gx8/YyonTymy6b4pKDpGeLQn4Lp03JCRTlWDwK1EAvhW4FFV3Swid4rINc5hPwLygcdEZL2ILHPObQG+T+SDYw1wp7PNNXpGeQ59WUEWPk8GB9p6Ru0aJvWpKr95cy8fvvdV9hzp4r4bTueeTy2wgJ+iYoO+28SV01fVZ4BnBmy7I+bxpcc49yHgoeNt4Ggb7RunMkSYPCGbA60W9M3gDh7t4cYHV7OzsZPq8kihtPaeYH+K0KSeaEqnPVWD/niW6FWzBjNtQg6r97YQCrtzoWSTHKrKY2sb+P6ft9AbDHHNaVNYVFXiutkeZvhSvqc/nvUGQ+SO8kDL1OIcAruUpk4rc2siGjt6uf2JTazY1sjCqhLOP2EipflZyW6WSRAL+i7W4w9Rkucb1WtMiQ7mWorHAP/y5Cb+uP4A/mCYq06ZzDmzS62MwjiTZUHfvcaiGNrE/Cx8mRkcaOse+mAzbrX3Bvjfyzbz5FsHmDohh0/WTKO8wF0zO0xieDIEX2YGHb3BZDflb6R10I/U0h/9nH6GCFMn5FhPP42t3tPCPz+ynkPtvVw8t5wPnlhuN+uNczleD+297uvpp/Xtfb2BMCHVMZkWN3VCDgeP9hIIubMIkxkd/mCYH/xlG59e8iaZHuGxW87h0nkVFvDTQLY3gw4XBv207ulHP4VHo8LmQFMn5BAMKzsPdzJ/SuGoX88k34+f38GjtfUcPNpLzcxirjp1MtsO2qI66SI700N7j6V3XCX6KTza6R2IzOAB2NDQZkF/nAuHlV++sZf7X6wjKzODz589k3mT7XuebrJdmt5J66B/1PkUHov0Tmmej/ysTFbvaeEzC63+znjV2N7LNx7bwKs7jzB3UgEfO30qBdneZDfLJEGOz0NLlz/ZzfgbaR30303vjH7QFxGqJuaxcnezK8utmpF7YethvvX4Rrr9Qe762Mmg2Pc5jWVlZlhP322it0iPRU4foGpiHpsOHKW+pYcZpbljck0z+n795l6Wv32I13c1M7kom8+fPRtBwOJ9WsvxeujoDbquk5feQd+ZQzsWOX2IBH2AlbubLeiPE/ubu/n5y7s50NbD2bNKufLkSXit5r0hkkEIhZVuf4g8F61yltY/nR1jmN4BKC/IoiTPx8o9zWNyPTO6/rzhHa766as0d/Xx2UUzuOa0KRbwTb9oXHFbisc9Hz9J0N4TxJMhY/aLKiIsqiph1W5XVZc2w9TZF+R7f9rME281cPqMCVw6t4LiUS7lYVJPNG3c0RtkclGSGxMjrbsl7b2BMa9XfvasUg609VDfYiUZUtHafa1cde+r/GFdA7ddUs1jXz7HAr4ZVDRt7Lbyymne0w+QM0aDuFGLZpUAsGpPC9NLLK+fKnoDIX78/A6WvLKbolwvXzxvFpMKs3m0tiHZTTMu5db0TlwRT0SuEJHtIlInIrcPsv8CEXlLRIIict2AfSFnNa3+FbXcor03OOY9/TnlBRTnelm52/L6qWJDfRsf+elr/PyV3dRUlvDVi6ttzVozpGhscVvRtSF7+iLiAe4HLgMagDUiskxVt8Qcth/4AvDNQV6iR1UXJKCtCdfeM/bpnYwMYWFVCatsMNf1/MEwP12xkwde2kVZfha/+vuFVjTPxC2a00/F9M5CoE5VdwOIyFLgWqA/6KvqXmdfSlUTO9oToCB77DNcZ88qZfnmwxxo6+lfON24yz3P7eCxtZG6OWfMmMBVp0yxgG+G5d30jrt6+vGkd6YC9THPG5xt8coWkVoRWSkiHx1W60ZZa7d/zObox1pUVQrAyl3W23ebcFj5xau7uf+lOtp7g3xu0UyuO3N6/5qnxsTL68nA58K7cseimztTVQ+IyCxghYhsUtVdsQeIyGJgMcCMGWNTlyYUVo72BMj1jX1Pf+6kAibm+3hpRxOfOHPamF/fDO7g0R6+8egG3tjVzLzJhXzs9Knku+imGpN6CrO9rqu0Gc9P9AFgeszzac62uKjqAeff3SLyEnA6sGvAMUuAJQA1NTVjsnp4R28AVUZ9fdyBHl61H4DK0jye23yIX7+5lxvPqRzTNpi/9eymg9z+5CYCoTA/+MQpBEPuunXepKbC7EzX9fTjSe+sAapFpEpEfMD1QFyzcESkWESynMcTgXOJGQtIptbuyDdirIN+1LzJhfQFw+w9YvP1k6mrL8gn/usN/uF3b1GQnck/XDibUNgKpZnEKMjxpt7sHVUNisitwHLAAzykqptF5E6gVlWXichZwB+AYuBqEfk3VT0JmAf83BngzQDuHjDrJ2nauiMlT5MV9GeX5eP1CFsOtifl+iYyFfOrS9exr7mbi04s45K5tqKVSazC7MyUnL2Dqj4DPDNg2x0xj9cQSfsMPO8N4JQRtnFUtPX39JOTs/VlZnCGt2fJAAAVsklEQVRCWT7bDra7rgrfeBcKKz97eRc/fn4H5QVZfPH8Wf3F8IxJpMJsL++0uWvWV9qWYWh1evrJnJUxb3IhbT0BttoSemOmvqWbzyxZyY+Wb+eKkyfx7NcusIBvRk1hTqbrpmym7dSEtiTn9AFOnFSAAM9vOWxLKI4yVeXbj29k2YZ3ALjuzGmcPn0CT288mOSWmfEsMnvHXemdtO3pt3X7ERm7ssqDKcj2Mr0kl79uPZy0NqSD1i4/tz68jsfWNjC5KJvbLq7mjBnFllIzo64gO5O+YJi+YCjZTemXtj391u4ARTleMpL8iz9vciHLNx/i4NEeJhfZ3bmJ9vKOJr712AZau/18aH4FF8wpS/r33KSPwpzI+sgdvUGy8t1xg1/69vR7AhTnJr8k7rxJBQD8dYv19hOp2x/kjj+9zU0PraYox8sf/vFcLjqx3AK+GVOF2ZGg76YUT9r29Nu6/RQ5n8LJVFaQxZyKfJ5cd4DP201aCXH3s9t4rLae5i4/584u5UMnTWJjw9FkN8ukoWhtLzfN1U/bnn5rt5/i3OQHfRHhUzXTWbe/je2HbBbPSPQFQ/zf5dv5+cu7CIWVm8+r4qpTbQlDkzzR9I6b7spN29+Gtm53pHcAPn7GNLwe4ZE19UMfbAa1qeEo1/z0de57sY7TZxRz2yXVzC7LT3azTJp7N71jPf2ka+sOUOSCnj5ASZ6PD500iSfXNbhqlD8V9AYivfuPPvA6bT1+HvpCDdedOS2ps7KMiXo3vWM9/aTyB8N09gVd09MHuP6s6bR1B3husw3oxmvl7mbO+8GL3PdiHadNK2Lx+bM5dLQv2c0ypp8b0ztpOZB71BlJd0NOP+rc2ROZOiGHR9bUc/VpU5LdHFc72h3g/zy7laVr6inO9fJ3H6ikuqIg2c0y5m/k+TxkiLsGctMy6EeLrRXl+uh0yTcjI0P49FnTuef5HdS3dNui6YNQVZ59+xB3/Gkzrd1+Fl8wiylFOfgy0/IPVpMCRIQCl92Vm5a/LW0u7OlDpDRAhsDDq/cnuymu87OXdvHhe1/jH3/3Fj6PcMuFs6kszbOAb1zPbfV30rKn39oV6ekX5/qob0l+BbzowioA86cU8eBreyjO9bH4gllJbJU7qCqPrW3gJy/sIBhSrjx5Eh+YPdFKIJuUUZjttYHcZIsWW3PDzVkDXTqvnEAwzCs7mpLdlKRraO3mC79cw7cf38ikwmxuu6Sa86vLLOCblFKQnemqKZtp2dNv63F6+nnumb0TVV6QzRkzilm5uzlt6/EEQ2F++fpe7nl+ByLwb9echCdDrISCSUmF2V72t7hnhby4evoicoWIbBeROhG5fZD9F4jIWyISFJHrBuy7SUR2Ol83JarhI9HaHcDrEfKSWFb5WC6eV44q3PtCXbKbMubW7G3h/B++yF3PbGVmaS63fvAEvJ4MC/gmZRW6bMnEIXv6IuIB7gcuAxqANSKybMCyh/uBLwDfHHBuCfA9oAZQYK1zbmtimn98InV3fK4trVuc62NhVQmP1tbz5QtmUZkGi3zsa+7i7me38ezbhyjMzuSGhTM4aUqha79HxsSrwGVLJsaT3lkI1KnqbgARWQpcS8wC56q619kXHnDu5cDzqtri7H8euAL4/YhbPgKREgzuy+fHuujEMtbXt/HvT2/lv288c9wGv5YuPz9dsZPfrtyH15PBNy6bQ0G212blmHGjMNtLR1+QUFhdMR4VT9CfCsQWhWkAFsX5+oOdOzXOc0dNa7efCS4P+gXZXr5+2RzuemYrv121n8+fPTPZTUqobn+Qry5dzys7mvAHw9RUFnPJvIr+WiXGjBfRu3I7+4KumDziioFcEVkMLAaYMWPGqF+vrTuQEjc/3XxeFa/VHeH7T22hZmYx8yan/pKKobDyxNoG/vP57Rxu72P+5EI+NL+C8sLsZDfNmFERrb/T3hNwRdCP52/oA8D0mOfTnG3xiOtcVV2iqjWqWlNWVhbnSx+/VEjvQOQu3f/81GkU5Xj5p9+vo9vvnsGg4VJVXtreyFX3vsq3n9jI5KIcFp8/i8+dPdMCvhnXon+9umUwN56e/hqgWkSqiATs64Eb4nz95cB/iEix8/xDwHeG3coEi6R33Dddc6DoTVtXnzqFX76+h+uXrORTNdP5XIqlev7j6a08t+UQe5u7Kc71cv1Z0zllatG4HacwJlZhjtPTd8kNWkMGfVUNisitRAK4B3hIVTeLyJ1AraouE5GzgD8AxcDVIvJvqnqSqraIyPeJfHAA3Bkd1E2WHn+IvmDY9Tn9WCeU53PZ/Aqe23IYVfj0WdNdvzCIqvLqziMseWU3r9UdoSArk2tOm0JNZTGZGe5uuzGJ5LYlE+PK6avqM8AzA7bdEfN4DZHUzWDnPgQ8NII2JlT/jVkp0NOPFV3f9S+bD/GPv3uL+244naxM991n0OMP8dTGd3jwtT1sO9RBeUEWV5w0ibNnldqMHJOWUjG9M660dkU+bSe4YEBluC6YU4Y3M4M/b3iHGx9czT2fXsDUCe64Y7eusYOHV9Xz8Op99AbClBdk8YkzpnHa9CLr2Zu0lnLpnfEm2tNPhZz+YM6ZVcrFc8v41z+8zRU/foU7rp7PdWdOS0p+vKM3wNMbD/JobT1v7W/D6xHmTS5kUVUplaW5lrM3BsjPis7esZ5+UkSLrRXnpV5PP6rHH+YfLzqBx9c28K3HN7Lkld384LpTOWNG8dAnj1BvIMSKbY3c/2Id2w91EAwr5QVZXHnyJE6fUdz/A26Micj0ZJCfldm/eFOypd1vaKuzgMqEnNTs6UeV5Pn44vlVvLmrmRXbGvn4A29w9qwSvnT+LM6rnpiwfL+qsre5m1d2NPHS9kbe3N1MbyBMflYmZ1WWsGD6BKYV51iv3phjmJjvo6nTHUt5pl3Qj/b0U2n2zvvJEOHcEyZSU1lMKKz896u7uflXteT6PJxfPZEL5pQxd1IBJ5QVxLUIvKpy8Ggv2w91sPVQO+v2t/HWvlaanfUHSvN8nD49cpPYrLI8K4JmTJzKC7M53N6b7GYAaRn0/eR4PWR73Tfz5XhFe/VfuegE6po62Xaog5W7W1ges8j6xHwfE/OzKMnzUZzrIyNDECCsSnOnn8aOXhpae+gLht9zzszSXM6rnsgJZfmU5meN9VszZlyoKMzm7QNHk90MIA2Dfmt3YFz08geT6clg7qRC5k4qRFVp7Q7Q2N5LY0cfRzr76PKHaGjtYcfhTlS1/7y8rEwKsjM5Y0YxZQVZVBRmM6kwmxyXlp42JtVUFGTxQnsvqpr0VGjaBf227kDKztwZDhGhJM9HSZ6PuZOT3Rpj0ltFYTbd/hCdfUEKklxUMO0mULd1+1Oi7o4xZvwoL4ykRg+3J38wN+2CfiqUVTbGjC/lBZGigo0uGMxNq6Cvqhxu76PCqjoaY8ZQRbSn32FBf0wd7QnQ2RdkWrH7a+kbY8aPaPnwRkvvjK2G1h4AphW7o16NMSY95Gdlkp+VaTn9sdbQ2g1Y0DfGjL3ywixL74y1d3v6lt4xxoyt8oIsG8gdaw2tPRRkZbpinUpjTHqpKMy29M5Ya2jtZqqldowxSVDh1N+JvRs+GeIK+iJyhYhsF5E6Ebl9kP1ZIvKIs3+ViFQ62ytFpEdE1jtfP0ts84enobXHUjvGmKQoL8iiLximPckraA0Z9EXEA9wPXAnMBz4jIvMHHHYz0KqqJwA/Bn4Qs2+Xqi5wvm5JULuHTVWdoG89fWPM2KsodMcNWvH09BcCdaq6W1X9wFLg2gHHXAv8ynn8OHCJJLuq0ADvztG3oG+MGXvlBe4oxRBP0J8K1Mc8b3C2DXqMqgaBo0Cps69KRNaJyMsicv5gFxCRxSJSKyK1TU1Nw3oD8bKZO8aYZIr29JNdV3+0B3IPAjNU9XTg68DDIlI48CBVXaKqNapaU1ZWNioNsTn6xphkKndJKYZ4gv4BYHrM82nOtkGPEZFMoAhoVtU+VW0GUNW1wC5gzkgbfTyiPf3p1tM3xiRBri+ybkWySzHEE/TXANUiUiUiPuB6YNmAY5YBNzmPrwNWqKqKSJkzEIyIzAKqgd2JafrwROfoF+ak3RICxhiXqCjMpjHJPf0hI6CqBkXkVmA54AEeUtXNInInUKuqy4AHgd+ISB3QQuSDAeAC4E4RCQBh4BZVbRmNNzKU6Bx9l40vG2PSSEVhVtIHcuPq9qrqM8AzA7bdEfO4F/jkIOc9ATwxwjYmhM3RN8YkW3lBNmv2JqXf2y8t7si1OfrGGDcoL8yisb0vqXflpkXQtzn6xhg3qCjIxh8K09YdSFob0iLo2xx9Y4wb9M/VT+JgbpoEfZujb4xJvuiyicmctpkmQd/m6Btjki+6QHoy78pNm6Bvc/SNMckWvSu3scN6+qPK5ugbY9wg2+uhKMfLoaPW0x9VNkffGOMW1eX5bGhoS9r1x33Q7/GH2NvcxcxSC/rGmOT74NxyNjYcTVo5hnEf9F/c3khvIMwlc8uT3RRjjOGDJ0Zi0UvbRqeM/FDGfdB/auM7TMz3sWhW6dAHG2PMKJs3uYDJRdms2NaYlOuP66Df1RdkxbZGrjx5Mp4MG8Q1xiSfiHDRieW8VncEfzA85tcf10H/hW2R1M5HTp2c7KYYY0y/i+eW09kXTErxtXEd9J/a8A7lBVmcVVmS7KYYY0y/c08oxZeZwQtbxz7FM26DfkdvgJd2NPHhUyaTYakdY4yL5PoyOWdWKS9ut6CfMH/dehh/MMzVp1lqxxjjPhfPLWfPkS72HOka0+vGFfRF5AoR2S4idSJy+yD7s0TkEWf/KhGpjNn3HWf7dhG5PHFNP7anNhxkSlE2p08vHqtLGmNM3C52ppGP9SyeIYO+s8bt/cCVwHzgMyIyf8BhNwOtqnoC8GPgB86584ksnXgScAXwQHTN3NHS7Q/y709tYcX2Rq4+bYqldowxrjS9JJe5kwq494Wd/HHdgTFbWCWenv5CoE5Vd6uqH1gKXDvgmGuBXzmPHwcukUihm2uBparap6p7gDrn9UbFG3VHuOInr/KL1/Zww8IZ3HZJ9WhdyhhjRuyBz57B7LI8vvbIer7067U0jkH1zXiC/lSgPuZ5g7Nt0GNUNQgcBUrjPDchdjd18tkHV+HJEJYuPpu7PnYKeVlWVdMY416zyvJ57JYP8N2r5vHqziZu+MUqwuHR7fG7IiqKyGJgsfO0U0S2H+9r7QXO+dawTpkIHDne66UAe3+pzd5fCvvsMI/fAXi+cdyXmxnPQfEE/QPA9Jjn05xtgx3TICKZQBHQHOe5qOoSYEk8DU40EalV1ZpkXHss2PtLbfb+TKLFk95ZA1SLSJWI+IgMzC4bcMwy4Cbn8XXACo2MSiwDrndm91QB1cDqxDTdGGPMcA3Z01fVoIjcCiwHPMBDqrpZRO4EalV1GfAg8BsRqQNaiHww4Bz3KLAFCAJfUdXQKL0XY4wxQ5CxmibkViKy2EkvjUv2/lKbvT+TaGkf9I0xJp2M2zIMxhhj/lbaBv2hSkukOhF5SEQaReTtZLcl0URkuoi8KCJbRGSziHw12W1KJBHJFpHVIrLBeX//luw2jQYR8YjIOhF5KtltSSdpGfTjLC2R6v6HSOmL8SgIfENV5wNnA18ZZ9+/PuBiVT0NWABcISJnJ7lNo+GrwNZkNyLdpGXQJ77SEilNVV8hMpNq3FHVg6r6lvO4g0jgGJU7vZNBIzqdp17na1wNvonINOAq4BfJbku6SdegP2blIczociq6ng6sSm5LEstJfawHGoHnVXVcvT/gJ8C3gbFfLzDNpWvQN+OAiOQDTwBfU9X2ZLcnkVQ1pKoLiNzFvlBETk52mxJFRD4CNKrq2mS3JR2la9CPqzyEcS8R8RIJ+L9T1SeT3Z7RoqptwIuMr/GZc4FrRGQvkdTqxSLy2+Q2KX2ka9CPp7SEcSmnbPeDwFZVvSfZ7Uk0ESkTkQnO4xzgMmBbcluVOKr6HVWdpqqVRH73Vqjq55LcrLSRlkHfKf8cLS2xFXhUVTcnt1WJJSK/B94EThSRBhG5OdltSqBzgc8T6SGud74+nOxGJdBk4EUR2Uikg/K8qtq0RpMQdkeuMcakkbTs6RtjTLqyoG+MMWnEgr4xxqQRC/rGGJNGLOgbY0wasaBvjDFpxIK+SToRCTlz7Tc75YS/ISIZzr4aEbn3GOdWisgNY9fav7l2j1MjxxVE5NNOuXCb128GZUHfuEGPqi5Q1ZOI3H16JfA9AFWtVdXbjnFuJZCUoO/Y5dTIiZtT2ntUqOojwBdH6/VN6rOgb1xFVRuBxcCtEnFRtNcqIhfG3IG7TkQKgLuB851t/+z0vl8Vkbecrw84514kIi+JyOMisk1EfueUc0BEzhKRN5y/MlaLSIFT5fJHIrJGRDaKyJfjab+I/FFE1jp/tSyO2d4pIv8pIhuAc97nmic5j9c716x2zv1czPafRz80nIWA3nJe44UEfhvMeKaq9mVfSf0COgfZ1gZUABcBTznb/gyc6zzOBzJj9zvbc4Fs53E1UOs8vgg4SqS4XgaREhXnAT5gN3CWc1yh87qLge8627KAWqBqQBsrgbcHbCtx/s0B3gZKnecKfMp5/H7X/Cnw2ZhjcoB5zvv2OtsfAG4EyoiUB6+KvW7Me31qsP9r+7KvzGF+RhiTTK8D94jI74AnVbXB6azH8gL3icgCIATMidm3WlUbAJw8fCWRD4KDqroGQJ0SzSLyIeBUEbnOObeIyIfIniHaeJuIfMx5PN05p9lpyxPO9hPf55pvAv/qLDDypKruFJFLgDOBNc57zSFSY/9s4BVV3eO8xrhcMMckngV94zoiMotIkGwk0tMFQFXvFpGngQ8Dr4vI5YOc/s/AYeA0Ij363ph9fTGPQxz751+Af1LV5cNo90XApcA5qtotIi8B2c7uXlUNHet8VX1YRFYRWVHqGSelJMCvVPU7A651dbztMiaW5fSNq4hIGfAz4D5V1QH7ZqvqJlX9AZHqk3OBDqAg5rAiIr3oMJFKnEMNmm4HJovIWc41CkQkk0gF1n9w6vYjInNEJG+I1yoCWp2AP5dIbzzuazofdrtV9V7gT8CpwAvAdSJS7hxbIiIzgZXABSJSFd0+RNuMAaynb9whx0m3eIksev4bYLA6+V8TkQ8SWWJvM/Cs8zjkDJD+D5Gc9xMiciPwF6DrWBdWVb+IfBr4qVO7vodIb/0XRNI/bzkDvk3AR4d4H38BbhGRrUQC+8phXvNTwOdFJAAcAv5DVVtE5LvAc8401gDwFVVd6QwUP+lsbyQy88mYY7LSysYcJ4msz/uUqrpqKUMnzfRNVf1Istti3MfSO8YcvxBQ5Labs4j8tdOa7LYYd7KevjHGpBHr6RtjTBqxoG+MMWnEgr4xxqQRC/rGGJNGLOgbY0wa+f/NKwOrhWVWlQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nb_merge_dist_plot(\n", " SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n", " SkyCoord(irac['irac_ra'], irac['irac_dec'])\n", ")" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Given the graph above, we use 1 arc-second radius\n", "master_catalogue = merge_catalogues(master_catalogue, irac, \"irac_ra\", \"irac_dec\", radius=1.*u.arcsec)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(np.sum(master_catalogue['flag_merged'])/len(master_catalogue)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cleaning\n", "\n", "When we merge the catalogues, astropy masks the non-existent values (e.g. when a row comes only from a catalogue and has no counterparts in the other, the columns from the latest are masked for that row). We indicate to use NaN for masked values for floats columns, False for flag columns and -1 for ID columns." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": true }, "outputs": [], "source": [ "for col in master_catalogue.colnames:\n", " #print(col)\n", " if (col.startswith(\"m_\") or col.startswith(\"merr_\") or col.startswith(\"f_\") or col.startswith(\"ferr_\")):\n", " master_catalogue[col] = master_catalogue[col].astype(float)\n", " master_catalogue[col].fill_value = np.nan\n", " elif \"flag\" in col:\n", " master_catalogue[col].fill_value = 0\n", " elif \"id\" in col:\n", " master_catalogue[col].fill_value = -1\n", " \n", "master_catalogue = master_catalogue.filled()" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=10\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
idxps1_idradecflag_mergedcfht_intidcandels_idcfhtls-wide_idcfhtls-deep_idsparcs_intidwirds_idvipers_idsxds_intidsxds_b_idsxds_v_idsxds_r_idsxds_i_idsxds_z_idhsc_intidhsc-wide_idhsc-deep_idhsc-udeep_iddecam_intiddecals_iddes_idukidss_intiddxs_iduds_idvircam_intidvhs_idvideo_idviking_idirac_intidservs_intidswire_intid
degdeg
010448034793855023934.7939347122776-2.9337278347732303False-1-10.00.0-1-1-1-10.00.00.00.00.00.0-1-1-10.0-1-10.0-1-10.0-1-1-1-1-1-1
110447034788859090934.7889409922776-2.94131771477323False-1-10.00.0-1-1-1-10.00.00.00.00.00.0-1-1-10.0-1-10.0-1-10.0-1-1-1-1-1-1
210438034869156621634.8690964322776-3.01188159477323False-1-10.00.0-1-1-1-10.00.00.00.00.00.0-1-1-10.0-1-10.0-1-10.0-1-1-1-1-1-1
310334034226446212134.2265307522776-3.88177687477323False-1-10.00.0-1-1-1-10.00.00.00.00.00.0-1-1-10.0-1-10.0-1-10.0-1-1-1-1-1-1
410215034626400630934.6264773122776-4.87011075477323False-1-10.00.0-1-1-1-10.00.00.00.00.00.0-1-1-10.0-1-10.0-1-10.0-1-1-1-1-1-1
510291036934859514736.9348618522776-4.23782747477323False-1-10.00.0-1-1-1-10.00.00.00.00.00.0-1-1-10.0-1-10.0-1-10.0-1-1-1-1-1-1
610410035464316923635.4645593722776-3.24264912477323False-1-10.00.0-1-1-1-10.00.00.00.00.00.0-1-1-10.0-1-10.0-1-10.0-1-1-1-1-1-1
710128036323894127036.3242755922776-5.59929363477323False-1-10.00.0-1-1-1-10.00.00.00.00.00.0-1-1-10.0-1-10.0-1-10.0-1-1-1-1-1-1
810442034823753762534.8238613522776-2.97741051477323False-1-10.00.0-1-1-1-10.00.00.00.00.00.0-1-1-10.0-1-10.0-1-10.0-1-1-1-1-1-1
910215034626398566734.6265266322776-4.8707222947732305False-1-10.00.0-1-1-1-10.00.00.00.00.00.0-1-1-10.0-1-10.0-1-10.0-1-1-1-1-1-1
\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "master_catalogue[:10].show_in_notebook()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## III - Merging flags and stellarity\n", "\n", "This all happens at the end now after the catalogue has been cut into strips." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## IV - Adding E(B-V) column" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": true }, "outputs": [], "source": [ "master_catalogue.add_column(\n", " ebv(master_catalogue['ra'], master_catalogue['dec'])\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## V - Adding HELP unique identifiers and field columns" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": true }, "outputs": [], "source": [ "master_catalogue.add_column(Column(gen_help_id(master_catalogue['ra'], master_catalogue['dec']),\n", " name=\"help_id\"))\n", "master_catalogue.add_column(Column(np.full(len(master_catalogue), \"XMM-LSS\", dtype='" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nb_merge_dist_plot(\n", " SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n", " SkyCoord(specz['ra'] * u.deg, specz['dec'] * u.deg)\n", ")" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": true }, "outputs": [], "source": [ "master_catalogue = specz_merge(master_catalogue, specz, radius=1. * u.arcsec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## VIII.a Wavelength domain coverage\n", "\n", "We add a binary `flag_optnir_obs` indicating that a source was observed in a given wavelength domain:\n", "\n", "- 1 for observation in optical;\n", "- 2 for observation in near-infrared;\n", "- 4 for observation in mid-infrared (IRAC).\n", "\n", "It's an integer binary flag, so a source observed both in optical and near-infrared by not in mid-infrared would have this flag at 1 + 2 = 3.\n", "\n", "*Note 1: The observation flag is based on the creation of multi-order coverage maps from the catalogues, this may not be accurate, especially on the edges of the coverage.*\n", "\n", "*Note 2: Being on the observation coverage does not mean having fluxes in that wavelength domain. For sources observed in one domain but having no flux in it, one must take into consideration the different depths in the catalogue we are using.*" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#candels_moc = MOC(filename=\"../../dmu0/dmu0_CANDELS-3D-HST/data/CANDELS-3D-HST_XMM-LSS_MOC.fits\") # 1.1\n", "candels_moc = MOC(filename=\"../../dmu0/dmu0_CANDELS-UDS/data/hlsp_candels_hst_wfc3_uds-tot-multiband_f160w_v1_MOC.fits\") # 1.2\n", "cfht_wirds_moc = MOC(filename=\"../../dmu0/dmu0_CFHT-WIRDS/data/XMM-LSS_Ks-priors_MOC.fits\") # 1.3\n", "cfhtls_wide_moc = MOC(filename=\"../../dmu0/dmu0_CFHTLS/data/CFHTLS-WIDE_XMM-LSS_MOC.fits\") # 1.4a\n", "cfhtls_deep_moc = MOC(filename=\"../../dmu0/dmu0_CFHTLS/data/CFHTLS-DEEP_XMM-LSS_MOC.fits\") # 1.4b\n", "cfhtls_moc = cfhtls_wide_moc + cfhtls_deep_moc\n", "#cfhtlens_moc = MOC(filename=\"../../dmu0/dmu0_CFHTLenS/data/CFHTLenS_XMM-LSS_MOC.fits\") # 1.5\n", "decals_moc = MOC(filename=\"../../dmu0/dmu0_DECaLS/data/DECaLS_XMM-LSS_MOC.fits\") # 1.6\n", "des_moc = MOC(filename=\"../../dmu0/dmu0_DES/data/DES-DR1_XMM-LSS_MOC.fits\") # 1.6\n", "decam_moc = decals_moc + des_moc\n", "servs_moc = MOC(filename=\"../../dmu0/dmu0_DataFusion-Spitzer/data/DF-SERVS_XMM-LSS_MOC.fits\") # 1.8\n", "swire_moc = MOC(filename=\"../../dmu0/dmu0_DataFusion-Spitzer/data/DF-SWIRE_XMM-LSS_MOC.fits\") # 1.7\n", "hsc_wide_moc = MOC(filename=\"../../dmu0/dmu0_HSC/data/HSC-PDR1_wide_XMM-LSS_MOC.fits\") # 1.9a\n", "hsc_deep_moc = MOC(filename=\"../../dmu0/dmu0_HSC/data/HSC-PDR1_deep_XMM-LSS_MOC.fits\") # 1.9b\n", "hsc_udeep_moc = MOC(filename=\"../../dmu0/dmu0_HSC/data/HSC-PDR1_uDeep_XMM-LSS_MOC.fits\") # 1.9c\n", "hsc_moc = hsc_wide_moc + hsc_deep_moc + hsc_udeep_moc\n", "ps1_moc = MOC(filename=\"../../dmu0/dmu0_PanSTARRS1-3SS/data/PanSTARRS1-3SS_XMM-LSS_MOC.fits\") # 1.10\n", "sxds_moc = MOC(filename=\"../../dmu0/dmu0_SXDS/data/dmu0_SXDS_MOC.fits\") # 1.11\n", "sparcs_moc = MOC(filename=\"../../dmu0/dmu0_SpARCS/data/SpARCS_HELP_XMM-LSS_MOC.fits\") # 1.12\n", "dxs_moc = MOC(filename=\"../../dmu0/dmu0_UKIDSS-DXS_DR10plus/data/UKIDSS-DR10plus_XMM-LSS_MOC.fits\") # 1.13\n", "uds_moc = MOC(filename=\"../../dmu0/dmu0_UKIDSS-UDS/data/UKIDSS-UDS_XMM-LSS_MOC.fits\") # 1.14\n", "vipers_moc = MOC(filename=\"../../dmu0/dmu0_VIPERS-MLS/data/VIPERS-MLS_20160502_MOC.fits\") # 1.15\n", "vhs_moc = MOC(filename=\"../../dmu0/dmu0_VISTA-VHS/data/VHS_XMM-LSS_MOC.fits\") # 1.16\n", "video_moc = MOC(filename=\"../../dmu0/dmu0_VISTA-VIDEO-private/data/VIDEO-all_2016-04-14_fullcat_errfix_XMM-LSS_MOC.fits\") # 1.17\n", "viking_moc = MOC(filename=\"../../dmu0/dmu0_VISTA-VIKING/data/VIKING_XMM-LSS_MOC.fits\") # 1.18" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": true }, "outputs": [], "source": [ "was_observed_optical = inMoc(\n", " master_catalogue['ra'], master_catalogue['dec'],\n", " (candels_moc + \n", " cfht_wirds_moc + \n", " cfhtls_moc + \n", " #cfhtlens_moc + \n", " sparcs_moc + \n", " decam_moc + \n", " hsc_moc + \n", " ps1_moc) )\n", "\n", "was_observed_nir = inMoc(\n", " master_catalogue['ra'], master_catalogue['dec'],\n", " dxs_moc + uds_moc + vhs_moc + video_moc + viking_moc\n", ")\n", "\n", "was_observed_mir = inMoc(\n", " master_catalogue['ra'], master_catalogue['dec'],\n", " servs_moc + swire_moc\n", ")" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": true }, "outputs": [], "source": [ "master_catalogue.add_column(\n", " Column(\n", " 1 * was_observed_optical + 2 * was_observed_nir + 4 * was_observed_mir,\n", " name=\"flag_optnir_obs\")\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## VIII.b Wavelength domain detection\n", "\n", "We add a binary `flag_optnir_det` indicating that a source was detected in a given wavelength domain:\n", "\n", "- 1 for detection in optical;\n", "- 2 for detection in near-infrared;\n", "- 4 for detection in mid-infrared (IRAC).\n", "\n", "It's an integer binary flag, so a source detected both in optical and near-infrared by not in mid-infrared would have this flag at 1 + 2 = 3.\n", "\n", "*Note 1: We use the total flux columns to know if the source has flux, in some catalogues, we may have aperture flux and no total flux.*\n", "\n", "To get rid of artefacts (chip edges, star flares, etc.) we consider that a source is detected in one wavelength domain when it has a flux value in **at least two bands**. That means that good sources will be excluded from this flag when they are on the coverage of only one band." ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# SpARCS is a catalogue of sources detected in r (with fluxes measured at \n", "# this prior position in the other bands). Thus, we are only using the r\n", "# CFHT band.\n", "# Check to use catalogue flags from HSC and PanSTARRS.\n", "#nb_optical_flux = np.zeros(len(master_catalogue), dtype=float)\n", "\n", "#nb_nir_flux = np.zeros(len(master_catalogue), dtype=float)\n", "\n", "#nb_mir_flux = np.zeros(len(master_catalogue), dtype=float)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#has_optical_flux = nb_optical_flux >= 2\n", "#has_nir_flux = nb_nir_flux >= 2\n", "#has_mir_flux = nb_mir_flux >= 2\n", "\n", "#master_catalogue.add_column(\n", "# Column(\n", "# 1 * has_optical_flux + 2 * has_nir_flux + 4 * has_mir_flux,\n", "# name=\"flag_optnir_det\")\n", "#)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## IX - Cross-identification table\n", "\n", "We are producing a table associating to each HELP identifier, the identifiers of the sources in the pristine catalogues. This can be used to easily get additional information from them.\n", "\n", "For convenience, we also cross-match the master list with the SDSS catalogue and add the objID associated with each source, if any. **TODO: should we correct the astrometry with respect to Gaia positions?**" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "315 master list rows had multiple associations.\n" ] } ], "source": [ "#\n", "# Addind SDSS ids\n", "#\n", "sdss = Table.read(\"../../dmu0/dmu0_SDSS-DR13/data/SDSS-DR13_XMM-LSS.fits\")['objID', 'ra', 'dec']\n", "sdss_coords = SkyCoord(sdss['ra'] * u.deg, sdss['dec'] * u.deg)\n", "idx_ml, d2d, _ = sdss_coords.match_to_catalog_sky(SkyCoord(master_catalogue['ra'], master_catalogue['dec']))\n", "idx_sdss = np.arange(len(sdss))\n", "\n", "# Limit the cross-match to 1 arcsec\n", "mask = d2d <= 1. * u.arcsec\n", "idx_ml = idx_ml[mask]\n", "idx_sdss = idx_sdss[mask]\n", "d2d = d2d[mask]\n", "nb_orig_matches = len(idx_ml)\n", "\n", "# In case of multiple associations of one master list object to an SDSS object, we keep only the\n", "# association to the nearest one.\n", "sort_idx = np.argsort(d2d)\n", "idx_ml = idx_ml[sort_idx]\n", "idx_sdss = idx_sdss[sort_idx]\n", "_, unique_idx = np.unique(idx_ml, return_index=True)\n", "idx_ml = idx_ml[unique_idx]\n", "idx_sdss = idx_sdss[unique_idx]\n", "print(\"{} master list rows had multiple associations.\".format(nb_orig_matches - len(idx_ml)))\n", "\n", "# Adding the ObjID to the master list\n", "master_catalogue.add_column(Column(data=np.full(len(master_catalogue), -1, dtype='>i8'), name=\"sdss_id\"))\n", "master_catalogue['sdss_id'][idx_ml] = sdss['objID'][idx_sdss]" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['ps1_id', 'cfht_intid', 'candels_id', 'cfhtls-wide_id', 'cfhtls-deep_id', 'sparcs_intid', 'wirds_id', 'vipers_id', 'sxds_intid', 'sxds_b_id', 'sxds_v_id', 'sxds_r_id', 'sxds_i_id', 'sxds_z_id', 'hsc_intid', 'hsc-wide_id', 'hsc-deep_id', 'hsc-udeep_id', 'decam_intid', 'decals_id', 'des_id', 'ukidss_intid', 'dxs_id', 'uds_id', 'vircam_intid', 'vhs_id', 'video_id', 'viking_id', 'irac_intid', 'servs_intid', 'swire_intid', 'help_id', 'specz_id', 'sdss_id']\n" ] } ], "source": [ "id_names = []\n", "for col in master_catalogue.colnames:\n", " if '_id' in col:\n", " id_names += [col]\n", " if '_intid' in col:\n", " id_names += [col]\n", " \n", "print(id_names)\n" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#master_catalogue[id_names].write(\n", "# \"{}/master_list_cross_ident_xmm-lss{}.fits\".format(OUT_DIR, SUFFIX), overwrite=True)\n", "id_names.remove('help_id')\n", "#master_catalogue.remove_columns(id_names)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## X - Adding HEALPix index\n", "\n", "We are adding a column with a HEALPix index at order 13 associated with each source." ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": true }, "outputs": [], "source": [ "master_catalogue.add_column(Column(\n", " data=coords_to_hpidx(master_catalogue['ra'], master_catalogue['dec'], order=13),\n", " name=\"hp_idx\"\n", "))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## XI - Saving the catalogue" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": true }, "outputs": [], "source": [ "columns = [\"help_id\", \"field\", \"ra\", \"dec\", \"hp_idx\"]\n", "\n", "bands = [column[5:] for column in master_catalogue.colnames if 'f_ap' in column]\n", "for band in bands:\n", " columns += [\"f_ap_{}\".format(band), \"ferr_ap_{}\".format(band),\n", " \"m_ap_{}\".format(band), \"merr_ap_{}\".format(band),\n", " \"f_{}\".format(band), \"ferr_{}\".format(band),\n", " \"m_{}\".format(band), \"merr_{}\".format(band),\n", " \"flag_{}\".format(band)] \n", " \n", "columns += [\"stellarity\", \"stellarity_origin\", \"flag_cleaned\", \"flag_merged\", \"flag_gaia\", \"flag_optnir_obs\", \"flag_optnir_det\", \n", " \"zspec\", \"zspec_qual\", \"zspec_association_flag\", \"ebv\"]" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Missing columns: {'vhs_id', 'irac_intid', 'candels_id', 'decals_id', 'sxds_intid', 'vipers_id', 'des_id', 'swire_intid', 'dxs_id', 'cfhtls-deep_id', 'sxds_z_id', 'vircam_intid', 'sxds_r_id', 'wirds_id', 'sxds_v_id', 'cfhtls-wide_id', 'hsc-wide_id', 'ukidss_intid', 'sxds_i_id', 'video_id', 'ps1_id', 'viking_id', 'servs_intid', 'hsc-deep_id', 'hsc_intid', 'sdss_id', 'sxds_b_id', 'hsc-udeep_id', 'decam_intid', 'cfht_intid', 'specz_id', 'uds_id', 'sparcs_intid'}\n" ] } ], "source": [ "# We check for columns in the master catalogue that we will not save to disk.\n", "print(\"Missing columns: {}\".format(set(master_catalogue.colnames) - set(columns)))" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#master_catalogue[columns].write(\"{}/master_catalogue_xmm-lss_low-memory{}.fits\".format(OUT_DIR, SUFFIX), overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## XII - folding in the photometry\n", "On XMM-LSS there is too much data to load all in to memory at once so we perform the cross matching without photometry columns. Only now do we fold in the photometry data by first cutting the catalogue up in to manageable sizes." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "88\n" ] } ], "source": [ "split_length = 100000 #number of sources to include in every sub catalogue\n", "num_files = int(np.ceil(len(master_catalogue)/split_length))\n", "print(num_files)\n" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#Sort catalogue by HELP id so that it is split up in RA strips\n", "master_catalogue.sort('help_id')\n", "#Remove all the old ids as they interfere with join\n", "old_id_names = [\n", " 'candels_id',\n", " 'cfhtls-wide_id',\n", " 'cfhtls-deep_id',\n", " 'sparcs_intid',\n", " 'wirds_id',\n", " 'vipers_id',\n", " 'sxds_intid',\n", " 'sxds_b_id',\n", " 'sxds_v_id',\n", " 'sxds_r_id',\n", " 'sxds_i_id',\n", " 'sxds_z_id',\n", " \n", " 'hsc-wide_id',\n", " 'hsc-deep_id',\n", " 'hsc-udeep_id',\n", " \n", " 'des_id',\n", " 'decals_id',\n", "\n", " 'dxs_id',\n", " 'uds_id',\n", " \n", " 'vhs_id',\n", " 'video_id',\n", " 'viking_id',\n", "\n", " 'servs_intid',\n", " 'swire_intid']\n", "master_catalogue.remove_columns(old_id_names)\n", "# The old id names will be added back in by \n", "# the join so I leave them in the list here to be \n", "# removed before the save of the sub catalogue\n", "# id_names.remove(old_id_names)" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#There is candels data on three seperate merged catalogues so duplicates must be removed from two\n", "irac_merged = Table.read(\"{}/irac_merged_catalogue_xmm-lss.fits\".format(TMP_DIR))\n", "irac_merged.remove_columns(['candels_id',\n", " 'candels_stellarity',\n", " 'candels_flag_cleaned',\n", " 'candels_flag_gaia'])\n", "\n", "ukidss_merged = Table.read(\"{}/ukidss_merged_catalogue_xmm-lss.fits\".format(TMP_DIR))\n", "ukidss_merged.remove_columns(['candels_id',\n", " 'candels_stellarity',\n", " 'candels_flag_cleaned',\n", " 'candels_flag_gaia'])" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "collapsed": true }, "outputs": [], "source": [ "surveys = [\n", " \n", "\n", " #['sxds', \"SXDS.fits\", Table.read(\"{}/SXDS.fits\".format(TMP_DIR)), 'sxds_intid'], \n", " #['candels', \"CANDELS.fits\", Table.read(\"{}/CANDELS.fits\".format(TMP_DIR)), 'candels_id'],\n", " #['cfhtls_wide', \"CFHTLS-WIDE.fits\", Table.read(\"{}/CFHTLS-WIDE.fits\".format(TMP_DIR)), 'cfhtls-wide_id'],\n", " #['cfhtls_deep', \"CFHTLS-DEEP.fits\", Table.read(\"{}/CFHTLS-DEEP.fits\".format(TMP_DIR)), 'cfhtls-deep_id'],\n", " #['sparcs', \"SpARCS.fits\", Table.read(\"{}/SpARCS.fits\".format(TMP_DIR)), 'sparcs_intid'], \n", " #['cfht_wirds', \"CFHT-WIRDS.fits\", Table.read(\"{}/CFHT-WIRDS.fits\".format(TMP_DIR)), 'wirds_intid' ],\n", " #['vipers', \"VIPERS.fits\", Table.read(\"{}/VIPERS.fits\".format(TMP_DIR)), 'vipers_id'],\n", " ['cfht', \n", " \"cfht_merged_catalogue_xmm-lss.fits\", \n", " Table.read(\"{}/cfht_merged_catalogue_xmm-lss.fits\".format(TMP_DIR)), \n", " 'cfht_intid'], \n", " \n", " #['decals', \"DECaLS.fits\", Table.read(\"{}/DECaLS.fits\".format(TMP_DIR)), 'decals_id'],\n", " ['decam', \n", " \"decam_merged_catalogue_xmm-lss.fits\", \n", " Table.read(\"{}/decam_merged_catalogue_xmm-lss.fits\".format(TMP_DIR)), \n", " 'decam_intid'],\n", " \n", " #['servs', \"SERVS.fits\", Table.read(\"{}/SERVS.fits\".format(TMP_DIR)), 'servs_intid'],\n", " #['swire', \"SWIRE.fits\", Table.read(\"{}/SWIRE.fits\".format(TMP_DIR)), 'swire_intid'], \n", " ['irac', \n", " \"irac_merged_catalogue_xmm-lss.fits\",\n", " irac_merged, \n", " 'irac_intid'], \n", " \n", " #['hsc_wide ', \"HSC-WIDE.fits\", Table.read(\"{}/HSC-WIDE.fits\".format(TMP_DIR)), 'hsc-wide_id' ],\n", " #['hsc_deep', \"HSC-DEEP.fits\", Table.read(\"{}/HSC-DEEP.fits\".format(TMP_DIR)), 'hsc-deep_id'],\n", " #['hsc_udeep', \"HSC-UDEEP.fits\", Table.read(\"{}/HSC-UDEEP.fits\".format(TMP_DIR)), 'hsc-udeep_id'], \n", " ['hsc', \n", " \"hsc_merged_catalogue_xmm-lss.fits\", \n", " Table.read(\"{}/hsc_merged_catalogue_xmm-lss.fits\".format(TMP_DIR)), \n", " 'hsc_intid'], \n", " \n", " ['ps1', \"PS1.fits\", Table.read(\"{}/PS1.fits\".format(TMP_DIR)), 'ps1_id'], \n", " \n", " \n", " \n", " #['dxs', \"UKIDSS-DXS.fits\", Table.read(\"{}/UKIDSS-DXS.fits\".format(TMP_DIR)), 'dxs_id'], \n", " #['uds', \"UKIDSS-UDS.fits\", Table.read(\"{}/UKIDSS-UDS.fits\".format(TMP_DIR)), 'uds_id'], \n", " ['ukidss', \n", " \"ukidss_merged_catalogue_xmm-lss.fits\", \n", " ukidss_merged, \n", " 'ukidss_intid'], \n", " \n", " \n", " \n", " #['vhs', \"VISTA-VHS.fits\", Table.read(\"{}/VISTA-VHS.fits\".format(TMP_DIR)), 'vhs_id'], \n", " #['video', \"VISTA-VIDEO.fits\", Table.read(\"{}/VISTA-VIDEO.fits\".format(TMP_DIR)), 'video_id'], \n", " #['viking', \"VISTA-VIKING.fits\", Table.read(\"{}/VISTA-VIKING.fits\".format(TMP_DIR)), 'viking_id'] \n", " ['vircam', \n", " \"vista_merged_catalogue_xmm-lss.fits\", \n", " Table.read(\"{}/vista_merged_catalogue_xmm-lss.fits\".format(TMP_DIR)), \n", " 'vircam_intid']\n", "]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(np.sum(master_catalogue['flag_merged'])/len(master_catalogue)) " ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "collapsed": true }, "outputs": [], "source": [ "n=0\n", "for sub_file in range(num_files):\n", " print(\"Run: {}\".format(sub_file))\n", " sub_catalogue = master_catalogue[n*split_length:(n+1)*split_length]\n", " #print(n)\n", " print(np.sum(sub_catalogue['flag_merged'])/len(sub_catalogue)) \n", " for survey in surveys:\n", " #print(survey[0])\n", " sub_catalogue = join(sub_catalogue, \n", " survey[2], #Table.read(\"{}/{}\".format(TMP_DIR, survey[1])),\n", " join_type='left',\n", " metadata_conflicts='silent',\n", " keys=survey[3]\n", " )\n", " \n", " \n", " \n", " #Each pristine catalogue contains a flag indicating if the source was associated to a another nearby source that was removed during the cleaning process. We merge these flags in a single one.\n", "\n", "\n", " flag_cleaned_columns = [column for column in sub_catalogue.colnames\n", " if 'flag_cleaned' in column]\n", "\n", " flag_column = np.zeros(len(sub_catalogue), dtype=bool)\n", " for column in flag_cleaned_columns:\n", " flag_column |= sub_catalogue[column]\n", " \n", " sub_catalogue.add_column(Column(data=flag_column, name=\"flag_cleaned\"))\n", " sub_catalogue.remove_columns(flag_cleaned_columns)\n", " #combining the flag_merged column which contains information regarding multiple associations\n", " \n", " print(np.sum(sub_catalogue['flag_merged'])/len(sub_catalogue)) \n", " sub_catalogue['flag_merged'].name = 'flag_merged_tmp'\n", " flag_merged_columns = [column for column in sub_catalogue.colnames\n", " if 'flag_merged' in column]\n", " \n", " flag_merged_column = np.zeros(len(sub_catalogue), dtype=bool)\n", " for column in flag_merged_columns:\n", " flag_merged_column |= sub_catalogue[column]\n", " print(\"flag_merged_count: {}\".format(np.sum(flag_merged_column)))\n", " print(np.sum(sub_catalogue['flag_merged'])/len(sub_catalogue)) \n", " \n", " sub_catalogue.add_column(Column(data=flag_merged_column, name=\"flag_merged\"))\n", " sub_catalogue.remove_columns(flag_merged_columns)\n", " print(np.sum(sub_catalogue['flag_merged'])/len(sub_catalogue)) \n", "\n", " #Each pristine catalogue contains a flag indicating the probability of a source being a Gaia object (0: not a Gaia object, 1: possibly, 2: probably, 3: definitely). We merge these flags taking the highest value.\n", "\n", "\n", " flag_gaia_columns = [column for column in sub_catalogue.colnames\n", " if 'flag_gaia' in column]\n", "\n", " sub_catalogue.add_column(Column(\n", " data=np.max([sub_catalogue[column] for column in flag_gaia_columns], axis=0),\n", " name=\"flag_gaia\"\n", " ))\n", " sub_catalogue.remove_columns(flag_gaia_columns)\n", " #Each prisitine catalogue may contain one or several stellarity columns indicating the #probability (0 to 1) of each source being a star. We merge these columns taking the highest value. We keep trace of the origin of the stellarity.\n", "\n", "\n", "\n", " stellarity_columns = [column for column in sub_catalogue.colnames\n", " if 'stellarity' in column]\n", "\n", " #print(\", \".join(stellarity_columns))\n", "\n", "\n", "\n", " # We create an masked array with all the stellarities and get the maximum value, as well as its\n", " # origin. Some sources may not have an associated stellarity.\n", " stellarity_array = np.array([sub_catalogue[column] for column in stellarity_columns])\n", " stellarity_array = np.ma.masked_array(stellarity_array, np.isnan(stellarity_array))\n", "\n", " max_stellarity = np.max(stellarity_array, axis=0)\n", " max_stellarity.fill_value = np.nan\n", "\n", " no_stellarity_mask = max_stellarity.mask\n", "\n", " sub_catalogue.add_column(Column(data=max_stellarity.filled(), name=\"stellarity\"))\n", "\n", " stellarity_origin = np.full(len(sub_catalogue), \"NO_INFORMATION\", dtype=\"S20\")\n", " stellarity_origin[~no_stellarity_mask] = np.array(stellarity_columns)[np.argmax(stellarity_array, axis=0)[~no_stellarity_mask]]\n", "\n", " sub_catalogue.add_column(Column(data=stellarity_origin, name=\"stellarity_origin\"))\n", "\n", " sub_catalogue.remove_columns(stellarity_columns)\n", "\n", " for col in sub_catalogue.colnames:\n", " #print(col)\n", " if (col.startswith(\"m_\") or col.startswith(\"merr_\") or col.startswith(\"f_\") or col.startswith(\"ferr_\")):\n", " sub_catalogue[col] = sub_catalogue[col].astype(float)\n", " sub_catalogue[col].fill_value = np.nan\n", " elif \"flag\" in col:\n", " sub_catalogue[col].fill_value = 0\n", " elif \"id\" in col:\n", " sub_catalogue[col].fill_value = -1\n", " \n", " sub_catalogue = sub_catalogue.filled()\n", " \n", " \n", " nb_optical_flux = (\n", " # CANDELS\n", " 1 * ~np.isnan(sub_catalogue['f_acs_f606w']) + \n", " 1 * ~np.isnan(sub_catalogue['f_acs_f814w']) +\n", " 1 * ~np.isnan(sub_catalogue['f_wfc3_f125w']) +\n", " 1 * ~np.isnan(sub_catalogue['f_wfc3_f160w']) +\n", "\n", " # DES and DECaLS\n", " 1 * ~np.isnan(sub_catalogue['f_decam_g']) + \n", " 1 * ~np.isnan(sub_catalogue['f_decam_r']) +\n", " 1 * ~np.isnan(sub_catalogue['f_decam_i']) +\n", " 1 * ~np.isnan(sub_catalogue['f_decam_z']) +\n", " 1 * ~np.isnan(sub_catalogue['f_decam_y']) +\n", " # HSC\n", " 1 * ~np.isnan(sub_catalogue['f_suprime_g']) + \n", " 1 * ~np.isnan(sub_catalogue['f_suprime_r']) +\n", " 1 * ~np.isnan(sub_catalogue['f_suprime_i']) +\n", " 1 * ~np.isnan(sub_catalogue['f_suprime_z']) +\n", " 1 * ~np.isnan(sub_catalogue['f_suprime_y']) +\n", " 1 * ~np.isnan(sub_catalogue['f_suprime_n921']) +\n", " 1 * ~np.isnan(sub_catalogue['f_suprime_n816']) +\n", " # PanSTARRS\n", " 1 * ~np.isnan(sub_catalogue['f_gpc1_g']) +\n", " 1 * ~np.isnan(sub_catalogue['f_gpc1_r']) +\n", " 1 * ~np.isnan(sub_catalogue['f_gpc1_i']) +\n", " 1 * ~np.isnan(sub_catalogue['f_gpc1_z']) +\n", " 1 * ~np.isnan(sub_catalogue['f_gpc1_y']) +\n", " # CFHT\n", " 1 * ~np.isnan(sub_catalogue['f_megacam_u']) +\n", " 1 * ~np.isnan(sub_catalogue['f_megacam_g']) +\n", " 1 * ~np.isnan(sub_catalogue['f_megacam_r']) +\n", " 1 * ~np.isnan(sub_catalogue['f_megacam_z']) +\n", " 1 * ~np.isnan(sub_catalogue['f_megacam_y']) \n", " )\n", "\n", " nb_nir_flux = (\n", " # CFHT-WIRCAM\n", " 1 * ~np.isnan(sub_catalogue['f_wircam_j']) +\n", " 1 * ~np.isnan(sub_catalogue['f_wircam_h']) +\n", " 1 * ~np.isnan(sub_catalogue['f_wircam_ks']) +\n", " # UKIDSS\n", " 1 * ~np.isnan(sub_catalogue['f_ukidss_j']) +\n", " 1 * ~np.isnan(sub_catalogue['f_ukidss_h']) +\n", " 1 * ~np.isnan(sub_catalogue['f_ukidss_k']) +\n", " # VISTA\n", " 1 * ~np.isnan(sub_catalogue['f_vista_z']) +\n", " 1 * ~np.isnan(sub_catalogue['f_vista_y']) +\n", " 1 * ~np.isnan(sub_catalogue['f_vista_j']) +\n", " 1 * ~np.isnan(sub_catalogue['f_vista_h']) +\n", " 1 * ~np.isnan(sub_catalogue['f_vista_ks'])\n", " )\n", "\n", " nb_mir_flux = (\n", " 1 * ~np.isnan(sub_catalogue['f_irac_i1']) +\n", " 1 * ~np.isnan(sub_catalogue['f_irac_i2']) +\n", " 1 * ~np.isnan(sub_catalogue['f_irac_i3']) +\n", " 1 * ~np.isnan(sub_catalogue['f_irac_i4']) \n", " )\n", " \n", " has_optical_flux = nb_optical_flux >= 2\n", " has_nir_flux = nb_nir_flux >= 2\n", " has_mir_flux = nb_mir_flux >= 2\n", "\n", " sub_catalogue.add_column(\n", " Column(\n", " 1 * has_optical_flux + 2 * has_nir_flux + 4 * has_mir_flux,\n", " name=\"flag_optnir_det\")\n", " )\n", " \n", " #print('Finished join')\n", " sub_catalogue.remove_columns(id_names)\n", " #sub_catalogue.remove_columns(['sxds_b_id', 'sxds_v_id', 'sxds_r_id', 'sxds_i_id', 'sxds_z_id'])\n", " \n", " for col in sub_catalogue.colnames:\n", " if col.endswith('_ra'):\n", " sub_catalogue.remove_column(col)\n", " if col.endswith('_dec'):\n", " sub_catalogue.remove_column(col)\n", " if '_vircam_' in col:\n", " sub_catalogue.rename_column(col, col.replace('_vircam_', '_vista_'))\n", " \n", "\n", " \n", " \n", " \n", " #sub_catalogue.write(\"{}/tiles/sub_catalogue_xmm-lss{}_{}.fits\".format(OUT_DIR, SUFFIX, n), overwrite=True)\n", " \n", " \n", " del sub_catalogue\n", " del flag_column\n", " del flag_merged_column\n", " gc.collect()\n", " time.sleep(10)\n", " \n", " n += 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How to generate final catalogue\n", "After this notebook has been run there will be a set of sub catalogues in data/tiles/\n", "\n", "These need to be stacked using stilts:\n", "\n", "
\n",
    "suffix=20180111\n",
    "\n",
    "ls ./data/tiles/sub_catalogue_xmm-lss_$suffix_*.fits > ./data/tiles/fits_list_$suffix.lis\n",
    "\n",
    "stilts tcat in=@./data/tiles/fits_list_$suffix.lis out=./data/master_catalogue_xmm-lss_$suffix.fits\n",
    "
\n", "\n", "For many purposes this file may be too large. In order to run checks and diagnostics we typically take a subset using something like:\n", "\n", "
\n",
    "stilts tpipe cmd='every 10' ./data/master_catalogue_xmm-lss_$suffix.fits omode=out out=./data/master_catalogue_xmm-lss_RANDOM10PCSAMPLE_$suffix.fits\n",
    "
" ] } ], "metadata": { "kernelspec": { "display_name": "Python (herschelhelp_internal)", "language": "python", "name": "helpint" }, "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.8" } }, "nbformat": 4, "nbformat_minor": 2 }