{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Selection Functions\n", "## Depth maps and selection functions\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", "This notebook should determine the correct field and filenames based on being placed in the correct folder.\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": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This notebook was run with herschelhelp_internal version: \n", "017bb1e (Mon Jun 18 14:58:59 2018 +0100) [with local modifications]\n", "This notebook was executed on: \n", "2020-12-01 11:55:20.600299\n" ] } ], "source": [ "from herschelhelp_internal import git_version\n", "print(\"This notebook was run with herschelhelp_internal version: \\n{}\".format(git_version()))\n", "import datetime\n", "print(\"This notebook was executed on: \\n{}\".format(datetime.datetime.now()))" ] }, { "cell_type": "code", "execution_count": 2, "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", "import glob\n", "import gc\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, flux_to_mag\n", "from herschelhelp_internal.masterlist import find_last_ml_suffix, nb_ccplots\n", "\n", "from astropy.io.votable import parse_single_table\n", "import yaml\n", "import time" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "TODAY = time.strftime(\"%Y%m%d\")" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'20210118'" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "TODAY" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "def make_depth_map(meta_yml):\n", " meta = yaml.load(open(meta_yml, 'r'))\n", " FIELD = meta['field']\n", " #FILTERS_DIR = \"/Users/rs548/GitHub/herschelhelp_python/database_builder/filters/\"\n", " FILTERS_DIR = \"/opt/herschelhelp_python/database_builder/filters/\"\n", "\n", " OUT_DIR = meta_yml.replace('meta_main.yml','data')\n", " #SUFFIX = find_last_ml_suffix()\n", " SUFFIX = meta['final'].split('/')[-1].split('_')[-1].strip('.fits')\n", "\n", " master_catalogue_filename = meta['final'].split('/')[-1]\n", " master_catalogue = Table.read(\"{}/{}\".format(OUT_DIR, master_catalogue_filename))\n", "\n", " for col in master_catalogue.colnames:\n", " if (\n", " col.startswith('m')\n", " or col.startswith('flag')\n", " or col =='redshift'\n", " or col=='zspec'\n", " \n", " ):\n", " master_catalogue.remove_column(col)\n", " \n", " print(\"Depth maps produced using: {}\".format(master_catalogue_filename))\n", "\n", " ORDER = 10\n", " \n", "\n", " field_moc = MOC(filename=\"../../dmu2/dmu2_field_coverages/{}_MOC.fits\".format(FIELD))\n", " # Remove sources whose signal to noise ratio is less than five as these will have been selected using forced \n", " # photometry and so the errors will not refelct the RMS of the map \n", " # For xid+ saource use s/n >2 as that is criteria we use for CIGALE\n", " for n,col in enumerate(master_catalogue.colnames):\n", " if col.startswith(\"f_spire\") or col.startswith(\"f_pacs\")or col.startswith(\"f_mips\"):\n", " err_col = \"ferr{}\".format(col[1:])\n", " errs = master_catalogue[err_col]\n", " fluxes = master_catalogue[col]\n", " mask = fluxes/errs < 2.0\n", " master_catalogue[col][mask] = np.nan\n", " master_catalogue[err_col][mask] = np.nan\n", " elif col.startswith(\"f_\"):\n", " err_col = \"ferr{}\".format(col[1:])\n", " errs = master_catalogue[err_col]\n", " fluxes = master_catalogue[col]\n", " mask = fluxes/errs < 5.0\n", " master_catalogue[col][mask] = np.nan\n", " master_catalogue[err_col][mask] = np.nan\n", " \n", " \n", " #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", " )\n", " # Convert catalogue to pandas and group by the order=ORDER pixel\n", " group = master_catalogue.group_by([\"hp_idx_O_{}\".format(str(ORDER))])\n", " depths = Table()\n", " depths['hp_idx_O_13'] = list(field_moc.flattened(13))\n", " depths.add_column(Column(data=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", " )\n", " 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", " standard = ['u', 'g', 'r', 'i', 'z', 'y', 'j', 'h', 'k', 'ks']\n", " tot_bands = [column[2:] for column in master_catalogue.colnames \n", " if (column.startswith('f_') & ~column.startswith('f_ap_'))]\n", " print('tot_bands', tot_bands)\n", " ap_bands = [column[5:] for column in master_catalogue.colnames \n", " if column.startswith('f_ap_') ]\n", " print('ap_bands', ap_bands)\n", " bands = set(tot_bands) | set(ap_bands)\n", " bands\n", " \n", " for col in depths.colnames:\n", " if depths[col].dtype == 'float64' or depths[col].dtype == 'float32':\n", " depths[col].fill_value = np.nan\n", " \n", " depths = depths.filled()\n", " \n", " for band in standard:\n", " print(band)\n", " #ap_bands\n", " list_of_ap_cols = np.array([ b if b.endswith('_{}'.format(band)) else None for b in ap_bands])\n", " list_of_ap_cols = list_of_ap_cols[list_of_ap_cols != None]\n", " list_of_ap_cols = ['ferr_ap_{}_mean'.format(b) for b in list_of_ap_cols]\n", " print(list_of_ap_cols)\n", " if len(list_of_ap_cols) >0:\n", " depths['ferr_ap_{}_mean'.format(band)] = np.nanmin([depths[t] for t in list_of_ap_cols], axis=0)\n", " #tot_bands\n", " list_of_tot_cols = np.array([ b if b.endswith('_{}'.format(band)) else None for b in tot_bands])\n", " list_of_tot_cols = list_of_tot_cols[list_of_tot_cols != None]\n", " list_of_tot_cols = ['ferr_{}_mean'.format(b) for b in list_of_tot_cols]\n", " print(list_of_tot_cols)\n", " if len(list_of_tot_cols) >0:\n", " depths['ferr_{}_mean'.format(band)] = np.nanmin([depths[t] for t in list_of_tot_cols], axis=0)\n", " \n", " depth_filename = \"../dmu32_{}/data/depths_{}_{}.fits\".format(\n", " FIELD.replace('HATLAS-',''), \n", " FIELD.lower(), \n", " SUFFIX )\n", " depths.write(depth_filename, overwrite=True)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "meta_files = glob.glob('../*/meta_main.yml')" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['../dmu32_COSMOS/meta_main.yml',\n", " '../dmu32_ELAIS-S1/meta_main.yml',\n", " '../dmu32_GAMA-09/meta_main.yml',\n", " '../dmu32_Lockman-SWIRE/meta_main.yml',\n", " '../dmu32_ELAIS-N2/meta_main.yml',\n", " '../dmu32_Herschel-Stripe-82/meta_main.yml',\n", " '../dmu32_SA13/meta_main.yml',\n", " '../dmu32_GAMA-12/meta_main.yml',\n", " '../dmu32_GAMA-15/meta_main.yml',\n", " '../dmu32_AKARI-SEP/meta_main.yml',\n", " '../dmu32_SGP/meta_main.yml',\n", " '../dmu32_xFLS/meta_main.yml',\n", " '../dmu32_CDFS-SWIRE/meta_main.yml',\n", " '../dmu32_SSDF/meta_main.yml',\n", " '../dmu32_XMM-LSS/meta_main.yml',\n", " '../dmu32_SPIRE-NEP/meta_main.yml',\n", " '../dmu32_AKARI-NEP/meta_main.yml',\n", " '../dmu32_ELAIS-N1/meta_main.yml',\n", " '../dmu32_NGP/meta_main.yml',\n", " '../dmu32_XMM-13hr/meta_main.yml',\n", " '../dmu32_Bootes/meta_main.yml',\n", " '../dmu32_HDF-N/meta_main.yml',\n", " '../dmu32_EGS/meta_main.yml']" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "meta_files" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "../dmu32_COSMOS/meta_main.yml failed\n", "../dmu32_ELAIS-S1/meta_main.yml failed\n", "../dmu32_GAMA-09/meta_main.yml failed\n", "../dmu32_Lockman-SWIRE/meta_main.yml failed\n" ] }, { "ename": "KeyboardInterrupt", "evalue": "", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0mfailures\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0mgc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcollect\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msleep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 10\u001b[0m \u001b[0mfailures\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mKeyboardInterrupt\u001b[0m: " ] } ], "source": [ "failures = set()\n", "for m in meta_files:\n", " try:\n", " make_depth_map(m)\n", " except:\n", " print(m, 'failed')\n", " failures.add(m)\n", " gc.collect()\n", " time.sleep(10)\n", "failures" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.10" } }, "nbformat": 4, "nbformat_minor": 2 }