Notebook

Spectroscopy Module¶

Lecturer: Robert Quimby
Jupyter Notebook Author: Robert Quimby
With Contributions from: Cameron Hummels, Matt Hankins, & Leo Singer

This is a Jupyter notebook lesson taken from the GROWTH Summer School 2019. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-school-2019.html

Objective¶

Extract a 1-D spectrum from a 2-D image of a long-slit spectrum, determine wavelength solution, and then measure the redshift of the target.

Key steps¶

determine the mapping between pixel coordinates and lines of constant wavelength
create and subtract a 2D model of the sky background from the science image
determine the mapping between pixel coordinates and lines of constant spatial position
extract (sum) the background subtracted target counts in each wavelength bin
determine the wavelengths for each wavelength bin
identify line features in the extracted spectrum and use the known rest wavelength to compute a redshift
normalize a flat field image wavelength bin by wavelength bin and correct the images for pixel-to-pixel efficiencies
extract a spectrum of a "standard star" with a known flux densities
use the standard star spectrum to calibrate the instrumental response function
use the response function to convert the target counts into flux densities
identify telluric absorption features in the standard star spectrum and correct the target spectrum for these absorption bands

*Note we may not get through all of these steps today. We will do our best to get as far down the list as we can and leave any remaining steps for you as homework.

Required dependencies¶

See GROWTH school webpage for detailed instructions on how to install these modules and packages. Nominally, you should be able to install the python modules with pip install <module>. The external astromatic packages are easiest installed using package managers (e.g., rpm, apt-get).

Python modules¶

python 3
astropy
numpy
scipy
matplotlib
pyds9

External packages¶

ds9 http://ds9.si.edu/site/Home.html

In [1]:

# python modules that we will use
import numpy as np
from astropy.io import fits
from astropy.stats import sigma_clip
from scipy.interpolate import interp1d
from scipy.interpolate import LSQBivariateSpline, LSQUnivariateSpline
from scipy.optimize import fmin
from scipy.optimize import least_squares
import pickle

%matplotlib inline
import matplotlib.pylab as plt
import os
import subprocess

In [2]:

# optional launch DS9 for displaying image data
import pyds9
ds9 = pyds9.DS9()

In [3]:

# change plotting defaults
plt.rc('axes', labelsize=14)
plt.rc('axes', labelweight='bold')
plt.rc('axes', titlesize=16)
plt.rc('axes', titleweight='bold')
plt.rc('font', family='sans-serif')
plt.rcParams['figure.figsize'] = (15, 7)

In [4]:

cwd = os.getcwd()
data_dir = os.path.join(cwd, 'data')

Test dependencies¶

In order for this jupyter notebook to function correctly, we must have some external software installed, as described above. The following step assures that these are installed properly before getting to the rest of the content of this lesson.

In [5]:

dependencies = ['ds9']

def test_dependency(dep):
    try:
        subprocess.Popen(dep, stderr=subprocess.PIPE).stderr.read()
        print("%s is installed properly. OK" % dep)
        return 1
    except ImportError:
        print("===%s IS NOT YET INSTALLED PROPERLY===" % dep)
        return 0
        
i = 0
for dep in dependencies:
    i += test_dependency(dep)
print("\n%i out of %i dependencies installed properly." % (i, len(dependencies)))
if i != len(dependencies):
    print("Please correctly install these programs before continuing.")
else:
    print("You are ready to continue.")

ds9 is installed properly. OK

1 out of 1 dependencies installed properly.
You are ready to continue.

Scipy bug check¶

In [6]:

# test if spline fitting works
x = np.array([0, 1, 2, 2, 3, 4])
y = np.array([0, 1, 2, 2, 3, 4])
try:
    spl = LSQUnivariateSpline(x, y, []) # new scipy crashes here because of duplicate x=2 values
    SPLINE_WORKS = True
    print('Congratulations! Your version of scipy has 1D spline fitting')
except ValueError:
    SPLINE_WORKS = False
    print('Nuts! Your version of scipy does not do proper 1D spline fitting.')

Nuts! Your version of scipy does not do proper 1D spline fitting.

In [7]:

# *** HACK! -- HACK! -- HACK! ***
if SPLINE_WORKS is False:
    from scipy.interpolate import LSQUnivariateSpline as real_LSQUnivariateSpline
    
    # re-define LSQUnivariateSpline
    def LSQUnivariateSpline(x, y, t, weights=None):
        # x must be **strictly** increasing
        # so remove (average over) duplicates
        uniqx, origi, counts = np.unique(x, return_inverse=True, return_counts=True)

        if weights is not None:
            newy = np.bincount(origi, weights=y*weights) / np.bincount(origi, weights=weights)
            newweights = np.sqrt( 1 / np.bincount(origi, weights=1/weights**2) )
        else:
            newy = np.bincount(origi, weights=y/counts[origi])
            newweights = None
        return real_LSQUnivariateSpline(uniqx, newy, t, newweights)   

Load the 2-D science Frame¶

In [8]:

# load the 2D science image data
image = fits.getdata(os.path.join(data_dir,'spec_sci.fits'))

In [9]:

# determine the image pixel distribution (used for displaying below)
sample = sigma_clip(image)
vmin = sample.mean() - 1 * sample.std()
vmax = sample.mean() + 3 * sample.std()

In [10]:

# show the image using matplotlib's `imshow` function
plt.figure(figsize=(15, 3))
plt.imshow(image, origin='lower', cmap='gray', aspect='auto', vmin=vmin, vmax=vmax)
plt.xlabel('Column Number')
plt.ylabel('Row Number');

In [ ]:

# optional: display image in DS9
ds9.set('frame 1')
ds9.set_np2arr(image)
ds9.set('scale zscale');

Take a look at the science image above. The supplied spec_sci.fits image is the red channel of a Keck/LRIS spectrum taken in 2014. The image has already been process to remove bias, and the image has been trimmed and rotated (what we will call "rows" are actually columns on the CCD). This science image is actually a combination of two individual exposures that has been filtered to remove most (but not quite all) cosmic ray events.

Notice the vertical bands of light. These are the sky lines. The spectrum runs from about 5000 Angstroms on the left to about 10000 Angstroms on the right (you will determine the exact range below!). Notice that there is a significantly higher density of sky lines in the near infrared. Many of these lines in the infrared are caused by different molecules in Earth's atmosphere. In particular O-H bonds are the primary culprit for sky lines beyond ~750 angstroms. Later we'll have a closer look at the sky lines to examine the wavelength solution for the spectra.

Notice the two horizontal bands of light. The lower band is the target (a small galaxy), and the upper band is a second object in the slit. If you look closely at the lower band (the target trace), you will see a number of emission lines. You will use these to determine the target redshift.

Set up some variables for later use¶

How big is the image? What are the pixel coordinates near the center of the image? Create 1-D and 2-D arrays with indices for each axis.

In [11]:

# get the a pixel coordinate near the image center
ny, nx = image.shape
cy, cx = ny//2, nx//2

# create 1d arays of the possible x and y values
xs = np.arange(nx)
ys = np.arange(ny)

# pixel coordinates for each pixel
yvals, xvals = np.indices(image.shape)

Find the rows dominated by sky light (exclude rows with object light, bad rows)¶

We can use the sky lines to map out lines of constant wavelength. To do this it is best to exclude portions of the image that are affected by light from the target trace, etc.

To find the rows affected by the object traces, we can just marginalize over all the columns to determine the average profile and then note where this is significantly above the background average.

In [12]:

# compute the row averages and normalize so that the background is near 0 and the peaks are near 1
rowaverage = image.mean(axis=1)
rowaverage -= np.median(rowaverage)
rowaverage /= rowaverage.max()

plt.plot(ys, rowaverage)
plt.xlabel('Row Number (y-coordinate)'), plt.ylabel('Normalized Row Average')
plt.grid();

Notes¶

the plot above should show peaks in the rows containing object light
there are some bad rows near the edges of the image
other rows also show significant deviations from the median and should not be used

Based on the plot above, record the row coordinates (ys) that are brighter that some threshold, say 20% of the profile peak. To be conservative, also include the 5 rows above and 5 rows below.

Then create a 2D boolean mask that is True for pixels that are dominated by sky light and False for other pixels.

In [13]:

# find the rows with object light
objrows = ys[rowaverage>0.1]

# add some margin to object rows
ngrow = 5 # number of rows to include above and below object rows
newobjrows = []
for row in objrows:
    newobjrows.extend([row + i for i in np.arange(-ngrow, ngrow + 1)])
objrows = np.unique(newobjrows)

# mask to mark sky rows
skymask = np.ones(image.shape, dtype=bool)
skymask[objrows, :] = False

# also exclude bad rows
badrows = ys[rowaverage < -0.05]
skymask[badrows, :] = False

# rows with mostly sky background light
skyrows = ys[skymask.mean(axis=1) == 1]

With the object traces and bad rows masked, we can also check for cosmic rays and mask these as well. To do this we calculate the median value and standard deviation along each column and then reject pixels that are more than a certain number (typically 3-5) of standard deviations away from the median. These deviant pixels are then noted on the skymask.

In [14]:

# median (unmasked) sky spectrum and standard deviation
medspec = np.median(image[skyrows],axis=0)
stdspec = np.std(image[skyrows],axis=0)

# exclude deviant pixels from the skymask
pull = (image - medspec) / stdspec
w = pull > 5
skymask[w] = False

In [15]:

# show the mask
plt.figure(figsize=(15, 3))
plt.imshow(skymask, origin='lower', aspect='auto');

In [16]:

# optional: display the mask in DS9
mask = np.zeros(image.shape)
mask[skymask] = 1

ds9.set('frame 1')
ds9.set_np2arr(image)
ds9.set('frame 2')
ds9.set_np2arr(mask)
ds9.set('lock frame image');

Let's look at a small section of image up close¶

In [17]:

# cut out a small image "stamp" near the center of the frame
row = cy
col = cx
hwidth = 50
stamp = image[:, col - hwidth : col + hwidth]
ys_stamp = yvals[:, col - hwidth : col + hwidth]
xs_stamp = xvals[:, col - hwidth : col + hwidth]

In [18]:

# show the image stamp
plt.figure(figsize=(15, 5))
extent = (xs_stamp.min(), xs_stamp.max(), ys_stamp.min(), ys_stamp.max())
plt.imshow(stamp, origin='lower', cmap='gray', aspect='auto', vmin=vmin, vmax=vmax, extent=extent)
plt.xlabel('Column Number')
plt.ylabel('Row Number');

Recall the vertical bands are sky lines that mark lines of constant wavelength. Notice that these do not run perfectly parallel to the columns. Rather, the x coordinate for a given wavelength will be a function of the row number.

As the plot should demonstrate, the lines of constant wavelength are slightly tilted with respect to the columns and there is also slight curvature. Thus we can approximate that the x coordinate for a given wavelength will be a quadratic function of the row number.

Note¶

Because wavelength varies along a given column, if we simply plot the counts in each pixel against each pixel's column number then we get a range of values in each column:

In [19]:

# plot stamp values against column numbers
plt.plot(xs_stamp.ravel(), stamp.ravel(), 'r.');
plt.xlabel('Column Number'), plt.ylabel('Counts');

Map out lines of constant wavelength (determine the wavelength for each pixel)¶

We can model the change in wavelength coordinate in arbitrary units, dl, from the wavelength at some reference pixel in terms of the offsets, dx and dy from the reference pixel.

We can then write down a function that takes our model parameters (the slope and curvature of the lines of constant wavelength with respect to the columns), the offsets in the x and y coordinates, dxs and dys, respectively, and returns the wavelength offset from the reference pixel (i.e. dx = dy = 0).

In [20]:

def get_dl_model(params, dxs, dys):
    return dxs + params[0] * dys + params[1] * dys ** 2

To see how this works, make a guess for the slope and curvature parameters and try plotting the wavelength offsets from our reference pixel (the one at x=col and y=row).

In [21]:

# pixel offsets from the refernece pixel
dxs = xs_stamp - col
dys = ys_stamp - row

# parameter guess
#guess = (-0.01, 0) #initial guess
guess = (-1.30240870e-02, -2.58517988e-05) #determined from best fit in the code below

# get the wavelength offsets and plot vs. counts
dls = get_dl_model(guess, dxs, dys)
plt.plot(dls.ravel(), stamp.ravel(), 'r.')
plt.xlabel('Wavelength Offset')
plt.ylabel('Counts');

You should notice that for the right choice of model parameters, the vertical scatter is significantly reduced. This demonstrates one way to decide on the best model parameters: find the model parameters that minimize the vertical scatter. But how do we measure this scatter?

Fit a spline to the counts spectrum above¶

Given the data above (the wavelength offsets, dls, and the stamp count values) we can fit a spline to the resulting curve above and then subtract this off from the curve and measure the scatter. Notice that for every column we get a lot of samples (ny to be precise). Thus we can set the spline knots every pixel (or even closer).

We can define a function that will return the best fit spline object (in the least-squares sense) for our data.

In [22]:

def get_profile_spl(dls, stamp):
    # need to sort the data (and weights) so that the x values are increasing
    x, y = dls.ravel(), stamp.ravel()
    weights = np.sqrt(np.abs(y)) # no