{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HDF-N master catalogue\n",
    "## Preparation of GOODS-ACS data\n",
    "\n",
    "The catalogue comes from GOODS-ACS but has been superseeded by CANDELS-GOODS-N so is not used here\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, join\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 = \"acs_ra\"\n",
    "DEC_COL = \"acs_dec\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I - Column selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "b_catalogue = Table.read(\"../../dmu0/dmu0_GOODS-ACS/data/h_goods_nb_r2.0z_cat.fits\") #[list(imported_columns)]\n",
    "i_catalogue = Table.read(\"../../dmu0/dmu0_GOODS-ACS/data/h_goods_ni_r2.0z_cat.fits\") #[list(imported_columns)]\n",
    "v_catalogue = Table.read(\"../../dmu0/dmu0_GOODS-ACS/data/h_goods_nv_r2.0z_cat.fits\") #[list(imported_columns)]\n",
    "z_catalogue = Table.read(\"../../dmu0/dmu0_GOODS-ACS/data/h_goods_nz_r2.0z_cat.fits\") #[list(imported_columns)]\n",
    "\n",
    "for column in b_catalogue.colnames:\n",
    "    if not b_catalogue[column].name == 'ID_IAU':\n",
    "        b_catalogue[column].name = 'b_' + b_catalogue[column].name \n",
    "\n",
    "for column in i_catalogue.colnames:\n",
    "    if not i_catalogue[column].name == 'ID_IAU':\n",
    "        i_catalogue[column].name = 'i_' + i_catalogue[column].name \n",
    "\n",
    "for column in v_catalogue.colnames:\n",
    "    if not v_catalogue[column].name == 'ID_IAU':\n",
    "        v_catalogue[column].name = 'v_' + v_catalogue[column].name \n",
    "    \n",
    "for column in z_catalogue.colnames:\n",
    "    if not z_catalogue[column].name == 'ID_IAU':\n",
    "        z_catalogue[column].name = 'z_' + z_catalogue[column].name \n",
    "    \n",
    "epoch = 2012 #Year of publication\n",
    "\n",
    "catalogue = b_catalogue\n",
    "\n",
    "catalogue = join(b_catalogue,i_catalogue, keys='ID_IAU')\n",
    "catalogue = join(catalogue,v_catalogue, keys='ID_IAU')\n",
    "catalogue = join(catalogue,z_catalogue, keys='ID_IAU')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "catalogue[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "imported_columns = OrderedDict({\n",
    "        'ID_IAU': \"acs_id\",\n",
    "        'b_ALPHA_J2000': \"acs_ra\",\n",
    "        'b_DELTA_J2000': \"acs_dec\",\n",
    "        'b_MAG_AUTO':  \"m_acs_b\",\n",
    "        'b_MAGERR_AUTO': \"merr_acs_b\",\n",
    "        'b_FLUX_AUTO': \"f_acs_b\",\n",
    "        'b_FLUXERR_AUTO': \"ferr_acs_b\",\n",
    "        'e_F140W': \"ferr_candels_f140w\",\n",
    "\n",
    "\n",
    "    })\n",
    "catalogue = catalogue[list(imported_columns)]\n",
    "\n",
    "for column in imported_columns:\n",
    "    catalogue[column].name = imported_columns[column]"
   ]
  },
  {
   "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_XMM-LSS.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
}