**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

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.

- 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.

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 3
- astropy
- numpy
- scipy
- matplotlib
- pyds9

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')
```

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.

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)
```

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.

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)
```

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();
```

- 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]
```

`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');
```

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.

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');
```

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
```

`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');
```

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
```