dmu32/1_seds_to_vo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write SEDs to the VO\n",
    "\n",
    "Loop over all the photo-z 'best' SEDs in dmu28 and write them all to the appropriate VO folders with the final published columns"
   ]
  },
  {
   "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)\n",
      "This notebook was executed on: \n",
      "2018-07-17 13:19:19.664728\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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.table import Table, Column\n",
    "from astropy import units as u\n",
    "from astropy.io import fits\n",
    "from astropy.coordinates import SkyCoord\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import re\n",
    "\n",
    "import yaml\n",
    "\n",
    "import glob\n",
    "\n",
    "from shutil import copyfile\n",
    "\n",
    "\n",
    "import os\n",
    "import time\n",
    "\n",
    "import humanfriendly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "TODAY = os.environ.get('SUFFIX', time.strftime(\"_%Y%m%d\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert the SEDs\n",
    "Here we define a function to take the SED fits file as output by CIGALE and change it to our final published format. This means replacing the columns thus:\n",
    "\n",
    "\n",
    "| Output column | relation to input columns                                     |\n",
    "|---------------|---------------------------------------------------------------|\n",
    "| WAVE          | wavelength                                                    |\n",
    "| FLUX          | Fnu                                                           |\n",
    "| LUMIN         | L_lambda_total                                                |\n",
    "| STELLAR       | stellar.old + stellar.young                                   |\n",
    "| DUST          | dust.Umin_Umin + dust.Umin_Umax                               |\n",
    "| AGN           | agn.fritz2006_therm + agn.fritz2006_scatt + agn.fritz2006_agn |\n",
    "| IGM           | igm                                                           |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def name_to_ra_dec(help_id):\n",
    "    \"\"\"Take a HELP id and return the ra dec in deg\n",
    "    \n",
    "    parameters\n",
    "    ----------\n",
    "    \n",
    "    help_id: string\n",
    "       The help ID\n",
    "       \n",
    "    returns\n",
    "    -------\n",
    "    \n",
    "    ra: flaot\n",
    "       Right Ascension in degrees\n",
    "       \n",
    "    dec: float\n",
    "        Declination in degrees\n",
    "    \"\"\"\n",
    "    if '-' in help_id:\n",
    "        ra, dec = help_id.strip('HELP_J').split('-')\n",
    "        dec = '-' + dec\n",
    "    if '+' in help_id:\n",
    "        ra, dec = help_id.strip('HELP_J').split('+')\n",
    "        dec = '+' + dec\n",
    "    \n",
    "    ra = ra[:2] + 'h' + ra[2:]\n",
    "    ra = ra[:5] + 'm' + ra[5:]\n",
    "    ra = ra + 's' \n",
    "    \n",
    "    \n",
    "    \n",
    "    dec = dec[:3] + 'd' + dec[3:]\n",
    "    dec = dec[:6] + 'm' + dec[6:]\n",
    "    dec = dec + 's'\n",
    "    \n",
    "    coords = SkyCoord('{} {}'.format(ra, dec))\n",
    "    \n",
    "    return coords.ra.value, coords.dec.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cigale_to_vo(in_sed_name, out_folder, help_id='None', ra=np.nan, dec=np.nan, z=np.nan):\n",
    "    \"\"\"Convert CIGALE output SED to final VO format\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    in_sed_name: string\n",
    "        Link to input SED fits file containg astropy Table\n",
    "        \n",
    "    out_folder: string\n",
    "        Folder to write outputs to\n",
    "        \n",
    "    ob_name: string\n",
    "        Name of object \n",
    "        \n",
    "    ra: float\n",
    "        Right Ascension [deg]\n",
    "        \n",
    "    dec: float\n",
    "        Declination [deg]\n",
    "        \n",
    "    z: float\n",
    "        Redshift\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    out_sed_name: string\n",
    "        Link to output SED fits file that was written during function run.\n",
    "    \"\"\"\n",
    "    try:\n",
    "        in_sed = Table.read(in_sed_name)\n",
    "        \n",
    "        z = in_sed.meta['universe.redshift']\n",
    "        help_id = in_sed_name.split('/')[-1].strip('_best_model.fits')\n",
    "        ra,dec = name_to_ra_dec(help_id)\n",
    "        \n",
    "        out_sed_name = out_folder + in_sed_name.split('/')[-1]\n",
    "        out_sed_name=out_sed_name.replace('+', 'p')\n",
    "        #print(out_sed_name)\n",
    "        out_sed = Table()\n",
    "    \n",
    "        out_sed.add_column(Column(data = in_sed['wavelength'], \n",
    "                                name = 'WAVE'))\n",
    "        out_sed.add_column(Column(data = in_sed['Fnu'], \n",
    "                                name = 'FLUX'))\n",
    "        out_sed.add_column(Column(data = in_sed['L_lambda_total'], \n",
    "                                name = 'LUMIN'))\n",
    "        out_sed.add_column(Column(data = in_sed['stellar.old'] + in_sed['stellar.young'], \n",
    "                                name = 'STELLAR'))\n",
    "        out_sed.add_column(Column(data = in_sed['dust.Umin_Umin'] + in_sed['dust.Umin_Umax'], \n",
    "                                name = 'DUST'))\n",
    "        out_sed.add_column(Column(data = in_sed['agn.fritz2006_therm'] + in_sed['agn.fritz2006_scatt'] + in_sed['agn.fritz2006_agn'], \n",
    "                                name = 'AGN'))\n",
    "        out_sed.add_column(Column(data = in_sed['igm'], \n",
    "                                name = 'IGM'))\n",
    "        \n",
    "    \n",
    "        #Add header\n",
    "\n",
    "        out_sed.meta['VOCLASS'] = 'SPECTRUM 1.0'\n",
    "        out_sed.meta['RADECSYS'] = 'ICRS'\n",
    "        out_sed.meta['EQUINOX'] = 2000.0\n",
    "        out_sed.meta['VOPUB'] = 'Herschel Extragalactic Legacy Programme'\n",
    "        out_sed.meta['OBJECT'] = help_id\n",
    "        out_sed.meta['REDSHIFT'] = z\n",
    "        out_sed.meta['RA_TARG'] = ra\n",
    "        out_sed.meta['DEC_TARG'] = dec\n",
    "        out_sed.meta['TITLE'] = '{} CIGALE spectrum'.format(help_id)\n",
    "        out_sed.write(out_sed_name, overwrite=True)\n",
    "    except FileNotFoundError:\n",
    "        return 'No file found!' \n",
    "    \n",
    "    \n",
    "    \n",
    "    #Also try to copy the png\n",
    "    try:\n",
    "        copyfile(in_sed_name.replace('fits', 'png'), \n",
    "                 out_sed_name.replace('fits', 'png'))\n",
    "\n",
    "    except FileNotFoundError:\n",
    "        pass\n",
    "        \n",
    "    \n",
    "    return out_sed_name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dmu32_ymls = glob.glob('../dmu32/*/meta_main.yml', recursive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dmu_products/dmu28/dmu28_GAMA-15/data/zphot/HELP_final_results.fits\n",
      "../dmu28/dmu28_GAMA-15/data/zphot/SEDs/fits/\n",
      "There are 117234 SEDs here\n",
      "dmu_products/dmu28/dmu28_GAMA-12/data/zphot/HELP_final_results.fits\n",
      "../dmu28/dmu28_GAMA-12/data/zphot/SEDs/fits/\n",
      "There are 108139 SEDs here\n",
      "dmu_products/dmu28/dmu28_COSMOS/data/zphot/final_results.fits\n",
      "../dmu28/dmu28_COSMOS/data/zphot/SEDs/fits/\n",
      "There are 15747 SEDs here\n",
      "dmu_products/dmu28/dmu28_SGP/data/zphot/HELP_final_results.fits\n",
      "../dmu28/dmu28_SGP/data/zphot/SEDs/fits/\n",
      "There are 352804 SEDs here\n",
      "dmu_products/dmu28/dmu28_Lockman-SWIRE/data/zphot/best_extcor/HELP_final_results.fits\n",
      "../dmu28/dmu28_Lockman-SWIRE/data/zphot/best_extcor/SEDs/fits/\n",
      "There are 46719 SEDs here\n",
      "dmu_products/dmu28/dmu28_ELAIS-N1/data/zphot/HELP_final_results.fits\n",
      "../dmu28/dmu28_ELAIS-N1/data/zphot/SEDs/fits/\n",
      "There are 49986 SEDs here\n",
      "dmu_products/dmu28/dmu28_GAMA-09/data/zphot/HELP_final_results.fits\n",
      "../dmu28/dmu28_GAMA-09/data/zphot/SEDs/fits/\n",
      "There are 130293 SEDs here\n",
      "dmu_products/dmu28/dmu28_ELAIS-S1/data/zphot/HELP_final_results.fits\n",
      "../dmu28/dmu28_ELAIS-S1/data/zphot/SEDs/fits/\n",
      "There are 25393 SEDs here\n",
      "dmu_products/dmu28/dmu28_CDFS-SWIRE/data/zphot/HELP_final_results.fits\n",
      "../dmu28/dmu28_CDFS-SWIRE/data/zphot/SEDs/fits/\n",
      "There are 9308 SEDs here\n"
     ]
    }
   ],
   "source": [
    "out_folder = '/mnt/hedam/data_vo/bestseds/{}/'\n",
    "#out_folder = '/Users/rs548/GitHub/gavo_inputs/bestseds/data/{}/'\n",
    "\n",
    "for yml_file in dmu32_ymls:\n",
    "\n",
    "    dmu32_yml = yaml.load(open(yml_file, 'r'))\n",
    "    field = dmu32_yml['field']\n",
    "    cigale = dmu32_yml['cigale']\n",
    "    if cigale != 'None':\n",
    "        print(cigale)\n",
    "        try:\n",
    "            os.mkdir(out_folder.format(field))\n",
    "        except FileExistsError:\n",
    "            pass\n",
    "        except FileNotFoundError:\n",
    "            continue\n",
    "        \n",
    "        sed_folder = re.sub('/*final_results.fits', '/SEDs/fits/', cigale)\n",
    "        sed_folder = re.sub('dmu_products', '..', sed_folder)\n",
    "        sed_folder = re.sub('/HELP_', '', sed_folder)\n",
    "        print(sed_folder)\n",
    "        \n",
    "        all_seds = glob.glob(sed_folder + '*.fits')\n",
    "        print(\"There are {} SEDs here\".format(len(all_seds)))\n",
    "        for in_sed_name in all_seds:\n",
    "            cigale_to_vo(in_sed_name, out_folder.format(field))\n",
    "\n",
    "  "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}