Historical Quasar Observations#

Learning Goals#

By the end of this tutorial, you will:

  • Understand how to query data on a target from the MAST archive.

  • Be able to plot a historic observation coverage plot of a target.

  • Learn to plot a spectrum from Hubble


Quasars are extremely luminous astronomical objects that can be found at the center of some galaxies. They are powered by gas spiraling at high velocity into a super-massive black hole. The brightest quasars are capable of outshining all the stars in their galaxy; they can be seen from billions of light-years away. The first quasar ever discovered is called 3c273. It is one of the most luminous quasars and therefore one of the most luminous objects in the observable universe. It is at a distance of 749 Megaparsecs [1 Megaparsec = 1 million parsecs = 3.26 million lightyears] with an absolute magnitude of −26.7, meaning that if it were at a distance of 10 parsecs, it would be as bright in our sky as the Sun.

Quasar 3c273 is the our target in this tutorial. We will first search the MAST Archive for all the observations of this quasar. Then, we will display those observations in a historic coverage plot; that is, a plot of the wavelengths in which 3c273 was observed for a given year. This will help us understand what the history of observations of this quasar looks over time. Finally, we will plot a spectrum from one of the observations.


  • Imports

  • Historic Observation Coverage

    • Query the MAST Archive

    • Create Plotting Variables

      • Time Data

      • Wavelength Data

      • Mission Names

    • Plotting Historical Observation Coverage

  • Plot a Spectrum

    • Query HST for Spectra

    • List Available Instrument and Filter Combinations

    • Select Desired Observations

    • Filter for Relevant Products

    • Download the Data

    • Plot the Spectrum

  • Exercises


The following cell holds the imported packages. These packages are necessary for running the rest of the cells in this notebook. A description of each import is as follows:

  • numpy to handle array functions

  • pandas to handle date conversions

  • fits from astropy.io for accessing FITS files

  • Table from astropy.table for creating tidy tables of the data

  • matplotlib.pyplot for plotting data

  • Mast and Observations from astroquery.mast for querying data and observations from the MAST archive

from astropy.io import fits
from astropy.table import Table
from astroquery.mast import Mast, Observations

import itertools
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

%matplotlib inline

Historical Observation Coverage#

The MAST Archive has data from as early as the 1970s. In this section, we’ll search for a target by name and examine its observational history. We’ll create a plot of this history, showing when the target was observed, in what wavelength, and by which mission.

Query the MAST Archive by Object Name#

We are going to use the astroquery.mast Observations package to gather our data from the MAST Archive. In this tutorial, we will use the Observations.query_object() function which takes the name of the target and an optional radius. If you don’t specify a radius, 0.2 degrees will be used by default. For more information about queries you can read the astroquery.mast readthedocs.

# Define target name
target_name = "3c273"

# We'll name the data "obs_table" for "observations table"
obs_table = Observations.query_object(target_name)
# Print out the first 10 entries with limited columns, if you want to see a preview
columns = ['intentType', 'obs_collection', 'wavelength_region', 'target_name', 'dataproduct_type', 'em_min']
WARNING: AstropyDeprecationWarning: show_in_notebook() is deprecated as of 6.1 and to create
         interactive tables it is recommended to use dedicated tools like:
         - https://github.com/bloomberg/ipydatagrid
         - https://docs.bokeh.org/en/latest/docs/user_guide/interaction/widgets.html#datatable
         - https://dash.plotly.com/datatable [warnings]
Table length=10
2scienceTESSOpticalTESS FFIimage600.0

Create Plotting Variables#

Before we start parsing our observations table, let’s recall what we want to do with it.

First, we want to plot a historical observation coverage plot, where the horizontal axis will be time and the vertical axis will be observed wavelength. We should label or color each observation according to what mission it corresponds to. So, we are going need variables for:

  • array of times of all observations = times

  • array of wavelengths of all observations = waves

  • array of mission names of all observations = mission

We will want to modify the queried data for easy visualization:

  1. We will want to convert the Modified Julian Date (MJD) to a calendar year so when we plot the timeline, it will be easy to tell when each observation was made.

  2. MAST only holds columns for the minimum and maximum wavelength of the observation, but for the historic observation coverage plot, we only want one wavelength per observation. We will calculate the average of the min/max values, and use this number instead.

Time Data#

We’ll use t_min for our time data, which corresponds to the beginning of the observation. This is “good enough” for our plot, since it will span multiple decades.

# Parse the observations table to get the Time data
obs_times = obs_table["t_min"]

# Convert MJD to Calendar Date:
# Initialize list for times as calendar dates
times = []
# Loop through times queried from MAST
for t in obs_times:
    # Convert MJD to Julian date
    t = t + 2400000.5
    # Convert Julian date to Calendar date 
    time = pd.to_datetime(t, unit = 'D', origin = 'julian')
    # Add converted date to times list

# Change list to numpy array for easy plotting
times = np.array(times)

Wavelength Data#

Again, we’ll compute the average value reported for the wavelength. One complication here is the presence of a database error in some legacy missions. The current standard is to report wavelengths in meters, but some older missions reported the value in nanometers. Work is ongoing to bring these legacy missions up date; the code in the cell below will correctly parse the values, whether or not the database has been fixed.

# Parse observations table to get the wavelength data
wavelength_min = obs_table["em_min"]
wavelength_max = obs_table["em_max"]

# Get the average wavelength
wavelength_avg = (wavelength_max+wavelength_min)/2

# Some older missions have wavelengths that are off by a factor of 10^9
waves = []
for wave in wavelength_avg:
    if wave/1e9 >= 1:
#change list to numpy array
waves = np.array(waves)
/tmp/ipykernel_1973/317441936.py:17: UserWarning: Warning: converting a masked element to nan.
  waves = np.array(waves)

Mission Names#

We’ll use the mission names to create labels for the data. It will make our plot a little noisier, but it’s good to know where the data is coming from.

#Parse the observations table to get the mission names data
mission = obs_table["obs_collection"]
mission = np.array(mission)

Plotting Historical Observation Coverage#

We want to visualize the history of observations of 3c273 according to the wavelength of the observation. We’ll change the color and shape of each point to indicate which mission made the observation.

# Initialize the figure
fig = plt.figure()
ax = fig.add_subplot()
# Set a color for every unique mission name
cm = plt.cm.get_cmap("plasma") 
num_colors = len(np.unique(mission))
ax.set_prop_cycle(color = [cm(1.*i/num_colors) for i in range(num_colors)])
marker = itertools.cycle(('^', 'o', 's','*')) 

# Loop through the mission names
for i in np.unique(mission):
    # Filter times and wavelengths by mission name
    ind = np.where(mission == i)
    # Plot the mission in its color as a scatterplot
    ax.scatter(times[ind], waves[ind], label = i, s = 100, marker = next(marker))

# Place the legend

# Set the label of the x and y axes
plt.xlabel("Time of Observation [Year]", fontsize = 15)
plt.ylabel("Wavelength of Observation [nm]", fontsize = 15)

# Show the plot
/tmp/ipykernel_1973/2250698704.py:7: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
  cm = plt.cm.get_cmap("plasma")

Whoops! Looks like our results are dominated by Spitzer; let’s drop that mission so we can see the other results.

# Initialize the figure
fig = plt.figure()
ax = fig.add_subplot()
# Set a color for every unique mission name
cm = plt.cm.get_cmap("plasma") 
num_colors = len(np.unique(mission))
ax.set_prop_cycle(color = [cm(1.*i/num_colors) for i in range(num_colors)])
marker = itertools.cycle(('^', 'o', 's','*')) 

# Loop through the mission names
for i in np.unique(mission):
    # Filter times and wavelengths by mission name
    ind = np.where(mission == i)
    # Plot the mission in its color as a scatterplot
    if i != "SPITZER_SHA":
        ax.scatter(times[ind], waves[ind], label = i, s = 100, marker = next(marker))

# Place the legend

# Set the label of the x and y axes
plt.xlabel("Time of Observation [Year]", fontsize = 15)
plt.ylabel("Wavelength of Observation [nm]", fontsize = 15)

# Show the plot
/tmp/ipykernel_1973/3164939797.py:7: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
  cm = plt.cm.get_cmap("plasma")

This gives us a much better view of the wavelength coverage. We have good coverage of wavelengths from the mid-90s, thanks to Hubble. Let’s take a look at some of this data and create a spectrum.

Plot a Spectrum#

Now, we want to plot a spectrum from one of our observations, where the x-axis will be wavelength and the y-axis will be flux (brightness).

We can use our historical observational coverage plot to choose which observation to plot. Let’s pick one from the Hubble Space Telescope (HST). Its high resolution coverage of the ultraviolet makes emission and absorption lines in this region clear; this can help us deduce the composition of 3c273.

Query HST for Spectra#

# Query all the observations of 3c273 from the Hubble Space Telescope
hst_table = Observations.query_criteria(objectname="3c273",radius="10 arcsec", 
                                        dataproduct_type="spectrum", obs_collection="HST")

# Let's print out some relevant columns of this table
columns = ["instrument_name","filters","target_name","obs_id","calib_level","t_exptime"]
Table length=10

List Available Instrument and Filter Combinations#

Most telescopes will have multiple instruments and observing modes. Here we’ll print a summary of the filters and instruments that are available for our search results. This is useful to us because the different instruments aboard Hubble will cover different wavelength ranges.

In our table, we’ll also create two new columns: average exposure time and maximum exposure time. This can help constrain searches for faint objects, or targets that need a longer exposure to be fully resolved.

hst_table['count'] = 1
columns = hst_table.group_by(["instrument_name","filters"])
summary_table = columns["instrument_name","filters","count"].groups.aggregate(np.sum)

# Create two new columns: the average exposure time, and the maximum
summary_table["avg_exptime"] = columns['t_exptime'].groups.aggregate(np.mean)
summary_table["max_exptime"] = columns['t_exptime'].groups.aggregate(np.max)
summary_table["avg_exptime"].format = ".1f"
summary_table["max_exptime"].format = ".1f"

#Take a look at the summary table
Table length=23

Select Observations with a Specific Instrument and Filter Combination#

We are interested in an ultraviolet observation that has an appropriate number of observations. Many of these instrument and filter combinations only have 1 or 2 observations. The COS G130M data looks like a good possibility. Let’s look at the observations for that mode.

g130m_table = hst_table['obsid','obs_id','target_name','calib_level',

# Print out the table of data for this specific filter configuration
Table length=9
# We'll take the longest exposure in this filter and plot the spectrum
sel_table = g130m_table[np.argmax(g130m_table['t_exptime'])]
#Take a look at the selected observation's data table
Row index=5

Filtering for Relevant Products#

Our table of selected observations includes not just the spectra but also RAW files, dark scans, bias scans, and others, all of which are used in the calibration of the data. These are not necessary for plotting the spectrum, so we’ll filter them out by selecting only the Minimum Recommended Products which are marked as intended for science.

# Query the observations from MAST to get a list of products for our selected observation
data_products = Observations.get_product_list(sel_table)

# Get the minimum required products
filtered = Observations.filter_products(data_products, productType='SCIENCE', 
                                        productGroupDescription='Minimum Recommended Products')

# Let's take a look at the products available for our selected observation
Table masked=True length=1
24842820HSTspectrumlbgl31030DADS XSM file - Calibrated combined extracted 1D spectrum COSDmast:HST/product/lbgl31030_x1dsum.fitsSCIENCEMinimum Recommended ProductsX1DSUM--CALCOS3.4.312038lbgl31030_x1dsum.fits181440024842820PUBLIC3G130M

Download the Data to Plot#

# Download the data for our selected observation
data = Observations.download_products(filtered)
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/lbgl31030_x1dsum.fits to ./mastDownload/HST/lbgl31030/lbgl31030_x1dsum.fits ...

We’ve downloaded a Flexible Image Transport System (FITS) file. This is a very common file type used in astronomy for holding data of multiple dimensions. FITS files can hold images but can also contain spectral and temporal information.

We can read the columns of our FITS file to see that it holds two segements of data for this observation, FUVA and FUVB. These are two different subsections of the far-ultraviolet spectrum that Hubble observes.

To plot a spectrum, we’ll need to get the data from the Wavelength and Flux columns.

#Take a peek at the FITS file we downloaded
filename = data['Local Path'][0]

#Read the table with the spectrum from the FITS file
tab = Table.read(filename)
Filename: ./mastDownload/HST/lbgl31030/lbgl31030_x1dsum.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     180   ()      
  1  SCI           1 BinTableHDU    318   2R x 16C   [4A, 1D, 1J, 16384D, 16384E, 16384E, 16384E, 16384E, 16384E, 16384E, 16384E, 16384E, 16384E, 16384E, 16384I, 16384E]   
WARNING: UnitsWarning: 'erg /s /cm**2 /angstrom' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
Table length=2
sAngstromerg / (Angstrom s cm2)erg / (Angstrom s cm2)ct / sct / sctctctctct / sct / s
FUVA1110.016163841297.2715719782857 .. 1460.57400135876260.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.0128 .. 1280.0 .. 0.0
FUVB1110.016163841143.983070446528 .. 1307.27681551990620.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.0128 .. 1280.0 .. 0.0

Plot the Spectrum#

# Make the figure and set the font size globally
plt.rcParams.update({"font.size": 14})

# Gather the arrays from our data table
waves = tab['WAVELENGTH']
fluxes = tab["FLUX"]
segment = tab['SEGMENT']

# You'll notice from our data table that there are two segments to this observation, FUV A and FUV B
# Let's parse the spectra by their segment and plot them separately
ind_A = np.squeeze(np.where(fluxes != 0) and np.where(segment == 'FUVA'))
waves_A = waves[ind_A]
fluxes_A = fluxes[ind_A]
ind_B = np.squeeze(np.where(fluxes != 0) and np.where(segment == 'FUVB'))
waves_B = waves[ind_B]
fluxes_B = fluxes[ind_B]

# Plot both segments
plt.plot(waves_B, fluxes_B, label = "FUV B", color = 'lightcoral')
plt.plot(waves_A, fluxes_A, label = "FUV A", color = 'turquoise')

# Set the x and y axes labels and the title
plt.xlabel('Wavelength [{}]'.format(tab['WAVELENGTH'].unit))
plt.ylabel('Flux [{}]'.format(tab['FLUX'].unit))
plt.title("HST Spectrum of 3c273")

# Plot the legend

# Give the figure a tight layout (optional)


Recognizing Familiar Emission Lines#

Look at the spectrum we plotted, does anything stand out to you?

In astronomy, spectral features at specific wavelengths are indicative of known elements. For example, an emission line at 1216 Angstroms is called Lyman Alpha. It is produced when an orbital electron of a hydrogen atom drops from the first excited state down to the ground state, emitting a photon.

Does 3c273 have a Lyman Alpha emission line? Plot a vertical line at 1216 Angstroms to find out.


This question is a bit disingenuous. We can add plot a vertical line at the Lyman Alpha wavelength, and we’ll actually find quite a nice fit for our data.

# appending these lines to the last cell above will give us the answer
xval = 1216 #Lyman Alpha in Angstroms
plt.axvline(xval, color = "black", linestyle = "--")
<matplotlib.lines.Line2D at 0x7f1fcd5c9450>

However, our target is an extragalactic quasar, so we should expect some redshifting. Indeed, the peak we’ve found here isn’t from the quasar at all! It’s contamination from a source along our line-of-sight.

The “true” Lyman Alpha emission is actually the enormous, broadened peak in FUV A. As an advanced exercise for the reader, you might try fitting a normal distribution to that peak, determining the center of the wavelength, then calculating the redshift using $\(z = \frac{\lambda_{obs}-\lambda_{actual}}{\lambda_{actual}}\)$


About this Notebook#

Author: Emma Lieb
Last updated: Apr 2023

If you have questions, comments, or other feedback, please contact the Archive HelpDesk at archive@stsci.edu.

Space Telescope Logo

Return to top of page