{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SPIRE Aperture Photometry for Point Sources\n",
    "\n",
    "The method is for photometry is outlined in [the SPIRE Handbook](http://herschel.esac.esa.int/Docs/SPIRE/spire_handbook.pdf), section 5.12.1. This notebook just shows one practical example.\n",
    "\n",
    "The workflow is the following:\n",
    "1. Using `ESASky` module of `astroquery` search for the source in _Herschel_ observations. Pick one observation for the example.\n",
    "2. Download the map using the [astroquery.esasky](https://astroquery.readthedocs.io/en/latest/esasky/esasky.html) results.\n",
    "3. Perform aperture photometry using [photutils](https://photutils.readthedocs.io/en/stable).\n",
    "4. Bonus: compare the aperture-derived flux density with the one in the SPIRE point source catalogue, also available in `astroquery.esasky`.\n",
    "\n",
    "For this exercise we are going to use a known stellar calibrator: $\\beta$ And, with RA=01:09:43.9, Dec= +35:37:14.\n",
    "\n",
    "For reference, the model flux densities of $\\beta$ And in the SPIRE bands ([Decin et al. 2007](https://ui.adsabs.harvard.edu/#abs/2007A&A...472.1041D)) are (430, 217, 105) mJy at (250, 350, 500) µm. \n",
    "\n",
    "**Requirements**: `python 3.x`, `matplotlib 2.2.2` (in `matplotlib 3.0.0` the WCS from `astropy` 3.0.4 is broken), `astropy`, `astroquery`, `photutils`.\n",
    "\n",
    "### Modification history:\n",
    "* Notebook updated to reflect backend changes as of August 2021."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Set up matplotlib and use a nicer set of plot parameters\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "from astropy.wcs import WCS\n",
    "from astropy.visualization import PercentileInterval, ImageNormalize\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.table import hstack\n",
    "from astropy.nddata import Cutout2D\n",
    "\n",
    "from photutils import SkyCircularAperture, SkyCircularAnnulus, aperture_photometry\n",
    "\n",
    "from astroquery.esasky import ESASky\n",
    "from astroquery.simbad import Simbad "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# set the download dir for ESASky products\n",
    "#\n",
    "#download_dir=\"/Users/nalvarez/Downloads\"\n",
    "#\n",
    "# define the target and get its coordinates from Simbad\n",
    "target_name = 'beta And'\n",
    "target = Simbad.query_object(target_name)\n",
    "target_coord = SkyCoord(ra=target[\"RA\"][0], dec=target[\"DEC\"][0], unit=(u.hourangle,u.deg), frame='icrs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Search in ESASky\n",
    "\n",
    "We search for _Herschel_ maps in ESASky, using as an inpu the taret name. The search radius is set to 6\" (that's one pixel of the SPIRE 250 µm map). This should be sufficient. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# We'll search ESASky for $\\beta$ And directly, using the Simbad name resolver\n",
    "maps = ESASky.query_region_maps(position=target_name, radius=\"6 arcsec\", missions=['Herschel'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check the ESASky results table\n",
    "\n",
    "The results are returned in a `astropy.Table.table`. It is useful to check it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# let's check the results table and what observations are there\n",
    "#\n",
    "print (maps[\"HERSCHEL\"].info)\n",
    "maps[\"HERSCHEL\"][\"observation_id\",\"filter\",\"duration\"].pprint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's pick up SPIRE observation, i.e. one that has 250, 350, and 500 µm\n",
    "# one example is 13422263815, this is index 3, we remove all the rest\n",
    "nher = len(maps[\"HERSCHEL\"])\n",
    "ikeep = 3\n",
    "maps[\"HERSCHEL\"].remove_rows(np.delete(range(nher),ikeep))\n",
    "maps[\"HERSCHEL\"][\"observation_id\",\"filter\",\"duration\"].pprint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we can download the map\n",
    "#maps_data = ESASky.get_maps(query_table_list=maps, download_dir=download_dir)\n",
    "maps_data = ESASky.get_maps(query_table_list=maps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# extract the HDU for the maps in a dictionary\n",
    "# there should be only one index [0]\n",
    "hdu = {}\n",
    "spire_bands = ['250','350','500']\n",
    "for band in spire_bands:\n",
    "    hdu[band] = maps_data[\"HERSCHEL\"][0][band]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualise the maps\n",
    "\n",
    "To visualise each SPIRE map we use `astropy.visualization.ImageNormalize` module, doing a colour cut at the 98% of the pixel distribution, using the `PercentileInterval`. There are many options for pixel cut ans stretch selection, so check the module documentation.\n",
    "\n",
    "It is necessary to ignore the NaNs in during the `ImageNormalization`. \n",
    "\n",
    "Alternatively one can use [APLpy](https://aplpy.github.io) for visualisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Display the three SPIRE maps, will use colour interval is 98%.\n",
    "# nan are messing the normalization, so will take care of it\n",
    "#\n",
    "fig = plt.figure(figsize=(15,5),dpi=100)\n",
    "pp = 98.0 # \n",
    "for k,band in enumerate(spire_bands):\n",
    "    wcs = WCS(hdu[band]['image'].header)\n",
    "    ax = fig.add_subplot(1,3,k+1,projection=wcs)\n",
    "    ax.set_title(f'{band} µm')\n",
    "    lon = ax.coords['ra']\n",
    "    lon.set_axislabel('RA (J2000.0)')\n",
    "    lon.set_major_formatter('hh:mm:ss.s')\n",
    "    lat = ax.coords['dec']\n",
    "    if (k == 0):\n",
    "        lat.set_axislabel('Dec (J2000.0)')\n",
    "    else:\n",
    "        lat.set_axislabel('')\n",
    "    lat.set_major_formatter('dd:mm')    \n",
    "    ximage = hdu[band]['image']\n",
    "    norm = ImageNormalize(ximage.data[~np.isnan(ximage.data)], interval=PercentileInterval(pp))\n",
    "    ax.imshow(ximage.data,norm=norm,cmap=plt.cm.gray,origin='lower',interpolation='nearest')\n",
    "    ax.grid(True)\n",
    "plt.tight_layout(pad=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Converting the input maps from MJy/sr to Jy/pixel\n",
    "\n",
    "To do aperture photmetry we need to convert the SPIRE maps to `Jy/pixel`. The pipeline produced SPIRE maps are either in units of `Jy/beam` (point-source calibrated maps) or in `MJy/sr` (extended-source calibrated maps). The maps available in `ESASky` are the extended-source calibrated ones, so they are in units of `MJy/sr`.\n",
    "\n",
    "Doing aperture photometry of a point source starting from extended-source calibrated map will require the following steps, as explained in the SPIRE Handbook:\n",
    "\n",
    "1. Convert to point-source calibration, dividing by `KPtoE`\n",
    "2. Divide by $\\Omega_\\mathrm{pipe}$ to convert from _per steradian_ to _per pixel_.\n",
    "\n",
    "After these two steps the maps should be in `Jy/pixel` and we have to set the proper unit in the header."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the map units, to make sure it's MJy/sr\n",
    "print (hdu['250']['image'].header['BUNIT'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Using the SPIRE Handbook we copy the necessary parameters\n",
    "# NOTE: we also apply the same scaling for the error image!\n",
    "#\n",
    "KPtoE = {'250': 90.646, '350': 51.181, '500': 23.580}\n",
    "# the beam solid angle in arcsec^2\n",
    "omega_pipe = {'250': 469.35, '350': 831.27, '500': 1804.31} # arcsec^2\n",
    "for band in spire_bands:\n",
    "    wcs = WCS(hdu[band]['image'].header)\n",
    "    pixscale = wcs.wcs.cdelt[1]*3600.0\n",
    "    hdu[band]['image'].data = pixscale**2*hdu[band]['image'].data/KPtoE[band]/omega_pipe[band]\n",
    "    hdu[band]['image'].header['BUNIT'] = 'Jy/pixel'\n",
    "    hdu[band]['error'].data = pixscale**2*hdu[band]['error'].data/KPtoE[band]/omega_pipe[band]\n",
    "    hdu[band]['error'].header['BUNIT'] = 'Jy/pixel'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the maps are in `Jy/pixel`. Let's define the coordinates of the target $\\beta$ Andromeda:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And define the apertures for the target, and the background. We use the default aperture sizes as described in the SPIRE Handbook, table 5.8."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aperture = dict()\n",
    "aperture['250'] = SkyCircularAperture(target_coord, r=22. * u.arcsec)\n",
    "aperture['350'] = SkyCircularAperture(target_coord, r=30. * u.arcsec)\n",
    "aperture['500'] = SkyCircularAperture(target_coord, r=40. * u.arcsec)\n",
    "# and the background annuli\n",
    "back = SkyCircularAnnulus(target_coord, r_in = 60.0*u.arcsec, r_out=90.0*u.arcsec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integrate the total flux\n",
    "\n",
    "Next is the actual measurement of the total flux within the target aperture and the background annulus. We do the background subtraction within the loop over the bands.\n",
    "\n",
    "We use the `error` extension of each SPIRE map and estimate the error on the derived flux within the input aperture.\n",
    "\n",
    "Few comments:\n",
    "* the aperture areas can only be calculated if the aperture is converted to pixels. To do this the method to_pixlel() needs a WCS.\n",
    "* The final result is the flux density and the error both `Jy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ap_flux = dict()\n",
    "ap_flux_err = dict()\n",
    "print (\"band,flux,error\")\n",
    "for band in spire_bands:\n",
    "    #photo[band] = aperture_photometry(hdu[band]['image'], aperture[band])\n",
    "    img = hdu[band]['image'].data\n",
    "    wcs = WCS(hdu[band]['image'].header)\n",
    "    err_img = hdu[band]['error'].data\n",
    "    photo = aperture_photometry(img, aperture[band],error=err_img,wcs=wcs)\n",
    "    bkg_photo = aperture_photometry(img, back,wcs=wcs)\n",
    "    wcs = WCS(hdu[band]['image'].header)\n",
    "    bkg_area = back.to_pixel(wcs).area\n",
    "    aper_area = aperture[band].to_pixel(wcs).area\n",
    "    ap_flux[band] = photo['aperture_sum'][0] - aper_area*bkg_photo['aperture_sum'][0]/bkg_area\n",
    "    ap_flux_err[band] = photo['aperture_sum_err'][0]\n",
    "    print (\"{},{:.4f},{:.4f} Jy\".format(band,ap_flux[band],ap_flux_err[band]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Corrections\n",
    "\n",
    "Fow all the subsequent corrections we need the source spectrum, as the pipeline assumes a source with $\\nu F_\\nu$ = const. Our assumption is that in the SPIRE bands $\\beta$ And has a powerlaw spectrum with an index 2 (i.e. Rayleigh-Jeans). \n",
    "\n",
    "We can use the tables in the SPIRE Handbook and pick up the beam correction (Table 5.5) and the colour-correction for a point source (Table 5.7).\n",
    "\n",
    "And finally, we apply the aperture correciton, which corrects for the point-like source flux lost outside the selected aperture. The corrections for the default apertures and sources with different spectral index, are listed in Table 5.9. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the beam correction\n",
    "kbeam = {'250': 1.0446, '350': 1.0434, '500': 1.0746} # powerlaw with alpha=2\n",
    "print (\"band,flux,error\")\n",
    "for band in spire_bands:\n",
    "    ap_flux[band] = ap_flux[band]*kbeam[band]\n",
    "    ap_flux_err[band] = ap_flux_err[band]*kbeam[band]\n",
    "    print (\"{},{:.4f},{:.4f} Jy\".format(band,ap_flux[band],ap_flux_err[band]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we need to apply the colour correction KColP (Table 5.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the colour correction\n",
    "\n",
    "kcolp = {'250': 0.9454, '350': 0.9481, '500': 0.9432} # powerlaw with alpha=2\n",
    "print (\"band,flux,error\")\n",
    "for band in spire_bands:\n",
    "    ap_flux[band] = ap_flux[band]*kcolp[band]\n",
    "    ap_flux_err[band] = ap_flux_err[band]*kcolp[band]\n",
    "    print (\"{},{:.4f},{:.4f} Jy\".format(band,ap_flux[band],ap_flux_err[band]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally the aperture correction, which is also source SED dependent (Table 5.8). Note that we use the standard apertures of (22,30,40) arcsec. For different ones the user has to use the beam profiles from the SPIRE calibration context and dervie the correction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the aperture correction\n",
    "\n",
    "kaper = {'250': 1.28, '350': 1.242, '500': 1.2610}\n",
    "print (\"band,flux,error\")\n",
    "for band in spire_bands:\n",
    "    ap_flux[band] = ap_flux[band]*kaper[band]\n",
    "    ap_flux_err[band] = ap_flux_err[band]*kaper[band]\n",
    "    print (\"{},{:.4f},{:.4f} Jy\".format(band,ap_flux[band],ap_flux_err[band]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison with models\n",
    "\n",
    "And finally let's compare with Decin et al. (2007) models:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = {'250': 430, '350': 217, '500': 105} # mJy\n",
    "print ('band,measured,error,model,model/measured')\n",
    "for band in spire_bands:\n",
    "    print (\"{},{:.2f}+/-{:.2f},{},{:.2f}\".format(band,ap_flux[band]*1000, ap_flux_err[band]*1000,model[band], model[band]/ap_flux[band]/1000.0))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualise the apertures and the background annuli just to make sure they make sense.\n",
    "\n",
    "We'll use Cutout2d to zoom in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,5),dpi=100)\n",
    "pp = 98.0 # colour cut percentage\n",
    "zoomSize = u.Quantity((4.0,4.0 ), u.arcmin)\n",
    "for k,band in enumerate(spire_bands):\n",
    "    wcs = WCS(hdu[band]['image'].header)\n",
    "    ax = fig.add_subplot(1,3,k+1,projection=wcs)\n",
    "    ax.set_title(f\"{band} µm\")\n",
    "    lon = ax.coords['ra']\n",
    "    lon.set_axislabel('RA (J2000.0)')\n",
    "    lon.set_major_formatter('hh:mm:ss.s')\n",
    "    lat = ax.coords['dec']\n",
    "    if (k == 0):\n",
    "        lat.set_axislabel('Dec (J2000.0)')\n",
    "    else:\n",
    "        lat.set_axislabel('')\n",
    "    #lat.set_major_formatter('dd:mm')    \n",
    "    ximage = hdu[band]['image']\n",
    "    norm = ImageNormalize(ximage.data[~np.isnan(ximage.data)], interval=PercentileInterval(pp))\n",
    "    cutout = Cutout2D(ximage.data, target_coord, zoomSize, wcs=wcs)\n",
    "    ax.imshow(cutout.data,norm=norm,cmap=plt.cm.gray,origin='lower',interpolation='nearest')\n",
    "    wcs_cut = cutout.wcs\n",
    "    aperture[band].to_pixel(wcs_cut).plot(color='r')\n",
    "    back.to_pixel(wcs_cut).plot(color='b')\n",
    "    ax.grid(True)\n",
    "plt.tight_layout(pad=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare with the SPIRE Point Source Catalogue\n",
    "\n",
    "The first release of the SPIRE Point Source Catalogue (SPSC) is available in `ESASky`, so we can extract the photometry in there and compare with the results from the aperture photometry.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# search ESASky for $\\beta$ And\n",
    "#\n",
    "# first, let's see what catalogues are available\n",
    "ESASky.list_catalogs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# will extract the flux, its error and signal-to-noise ratio\n",
    "#\n",
    "# note the flux density is for a source with power-law SED with index $\\alpha$ = -1.\n",
    "#\n",
    "# The SPSC is based on the timeline fitter, so it only needs the colour correction, see the Handbook, flowchart 5.21.\n",
    "#\n",
    "spsc_f = {}\n",
    "spsc_f_err = {}\n",
    "print (\"band, Flux (mJy), Flux_err (mJy), SNR\")\n",
    "for band in ['250','350','500']:\n",
    "    output = ESASky.query_object_catalogs(position=target_name, catalogs=f'Herschel-SPSC-{band}')\n",
    "    spsc_f[band] = output[0]['flux'][0]*kcolp[band]\n",
    "    spsc_f_err[band] = output[0]['flux_err'][0]*kcolp[band]\n",
    "    print (\"{}, {:.2f}, {:.2f}, {}\".format(band,spsc_f[band],spsc_f_err[band],output[0]['snr'][0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# let's see them graphically, we use the model as the baseline and do the ratio\n",
    "#\n",
    "spire_bands = ['250','350','500']\n",
    "r1 = [] # aperture to model\n",
    "r2 = [] # spsc to model\n",
    "r1_err = [] # error aperture to model\n",
    "r2_err = [] # error spsc to model\n",
    "#\n",
    "# will propagate the errors, assuming the model has zero error\n",
    "#\n",
    "for band in spire_bands:\n",
    "    r1.append(1000*ap_flux[band]/model[band])\n",
    "    r2.append(spsc_f[band]/model[band])\n",
    "    r1_err.append(1000*ap_flux_err[band]/model[band])\n",
    "    r2_err.append(spsc_f_err[band]/model[band])\n",
    "#\n",
    "#\n",
    "fig = plt.figure(figsize=(8,5),dpi=100)\n",
    "fig, ax = plt.subplots()\n",
    "ax.axhline(y=1,color='k',ls='dashed')\n",
    "ax.errorbar(spire_bands,r1,yerr=r1_err,marker='s',label='Aperture/model')\n",
    "ax.errorbar(spire_bands,r2,yerr=r2_err,marker='s',label='SPSC flux/model')\n",
    "ax.set_xlabel('SPIRE band (µm)')\n",
    "ax.set_ylabel('Ratio to model')\n",
    "ax.set_title(r'{} photmetry result'.format(target_name))\n",
    "ax.grid(True)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## You can also try:\n",
    "\n",
    "* Another target\n",
    "* Center the aperture on the peak (note, the aperture centre should be the same at 250, 350 and 500 µm).\n",
    "* Overplot all SPSC source in the field.\n",
    "* Extract also PACS 70|100 and 160 µm and build an SED for a source. Note: PACS maps are already in Jy/pixel but they suffer from correlated noise, so obtaining uncertainy on flux is tricky.\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}