B_CRTS_SDSS_matching.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from astropy.table import join\n",
    "from astropy.table import Table\n",
    "from astropy.table import  hstack\n",
    "from astropy.table import Column\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from CRTS_paper_modules import update_progress\n",
    "import CRTS_paper_modules as mod\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "from os import listdir\n",
    "import pandas as pd \n",
    "import imp\n",
    "imp.reload(mod)\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CRTS stars "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aim : want to figure out WHY  all SDSS and CRTS stars can be matched to matching radius of 0.0 ! How is this possible? Are CRTS RA even different ? \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# CRTS stars radec\n",
    "address = '../raw_LC_CRTS/stars/radec.00'\n",
    "crts_stars_radec = Table.read(address, format='ascii', names=['CRTS_name','ra','dec'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# SDSS stars\n",
    "address = '../catalogs_SDSS/stripe82calibStars_v2.6.dat' \n",
    "colnames = ['calib_fla', 'ra', 'dec', 'raRMS', 'decRMS', 'nEpochs', 'AR_val', \n",
    "                'u_Nobs', 'u_mMed', 'u_mMean', 'u_mErr', 'u_rms_scatt', 'u_chi2',\n",
    "                'g_Nobs', 'g_mMed', 'g_mMean', 'g_mErr', 'g_rms_scatt', 'g_chi2',\n",
    "                'r_Nobs', 'r_mMed', 'r_mMean', 'r_mErr', 'r_rms_scatt', 'r_chi2',\n",
    "                'i_Nobs', 'i_mMed', 'i_mMean', 'i_mErr', 'i_rms_scatt', 'i_chi2',\n",
    "                'z_Nobs', 'z_mMed', 'z_mMean', 'z_mErr', 'z_rms_scatt', 'z_chi2']\n",
    "    \n",
    "sdss_stars = Table.read(address, format='ascii', names = colnames)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1006849"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(sdss_stars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " We combine aggregate info of raw LC and proc LC together with SDSS info, using the ra, dec which are for SDSS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# read in CRTS aggregate info  for stars \n",
    "address = '../data_products/CRTS_stars_raw_and_proc_combined_agg.dat'\n",
    "crts_stars = Table.read(address, format='ascii')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Add the SDSS ra, dec information.... \n",
    "crts_table = join(left=crts_stars, right = crts_stars_radec, keys='CRTS_name')\n",
    "\n",
    "# It's length is the same as the crts_stars :  48250 / 100000  stars, that fulfill the criteria \n",
    "# of being longer than 10 day-averaged epochs "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "48250"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check that indeed all ra are in the SDSS (because they are the same )\n",
    "\n",
    "np.sum(np.in1d(crts_table['ra'], sdss_stars['ra']))\n",
    "\n",
    "# Indeed they are! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "47748"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(np.unique(crts_table['dec']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Just join the SDSS info to CRTS table \n",
    "#crts_sdss_table = join(left=crts_table, right=sdss_stars )\n",
    "\n",
    "# Interestingly , for 48250 stars, there are only 47137 unique ra, and 47748 unique dec, i.e. some stars have either \n",
    "# the same ra, or the same dec. \n",
    "\n",
    "# So it seems that there are 820150 unique dec  , and 656841  unique ra  in sdss_stars.\n",
    "# even if we join on both 'ra' and 'dec', it seems that still we are getting more than one match from SDSS for\n",
    "# few stars . Just try using astropy catalog match, and see if it works perhaps better.... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Astropy  SkyCoord.match_to_catalog_sky \n",
    "crts_coord = SkyCoord(ra=crts_table['ra']*u.degree, \n",
    "                      dec=crts_table['dec']*u.degree)\n",
    "sdss_coord = SkyCoord(ra=sdss_stars['ra']*u.degree, \n",
    "                      dec=sdss_stars['dec']*u.degree)\n",
    "idx, sep2d, dist3d = crts_coord.match_to_catalog_sky(sdss_coord)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sdss_crts_table = hstack([crts_table, sdss_stars[idx]], table_names = ['crts','sdss'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Save the combined CRTS - SDSS table for stars... \n",
    "address = '../data_products/CRTS_SDSS_combined_stars_catalog.dat'\n",
    "sdss_crts_table.write(address, format='ascii', overwrite='True')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CRTS quasars  matching to SDSS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here I use CRTS aggregates on raw and processed lightcurves, and obtain ra,dec information from the filename.  I convert the hms to deg in two  ways:  using a handwritten routine, and using AstroPy SkyCoord.  Both yield exactly the same results : 15 QSO without a counterpart within 50 arcsec... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Use the aggregate info on raw and proc CRTS lightcurves, obtain the CRTS ra, dec from filename,\n",
    "# Read in the raw and processed CRTS LC aggregates \n",
    "\n",
    "address = '../data_products/CRTS_qso_raw_and_proc_combined_agg.dat'\n",
    "crts_qso = Table.read(address, format='ascii')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of those CRTS QSO without an SDSS match within 15 arcsec : 15\n",
      "The 15 outliers.... \n",
      "raw_timespan_obs  raw_mean_err  raw_lc_length ...   Au  match_angle_arcsec\n",
      "---------------- -------------- ------------- ... ----- ------------------\n",
      "      2996.86109 0.209255319149         282.0 ... 0.137      351.674559391\n",
      "      2996.86449        0.24195         200.0 ... 0.207      186.993884728\n",
      "      2933.01351 0.145868852459         305.0 ... 0.153      271.525022433\n",
      "      2917.05774 0.253652173913         115.0 ... 0.142      87.7045721106\n",
      "      2938.08685 0.114838709677         279.0 ... 0.143      450.192877296\n",
      "      2914.09165 0.343448275862          58.0 ... 0.371      247.422614573\n",
      "      2849.23653 0.364444444444          45.0 ... 0.484      144.748025985\n",
      "      2938.09008 0.337542372881         118.0 ... 0.585      551.280989912\n",
      "      2920.01124 0.246804979253         241.0 ... 0.637      213.254624003\n",
      "      1924.83354 0.182708333333          48.0 ... 0.468      142.664246028\n",
      "      3055.87807 0.200588235294         340.0 ... 0.348      390.143285679\n",
      "      2878.11108 0.383571428571          14.0 ... 0.313      470.316607554\n",
      "      3025.83213 0.177466666667         375.0 ... 0.187      242.141357687\n",
      "      3025.82255 0.225971223022         139.0 ... 0.225      410.624632156\n",
      "      3005.89569 0.252040816327         196.0 ...  0.16      450.966870305\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "# extract ra, dec from  CRTS_name   : \n",
    "ra_string = [name[:-9] for name in  crts_qso['CRTS_name']]\n",
    "dec_string  = [name[9:] for name in crts_qso['CRTS_name']]\n",
    "\n",
    "\n",
    "# Split CRTS  ra, dec from hms to h m s \n",
    "ra_hms_split, dec_hms_split = mod.get_ra_dec_CRTS(ra_string, dec_string)\n",
    "\n",
    "# Convert CRTS  ra, dec from hms to deg  \n",
    "ra_deg_CRTS, dec_deg_CRTS = mod.convert_to_deg(ra_hms_split, dec_hms_split)\n",
    "\n",
    "\n",
    "# add that as two columns :\n",
    "ra = Column(name='ra_CRTS', data=ra_deg_CRTS)\n",
    "crts_qso.add_column(ra)\n",
    "\n",
    "dec = Column(name='dec_CRTS', data=dec_deg_CRTS)\n",
    "crts_qso.add_column(dec)\n",
    "\n",
    "# Load SDSS Quasars \n",
    "\n",
    "File = '../catalogs_SDSS/DB_QSO_S82.dat'\n",
    "colnames = ['dbID', 'ra', 'dec', 'SDR5ID', 'M_i', 'M_i_corr', 'redshift',\n",
    "                   'mass_BH', 'Lbol', 'u', 'g', 'r', 'i', 'z', 'Au']\n",
    "sdss_qso = Table.read(File, format='ascii', names = colnames)\n",
    "\n",
    "use_astropy_angles = False \n",
    "if use_astropy_angles == True : \n",
    "    CRTS_coords = np.zeros_like(ra_hms_split, dtype='<U31')\n",
    "    for  i in range(len(comb)):\n",
    "        CRTS_coords[i] = ra_hms_split[i]+' '+dec_hms_split[i]\n",
    "    CRTS_SkyCoords = SkyCoord(CRTS_coords, unit=(u.hourangle, u.deg))\n",
    "\n",
    "    cat1 = CRTS_SkyCoords\n",
    "else : \n",
    "    cat1_ra  = crts_qso['ra_CRTS']\n",
    "    cat1_dec = crts_qso['dec_CRTS']\n",
    "    cat1 = SkyCoord(ra=cat1_ra*u.degree, dec=cat1_dec*u.degree)\n",
    "    \n",
    "# Merge positionally the two catalogs with AStropy \n",
    "cat2_ra  = sdss_qso['ra']\n",
    "cat2_dec = sdss_qso['dec']\n",
    "\n",
    "cat2 = SkyCoord(ra=cat2_ra*u.degree, dec=cat2_dec*u.degree)\n",
    "idx, sep2d, dist3d = cat1.match_to_catalog_sky(cat2) \n",
    "\n",
    "match_angle_deg = sep2d.value\n",
    "match_angle_arcsec = match_angle_deg * 3600\n",
    "\n",
    "# The fifteen outliers...\n",
    "print('Number of those CRTS QSO without an SDSS match within 15 arcsec : %d' %np.sum(match_angle_arcsec> 15))\n",
    "\n",
    "# take only those rows that are a perfect positional match to our catalog ... \n",
    "# even though we didn't have to do this, the fact that some stars share ra or dec made it hard \n",
    "# to use join() \n",
    "sdss_qso_table = sdss_qso[idx]\n",
    "\n",
    "match_angle = Column(name='match_angle_arcsec',data= match_angle_arcsec)\n",
    "sdss_qso_table.add_column(match_angle)\n",
    "\n",
    "# add the sdss info to the crts info : stack two tables horizontally \n",
    "sdss_crts_table = hstack([crts_qso, sdss_qso_table])\n",
    "\n",
    "# Display the fifteen outliers...\n",
    "m = sdss_crts_table['match_angle_arcsec'] > 1\n",
    "print('The 15 outliers.... ')\n",
    "print(sdss_crts_table[m])\n",
    "\n",
    "# Save the combined CRTS - SDSS table for quasars ... \n",
    "address = '../data_products/CRTS_SDSS_combined_QSO_catalog.dat'\n",
    "sdss_crts_table.write(address, format='ascii', overwrite='True')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['raw_timespan_obs', 'raw_mean_err', 'raw_lc_length', 'raw_mean_mag',\n",
       "       'CRTS_name', 'proc_mean_N_day', 'proc_lc_length', 'proc_mean_mag',\n",
       "       'proc_mjd_span', 'proc_mean_err', 'proc_median_dt', 'ra_CRTS',\n",
       "       'dec_CRTS', 'dbID', 'ra', 'dec', 'SDR5ID', 'M_i', 'M_i_corr',\n",
       "       'redshift', 'mass_BH', 'Lbol', 'u', 'g', 'r', 'i', 'z', 'Au',\n",
       "       'match_angle_arcsec'], \n",
       "      dtype='<U18')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.ravel(sdss_crts_table.colnames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Make an ascii file with list of ra, dec in degrees, separated by white space, commas (,), or pipes (|). \n",
    "# RA and Dec must be in units decimal degrees.\n",
    "\n",
    "# first, feed the SDSS DR12 with the SDSS coords of mismatched QSO : \n",
    "radec_SDSS_mismatch = sdss_crts_table1[sdss_crts_table['match_angle_arcsec']>2][['ra','dec']]\n",
    "address = '../data_products/CRTS_SDSS_mismatch_QSO_SDSS_radec.dat'\n",
    "radec_SDSS_mismatch.write(address, format='ascii', overwrite='True')\n",
    "\n",
    "# Feed the SDSS with the CRTS coords of mismatched objects .... \n",
    "\n",
    "radec_CRTS_mismatch = sdss_crts_table1[sdss_crts_table1['match_angle_arcsec']>2][['ra_CRTS','dec_CRTS']]\n",
    "address = '../data_products/CRTS_SDSS_mismatch_QSO_CRTS_radec.dat'\n",
    "radec_CRTS_mismatch.write(address, format='ascii', overwrite='True')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    " Could do https://dr12.sdss.org/bulkFields/raDecs ,but it does not give much useful info...\n",
    "\n",
    " Instead, by hand try http://skyserver.sdss.org/dr7/en/tools/explore/obj.asp\n",
    "\n",
    " NOTES on individual objects : \n",
    "\n",
    " First, I use  SDSS coords from DB9  to query by hand these 15 objects...\n",
    "\n",
    "1) SDSS STAR ra=14.16039248, dec=0.61147335,   ObjId = 588015509812412634 : \n",
    "    \n",
    "2) SDSS GALAXY ra=16.57797276, dec=-0.36779433,   ObjId = 588015508739719362\n",
    "\n",
    "3) SDSS STAR    ra=21.79044991, dec=-0.82312779,   ObjId = 588015508205142178\n",
    "\n",
    "4) SDSS STAR    ra=23.57580185, dec=0.26019066,   ObjId = 587731513144770582\n",
    "\n",
    "5) SDSS STAR    ra=30.27577764, dec=-0.9485065,   ObjId = 587731511537041576\n",
    "\n",
    "6) SDSS STAR    ra=43.79397299, dec=-0.12030369,   ObjId = 587731512616747105\n",
    "\n",
    "7) SDSS STAR    ra=44.02973877, dec=1.07722393,   ObjId = 587731514227425430\n",
    "\n",
    "8) SDSS STAR    ra=51.85132224, dec=-0.85782689,   ObjId = 587731511546478754\n",
    "\n",
    "9) SDSS STAR    ra=55.73044693, dec=0.8231733,   ObjId = 587731513695666438\n",
    "\n",
    "10) SDSS STAR    ra=311.17075564, dec=0.94198223,   ObjId = 587730848497206143\n",
    "--> crowded field \n",
    "\n",
    "11) SDSS STAR    ra=335.40398149, dec=0.19558045,   ObjId = 587734304878690655\n",
    "\n",
    "12) SDSS STAR    ra=343.92565411, dec=-1.12342469,   ObjId = 587734303271813146\n",
    "\n",
    "13) SDSS STAR    ra=348.04136769, dec=-0.94963622,   ObjId = 587731185124966472\n",
    "\n",
    "14) SDSS STAR    ra=350.34892733, dec=-0.02322992,   ObjId = 587731186199691464\n",
    "\n",
    "15) SDSS STAR    ra=355.04130819, dec=0.93886394,   ObjId = 588015510340960369\n",
    "\n",
    "--> some looked like QSO, but then got spectra, but not all have spectra! \n",
    "--> this is surprising... well, I was looking at DR7, is there anything more possibly in DR9 ? \n",
    "\n",
    "--> DR12 catalog is not easy to browse... \n",
    "\n",
    "--> Here in DR13 I can paste ra, dec and get poststamps that can be further inspected ... \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exec. summary on the CRTS-SDSS mismatched QSO ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These 15 objects according to DB9 are QSO, and indeed they are. The CRTS objects that are mismatched sometimes are quasars, sometimes they are not.... I wonder where Brani took the QSO from. He said they were spectroscopically confirmed. All mismatches are 'catastrophic', i.e. more than 50 arcsec separation. \n",
    "\n",
    "Possible problems: \n",
    "   - are we using different coordinate systems between CRTS and SDSS ? The default frame in SkyCoord is ICRS   (\"The International Celestial Reference System (ICRS) is the current standard celestial reference system adopted by the International Astronomical Union (IAU). Its origin is at the barycenter of the Solar System, with axes that are intended to be \"fixed\" with respect to space.\")\n",
    "   - problems with hms --> deg conversion : unlikely. I checked the matching in two ways, and results still persist, with exactly the same 15 quasars being mismatched (without an SDSS counterpart within 2 arcsec).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}