build_kepler_catalogue.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# IMPORT STATEMENTS AND PRELIMINARIES\n",
    "\n",
    "import numpy as np\n",
    "import csv\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "\n",
    "pi = np.pi\n",
    "\n",
    "MAINPATH = '/Users/research/projects/maiasaurus/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in CKS and DR25 catalogues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# file names\n",
    "cks_file = MAINPATH + 'catalogues/cks_physical_merged.csv'\n",
    "dr25_file = MAINPATH + 'catalogues/kepler_dr25.csv'\n",
    "\n",
    "# convenience function to read in csv file\n",
    "def read_csv_file(filename):\n",
    "    data = []\n",
    "    with open(filename) as infile:\n",
    "        reader = csv.reader(infile)\n",
    "\n",
    "        for row in reader:\n",
    "            data.append(row)\n",
    "\n",
    "        if filename == cks_file:\n",
    "            keys   = data[0]\n",
    "            values = data[1:]\n",
    "\n",
    "        if filename == dr25_file:\n",
    "            keys   = data[0]\n",
    "            values = data[1:]\n",
    "\n",
    "            \n",
    "        return keys, values\n",
    "\n",
    "\n",
    "# READ IN CKS DATABASE --  http://www.astro.caltech.edu/~howard/cks/column-definitions.txt\n",
    "cks_keys, cks_data = read_csv_file(cks_file)\n",
    "\n",
    "# READ IN DR25 DATABASE -- https://exoplanetarchive.ipac.caltech.edu\n",
    "dr25_keys, dr25_data = read_csv_file(dr25_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convenience functions to pull data from csv files\n",
    "def getcks(keyname,keys=cks_keys,data=cks_data):\n",
    "    '''\n",
    "    keyname = (string) of column definition, see CKS documentation\n",
    "    '''\n",
    "    kid = keys.index(keyname)\n",
    "    \n",
    "    outdata = []\n",
    "    for row in data:\n",
    "        outdata.append(row[kid])\n",
    "    \n",
    "    return outdata\n",
    "\n",
    "def getdr25(keyname,keys=dr25_keys,data=dr25_data):\n",
    "    '''\n",
    "    keyname = (string) of column definition, see CKS documentation\n",
    "    '''\n",
    "    kid = keys.index(keyname)\n",
    "    \n",
    "    outdata = []\n",
    "    for row in data:\n",
    "        outdata.append(row[kid])\n",
    "    \n",
    "    return outdata\n",
    "\n",
    "def printkeys(data):\n",
    "    for k in data.keys():\n",
    "        print(k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read cks data into a dictionary\n",
    "cks = {}\n",
    "for k in cks_keys: \n",
    "    cks[k] = getcks(k)\n",
    "    \n",
    "# read dr25 data into a dictionary\n",
    "dr25 = {}\n",
    "for k in dr25_keys:\n",
    "    dr25[k] = getdr25(k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total number of queried CKS objects = 2025\n",
      "total number of queried DR25 objects = 9564\n"
     ]
    }
   ],
   "source": [
    "def check_lengths(data):\n",
    "    keys = data.keys()\n",
    "    k0   = list(keys)[0]\n",
    "    L0   = len(data[k0])\n",
    "    \n",
    "    for k in keys:\n",
    "        if len(data[k]) != L0:\n",
    "            raise ValueError('inconsistent array lengths')\n",
    "            \n",
    "    return None\n",
    "\n",
    "\n",
    "def convert_to_arrays(data):\n",
    "    keys = data.keys()\n",
    "    dnew = {}\n",
    "    \n",
    "    for k in keys:\n",
    "        dnew[k] = np.asarray(data[k])\n",
    "        \n",
    "    return dnew       \n",
    "\n",
    "\n",
    "# do some quick cleanup\n",
    "check_lengths(cks)\n",
    "check_lengths(dr25)\n",
    "\n",
    "cks  = convert_to_arrays(cks)\n",
    "dr25 = convert_to_arrays(dr25)\n",
    "\n",
    "# grab a reference key\n",
    "cks_k0  = list(cks.keys())[0]\n",
    "dr25_k0 = list(dr25.keys())[0]\n",
    "\n",
    "print('total number of queried CKS objects =', len(cks[cks_k0]))\n",
    "print('total number of queried DR25 objects =', len(dr25[dr25_k0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in Gaia DR2 constrained values (from Andrew)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the number of stars in the Gaia/CKS sample pre-cut is 1189\n",
      "the number of planets in the Gaia/CKS sample pre-cut is 1901\n"
     ]
    }
   ],
   "source": [
    "# load in Andrew's catalogues\n",
    "starsfile  = MAINPATH + 'catalogues/gaia_cks_stars_may14.tsv'\n",
    "planetfile = MAINPATH + 'catalogues/gaia_cks_planets_may14.tsv'\n",
    "\n",
    "gaia_stars    = np.genfromtxt(starsfile, skip_header=89, delimiter=';', usecols=np.arange(1,31))\n",
    "gaia_stars_id = np.genfromtxt(starsfile, skip_header=89, delimiter=';', usecols=[0], dtype='str')\n",
    "\n",
    "gaia_planets    = np.genfromtxt(planetfile, skip_header=64, delimiter=';', usecols=np.arange(1,16))\n",
    "gaia_planets_id = np.genfromtxt(planetfile, skip_header=64, delimiter=';', usecols=[0], dtype='str')\n",
    "\n",
    "\n",
    "print('the number of stars in the Gaia/CKS sample pre-cut is %d' %len(gaia_stars[:,0]))\n",
    "print('the number of planets in the Gaia/CKS sample pre-cut is %d' %len(gaia_planets[:,0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1189 1189\n",
      "1901 1901\n",
      "1901 2025\n"
     ]
    }
   ],
   "source": [
    "# ONLY WORRY ABOUT OVERLAPPING SYSTEMS\n",
    "\n",
    "# first clean up Andrew's Gaia catalogue\n",
    "keep = np.isin(gaia_stars_id, np.unique(cks['id_starname']))\n",
    "print(np.sum(keep), len(keep))\n",
    "\n",
    "gaia_stars    = gaia_stars[keep]\n",
    "gaia_stars_id = gaia_stars_id[keep]\n",
    "\n",
    "keep = np.isin(gaia_planets_id, cks['id_koicand'])\n",
    "print(np.sum(keep), len(keep))\n",
    "\n",
    "gaia_planets    = gaia_planets[keep]\n",
    "gaia_planets_id = gaia_planets_id[keep]\n",
    "\n",
    "\n",
    "# then clean up my CKS catalogue\n",
    "keep = np.isin(cks['id_koicand'], gaia_planets_id)\n",
    "print(np.sum(keep), len(keep))\n",
    "\n",
    "for k in cks.keys():\n",
    "    cks[k] = cks[k][keep]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rename variables (it's sloppy, but I'm too lazy to change it right now)\n",
    "cks_stars = gaia_stars\n",
    "cks_stars_id = gaia_stars_id\n",
    "cks_planets = gaia_planets\n",
    "cks_planets_id = gaia_planets_id\n",
    "n_cks = len(cks_planets_id)\n",
    "\n",
    "\n",
    "# make arrays to hold stellar properties\n",
    "cks_srad_raw = cks_stars[:,10]\n",
    "cks_srad_err1_raw = cks_stars[:,11]\n",
    "cks_srad_err2_raw = cks_stars[:,12]\n",
    "cks_smass_raw = cks_stars[:,13]\n",
    "cks_smass_err1_raw = cks_stars[:,14]\n",
    "cks_smass_err2_raw = cks_stars[:,15]\n",
    "cks_steff_raw = cks_stars[:,0]\n",
    "cks_steff_err1_raw = cks_stars[:,1]\n",
    "cks_steff_err2_raw = cks_stars[:,2]\n",
    "\n",
    "cks_dilution_raw = cks_stars[:,28]\n",
    "cks_rcf_raw = cks_stars[:,29]\n",
    "\n",
    "cks_parallax1_raw = cks_stars[:,8]\n",
    "cks_parallax1_err0_raw = cks_stars[:,9]\n",
    "cks_parallax2_raw = cks_stars[:,25]\n",
    "cks_parallax2_err1_raw = cks_stars[:,26]\n",
    "cks_parallax2_err2_raw = cks_stars[:,27]\n",
    "\n",
    "cks_srad = np.zeros(n_cks)\n",
    "cks_srad_err1 = np.zeros(n_cks)\n",
    "cks_srad_err2 = np.zeros(n_cks)\n",
    "cks_smass = np.zeros(n_cks)\n",
    "cks_smass_err1 = np.zeros(n_cks)\n",
    "cks_smass_err2 = np.zeros(n_cks)\n",
    "cks_steff = np.zeros(n_cks)\n",
    "cks_steff_err1 = np.zeros(n_cks)\n",
    "cks_steff_err2 = np.zeros(n_cks)\n",
    "\n",
    "cks_dilution = np.zeros(n_cks)\n",
    "cks_rcf = np.zeros(n_cks)\n",
    "cks_parallax1 = np.zeros(n_cks)\n",
    "cks_parallax1_err0 = np.zeros(n_cks)\n",
    "cks_parallax2 = np.zeros(n_cks)\n",
    "cks_parallax2_err1 = np.zeros(n_cks)\n",
    "cks_parallax2_err2 = np.zeros(n_cks)\n",
    "\n",
    "\n",
    "# get stellar properties for planets from CKS stellar catalog\n",
    "for i in range(n_cks):\n",
    "    idx = np.where(cks_stars_id == cks_planets_id[i][:-3])[0]\n",
    "    \n",
    "    cks_srad[i] = cks_srad_raw[idx]\n",
    "    cks_srad_err1[i] = cks_srad_err1_raw[idx]\n",
    "    cks_srad_err2[i] = cks_srad_err2_raw[idx]\n",
    "    cks_smass[i] = cks_smass_raw[idx]\n",
    "    cks_smass_err1[i] = cks_smass_err1_raw[idx]\n",
    "    cks_smass_err2[i] = cks_smass_err2_raw[idx]\n",
    "    cks_steff[i] = cks_steff_raw[idx]\n",
    "    cks_steff_err1[i] = cks_steff_err1_raw[idx]\n",
    "    cks_steff_err2[i] = cks_steff_err2_raw[idx]\n",
    "   \n",
    "    cks_dilution[i] = cks_dilution_raw[idx]\n",
    "    cks_rcf[i] = cks_rcf_raw[idx]\n",
    "    if np.isnan(cks_rcf[i]): cks_rcf[i] = 1.0\n",
    "    \n",
    "    cks_parallax1[i] = cks_parallax1_raw[idx]\n",
    "    cks_parallax1_err0[i] = cks_parallax1_err0_raw[idx]\n",
    "    cks_parallax2[i] = cks_parallax2_raw[idx]\n",
    "    cks_parallax2_err1[i] = cks_parallax2_err1_raw[idx]\n",
    "    cks_parallax2_err2[i] = cks_parallax2_err2_raw[idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reorder my CKS catalogue to match Andrew's\n",
    "pos = [list(cks['id_koicand']).index(pid) for pid in cks_planets_id]\n",
    "\n",
    "for k in cks.keys():\n",
    "    cks[k] = cks[k][pos]\n",
    "    \n",
    "    \n",
    "# and check that it worked properly\n",
    "for i in range(1901):\n",
    "    if cks_planets_id[i] != cks['id_koicand'][i]:\n",
    "        print(cks_planets_id[i], cks['id_koicand'][i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# update stellar properties\n",
    "cks['iso_srad']  = cks_srad\n",
    "cks['iso_srad_err1']  = cks_srad_err1\n",
    "cks['iso_srad_err2']  = cks_srad_err2\n",
    "cks['iso_smass'] = cks_smass\n",
    "cks['iso_smass_err1'] = cks_smass_err1\n",
    "cks['iso_smass_err2'] = cks_smass_err2\n",
    "cks['iso_steff'] = cks_steff\n",
    "cks['iso_steff_err1'] = cks_steff_err1\n",
    "cks['iso_steff_err2'] = cks_steff_err2\n",
    "\n",
    "\n",
    "# track dilution and radius correlation factor (RCF)\n",
    "cks['cks_dilution'] = np.array([1.0]*len(cks[cks_k0]))\n",
    "cks['cks_rcf'] = np.array([1.0]*len(cks[cks_k0]))\n",
    "\n",
    "cks['cks_dilution'] = cks_dilution\n",
    "cks['cks_rcf'] = cks_rcf\n",
    "\n",
    "\n",
    "# record paralax info\n",
    "cks['parallax1'] = np.array([None]*len(cks[cks_k0]))\n",
    "cks['parallax1_err0'] = np.array([None]*len(cks[cks_k0]))\n",
    "cks['parallax2'] = np.array([None]*len(cks[cks_k0]))\n",
    "cks['parallax2_err1'] = np.array([None]*len(cks[cks_k0]))\n",
    "cks['parallax2_err2'] = np.array([None]*len(cks[cks_k0]))\n",
    "\n",
    "cks['parallax1'] = cks_parallax1\n",
    "cks['parallax1_err0'] = cks_parallax1_err0\n",
    "cks['parallax2'] = cks_parallax2\n",
    "cks['parallax2_err1'] = cks_parallax2_err1\n",
    "cks['parallax2_err2'] = cks_parallax2_err2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# update planet radii\n",
    "cks['iso_prad'] = cks_planets[:,6]\n",
    "cks['iso_prad_err1'] = cks_planets[:,7]\n",
    "cks['iso_prad_err2'] = cks_planets[:,8]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Remove false-positives, low SNR objects, and grazing transits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "after removing false positives, CKS objects = 1752\n",
      "after removing false positives, DR25 objects = 4541\n"
     ]
    }
   ],
   "source": [
    "# remove false positives (and not dispositioned objects) from CKS\n",
    "fp = cks['koi_disposition'] == 'FALSE POSITIVE'\n",
    "fp += cks['cks_fp'] == 'True'\n",
    "\n",
    "for k in cks.keys():\n",
    "    cks[k] = cks[k][~fp]\n",
    "\n",
    "# remove false positives from DR25\n",
    "fp = dr25['koi_disposition'] == 'FALSE POSITIVE'\n",
    "fp += (dr25['koi_disposition'] != 'CONFIRMED')*(dr25['koi_pdisposition'] == 'FALSE POSITIVE')\n",
    "\n",
    "for k in dr25.keys():\n",
    "    dr25[k] = dr25[k][~fp]\n",
    "\n",
    "print('after removing false positives, CKS objects =', len(cks[cks_k0]))\n",
    "print('after removing false positives, DR25 objects =', len(dr25[dr25_k0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "after cutting on SNR > 7.1, CKS objects = 1718\n",
      "after cutting on SNR > 7.1, DR25 objects = 4341\n"
     ]
    }
   ],
   "source": [
    "# cut CKS low signal-to-noise objects (SNR < 7.1)\n",
    "cks_snr = cks['koi_model_snr']\n",
    "cks_snr[cks_snr == ''] = 0\n",
    "cks_snr = np.array(cks_snr, dtype='float')\n",
    "\n",
    "for k in cks.keys():\n",
    "    cks[k] = cks[k][cks_snr >= 7.1]\n",
    "    \n",
    "\n",
    "# cut DR25 low signal-to-noise objects (SNR < 7.1)\n",
    "dr25_snr = dr25['koi_model_snr']\n",
    "dr25_snr[dr25_snr == ''] = 0\n",
    "dr25_snr = np.array(dr25_snr, dtype='float')\n",
    "\n",
    "for k in dr25.keys():\n",
    "    dr25[k] = dr25[k][dr25_snr >= 7.1]\n",
    "    \n",
    "    \n",
    "print('after cutting on SNR > 7.1, CKS objects =', len(cks[cks_k0]))\n",
    "print('after cutting on SNR > 7.1, DR25 objects =', len(dr25[dr25_k0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "after removing grazing transits, objects = 1673\n"
     ]
    }
   ],
   "source": [
    "# remove high impact parameter systems\n",
    "b    = np.array(cks['koi_impact'], dtype='float')\n",
    "rprs = np.array(cks['koi_ror'], dtype='float')\n",
    "\n",
    "for k in cks.keys():\n",
    "    cks[k] = cks[k][b < 1-rprs]\n",
    "\n",
    "    \n",
    "print('after removing grazing transits, objects =', len(cks[cks_k0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make cuts based on stellar properties (Fulton 2018 - CKS VII)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "after cutting on stellar temperature and radius, CKS objects = 1563\n"
     ]
    }
   ],
   "source": [
    "cks_steff = np.array(cks['iso_steff'], dtype='float')\n",
    "cks_srad  = np.array(cks['iso_srad'], dtype='float')\n",
    "\n",
    "keep = (cks_steff > 4700) * (cks_steff < 6500) * (np.log10(cks_srad) < ((cks_steff - 5500.0)/4000.0 + 0.2))\n",
    "\n",
    "for k in cks.keys():\n",
    "    cks[k] = cks[k][keep]\n",
    "\n",
    "print('after cutting on stellar temperature and radius, CKS objects =', len(cks[cks_k0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "after accounting for dilution factors, CKS objects = 1419\n"
     ]
    }
   ],
   "source": [
    "dilution = np.array(cks['cks_dilution'], dtype='float')\n",
    "dilution[np.isnan(dilution)] = 99\n",
    "\n",
    "rcf = np.array(cks['cks_rcf'], dtype='float')\n",
    "rcf[np.isnan(rcf)] = 99\n",
    "\n",
    "keep = (dilution < 1.1) * (rcf < 1.05)\n",
    "\n",
    "for k in cks.keys():\n",
    "    cks[k] = cks[k][keep]\n",
    "\n",
    "print('after accounting for dilution factors, CKS objects =', len(cks[cks_k0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "after cutting stars with parallax discrepancies = CKS objects 1389\n"
     ]
    }
   ],
   "source": [
    "parallax1 = np.array(cks['parallax1'], dtype='float')\n",
    "parallax1_err0 = np.array(cks['parallax1_err0'], dtype='float')\n",
    "parallax2 = np.array(cks['parallax2'], dtype='float')\n",
    "parallax2_err1 = np.array(cks['parallax2_err1'], dtype='float')\n",
    "parallax2_err2 = np.array(cks['parallax2_err2'], dtype='float')\n",
    "\n",
    "numerator   = parallax1-parallax2\n",
    "denomenator = np.sqrt(parallax1_err0**2 + (0.5*(parallax2_err1+parallax2_err2))**2)\n",
    "\n",
    "keep = np.abs(numerator/denomenator) < 4\n",
    "keep[np.isnan(parallax1)] = False\n",
    "\n",
    "for k in cks.keys():\n",
    "    cks[k] = cks[k][keep]\n",
    "\n",
    "print('after cutting stars with parallax discrepancies = CKS objects', len(cks[cks_k0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adjust number of planets in each CKS system to reflect removals \n",
    "cks_koi = cks['id_starname']\n",
    "unique_cks_systems = np.unique(cks_koi)\n",
    "\n",
    "for us in unique_cks_systems:\n",
    "    reduced_npl = np.sum(cks_koi == us)\n",
    "    cks['koi_count'][cks_koi == us] = reduced_npl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([517., 414., 240., 144.,  55.,  19.]),\n",
       " array([1, 2, 3, 4, 5, 6, 7]),\n",
       " <a list of 6 Patch objects>)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADvpJREFUeJzt3VGsXVWdx/HvzxbRQcdqudOU3kpJJE7IJAJpGAzGODQaUGN5UIOZ0YZAOg91gnESB30xJvOgL+KYEJKmdabMqMigBEKIIwGM4wNoiwgCGispae8ArQgoY8QU//NwF5krFu6595zrvmf5/SQnZ6211zn7v9vkd/ddd599UlVIkvr1iqELkCStLINekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1Lm1QxcAcOqpp9aWLVuGLkOSpsqBAwd+XlUzi81bFUG/ZcsW9u/fP3QZkjRVkjw6yjyXbiSpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXNTH/SbNm8hSRePTZu3DP3PKalDq+IWCOP4nyOPcsU1Tw5dxkTs2bV+6BIkdWjqz+glSS9vpKBPcijJA0nuS7K/jb0hye1JftqeX9/Gk+SLSQ4muT/JuSt5AJKkl7eUM/q/qaqzq2pr618F3FFVZwJ3tD7AxcCZ7bETuHZSxUqSlm6cpZvtwL7W3gdcsmD8upp3N7AuycYx9iNJGsOoQV/At5IcSLKzjW2oqsda+3FgQ2tvAg4veO2RNiZJGsCoV928rarmkvwFcHuSHy/cWFWVpJay4/YDYyfAG9/4xqW8VJK0BCOd0VfVXHs+CtwEnAc88cKSTHs+2qbPAZsXvHy2jb34PXdX1daq2jozs+g3YUmSlmnRoE9ySpLXvtAG3gX8CLgF2NGm7QBubu1bgI+0q2/OB55ZsMQjSfojG2XpZgNwU5IX5n+lqr6Z5PvADUkuBx4FPtjm3wa8GzgI/Bq4bOJVS5JGtmjQV9UjwFtOMP4ksO0E4wXsmkh1kqSx+clYSeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOrR26AP2/NWtPJsnQZUzEabOnM3f40NBlSMKgX1WeP/4cV1zz5NBlTMSeXeuHLkFSM/LSTZI1SX6Q5NbWPyPJPUkOJvlakle28ZNb/2DbvmVlSpckjWIpa/RXAg8v6H8OuLqq3gQ8BVzexi8HnmrjV7d5kqSBjBT0SWaB9wB7Wj/AhcCNbco+4JLW3t76tO3b0svCsyRNoVHP6L8AfAL4XeuvB56uquOtfwTY1NqbgMMAbfszbb4kaQCLBn2S9wJHq+rAJHecZGeS/Un2Hzt2bJJvLUlaYJQz+guA9yU5BFzP/JLNvwDrkrxw1c4sMNfac8BmgLb9dcAfXEpSVburamtVbZ2ZmRnrICRJL23RoK+qT1bVbFVtAS4F7qyqvwXuAt7fpu0Abm7tW1qftv3OqqqJVi1JGtk4n4z9J+DjSQ4yvwa/t43vBda38Y8DV41XoiRpHEv6wFRVfRv4dms/Apx3gjm/AT4wgdokSRPgvW4kqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SercokGf5FVJvpfkh0keTPKZNn5GknuSHEzytSSvbOMnt/7Btn3Lyh6CJOnljHJG/xxwYVW9BTgbuCjJ+cDngKur6k3AU8Dlbf7lwFNt/Oo2T5I0kEWDvuY927ontUcBFwI3tvF9wCWtvb31adu3JcnEKpYkLclIa/RJ1iS5DzgK3A78DHi6qo63KUeATa29CTgM0LY/A6yfZNGSpNGNFPRV9XxVnQ3MAucBfznujpPsTLI/yf5jx46N+3aSpJewpKtuqupp4C7grcC6JGvbpllgrrXngM0AbfvrgCdP8F67q2prVW2dmZlZZvmSpMWMctXNTJJ1rf1q4J3Aw8wH/vvbtB3Aza19S+vTtt9ZVTXJoiVJo1u7+BQ2AvuSrGH+B8MNVXVrkoeA65P8M/ADYG+bvxf49yQHgV8Al65A3ZKkES0a9FV1P3DOCcYfYX69/sXjvwE+MJHqJElj85OxktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUudG+XJwacnWrD2ZJEOXMbbTZk9n7vChocuQxmLQa0U8f/w5rrjmyaHLGNueXeuHLkEam0s3ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjq3aNAn2ZzkriQPJXkwyZVt/A1Jbk/y0/b8+jaeJF9McjDJ/UnOXemDkCS9tFHO6I8D/1hVZwHnA7uSnAVcBdxRVWcCd7Q+wMXAme2xE7h24lVLkka2aNBX1WNVdW9r/wp4GNgEbAf2tWn7gEtaeztwXc27G1iXZOPEK5ckjWRJa/RJtgDnAPcAG6rqsbbpcWBDa28CDi942ZE2JkkawMhBn+Q1wNeBj1XVLxduq6oCaik7TrIzyf4k+48dO7aUl0qSlmCkoE9yEvMh/+Wq+kYbfuKFJZn2fLSNzwGbF7x8to39nqraXVVbq2rrzMzMcuuXJC1ilKtuAuwFHq6qzy/YdAuwo7V3ADcvGP9Iu/rmfOCZBUs8kqQ/slG+YeoC4MPAA0nua2OfAj4L3JDkcuBR4INt223Au4GDwK+ByyZasSRpSRYN+qr6LvBSX/657QTzC9g1Zl2SpAnxk7GS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUuVG+SlD6k7Vm7cnMf23y9Dtt9nTmDh8augwNwKCXXsbzx5/jimueHLqMidiza/3QJWggLt1IUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXOLBn2SLyU5muRHC8bekOT2JD9tz69v40nyxSQHk9yf5NyVLF6StLhRzuj/DbjoRWNXAXdU1ZnAHa0PcDFwZnvsBK6dTJmSpOVaNOir6jvAL140vB3Y19r7gEsWjF9X8+4G1iXZOKliJUlLt9w1+g1V9VhrPw5saO1NwOEF8460MUnSQMb+Y2xVFVBLfV2SnUn2J9l/7NixccuQJL2E5Qb9Ey8sybTno218Dti8YN5sG/sDVbW7qrZW1daZmZllliFJWsxyg/4WYEdr7wBuXjD+kXb1zfnAMwuWeCRJA1j0O2OTfBV4B3BqkiPAp4HPAjckuRx4FPhgm34b8G7gIPBr4LIVqFmStASLBn1VfeglNm07wdwCdo1blCRpcvxkrCR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucWvXulpD6sWXsySYYuYyJOmz2ducOHhi5jahj00p+I548/xxXXPDl0GROxZ9f6oUuYKi7dSFLnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalz3tRM0tTxTpxLY9BLmjreiXNpXLqRpM4Z9JLUuRUJ+iQXJflJkoNJrlqJfUiSRjPxoE+yBrgGuBg4C/hQkrMmvR9J0mhW4oz+POBgVT1SVb8Frge2r8B+JEkjWImg3wQcXtA/0sYkSQNIVU32DZP3AxdV1RWt/2Hgr6vqoy+atxPY2bpvBn6yzF2eCvx8ma9dbTyW1aeX4wCPZbUa51hOr6qZxSatxHX0c8DmBf3ZNvZ7qmo3sHvcnSXZX1Vbx32f1cBjWX16OQ7wWFarP8axrMTSzfeBM5OckeSVwKXALSuwH0nSCCZ+Rl9Vx5N8FPgvYA3wpap6cNL7kSSNZkVugVBVtwG3rcR7n8DYyz+riMey+vRyHOCxrFYrfiwT/2OsJGl18RYIktS5qQ36JF9KcjTJj4auZVxJNie5K8lDSR5McuXQNS1Hklcl+V6SH7bj+MzQNY0ryZokP0hy69C1jCPJoSQPJLkvyf6h61muJOuS3Jjkx0keTvLWoWtajiRvbv8XLzx+meRjK7a/aV26SfJ24Fnguqr6q6HrGUeSjcDGqro3yWuBA8AlVfXQwKUtSeZvEH5KVT2b5CTgu8CVVXX3wKUtW5KPA1uBP6+q9w5dz3IlOQRsraqpvvY8yT7gv6tqT7uq78+q6umh6xpHu23MHPOfN3p0JfYxtWf0VfUd4BdD1zEJVfVYVd3b2r8CHmYKP01c855t3ZPaYzrPJIAks8B7gD1D1yJI8jrg7cBegKr67bSHfLMN+NlKhTxMcdD3KskW4BzgnmErWZ621HEfcBS4vaqm8jiaLwCfAH43dCETUMC3khxon0qfRmcAx4B/bctpe5KcMnRRE3Ap8NWV3IFBv4okeQ3wdeBjVfXLoetZjqp6vqrOZv4T0eclmcpltSTvBY5W1YGha5mQt1XVuczfVXZXW/qcNmuBc4Frq+oc4H+Bqb4Nelt+eh/wnyu5H4N+lWhr2l8HvlxV3xi6nnG1X6nvAi4aupZlugB4X1vbvh64MMl/DFvS8lXVXHs+CtzE/F1mp80R4MiC3xJvZD74p9nFwL1V9cRK7sSgXwXaHzH3Ag9X1eeHrme5kswkWdfarwbeCfx42KqWp6o+WVWzVbWF+V+t76yqvxu4rGVJckr7Iz9tqeNdwNRdrVZVjwOHk7y5DW0DpuqChRP4ECu8bANT/OXgSb4KvAM4NckR4NNVtXfYqpbtAuDDwANtfRvgU+0TxtNkI7CvXUXwCuCGqprqyxI7sQG4af58grXAV6rqm8OWtGz/AHy5LXk8Alw2cD3L1n7ovhP4+xXf17ReXilJGo1LN5LUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TO/R/k0rIVTkdbyQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(np.array(cks['koi_count'], dtype='float'), bins=np.arange(1,8), color='cornflowerblue', edgecolor='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adjust number of planets in each system to reflect removals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adjust number of planets in each CKS system to reflect removals \n",
    "cks_koi = cks['id_starname']\n",
    "unique_cks_systems = np.unique(cks_koi)\n",
    "\n",
    "for us in unique_cks_systems:\n",
    "    reduced_npl = np.sum(cks_koi == us)\n",
    "    cks['koi_count'][cks_koi == us] = reduced_npl\n",
    "\n",
    "    \n",
    "# adjust number of planets in each DR25 system to reflect removals \n",
    "dr25_koi = dr25['kepoi_name']\n",
    "dr25['id_starname'] = np.array([name[:6] for name in dr25_koi])\n",
    "unique_dr25_systems = np.unique(dr25['id_starname'])\n",
    "\n",
    "for us in unique_dr25_systems:\n",
    "    reduced_npl = np.sum(dr25_koi == us)\n",
    "    dr25['koi_count'][dr25_koi == us] = reduced_npl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure out which objects overlap between the two databases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here are the systems which are only found in one catalogue\n",
    "\n",
    "cks_koi  = cks['id_koicand']\n",
    "dr25_koi = dr25['kepoi_name']\n",
    "\n",
    "\n",
    "inboth = dr25_koi[np.isin(dr25_koi, cks_koi)]\n",
    "\n",
    "only_in_cks  = cks_koi[~np.isin(cks_koi, dr25_koi)]\n",
    "only_in_dr25 = dr25_koi[~np.isin(dr25_koi, cks_koi)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1373 objects in both catalogues\n",
      "16 objects only in CKS\n",
      "2968 objects only in dr25\n"
     ]
    }
   ],
   "source": [
    "print(len(inboth), 'objects in both catalogues')\n",
    "print(len(only_in_cks), 'objects only in CKS')\n",
    "print(len(only_in_dr25), 'objects only in dr25')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# identify and extract the high purity objects\n",
    "good  = np.isin(cks_koi,inboth)                    # boolean array, True for good systems\n",
    "good += cks['koi_disposition'] == 'CONFIRMED'\n",
    "\n",
    "my_catalogue = {}\n",
    "for k in cks.keys():\n",
    "    my_catalogue[k] = cks[k][good]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "864 objects in 335 unique systems\n",
      "\n",
      "412 objects in 2-planet systems\n",
      "243 objects in 3-planet systems\n",
      "136 objects in 4-planet systems\n",
      "60 objects in 5-planet systems\n",
      "6 objects in 6-planet systems\n",
      "7 objects in 7-planet systems\n"
     ]
    }
   ],
   "source": [
    "# readjust number of planets in each system to reflect removals \n",
    "hq_objects_starname = my_catalogue['id_starname']\n",
    "unique_systems      = np.unique(my_catalogue['id_starname'])\n",
    "\n",
    "for us in unique_systems:\n",
    "    reduced_npl = np.sum(hq_objects_starname == us)\n",
    "    my_catalogue['koi_count'][hq_objects_starname == us] = reduced_npl\n",
    "\n",
    "\n",
    "# remove systems reduced to a single planet\n",
    "npl = np.array(my_catalogue['koi_count'], dtype='float')\n",
    "\n",
    "for k in my_catalogue.keys():\n",
    "    my_catalogue[k] = my_catalogue[k][npl > 1]\n",
    "\n",
    "# grab a reference key\n",
    "k0 = list(my_catalogue.keys())[1]\n",
    "\n",
    "\n",
    "print('%d objects in %d unique systems' %(len(my_catalogue[k0]), len(np.unique(my_catalogue[k0]))))\n",
    "print('')\n",
    "print('%d objects in 2-planet systems' %np.sum(npl==2))\n",
    "print('%d objects in 3-planet systems' %np.sum(npl==3))\n",
    "print('%d objects in 4-planet systems' %np.sum(npl==4))\n",
    "print('%d objects in 5-planet systems' %np.sum(npl==5))\n",
    "print('%d objects in 6-planet systems' %np.sum(npl==6))\n",
    "print('%d objects in 7-planet systems' %np.sum(npl==7))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write out catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/clean_kepler_catalogue.csv'\n",
    "\n",
    "    with open(filepath, \"w\") as outfile:\n",
    "        writer = csv.writer(outfile)\n",
    "        writer.writerow(my_catalogue.keys())\n",
    "        writer.writerows(zip(*my_catalogue.values()))"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}