{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HDF-N master catalogue\n",
    "## Preparation of HST CANDELS-3D data\n",
    "\n",
    "The catalogue comes from `dmu0_CANDELS-3D-HST`. We no longer use this catalogue since it is superseeded by CANDELS-HDFN\n",
    "\n",
    "In the catalogue, we keep:\n",
    "\n",
    "- The identifier (it's unique in the catalogue);\n",
    "- The position;\n",
    "- The stellarity;\n",
    "- The kron magnitude, there doesn't appear to be aperture magnitudes. This may mean the survey is unusable.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "from herschelhelp_internal import git_version\n",
    "print(\"This notebook was run with herschelhelp_internal version: \\n{}\".format(git_version()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
    "from collections import OrderedDict\n",
    "import os\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.table import Column, Table\n",
    "import numpy as np\n",
    "\n",
    "from herschelhelp_internal.flagging import  gaia_flag_column\n",
    "from herschelhelp_internal.masterlist import nb_astcor_diag_plot, remove_duplicates\n",
    "from herschelhelp_internal.utils import astrometric_correction, flux_to_mag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "OUT_DIR =  os.environ.get('TMP_DIR', \"./data_tmp\")\n",
    "try:\n",
    "    os.makedirs(OUT_DIR)\n",
    "except FileExistsError:\n",
    "    pass\n",
    "\n",
    "RA_COL = \"candels_ra\"\n",
    "DEC_COL = \"candels_dec\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I - Column selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "imported_columns = OrderedDict({\n",
    "        'ID': \"candels_id\",\n",
    "        'RAJ2000': \"candels_ra\",\n",
    "        'DEJ2000': \"candels_dec\",\n",
    "        'S/G':  \"candels_stellarity\",\n",
    "        'F140Wap': \"f_ap_candels_f140w\",\n",
    "        'e_F140Wap': \"ferr_ap_candels_f140w\",\n",
    "        'F140W': \"f_candels_f140w\",\n",
    "        'e_F140W': \"ferr_candels_f140w\",\n",
    "        'F160Wap': \"f_ap_candels_f160w\",\n",
    "        'e_F160Wap': \"ferr_ap_candels_f160w\",\n",
    "        'F160W': \"f_candels_f160w\",\n",
    "        'e_F160W': \"ferr_candels_f160w\",\n",
    "        'F606W': \"f_candels_f606w\",\n",
    "        'e_F606W': \"ferr_candels_f606w\",\n",
    "        'F814W': \"f_candels_f814w\",\n",
    "        'e_F814W': \"ferr_candels_f814w\",\n",
    "        'F125W': \"f_candels_f125w\",\n",
    "        'e_F125W': \"ferr_candels_f125w\"\n",
    "\n",
    "    })\n",
    "\n",
    "\n",
    "catalogue = Table.read(\"../../dmu0/dmu0_CANDELS-3D-HST/data/CANDELS-3D-HST_HDF-N.fits\")[list(imported_columns)]\n",
    "for column in imported_columns:\n",
    "    catalogue[column].name = imported_columns[column]\n",
    "\n",
    "epoch = 2012 #Year of publication\n",
    "\n",
    "#Clean table metadata\n",
    "catalogue.meta = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Adding flux and band-flag columns\n",
    "for col in catalogue.colnames:\n",
    "    if col.startswith('f_'):\n",
    "        \n",
    "        errcol = \"ferr{}\".format(col[1:])\n",
    "        \n",
    "        #Calculate mags, errors including the fact that fluxes are in units of 0.3631 uJy\n",
    "        mag, error = flux_to_mag(np.array(catalogue[col]) * 0.3631e-6, np.array(catalogue[errcol] * 0.3631e-6))\n",
    "        \n",
    "        \n",
    "        # magnitudes are added\n",
    "        catalogue.add_column(Column(mag, name=\"m{}\".format(col[1:])))\n",
    "        catalogue.add_column(Column(error, name=\"m{}\".format(errcol[1:])))\n",
    "        \n",
    "        \n",
    "        #Correct flux units to uJy\n",
    "        catalogue[col] = catalogue[col] * 0.3631\n",
    "        catalogue[col].unit = u.microjansky\n",
    "        catalogue[errcol] = catalogue[errcol] * 0.3631\n",
    "        catalogue[errcol].unit = u.microjansky        \n",
    "\n",
    "        if ('125' in col) or ('814' in col) or ('606' in col) :\n",
    "            # We add nan filled aperture photometry for consistency\n",
    "            catalogue.add_column(Column(np.full(len(catalogue), np.nan), name=\"m_ap{}\".format(col[1:])))\n",
    "            catalogue.add_column(Column(np.full(len(catalogue), np.nan), name=\"merr_ap{}\".format(col[1:])))\n",
    "            catalogue.add_column(Column(np.full(len(catalogue), np.nan), name=\"f_ap{}\".format(col[1:])))\n",
    "            catalogue.add_column(Column(np.full(len(catalogue), np.nan), name=\"ferr_ap{}\".format(col[1:])))\n",
    "\n",
    "        \n",
    "\n",
    "        \n",
    "        # Band-flag column\n",
    "        if \"ap\" not in col:\n",
    "            catalogue.add_column(Column(np.zeros(len(catalogue), dtype=bool), name=\"flag{}\".format(col[1:])))\n",
    "        \n",
    "#TODO: Set to True the flag columns for fluxes that should not be used for SED fitting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "catalogue[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II - Removal of duplicated sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We remove duplicated objects from the input catalogues."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "SORT_COLS = ['merr_candels_f140w', \n",
    "             'merr_candels_f160w', \n",
    "             'merr_candels_f606w', \n",
    "             'merr_candels_f814w',\n",
    "             'merr_candels_f125w']\n",
    "FLAG_NAME = 'candels_flag_cleaned'\n",
    "\n",
    "nb_orig_sources = len(catalogue)\n",
    "\n",
    "catalogue = remove_duplicates(catalogue, RA_COL, DEC_COL, sort_col=SORT_COLS,flag_name=FLAG_NAME)\n",
    "\n",
    "nb_sources = len(catalogue)\n",
    "\n",
    "print(\"The initial catalogue had {} sources.\".format(nb_orig_sources))\n",
    "print(\"The cleaned catalogue has {} sources ({} removed).\".format(nb_sources, nb_orig_sources - nb_sources))\n",
    "print(\"The cleaned catalogue has {} sources flagged as having been cleaned\".format(np.sum(catalogue[FLAG_NAME])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III - Astrometry correction\n",
    "\n",
    "We match the astrometry to the Gaia one. We limit the Gaia catalogue to sources with a g band flux between the 30th and the 70th percentile. Some quick tests show that this give the lower dispersion in the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "gaia = Table.read(\"../../dmu0/dmu0_GAIA/data/GAIA_HDF-N.fits\")\n",
    "gaia_coords = SkyCoord(gaia['ra'], gaia['dec'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL], \n",
    "                    gaia_coords.ra, gaia_coords.dec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "delta_ra, delta_dec =  astrometric_correction(\n",
    "    SkyCoord(catalogue[RA_COL], catalogue[DEC_COL]),\n",
    "    gaia_coords\n",
    ")\n",
    "\n",
    "print(\"RA correction: {}\".format(delta_ra))\n",
    "print(\"Dec correction: {}\".format(delta_dec))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "catalogue[RA_COL] +=  delta_ra.to(u.deg)\n",
    "catalogue[DEC_COL] += delta_dec.to(u.deg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL], \n",
    "                    gaia_coords.ra, gaia_coords.dec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IV - Flagging Gaia objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "catalogue.add_column(\n",
    "    gaia_flag_column(SkyCoord(catalogue[RA_COL], catalogue[DEC_COL]), epoch, gaia)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GAIA_FLAG_NAME = \"candels_flag_gaia\"\n",
    "\n",
    "catalogue['flag_gaia'].name = GAIA_FLAG_NAME\n",
    "print(\"{} sources flagged.\".format(np.sum(catalogue[GAIA_FLAG_NAME] > 0)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## V - Flagging objects near bright stars"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# VI - Saving to disk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "catalogue.write(\"{}/CANDELS.fits\".format(OUT_DIR), overwrite=True)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}