{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# ELAIS-N1 - Merging HELP data products\n", "\n", "This notebook merges the various HELP data products on ELAIS-N1.\n", "\n", "It is first used to create a catalogue that will be used for SED fitting by CIGALE by merging the optical master list, the photo-z and the XID+ far infrared fluxes. Then, this notebook is used to incorporate the CIGALE physical parameter estimations and generate the final HELP data product on the field." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Set this to true to produce only the catalogue for CIGALE and to false \n", "# to continue and merge the CIGALE results too.\n", "FIRST_RUN_FOR_CIGALE = False" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "\n", "from astropy.table import Column, MaskedColumn, Table, join, vstack\n", "\n", "from herschelhelp.filters import get_filter_meta_table\n", "\n", "from herschelhelp_internal.utils import add_column_meta\n", "\n", "filter_mean_lambda = {\n", " item['filter_id']: item['mean_wavelength'] for item in\n", " get_filter_meta_table()\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Reading the masterlist, XID+, and photo-z catalogues" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Master list\n", "\n", "SUFFIX = \"20171016\"\n", "ml = Table.read(\n", " \"../../dmu1/dmu1_ml_ELAIS-N1/data/master_catalogue_elais-n1_{}.fits\".format(SUFFIX))\n", "ml.meta = None" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# XID+ MIPS24\n", "# There are two catalogues on SWIRE and SERVS coverage that must be \n", "# merged (they don't overlap).\n", "\n", "mips24_servs = Table.read(\"../../dmu26/dmu26_XID+MIPS_ELAIS-N1/data/\"\n", " \"dmu26_XID+MIPS_ELAIS-N1_SERVS_cat_20170725.fits\")\n", "mips24_servs.meta = None\n", "mips24_swire = Table.read(\"../../dmu26/dmu26_XID+MIPS_ELAIS-N1/data/\"\n", " \"dmu26_XID+MIPS_ELAIS-N1_SWIRE_cat_20170725.fits\")\n", "mips24_swire.meta = None\n", "mips24_swire.remove_column(\"flag_mips24\") # duplicated column\n", "\n", "xid_mips24 = vstack([mips24_servs, mips24_swire])\n", "\n", "# Adding the error column\n", "xid_mips24.add_column(Column(\n", " data=np.max([xid_mips24['FErr_MIPS_24_u'] - xid_mips24['F_MIPS_24'],\n", " xid_mips24['F_MIPS_24'] - xid_mips24['FErr_MIPS_24_l']],\n", " axis=0),\n", " name=\"ferr_mips_24\"\n", "))\n", "xid_mips24['F_MIPS_24'].name = \"f_mips_24\"\n", "xid_mips24 = xid_mips24['help_id', 'f_mips_24', 'ferr_mips_24', 'flag_mips_24']" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# XID+ PACS\n", "# There are two catalogues on SWIRE and SERVS coverage that must be \n", "# merged (they don't overlap).\n", "\n", "pacs_swire = Table.read(\"../../dmu26/dmu26_XID+PACS_ELAIS-N1/data/\"\n", " \"dmu26_XID+PACS_ELAIS-N1_SWIRE_cat_20170808.fits\")\n", "pacs_swire.meta = None\n", "pacs_servs = Table.read(\"../../dmu26/dmu26_XID+PACS_ELAIS-N1/data/\"\n", " \"dmu26_XID+PACS_ELAIS-N1_SERVS_cat_20170808.fits\")\n", "pacs_servs.meta = None\n", "\n", "xid_pacs = vstack([pacs_servs, pacs_swire])\n", "\n", "# Convert from mJy to μJy\n", "for col in [\"F_PACS_100\", \"FErr_PACS_100_u\", \"FErr_PACS_100_l\",\n", " \"F_PACS_160\", \"FErr_PACS_160_u\", \"FErr_PACS_160_l\"]:\n", " xid_pacs[col] *= 1000\n", "\n", "xid_pacs.add_column(Column(\n", " data=np.max([xid_pacs['FErr_PACS_100_u'] - xid_pacs['F_PACS_100'],\n", " xid_pacs['F_PACS_100'] - xid_pacs['FErr_PACS_100_l']],\n", " axis=0),\n", " name=\"ferr_pacs_green\"\n", "))\n", "xid_pacs['F_PACS_100'].name = \"f_pacs_green\"\n", "xid_pacs['flag_PACS_100'].name = \"flag_pacs_green\"\n", "\n", "xid_pacs.add_column(Column(\n", " data=np.max([xid_pacs['FErr_PACS_160_u'] - xid_pacs['F_PACS_160'],\n", " xid_pacs['F_PACS_160'] - xid_pacs['FErr_PACS_160_l']],\n", " axis=0),\n", " name=\"ferr_pacs_red\"\n", "))\n", "xid_pacs['F_PACS_160'].name = \"f_pacs_red\"\n", "xid_pacs['flag_PACS_160'].name = \"flag_pacs_red\"\n", "\n", "xid_pacs = xid_pacs['help_id', 'f_pacs_green', 'ferr_pacs_green',\n", " 'flag_pacs_green', 'f_pacs_red', 'ferr_pacs_red',\n", " 'flag_pacs_red']\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# XID+ SPIRE\n", "# There are two catalogues on SWIRE and SERVS coverage that must be \n", "# merged (they don't overlap).\n", "\n", "spire_servs = Table.read(\"../../dmu26/dmu26_XID+SPIRE_ELAIS-N1/data/\"\n", " \"dmu26_XID+SPIRE_ELAIS-N1_SERVS_cat_20170725.fits\")\n", "spire_servs.meta = None\n", "spire_swire = Table.read(\"../../dmu26/dmu26_XID+SPIRE_ELAIS-N1/data/\"\n", " \"dmu26_XID+SPIRE_ELAIS-N1_SWIRE_cat_20170808.fits\")\n", "spire_swire.meta = None\n", "\n", "xid_spire = vstack([spire_servs, spire_swire])\n", "\n", "# Convert from mJy to μJy\n", "for col in [\"F_SPIRE_250\", \"FErr_SPIRE_250_u\", \"FErr_SPIRE_250_l\",\n", " \"F_SPIRE_350\", \"FErr_SPIRE_350_u\", \"FErr_SPIRE_350_l\",\n", " \"F_SPIRE_500\", \"FErr_SPIRE_500_u\", \"FErr_SPIRE_500_l\"]:\n", " xid_spire[col] *= 1000\n", "\n", "xid_spire.add_column(Column(\n", " data=np.max([xid_spire['FErr_SPIRE_250_u'] - xid_spire['F_SPIRE_250'],\n", " xid_spire['F_SPIRE_250'] - xid_spire['FErr_SPIRE_250_l']],\n", " axis=0),\n", " name=\"ferr_spire_250\"\n", "))\n", "xid_spire['F_SPIRE_250'].name = \"f_spire_250\"\n", "xid_spire.add_column(Column(\n", " data=np.max([xid_spire['FErr_SPIRE_350_u'] - xid_spire['F_SPIRE_350'],\n", " xid_spire['F_SPIRE_350'] - xid_spire['FErr_SPIRE_350_l']],\n", " axis=0),\n", " name=\"ferr_spire_350\"\n", "))\n", "xid_spire['F_SPIRE_350'].name = \"f_spire_350\"\n", "xid_spire.add_column(Column(\n", " data=np.max([xid_spire['FErr_SPIRE_500_u'] - xid_spire['F_SPIRE_500'],\n", " xid_spire['F_SPIRE_500'] - xid_spire['FErr_SPIRE_500_l']],\n", " axis=0),\n", " name=\"ferr_spire_500\"\n", "))\n", "xid_spire['F_SPIRE_500'].name = \"f_spire_500\"\n", "\n", "xid_spire = xid_spire['help_id',\n", " 'f_spire_250', 'ferr_spire_250', 'flag_spire_250',\n", " 'f_spire_350', 'ferr_spire_350', 'flag_spire_350',\n", " 'f_spire_500', 'ferr_spire_500', 'flag_spire_500']" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Photo-z\n", "\n", "photoz = Table.read(\"../../dmu24/dmu24_ELAIS-N1/data/\"\n", " \"master_catalogue_elais-n1_20170706_photoz_20170725_irac1_optimised.fits\")\n", "photoz.meta = None\n", "\n", "photoz = photoz['help_id', 'z1_median']\n", "photoz['z1_median'].name = 'redshift'\n", "\n", "photoz['redshift'][photoz['redshift'] < 0] = np.nan # -99 used for missing values" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/astropy/table/column.py:965: RuntimeWarning: invalid value encountered in less\n", " return getattr(self.data, op)(other)\n" ] } ], "source": [ "# Spec-z\n", "ml['zspec'][ml['zspec'] < 0] = np.nan # -99 used for missing values" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Flags\n", "flags = Table.read(\"../../dmu6/dmu6_v_ELAIS-N1/data/elais-n1_20171016_flags.fits\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Merging" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "merged_table = join(ml, xid_mips24, join_type='left')\n", "\n", "# Fill values\n", "for col in xid_mips24.colnames:\n", " if col.startswith(\"f_\") or col.startswith(\"ferr_\"):\n", " merged_table[col].fill_value = np.nan\n", " elif col.startswith(\"flag_\"):\n", " merged_table[col].fill_value = False\n", "merged_table = merged_table.filled()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "merged_table = join(merged_table, xid_pacs, join_type='left')\n", " \n", "# Fill values\n", "for col in xid_pacs.colnames:\n", " if col.startswith(\"f_\") or col.startswith(\"ferr_\"):\n", " merged_table[col].fill_value = np.nan\n", " elif col.startswith(\"flag_\"):\n", " merged_table[col].fill_value = False\n", "merged_table = merged_table.filled()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "merged_table = join(merged_table, xid_spire, join_type='left')\n", " \n", "# Fill values\n", "for col in xid_spire.colnames:\n", " if col.startswith(\"f_\") or col.startswith(\"ferr_\"):\n", " merged_table[col].fill_value = np.nan\n", " elif col.startswith(\"flag_\"):\n", " merged_table[col].fill_value = False\n", "merged_table = merged_table.filled()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true }, "outputs": [], "source": [ "merged_table = join(merged_table, photoz, join_type='left')\n", "\n", "# Fill values\n", "merged_table['redshift'].fill_value = np.nan\n", "merged_table = merged_table.filled()" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "for col in flags.colnames:\n", " if 'flag' in col:\n", " try:\n", " merged_table.remove_column(col)\n", " except KeyError:\n", " print(\"Column: {} not in masterlist.\".format(col))\n", " \n", "merged_table = join(merged_table, flags, join_type='left')\n", "# Fill values\n", "for col in merged_table.colnames:\n", " if 'flag' in col:\n", " merged_table[col].fill_value = False\n", "merged_table = merged_table.filled()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Saving the catalogue for CIGALE (first run)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ "if FIRST_RUN_FOR_CIGALE:\n", " \n", " # Sorting the columns\n", " bands_tot = [col[2:] for col in merged_table.colnames\n", " if col.startswith('f_') and not col.startswith('f_ap')]\n", " bands_ap = [col[5:] for col in merged_table.colnames\n", " if col.startswith('f_ap_') ]\n", " bands = list(set(bands_tot) | set(bands_ap))\n", " bands.sort(key=lambda x: filter_mean_lambda[x])\n", " \n", " columns = ['help_id', 'field', 'ra', 'dec', 'hp_idx', 'ebv', 'redshift', \n", " 'zspec']\n", " for band in bands:\n", " for col_tpl in ['f_{}', 'ferr_{}', 'f_ap_{}', 'ferr_ap_{}',\n", " 'm_{}', 'merr_{}', 'm_ap_{}', 'merr_ap_{}',\n", " 'flag_{}']:\n", " colname = col_tpl.format(band)\n", " if colname in merged_table.colnames:\n", " columns.append(colname)\n", " columns += ['stellarity', 'stellarity_origin', 'flag_cleaned',\n", " 'flag_merged', 'flag_gaia', 'flag_optnir_obs',\n", " 'flag_optnir_det', 'zspec_qual', 'zspec_association_flag']\n", " \n", "\n", " # Check that we did not forget any column\n", " assert set(columns) == set(merged_table.colnames)\n", " \n", " merged_table = add_column_meta(merged_table, '../columns.yml')\n", " merged_table[columns].write(\"data/ELAIS-N1_{}_cigale.fits\".format(SUFFIX))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Merging CIGALE outputs\n", "\n", "We merge the CIGALE outputs to the main catalogue. The CIGALE products provides several χ² with associated thresholds. For simplicity, we convert these two values to flags." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Cigale outputs\n", "cigale = Table.read(\"../../dmu28/dmu28_ELAIS-N1/data/zphot/HELP_final_results.fits\")\n", "cigale['id'].name = \"help_id\"\n", "\n", "# We convert the various Chi2 and threshold to flags\n", "flag_cigale_opt = cigale[\"UVoptIR_OPTchi2\"] <= cigale[\"UVoptIR_OPTchi2_threshold\"]\n", "flag_cigale_ir = cigale[\"UVoptIR_IRchi2\"] <= cigale[\"UVoptIR_IRchi2_threshold\"]\n", "flag_cigale = (\n", " (cigale[\"UVoptIR_best.reduced_chi_square\"] \n", " <= cigale[\"UVoptIR_best.reduced_chi_square_threshold\"]) &\n", " flag_cigale_opt & flag_cigale_ir)\n", "flag_cigale_ironly = cigale[\"IRonly_IRchi2\"] <= cigale[\"IRonly_IRchi2_threshold\"]\n", "\n", "cigale.add_columns([\n", " MaskedColumn(flag_cigale, \"flag_cigale\", \n", " dtype=int, fill_value=-1),\n", " MaskedColumn(flag_cigale_opt, \"flag_cigale_opt\", \n", " dtype=int, fill_value=-1),\n", " MaskedColumn(flag_cigale_ir, \"flag_cigale_ir\", \n", " dtype=int, fill_value=-1),\n", " MaskedColumn(flag_cigale_ironly, \"flag_cigale_ironly\", \n", " dtype=int, fill_value=-1)\n", "])\n", "\n", "cigale['UVoptIR_bayes.stellar.m_star'].name = \"cigale_mstar\"\n", "cigale['UVoptIR_bayes.stellar.m_star_err'].name = \"cigale_mstar_err\"\n", "cigale['UVoptIR_bayes.sfh.sfr10Myrs'].name = \"cigale_sfr\"\n", "cigale['UVoptIR_bayes.sfh.sfr10Myrs_err'].name = \"cigale_sfr_err\"\n", "cigale['UVoptIR_bayes.dust.luminosity'].name = \"cigale_dustlumin\"\n", "cigale['UVoptIR_bayes.dust.luminosity_err'].name = \"cigale_dustlumin_err\"\n", "cigale['IR_bayes.dust.luminosity'].name = \"cigale_dustlumin_ironly\"\n", "cigale['IR_bayes.dust.luminosity_err'].name = \"cigale_dustlumin_ironly_err\"\n", "\n", "cigale = cigale['help_id', 'cigale_mstar', 'cigale_mstar_err', 'cigale_sfr',\n", " 'cigale_sfr_err', 'cigale_dustlumin', 'cigale_dustlumin_err', \n", " 'cigale_dustlumin_ironly', 'cigale_dustlumin_ironly_err', \n", " 'flag_cigale', 'flag_cigale_opt', 'flag_cigale_ir', \n", " 'flag_cigale_ironly']" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true }, "outputs": [], "source": [ "merged_table = join(merged_table, cigale, join_type='left')\n", "\n", "# Fill values\n", "for col in cigale.colnames:\n", " if col.startswith(\"cigale_\"):\n", " merged_table[col].fill_value = np.nan\n", " elif col.startswith(\"flag_\"):\n", " merged_table[col].fill_value = -1\n", "merged_table = merged_table.filled()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Sorting columns\n", "\n", "We sort the columns by increasing band wavelength." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": true }, "outputs": [], "source": [ "bands = [col[2:] for col in merged_table.colnames\n", " if col.startswith('f_') and not col.startswith('f_ap')]\n", "bands.sort(key=lambda x: filter_mean_lambda[x])" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": true }, "outputs": [], "source": [ "columns = ['help_id', 'field', 'ra', 'dec', 'hp_idx', 'ebv', 'redshift', 'zspec']\n", "for band in bands:\n", " for col_tpl in ['f_{}', 'ferr_{}', 'f_ap_{}', 'ferr_ap_{}',\n", " 'm_{}', 'merr_{}', 'm_ap_{}', 'merr_ap_{}',\n", " 'flag_{}']:\n", " colname = col_tpl.format(band)\n", " if colname in merged_table.colnames:\n", " columns.append(colname)\n", "columns += ['cigale_mstar', 'cigale_mstar_err', 'cigale_sfr', 'cigale_sfr_err',\n", " 'cigale_dustlumin', 'cigale_dustlumin_err', 'cigale_dustlumin_ironly', \n", " 'cigale_dustlumin_ironly_err', 'flag_cigale', 'flag_cigale_opt', \n", " 'flag_cigale_ir', 'flag_cigale_ironly', 'stellarity', \n", " 'stellarity_origin', 'flag_cleaned', 'flag_merged', 'flag_gaia', \n", " 'flag_optnir_obs', 'flag_optnir_det', 'zspec_qual', \n", " 'zspec_association_flag']" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "set()\n", "set()\n" ] } ], "source": [ "# Check that we did not forget any column\n", "print( set(columns) - set(merged_table.colnames))\n", "print(set(merged_table.colnames) - set(columns))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Saving" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": true }, "outputs": [], "source": [ "merged_table['help_id'] = merged_table['help_id'].astype('S27')\n", "merged_table = add_column_meta(merged_table, '../columns.yml')\n", "merged_table[columns].write(\"data/ELAIS-N1_{}.fits\".format(SUFFIX), overwrite=True)" ] } ], "metadata": { "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.4" } }, "nbformat": 4, "nbformat_minor": 2 }