{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "The PACS maps for HELP have been made by Matt Smith and Kyle Penner. Unfortunetly, the formats differ between the two. For instance Kyle's noise maps are variance rather than standard deviation maps and units differ. This script creates an homogenised format, and attaches other relevant meta data and the ObsIDs." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Steps:\n", "for each field:\n", "* read in coverage, image and noise\n", "* if Kyle's convert units and correct noise\n", "* attach obs IDs\n", "* attach meta data\n", "* create PSF extension" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from datetime import datetime\n", "from itertools import product\n", "\n", "import numpy as np\n", "\n", "from astropy.io import fits\n", "from astropy.table import Table\n", "VERSION = \"0.9\"" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "filelist=Table.read('./pacs_maps.csv', format='csv')" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "all_obsids=Table.read('./pacs_obsid.csv', format='csv')" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'np' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfield\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munique\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilelist\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'field'\u001b[0m\u001b[0;34m]\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[0m\u001b[1;32m 2\u001b[0m \u001b[0;31m#print(filelist[filelist['field']==field])\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'np' is not defined" ] } ], "source": [ "for i, field in enumerate(np.unique(filelist['field'])):\n", " #print(filelist[filelist['field']==field])\n", " \n", " \n", " \n", " obsids = sorted(list(all_obsids['ObsID'][\n", " all_obsids['field'] == field]))\n", " \n", " #loop over band\n", " for b in [100, 160]:\n", " \n", " # to get around Herve's annoying change in number of extensions for cosmos\n", " if field == 'COSMOS':\n", " ext=1\n", " else:\n", " ext=0\n", " \n", " \n", " # Load in image map\n", " image_map=fits.open(f\"/Volumes/Untitled 2/HELP/hedam.lam.fr/HELP/dataproducts/dmu18/dmu18_{field}/data/{filelist[(filelist['field']==field) & (filelist['band']==b) & (filelist['type']=='image')]['file'][0]}\")\n", " image_hdu = fits.ImageHDU(header=image_map[ext].header,\n", " data=image_map[ext].data)\n", " image_hdu.header['EXTNAME'] = \"IMAGE\"\n", " #convert unit if required\n", "\n", " if image_hdu.header['BUNIT'] =='Jy/pix':\n", " image_hdu.data=image_hdu.data*1.0/((3.0462E2)*image_hdu.header['CDELT1']**2)\n", " image_hdu.header['BUNIT']='MJy/sr'\n", " \n", " #to get around the fact Herve didnt bother adding units for COSMOS\n", " if image_hdu.header['BUNIT'] =='1':\n", " print('COSMOS now has units')\n", " image_hdu.header['BUNIT']='MJy/sr'\n", " \n", " #Load in noise map\n", " noise_map=fits.open(f\"/Volumes/Untitled 2/HELP/hedam.lam.fr/HELP/dataproducts/dmu18/dmu18_{field}/data/{filelist[(filelist['field']==field) & (filelist['band']==b) & (filelist['type']=='noise')]['file'][0]}\")\n", " noise_hdu = fits.ImageHDU(header=noise_map[ext].header,\n", " data=noise_map[ext].data)\n", " noise_hdu.header['EXTNAME'] = \"ERROR\"\n", " #convert unit if required\n", " if noise_hdu.header['BUNIT'] =='Jy/pix':\n", " noise_hdu.data=noise_hdu.data*1.0/((3.0462E2)*noise_hdu.header['CDELT1']**2)\n", " noise_hdu.header['BUNIT']='MJy/sr'\n", " #to get around the fact Herve didnt bother adding units for COSMOS\n", " if noise_hdu.header['BUNIT'] =='1':\n", " print('COSMOS now has units')\n", " noise_hdu.header['BUNIT']='MJy/sr'\n", " \n", " #Check if noise is standard deviation by comparing median of noise to standard deviation of image.\n", " #If not, then correct it\n", " if np.std(image_hdu.data[np.isfinite(image_hdu.data)]) < np.median(noise_hdu.data[np.isfinite(noise_hdu.data)]):\n", " noise_hdu.data=np.sqrt(noise_hdu.data)\n", " \n", " \n", " #Load in coverage map. Not all maps have coverage maps so check\n", " try:\n", " cov_map=fits.open(f\"/Volumes/Untitled 2/HELP/hedam.lam.fr/HELP/dataproducts/dmu18/dmu18_{field}/data/{filelist[(filelist['field']==field) & (filelist['band']==b) & (filelist['type']=='cov')]['file'][0]}\")\n", " cov_hdu = fits.ImageHDU(header=cov_map[0].header,\n", " data=cov_map[0].data)\n", " cov_hdu.header['EXTNAME'] = \"COVERAGE\"\n", " #convert unit if required\n", " if cov_hdu.header['BUNIT'] =='Jy/pix':\n", " cov_hdu.data=cov_hdu.data*1.0/((3.0462E2)*cov_hdu.header['CDELT1']**2)\n", " cov_hdu.header['BUNIT']='MJy/sr'\n", " except(IndexError):\n", " cov_hdu = fits.ImageHDU()\n", " cov_hdu.header['EXTNAME'] = \"COVERAGE\"\n", " cov_hdu.header.add_comment(\"The coverage map is not available.\")\n", " \n", " \n", " \n", " \n", " #Create primary header\n", " primary_hdu = fits.PrimaryHDU()\n", " primary_hdu.header.append((\n", " \"CREATOR\", \"Herschel Extragalactic Legacy Project\"\n", " ))\n", " primary_hdu.header.append((\n", " \"TIMESYS\", \"UTC\", \"All dates are in UTC time\"\n", " ))\n", " primary_hdu.header.append((\n", " \"DATE\", datetime.now().replace(microsecond=0).isoformat(),\n", " \"Date of file creation\"\n", " ))\n", " primary_hdu.header.append((\n", " \"VERSION\", VERSION, \"HELP product version\"\n", " ))\n", " primary_hdu.header.append((\n", " \"TELESCOP\", \"Herschel\", \"Name of the telescope\"\n", " ))\n", " primary_hdu.header.append((\n", " \"INSTRUME\", \"PACS\", \"Name of the instrument\"\n", " ))\n", " primary_hdu.header.append((\n", " \"FILTER\", f\"PACS-{str(b)}\", \"Name of the filter\"\n", " ))\n", " primary_hdu.header.append((\n", " \"FIELD\", field, \"Name of the HELP field\"\n", " ))\n", " for idx, obs_id in enumerate(obsids):\n", " keyword = \"OBSID\" + str(idx).zfill(3)\n", " primary_hdu.header.append((keyword, obs_id))\n", " \n", " \n", " #Add empty extension for PSF\n", " psf_hdu = fits.ImageHDU()\n", " psf_hdu.header['EXTNAME'] = \"PSF\"\n", " psf_hdu.header.add_comment(\"The PSF is not available.\")\n", " \n", " hdu_list = fits.HDUList([primary_hdu, image_hdu, noise_hdu, cov_hdu, psf_hdu])\n", " hdu_list.writeto(f\"data/{field}_PACS{str(b)}_v{VERSION}.fits\",\n", " checksum=True, overwrite=True)\n", " # print(f\"{field} / {b} processed...\")\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "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.4" } }, "nbformat": 4, "nbformat_minor": 2 }