{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction\n", "\n", "The Herschel/SPIRE Spectral Feature Finder (FF) is an automated tool which detects significant spectral features within SPIRE Fourier Transform Spectrometer (FTS) data. The FF is designed to facilitate data-mining and preliminary analysis of SPIRE spectral data and is not a replacement for detailed scientific analysis. A full introduction to the FF routine and products is presented in the project's [readme](http://archives.esac.esa.int/hsa/legacy/HPDP/SPIRE/SPIRE-S/spectral_feature_catalogue/SpectralFeatureFinder_releaseNote.html) page. The FF products can be accessed from this [archive](https://www.cosmos.esa.int/web/herschel/spire-spectral-feature-catalogue). In this notebook we show the structure of the FF products and provide examples of how to search FF products using typical search criteria. We also include examples for visually inspecting FF products and SPIRE data. A brief section on estimating integrated line flux for features detected by the FF is also presented. \n", "\n", "The full set of FF products can be downloaded from the link above, however, in this notebook we will only download select products when necessary. Products will be downloaded from the archive and will be placed in a chosen directory which we will assigned to the variable `data_dir`. The archive has a directory structure identical to the FF products you can download from the link above.\n", "\n", "The standard Anaconda distribution is recommended as it includes most of the packages required to run this notebook by default. However, Anaconda does not include the `lmfit` package and likely needs to be installed by the user. This notebook makes use of a custom package called `spireff` developed to aid in the acquisition, manipulation, and visualization of FF and SPIRE Herschel Science Archive (HSA) products. `spireff` is available in the `archive` and has the dependencies listed below.\n", "\n", "## Acknowledging or Citing the Feature Finder Software or Products If you plan to use this notebook\n", "for research or work publications we ask that you please cite the following Feature Finder Paper:\n", "\n", "* R Hopwood, I Valtchanov, L D Spencer, J Scott, C Benson, N Marchili, N Hładczuk, E T Polehampton, N Lu, G Makiwa, D A Naylor, B G Gom, G Noble, M J Griffin, The Herschel SPIRE Fourier Transform Spectrometer Spectral Feature Finder I. The Spectral Feature Finder and Catalogue, Monthly Notices of the Royal Astronomical Society, Volume 496, Issue 4, August 2020, Pages 4874–4893, https://doi.org/10.1093/mnras/staa1612\n", "\n", "* and provide link and refer to DOI [10.5270/esa-lysf2yi](https://doi.org/10.5270/esa-lysf2yi)\n", "\n", "This notebook was developed with Python 3.7 and the following packages:\n", "* Numpy 1.18.1\n", "* Matplotlib 3.1.3\n", "* Astropy 4.0.1\n", "* Scipy 1.4.1\n", "* Pandas 1.0.4\n", "* lmfit 0.9.13\n", "* Requests 2.23.0\n", "* Other standard packages" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "import sys\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import astropy.units as u\n", "from astropy.table import Table\n", "from astropy.coordinates import SkyCoord\n", "from scipy.interpolate import interp1d\n", "%matplotlib inline\n", "\n", "spireff_path = '../'\n", "sys.path.insert(1, spireff_path)\n", "import spireff.io as spio\n", "import spireff.figtools as tools\n", "#from spireff.gui import GUI\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "#All saved data will be stored in this directory\n", "save_dir = './savedata'\n", "if (not os.path.isdir(save_dir)):\n", " os.mkdir(save_dir)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# SPIRE Automated Feature Extraction Catalogue\n", "\n", "The Automated Feature Extraction Catalogue (SAFECAT) collates all of the features detected by the FF and their associated parameters. It can be used as a convenient resource for searching FF products. A version of SAFECAT which shifts detected features to their rest frame based on the radial velocity estimate obtained by the FF also exists. SAFECAT can be accessed with the following function." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "SAFECAT = spio.safecat(save_dir=save_dir, save=True)\n", "SAFECAT = Table.read(SAFECAT)\n", "SAFECAT.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "FF products are natively stored as FITS files. This file format is quite useful for storing general data and is easy to work with, but not ideal for visual inspection. To improve visualization, we convert SAFECAT to an Astropy table." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Searching the Catalogue\n", "\n", "In this section we demonstrate how to search SAFECAT using a select few search criteria. The searches shown generalize naturally to other SAFECAT columns show above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Single Column\n", "\n", "To search for all rows with a specific observation identification number (obsid) one can use the code below. Note this works for any valid column name, not just \"obsid\"." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "column, val = 'obsid', 1342188190\n", "ind = np.where(SAFECAT[column] == val)[0]\n", "\n", "newSAFECAT = SAFECAT.copy() #Creates a copy of the SAFECAT table that has reduced significant figures\n", "\n", "colTitles = ['frequency', 'frequencyError', 'SNR', 'ra', 'dec', 'velocity', 'velocityError', 'velocityCorrFreq', 'templateFrequency']\n", "\n", "for i in range(len(colTitles)):\n", " if colTitles[i] == 'ra' or colTitles[i] == 'dec':\n", " newSAFECAT[colTitles[i]].info.format = \"0.5f\"\n", " else:\n", " newSAFECAT[colTitles[i]].info.format = \"0.2f\"\n", "\n", "print(newSAFECAT[ind])#.show_in_notebook() *replace print() with .show_in_notebook() for high quality tables" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Multiple Columns\n", "\n", "The above observation targets an extended source. In the FF extended and semi-extended sources are processed using both their point-source and extended-source calibrated spectra. If only the extended calibration products are desired, one can use the following code. Note this extends naturally to any number of columns and values." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "column1, val1 = 'obsid', 1342188190\n", "column2, val2 = 'calibration', 'extendedSource'\n", "ind = np.where((SAFECAT[column1] == val1) & \\\n", " (SAFECAT[column2] == val2))[0]\n", "\n", "print(newSAFECAT[ind])#.show_in_notebook()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Range of Values\n", "\n", "If the user is interested in only a range of values, a range of frequency values for example, one can use the following code. Again, this extends to any number of columns and values." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "column1, val1 = 'obsid', 1342188190\n", "column2, val2 = 'calibration', 'extendedSource'\n", "column3, v_range = 'frequency', (500.0, 700.0)\n", "ind = np.where((SAFECAT[column1] == val1) & \\\n", " (SAFECAT[column2] == val2) & \\\n", " (SAFECAT[column3] >= v_range[0]) & \\\n", " (SAFECAT[column3] <= v_range[1]))[0]\n", "\n", "print(newSAFECAT[ind])#.show_in_notebook()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Coordinate Neighborhood\n", "\n", "When searching sky coordinates values, the use of the `astropy.coordinates` module is recommended. This module, along with the `astropy.wcs` module, provide convenient high-level functionality for working with sky coordinates and astronomy data.\n", "\n", "As a demonstration, we will search SAFECAT for columns within 1 arcmin of ra = 316.75 deg, dec = 42.23 deg." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ra, dec = [316.75], [42.23]\n", "radius = 1*u.arcmin\n", "coord = SkyCoord(ra, dec, unit='deg')\n", "search_coords = SkyCoord(SAFECAT['ra'], SAFECAT['dec'], unit='deg')\n", "ind, _, sep, _ = coord.search_around_sky(search_coords, radius)\n", "\n", "print('Maximum separation :', max(sep))\n", "\n", "print(newSAFECAT[ind])#.show_in_notebook()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Graphical Interface\n", "\n", "A custom graphical interface has been designed to provide a quick interactive way to search FF products and find relevant information. The interface uses `SAFECAT` to facilitate its functionality and runs on `tkinter`. The graphical interface can be used by executing the following cell." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#GUI(SAFECAT, save_dir=save_dir)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Observation Specific Products\n", "\n", "A set of observation specific products are generated for each observation processed by the FF. These products consist of a fits table for the set of detected features and associated parameters found within the spectrum, a fits table for the fitted continuum parameters, and a quick inspection postcard. The `spireff` module provides a convenience function for retrieving these products from the archive." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sparse Products\n", "\n", "In this section we focus on sparse products." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "spio.set_verbose(True)\n", "sparse_obsid = 1342270193\n", "result = spio.get_observation(sparse_obsid, useFF=True, useHsa=True, \n", " save_dir=save_dir, save=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `get_observation` function can print useful information about the products it retrieves. These print statements can be enabled or disabled (the default) using the `set_verbose` function. The print statements above show the `spec_type`, `obs_mode`, and `cal` used when retrieving products, which can be set manually in the `get_observation` function. However, if `useFF=True` and `spec_type=None` (the default) these values are overwritten by their nominal values as determined by the FF. If `useHsa=True`, the function will also return level 2 HSA (spectrometer or photometer) products. If `spec_type` is 'calHpdp' or 'calBgs', the corresponding [HPDP](http://archives.esac.esa.int/hsa/legacy/HPDP/SPIRE/SPIRE-S/cal_targets) or [BGS](http://archives.esac.esa.int/hsa/legacy/HPDP/SPIRE/SPIRE-S/BKGS/) spectra will be returned. To get the standard HSA products for these obsids, set `spec_type='spg'`. \n", "\n", "The `result` of `get_observation` is a dictionary of dictionaries containing the FF products for the input obsid. The key structure of `result` can get complicated but calling `spio._keys_summary(result)` will print a summary of the key structure of `result`. This is done automatically in the `get_observation` function and the output is shown above.\n", "\n", "Let's inspect the content of the FF products in the `result` variable. As previously mentioned, FF products are natively stored as fits files. The 'catalogue' field contains an HDUList and its `data` table contains the features catalogue." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sp_prod = result['ff']['sparse_point']\n", "\n", "sp_prod_cat = sp_prod['catalogue']\n", "print('Type :', type(sp_prod_cat), '\\n')\n", "sp_prod_cat.info(); print()\n", "\n", "sp_cat_table = Table.read(sp_prod_cat, hdu='data')\n", "#Reduce the columns to 2 significant figures\n", "sp_cat_table['frequency'].info.format = \"0.2f\"\n", "sp_cat_table['frequencyError'].info.format = \"0.4f\"\n", "sp_cat_table['SNR'].info.format = \"0.2f\"\n", "\n", "print(sp_cat_table)#.show_in_notebook()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The interesting meta data is found in the catalogue table header." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sp_prod_cat['data'].header" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The 'continuum' field contains an HDUList and its `data` table contains the fitted continuum parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sp_prod_cont = sp_prod['continuum']\n", "print('Type :', type(sp_prod_cont), '\\n')\n", "sp_prod_cont.info(); print()\n", "\n", "sp_cont_table = Table.read(sp_prod_cont, hdu='data')\n", "\n", "colTitles = ['p0', 'p0err', 'p1', 'p1err', 'p2', 'p2err', 'p3', 'p3err']\n", "\n", "for i in range(len(colTitles)):\n", " sp_cont_table[colTitles[i]].info.format = \"0.3g\"\n", "\n", "print(sp_cont_table)#.show_in_notebook()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The remaining 'postcard' field contains the quick inspection postcard for this obsid." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(dpi=150)\n", "plt.imshow(sp_prod['postcard'])\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The 'spectrum' field in `result` contains the HPDP spectrum for this observation. The data is again stored as an HDUList and it contains detector specific tables. It is worth noting that if FF products do not use HPDP or BGS spectra, the 'spectrum' variable may be either 'HR_point_pnt' or 'HR_point_ext' for point-source calibrated or extended-source calibrated SPG spectra, respectively." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hsa_prod = result['spectrum']\n", "\n", "print('Type :', type(hsa_prod), '\\n')\n", "hsa_prod.info(); print()\n", "\n", "spectrum_table = Table.read(hsa_prod, hdu='SLWC3')\n", "\n", "spectrum_table['wave'].info.format = \"0.2f\"\n", "spectrum_table['flux'].info.format = \"0.2f\"\n", "spectrum_table['error'].info.format = \"0.2f\"\n", "\n", "print(spectrum_table)#.show_in_notebook()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A convenience function called `parse_sparse` exists to parse the `result` variable of sparse observations into something more convenient." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sparse_prod = spio.parse_sparse(sp_prod['catalogue'], \n", " sp_prod['continuum'], \n", " spec_hdu=hsa_prod)\n", "spio._keys_summary(sparse_prod)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now each detector contains its own dictionary with fields for detected lines, fitted continuum parameters, spectral frequency axis ('wave'), and spectral flux. Note there is a `detectors` variable in `parse_sparse` which allows you to select a subset of detectors. The default value returns all detectors in the FF features catalogue, and the subset should only included detectors in the default list." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Off-Axis Sparse Products\n", "\n", "The sparse catalogue focuses on the central SLWC3 and SSWD4 detectors from the SPIRE detector arrays. However, a complementary catalogue exists for off-axis detectors. These products can be accessed using the argument `off_axis=True` in the `get_observation` function. The off-axis catalogue is generated by processing exclusively extended calibration spectra to give access to the full detector array including vignetted detectors." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "spio.set_verbose(False) # Trun off print statements\n", "offax_obsid = sparse_obsid # 1342191234 obsid currently missing in archive\n", "result = spio.get_observation(offax_obsid, useFF=True, useHsa=True, \n", " cal='ext', spec_type='spg', off_axis=True,\n", " save_dir=save_dir, save=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The off-axis products have the same structure as the sparse products, but include an additional table containing the radial velocity estimates for each detector. Velocity estimates are included in the meta data of FF sparse products." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "offax_prod = result['ff']['sparse_offAx']\n", "\n", "offax_prod_cat = offax_prod['catalogue']\n", "print('Type :', type(offax_prod_cat), '\\n')\n", "offax_prod_cat.info(); print()\n", "\n", "offax_velocity_table = Table.read(offax_prod_cat, hdu='velocity')\n", "\n", "offax_velocity_table['velocity'].info.format = \"0.2f\"\n", "offax_velocity_table['velocityError'].info.format = \"0.2f\"\n", "\n", "print(offax_velocity_table)#.show_in_notebook()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The off-axis postcards are also different." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(dpi=200)\n", "plt.imshow(offax_prod['postcard'])\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mapping Products\n", "\n", "FF mapping products are generated by processing the spectrum in each pixel of an observation's HSA convolved spectral cube. Spectral cubes exist for both the SLW and SSW band but the results for both cubes are combined into the same FF product. One can access FF and HSA products using the `get_observation` function as before.\n", "\n", "Note that mapping products are large and this download may take long time to finish." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "map_obsid = 1342192173\n", "result = spio.get_observation(map_obsid, useFF=True, useHsa=True,\n", " save_dir=save_dir, save=True)\n", "\n", "spio._keys_summary(result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The features catalogue is again an HDUList containing multiple tables. The `data` table contains the detected features with associated parameters. For mapping products, we project the SLW lines into an SSW equivalent grid in order to obtain line lists for the full SPIRE band before attempting to estimate radial velocities. The extra tables give radial `velocity` estimates (and associated parameters) for this SSW equivalent grid, and the number of lines (`nLines`) in each pixel of the combined grid. Since the sky coverage of the SLW and SSW cubes is not identical, the `array` table indicates which bands were used when populating the combined grid." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "map_prod = result['ff']['mapping']\n", "\n", "map_prod_cat = map_prod['catalogue']\n", "print('Type :', type(map_prod_cat), '\\n')\n", "map_prod_cat.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The mapping features catalogue table has additional columns to identify the pixel (`column`, `row`) and coordinate (`ra`, `dec`) the feature was found in, and its associated radial `velocity` estimate." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "map_cat_table = Table.read(map_prod_cat, hdu='data')\n", "map_cat_table.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The continuum table is modified in a similar way." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "map_prod_cont = map_prod['continuum']\n", "map_cont_table = Table.read(map_prod_cont, hdu=1)\n", "map_cont_table.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The mapping product postcard is shown below." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(dpi=200)\n", "plt.imshow(map_prod['postcard'])\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we consider one of the spectral cubes in the `result` variable." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cube_slw = result['HR_SLW_cube_convol']\n", "cube_ssw = result['HR_SSW_cube_convol']\n", "print('Type :', type(cube_slw), '\\n')\n", "cube_slw.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Spectra are stored in the `image` table. It is worth noting that the image is indexed using image[freq][row][col]. A convenience function in the `spireff` package exists to generate the frequency axis for spectral cubes." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "header = cube_slw['image'].header\n", "slw_spec = cube_slw['image'].data\n", "\n", "row, col = 6,5\n", "flux = slw_spec[:, row, col]\n", "freq = spio.gen_frequency(cube_slw)\n", "\n", "plt.figure(figsize=(8,4.5))\n", "plt.plot(freq, flux)\n", "plt.xlabel('Frequency [GHz]')\n", "plt.ylabel(f'Brightness\\n[{header[\"BUNIT\"]}]')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Parse Mapping\n", "\n", "The `spireff` package also has a function for parsing mapping products into a more convenient dictionary." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Search by Pixels\n", "\n", "By default the `parse_mapping` function will return outputs for all pixels in both bands of the spectral cubes. If the `pixels` argument is provided, it will only return products associated with that pixel." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# for multiple pixels\n", "# pixels = {'SLW': [[cols], [rows]]: SSW: [[cols], [rows]]}\n", "\n", "pixels = {'SLW': [col, row]}\n", "cube_hdu = {'SLW': cube_slw, 'SSW': cube_ssw}\n", "map_prod = spio.parse_mapping(map_prod_cat, map_prod_cont, \n", " pixels=pixels, cube_hdu=cube_hdu)\n", "spio._keys_summary(map_prod)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the above implementation, a dictionary is produced for each pixel which contains the main FF products and associated spectrum (if the optional `cube_hdu` argument is defined). Note the notation uses \"ary[col][row]\", indicating right ascension then declination, to be consistent with the astropy WCS conventions. In some cases, a cascading key structure may be more desirable, which is obtained by setting `flatten=False`. Note the cascading key structure is not supported by the plotting functions in the `spireff.figtools` module." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "map_prod = spio.parse_mapping(map_prod_cat, map_prod_cont, \n", " pixels=pixels, cube_hdu=cube_hdu,\n", " flatten=False)\n", "spio._keys_summary(map_prod)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If `nearest=True`, the output will contain the nearest pixel in the complementary spectral cube." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "map_prod = spio.parse_mapping(map_prod_cat, map_prod_cont, \n", " pixels=pixels, cube_hdu=cube_hdu,\n", " nearest=True)\n", "spio._keys_summary(map_prod)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Search by Coordinates\n", "\n", "It is also possible to search for pixels by coordinates (assumed to use the ICRS frame in degrees) using the `coord` argument. The `radius` argument must be an astropy Quanity and defines a circular region about the coordinates to search for pixels who's centers lie within this region." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# for multiple coords\n", "# coord = [[ras], [decs]]\n", "\n", "ra, dec = 83.841, -5.411\n", "coord = [ra, dec]\n", "\n", "map_prod = spio.parse_mapping(map_prod_cat, map_prod_cont, \n", " cube_hdu=cube_hdu, coords=coord, radius=10*u.arcsec)\n", "spio._keys_summary(map_prod)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To recover a pixels list from the outputs of `parse_mapping`, one can use the `_pixels_in_mapping_products` function." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pixels_ = spio._pixels_in_mapping_products(map_prod)\n", "print(pixels_)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Velocity Maps\n", "\n", "The FF mapping products contain a table with radial velocity estimates for each pixel of the SSW equivalent gird. We can access this table and create a 2d-numpy array from the data in the following way." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "velocity_table = map_prod_cat['velocity'].data\n", "velocity_map = np.array([list(row) for row in velocity_table])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To display these maps on sky, we require an astropy WCS object. The `spireff` module can generate these objects from the FF continuum products using the `_make_wcs` function." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wcs = spio._make_wcs(map_prod_cont)\n", "print(wcs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the velocity maps use the SSW cube griding, we only use the SSW WCS object." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "v_unit = map_prod_cat['data'].header['TUNIT11']\n", "\n", "fig = plt.figure(figsize=(6, 5.5))\n", "ax = fig.add_subplot(111, projection=wcs['SSW'])\n", "img = ax.imshow(velocity_map, origin='lower')\n", "fig.colorbar(img, label=f\"Velocity [{v_unit}]\")\n", "ax.set_xlabel('RA')\n", "ax.set_ylabel('Dec')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Visualization\n", "\n", "The `spireff` module contains a few functions to help quickly visualize FF and HSA products." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sparse Postcards\n", "\n", "The `sparse_postcard` function takes FF and HSA products as inputs and produces a plot that follows the style of the standard FF sparse postcards. Although not strictly required, it is recommended that only one set of SLW and SSW products be passed to the function at a time." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "obsid = 1342204037\n", "result = spio.get_observation(obsid, useFF=True, useHsa=True,\n", " save_dir=save_dir, save=True)\n", "cat = result['ff']['sparse_ext']['catalogue']\n", "cont = result['ff']['sparse_ext']['continuum']\n", "spec_hdu = result['HR_spectrum_ext']\n", "spec_prod = spio.parse_sparse(cat, cont, spec_hdu=spec_hdu)\n", "\n", "fig_kwargs = {'figsize': (8,4.5)}\n", "fig, ax = tools.sparse_postcard(spec_prod, fig_kwargs=fig_kwargs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Photometer Images\n", "\n", "The discontinuity between the SLW and SSW continua is due to the semi-extended nature of the source. We can get a broader view of the source by looking at its SPIRE photometer image. Note that the photometer bands are centered at 250$\\mu$m (PSW), 350$\\mu$m (PMW), and 500$\\mu$m (PLW) while the spectrometer covers a range from 191–671$\\mu$m.\n", "\n", "We can locate a complementary photometer obsid in the HSA product meta data. The keyword is typically \"POBS_XXX\" where \"XXX\" is a three digit number starting at 000 indexing a list of photometer obsids. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "phot_obsid = spec_hdu[0].header['POBS_000']\n", "result = spio.get_observation(phot_obsid, useHsa=True, \n", " photo_band='PSW', save_dir=save_dir, save=True)\n", "phot_hdu = result['psrcPSW']\n", "\n", "detectors = ['SLWC3', 'SSWD4'] # try : None, ['SLWC3', 'SSWD4']\n", "annotate = True # try : True, False\n", "fig_kwargs = {'figsize': (6.8, 6)}\n", "foot_color = {'SLW': 'limegreen', 'SSW': 'blue'}\n", "fig, ax, img, wcs = tools.sky_image(phot_hdu, foot_data=spec_hdu, \n", " annotate=annotate, detectors=detectors,\n", " foot_color=foot_color, fig_kwargs=fig_kwargs)\n", "_ = ax.legend(loc=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Virtually all elements of the plot can be controlled by additional kwargs in the `sky_image` function (see doc string), or by manipulating the function output variables." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Additional Example\n", "\n", "As an additional example of how to leverage the `spireff` package, we generate a figure with a photometer image and sparse postcard as subplots for a mapping observation. In doing so we only need specify the mapping obsid and the pixel(s) from a spectral cube we are interested in." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "map_obsid = 1342192173\n", "pixels = {'SLW': [3, 9]}\n", "\n", "# Retrieve FF and HSA products\n", "result = spio.get_observation(map_obsid, useFF=True, useHsa=True,\n", " save_dir=save_dir, save=True)\n", "map_prod_cat = result['ff']['mapping']['catalogue']\n", "map_prod_cont = result['ff']['mapping']['continuum']\n", "cube_slw = result['HR_SLW_cube_convol']\n", "cube_ssw = result['HR_SSW_cube_convol']\n", "cube_hdu = {'SLW': cube_slw, 'SSW': cube_ssw}\n", "phot_obsid = cube_slw[0].header['POBS_000']\n", "\n", "# Parse the data\n", "map_prod = spio.parse_mapping(map_prod_cat, map_prod_cont, \n", " pixels=pixels, cube_hdu=cube_hdu,\n", " nearest=True)\n", "pixels_ = spio._pixels_in_mapping_products(map_prod)\n", "\n", "# Retrieve Photometer Image\n", "result = spio.get_observation(phot_obsid, useHsa=True, \n", " photo_band='PSW', save_dir=save_dir, save=True)\n", "phot_hdu = result['psrcPSW']\n", "\n", "# Plotting\n", "fig = plt.figure(figsize=(10,4), dpi=300)\n", "gs = fig.add_gridspec(1, 5)\n", "fig, ax = tools.sparse_postcard(map_prod, fig=fig, index=gs[0,2:])\n", "fig, ax, img, wcs = tools.sky_image(phot_hdu, foot_data=[cube_slw, cube_ssw], \n", " pixels=pixels_, fig=fig, index=gs[0,:2])\n", "_ = ax.legend(loc=2)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Integrated Flux\n", "\n", "When analyzing a spectrum, one is often interested in the integrated flux of spectral features. Integrated flux of spectral features is not provided in the FF products because obtaining accurate values typically requires careful analysis of the associated spectra. With this consideration, we demonstrate below a simple approach for obtaining rough estimates of the integrated flux for features detected by the FF. \n", "\n", "The proposed procedure involves fitting detected features using a Sinc model, which closely approximates the empirical line shape of the SPIRE FTS, and extracting the integrated flux analytically from the fitted line parameters. This process is similar to the method employed by the the standard spectral line fitting analysis script used by the Herschel Interactive Processing Environment (HIPE). Typically one would include the continuum in the fitted model, but in our case, we will simply use continuum subtracted spectral data where the fitted continuum is provided by the FF. \n", "\n", "The Sinc function is given by\n", "\n", "\\begin{equation}\n", "\\label{eq:sinc}\n", "f(\\nu) = A \\frac{\\sin(\\mathrm{arg})}{\\mathrm{arg}} \\\\\n", "\\mathrm{arg} = \\frac{\\pi(\\nu - C)}{W}\n", "\\end{equation}\n", "\n", "where $A$ and $C$ define the amplitude and centroid of the spectral feature, respectively. The width of the spectral feature, $W$, corresponds the the spectral resolution of the FTS ($\\approx$ 1.18 GHz for SPIRE). With these parameters, the integrated line flux of an unresolved spectral feature is\n", "\n", "\\begin{equation}\n", "\\label{eq:integratedFlux}\n", "\\int_{-\\infty}^{\\infty} f(\\nu) \\; d\\nu = AW\n", "\\end{equation}\n", "\n", "This analytic expression is used to calculate integrated line flux in the code below.\n", "\n", "Note: The Sinc function expresses decaying side-lobes which produce non-local effects in the spectrum. For this reasons, fitting a single spectral feature is not recommended. The best most physically real results are obtained when simultaneously fitting all detected features in the spectrum (ie, either the entire SLW or SSW band). We will perform our analysis on the already loaded mapping products from above." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model_prod = tools.model_spectrum(map_prod, vary_width=False)\n", "colors = {'SLW': 'red', 'SSW': 'blue'}\n", "flux_unit = cube_slw[1].header['BUNIT']\n", "freq_unit = cube_slw[1].header['CUNIT3']\n", "\n", "fig = plt.figure(figsize=(6,4), dpi=400)\n", "for key, val in model_prod.items():\n", " if key == 'header':\n", " continue\n", " ary = 'SLW' if 'SLW' in key else 'SSW'\n", " flux_sub = val['flux'] - val['cont']\n", " plt.plot(val['wave'], flux_sub, color=colors[ary], linewidth=4.0, label=key)\n", " (axm,) = plt.plot(val['wave'], val['model'].best_fit, color='orange')\n", " (axr,) = plt.plot(val['wave'], val['model'].residual, color='green')\n", " if ary == 'SSW':\n", " axm.set_label('Model')\n", " axr.set_label('Residual')\n", "plt.xlabel(f'Frequency [{freq_unit}]')\n", "plt.ylabel(f'Brightness [{flux_unit}]')\n", "plt.legend(loc=0)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the above example we have assumed the spectral features are not resolved. When analyzing partially or fully resolved spectral features, best results are obtained when modeling the spectral features as Sinc convolved Gaussians. However, for a quick estimate of integrated flux, one can simply set the `vary_width` keyword variable in the `model_spectrum` function to `True`. The large majority of spectral features in the SPIRE FTS dataset are unresolved and the above script should be sufficient in most cases. \n", "\n", "Below we calculate the integrated flux of spectral features using the fitted model parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for key, val in model_prod.items():\n", " if key == 'header':\n", " continue\n", " fit_pars = val['model'].params\n", " print(f'\\n{key:-^65}')\n", " print('Frequency [GHz] Flux [W m-2 sr-1] Flux Error [W m-2 sr-1]')\n", " print('-'*65)\n", " for comp in val['model'].components:\n", " pref = comp.prefix\n", " amp = fit_pars[f'{pref}amp']\n", " freq = fit_pars[f'{pref}cent']\n", " width = fit_pars[f'{pref}width']\n", " flux = amp.value*width.value*1e9\n", " fluxErr = flux*np.sqrt((amp.stderr/amp.value)**2 + (width.stderr/width.value)**2)\n", " print(f\"{freq.value:10.2f} {flux:21.3e} {fluxErr:25.3e}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Unit Conversion for Integrated Flux\n", "\n", "It is often convenient to express integrated flux in alternative units. Here we summarize how to convert from standard SPIRE FTS integrated line flux units to a few other commonly encountered units. For a more comprehensive overview, please see \"SPIRE FTS: converting integrated line fluxes to alternative units\"." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Point Source Calibration\n", "\n", "Since point source calibrated spectra in SPIRE are presented in units of Janskys, employing the analytical equation for integrated line flux returns units of [Jy GHz]. To obtain integrated line flux in SI units use\n", "\n", "\\begin{equation}\n", "\\label{eq:flux_Wm2}\n", "F_{\\mathrm{line}} = AW \\times 10^{-17} \\qquad [\\mathrm{W \\; m^{-2}}]\n", "\\end{equation}\n", "\n", "where $A$ is in units of Jy and $W$ is in units of GHz. All the following unit conversions for point source calibrated spectra will be done with respect to $F_{\\mathrm{line}}$. \n", "\n", "#### $[\\mathrm{W \\; m^{-2}}] \\rightarrow [\\mathrm{erg \\; cm^{-2} \\; s^{-1}}]$ \n", "\n", "This conversion requires only a multiplicative factor,\n", "\n", "\\begin{equation}\n", "\\label{eq:flux_ergcm2s}\n", "F' = F_{\\mathrm{line}} \\times 10^3 \\qquad [\\mathrm{erg \\; cm^{-2} \\; s^{-1}}]\n", "\\end{equation}\n", "\n", "#### $[\\mathrm{W \\; m^{-2}}] \\rightarrow [\\mathrm{Jy \\; km \\; s^{-1}}]$\n", "\n", "This conversion requires a multiplicative factor and the rest frame emission frequency, $\\nu_0$, of the line in question using units of GHz.\n", "\n", "\\begin{equation}\n", "\\label{eq:flux_jykms1}\n", "F' = \\frac{F_{\\mathrm{line}} \\; (9.418 \\times 10^{19})}{\\nu_0} \\qquad [\\mathrm{Jy \\; km \\; s^{-1}}]\n", "\\end{equation}\n", "\n", "Alternatively, one can use the rest frame emission wavelength, $\\lambda_0$, of the line in question using units of $\\mu\\mathrm{m}$.\n", "\n", "\\begin{equation}\n", "\\label{eq:flux_jykms2}\n", "F' = F_{\\mathrm{line}} \\; \\lambda_0 \\; (\\pi \\times 10^{7}) \\qquad [\\mathrm{Jy \\; km \\; s^{-1}}]\n", "\\end{equation}\n", "\n", "#### $[\\mathrm{W \\; m^{-2}}] \\rightarrow [\\mathrm{K \\; km \\; s^{-1}}]$\n", "\n", "The temperature unit Kelvins in the desired units is in reference to the so called \"brightness temperature\" and is a useful quantity for FIR/submm and radio astronomy. This conversion is somewhat more complicated in that it requires knowledge of the frequency dependent SPIRE FTS beam width. We use $\\theta_b(\\nu)$, in units of arcseconds, to denote the frequency dependent Full Width at Half Maximum (FWHM) of the beam. For SPIRE, this empirically derived function can be obtained from the following [link](http://archives.esac.esa.int/hsa/legacy/ADP/PSF/SPIRE/SPIRE-S/), but can also be accessed using the `frequency_beam` function in the `spireff.io` library. Both the SPIRE SLW and SSW detector arrays have their own functions. The conversion is then\n", "\n", "\\begin{equation}\n", "\\label{eq:flux_kkms}\n", "F' = \\frac{F_{\\mathrm{line}} \\; (1.149 \\times 10^{29})}{\\theta_b^2(\\nu) \\; \\nu_0^3} \\qquad [\\mathrm{K \\; km \\; s^{-1}}]\n", "\\end{equation}\n", "\n", "Note that $\\nu_0$ is again the rest frame emission frequency of the line in question in units of GHz, and $\\theta_b(\\nu)$ is the beam FWHM at the detected frequency of the line in question. Since this conversion is nontrivial, we provide a sample conversion below." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "freq_det = 482.1 # feature detection frequency [GHz]\n", "freq_0 = 500.0 # rest frame frequency for detected feature [GHz]\n", "amp = 1 # amplitude of detected feature [Jy]\n", "width = 1.18448225 # width of sinc function [GHz]\n", "\n", "slw_beam = spio.frequency_beam(save_dir=save_dir, save=True)['SLW']\n", "\n", "# beam FWHM at detection frequency\n", "interp = interp1d(slw_beam['frequency'], slw_beam['FWHM'])\n", "theta = interp(freq_det)\n", "\n", "# calculated integrated line flux [W m-2]\n", "F_line = amp*width*1e-17\n", "# conversion to [K km s-1]\n", "F_ = F_line*1.149e29/theta**2/freq_0**3\n", "\n", "print(f\"Integrated line flux:\\n F_line = {F_line:.2e} \"\n", " f\"[W m-2]\\n F_ = {F_:.2f} [K km s-1]\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Extended Source Calibration\n", "\n", "Since extended source calibrated spectra in SPIRE are presented in units of [W m-2 Hz-1 sr-1], employing the analytical equation for integrated line flux returns units of [W GHz m-2 Hz-1 sr-1]. To obtain integrated line flux in SI units use\n", "\n", "\\begin{equation}\n", "\\label{eq:flux_Wm2sr1}\n", "F_{\\mathrm{line}} = AW \\times 10^{9} \\qquad [\\mathrm{W \\; m^{-2} \\; sr^{-1}}]\n", "\\end{equation}\n", "\n", "where $A$ is the amplitude of the feature in units of [W m-2 Hz-1 sr-1], and $W$ is the width of the feature in units of GHz.\n", "\n", "#### $[\\mathrm{W \\; m^{-2} \\; sr^{-1}}] \\rightarrow [K]$\n", "\n", "The desired temperature units of Kelvins is again in relation to the so called \"brightness temperature\". For extended source calibrated spectra, the conversion requires simply a multiplicative factor and the emission rest frequency, $\\nu_0$, of the spectral feature in question in units of GHz. The conversion is given by\n", "\n", "\\begin{equation}\n", "\\label{eq:flux_K}\n", "F' = \\frac{F_{\\mathrm{line}} (1.04 \\times 10^{12})}{\\nu_0^2} \\qquad [K]\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Conclusion\n", "\n", "The Herschel/SPIRE Spectral Feature Finder is an automated tool which attempts to detect significant spectral features within SPIRE FTS data. This tool and the resulting dataset is designed to help with preliminary analysis and data mining of SPIRE spectral data and is not intended to replaced detailed analysis. In this document we have demonstrated how to retrieve FF products through a Python interface and how to query these products for features meeting specific search criteria. We also outlined simple methods for displaying FF results including plotting FF features and fitted continua over the SPIRE spectra obtained from the HSA, and displaying velocity maps for mapping observations. Finally, we provided sample code for how to extract integrated line flux for detected spectral features and outlined conversion methods for obtaining the integrated flux in other common units." ] }, { "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.9.2" } }, "nbformat": 4, "nbformat_minor": 4 }