Abell 1656: the Coma Cluster of Galaxies __ Simple Spectral Access Protocol (SSA)#

Stefania Amodeo¹, Thomas Boch¹, Caroline Bot¹, Giulia Iafrate², Katharina A. Lutz¹, Manon Marchand¹, Massimo Ramella², and Jenny G.Sorce³⁴⁵

  1. Université de Strasbourg, CNRS, Observatoire Astronomique de Strasbourg, UMR 7550, F-67000, Strasbourg, France

  2. INAF - Osservatorio Astronomico di Trieste

  3. Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRIStAL, F-59000 Lille, France

  4. Université Paris-Saclay, CNRS, Institut d’Astrophysique Spatiale, 91405, Orsay, France

  5. Leibniz-Institut für Astrophysik (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany

The original form of this tutorial authored by Massimo Ramella & Giulia Iafrate can be found on the EURO-VO tutorials page. The version here is an adaptation into a jupyter notebook by the Strasbourg astronomical Data Center (CDS) team.


Introduction#

The Coma Cluster is an ensemble of over 1000 galaxies that takes its name from the constellation Coma Berenices.

The goals of this notebook tutorial are to:

  1. examine the Coma cluster of galaxies (Abell 1656) using services and data from the virtual observatory within a jupyter notebook. This allows a quick evaluation of the mean redshift and velocity dispersion of the cluster. Both measurements are important to study the evolution of galaxy clusters. To do so we use redshifts and photometry (Petrosian r magnitude) from the Sloan Digital Sky Survey (SDSS) and then add redshifts from the Cluster And Infall Region Nearby Survey (CAIRNS) (Rines et al. 2003). This improves the completeness of the redshift sample,

  2. research the Mikulski Archive for Space Telescopes (MAST) for Hubble Space Telescope (HST) spectra in the region of the Coma cluster,

  3. download a spectrum from MAST and do a quick analysis of the redshift of the emission lines in the spectrum.

# Astronomy tools
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.io import fits
from astropy.table import vstack

# Access astronomical databases
import pyvo
from astroquery.simbad import Simbad
from astroquery.vizier import Vizier
from astroquery.xmatch import XMatch

# Sky visualization
from ipyaladin import Aladin

# For plots
import matplotlib.pyplot as plt

# Data handling
import numpy as np
import specutils.analysis as spec_ana
from scipy.optimize import curve_fit
from specutils import SpectralRegion, Spectrum1D

Step #1: Exploration of the Coma cluster#

Display the region of Abell 1656 in Aladin lite#

Aladin

We start by displaying the Coma Cluster in an Aladin Lite widget. We will display DSS2 color images (survey='P/DSS2/color') centered on Abell 1656 (target='A1656') and set the field of view to \(0.7^\circ\) (fov=0.7). At this distance, this field of view corresponds to approximately 1Mpc: this is large enough to look at the cluster.

If you are using Jupyter lab you can open a second python3 notebook. On this new notebook click the “Python 3” button in the top-right corner of this notebook to switch kernel. A new window will pop up and you can then select this notebook (i.e. “Abel1656…”) as a kernel. This will link the two notebooks such that they see the same variables. You may use the second notebook to look at Aladin Lite widget while doing the analysis in this notebook.

aladin = Aladin(target="A1656", fov=0.7, survey="P/DSS2/color")
aladin

You can already see galaxies composing the cluster, and a few stars, like for example the bright HD112887 star (click for a surprise Simbad query). But there is more! As with any Aladin Lite implementation, you can interact with this widget.

  • to zoom in and out, place your mouse pointer on the image and scroll,

  • you can right-click on a point to grab it and center the reticule on it,

  • with layers you can select other image surveys,

  • if you would like to look at another target, you can use the search field search to get there.

We can also interact with the variable aladin. If, for example, after zooming in and out, you wanted to set the FoV (field of view) again to \(0.7^\circ\), run:

aladin.fov = 0.4

and look back up … or call the variable aladin again to pop a new widget window. You’ll remark that the field of view and reticule position are synchronized, but not the survey if you change if manually with layers. It can allow the exploration of the same region in two different wavelengths.

aladin

Load the SDSS catalog and select galaxies#

Vizier

In this section, we will access the SDSS DR12 catalog from the VizieR catalog service. We will download all entries that are located within 40arcmin of the center of the cluster A1656. We start by querying VizieR for all catalogs that match the term SDSS DR12.

Since this query takes a few seconds, we will store its results in the variable catalog_list_sdss for further analysis.

catalog_list_sdss = Vizier.find_catalogs("SDSS DR12")
catalog_list_sdss
OrderedDict([('V/147', </>),
             ('J/ApJ/807/178', </>),
             ('J/ApJ/835/161', </>),
             ('J/ApJ/888/85', </>),
             ('J/ApJ/901/93', </>),
             ('J/ApJS/225/23', </>),
             ('J/ApJS/228/9', </>),
             ('J/ApJS/229/39', </>),
             ('J/A+A/583/A86', </>),
             ('J/A+A/596/A14', </>),
             ('J/PASP/130/H4203', </>),
             ('J/MNRAS/452/4153', </>),
             ('J/MNRAS/455/3413', </>),
             ('J/MNRAS/469/2102', </>)])

The output is not easily readable. Let’s inspect a bit more this object.

print(f"The catalog is returned and stored in a {type(catalog_list_sdss)}")
The catalog is returned and stored in a <class 'collections.OrderedDict'>

To know how to deal with this ordered dictionary, we look at its first entry and list the methods applicable to the associated value.

# just selecting the value of the first item in the dictionnary
first_entry_value = catalog_list_sdss[next(iter(catalog_list_sdss))]
# the dir() function is built in python > 3. It lists all available methods for an object
# it is a really useful tool when discovering a new python library
print(dir(first_entry_value))
['ID', '_ID', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_add_coosys', '_add_definitions', '_add_group', '_add_info', '_add_link', '_add_param', '_add_resource', '_add_table', '_add_timesys', '_add_unknown_tag', '_attr_list', '_config', '_coordinate_systems', '_description', '_element_name', '_extra_attributes', '_groups', '_ignore_add', '_infos', '_links', '_name', '_params', '_pos', '_resources', '_tables', '_time_systems', '_type', '_utype', '_utype_in_v1_2', 'coordinate_systems', 'description', 'extra_attributes', 'groups', 'infos', 'iter_coosys', 'iter_fields_and_params', 'iter_info', 'iter_tables', 'iter_timesys', 'links', 'name', 'params', 'parse', 'resources', 'tables', 'time_systems', 'to_xml', 'type', 'utype']

Among the possible methods, we chose to print the description, but feel free to explore the catalog list in any other way :)

# Let's display a list of the catalogs' names and descriptions
for key, catalog in catalog_list_sdss.items():
    print(key, ": ", catalog.description)
V/147 :  The SDSS Photometric Catalogue, Release 12 (Alam+, 2015)
J/ApJ/807/178 :  Newly rich galaxy clusters identified in SDSS-DR12 (Wen+, 2015)
J/ApJ/835/161 :  A cosmic void catalog of SDSS DR12 BOSS galaxies (Mao+, 2017)
J/ApJ/888/85 :  Ghostly strong Lya absorbers in SDSS DR12 (Fathivavsari, 2020)
J/ApJ/901/93 :  Model atmosphere analysis of hot WDs from SDSS DR12 (Bedard+, 2020)
J/ApJS/225/23 :  Compact groups of galaxies from SDSS-DR12 (MLCG) (Sohn+, 2016)
J/ApJS/228/9 :  Physical parameters of ~300000 SDSS-DR12 QSOs (Kozlowski, 2017)
J/ApJS/229/39 :  Narrow line Seyfert 1 galaxies from SDSS-DR12 (Rakshit+, 2017)
J/A+A/583/A86 :  DB white dwarfs in SDSS DR10 and DR12 (Koester+, 2015)
J/A+A/596/A14 :  Group catalogues of the local universe (Saulder+, 2016)
J/PASP/130/H4203 :  Newly spectroscopically confirmed DB white dwarfs (Kong+, 2018)
J/MNRAS/452/4153 :  SDSS DR12 QSOs [OIII] doublet (Albareti+, 2015)
J/MNRAS/455/3413 :  New white dwarf and subdwarf stars in SDSS DR12 (Kepler+, 2016)
J/MNRAS/469/2102 :  White dwarf population from the SDSS DR12 (Anguiano+, 2017)

In this list, the first item is the complete SDSS catalog, while the others are sub-catalogs produced from the full one that often adds information. These other catalogs are published in astronomical journals (ex: ApJ stands for the Astrophysical Journal) and are curated, formated, and included in the Vizier catalogs list by a team of documentalists and astronomes in the CDS.

We are interested in data from the main SDSS DR12 photometric catalog: “The SDSS Photometric Catalogue, Release 12 (Alam+, 2015)”. We query its content in 10arcmin around the center of the Coma Cluster.

results_test_sdss = Vizier.query_region("A1656", radius="0d10m0s", catalog="V/147")
print(results_test_sdss)
TableList with 1 tables:
	'0:V/147/sdss12' with 23 column(s) and 50 row(s) 

As you see, there is only one table in our resulting list of tables. Moreover, this table has only 50 rows, an extremely small number! It is actually due to the default limit in the number of rows of the Vizier.query_region function. For now, we will use this small table to see what columns could be interesting. We will query without a row limit once we know what we need. This is a general good practice and allows to avoid unneccesary manipulation of huge tables.

To access the unique table in the list, we select the first table of the catalogue with results_test_sdss[0]. And we can again print the list of available methods for this object:

print(dir(results_test_sdss[0]))
['Column', 'ColumnClass', 'MaskedColumn', 'Row', 'TableColumns', 'TableFormatter', '__array__', '__bytes__', '__class__', '__copy__', '__deepcopy__', '__delattr__', '__delitem__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__ior__', '__le__', '__len__', '__lt__', '__module__', '__ne__', '__new__', '__or__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setitem__', '__setstate__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_base_repr_', '_check_names_dtype', '_column_class', '_convert_col_for_table', '_convert_data_to_col', '_convert_string_dtype', '_copy_indices', '_first_colname', '_get_col_cls_for_table', '_init_from_cols', '_init_from_dict', '_init_from_list', '_init_from_list_of_dicts', '_init_from_ndarray', '_init_indices', '_ipython_key_completions_', '_is_list_or_tuple_of_str', '_is_mixin_for_table', '_make_index_row_display_table', '_make_table_from_cols', '_mask', '_masked', '_meta', '_new_from_slice', '_replace_cols', '_replace_column_warnings', '_repr_html_', '_rows_equal', '_set_col_parent_table_and_mask', '_set_column_attribute', '_set_masked', '_set_of_names_in_colnames', '_set_row', 'add_column', 'add_columns', 'add_index', 'add_row', 'argsort', 'as_array', 'colnames', 'columns', 'convert_bytestring_to_unicode', 'convert_unicode_to_bytestring', 'copy', 'dtype', 'field', 'filled', 'formatter', 'from_pandas', 'group_by', 'groups', 'has_masked_columns', 'has_masked_values', 'has_mixin_columns', 'iloc', 'index_column', 'index_mode', 'indices', 'info', 'insert_row', 'items', 'itercols', 'iterrows', 'keep_columns', 'keys', 'loc', 'loc_indices', 'mask', 'masked', 'meta', 'more', 'pformat', 'pformat_all', 'pprint', 'pprint_all', 'pprint_exclude_names', 'pprint_include_names', 'primary_key', 'read', 'remove_column', 'remove_columns', 'remove_indices', 'remove_row', 'remove_rows', 'rename_column', 'rename_columns', 'replace_column', 'reverse', 'round', 'show_in_browser', 'show_in_notebook', 'sort', 'to_pandas', 'update', 'values', 'values_equal', 'write']

Let’s see its columns:

print(results_test_sdss[0].colnames)
['RA_ICRS', 'DE_ICRS', 'mode', 'q_mode', 'class', 'SDSS12', 'm_SDSS12', '_tab1_18', 'Q', 'umag', 'e_umag', 'gmag', 'e_gmag', 'rmag', 'e_rmag', 'imag', 'e_imag', 'zmag', 'e_zmag', 'zsp', 'zph', 'e_zph', '__zph_']

But since the names don’t mean a lot by themselves, let’s also look at another method, the info one that gives more detailed information about each column. If these descriptions are still not clear, you can also have a look directly at the SDSS DR12 webpage.

print(results_test_sdss[0].info)
<Table length=50>
  name    dtype  unit  format                                                                   description                                                                  n_bad
-------- ------- ---- -------- --------------------------------------------------------------------------------------------------------------------------------------------- -----
 RA_ICRS float64  deg {:10.6f}                                                                                  Right Ascension of the object (ICRS) at Epoch="ObsDate" (ra)     0
 DE_ICRS float64  deg {:10.6f}                                                                                     Declination of the object (ICRS) at Epoch="ObsDate" (dec)     0
    mode   uint8                                                                                  [1/2] 1: primary (469,053,874 sources), 2: secondary (324,960,094 sources)     0
  q_mode    str1                                                                               [+] '+' indicates clean photometry (310,190,812 sources with mode 1+) (clean)     0
   class   uint8                                                                                                                       Type of object (3=galaxy, 6=star) (1)     0
  SDSS12   str19                                                                                                                     SDSS-DR12 name, based on J2000 position     0
m_SDSS12    str1                                                                         [*] The asterisk indicates that 2 different SDSS objects share the same SDSS12 name     0
_tab1_18 float64   yr  {:9.4f}                                                                                                                     Mean Observation date (6)     0
       Q   uint8                                                                                                 [1/3] Quality of the observation: 1=bad 2=acceptable 3=good     0
    umag float32  mag  {:6.3f}                                                                                         [4/38]? Model magnitude in u filter, AB scale (u) (5)     0
  e_umag float32  mag  {:6.3f}                                                                                                                  ? Mean error on umag (err_u)     0
    gmag float32  mag  {:6.3f}                                                                                         [5/40]? Model magnitude in g filter, AB scale (g) (5)     0
  e_gmag float32  mag  {:6.3f}                                                                                                                  ? Mean error on gmag (err_g)     0
    rmag float32  mag  {:6.3f}                                                                                         [4/39]? Model magnitude in r filter, AB scale (r) (5)     0
  e_rmag float32  mag  {:6.3f}                                                                                                                  ? Mean error on rmag (err_r)     0
    imag float32  mag  {:6.3f}                                                                                         [3/40]? Model magnitude in i filter, AB scale (i) (5)     0
  e_imag float32  mag  {:6.3f}                                                                                                                  ? Mean error on imag (err_i)     0
    zmag float32  mag  {:6.3f}                                                                                         [3/38]? Model magnitude in z filter, AB scale (z) (5)     0
  e_zmag float32  mag  {:6.3f}                                                                                                                  ? Mean error on zmag (err_z)     0
     zsp float64       {:8.5f}                                                                                      [-0.02/7.1]? Spectroscopic redshift (when SpObjID>0) (7)    47
     zph float64      {:10.4f}                                          [-9999/]? Photometric redshift; estimated by robust fit to nearest neighbors in a reference set (12)    17
   e_zph float64      {:10.4f}                                                                             [-9999/]? Estimated error of the photometric redshift (zErr) (12)    17
  __zph_ float32       {:7.4f} [0.009/0.9]? average redshift of the nearest neighbors; if significantly different from zph this might be a better estimate than zph (nnAvgZ)    17

For our current goal, we will need the coordinates (RA_ICRS and DE_ICRS), r-band magnitudes (rmag),and redshifts (zsp) of the galaxies in the Coma cluster. Hence, we only keep these columns.

SDSS furthermore provides information on the type of source (galaxies correspond to class = 3, see the info output above) so we’ll want to keep these.

Additionally, and this information is specific to the SDSS survey (see their website), some objects are observed more than once because of the overlaps between plates. To avoid these duplicates, we select only the primary observation of each object (corresponding to mode = 1).

A first way of obtaining unique records (mode = 1) for the galaxies (class = 3) would be to apply masks on the precedent test table:

masked_sdss = results_test_sdss[0]
masked_sdss = masked_sdss[(masked_sdss["mode"] == 1) & (masked_sdss["class"] == 3)]
# note the important double [] here
masked_sdss[["RA_ICRS", "DE_ICRS", "rmag", "zsp"]].pprint()
 RA_ICRS    DE_ICRS    rmag    zsp   
   deg        deg      mag           
---------- ---------- ------ --------
194.898322  27.824928 25.666       --
194.866758  27.836941 18.788       --
194.880537  27.835338 21.533       --
194.882785  27.830918 20.872       --
194.893763  27.823545 19.719       --
194.893826  27.829963 23.613       --
194.868513  27.834005 22.168       --
194.883249  27.839735 22.061       --
194.872526  27.850157 14.436  0.02264
194.860237  27.856882 17.548       --
194.864856  27.839634 21.783       --
194.849198  27.844535 19.528       --
194.861678  27.848732 20.482       --
194.838340  27.853454 22.103       --
194.850056  27.862254 22.681       --
194.850795  27.841235 22.446       --
194.862376  27.853645 22.085       --
194.860047  27.859902 23.805       --
194.843658  27.851190 25.908       --
194.831179  27.861593 24.187       --
194.834057  27.885992 15.178  0.02162
194.821502  27.877073 20.935       --
194.828192  27.878889 21.441       --
194.832127  27.880659 21.086       --
194.818875  27.873814 22.804       --
194.821035  27.872558 22.729       --
194.826019  27.882583 22.282       --
194.839031  27.866565 21.038  0.51578
194.837992  27.862871 22.196       --
194.836993  27.865032 24.383       --
194.838772  27.871596 18.756       --
194.844106  27.868562 21.546       --
194.844697  27.869371 22.301       --

But here is a second option using a specialized instance of the Vizier class in which we will lift the restriction on the number of rows.

To know how to do so, always use the magic help in jupyter notebooks:

help(Vizier)
Help on VizierClass in module astroquery.vizier.core object:

class VizierClass(astroquery.query.BaseQuery)
 |  VizierClass(columns=['*'], column_filters={}, catalog=None, keywords=None, ucd='', timeout=60, vizier_server='vizier.u-strasbg.fr', row_limit=50)
 |  
 |  Method resolution order:
 |      VizierClass
 |      astroquery.query.BaseQuery
 |      builtins.object
 |  
 |  Methods defined here:
 |  
 |  __init__(self, columns=['*'], column_filters={}, catalog=None, keywords=None, ucd='', timeout=60, vizier_server='vizier.u-strasbg.fr', row_limit=50)
 |      Parameters
 |      ----------
 |      columns : list
 |          List of strings
 |      column_filters : dict
 |      catalog : str or None
 |      keywords : str or None
 |      ucd : string
 |          "Unified Content Description" column descriptions.  Specifying
 |          these will select only catalogs that have columns matching the
 |          column descriptions defined on the Vizier web pages.
 |          See http://vizier.u-strasbg.fr/vizier/vizHelp/1.htx#ucd and
 |          http://cds.u-strasbg.fr/w/doc/UCD/
 |      vizier_server : string
 |          Name of the VizieR mirror to use.
 |          (This parameter's default is set from a configuration object.)
 |      timeout : number
 |          timeout for connecting to server
 |          (This parameter's default is set from a configuration object.)
 |      row_limit : int
 |          Maximum number of rows that will be fetched from the result
 |          (set to -1 for unlimited).
 |          (This parameter's default is set from a configuration object.)
 |  
 |  find_catalogs(self, keywords, include_obsolete=False, verbose=False, max_catalogs=None, return_type='votable')
 |      Search Vizier for catalogs based on a set of keywords, e.g. author name
 |      
 |      Parameters
 |      ----------
 |      keywords : list or string
 |          List of keywords, or space-separated set of keywords.
 |          From `Vizier <http://vizier.u-strasbg.fr/doc/asu-summary.htx>`_:
 |          "names or words of title of catalog. The words are and'ed, i.e.
 |          only the catalogues characterized by all the words are selected."
 |      include_obsolete : bool, optional
 |          If set to True, catalogs marked obsolete will also be returned.
 |      max_catalogs : int or None
 |          The maximum number of catalogs to return.  If ``None``, all
 |          catalogs will be returned.
 |      
 |      Returns
 |      -------
 |      resource_dict : dict
 |          Dictionary of the "Resource" name and the VOTable resource object.
 |          "Resources" are generally publications; one publication may contain
 |          many tables.
 |      
 |      Examples
 |      --------
 |      >>> from astroquery.vizier import Vizier
 |      >>> catalog_list = Vizier.find_catalogs('Kang W51')
 |      >>> print(catalog_list)
 |      {u'J/ApJ/706/83': <astropy.io.votable.tree.Resource at 0x108d4d490>,
 |       u'J/ApJS/191/232': <astropy.io.votable.tree.Resource at 0x108d50490>}
 |      >>> print({k:v.description for k,v in catalog_list.items()})
 |      {u'J/ApJ/706/83': u'Embedded YSO candidates in W51 (Kang+, 2009)',
 |       u'J/ApJS/191/232': u'CO survey of W51 molecular cloud (Bieging+, 2010)'}
 |  
 |  get_catalogs(self, *args, **kwargs)
 |      Queries the service and returns a table object.
 |      
 |      Query the Vizier service for a specific catalog
 |      
 |      Parameters
 |      ----------
 |      catalog : str, Resource, or list, optional
 |          The catalog(s) that will be retrieved
 |      
 |      Returns
 |      -------
 |      table : A `~astropy.table.Table` object.
 |  
 |  get_catalogs_async(self, catalog, verbose=False, return_type='votable', get_query_payload=False)
 |      Query the Vizier service for a specific catalog
 |      
 |      Parameters
 |      ----------
 |      catalog : str, Resource, or list, optional
 |          The catalog(s) that will be retrieved
 |      
 |      Returns
 |      -------
 |      response : `~requests.Response`
 |          Returned if asynchronous method used
 |  
 |  query_constraints(self, *args, **kwargs)
 |      Queries the service and returns a table object.
 |      
 |      Send a query to Vizier in which you specify constraints with
 |      keyword/value pairs.
 |      
 |      See `the vizier constraints page
 |      <http://vizier.cfa.harvard.edu/vizier/vizHelp/cst.htx>`_ for details.
 |      
 |      Parameters
 |      ----------
 |      catalog : str or list, optional
 |          The catalog(s) which must be searched for this identifier.
 |          If not specified, all matching catalogs will be searched.
 |      kwargs : dict
 |          Any key/value pairs besides "catalog" will be parsed
 |          as additional column filters.
 |      
 |      Examples
 |      --------
 |      >>> from astroquery.vizier import Vizier
 |      >>> # note that glon/glat constraints here *must* be floats
 |      >>> result = Vizier.query_constraints(catalog='J/ApJ/723/492/table1',
 |      ...                                   GLON='>49.0 & <51.0', GLAT='<0')
 |      >>> result[result.keys()[0]].pprint()
 |          GRSMC      GLON   GLAT   VLSR  ... RD09 _RA.icrs _DE.icrs
 |      ------------- ------ ------ ------ ... ---- -------- --------
 |      G049.49-00.41  49.49  -0.41  56.90 ... RD09   290.95    14.50
 |      G049.39-00.26  49.39  -0.26  50.94 ... RD09   290.77    14.48
 |      G049.44-00.06  49.44  -0.06  62.00 ... RD09   290.61    14.62
 |      G049.04-00.31  49.04  -0.31  66.25 ... RD09   290.64    14.15
 |      G049.74-00.56  49.74  -0.56  67.95 ... RD09   291.21    14.65
 |      G050.39-00.41  50.39  -0.41  41.17 ... RD09   291.39    15.29
 |      G050.24-00.61  50.24  -0.61  41.17 ... RD09   291.50    15.06
 |      G050.94-00.61  50.94  -0.61  40.32 ... RD09   291.85    15.68
 |      G049.99-00.16  49.99  -0.16  46.27 ... RD09   290.97    15.06
 |      G049.44-00.06  49.44  -0.06  46.27 ... RD09   290.61    14.62
 |      G049.54-00.01  49.54  -0.01  56.05 ... RD09   290.61    14.73
 |      G049.74-00.01  49.74  -0.01  48.39 ... RD09   290.71    14.91
 |      G049.54-00.91  49.54  -0.91  43.29 ... RD09   291.43    14.31
 |      G049.04-00.46  49.04  -0.46  58.60 ... RD09   290.78    14.08
 |      G049.09-00.06  49.09  -0.06  46.69 ... RD09   290.44    14.31
 |      G050.84-00.11  50.84  -0.11  50.52 ... RD09   291.34    15.83
 |      G050.89-00.11  50.89  -0.11  59.45 ... RD09   291.37    15.87
 |      G050.44-00.41  50.44  -0.41  64.12 ... RD09   291.42    15.34
 |      G050.84-00.76  50.84  -0.76  61.15 ... RD09   291.94    15.52
 |      G050.29-00.46  50.29  -0.46  14.81 ... RD09   291.39    15.18
 |      
 |      Returns
 |      -------
 |      table : A `~astropy.table.Table` object.
 |  
 |  query_constraints_async(self, catalog=None, return_type='votable', cache=True, get_query_payload=False, **kwargs)
 |      Send a query to Vizier in which you specify constraints with
 |      keyword/value pairs.
 |      
 |      See `the vizier constraints page
 |      <http://vizier.cfa.harvard.edu/vizier/vizHelp/cst.htx>`_ for details.
 |      
 |      Parameters
 |      ----------
 |      catalog : str or list, optional
 |          The catalog(s) which must be searched for this identifier.
 |          If not specified, all matching catalogs will be searched.
 |      kwargs : dict
 |          Any key/value pairs besides "catalog" will be parsed
 |          as additional column filters.
 |      
 |      Returns
 |      -------
 |      response : `requests.Response`
 |          The response of the HTTP request.
 |      
 |      Examples
 |      --------
 |      >>> from astroquery.vizier import Vizier
 |      >>> # note that glon/glat constraints here *must* be floats
 |      >>> result = Vizier.query_constraints(catalog='J/ApJ/723/492/table1',
 |      ...                                   GLON='>49.0 & <51.0', GLAT='<0')
 |      >>> result[result.keys()[0]].pprint()
 |          GRSMC      GLON   GLAT   VLSR  ... RD09 _RA.icrs _DE.icrs
 |      ------------- ------ ------ ------ ... ---- -------- --------
 |      G049.49-00.41  49.49  -0.41  56.90 ... RD09   290.95    14.50
 |      G049.39-00.26  49.39  -0.26  50.94 ... RD09   290.77    14.48
 |      G049.44-00.06  49.44  -0.06  62.00 ... RD09   290.61    14.62
 |      G049.04-00.31  49.04  -0.31  66.25 ... RD09   290.64    14.15
 |      G049.74-00.56  49.74  -0.56  67.95 ... RD09   291.21    14.65
 |      G050.39-00.41  50.39  -0.41  41.17 ... RD09   291.39    15.29
 |      G050.24-00.61  50.24  -0.61  41.17 ... RD09   291.50    15.06
 |      G050.94-00.61  50.94  -0.61  40.32 ... RD09   291.85    15.68
 |      G049.99-00.16  49.99  -0.16  46.27 ... RD09   290.97    15.06
 |      G049.44-00.06  49.44  -0.06  46.27 ... RD09   290.61    14.62
 |      G049.54-00.01  49.54  -0.01  56.05 ... RD09   290.61    14.73
 |      G049.74-00.01  49.74  -0.01  48.39 ... RD09   290.71    14.91
 |      G049.54-00.91  49.54  -0.91  43.29 ... RD09   291.43    14.31
 |      G049.04-00.46  49.04  -0.46  58.60 ... RD09   290.78    14.08
 |      G049.09-00.06  49.09  -0.06  46.69 ... RD09   290.44    14.31
 |      G050.84-00.11  50.84  -0.11  50.52 ... RD09   291.34    15.83
 |      G050.89-00.11  50.89  -0.11  59.45 ... RD09   291.37    15.87
 |      G050.44-00.41  50.44  -0.41  64.12 ... RD09   291.42    15.34
 |      G050.84-00.76  50.84  -0.76  61.15 ... RD09   291.94    15.52
 |      G050.29-00.46  50.29  -0.46  14.81 ... RD09   291.39    15.18
 |  
 |  query_object(self, *args, **kwargs)
 |      Queries the service and returns a table object.
 |      
 |      Serves the same purpose as `query_object` but only
 |      returns the HTTP response rather than the parsed result.
 |      
 |      Parameters
 |      ----------
 |      object_name : str
 |          The name of the identifier.
 |      catalog : str or list, optional
 |          The catalog(s) which must be searched for this identifier.
 |          If not specified, all matching catalogs will be searched.
 |      radius : `~astropy.units.Quantity` or None
 |          A degree-equivalent radius (optional).
 |      coordinate_system : str or None
 |          If the object name is given as a coordinate, you *should* use
 |          `~astroquery.vizier.VizierClass.query_region`, but you can
 |          specify a coordinate frame here instead (today, J2000, B1975,
 |          B1950, B1900, B1875, B1855, Galactic, Supergal., Ecl.J2000, )
 |      
 |      
 |      Returns
 |      -------
 |      table : A `~astropy.table.Table` object.
 |  
 |  query_object_async(self, object_name, catalog=None, radius=None, coordinate_frame=None, get_query_payload=False, return_type='votable', cache=True)
 |      Serves the same purpose as `query_object` but only
 |      returns the HTTP response rather than the parsed result.
 |      
 |      Parameters
 |      ----------
 |      object_name : str
 |          The name of the identifier.
 |      catalog : str or list, optional
 |          The catalog(s) which must be searched for this identifier.
 |          If not specified, all matching catalogs will be searched.
 |      radius : `~astropy.units.Quantity` or None
 |          A degree-equivalent radius (optional).
 |      coordinate_system : str or None
 |          If the object name is given as a coordinate, you *should* use
 |          `~astroquery.vizier.VizierClass.query_region`, but you can
 |          specify a coordinate frame here instead (today, J2000, B1975,
 |          B1950, B1900, B1875, B1855, Galactic, Supergal., Ecl.J2000, )
 |      
 |      Returns
 |      -------
 |      response : `~requests.Response`
 |          The response of the HTTP request.
 |  
 |  query_region(self, *args, **kwargs)
 |      Queries the service and returns a table object.
 |      
 |      Serves the same purpose as `query_region` but only
 |      returns the HTTP response rather than the parsed result.
 |      
 |      Parameters
 |      ----------
 |      coordinates : str, `astropy.coordinates` object, or `~astropy.table.Table`
 |          The target around which to search. It may be specified as a
 |          string in which case it is resolved using online services or as
 |          the appropriate `astropy.coordinates` object. ICRS coordinates
 |          may also be entered as a string.  If a table is used, each of
 |          its rows will be queried, as long as it contains two columns
 |          named ``_RAJ2000`` and ``_DEJ2000`` with proper angular units.
 |      radius : convertible to `~astropy.coordinates.Angle`
 |          The radius of the circular region to query.
 |      inner_radius : convertible to `~astropy.coordinates.Angle`
 |          When set in addition to ``radius``, the queried region becomes
 |          annular, with outer radius ``radius`` and inner radius
 |          ``inner_radius``.
 |      width : convertible to `~astropy.coordinates.Angle`
 |          The width of the square region to query.
 |      height : convertible to `~astropy.coordinates.Angle`
 |          When set in addition to ``width``, the queried region becomes
 |          rectangular, with the specified ``width`` and ``height``.
 |      catalog : str or list, optional
 |          The catalog(s) which must be searched for this identifier.
 |          If not specified, all matching catalogs will be searched.
 |      column_filters: dict, optional
 |          Constraints on columns of the result. The dictionary contains
 |          the column name as keys, and the constraints as values.
 |      frame : str, optional
 |          The frame to use for the request. It should be 'fk5', 'icrs',
 |          or 'galactic'. This choice influences the the orientation of
 |          box requests.
 |      
 |      
 |      Returns
 |      -------
 |      table : A `~astropy.table.Table` object.
 |  
 |  query_region_async(self, coordinates, radius=None, inner_radius=None, width=None, height=None, catalog=None, get_query_payload=False, cache=True, return_type='votable', column_filters={}, frame='fk5')
 |      Serves the same purpose as `query_region` but only
 |      returns the HTTP response rather than the parsed result.
 |      
 |      Parameters
 |      ----------
 |      coordinates : str, `astropy.coordinates` object, or `~astropy.table.Table`
 |          The target around which to search. It may be specified as a
 |          string in which case it is resolved using online services or as
 |          the appropriate `astropy.coordinates` object. ICRS coordinates
 |          may also be entered as a string.  If a table is used, each of
 |          its rows will be queried, as long as it contains two columns
 |          named ``_RAJ2000`` and ``_DEJ2000`` with proper angular units.
 |      radius : convertible to `~astropy.coordinates.Angle`
 |          The radius of the circular region to query.
 |      inner_radius : convertible to `~astropy.coordinates.Angle`
 |          When set in addition to ``radius``, the queried region becomes
 |          annular, with outer radius ``radius`` and inner radius
 |          ``inner_radius``.
 |      width : convertible to `~astropy.coordinates.Angle`
 |          The width of the square region to query.
 |      height : convertible to `~astropy.coordinates.Angle`
 |          When set in addition to ``width``, the queried region becomes
 |          rectangular, with the specified ``width`` and ``height``.
 |      catalog : str or list, optional
 |          The catalog(s) which must be searched for this identifier.
 |          If not specified, all matching catalogs will be searched.
 |      column_filters: dict, optional
 |          Constraints on columns of the result. The dictionary contains
 |          the column name as keys, and the constraints as values.
 |      frame : str, optional
 |          The frame to use for the request. It should be 'fk5', 'icrs',
 |          or 'galactic'. This choice influences the the orientation of
 |          box requests.
 |      
 |      Returns
 |      -------
 |      response : `requests.Response`
 |          The response of the HTTP request.
 |  
 |  ----------------------------------------------------------------------
 |  Readonly properties defined here:
 |  
 |  valid_keywords
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |  
 |  catalog
 |      The default catalog to search.  If left empty, will search all
 |      catalogs.
 |  
 |  column_filters
 |      Filters to run on the individual columns. See the Vizier website
 |      for details.
 |  
 |  columns
 |      Columns to include.  The special keyword 'all' will return ALL
 |      columns from ALL retrieved tables.
 |  
 |  keywords
 |      The set of keywords to filter the Vizier search
 |  
 |  ucd
 |      UCD criteria: see http://vizier.u-strasbg.fr/vizier/vizHelp/1.htx#ucd
 |      
 |      Examples
 |      --------
 |      >>> Vizier.ucd = '(spect.dopplerVeloc*|phys.veloc*)'
 |  
 |  ----------------------------------------------------------------------
 |  Data and other attributes defined here:
 |  
 |  __abstractmethods__ = frozenset()
 |  
 |  ----------------------------------------------------------------------
 |  Methods inherited from astroquery.query.BaseQuery:
 |  
 |  __call__(self, *args, **kwargs)
 |      init a fresh copy of self
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from astroquery.query.BaseQuery:
 |  
 |  __dict__
 |      dictionary for instance variables (if defined)
 |  
 |  __weakref__
 |      list of weak references to the object (if defined)

Where we can directly read that we should give the columns names as a list, the column_filters as a dictionary, and the row_limit as an int. It is even specified that row_limit = -1 puts no limit on the number of rows. So let us do this:

vizier_instance = Vizier(
    columns=["RA_ICRS", "DE_ICRS", "rmag", "zsp"],
    column_filters={"class": "==3", "mode": "==1"},
    row_limit=-1,
)

Now that we have prepared everything, we can query VizieR for all galaxies, which are primary sources and within 40arcmin of the center of the Coma cluster. Note that VizieR knows that A1656 is the center of the Coma cluster. If you wanted to search at some other place in the sky, you could also give coordinates (in the form of Astropy’s SkyCoord) instead of 'A1656'.

results_coma_sdss = vizier_instance.query_region(
    "A1656",
    radius="0d40m0s",
    catalog="V/147",
)
print(results_coma_sdss)
TableList with 1 tables:
	'0:V/147/sdss12' with 4 column(s) and 23770 row(s) 

As you can see, the result of our query includes data from one catalog for 23770 objects (i.e. the table has 23770 rows). The output from this query is again a list of Astropy Table objects. We assign the table to the variable sdss and work with this from here on.

sdss = results_coma_sdss[0]
sdss
Table length=23770
RA_ICRSDE_ICRSrmagzsp
degdegmag
float64float64float32float64
195.20810327.36213622.487--
195.21148827.36813221.056--
195.23487027.36629220.244--
195.21753227.36635422.695--
195.22089427.37538122.440--
195.22125327.37226222.266--
195.24037127.36685121.962--
195.24266227.36937622.027--
195.22176627.36266224.560--
195.23426127.38115521.036--
195.24207927.37813420.934--
195.24120927.37705720.794--
195.23818627.37735721.092--
195.24513327.38984119.907--
195.25197627.40214521.357--
195.25695227.39001721.258--
195.26546027.38060121.522--
195.24643127.39735022.648--
195.25009727.40102522.637--
195.25392627.39940823.011--
195.26249127.38563522.951--
195.28657727.38824820.959--
195.28856127.39200519.016--
195.29003727.39643621.412--
195.28867227.39527622.572--
195.29154127.38537121.636--
195.29656427.38835121.005--
195.30202727.39381622.420--
195.30089527.40290123.235--
195.30931827.39658923.431--
195.31769627.41903822.683--
195.33056627.40838323.092--
195.32945227.40903721.482--
195.31354327.40823921.673--
195.32008627.41069522.400--
............
194.70800928.58121120.669--
194.71008028.57981423.138--
194.71391728.58665420.766--
194.68811928.58111122.428--
194.69110028.57699722.265--
194.69352828.58678223.271--
194.69384128.60403918.491--
194.69732028.59363123.191--
194.69185028.59734922.257--
194.69814128.59632820.575--
194.69914128.60326820.506--
194.70227828.60264221.307--
194.71025128.59274121.945--
194.70795328.59293821.489--
194.71425628.59025122.083--
194.69660828.58983722.693--
194.64757528.59030722.930--
194.65536628.58996123.437--
194.67001228.58470721.504--
194.64238128.58283520.433--
194.65221628.58435420.563--
194.66166928.59055721.955--
194.66505928.57847421.662--
194.64483528.58047923.652--
194.66722628.59231821.988--
194.66643628.59384823.527--
194.67490228.58984121.980--
194.71671928.61199020.690--
194.72502128.61468219.639--
194.72894428.61149221.789--
194.73944128.61139221.153--
194.70883028.60383322.000--
194.72613828.59585122.754--
194.72273628.60698323.570--
194.74290728.61141522.426--

By navigating this table, you’ll remark that very few entries have a spectroscopic redshift zsp. Thus we will complete the results later with information from another table. But we’ll first ensure the validity of the entries here.

Identify the brightest sources as stars contaminating the sample#

aladin's logo We have already restricted our sample to sources that are classified as galaxies (`class` = 3). However, for very bright sources, stars might be confused with galaxies. To test this and exclude any contamination from stars, we now take a closer look at the brightest sources. To do so we select sources brighter than 11.5 mag in r-band (`rmag < 11.5`) and check with the Aladin widget what these sources look like.
aladin
limit_magnitude = 11.5
stars = sdss[sdss["rmag"] < limit_magnitude]
aladin.add_table(stars["RA_ICRS", "DE_ICRS"])
print(
    f"Our sample contains {len(stars)} really bright sources with a magnitude below {limit_magnitude}.",
)
aladin.fov = 9
Our sample contains 11 really bright sources with a magnitude below 11.5.

With aladin.add_table(stars), we have added symbols to the Aladin Lite widget at the location of the brightest sources (i.e. stars). You will be able to find all the brightest sources. By looking at each source you will find that these are indeed stars. Feel free to change the limit magnitude or other parameters to see if you can find better filters, this one is pretty rudimentary.

You can choose to hide or show the symbols for the table entries in the layer options, the tableappears as catalog

Build a subset of galaxies with photometry and redshift in SDSS#

Now on to exploring the galaxies in our sample. We again build a subset, this time for all sources fainter than 11.5mag (to leave out the stars identified in the section above) but brighter than 17.77mag, which is the completeness limit of the SDSS spectroscopic sample.

sdss_sample = sdss[(sdss["rmag"] > 11.5) & (sdss["rmag"] < 17.77)]
sdss_sample
Table length=632
RA_ICRSDE_ICRSrmagzsp
degdegmag
float64float64float32float64
195.27580027.39782615.8420.02706
194.90563127.33604716.0110.02326
195.12170627.33321616.0151.33693
195.07422627.38751615.3670.03678
195.07495427.38596617.653--
195.21991527.38436817.561--
195.20801027.40578213.9260.02198
195.05383227.39651317.5530.29843
195.04666827.46008417.6250.02282
195.04602827.48077617.3380.15689
195.01688327.51268817.2390.12721
195.13978427.50412416.1170.01861
195.11189727.51559115.0340.02595
195.12180427.51484215.6621.32939
195.07959327.55372114.6780.01964
195.10491527.55226016.4430.02424
195.14864827.57422515.1290.01700
195.15253127.57584517.521--
194.58554427.42938316.9250.02516
194.56361427.46469717.736--
194.50313927.45398816.6400.02562
194.50649527.48966614.4160.02543
194.46148227.49086417.4160.02452
194.47651127.49063314.9980.01661
194.48206227.49009917.555--
194.43013227.57759416.2570.02411
194.43726827.57612017.760--
194.41420827.60607716.7000.07007
194.78583527.36713617.3770.13948
194.80710227.40256914.3840.01886
194.72412327.38863617.521--
194.83459027.43523717.2980.02567
194.85245127.51733817.2120.01965
194.68955127.44695417.474--
194.65762427.46395415.8400.02091
............
195.09232728.40072616.3780.09117
195.12231228.45564317.3800.02412
195.19894128.40614617.2850.06196
195.20240428.41431716.8360.14366
195.21150428.41994017.7350.22162
195.18279528.41633316.1030.02093
195.17193328.45500817.4740.14457
195.17049028.51194317.620--
195.17815528.51042316.563--
195.01706328.41408915.433--
195.00764128.43481616.2090.02699
195.05924128.45806617.7470.02610
195.16950628.51989315.869--
195.09458928.57447515.8500.02159
195.08081228.59127517.737--
195.31231928.52175517.1170.02813
194.92210628.50719515.0100.02748
194.89493728.53728217.397--
194.89863028.55136514.9600.02530
194.95276228.59333317.666--
194.94571128.62198817.545--
195.01829528.60348215.4490.02274
195.06527628.60489316.6380.17906
195.08483628.63667416.414--
195.04358328.62524617.696--
194.97372128.63273914.688--
194.97757228.63338515.845--
194.98097128.62917216.499--
194.91720128.63080215.0660.01784
194.91688728.63382517.413--
194.82469928.62400116.5500.03470
194.34604928.36914517.7290.06847
194.50832728.45636116.7380.16127
194.61924728.52682515.4380.06149
194.56276828.52178316.3580.02304

Improve the completeness with other sources of redshifts in Vizier#

Vizier

As you can see in the table above, not all galaxies in our zsp17 sample have redshift measurements (some rows have ‘–’ in the ‘zsp’ column, i.e. they are masked). So to improve the completeness of our sample we will now use Vizier to search for redshifts in the Rines et al (2003) catalog. First, find all catalogs that match the search terms ‘redshifts Rines 2003’:

catalog_list_rines = Vizier.find_catalogs("Rines 2003")
for k, v in catalog_list_rines.items():
    print(k, ": ", v.description)
J/ApJ/700/331 :  Light curves of type Ia supernovae (CfA3) (Hicken+, 2009)
J/ApJ/757/22 :  Strong and weak lensing analysis of A2261 (Coe+, 2012)
J/ApJ/767/15 :  Hectospec Cluster Survey (HeCS) (Rines+, 2013)
J/ApJ/783/52 :  Redshifts in the field of A383 (Geller+, 2014)
J/ApJ/797/106 :  Redshifts in nine galaxy cluster fields (Hwang+, 2014)
J/ApJ/819/63 :  Hectospec survey of SZ clusters (HeCS-SZ) (Rines+, 2016)
J/ApJ/855/100 :  The HectoMAP cluster survey. II. X-ray clusters (Sohn+, 2018)
J/ApJ/856/172 :  The HectoMAP cluster survey. I. (Sohn+, 2018)
J/ApJ/862/172 :  HeCS-red: Hectospec surveys of redMaPPer clusters (Rines+, 2018)
J/ApJ/871/129 :  A redshift catalog of the galaxy cluster A2029 (Sohn+, 2019)
J/ApJ/891/129 :  The HeCS-omnibus catalog: SDSS & MMT sp. obs. (Sohn+, 2020)
J/ApJS/229/20 :  MMT/Hectospec redshift survey for Abell 2029 (Sohn+, 2017)
J/AJ/120/2338 :  Abell 576 redshifts (Rines+, 2000)
J/AJ/124/1266 :  Redshift survey around Abell 2199 (Rines+, 2002)
J/AJ/126/2152 :  Cluster And Infall Region Nearby Survey. I (Rines+, 2003)
J/AJ/128/1078 :  Cluster and Infall Region Nearby Survey. II (Rines+, 2004)
J/AJ/131/527 :  UBVRI light curves of 44 type Ia supernovae (Jha+, 2006)
J/AJ/132/1275 :  CIRS (Cluster Infall Regions in the SDSS). I. (Rines+, 2006)
J/AJ/135/1598 :  Optical spectroscopy of type Ia supernovae (Matheson+, 2008)
J/AJ/135/1837 :  Spectroscopy in A2199 and Virgo clusters (Rines+, 2008)

Among these, we find the ‘J/AJ/126/2152’ catalog: Cluster And Infall Region Nearby Survey. I (Rines+, 2003): the one we’d like to read.

So again, let’s take a quick look at this catalog.

results_test_rines = Vizier.query_region(
    "A1656",
    radius="0d10m0s",
    catalog="J/AJ/126/2152",
)
print(results_test_rines)
TableList with 2 tables:
	'0:J/AJ/126/2152/clusters' with 10 column(s) and 1 row(s) 
	'1:J/AJ/126/2152/galaxies' with 6 column(s) and 50 row(s) 
results_test_rines[0]
Table length=1
Clustern_ClusterRAJ2000DEJ2000czsigmap_3s_sigmap_ca_LXTXR
km / skm / skm / s1e+36 WkeV
str5str1str10str9int32int16int16float32float32uint8
A1656g12 59 31.9+27 54 106973104295718.08.02
results_test_rines[1]
Table length=50
RAJ2000DEJ2000cze_czr_czCluster
km / skm / s
str11str11int32int32uint8str5
12 59 03.85+27 57 32.6697862A1656
12 59 05.83+27 59 49.57699271A1656
12 59 09.43+28 02 27.07220312A1656
12 59 11.51+28 00 33.16942302A1656
12 59 13.00+27 58 39.06747641A1656
12 59 13.83+28 04 36.07805281A1656
12 59 14.36+27 55 58.258669752A1656
12 59 14.61+27 53 43.96450212A1656
12 59 15.20+27 58 14.04849241A1656
12 59 15.61+27 59 48.91208292A1656
12 59 16.00+27 58 14.081993672A1656
12 59 16.22+27 59 10.747487752A1656
12 59 16.68+27 58 56.97195752A1656
12 59 17.47+27 57 16.947607752A1656
12 59 17.97+28 05 27.982317592A1656
12 59 19.71+27 58 24.47140692A1656
12 59 19.90+28 05 03.04642221A1656
12 59 20.16+28 04 27.27195752A1656
12 59 20.20+27 53 09.0645292A1656
12 59 21.30+27 52 17.485471752A1656
12 59 21.40+27 58 24.764111002A1656
12 59 21.94+28 06 17.025782752A1656
12 59 22.82+27 53 49.0508782A1656
12 59 23.20+27 54 43.06899311A1656
12 59 23.21+27 58 07.874768752A1656
12 59 23.24+27 58 54.98874752A1656
12 59 23.41+27 55 10.06824482A1656
12 59 23.83+27 53 49.01838882A1656
12 59 25.02+27 59 48.15507372A1656
12 59 25.31+27 58 04.47683452A1656
12 59 25.34+27 56 04.27604442A1656
12 59 25.69+27 58 32.16082452A1656
12 59 26.40+27 59 54.3668462A1656
12 59 26.44+27 51 24.65007752A1656
12 59 27.74+27 50 12.169256772A1656
12 59 27.76+28 03 17.182309642A1656
12 59 28.40+28 05 06.636501002A1656
12 59 28.80+28 02 25.05569112A1656
12 59 29.20+27 51 03.06814401A1656
12 59 29.21+27 56 30.85613582A1656
12 59 30.00+27 57 22.06800271A1656
12 59 30.10+28 05 44.0434352A1656
12 59 30.27+28 01 14.671661732A1656
12 59 30.50+27 56 24.0123282A1656
12 59 30.80+27 53 03.04741251A1656
12 59 31.32+28 02 49.46935211A1656
12 59 31.57+28 06 02.17420362A1656
12 59 31.87+27 51 40.74640342A1656
12 59 32.05+27 55 15.268151072A1656
12 59 32.32+27 51 02.28184902A1656
results_test_rines[1].info()
<Table length=50>
  name  dtype  unit                        description                      
------- ----- ------ -------------------------------------------------------
RAJ2000 str11                                Hour of right ascension (J2000)
DEJ2000 str11                                  Degree of declination (J2000)
     cz int32 km / s                                                Redshift
   e_cz int32 km / s ? Error in cz [NULL integer written as an empty string]
   r_cz uint8                                   [1/2] Redshift reference (1)
Cluster  str5                                  Abell cluster name, or "Both"

Xmatch

After inspecting the result of the test query, we see that the first table describes the cluster as an ensemble. The second one describes individual galaxies in the cluster. The later, named 1:J/AJ/126/2152/galaxies contains the information we want.

To see which of the galaxies in the 1:J/AJ/126/2152/galaxies table could fill in the gaps in our SDSS table, we first isolate the galaxies without redshifts in zsp17. Then we spatially crossmatch the two tables using the CDS XMatch service via the astroquery.XMatch.query module.

# mask zsp17 for entries with a nan value in the zsp column
sdss_sample_without_zsp = sdss_sample[np.isnan(sdss_sample["zsp"])]
sdss_sample_with_zsp = sdss_sample[~np.isnan(sdss_sample["zsp"])]
sdss_sample_without_zsp
Table length=143
RA_ICRSDE_ICRSrmagzsp
degdegmag
float64float64float32float64
195.07495427.38596617.653--
195.21991527.38436817.561--
195.15253127.57584517.521--
194.56361427.46469717.736--
194.48206227.49009917.555--
194.43726827.57612017.760--
194.72412327.38863617.521--
194.68955127.44695417.474--
194.64944427.46614815.715--
194.66766727.47290915.452--
195.01113127.71133717.675--
194.92457027.72366817.388--
194.62757527.56440117.750--
194.50398427.55912416.508--
194.86023727.85688217.548--
194.73490027.80356817.533--
194.81493227.89915217.438--
194.73442727.90813017.753--
195.53511427.63695917.496--
195.40198027.62465616.385--
195.21377327.68915617.643--
195.45929527.62571615.971--
195.46688327.63619417.703--
195.45894727.67597416.417--
195.37507127.85272616.690--
195.65864528.11366216.557--
195.05384527.71887217.576--
195.11840227.76122117.495--
195.21289127.77646017.229--
195.26623827.76248217.335--
195.24344027.76120917.208--
195.25996627.77531917.583--
195.26253527.78811617.240--
195.25952227.79012717.686--
195.26315227.78456617.198--
............
194.86957228.08577916.586--
194.83587328.08127016.341--
194.91736028.17934717.323--
194.98150828.24319017.449--
194.88037028.21731417.509--
194.88841228.23894711.939--
194.89432428.24145315.136--
194.90534928.23776116.443--
194.89899228.23741716.831--
194.89352728.22980516.541--
194.88604828.24316915.754--
194.87308528.24002616.263--
194.97571428.40250615.953--
194.97635328.40249115.579--
194.63815928.32487514.218--
194.64502928.32144314.610--
194.63591728.31751316.208--
194.65798128.31576417.002--
194.52155128.24997117.705--
194.62513828.42663117.051--
194.75143128.37425617.215--
195.17049028.51194317.620--
195.17815528.51042316.563--
195.01706328.41408915.433--
195.16950628.51989315.869--
195.08081228.59127517.737--
194.89493728.53728217.397--
194.95276228.59333317.666--
194.94571128.62198817.545--
195.08483628.63667416.414--
195.04358328.62524617.696--
194.97372128.63273914.688--
194.97757228.63338515.845--
194.98097128.62917216.499--
194.91688728.63382517.413--
xmatch_sdss_rines = XMatch.query(
    cat1=sdss_sample_without_zsp,
    cat2="vizier:J/AJ/126/2152/galaxies",
    max_distance=5 * u.arcsec,
    colRA1="RA_ICRS",
    colDec1="DE_ICRS",
)
xmatch_sdss_rines
Table length=25
angDistRA_ICRSDE_ICRSrmagzsp_RAJ2000_DEJ2000RAJ2000DEJ2000cze_czr_czCluster
float64float64float64float64int64float64float64str11str11int64int64int64str5
4.756129195.07495427.38596617.653--195.074333327.387166713 00 17.84+27 23 13.81111492A1656
0.268487194.56361427.46469717.736--194.563666727.464638912 58 15.28+27 27 52.77625272A1656
0.089868194.72412327.38863617.521--194.72412527.388611112 58 53.79+27 23 19.07654502A1656
0.284257194.86023727.85688217.548--194.860166727.856833312 59 26.44+27 51 24.65007752A1656
1.346011194.734927.80356817.533--194.734916727.803194412 58 56.38+27 48 11.57864642A1656
0.271796195.21377327.68915617.643--195.2137527.689083313 00 51.30+27 41 20.78184752A1656
0.993647195.46688327.63619417.703--195.4667527.635944413 01 52.02+27 38 09.417437272A1656
0.503127195.11840227.76122117.495--195.11837527.761083313 00 28.41+27 45 39.96835752A1656
2.31869195.44246127.957815.158--195.44187527.957416713 01 46.05+27 57 26.74706132A1656
0.745055195.24513628.02685117.642--195.244916728.026777813 00 58.78+28 01 36.442190522A1656
1.981451195.21475528.04287313.603--195.214166728.043055613 00 51.40+28 02 35.08748231A1656
0.822759194.92607527.85886716.959--194.92587527.858722212 59 42.21+27 51 31.41691182A1656
0.650358194.89986827.90611217.639--194.899791727.905944412 59 35.95+27 54 21.49054752A1656
2.615351195.07366927.95526213.772--195.072916727.955555613 00 17.50+27 57 20.06899231A1656
0.649338195.04483227.99094314.771--195.0447527.990777813 00 10.74+27 59 26.823078412A1656
2.674574194.84741427.91160413.451--194.846666727.911944412 59 23.20+27 54 43.06899311A1656
3.582934195.16057528.0160616.203--195.160791728.015083313 00 38.59+28 00 54.37594321A1656
0.815348195.17734128.11630517.018--195.177208328.116111113 00 42.53+28 06 58.06206752A1656
0.843958195.08209928.12152117.696--195.081833328.121527813 00 19.64+28 07 17.56659652A1656
4.3169195.54976928.17322316.363--195.548791728.172388913 02 11.71+28 10 20.68950112A1656
0.546321194.62947927.89550317.676--194.629333327.895583312 58 31.04+27 53 44.16804832A1656
0.518468194.56269328.12590714.632--194.562833328.125833312 58 15.08+28 07 33.07428231A1656
4.509648194.86957228.08577916.586--194.868333328.085166712 59 28.40+28 05 06.636501002A1656
2.489322194.9173628.17934717.323--194.916583328.1792512 59 39.98+28 10 45.35437--2A1656
1.210618195.16950628.51989315.869--195.16912528.519861113 00 40.59+28 31 11.5890142A1656

Build the final catalog including the Rines et al. (2003) redshifts#

The resulting table of the cross-match above contains 25 rows, so we have found recession velocity (‘cz’) measurements for 25 galaxies. Now let’s add these data to the zsp17_with table.

# converts redshift into velocity
c = 2.998e5  # km/s speed of light
sdss_sample_with_zsp["cz"] = sdss_sample_with_zsp["zsp"] * c
# from the cross math result, only keep columns that are in sdss_sample_with_zsp
columns = ["RA_ICRS", "DE_ICRS", "rmag", "cz"]
# put together the two tables
complete_sample = vstack(
    [sdss_sample_with_zsp[columns], xmatch_sdss_rines[columns]],
    metadata_conflicts="silent",
)
# don't forget to set the unit of the newly created column
complete_sample["cz"].unit = u.km / u.s
complete_sample
Table length=514
RA_ICRSDE_ICRSrmagcz
degdegmagkm / s
float64float64float64float64
195.27580027.39782615.8428112.588000000001
194.90563127.33604716.0116973.348
195.12170627.33321616.015400811.614
195.07422627.38751615.36711026.644
195.20801027.40578213.9266589.604
195.05383227.39651317.55389469.314
195.04666827.46008417.6256841.436
195.04602827.48077617.33847035.622
195.01688327.51268817.23938137.558
195.13978427.50412416.1175579.278
195.11189727.51559115.0347779.81
195.12180427.51484215.662398551.12200000003
195.07959327.55372114.6785888.072
195.10491527.55226016.4437267.152
195.14864827.57422515.1295096.6
194.58554427.42938316.9257542.968
194.50313927.45398816.6407680.876
194.50649527.48966614.4167623.914000000001
194.46148227.49086417.4167351.0960000000005
194.47651127.49063314.9984979.678
194.43013227.57759416.2577228.178
194.41420827.60607716.70021006.985999999997
194.78583527.36713617.37741816.104
194.80710227.40256914.3845654.227999999999
194.83459027.43523717.2987695.865999999999
194.85245127.51733817.2125891.070000000001
194.65762427.46395415.8406268.818
194.63360127.45633415.0597021.316
194.77631827.49363117.6678163.554
194.82687227.51351116.1106040.97
194.75677327.53692915.576934530.5639999999
194.93845127.45104817.32836284.794
194.91587127.57646815.7465009.657999999999
195.04337927.59499815.5855594.268
195.08384227.61570617.44643468.002
............
194.92210628.50719515.0108238.504
194.89863028.55136514.9607584.94
195.01829528.60348215.4496817.452
195.06527628.60489316.63853682.188
194.91720128.63080215.0665348.432000000001
194.82469928.62400116.55010403.060000000001
194.34604928.36914517.72920527.306
194.50832728.45636116.73848348.746
194.61924728.52682515.43818434.702
194.56276828.52178316.3586907.392000000001
195.07495427.38596617.65311114.0
194.56361427.46469717.7367625.0
194.72412327.38863617.5217654.0
194.86023727.85688217.5485007.0
194.73490027.80356817.5337864.0
195.21377327.68915617.6438184.0
195.46688327.63619417.70317437.0
195.11840227.76122117.4956835.0
195.44246127.95780015.1584706.0
195.24513628.02685117.64242190.0
195.21475528.04287313.6038748.0
194.92607527.85886716.95916911.0
194.89986827.90611217.6399054.0
195.07366927.95526213.7726899.0
195.04483227.99094314.77123078.0
194.84741427.91160413.4516899.0
195.16057528.01606016.2037594.0
195.17734128.11630517.0186206.0
195.08209928.12152117.6966659.0
195.54976928.17322316.3638950.0
194.62947927.89550317.6766804.0
194.56269328.12590714.6327428.0
194.86957228.08577916.5863650.0
194.91736028.17934717.3235437.0
195.16950628.51989315.8698901.0

Now we have a table with all galaxies that either have a redshift measurement in SDSS or a velocity value obtained by (Rines et al. 2003).

Overall, this sample within 40arcmin of the center of the Coma cluster contains 514 galaxies.

Before we start the analysis of the data, let’s look at the galaxies in the sample by loading the table into the Aladin Lite widget. They will appear as catalog_1 and with a different color than the bright stars we added before.

aladin.add_table(complete_sample)

scroll up to the aladin widget to see the table

Determine velocity distribution, cluster average velocity, and velocity dispersion#

Based on the 514 galaxies, we can now analyze the recession velocity and velocity dispersion of the Abell 1656 galaxy cluster. First, we visualize the recession velocity distribution of the entire sample:

# bins is the number of bars in the histogram
plt.hist(data=complete_sample, x="cz", bins=50);
_images/efc106939adbb83cf52a4ff35625ffdf9247cce87294397b1f20d44c8e1b419b.png

Note how there is a large range of recession velocities in our sample. We are only interested in the range of recession velocities of the Coma cluster. These are around the peak at low velocities. Thus, we restrict our sample to a subset df_Coma to recession velocities between 3000 and 11000 km/s:

Coma = complete_sample[
    (complete_sample["cz"] > 3000.0) & (complete_sample["cz"] < 11000.0)
]  # here we apply a mask with two conditions to the sample based on criteria about the column "cz"

plt.hist(data=Coma, x="cz")
plt.xlabel(
    r"$cz~(\mathrm{km/s})$"
);  # the r before the string allows to write the label in LaTex
_images/e339c9e42abb3c61fdaf6e3ced3159c4297a3c7d92a2b96af2eac2cbcc6a8b47.png

This subset corresponds to galaxies in the vicinity of the cluster (both spatially and in recession velocity). Let’s calculate the mean recession velocity of the cluster and its velocity dispersion:

print(
    f"The mean velocity in Coma is {Coma['cz'].mean()} km/s."
    f" Its velocity dispersion (i.e. standard deviation) is {Coma['cz'].std()} km/s.",
)
The mean velocity in Coma is 6974.036149625936 km/s. Its velocity dispersion (i.e. standard deviation) is 1138.227427194869 km/s.

It is in agreement with more refined analyses (e.g. Sohn et al. 2017, ApJS, 229, 20).

When looking back at the query results for the Rines et al. (2003) catalog, you can check again the table that describes the full cluster. The mean recession velocity cz = 6973km/s and dispersion sigmap_3s_ = 1042km/s for the Coma cluster are is also in good agreement with our results.

results_test_rines[0]
Table length=1
Clustern_ClusterRAJ2000DEJ2000czsigmap_3s_sigmap_ca_LXTXR
km / skm / skm / s1e+36 WkeV
str5str1str10str9int32int16int16float32float32uint8
A1656g12 59 31.9+27 54 106973104295718.08.02

Search for Hubble Space Telescope (HST) spectra from the Coma Cluster#

We now want to find out whether there are HST spectra available for the galaxies that had neither a redshift in SDSS nor a velocity in the catalog Rines et al. (2003).

We use the Simple Spectral Access (SSA) protocol from the IVOA to query the Mikulski Archive for Space Telescopes (MAST). Once again, we look at an area of 40arcmin around the center of the Coma Cluster.

mast_ssa_service = pyvo.dal.SSAService(
    "https://archive.stsci.edu/ssap/search2.php?id=HST&",
)
diameter = u.Quantity(2 * 40.0, unit="arcmin")
position = SkyCoord.from_name("A1656")
mast_hst_results = mast_ssa_service.search(pos=position, diameter=diameter)
mast_hst_results
<Table length=21>
   ra_obs    dec_obs           coord_targ                 coord_obs         ...       fluxunit                 fluxucd             fluxcal    ang_sep
    deg        deg                deg                        deg            ...                                                                      
  float64    float64           float64[2]                 float64[2]        ...        object                   object              object    float64
----------- ---------- -------------------------- ------------------------- ... -------------------- --------------------------- ------------ -------
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ...              count/s arith.rate;phot.count;em.wl UNCALIBRATED  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ...              count/s arith.rate;phot.count;em.wl UNCALIBRATED  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ...              count/s arith.rate;phot.count;em.wl UNCALIBRATED  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ...              count/s arith.rate;phot.count;em.wl UNCALIBRATED  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
195.2537375 28.3289472 195.2537375 .. 28.32894722 195.2537375 .. 28.3289472 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.271
   195.0925    28.4005        195.0925 .. 28.4005       195.0925 .. 28.4005 ... erg/cm**2/s/angstrom     phot.flux.density;em.wl     ABSOLUTE  26.255

Note that mast_hst_results is not a list of tables as we had for astroquery queries. This time, we get a pyvo resultset. Thus the methods to handle the resultset are slightly different but can still be printed out with the dir() function which is generic in python. Let’s find out which columns are available:

print("--- Available methods: ", dir(mast_hst_results))
print("--- Name of columns: ", mast_hst_results.fieldnames)
--- Available methods:  ['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__iter__', '__le__', '__len__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_adhocservices', '_findinfos', '_findresultsresource', '_findresultstable', '_findstatus', '_findstatusinfo', '_fldnames', '_from_result_url', '_infos', '_resultstable', '_session', '_status', '_url', '_votable', 'broadcast_samp', 'cursor', 'fielddescs', 'fieldname_with_ucd', 'fieldname_with_utype', 'fieldnames', 'from_result_url', 'get_adhocservice_by_id', 'get_adhocservice_by_ivoid', 'getcolumn', 'getdesc', 'getrecord', 'getvalue', 'iter_adhocservices', 'iter_datalinks', 'queryurl', 'resultstable', 'status', 'table', 'to_table', 'votable']
--- Name of columns:  ('ra_obs', 'dec_obs', 'coord_targ', 'coord_obs', 'url', 'object', 'datalen', 'exposure', 'date_obs', 'tstart', 'tmid', 'tstop', 'format', 'radecsys', 'equinox', 'preview', 'representative', 'min_wavelength', 'max_wavelength', 'title', 'timesys', 'specsys', 'vover', 'vodate', 'author', 'collection', 'ds_ident', 'cr_ident', 'date', 'version', 'instrume', 'dssource', 'specband', 'der_snr', 'spec_val', 'spec_bw', 'spec_fil', 'spec_rp', 's_fov', 'aperture', 'telescop', 'fluxavg', 'fluxmax2', 'min_flux', 'max_flux', 'min_error', 'max_error', 'filesize', 'spectralaxisname', 'fluxaxisname', 'spectralsi', 'fluxsi', 'spectralunit', 'fluxunit', 'fluxucd', 'fluxcal', 'ang_sep')
for observation in mast_hst_results:
    print(observation["object"])
WAVE
QSO-1258+285
QSO-1258+285
QSO-1258+285
QSO-1258+285
WAVE
QSO-1258+285
QSO-1258+285
QSO-1258+285
QSO-1258+285
WAVE
QSO-1258+285
QSO-1258+285
QSO-1258+285
QSO-1258+285
WAVE
QSO-1258+285
QSO-1258+285
QSO-1258+285
QSO-1258+285
1257+2840

Simbad

Often Quasars are further away than the Coma cluster, so let’s check quickly on Simbad whether this source is interesting for further analysis. Usually a Simbad query would only return information on the object’s identifier and coordinates. We are, however, also interested in its redshift, so we first create a customised Simbad query (as we did above for VizieR, for more details see here) and then submit the query.

custom_Simbad = Simbad()
custom_Simbad.add_votable_fields("z_value")
qso_table = custom_Simbad.query_object("QSO 1258+285")
qso_table
Table length=1
MAIN_IDRADECRA_PRECDEC_PRECCOO_ERR_MAJACOO_ERR_MINACOO_ERR_ANGLECOO_QUALCOO_WAVELENGTHCOO_BIBCODEZ_VALUESCRIPT_NUMBER_ID
"h:m:s""d:m:s"masmasdeg
objectstr13str13int16int16float32float32int16str1str1objectfloat64int32
QSO B1258+283513 01 00.8675+28 19 44.74114140.0490.04090AO2020yCat.1350....0G1.36102001

As you can see in the column before last, the Quasar is at a redshift of 1.36. This is far beyond the Coma Cluster. Therefore, we focus on the source ‘1257+2840’ for now. 1257+2840 is the last source in the list: we assign it to a new variable (interesting_observation). Then we exploit the functionalities of resultset to find out where the data is and what kind of file it is.

interesting_observation = mast_hst_results[-1]
observation_url = interesting_observation.getdataurl()
print(observation_url)
http://archive.stsci.edu/pub/vospectra/fos2/y1hi1402t_vop.fits

A quick analysis of the discovered spectrum#

With the previous step, we obtained a link to a fits file which we can download and open with astropy.

spectrum_fits = fits.open(observation_url)
spectrum_fits
[<astropy.io.fits.hdu.image.PrimaryHDU object at 0x7fe0c26b11d0>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x7fe0c2386710>]

we see that spectrum_fits is a list of two astropy HDU (Header Data Unit) objects. Let’s look at the second one (with index 1)

spectrum_fits[1].header
XTENSION= 'BINTABLE'           /Binary table written by MWRFITS v1.4b           
BITPIX  =                    8 /Required value                                  
NAXIS   =                    2 /Required value                                  
NAXIS1  =                24768 /Number of bytes per row                         
NAXIS2  =                    1 /Number of rows                                  
PCOUNT  =                    0 /Normally 0 (no varying arrays)                  
GCOUNT  =                    1 /Required value                                  
TFIELDS =                    3 /Number of columns in table                      
COMMENT                                                                         
COMMENT  *** Column names ***                                                   
COMMENT                                                                         
TTYPE1  = 'WAVE    '           /                                                
TTYPE2  = 'FLUX    '           /                                                
TTYPE3  = 'SIGMA   '           /                                                
COMMENT                                                                         
COMMENT  *** Column formats ***                                                 
COMMENT                                                                         
TFORM1  = '2064E   '           /                                                
TFORM2  = '2064E   '           /                                                
TFORM3  = '2064E   '           /                                                
COMMENT                                                                         
COMMENT   *** Column units ***                                                  
COMMENT                                                                         
TUNIT1  = 'angstrom'           / wavelength unit is Angstrom                    
TUNIT2  = 'erg/cm**2/s/angstrom' / flux units                                   
TUNIT3  = 'erg/cm**2/s/angstrom' / sigma units                                  
COMMENT                                                                         
COMMENT   *** Column limits ***                                                 
COMMENT                                                                         
TDMIN1  =              1087.14 / [angstrom] Min Value for Field 1               
TDMAX1  =              1605.76 / [angstrom] Max Value for Field 1               
TDMIN2  =         -3.26402E-14 / [erg/cm**2/s/angstrom] Min Value for Field 2   
TDMAX2  =          1.39739E-13 / [erg/cm**2/s/angstrom] Max Value for Field 2   
TDMIN3  =          0.00000E+00 / [erg/cm**2/s/angstrom] Min Value for Field 3   
TDMAX3  =          5.63300E-14 / [erg/cm**2/s/angstrom] Max Value for Field 3   
COMMENT                                                                         
COMMENT   *** Column UCDs, UTYPEs ***                                           
COMMENT                                                                         
TUCD1   = 'em.wl'              / Wavelength UCD                                 
TUCD2   = 'phot.flux.density;em.wl' / Flux UCD                                  
TUCD3   = 'stat.error;em.wl'   / Sigma UCD                                      
TUTYP1  = 'Spectrum.Data.SpectralAxis.Value' / Wavelength UTYPE                 
TUTYP2  = 'Spectrum.Data.FluxAxis.Value' / Flux UTYPE                           
TUTYP3  = 'Spectrum.Data.FluxAxis.Accuracy.StatError' / Sigma UTYPE             
                                                                                
          GENERAL KEYWORDS                                                      
                                                                                
EXTNAME = 'Spectral Container' / HST VO spectral container format               
EQUINOX =              2000.00 / Coordinates precessed to J2000                 
RADECSYS= 'FK4'                /                                                
TIMESYS = 'UTC'                / Time system                                    
MJDREF  =              0.00000 / [d] MJD zero point for times                   
SPECSYS = 'TOPOCENT'           / no velocity corrections applied                
VOCLASS = 'SPECTRUM V1.0'      / VO Data Model                                  
VOSEGT  = 'SPECTRUM'           / Segment type                                   
DATALEN =                 2064 / [ ] Number of points in spectrum               
DSSOURCE= 'POINTED'            / Survey or Pointed                              
                                                                                
          DATA ID KEYWORDS                                                      
                                                                                
DATE    = '2008-06-28'         / Processing-Creation Date                       
VODATE  = '2010-08-12'         / Curation Date                                  
VOREF   = 'http://archive.stsci.edu/vodocs/hst/fos' / URL for Documentation     
DS_IDENT= 'ads/sa.hst#y1hi1402t' / Publishers Data Set ID                       
VERSION = '1.0'                / Processing Version                             
COLLECT = 'HST/FOS'            / Data Set Collection                            
INSTRUME= 'FOS'                / Instrument                                     
TELESCOP= 'HST'                / Telescope                                      
FILENAME= 'y1hi1402t_vop.fits' / File Name                                      
CR_IDENT= 'y1hi1402t'          / Internal dataset ID                            
CRETYPE = 'Archival'           / Not an on-the-fly dataset                      
                                                                                
          CURATION KEYWORDS                                                     
                                                                                
VOPUB   = 'MAST'               / VO Publishing Authority                        
AUTHOR  = 'HST project'        / Creator                                        
VOPUBID = 'ivo://mast.stsci.edu' / URI for VO Publisher                         
VOLOGO  = 'http://archive.stsci.edu/images/100.mastlogo.gif' / MAST logo        
VOVER   = '1.0'                / VO Curation version                            
VORIGHTS= 'PUBLIC'             / Data is public                                 
CONTRIB1=                 4952 / PEP proposal identifier                        
CONTACT = 'archive help desk'  /                                                
EMAIL   = 'archive@stsci.edu'  /                                                
                                                                                
          TARGET KEYWORDS                                                       
                                                                                
OBJECT  = '1257+2840'          / proposer's target name                         
TITLE   = 'y1hi1402t, 1257+2840' / filename, object                             
RA_TARG =       195.0925000000 / right ascension of target                      
DEC_TARG=        28.4005000000 / declination of target                          
                                                                                
          WCS Paper 3 Keywords                                                  
                                                                                
1S2_1   = 'WAVE'               / Column name with spectral coordinates          
1CYTP2  = 'WAVE-TAB'           / Spectral coordinate is wavelength              
1S3_1   = 'WAVE'               / Column name with spectral coordinates          
1CYTP3  = 'WAVE-TAB'           / Spectral coordinate is wavelength              
                                                                                
          COVERAGE: SPATIAL                                                     
                                                                                
RA      =       195.0925000000 / right ascension of aperture                    
DEC     =        28.4005000000 / declination of aperture                        
APERTURE= '4.3x4.3'            / [arcsec] Aperture (width or lengthxwidth)      
LONGSTRN= 'OGIP 1.0'           / The OGIP long string convention may be used.'  
COMMENT   This FITS file may contain long string keyword values that are        
COMMENT   continued over multiple keywords.  This convention uses the  '&'      
COMMENT   character at the end of a string which is then continued              
COMMENT   on subsequent keywords whose name = 'CONTINUE'.                       
FTPRT_AP= 'Polygon J2000 195.09191171 28.40116854 195.09174013 28.39998032 &' / 
CONTINUE= '195.09308871 28.39983280 195.09325945 28.40101834' / aperture footpri
                                                                                
          COVERAGE: TEMPORAL                                                    
                                                                                
DATE-OBS= '1993-07-15T04:07:11' / UT observation start time                     
EXPOSURE=             759.9844 / [s] exposure duration                          
TSTART  =       49183.16861829 / [d] MJD exposure start time                    
TSTOP   =       49183.17779164 / [d] MJD exposure stop time                     
TMID    =       49183.17320497 / [d] MJD exposure mid time                      
                                                                                
          COVERAGE: SPECTRAL                                                    
                                                                                
SPECSDIM= '1.0E-10 L'          /                                                
SPEC_CAL= 'ABSOLUTE'           / Calibrated spectral coordinate                 
SPEC_ERR=                  0.2 / [angstrom] Wavelength statistical error        
SPECBAND= 'Far-UV'             / Spectral Band(s)                               
SPEC_VAL=             1346.450 / [angstrom] Mean Wavelength                     
SPEC_BW =              518.618 / [angstrom] Bandpass Width Wmax - Wmin          
SPEC_FIL=                 1.00 / [ ] No gaps between points                     
SPEC_RP =               1300.0 / [ ] average resolving power W/dW               
                                                                                
          COVERAGE: OBSERVABLE                                                  
                                                                                
FLUXSDIM= '1.0E+7 ML-1T-3'     / Flux SI Dimensions                             
FLUX_CAL= 'ABSOLUTE'           / Calibrated or uncalibrated Fluxes              
STAT_ERR=                    3 / [ ] flux error in percent                      
DER_SNR =                 1.15 / [ ] Mean signal-to-noise ratio                 
FLUXAVG =           4.6306E-15 / [erg/cm**2/s/angstrom] Average Flux            
FLUXMAX2=           8.7484E-14 / [erg/cm**2/s/angstrom] 98% Maximum Flux        
SPECDISP=                0.251 / [angstrom/pixel] Dispersion                    
                                                                                

From the fits information and the header, it appears that we have three columns (listed in one axis though):

  • wavelength in Angstrom,

  • flux and flux error in \(\mathrm{erg \cdot cm}^{-2} \mathrm{s}^{-1} \mathrm{\r{A}}^{-1}\).

For a first quick look we can plot the spectrum:

fig = plt.figure(figsize=(10.0, 8.0))
ax = fig.add_axes([0.17, 0.17, 0.77, 0.77])
ax.plot(spectrum_fits[1].data[0][0], spectrum_fits[1].data[0][1])
ax.set_xlabel(r"Wavelength ($\mathrm{\AA}$)", fontsize=14)
ax.set_ylabel(
    r"Flux ($\mathrm{erg \cdot cm}^{-2} \mathrm{s}^{-1} \mathrm{\AA}^{-1}$)",
    fontsize=14,
)
Text(0, 0.5, 'Flux ($\\mathrm{erg \\cdot cm}^{-2} \\mathrm{s}^{-1} \\mathrm{\\AA}^{-1}$)')
_images/1a18f1c3cfa8b0943c8e918079197c6251447186a76fd449fdf59da3449831a0.png

It is a spectrum in the ultraviolet with two visible emission lines, one around \(1220\mathrm{\r{A}}\) and one around \(1330\mathrm{\r{A}}\). We know that wavelength at rest of the Lyman \(\alpha\) line is at \(1216\mathrm{\r{A}}\). This spectrum might thus show Ly\(\alpha\) (atomic hydrogen, HI) emission of the Milky Way (hardly redshifted) and a redshifted extragalactic source.

To investigate this further, we use the specutils package. First, we define a 1D spectrum: the data format that specutils accepts.

flux_unit = u.erg / u.cm**2 / u.s / u.Angstrom
spectrum = Spectrum1D(
    spectral_axis=spectrum_fits[1].data[0][0] * u.Angstrom,
    flux=spectrum_fits[1].data[0][1] * flux_unit,
)
spectrum
<Spectrum1D(flux=<Quantity [ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00, ...,
           -5.2580507e-16,  0.0000000e+00,  0.0000000e+00] erg / (Angstrom s cm2)>, spectral_axis=<SpectralAxis [1087.1409, 1087.3932, 1087.6455, ..., 1605.2651, 1605.5122, 1605.7592] Angstrom>)>

Now we can use specutils_analysis functions to analyze the spectrum. Let’s find the centroid of the two lines.

# find the first peak, between 1200 and 1260 Angstrom
centroid_Milky_Way = spec_ana.centroid(
    spectrum,
    SpectralRegion(1200 * u.Angstrom, 1260 * u.Angstrom),
)
# find the second peak, between 1300 and 1370 Angstrom
centroid_second_peak = spec_ana.centroid(
    spectrum,
    SpectralRegion(1300 * u.Angstrom, 1370 * u.Angstrom),
)
print("The centroid of the first peak is located at: ", centroid_Milky_Way)
print("The centroid of the second peak is located at: ", centroid_second_peak)
The centroid of the first peak is located at:  1216.476806640625 Angstrom
The centroid of the second peak is located at:  1332.959228515625 Angstrom

Indeed the first peak is centered around the rest wavelength of the Ly\(\alpha\) line. We may thus assume that this is the HI emission from the Milky Way in the foreground. Now assuming that the second line is also Ly\(\alpha\) emission, let’s calculate the redshift and recession velocity:

rest_Ly_alpha = 1216.0
redshift_z = (centroid_second_peak.value - rest_Ly_alpha) / rest_Ly_alpha
cz = redshift_z * c  # speed of light in km/s, defined before
print(
    f"The source has a redshift of {round(redshift_z, 3)}"
    f" and a recession velocity of {round(cz, 2)} km/s",
)
The source has a redshift of 0.096 and a recession velocity of 28835.84 km/s

Although this source is much closer than the Quasar, it is still further away than the Coma Cluster and thus not a member of the Cluster.

An alternative to using specutils we can also use more generic python packages and fit the emission lines with simple Gaussians. Let’s define the Gaussian:

def gauss(x, height, peak_position, sigma):
    """Gaussian 1d function.

    Parameters
    ----------
    x : numpy 1d array or a list
        list or array of wavelengths
    height : float
        multiplicative factor for the height of the gaussian function
    peak_position : float
        wavelength corresponding to the peak position
    sigma : float
        standard deviation that controlls with of the peak

    Returns
    -------
    numpy array
    """
    return height * np.exp(
        -((np.asarray(x) - peak_position) ** 2.0) / (2 * sigma**2.0),
    )

Next, we select the two parts of the spectrum where the emission lines are:

# separate wavelengths and fluxes in different objects
spectrum_wavelengths = spectrum_fits[1].data[0][0]
spectrum_flux = spectrum_fits[1].data[0][1]

# make a mask to select wavelengths between 1190 and 1240 Angstrom
mask_Milky_Way = np.where(
    (spectrum_wavelengths < 1240.0) & (spectrum_wavelengths > 1190),
)
wavelengths_Milky_Way = spectrum_wavelengths[mask_Milky_Way]
flux_Milky_Way = spectrum_flux[mask_Milky_Way]

# make a mask to select wavelengths between 1300 and 1380 Angstrom
mask_second_peak = np.where(
    (spectrum_wavelengths < 1380) & (spectrum_wavelengths > 1300),
)
wavelengths_second_peak = spectrum_wavelengths[mask_second_peak]
flux_second_peak = spectrum_flux[mask_second_peak]

The Gaussian fit is done with the curve_fit function of the scipy.optimize library.

popt_Milky_Way_line, pcov_Milky_Way_line = curve_fit(
    gauss,
    wavelengths_Milky_Way,
    flux_Milky_Way,
    p0=[1.25e-13, 1220.0, 10.0],
)
print(
    "The first peak, attributed to the Milky Way, has a",
    f" central wavelength of {round(popt_Milky_Way_line[1], 2)} ",
    f"+/- {round(np.sqrt(np.diag(pcov_Milky_Way_line))[1], 2)} km/s.",
)


popt_second_line, pcov_second_line = curve_fit(
    gauss,
    wavelengths_second_peak,
    flux_second_peak,
    p0=[0.2e-13, 1330.0, 10.0],
)
print(
    "The first peak for the other object has a",
    f"central wavelength of {round(popt_second_line[1], 2)}",
    f" +/- {round(np.sqrt(np.diag(pcov_second_line))[1], 2)} km/s.",
)
The first peak, attributed to the Milky Way, has a  central wavelength of 1215.38  +/- 0.13 km/s.
The first peak for the other object has a central wavelength of 1329.18  +/- 0.6 km/s.

These values are close to the results from the specutils package.

To further check the fitting results, we plot the data, our model, and residuals:

fig, axes = plt.subplots(2, 1)

# First plot is the data and models
axes[0].plot(wavelengths_Milky_Way, flux_Milky_Way)
axes[0].plot(
    wavelengths_Milky_Way,
    gauss(wavelengths_Milky_Way, *popt_Milky_Way_line),
    ls="-",
    color="k",
)
axes[0].plot(wavelengths_second_peak, flux_second_peak)
axes[0].plot(
    wavelengths_second_peak,
    gauss(wavelengths_second_peak, *popt_second_line),
    ls="-",
    color="k",
    label="gaussian fits",
)
axes[0].legend()
axes[0].set_ylabel(r"Flux [erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$]")


# Second plot is the residuals = data - model
axes[1].plot(
    wavelengths_Milky_Way,
    flux_Milky_Way - gauss(wavelengths_Milky_Way, *popt_Milky_Way_line),
)
axes[1].plot(
    wavelengths_second_peak,
    flux_second_peak - gauss(wavelengths_second_peak, *popt_second_line),
)
axes[1].axhline(0, ls="--", c="k")
axes[1].set_xlabel(r"Wavelength [$\mathrm{\AA}$]")
axes[1].set_ylabel("Residual")
Text(0, 0.5, 'Residual')
_images/852fe4bcbafc3e14fc8e0d0d09c7ba420a99a9b0da39f2185d8772b5482cf0fd.png

While fitting a Gaussian to the emission lines provides very good results with regard to the central wavelength (and thus redshift) of the observed objects, the residuals show that the emission from the Milky Way is much more complicated.

# End of the tutorial
spectrum_fits.close()