Started work on selection func on Lockman-SWIRE

These mean errors can be used to calculate probabilities of a given object being detected and possibilty generate a pdf of true flux given measured flux and a model of number counts #9

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 <i>per object</i>."
+   ]
+  },
+  {
+   "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
+}