{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# XID+ prior model for blind catalogues\n", "\n", "![HELP LOGO](https://avatars1.githubusercontent.com/u/7880370?s=100&v=4>)\n", "\n", "\n", "The final processing stage requires:\n", "1. Create xid+ prior files for all fields\n", "2. Description on runnning xid+\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/Steven/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: The mpl_toolkits.axes_grid module was deprecated in version 2.1. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist provies the same functionality instead.\n", " warnings.warn(message, mplDeprecation, stacklevel=1)\n" ] } ], "source": [ "import pylab\n", "import pymoc\n", "import xidplus\n", "\n", "import numpy as np\n", "%matplotlib inline\n", "from astropy.table import Table\n", "import os\n", "import seaborn as sns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook uses all the raw data from the Blind SPIRE catalogue to create XID+ prior object and relevant tiling scheme" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in XID+SPIRE catalogue" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SA13_SPIRE\n" ] } ], "source": [ "name_field= ['AKARI-SEP_SPIRE','EGS_SPIRE','ELAIS-S1_SPIRE','SA13_SPIRE','XMM-13hr_SPIRE',\\\n", " 'COSMOS_SPIRE','ELAIS-N2_SPIRE','HDF-N_SPIRE','SPIRE-NEP_SPIRE','xFLS_SPIRE','ELAIS-N1_SPIRE',\\\n", " 'XMM-LSS_SPIRE','Lockman-SWIRE_SPIRE','CDFS-SWIRE_SPIRE', 'GAMA-09_SPIRE','GAMA-12_SPIRE',\\\n", " 'GAMA-15_SPIRE','SSDF_SPIRE','Bootes_SPIRE','HATLAS-NGP_SPIRE','HATLAS-SGP_SPIRE','AKARI-NEP_SPIRE']\n", "which_field = 3\n", "name_field = name_field[which_field] # select a field, in this case SA13\n", "print(name_field)\n", "loc = 'dmu22_'+name_field[0:-6]+'/data/'\n", "XID_table=Table.read(loc+name_field+'_all.fits')" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=10\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
RADecF_BLIND_MF_SPIRE_250FErr_BLIND_MF_SPIRE_250F_BLIND_MF_SPIRE_350FErr_BLIND_MF_SPIRE_350F_BLIND_MF_SPIRE_500FErr_BLIND_MF_SPIRE_500rPRA_pixDec_pixF_BLIND_pix_SPIREFErr_BLIND_pix_SPIREflag
degdegJyJyJyJyJyJydegdegJyJy
float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64
198.0633796142.7772856940.1749214630620.004116980692740.08483415681040.003837016298640.02662200195380.004422679457670.8803060193811.0198.062244442.77811912160.1523224962080.004032605971640.0
198.32940754342.7733430720.1189758700060.005403196408650.05458381472230.005659005076260.01754280125070.0055632664140.7846320107951.0198.33016933342.77445234740.1309948257070.006050677334620.0
198.0063990442.55922147740.1132737837980.004066338841280.07932375875660.004103916730760.04440536024420.00426889525430.9118685421591.0198.00564438142.55977671960.1119217837350.003987072402690.0
197.99179184342.83976958350.08967807580990.003979244982750.05726632339010.003898969815980.03512359684330.004739984902230.8477153052861.0197.99179184342.83976958350.08938605838270.003905327632680.0
198.34259837742.59969984430.07620235762270.00565829876140.03804125049120.006853150108160.02172573309370.006529343549660.8095718762681.0198.34297443842.59942112080.08734100275860.005550155943390.0
198.26716176842.82930898820.08177526156380.005319351985850.02233105622850.005821530191350.003087966382230.006173190543220.6921363387871.0198.26678396842.82958746760.0794898393120.005217857956320.0
198.02515237242.78700494360.07552919995020.003828648268690.0442331466780.003743778642380.01746502905690.003927393009450.836727502941.0198.02590962342.7864495710.07412386677370.003756105478060.0
197.85514999342.56711602120.06856850280670.00515664573930.06083904121020.005742533148920.05445414635140.006631577065650.8933208806451.0197.8540158442.56794739040.0737391952980.0050562206990.0
198.26473341342.8945907820.07552865428110.007538084227570.02043141749330.007605852873130.003891193557370.006734420647520.7078056674341.0198.26473341342.8945907820.07187014873370.007368869401020.0
197.93872088642.65389590620.07475329153570.003888329959030.02407165695640.003767503364670.003746431736610.004040928101270.7176607392071.0197.93758626242.65472810360.07174016249250.003816528263730.0
" ], "text/plain": [ "\n", " RA Dec ... FErr_BLIND_pix_SPIRE flag \n", " deg deg ... Jy \n", " float64 float64 ... float64 float64\n", "------------- ------------- ... -------------------- -------\n", " 198.06337961 42.777285694 ... 0.00403260597164 0.0\n", "198.329407543 42.773343072 ... 0.00605067733462 0.0\n", " 198.00639904 42.5592214774 ... 0.00398707240269 0.0\n", "197.991791843 42.8397695835 ... 0.00390532763268 0.0\n", "198.342598377 42.5996998443 ... 0.00555015594339 0.0\n", "198.267161768 42.8293089882 ... 0.00521785795632 0.0\n", "198.025152372 42.7870049436 ... 0.00375610547806 0.0\n", "197.855149993 42.5671160212 ... 0.005056220699 0.0\n", "198.264733413 42.894590782 ... 0.00736886940102 0.0\n", "197.938720886 42.6538959062 ... 0.00381652826373 0.0" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "XID_table[0:10]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The uncertianties become Gaussian by $\\sim 20 \\mathrm{\\mu Jy}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in Maps" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "loc = '../dmu19/dmu19_HELP-SPIRE-maps/data/'\n", "name = [name_field+'250_v1.0.fits',name_field+'350_v1.0.fits',name_field+'500_v1.0.fits']\n", "pswfits=loc+name[0]\n", "pmwfits=loc+name[1]\n", "plwfits=loc+name[2]\n", "os.mkdir('./OUT_'+name_field[0:-6])\n", "output_folder ='./OUT_'+name_field[0:-6]+'/'\n" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "from astropy.io import fits\n", "from astropy import wcs\n", "\n", "#-----250-------------\n", "hdulist = fits.open(pswfits)\n", "im250phdu=hdulist[0].header\n", "im250hdu=hdulist[1].header\n", "\n", "im250=hdulist['IMAGE'].data*1.0E3 #convert to mJy\n", "nim250=hdulist['ERROR'].data*1.0E3 #convert to mJy\n", "w_250 = wcs.WCS(hdulist[1].header)\n", "pixsize250 = hdulist[\"Matchedfilter\"].header[\"PIXSIZE\"] \n", "\n", "hdulist.close()\n", "#-----350-------------\n", "hdulist = fits.open(pmwfits)\n", "im350phdu=hdulist[0].header\n", "im350hdu=hdulist[1].header\n", "\n", "im350=hdulist['IMAGE'].data*1.0E3 #convert to mJy\n", "nim350=hdulist['ERROR'].data*1.0E3 #convert to mJy\n", "w_350 = wcs.WCS(hdulist[1].header)\n", "pixsize350 = hdulist[\"Matchedfilter\"].header[\"PIXSIZE\"] \n", "hdulist.close()\n", "#-----500-------------\n", "hdulist = fits.open(plwfits)\n", "im500phdu=hdulist[0].header\n", "im500hdu=hdulist[1].header\n", "im500=hdulist['IMAGE'].data*1.0E3 #convert to mJy\n", "nim500=hdulist['ERROR'].data*1.0E3 #convert to mJy\n", "w_500 = wcs.WCS(hdulist[1].header)\n", "pixsize500 = hdulist[\"Matchedfilter\"].header[\"PIXSIZE\"] \n", "\n", "hdulist.close()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "## Set XID+ prior class\n", "ID = np.arange(1,np.size(XID_table['RA'])+1)\n", "ID = ID.astype(str)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "#---prior250--------\n", "prior250=xidplus.prior(im250,nim250,im250phdu,im250hdu)#Initialise with map, uncertianty map, wcs info and primary header\n", "prior250.prior_cat(XID_table['RA'],XID_table['Dec'],name_field+'_250_XID.fits',ID=ID)#Set input catalogue\n", "prior250.prior_bkg(-5.0,5)#Set prior on background (assumes Gaussian pdf with mu and sigma)\n", "#---prior350--------\n", "prior350=xidplus.prior(im350,nim350,im350phdu,im350hdu)\n", "prior350.prior_cat(XID_table['RA'],XID_table['Dec'],name_field+'_350_XID.fits',ID=ID)#Set input catalogue\n", "prior350.prior_bkg(-5.0,5)\n", "\n", "#---prior500--------\n", "prior500=xidplus.prior(im500,nim500,im500phdu,im500hdu)\n", "prior500.prior_cat(XID_table['RA'],XID_table['Dec'],name_field+'_500_XID.fits',ID=ID)#Set input catalogue\n", "prior500.prior_bkg(-5.0,5)\n" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "#pixsize array (size of pixels in arcseconds)\n", "pixsize=np.array([pixsize250,pixsize350,pixsize500])\n", "#point response function for the three bands\n", "prfsize=np.array([18.15,25.15,36.3])\n", "#use Gaussian2DKernel to create prf (requires stddev rather than fwhm hence pfwhm/2.355)\n", "from astropy.convolution import Gaussian2DKernel\n", "\n", "##---------fit using Gaussian beam-----------------------\n", "prf250=Gaussian2DKernel(prfsize[0]/2.355,x_size=101,y_size=101)\n", "prf250.normalize(mode='peak')\n", "prf350=Gaussian2DKernel(prfsize[1]/2.355,x_size=101,y_size=101)\n", "prf350.normalize(mode='peak')\n", "prf500=Gaussian2DKernel(prfsize[2]/2.355,x_size=101,y_size=101)\n", "prf500.normalize(mode='peak')\n", "\n", "pind250=np.arange(0,101,1)*1.0/pixsize[0] #get 250 scale in terms of pixel scale of map\n", "pind350=np.arange(0,101,1)*1.0/pixsize[1] #get 350 scale in terms of pixel scale of map\n", "pind500=np.arange(0,101,1)*1.0/pixsize[2] #get 500 scale in terms of pixel scale of map\n", "\n", "prior250.set_prf(prf250.array,pind250,pind250)#requires psf as 2d grid, and x and y bins for grid (in pixel scale)\n", "prior350.set_prf(prf350.array,pind350,pind350)\n", "prior500.set_prf(prf500.array,pind500,pind500)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "----- There are 5 tiles required for input catalogue and 2 large tiles\n" ] } ], "source": [ "import pickle\n", "#from moc, get healpix pixels at a given order\n", "from xidplus import moc_routines\n", "\n", "order=7 \n", "tiles=moc_routines.get_HEALPix_pixels(order,prior250.sra,prior250.sdec,unique=True)\n", "order_large=6\n", "tiles_large=moc_routines.get_HEALPix_pixels(order_large,prior250.sra,prior250.sdec,unique=True)\n", "print('----- There are '+str(len(tiles))+' tiles required for input catalogue and '+str(len(tiles_large))+' large tiles')\n", "#output_folder='./'\n", "outfile=output_folder+'Master_prior.pkl'\n", "\n", "with open(outfile, 'wb') as f:\n", " xidplus.io.pickle.dump({'priors':[prior250,prior350,prior500],'tiles':tiles,'order':order,'version':xidplus.io.git_version()},f)\n", " \n", "outfile=output_folder+'Tiles.pkl'\n", "with open(outfile, 'wb') as f:\n", " pickle.dump({'tiles':tiles,'order':order,'tiles_large':tiles_large,'order_large':order_large,'version':xidplus.io.git_version()},f)\n" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "./OUT_SA13/\n" ] } ], "source": [ "print(output_folder)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## The Master_prior.pkl and Tiles.pkl are used to run XID+\n", "# Running on Apollo\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "% mv XID\\_plus\\_hier.sh OUT\\_\"name_field\"\n", "\n", "% cd OUT_\"name_field\"\n", "\n", "% module load sge\n", "\n", "% qsub -t 1-$n_hier -q mps.q -jc mps.short XID_plus_hier.sh \n", "\n", "### n_hier is the number of large tiles\n", "\n", "% cd ..\n", "\n", "% qsub -t 1-$n_tiles -pe openmp 4 -l h_rt=6:00:00 -l m_mem_free=13G -q mps.q XID_plus_tile.sh\n", "### n_hier is the number of small tiles\n", "### Then combine the Bayesian maps into one:\n", "\n", "% python make_combined_map.py\n", "\n", "% cd output\n", "\n", "% ls *cat.fits > cat_files \n", "\n", "% module load starlink/hikianalia-64bit\n", "\n", "% stilts tcat ifmt=fits in=@cat_files out=dmu22_XID+SPIRE_\"name_field\"_BLIND.fits\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## The data products are:\n", "### dmu22XID+SPIRE\"name_field\"_BLIND.fits\n", "#### dmu22_XID+SPIRE\\_psw_\"name_field\"_Bayes_Pval\n", "#### dmu22_XID+SPIRE\\_pmw_\"name_field\"_Bayes_Pval\n", "#### dmu22_XID+SPIRE\\_plw_\"name_field\"_Bayes_Pval" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Final validation of the data can be found at: \n", "# http://hedam.lam.fr/HELP/dataproducts/dmu22/\n", "# or\n", "# https://github.com/H-E-L-P/dmu_products/tree/master/dmu22\n", "# dmu22_\"name\\_field\"/XID+BLIND_\"name_field\"_final_processing.ipynb" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "*This is a default HELP jupyter notebook *\n", "\n", " ![HELP LOGO](https://avatars1.githubusercontent.com/u/7880370?s=75&v=4)\n", "\n", "**Authors**: S. Duivenvoorden, P.D. Hurely\n", "\n", " \n", "For a full description of the database and how it is organised in to `dmu_products` please the top level [readme](../readme.md).\n", " \n", "The Herschel Extragalactic Legacy Project, ([HELP](http://herschel.sussex.ac.uk/)), is a [European Commission Research Executive Agency](https://ec.europa.eu/info/departments/research-executive-agency_en)\n", "funded project under the SP1-Cooperation, Collaborative project, Small or medium-scale focused research project, FP7-SPACE-2013-1 scheme, Grant Agreement\n", "Number 607254.\n", "\n", "[Acknowledgements](http://herschel.sussex.ac.uk/acknowledgements)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "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.5" } }, "nbformat": 4, "nbformat_minor": 2 }