diff --git a/dmu1/dmu1_ml_Lockman-SWIRE/4_Selection_function.ipynb b/dmu1/dmu1_ml_Lockman-SWIRE/4_Selection_function.ipynb new file mode 100644 index 00000000..2e8cba6c --- /dev/null +++ b/dmu1/dmu1_ml_Lockman-SWIRE/4_Selection_function.ipynb @@ -0,0 +1,336 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Lockman-SWIRE Selection Functions\n", + "## Depth maps and selection functions for Lockman-SWIRE\n", + "\n", + "The simplest selection function available is the field MOC which specifies the area for which there is Herschel data. Each pristine catalogue also has a MOC defining the area for which that data is available.\n", + "\n", + "The next stage is to provide mean flux standard deviations which act as a proxy for the catalogue's 5$\\sigma$ depth" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "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", + "import os\n", + "import time\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", + "from pymoc import MOC\n", + "import healpy as hp\n", + "#import pandas as pd #Astropy has group_by function so apandas isn't required.\n", + "import seaborn as sns\n", + "\n", + "import warnings\n", + "#We ignore warnings - this is a little dangerous but a huge number of warnings are generated by empty cells later\n", + "warnings.filterwarnings('ignore')\n", + "\n", + "from herschelhelp_internal.utils import inMoc, coords_to_hpidx\n", + "from herschelhelp_internal.masterlist import find_last_ml_suffix, nb_ccplots" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "TMP_DIR = os.environ.get('TMP_DIR', \"./data_tmp\")\n", + "OUT_DIR = os.environ.get('OUT_DIR', \"./data\")\n", + "SUFFIX = find_last_ml_suffix()\n", + "#SUFFIX = \"20171016\"\n", + "\n", + "master_catalogue_filename = \"master_catalogue_lockman-swire_{}.fits\".format(SUFFIX)\n", + "master_catalogue = Table.read(\"{}/{}\".format(OUT_DIR, master_catalogue_filename))\n", + "\n", + "print(\"Diagnostics done using: {}\".format(master_catalogue_filename))\n", + "\n", + "ORDER = 10\n", + "#TODO write code to decide on appropriate order\n", + "\n", + "field_moc = MOC(filename=\"../../dmu2/dmu2_field_coverages/Lockman-SWIRE_MOC.fits\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plotting a histogram of magnitudes gives a quick idea of depth. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "sns.distplot(master_catalogue['m_wfc_g'][(master_catalogue['m_wfc_g'] > 0) & (master_catalogue['m_wfc_g'] < 90.)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and for flux" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sns.distplot(master_catalogue['f_wfc_g'][(master_catalogue['f_wfc_g'] < 10.) & (master_catalogue['f_wfc_g'] > 0)] )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Distribution of flux errors:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sns.distplot(master_catalogue['ferr_wfc_g'][(master_catalogue['ferr_wfc_g'] > 0.) &(master_catalogue['ferr_wfc_g'] < 3.)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Magnitude error as function of magnitude" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ax = sns.jointplot((np.log10(master_catalogue['m_wfc_g'])),\n", + " (master_catalogue['merr_wfc_g'] ))\n", + "#ax.set(xticklabels=)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## I - Group masterlist objects by healpix cell and calculate depths\n", + "We add a column to the masterlist catalogue for the target order healpix cell per object." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "#Add a column to the catalogue with the order=ORDER hp_idx\n", + "master_catalogue.add_column(Column(data=coords_to_hpidx(master_catalogue['ra'],\n", + " master_catalogue['dec'],\n", + " ORDER), \n", + " name=\"hp_idx_O_{}\".format(str(ORDER))\n", + " )\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Convert catalogue to pandas and group by the order=ORDER pixel\n", + "\n", + "group = master_catalogue.group_by([\"hp_idx_O_{}\".format(str(ORDER))])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "f_wfc_g_p90 = group[\"hp_idx_O_{}\".format(str(ORDER)), 'f_wfc_g'].groups.aggregate(lambda x: np.nanpercentile(x, 10.))\n", + "f_wfc_g_p90['f_wfc_g'].name = 'f_wfc_g_p90'\n", + "f_wfc_g_p90[:10].show_in_notebook()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ferr_wfc_g_mean = group[\"hp_idx_O_{}\".format(str(ORDER)), 'ferr_wfc_g'].groups.aggregate(np.nanmean)\n", + "ferr_wfc_g_mean['ferr_wfc_g'].name = 'ferr_wfc_g_mean'\n", + "ferr_wfc_g_mean[:10].show_in_notebook()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "#Downgrade the groups from order=ORDER to order=13 and then fill out the appropriate cells\n", + "#hp.pixelfunc.ud_grade([2599293, 2599294], nside_out=hp.order2nside(13))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## II Create a table of all Order=13 healpix cells in the field and populate it\n", + "We create a table with every order=13 healpix cell in the field MOC. We then calculate the healpix cell at lower order that the order=13 cell is in. We then fill in the depth at every order=13 cell as calculated for the lower order cell that that the order=13 cell is inside." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "depths = Table()\n", + "depths['hp_idx_O_13'] = list(field_moc.flattened(13))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "depths[:10].show_in_notebook()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "depths.add_column(hp.pixelfunc.ang2pix(2**ORDER,\n", + " hp.pixelfunc.pix2ang(2**13, depths['hp_idx_O_13'], nest=True)[0],\n", + " hp.pixelfunc.pix2ang(2**13, depths['hp_idx_O_13'], nest=True)[1],\n", + " nest = True),\n", + " name=\"hp_idx_O_{}\".format(str(ORDER))\n", + " )\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "depths[:10].show_in_notebook()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "join(depths, ferr_wfc_g_mean)[:10].show_in_notebook()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for col in master_catalogue.colnames:\n", + " if col.startswith(\"f_\"):\n", + " errcol = \"ferr{}\".format(col[1:])\n", + " depths = join(depths, \n", + " group[\"hp_idx_O_{}\".format(str(ORDER)), errcol].groups.aggregate(np.nanmean),\n", + " join_type='left')\n", + " depths[errcol].name = errcol + \"_mean\"\n", + " depths = join(depths, \n", + " group[\"hp_idx_O_{}\".format(str(ORDER)), col].groups.aggregate(lambda x: np.nanpercentile(x, 90.)),\n", + " join_type='left')\n", + " depths[col].name = col + \"_p90\"\n", + "\n", + "depths[:10].show_in_notebook()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## III - Save the table of depth maps" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "depths.write(\"{}/depths_lockman-swire{}.fits\".format(OUT_DIR, SUFFIX))" + ] + } + ], + "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.1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}