PanSTARRS_WISE_ML.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  PanSTARRS - WISE cross-match: Maximum Likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy.table import Table, join\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord, search_around_sky\n",
    "import pickle\n",
    "from tqdm import tnrange, tqdm_notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import clear_output\n",
    "from mltier1 import parallel_process, SingleMLEstimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%autoreload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data\n",
    "\n",
    "### Load the catalogues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "panstarrs = Table.read(\"panstarrs_u2.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wise = Table.read(\"wise_u2.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21650894"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(panstarrs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13219700"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(wise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coordinates\n",
    "\n",
    "As we will use the coordinates to make a cross-match we to load them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "coords_panstarrs = SkyCoord(panstarrs['raMean'], panstarrs['decMean'], unit=(u.deg, u.deg), frame='icrs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "coords_wise = SkyCoord(wise['raWise'], wise['decWise'], unit=(u.deg, u.deg), frame='icrs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import the ML parameeters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bin_list, center, q0, n_m, q_m = pickle.load(open(\"pw_params.pckl\", \"rb\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Define the main functions\n",
    "\n",
    "We define the main functions used in the ML computing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "likelihood_ratio = SingleMLEstimator(q0, n_m, q_m, center)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function could be used to get the sigma using the errors in the two optical catalogues but it may underestimate the error (not used)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_sigma(ra1, dec1, ra1err, dec1err, ra2, dec2, ra2err, dec2err):\n",
    "    \"\"\"Input positions in degrees.\n",
    "    Errors in arcsecs\n",
    "    Output in arcsecs\n",
    "    \"\"\"\n",
    "    cosadj = np.cos(np.deg2rad(0.5*(dec1 + dec2)))\n",
    "    phi = np.arctan2((dec2 - dec1), ((ra2 - ra1)*cosadj))\n",
    "    sigma = np.pi - phi\n",
    "    err1squared  = (ra1err * np.cos(sigma))**2 + (dec1err * np.sin(sigma))**2\n",
    "    err2squared  = (ra2err * np.cos(phi))**2 + (dec2err * np.sin(phi))**2\n",
    "    return np.sqrt(err1squared + err2squared)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the ML ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "radius = 15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "idx_wise, idx_panstarrs, d2d, d3d = search_around_sky(\n",
    "    coords_wise, coords_panstarrs, radius*u.arcsec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "idx_wise_unique = np.unique(idx_wise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wise[\"lr\"] = np.nan                   # Likelihood ratio\n",
    "wise[\"lr_dist\"] = np.nan              # Distance to the selected source\n",
    "wise[\"lr_panstarrs_index\"] = np.nan   # Index of the PanSTARRS source"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "total_sources = len(idx_wise_unique)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11539194"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_sources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "panstarrs_aux_index = np.arange(len(panstarrs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ml(i):\n",
    "    idx_0 = idx_panstarrs[idx_wise == i]\n",
    "    d2d_0 = d2d[idx_wise == i]\n",
    "    i_mag = panstarrs[\"i\"][idx_0]\n",
    "#     sigma = get_sigma(wise[\"ra\"][i], \n",
    "#                       wise[\"dec\"][i], \n",
    "#                       wise[\"sigra\"][i], \n",
    "#                       wise[\"sigdec\"][i],\n",
    "#                       panstarrs[idx_0][\"raMean\"][idx_0],\n",
    "#                       panstarrs[idx_0][\"decMean\"][idx_0],\n",
    "#                       panstarrs[idx_0][\"raMeanErr\"][idx_0],\n",
    "#                       panstarrs[idx_0][\"decMeanErr\"][idx_0],\n",
    "#              )\n",
    "    sigma = 1.\n",
    "    lr_0 = likelihood_ratio(i_mag, d2d_0.arcsec, sigma)\n",
    "    chosen_index = np.argmax(lr_0)\n",
    "    \n",
    "    result = [panstarrs_aux_index[idx_0[chosen_index]], # Index\n",
    "              (d2d_0.arcsec)[chosen_index],             # distance\n",
    "              lr_0[chosen_index]]                       # LR\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "step_size = 100000\n",
    "nsteps = total_sources//step_size + 1\n",
    "res = []\n",
    "for k in tnrange(nsteps, desc=\"Blocks\"):\n",
    "    low_limit = k*step_size\n",
    "    high_limit = (k+1)*step_size\n",
    "    res += parallel_process(idx_wise_unique[low_limit:high_limit], \n",
    "                            ml, \n",
    "                            n_jobs=8,\n",
    "                            notebook=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11539194"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "(wise[\"lr_panstarrs_index\"][idx_wise_unique], \n",
    " wise[\"lr_dist\"][idx_wise_unique], \n",
    " wise[\"lr\"][idx_wise_unique]) = list(map(list, zip(*res)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the ML results\n",
    "\n",
    "The auxiliary column \"lr_pc\" will have a 0 value for all non-matches and will be used to get the threshold in the ML parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wise[\"lr_pc\"] = wise[\"lr\"]\n",
    "wise[\"lr_pc\"][np.isnan(wise[\"lr_pc\"])] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "threshold = np.percentile(wise[\"lr_pc\"], 100*(1 - q0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.4306303253655415"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAFpCAYAAADEC+ZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+wbeV5H/bvY66FkW0hEIQQfvSQiFEHkXgi3cHEzqS0\nZAQxjlFnJM3VxBZJqJhUiuukmboXdyZ07GEG2owVK63UYQwBKbIQJXbF+JrIFKLxNAlIV7YcBBhz\nK5C4NyAQYIibCuvSp3/sdZTN4f4495x9zt77rM9nZs9Z+91rrfPudddd7/qu913rVHcHAACA8fm+\neVcAAACA+RAIAQAARkogBAAAGCmBEAAAYKQEQgAAgJESCAEAAEZKIAQAABgpgRAAAGCkBEIAAICR\nEggBAABGate8KzBrZ5xxRq+srGx6PQ8fejl//pxTN1+hweOPP54kecc73jGzdQKM3Ve+8pVvd/eZ\n867HsphFG/nwoZeTZKZtJOPl/Ih5Wj2eLZpZHF9PpH3ccYFwZWUl+/fv3/x69u7L/puumkGNJi67\n7LIkyRe/+MWZrRNg7KrqG/OuwzKZRRu5sndfksy0jWS8nB8xT6vHs0Uzi+PribSPhowCAACMlEAI\nAAAwUjtuyCgAAMDRLOpQ0XnRQwgAADBSAiEAAMBICYQAAAAjJRACAACMlEAIAAAwUgIhAADASAmE\nAAAAIyUQAgAAjJRACAAAMFICIQAAwEgJhAAAACMlEAIAAIzUcQNhVd1WVc9V1deO8Nk/qKquqjOm\nyq6vqgNV9XhVXTFV/u6qenj47ONVVUP5yVX1uaH8oapamVrmmqp6Ynhds9kvCwCzpI0EWC4re/fN\nuwoLZz09hLcnuXJtYVWdl+Q9Sb45VXZRkj1J3jks84mqOmn4+JNJPpzkwuG1us5rk7zU3W9P8rEk\nNw/rOj3JDUl+NMklSW6oqtNO7OsBwJa6PdpIAJbYcQNhd/9OkheP8NHHkvx8kp4quzrJnd39anc/\nmeRAkkuq6uwkb+nuB7u7k3wqyXunlrljmL47yeXDldErktzX3S9290tJ7ssRGl0AmBdtJADLbkP3\nEFbV1UkOdffvr/nonCRPT70/OJSdM0yvLX/dMt19OMnLSd52jHUBwMLSRgKwTHad6AJV9eYkv5DJ\nUJiFUFXXJbkuSc4///w51waAsdJGArBsNtJD+OeSXJDk96vqqSTnJvndqvrTSQ4lOW9q3nOHskPD\n9NryTC9TVbuSnJrkhWOs6w26+5bu3t3du88888wNfCUAmAltJABL5YQDYXc/3N1/qrtXunslk2Eq\n7+ruZ5Pck2TP8FS0CzK5Mf5L3f1Mkleq6tLh3ocPJfn8sMp7kqw+He19SR4Y7qH4QpL3VNVpw43y\n7xnKAGAhaSMBWDbHHTJaVZ9NclmSM6rqYJIbuvvWI83b3Y9U1V1JHk1yOMlHu/u14eOPZPI0tlOS\n3Du8kuTWJJ+uqgOZ3Ji/Z1jXi1X1S0m+PMz3i919pBv3AWAutJEALLvjBsLu/uBxPl9Z8/7GJDce\nYb79SS4+Qvl3krz/KOu+Lcltx6sjAMyDNhKAZbehp4wCAACw/ARCAACAkRIIAQAARkogBAAAGCmB\nEAAAYKQEQgAAgJESCAEAAEbquH+HEAAAYJmt7N037yosLD2EAAAAIyUQAgAAjJRACAAAMFICIQAA\nwEgJhAAAACMlEAIAAIyUQAgAADBSAiEAAMBICYQAAAAjJRACAACMlEAIAAAwUgIhAADASAmEAAAA\nIyUQAgAAjJRACAAAMFK75l0BAACArbCyd9+8q7Dw9BACAACMlEAIAAAwUgIhAADASAmEAAAAIyUQ\nAgAAjJRACAAAMFICIQAAwEgJhAAAACMlEAIAAIyUQAgAADBSAiEAAMBIHTcQVtVtVfVcVX1tqux/\nrqo/qKp/W1W/UVVvnfrs+qo6UFWPV9UVU+XvrqqHh88+XlU1lJ9cVZ8byh+qqpWpZa6pqieG1zWz\n+tIAMAvaSACW3Xp6CG9PcuWasvuSXNzdfyHJHya5Pkmq6qIke5K8c1jmE1V10rDMJ5N8OMmFw2t1\nndcmeam7357kY0luHtZ1epIbkvxokkuS3FBVp534VwSALXN7tJEALLHjBsLu/p0kL64p++3uPjy8\nfTDJucP01Unu7O5Xu/vJJAeSXFJVZyd5S3c/2N2d5FNJ3ju1zB3D9N1JLh+ujF6R5L7ufrG7X8qk\ngV3b6ALA3GgjAVh2s7iH8G8nuXeYPifJ01OfHRzKzhmm15a/bpmhAX05yduOsS4AWBbaSAAW2qYC\nYVX9D0kOJ/nMbKqz4XpcV1X7q2r/888/P8+qAEASbSQAy2HDgbCq/maSn0zyN4YhLklyKMl5U7Od\nO5Qdyn8cMjNd/rplqmpXklOTvHCMdb1Bd9/S3bu7e/eZZ5650a8EADOhjQRgWWwoEFbVlUl+PslP\ndfd/mProniR7hqeiXZDJjfFf6u5nkrxSVZcO9z58KMnnp5ZZfTra+5I8MDSeX0jynqo6bbhR/j1D\nGQAsLG0kwPyt7N2Xlb375l2NpbDreDNU1WeTXJbkjKo6mMlTza5PcnKS+4YnYz/Y3X+nux+pqruS\nPJrJMJmPdvdrw6o+ksnT2E7J5H6K1Xsqbk3y6ao6kMmN+XuSpLtfrKpfSvLlYb5f7O7X3bgPAPOk\njQRg2R03EHb3B49QfOsx5r8xyY1HKN+f5OIjlH8nyfuPsq7bktx2vDoCwDxoIwFYdrN4yigAAABL\nSCAEAAAYKYEQAABgpARCAACAkRIIAQAARkogBAAAGCmBEAAAYKQEQgAAgJESCAEAAEZKIAQAABgp\ngRAAAGCkBEIAZmZl776s7N0372oAAOskEAIAADuGC5MnRiAEAAAYKYEQAABgpARCAACAkRIIAZg5\nD5cBgOUgEAIAAIyUQAgAADBSAiEAAMBICYQAAAAjJRACAACMlEAIAAAwUgIhAADASAmEAAAAIyUQ\nAgAAjJRACAAAMFICIQAAwEjtmncFAAAANmtl7755V2Ep6SEEAAAYKYEQAABgpARCAACAkRIIAQAA\nRkogBGAm3MwPAMvnuIGwqm6rqueq6mtTZadX1X1V9cTw87Spz66vqgNV9XhVXTFV/u6qenj47ONV\nVUP5yVX1uaH8oapamVrmmuF3PFFV18zqSwPALGgjAVh26+khvD3JlWvK9ia5v7svTHL/8D5VdVGS\nPUneOSzziao6aVjmk0k+nOTC4bW6zmuTvNTdb0/ysSQ3D+s6PckNSX40ySVJbphuVAFgAdwebSQA\nS+y4gbC7fyfJi2uKr05yxzB9R5L3TpXf2d2vdveTSQ4kuaSqzk7ylu5+sLs7yafWLLO6rruTXD5c\nGb0iyX3d/WJ3v5Tkvryx0QWAudFGArDsNnoP4Vnd/cww/WySs4bpc5I8PTXfwaHsnGF6bfnrlunu\nw0leTvK2Y6wLABaZNhKApbHph8oMVzN7BnXZsKq6rqr2V9X+559/fp5VAYDv0UYCsOg2Ggi/NQxx\nyfDzuaH8UJLzpuY7dyg7NEyvLX/dMlW1K8mpSV44xrreoLtv6e7d3b37zDPP3OBXAoCZ0EYCsDQ2\nGgjvSbL6RLNrknx+qnzP8FS0CzK5Mf5Lw9CZV6rq0uHehw+tWWZ1Xe9L8sBwRfULSd5TVacNN8q/\nZygDgEWmjQRgaew63gxV9dkklyU5o6oOZvJUs5uS3FVV1yb5RpIPJEl3P1JVdyV5NMnhJB/t7teG\nVX0kk6exnZLk3uGVJLcm+XRVHcjkxvw9w7perKpfSvLlYb5f7O61N+4DwNxoIwFYdscNhN39waN8\ndPlR5r8xyY1HKN+f5OIjlH8nyfuPsq7bktx2vDoCwDxoIwFYdscNhAAAAItqZe++eVdhqW36KaMA\nAAAsJ4EQAABgpARCAACAkRIIAQAARkogBAAAGCmBEAAAYKQEQgAAgJESCAEAAEZKIAQAABipXfOu\nAADLbWXvvnlXAQDYID2EAAAAIyUQAgAAjJRACAAAMFICIQAAwEgJhAAAACMlEAIAAIyUQAgAACwl\nf/po8wRCAACAkRIIAQAARkogBAAAGCmBEAAAYKR2zbsCACwnN/IDwPLTQwgAADBSAiEAAMBICYQA\nAAAjJRACsGXcZwgAi00gBAAAGCmBEAAAYKQEQgAAgJESCAEAAEZKIAQAABipXfOuAADLxZNDAZg3\nbdHs6CEEAAAYKYEQAABgpDYVCKvq71fVI1X1tar6bFX9QFWdXlX3VdUTw8/Tpua/vqoOVNXjVXXF\nVPm7q+rh4bOPV1UN5SdX1eeG8oeqamUz9QWA7aKNBGAZbDgQVtU5Sf6bJLu7++IkJyXZk2Rvkvu7\n+8Ik9w/vU1UXDZ+/M8mVST5RVScNq/tkkg8nuXB4XTmUX5vkpe5+e5KPJbl5o/UFgO2ijQRgWWx2\nyOiuJKdU1a4kb07y75JcneSO4fM7krx3mL46yZ3d/Wp3P5nkQJJLqursJG/p7ge7u5N8as0yq+u6\nO8nlq1dGAWDBaSMBWHgbDoTdfSjJP0ryzSTPJHm5u387yVnd/cww27NJzhqmz0ny9NQqDg5l5wzT\na8tft0x3H07ycpK3bbTOALAdtJEALIvNDBk9LZOrkxck+TNJfrCqfnp6nuFqZm+qhuury3VVtb+q\n9j///PNb/esA4Ji0kQAsi80MGf2rSZ7s7ue7+7tJfj3JjyX51jDEJcPP54b5DyU5b2r5c4eyQ8P0\n2vLXLTMMuTk1yQtrK9Ldt3T37u7efeaZZ27iKwHATGgjAVgKmwmE30xyaVW9ebhn4fIkjyW5J8k1\nwzzXJPn8MH1Pkj3DU9EuyOTG+C8NQ2deqapLh/V8aM0yq+t6X5IHhiuqAMyBPwS8btpIAJbCro0u\n2N0PVdXdSX43yeEkv5fkliQ/lOSuqro2yTeSfGCY/5GquivJo8P8H+3u14bVfSTJ7UlOSXLv8EqS\nW5N8uqoOJHkxkyewAcBC00YCsCw2HAiTpLtvSHLDmuJXM7kSeqT5b0xy4xHK9ye5+Ajl30ny/s3U\nEQDmQRsJwDLY7J+dAAAAYEkJhAAAACO1qSGjAAAA28XDzWZPDyEAAMBICYQAAAAjJRACAACMlEAI\nAAAwUgIhAADASAmEAAAAIyUQAgAAjJRACLAk5vm3l1b27vO3nwBgB/KH6QEW3HQQWxvKnrrpqu2u\nDgCwgwiEAAtqPT1yq/M8ddNVwiIAcMIEQoAFdKLDMw3nBGAn085tHYEQYIEsWoO3aPUBAGZLIARY\nAIIXADAPAiHAnAiBAMC8+bMTADvUovypiEWpBwDwRgIhwA4njAEAR2PIKMA2EcwAgEUjEAJssUUI\ngtN/rxAAYJVACDAi6w2GixBiAYCt5x5CgC0kWAEAi0wgBNgCi/5kzUWuGwCwfQwZBZihZQ5ay1x3\nAGBj9BACzMiyBapF78UEgGT52tdlo4cQYJM0VADAshIIAU7QTguAO+37AADrZ8goAADASAmEACdA\nbxoAsJMYMgpwDAIgALCTCYQAawiBAMBYGDIKAAAwUnoIgdHTIwgAi0f7vD021UNYVW+tqrur6g+q\n6rGq+ktVdXpV3VdVTww/T5ua//qqOlBVj1fVFVPl766qh4fPPl5VNZSfXFWfG8ofqqqVzdQXALaL\nNhKAZbDZIaO/kuRfdPd/muRHkjyWZG+S+7v7wiT3D+9TVRcl2ZPknUmuTPKJqjppWM8nk3w4yYXD\n68qh/NokL3X325N8LMnNm6wvQJLJVcfVF2wRbSQAC2/DgbCqTk3yV5LcmiTd/Sfd/UdJrk5yxzDb\nHUneO0xfneTO7n61u59MciDJJVV1dpK3dPeD3d1JPrVmmdV13Z3k8tUro9vBySIAGzGGNhKAnWEz\n9xBekOT5JP+0qn4kyVeS/FySs7r7mWGeZ5OcNUyfk+TBqeUPDmXfHabXlq8u83SSdPfhqno5yduS\nfHsT9QZGyMUdtpk2EoClsJlAuCvJu5L8bHc/VFW/kmHoy6ru7qrqzVRwParquiTXJcn555+/1b8O\nWBJCIHOkjQRgKWwmEB5McrC7Hxre351JY/etqjq7u58Zhro8N3x+KMl5U8ufO5QdGqbXlk8vc7Cq\ndiU5NckLayvS3bckuSVJdu/eveWNK7B4hD8WjDYSgKWw4XsIu/vZJE9X1TuGosuTPJrkniTXDGXX\nJPn8MH1Pkj3DU9EuyOTG+C8NQ2deqapLh3sfPrRmmdV1vS/JA8M9FACwsLSRACyLzf4dwp9N8pmq\nelOSryf5W5mEzLuq6tok30jygSTp7keq6q5MGsTDST7a3a8N6/lIktuTnJLk3uGVTG7G/3RVHUjy\nYiZPYANGbLUn8Kmbrnrde1hA2kgAFt6mAmF3fzXJ7iN8dPlR5r8xyY1HKN+f5OIjlH8nyfs3U0dg\n+a0NgdNlsKi0kQAbo43fXpvtIQTYMmsbBA0EAMBsCYTAQhH6AAC2j0AIzI3wBwAwXxt+yigAAADL\nTQ8hsOX0BAIALCY9hMCWEgYBABaXHkJgpgRAAIDlIRACMyMMAgAb5TxiPgwZBQAAGCk9hMCmuaIH\nALCcBEJgwwRBAIDlJhACJ0wQBADYGdxDCJwQYRAAYOfQQwgckwAIALBz6SEEAADmygXo+dFDCByR\nAzMAwM6nhxB4A2EQAGAcBEIAAICRMmQU+B49gwAA4yIQAoIgAMBICYQAAMBcuCg9f+4hhJFzIAYA\nGC+BEAAAYKQMGYWR0jMIAIBACCMjCAIAsEoghJEQBAEAWMs9hAAAACOlhxBGQO8gALBInJssDj2E\nAAAAI6WHEHYwV98AADgWgRB2IEEQAID1MGQUAABgpARCAABg2xjJtFgMGV2H1Z32qZuumnNN4Ngc\nYAEAOBF6CAEAAEZq04Gwqk6qqt+rqt8c3p9eVfdV1RPDz9Om5r2+qg5U1eNVdcVU+bur6uHhs49X\nVQ3lJ1fV54byh6pqZbP1hZ1K7yAsHm0kAItuFj2EP5fksan3e5Pc390XJrl/eJ+quijJniTvTHJl\nkk9U1UnDMp9M8uEkFw6vK4fya5O81N1vT/KxJDfPoL6wo6zs3ScMwuLSRgIMnLMspk0Fwqo6N8lV\nSX51qvjqJHcM03ckee9U+Z3d/Wp3P5nkQJJLqursJG/p7ge7u5N8as0yq+u6O8nlq1dGAWCRaSMB\nWAabfajMP07y80l+eKrsrO5+Zph+NslZw/Q5SR6cmu/gUPbdYXpt+eoyTydJdx+uqpeTvC3JtzdZ\nb1h6rrDBwtNGArDwNtxDWFU/meS57v7K0eYZrmb2Rn/HCdTluqraX1X7n3/++a3+dQBwTNpIAJbF\nZnoIfzzJT1XVTyT5gSRvqap/luRbVXV2dz8zDHV5bpj/UJLzppY/dyg7NEyvLZ9e5mBV7UpyapIX\n1laku29JckuS7N69e8sbV5gnPYOwFLSRACyFDfcQdvf13X1ud69kciP8A93900nuSXLNMNs1ST4/\nTN+TZM/wVLQLMrkx/kvD0JlXqurS4d6HD61ZZnVd7xt+h8YMgIWmjQT4jzxMZrFtxR+mvynJXVV1\nbZJvJPlAknT3I1V1V5JHkxxO8tHufm1Y5iNJbk9ySpJ7h1eS3Jrk01V1IMmLmTSqALCstJEALJSZ\nBMLu/mKSLw7TLyS5/Cjz3ZjkxiOU709y8RHKv5Pk/bOoI+wErq7B8tFGArDIZvF3CAEAAFhCWzFk\nFJgxPYMAwLJx/rIc9BACAACMlB5CWGCurAEAsJX0EMKCEgYBANhqAiEAADBTLmwvD0NGT8Dqjv3U\nTVfNuSbsZA6gAABsFz2EAMAJW9m7zwUsgB1ADyEsCCdWAABsN4EQAACYCRe4l49ACHPmwAkAwLy4\nhxDmSBgEAGCe9BACAACb4iL38hIIYQ4cNAEAWASGjAIAAIyUHkLYRnoGAQBYJALhBqzs3Zenbrpq\n3tVgiQiCAMBO5Bxn+RkyClvMgRIAgEWlhxC2iCAIAOxkznV2Bj2EAAAAI6WHEGbM1TIAAJaFQAgz\nJAwCADud852dRSCEGXBgBABgGQmEG7QaAPz5CQAAxsAF8J1JIIQNcEAEAGAn8JRROEHCIAAwNs5/\ndi6BEAAAYKQMGYV1cmUMAICdRiDcJA+X2bkEQABg7JwP7XyGjMIROPgBADAGeghhiiAIAOCcaEwE\nwhkxdHS5OegBADgnGiNDRhk9Bz4AAMZKDyGjJQgCADB2G+4hrKrzqupfVtWjVfVIVf3cUH56Vd1X\nVU8MP0+bWub6qjpQVY9X1RVT5e+uqoeHzz5eVTWUn1xVnxvKH6qqlY1/1e0hZCyulb37XvcC2Cra\nSGAZOT8ap80MGT2c5B9090VJLk3y0aq6KMneJPd394VJ7h/eZ/hsT5J3JrkyySeq6qRhXZ9M8uEk\nFw6vK4fya5O81N1vT/KxJDdvor6MmAMcsM1G00Y6vsLyc7F83DY8ZLS7n0nyzDD976vqsSTnJLk6\nyWXDbHck+WKS/34ov7O7X03yZFUdSHJJVT2V5C3d/WCSVNWnkrw3yb3DMv/jsK67k/wvVVXd3Rut\n93bwgJnF4eAGzIM2ElgGzpNIZnQP4TBM5S8meSjJWUNDmCTPJjlrmD4nyYNTix0cyr47TK8tX13m\n6STp7sNV9XKStyX59izqzc7jwAYsGm0kAIts04Gwqn4oyT9P8ve6+5Xh1oYkSXd3VW35lcqqui7J\ndUly/vnnb/WvWzc9hdtLGAQWjTYSWETOmZi2qUBYVd+fSUP3me7+9aH4W1V1dnc/U1VnJ3luKD+U\n5Lypxc8dyg4N02vLp5c5WFW7kpya5IW19ejuW5LckiS7d+9euKEyK3v3ZWXeldhhHMiARaeNBBaR\ncyjW2nAgHJ5ydmuSx7r7l6c+uifJNUluGn5+fqr816rql5P8mUxujP9Sd79WVa9U1aWZDKf5UJJ/\nsmZd/ybJ+5I84N6IcXLwApaJNhKAZbGZHsIfT/IzSR6uqq8OZb+QSSN3V1Vdm+QbST6QJN39SFXd\nleTRTJ6+9tHufm1Y7iNJbk9ySiY3yt87lN+a5NPDzfUvZvIEtqX04NffcNGW4xACgSWmjQQWivMq\njmYzTxn9v5LUUT6+/CjL3JjkxiOU709y8RHKv5Pk/Rut4yJyX+HROVABO4U2ElgUzq84npk8ZZQT\nN/2fc4zh0MEJAGDrONdivQTCBbDTew0dkAAAto9zL06EQLhA1gbDlb37li4kOgABAMDyEAgX0HSo\nOlrAmkdQFPYAOJKdPtIFFp1zNDZDIFxS/uMDAOCckM0SCAEAYMkIgsyKQAjAtjCsEGDzBEFmTSAE\nAIAFJwiyVb5v3hUAAACOThhkK+khBACABSQIsh0EQgBgJtwnCpsjADIPAiEAAMyJEMi8CYQAALDN\nBEEWhUAIAADbQAhkEQmEAMBMuZcQJgRAloE/OwEAADMmDLIs9BACAMAGCX4sO4EQANgSK3v3GTbK\njiQEspMIhAAAcBTT98QKguxEAiEAsGU8YIZlc7TQJwyyUwmEAACMkpAHAiEAsA30FDJvhn7CkQmE\nAMC2EQzZSusJesIgvJ5ACABsO08gZbMEO5gNgRAAmIu1J/QCImsJfbD1BEIAYCHoNRyP1aD37Ndf\neN17YPsJhADAwtBruPyEO1guAiEAsLCOFS6Exa2x9sE/Ah7sbAIhALCUhMU3OtI22WiwEwRhHARC\nAGDHmUWYOdEgtZ7559HrJtgBxyIQAgAcwVb0qAlnwKL5vnlXAAAAgPkQCAEAAEZKIARgW63s3WfY\nHAAsCIEQAABgpJYiEFbVlVX1eFUdqKq9864PACwC7SMAm7XwgbCqTkryvyb5a0kuSvLBqrpovrUC\ngPnSPgIwCwsfCJNckuRAd3+9u/8kyZ1Jrp5znQBg3rSPAGzaMgTCc5I8PfX+4FAGAGOmfQRg03bE\nH6avquuSXDe8/eOqenwGqz0jybdnsJ7X+cbNPznrVW6lLdkGS8Y2sA1W2Q4z3gZ180xW85/MZC07\n2Ba0kcv0f2FZ6ros9UyW6/xoWbbrstQzUdet8IZ6bnf7uAyB8FCS86benzuUfU9335Lklln+0qra\n3927Z7nOZWMb2AaJbbDKdrANFtBx28dk9m3kMu0Hy1LXZalnoq5bYVnqmajrVliEei7DkNEvJ7mw\nqi6oqjcl2ZPknjnXCQDmTfsIwKYtfA9hdx+uqr+b5AtJTkpyW3c/MudqAcBcaR8BmIWFD4RJ0t2/\nleS3tvnXznQI6pKyDWyDxDZYZTvYBgtH+3hcy1LXZalnoq5bYVnqmajrVph7Pau7510HAAAA5mAZ\n7iEEAABgCwiEa1TVlVX1eFUdqKq9867Pdqiq86rqX1bVo1X1SFX93FB+elXdV1VPDD9Pm3ddt1pV\nnVRVv1dVvzm8H+M2eGtV3V1Vf1BVj1XVXxrbdqiqvz/8X/haVX22qn5gp2+Dqrqtqp6rqq9NlR31\nO1fV9cNx8vGqumI+tWbWjtcG1sTHh8//bVW9a73LbnM9/8ZQv4er6l9X1Y9MffbUUP7Vqtq/lfVc\nZ10vq6qXh/p8tar+4XqXnUNd/7upen6tql6rqtOHz7Ztux7peLXm80XZT49Xz0XaT49X14XYT9dR\nz4XYR4ffd8Tz6zXzLMS+mu72Gl6Z3JT/fyf5s0nelOT3k1w073ptw/c+O8m7hukfTvKHSS5K8j8l\n2TuU701y87zrug3b4r9N8mtJfnN4P8ZtcEeS/2qYflOSt45pO2Tyh72fTHLK8P6uJH9zp2+DJH8l\nybuSfG2q7IjfeTg+/H6Sk5NcMBw3T5r3d/Da9D5w3DYwyU8kuTdJJbk0yUPrXXab6/ljSU4bpv/a\naj2H908lOWOBtullq23OiS673XVdM/9fT/LAnLbrG45Xi7afrrOeC7GfrrOui7KfHrOei7KPDr/v\niOfXi7iv6iF8vUuSHOjur3f3nyS5M8nVc67TluvuZ7r7d4fpf5/ksUxOiq/OJBxk+Pne+dRwe1TV\nuUmuSvKrU8Vj2wanZnKwvTVJuvtPuvuPMrLtkMkDt06pql1J3pzk32WHb4Pu/p0kL64pPtp3vjrJ\nnd39anfgPRmcAAAEuElEQVQ/meRAJsdPltt62sCrk3yqJx5M8taqOnudy25bPbv7X3f3S8PbBzP5\nG43zsJntst3nJCf6+z6Y5LNbWJ+jOsrxatoi7KfHrecC7afr2aZHs1DbdI257aPJMc+vpy3EvioQ\nvt45SZ6een8wb/yH29GqaiXJX0zyUJKzuvuZ4aNnk5w1p2ptl3+c5OeT/H9TZWPbBhckeT7JP63J\n0NlfraofzIi2Q3cfSvKPknwzyTNJXu7u386ItsGUo33n0R8rd6j1/LsebZ7t3CdO9Hddm8kV+FWd\n5P+sqq9U1XVbUL9p663rjw3Dxe6tqnee4LKzsu7fV1VvTnJlkn8+Vbyd2/V4FmE/PVHz3E/XaxH2\n03VZtH10zfn1tIXYV5fiz06wParqhzL5j/P3uvuVqvreZ93dVbVjH0lbVT+Z5Lnu/kpVXXakeXb6\nNhjsymQoxs9290NV9SuZDBX8np2+HWpyn9zVmYTjP0ryv1fVT0/Ps9O3wZGM8Tuz/KrqP8/kRPsv\nTxX/5e4+VFV/Ksl9VfUHQ6/DvPxukvO7+4+r6ieS/B9JLpxjfdbjryf5V9093VOzaNt1adhPt8TC\n7KNrz6+3+vdthB7C1zuU5Lyp9+cOZTteVX1/JjvrZ7r714fibw3d1hl+Pjev+m2DH0/yU1X1VCbd\n8v9FVf2zjGsbJJMrUAe7e/UK1t2ZBMQxbYe/muTJ7n6+u7+b5Nczuc9jTNtg1dG+82iPlTvcev5d\njzbPdu4T6/pdVfUXMrkF4OrufmG1fBgFkO5+LslvZGuHOx+3rt39Snf/8TD9W0m+v6rOWM+y213X\nKXuyZijeNm/X41mE/XRdFmQ/Pa4F2k/XayH20aOcX09biH1VIHy9Lye5sKouqKo3ZbIz3TPnOm25\nmnQF3prkse7+5amP7klyzTB9TZLPb3fdtkt3X9/d53b3Sib/7g90909nRNsgSbr72SRPV9U7hqLL\nkzyacW2Hbya5tKrePPzfuDyTcf9j2garjvad70myp6pOrqoLMrlK/KU51I/ZWk8beE+SDw1Pxrs0\nkyHVz6xz2W2rZ1Wdn8nFnJ/p7j+cKv/Bqvrh1ekk70lyxKcVbmNd//RwrElVXZLJudkL61l2u+s6\n1PHUJP9Zpo6Bc9iux7MI++lxLdB+elwLtJ8e16Lso8c4v562EPuqIaNTuvtwVf3dJF/I5Ok+t3X3\nI3Ou1nb48SQ/k+ThqvrqUPYLSW5KcldVXZvkG0k+MKf6zdMYt8HPJvnMcAD6epK/lcmBfxTbYRgq\ne3cmw2MOJ/m9JLck+aHs4G1QVZ/N5ClyZ1TVwSQ35Cj7f3c/UlV3ZXKx4HCSj3b3a3OpODNztDaw\nqv7O8Pn/luS3Mnkq3oEk/yGT48O2tp/rrOc/TPK2JJ8YzmEPd/fuTO6D/Y2hbFeSX+vuf7EV9TyB\nur4vyX9dVYeT/L9J9nR3J9nWc5J11jVJ/sskv93d/8/U4tu6XY9yvPr+qXrOfT9dZz0XYj9dZ10X\nYj9dRz2TBdhHB0c7vz5/qr6Lsa9O/i0BAAAYG0NGAQAARkogBAAAGCmBEAAAYKQEQgAAgJESCAEA\nAEZKIAQAABgpgRAAAGCkBEIAAICR+v8Bq/V1hQZJT60AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f28a861d358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (15,6)\n",
    "subplot(1,2,1)\n",
    "hist(wise[~np.isnan(wise[\"lr_pc\"])][\"lr_pc\"], bins=200)\n",
    "vlines([threshold], 0, 150000)\n",
    "ylim([0,150000])\n",
    "subplot(1,2,2)\n",
    "hist(np.log10(wise[~np.isnan(wise[\"lr_pc\"])][\"lr_pc\"]+1), bins=200)\n",
    "vlines(np.log10(threshold+1), 0, 150000)\n",
    "ylim([0,150000]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combine data\n",
    "\n",
    "First we discard the matches below the threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wise[\"lr_index\"] = wise[\"lr_panstarrs_index\"]\n",
    "wise[\"lr_index\"][wise[\"lr_pc\"] < threshold] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We combine the two catalogues using an outer join wich maintains all the data for the two catalogues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "panstarrs[\"lr_index\"] = np.arange(len(panstarrs)).astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "combined = join(wise, panstarrs, join_type='outer', keys='lr_index')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "combined['ra'] = combined['raMean']\n",
    "combined['dec'] = combined['decMean']\n",
    "combined['raErr'] = combined['raMeanErr']\n",
    "combined['decErr'] = combined['decMeanErr']\n",
    "\n",
    "combined['ra'][np.isnan(combined['raMean'])] = combined['raWise'][np.isnan(combined['raMean'])]\n",
    "combined['dec'][np.isnan(combined['decMean'])] = combined['decWise'][np.isnan(combined['decMean'])]\n",
    "combined['raErr'][np.isnan(combined['raMean'])] = combined['raWiseErr'][np.isnan(combined['raMean'])]\n",
    "combined['decErr'][np.isnan(combined['decMean'])] = combined['decWiseErr'][np.isnan(combined['decMean'])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Important step to solve problems with the default values for the columns. Check later why we have to do that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for col in [\"raMean\", \"decMean\", \"raMeanErr\", \"decMeanErr\",\n",
    "            \"raWise\", \"decWise\", \"raWiseErr\", \"decWiseErr\",\n",
    "            \"ra\", \"dec\", \"raErr\", \"decErr\"]:\n",
    "    combined[col].fill_value = 1e+20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "columns_save = ['AllWISE', 'objID', 'ra', 'dec', 'raErr', 'decErr', \n",
    "                'W1mag', 'W1magErr', 'i', 'iErr']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "combined[columns_save].write('pw.fits', format=\"fits\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also save a version with all the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['AllWISE', 'raWise', 'decWise', 'raWiseErr', 'decWiseErr', 'W1mag',\n",
       "       'W1magErr', 'lr', 'lr_dist', 'lr_panstarrs_index', 'lr_pc',\n",
       "       'lr_index', 'objID', 'raMean', 'decMean', 'raMeanErr', 'decMeanErr',\n",
       "       'i', 'iErr', 'ra', 'dec', 'raErr', 'decErr'],\n",
       "      dtype='<U18')"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(combined.colnames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "combined.write('pw_lrdata.fits', format=\"fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}