Spectroscopy Module

Lecturer: Robert Quimby
Jupyter Notebook Authors: Robert Quimby & Cameron Hummels

This is a Jupyter notebook lesson taken from the GROWTH Winter School 2018. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-astro-school-2018-resources.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

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

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; this runs VERY slow!
    def LSQUnivariateSpline(x, y, t, weights=None):
        # x must be **strictly** increasing
        # so remove (average over) duplicates
        uniqx = np.unique(x)
        newy = np.zeros(uniqx.shape)
        if weights is not None:
            newweights = np.zeros(uniqx.shape)
        else:
            newweights = None
            
        for i in range(uniqx.size):
            w = (x == uniqx[i])
            newy[i] = np.mean(y[w])
            if weights is not None:
                newweights[i] = np.mean(weights[w]) # probably not right...
                
        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 [11]:
# 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.

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 [12]:
# 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 [13]:
# 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 [14]:
# find the rows with object light
objrows = ys[rowaverage > 0.2]

# 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 [15]:
# median (unmasked) sky spectrum and standard deviation
medspec = np.median(image[skyrows, :], axis=0)
stdspec = np.std(image[skyrows, :], axis=0, ddof=1)

# exclude deviant pixels from the skymask
pull = (image - medspec) / stdspec
w = pull > 5
skymask[w] = False
In [16]:
# show the mask
plt.figure(figsize=(15, 3))
plt.imshow(skymask, origin='lower', aspect='auto');
In [17]:
# optional: display the mask in DS9
mask = np.zeros(image.shape)
mask[skymask] = 1

#optional use of ds9
#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 [18]:
# 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 [19]:
# 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 [20]:
# 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 [21]:
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 [22]:
# pixel offsets from the refernece pixel
dxs = xs_stamp - col
dys = ys_stamp - row

# parameter guess
guess = (-0.01, 0)

# 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 [23]:
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)) # not technically optimal for coadded data, but ok
    wsort = x.argsort()
    x, y, weights = x[wsort], y[wsort], weights[wsort]

    # set locations for spline knots
    t = np.linspace(x.min() + 1, x.max() - 1, np.int(x.max() - x.min()))
    spl = LSQUnivariateSpline(x, y, t, weights)
    return x, y, spl

Lets see how good this looks.

In [24]:
# fit a spline to the data and plot
x, y, spl = get_profile_spl(dls, stamp)

fig, axarr = plt.subplots(2, sharex=True)
axarr[0].plot(x, y, 'r.')
axarr[0].plot(x, spl(x))
axarr[1].plot(x, y - spl(x), 'r.')
plt.ylim(-200, 200)
plt.xlabel('Wavelength Offset');

Notice

  • the counts may vary along the slit (both because of real spatial variations and, more importantly, we have not accounted for efficiency variations along the slit), so there may be some biased scatter
  • the parameter guess above may be imperfect resulting in systematic residuals (we can minimize these below)

Minimize residuals from the spline model to determine the best tilt/curvature model

We need a metric to determine how well the spline fits the data. We could take a simple sum of the squares of the residuals or weight these by the expected poison noise to determine a $\chi^2$ value.

In [25]:
def check_dl_model(params, dxs, dys, stamp):
    dls = get_dl_model(params, dxs, dys)
    x, y, spl = get_profile_spl(dls, stamp)
    chisq = np.sum((stamp - spl(dls)) ** 2 / np.sqrt(np.abs(stamp)))
    return chisq / (stamp.size - len(params))
In [26]:
# see how good our guess is
check_dl_model(guess, dxs, dys, stamp)
Out[26]:
111.18612127340944

Now we just need to change the model parameters to minimize the residuals. We can use fmin for this.

In [27]:
if SPLINE_WORKS is False:
    print("the hacked version of LSQUnivariateSpline is VERY slow...this cell might take a minute or more to run!")

# get the best model parameters for this stamp
params = fmin(check_dl_model, guess, args=(dxs, dys, stamp))
print("best model parameters are", params)
the hacked version of LSQUnivariateSpline is VERY slow...this cell might take a minute or more to run!
Optimization terminated successfully.
         Current function value: 50.174574
         Iterations: 27
         Function evaluations: 54
best model parameters are [-1.30240870e-02 -2.58517988e-05]

You can plug in these parameters for your guess above and re-run the cells to check the residuals if you want.

Find the lines of constant wavelength at other parts of the 2-D image

The tilt/curvature between lines of constant wavelength and the columns can vary across the spectrum. So far we have only determined it for a small portion of the image. Now we can apply the same techniques to map out the relative wavelengths shifts across the image.

To do this, we need to pick a width (number of columns) within which we can expect to get enough sky lines to do the mapping. Considering the left side (shorter wavelength or "blue" portion) of the spectrum, we can set the width at about 400 columns. (Note it is possible to vary the width across the image to take advantage of the extra wavelength information in the near-IR, but this is not required).

In [28]:
# define the column centers for the stamps
hwidth = 200 # width = 2 * hwidth + 1
cols = np.arange(hwidth, nx, 2 * hwidth)
cols
Out[28]:
array([ 200,  600, 1000, 1400, 1800, 2200, 2600, 3000, 3400, 3800])

As above, map out the wavelength shits along columns, one stamp (image section) at a time. It is best to use the sky mask we defined above and possibly some outlier rejection to fend off deviant pixels.

Notice that the code cell below will keep track of the wavelength offsets for all columns, but really we only need the central columns defined above.

Note: the code below might take some time to run. If you are using the hacked version of LSQUnivariateSpline, you should probably skip this cell and load the solution below

In [ ]:
# reference wavelength offsets to the center row
row = cy

# define a 2D array to hold the wavelength offsets for each pixel
lambdas = np.zeros(image.shape) 

# loop over each central column
for col in cols:
    print('col = ', col)
    
    # slice the data
    inds = np.s_[:, col - hwidth : col + hwidth]
    stamp = image[inds]
    mask = skymask[inds]
    dys = yvals[inds] - row
    dxs = xvals[inds] - col

    # initial fit
    params = fmin(check_dl_model, guess, args=(dxs[mask], dys[mask], stamp[mask]))
    
    # check for outliers
    dls = get_dl_model(guess, dxs, dys)
    x, y, spl = get_profile_spl(dls, stamp)
    model = spl(dls)
    pull = (stamp - model) / np.sqrt(np.abs(model))
    w = (pull < 5) & mask
    params2 = fmin(check_dl_model, params, args=(dxs[w], dys[w], stamp[w]))

    # record
    lambdas[inds] = get_dl_model(params2, dxs, dys) + col
In [29]:
if SPLINE_WORKS is False:
    lambdas = pickle.load(open(os.path.join(data_dir,'lambdas.dat'), 'rb'))

# **********************
# save the lambdas
#pickle.dump(lambdas, open('lambdas.dat', 'wb'))
# **********************

Look at a few rows and see how the wavelength offsets vary with column number. We can fit a low-order polynomial to these.

In [30]:
# just plot offsets for a few of the rows across the image
order = 3
for y in range(10, ny, 50):
    p = plt.plot(cols, lambdas[y, cols] - xs[cols], 'o')
    c = np.polyfit(cols, lambdas[y, cols] - xs[cols], order)
    plt.plot(xs, np.polyval(c, xs), c=p[0].get_color(), label='row {}'.format(y))
plt.legend()
plt.xlabel('Column Number')
plt.ylabel('Wavelength Offset from Middle Row');

You may notice that the tilt (the wavelength difference in a given column between the first row and the last row) increases as we move to larger column numbers. Make sure the order is large enough to follow the trends without over fitting the data (i.e. there should not be noticeable wiggles between data points).

Now we can fit for every row. We could do a 2D fit to all the data at once, but simply fitting row by row works here.

In [31]:
# get the lambda values for the entire image (fit)
lambdafit = np.zeros(image.shape)
for y in range(ny):
    c = np.polyfit(cols, lambdas[y, cols] - xs[cols], order)
    lambdafit[y, :] = np.polyval(c, xs) + xs

Now for the fun part! We have wavelength values relative to the central row (in arbitrary units) for every pixel on the image, which means we can do some cool tricks.

Model the 2D sky

We want to measure the light from our target, not the night sky, so lets make those sky lines disappear! We can do a 2-D spline fit to the sky background pixels to create a sky model. As above, we have about ny measurements in each wavelength bin, so we can set the knot points along the wavelength direction about every pixel (or even less!). We do not expect strong variations in the sky brightness along the slit, so we can simply fit a low order polynomial along this direction (which we will approximate for now using the y-coordinates).

In [32]:
# function to fit a 2D spline
def fit_sky(xvals, yvals, image, ky=1, dx=0.5):
    # select knot points in the x (wavelength) direction
    tx = np.arange(xvals.min() + 2, xvals.max() - 2, dx)
    
    # select knot points in the y (spatial) direction
    ty = [] # if there are no knots, the fit will be a poly nomial of degree ky
    
    # fit the 2D spline
    return LSQBivariateSpline(xvals.ravel(), yvals.ravel(), image.ravel(), tx, ty, ky=ky)
In [33]:
# use the (unmasked) sky background pixels and fit the 2D spline
skyfit = fit_sky(lambdafit[skymask], yvals[skymask], image[skymask])

# evaluate the 2D sky at every pixel
sky = skyfit.ev(lambdafit, yvals)
In [35]:
# plot the image, sky model, and differece (and/or display in DS9)
#ds9.set('frame 1')
#ds9.set_np2arr(image)
#ds9.set('frame 2')
#ds9.set_np2arr(sky)
#ds9.set('frame 3')
#ds9.set_np2arr(image - sky);

Notice how well the object trace stands out now. Analogous to the sky lines, the object trace is slightly tilted with respect to the rows. Thus to extract our spectrum we need to map out the y-values where the trace is centered as a function of column number.

Locate object trace

There are actually two object traces in our data, so let's simultaneously find the centers of each of these. The object that makes the lower trace is actually a galaxy, but we can approximate its distribution of light along the slit with a Gaussian function, and we can do the same for the other object. Thus our profile model will be the sum of two Gaussians.

In [36]:
def get_profile_model(params, ys):
    a1, cy1, sigma1, a2, cy2, sigma2 = params
    
    p1 = np.exp(-(ys - cy1)**2 / 2 / sigma1**2) 
    p1 /= p1.max()
    
    p2 = np.exp(-(ys - cy2)**2 / 2 / sigma2**2) 
    p2 /= p2.max()

    return a1 * p1 + a2 * p2

Now that we have our model, we can make a guess at the parameters and check by eye how well these represent the data. Change the guess until the model is a reasonably close match to the data (this need not be perfect).

In [37]:
# get the median for each row
profile = np.median(image - sky, axis=1)

# starting guess for the profile model
guess = (80, 35, 3, 80, 155, 3)
model = get_profile_model(guess, ys)

plt.plot(ys, profile)
plt.plot(ys, model);
plt.xlabel('Row Number')
plt.ylabel('Median Row Counts');

Of course, we can improve the model by minimizing the model residuals (or $\chi^2$/dof).

In [38]:
def get_profile_chisq(params, ys, profile):
    model = get_profile_model(params, ys)
    return np.sum( (profile - model)**2 / np.sqrt(np.abs(profile)) ) / (profile.size - len(params))
In [39]:
# fit for the best model
params = fmin(get_profile_chisq, guess, args=(ys, profile))
print("best fit parameters are", params)

model = get_profile_model(params, ys)
plt.plot(ys, profile)
plt.plot(ys, model)
plt.xlabel('Row Number')
plt.ylabel('Median Row Counts')
plt.grid();
Optimization terminated successfully.
         Current function value: 2.582448
         Iterations: 294
         Function evaluations: 460
best fit parameters are [ 79.44334661  36.23045029   3.53395051  89.76771819 154.95509459
   2.68461778]

We could fit a more complex model that accounts for the vertical shift of the traces with column numbers, but here we simply divide the image up into sections by wavelength and find the Gaussian centroids in each of these.

In [40]:
# fit the profile centered at these columns
hwidth = 50
cols = np.arange(hwidth, nx + 1, 2 * hwidth)

ycenter = np.zeros( (len(cols), 2) )
for icol, col in enumerate(cols):
    stamp = (image - sky)[:, col - hwidth : col + hwidth]
    profile = np.mean(stamp, axis=1)
    params = fmin(get_profile_chisq, guess, args=(ys, profile))
    ycenter[icol, :] = params[[1, 4]]
Optimization terminated successfully.
         Current function value: 0.366600
         Iterations: 546
         Function evaluations: 855
Optimization terminated successfully.
         Current function value: 1.685744
         Iterations: 705
         Function evaluations: 1100
Optimization terminated successfully.
         Current function value: 1.233543
         Iterations: 366
         Function evaluations: 581
Optimization terminated successfully.
         Current function value: 1.765079
         Iterations: 344
         Function evaluations: 540
Optimization terminated successfully.
         Current function value: 3.269195
         Iterations: 313
         Function evaluations: 505
Optimization terminated successfully.
         Current function value: 2.536442
         Iterations: 382
         Function evaluations: 600
Optimization terminated successfully.
         Current function value: 2.471912
         Iterations: 315
         Function evaluations: 505
Optimization terminated successfully.
         Current function value: 2.824128
         Iterations: 314
         Function evaluations: 508
Optimization terminated successfully.
         Current function value: 3.543920
         Iterations: 256
         Function evaluations: 417
Optimization terminated successfully.
         Current function value: 3.341981
         Iterations: 257
         Function evaluations: 421
Optimization terminated successfully.
         Current function value: 3.224893
         Iterations: 252
         Function evaluations: 407
Optimization terminated successfully.
         Current function value: 3.017439
         Iterations: 325
         Function evaluations: 530
Optimization terminated successfully.
         Current function value: 3.589273
         Iterations: 288
         Function evaluations: 466
Optimization terminated successfully.
         Current function value: 4.645856
         Iterations: 305
         Function evaluations: 500
Optimization terminated successfully.
         Current function value: 5.627596
         Iterations: 317
         Function evaluations: 507
Optimization terminated successfully.
         Current function value: 4.259348
         Iterations: 281
         Function evaluations: 462
/home/chummels/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in true_divide
  
Warning: Maximum number of function evaluations has been exceeded.
Optimization terminated successfully.
         Current function value: 7.197677
         Iterations: 372
         Function evaluations: 590
Optimization terminated successfully.
         Current function value: 5.323692
         Iterations: 391
         Function evaluations: 609
Optimization terminated successfully.
         Current function value: 4.998407
         Iterations: 255
         Function evaluations: 419
Optimization terminated successfully.
         Current function value: 9.726273
         Iterations: 344
         Function evaluations: 559
Optimization terminated successfully.
         Current function value: 10.494806
         Iterations: 318
         Function evaluations: 513
Optimization terminated successfully.
         Current function value: 7.791675
         Iterations: 354
         Function evaluations: 581
Optimization terminated successfully.
         Current function value: 4.431137
         Iterations: 342
         Function evaluations: 550
Optimization terminated successfully.
         Current function value: 7.053048
         Iterations: 312
         Function evaluations: 497
Optimization terminated successfully.
         Current function value: 12.994357
         Iterations: 368
         Function evaluations: 583
Optimization terminated successfully.
         Current function value: 7.576511
         Iterations: 356
         Function evaluations: 573
Optimization terminated successfully.
         Current function value: 14.324226
         Iterations: 291
         Function evaluations: 470
Optimization terminated successfully.
         Current function value: 11.501018
         Iterations: 414
         Function evaluations: 662
Optimization terminated successfully.
         Current function value: 18.868989
         Iterations: 411
         Function evaluations: 664
Optimization terminated successfully.
         Current function value: 13.738735
         Iterations: 359
         Function evaluations: 572
Optimization terminated successfully.
         Current function value: 6.686031
         Iterations: 384
         Function evaluations: 622
Optimization terminated successfully.
         Current function value: 5.095648
         Iterations: 490
         Function evaluations: 795
Optimization terminated successfully.
         Current function value: 25.420442
         Iterations: 372
         Function evaluations: 591
Optimization terminated successfully.
         Current function value: 22.220226
         Iterations: 453
         Function evaluations: 702
Optimization terminated successfully.
         Current function value: 11.437006
         Iterations: 491
         Function evaluations: 774
Optimization terminated successfully.
         Current function value: 10.264501
         Iterations: 343
         Function evaluations: 552
Optimization terminated successfully.
         Current function value: 9.764020
         Iterations: 395
         Function evaluations: 630
Optimization terminated successfully.
         Current function value: 15.026831
         Iterations: 373
         Function evaluations: 602
Optimization terminated successfully.
         Current function value: 12.935445
         Iterations: 455
         Function evaluations: 716

Lets see how the y-center of our target's trace varies with column number. We can fit a polynomial to the trend and use this to estimate the y-center at any column.

In [41]:
# fit the relation with a polynomial
ind = 0 # which trace 0 or 1?
t_order = 3
trace_c = np.polyfit(cols, ycenter[:, ind], t_order)
fig, axarr = plt.subplots(2, sharex=True)
axarr[0].plot(cols, ycenter[:, ind], 'ro')
axarr[0].plot(xs, np.polyval(trace_c, xs), 'r')
axarr[0].axes.set_ylabel('y-coordinate')
axarr[1].plot(cols, ycenter[:, ind] - np.polyval(trace_c, cols), 'ro')
axarr[1].axes.set_ylim(-0.5, 0.5)
axarr[1].axes.set_ylabel('Fit Residual (pixels)')
plt.xlabel('Column Number');

The typical residuals should be a small fraction of a pixel.

Note we are fitting a Gaussian model to a non-Gaussian (galaxy) profile, so the center value may be systematically offset.

Now we can get the pixel offsets from the object trace for every pixel in the image. To properly account for non-linearities along the trace (i.e. the pixel scale may be 0.1 arcsec per y-pixel on the blue side but 0.12 arcsec per pixel on the red side) we could use the second object trace and the relative wavelength information, but here we will simply assume a simple linear relation.

In [42]:
# position offsets from the object trace (defined to be at slitpos = 0)
slitpos = yvals - np.polyval(trace_c, yvals)

Spatial profile of the target

Next, determine the spatial profile of the target and stretch this over the whole 2D image to make a weight map for later use.

The brightness of the target will vary as a function of wavelength, so to determine the profile we can normalize the counts in each wavelength bin using the counts at the trace center (or our approximation there of).

In [43]:
# subtract the sky
nosky = image - sky

# normalize to the pixel brightness at the trace center
yinds = (np.round(np.polyval(trace_c, xs))).astype(int)
normed = nosky / nosky[yinds, xs]

Fit a spline to the normalized profile

In [44]:
# get 1D arrays with the positions along the slit and the normalized counts
pos = slitpos.flatten()
counts = normed.flatten()

# sort by slit position
sort_inds = pos.argsort()
pos, counts = pos[sort_inds], counts[sort_inds]

# fit a spline to model the spatial profile
t = np.linspace(pos.min() + 2, pos.max() - 2, ny // 2) # spline knot points
profile_spl = LSQUnivariateSpline(pos, counts, t)

# remove outliers and re-fit
diff = counts - profile_spl(pos)
sample = sigma_clip(diff)
w = ((np.abs(diff) / sample.std()) < 5) & np.isfinite(diff)
profile_spl = LSQUnivariateSpline(pos[w], counts[w], t)
In [45]:
# plot the target profile
#plt.plot(pos, counts, 'r.')
plt.plot(pos, profile_spl(pos) )
plt.xlim(-40, 50)
#plt.ylim(-2, 4)
plt.grid();

Calculate the profile for every column to make a 2D profile image.

In [46]:
# create the profile image
profile_image = profile_spl(slitpos)

# de-weight negative values in provile_image
profile_image[profile_image < 0] = 0
In [47]:
# display the result in ds9
#ds9.set('frame 4')
#ds9.set_np2arr(profile_image);

Extract the target spectrum (basic version)

We have located the target trace, the sky background has been removed, so all that is left to do is sum up the counts in each wavelength bin in some aperture. We can pick any size wavelength bin we like, but at this phase it is common to use the native pixel scale to set the binning. If the target extends over many pixels in the spatial direction or the tilt/curvature of the lines of constant wavelength is large, you may have to take into account the changes in wavelength with row number. But for the supplied data the change in wavelength is not significant along the short extraction window, so we can simply sum along columns.

In [48]:
# select which rows to sum
w = (slitpos > -10) & (slitpos < 10)
ymin, ymax = yvals[w].min(), yvals[w].max()

# calculate the sum
spec_basic = nosky[ymin:ymax, :].sum(axis=0)

# sky background
skybg_basic = sky[ymin:ymax, :].sum(axis=0)

# plot the extracted spectrum
plt.plot(xs, spec_basic);

Optimal extraction

Analogous to PSF-fitting in 2D image data, we can weight the extraction using the spatial profile we calculated above to limit the noise from pixels that are dominated by the sky background. Instead of fitting a model to the data however, we will simply calculate a weighted average in each wavelength bin. Using our basic extraction above we can determine the bias between the sum and the weighted average and then use this to re-scale the weighted average, because ultimately we want the sum of the counts from the target, not the average in each wavelength bin.

In [49]:
# calculate the weighted average (for each column)
spec_opt = (nosky * profile_image)[ymin:ymax, :].sum(axis=0) / profile_image.sum(axis=0)

# calculate the bias factor needed to scale the average to a sum
bias_factor = np.median(spec_basic / spec_opt)
spec_opt *= bias_factor

# same for the sky background
skybg_opt = (sky * profile_image)[ymin:ymax, :].sum(axis=0) / profile_image.sum(axis=0)
bias_factor_sky = np.median(skybg_basic / skybg_opt)
skybg_opt *= bias_factor_sky

# plot the extracted spectrum
plt.plot(xs, spec_basic, label='basic extraction')
plt.plot(xs, spec_opt, label='optimal extraction')
plt.legend()
#plt.ylim(0, 2000);
#plt.yscale('log')
Out[49]:
<matplotlib.legend.Legend at 0x7f46700b1dd8>

If you look closely you should see that the basic extraction is noisier than the optimal.

Wavelength solution

Up to now we have been working with relative wavelengths. Now we will use the emission lines to calibrate these into absolute wavelengths.

There are several ways to do this. One method is to match the sky lines to a published atlas to determine the identity of several lines and then use the known wavelengths of these to determine the relative to absolute mapping function.

Instead of working with individual lines, we can instead find a mapping that matches our extracted sky spectrum to a reference spectrum with a calibrated wavelength axis. Ideally this reference should be matched in spectral resolution (and spectral response) to the data at hand. This module includes such a reference spectrum, which is a version of the UVES high-resolution night sky spectrum convolved to roughly match the resoultion of the Keck/LRIS data.

In [50]:
# re-scale the sky spectrum so that the max counts are around 1 for convience
bgspec = skybg_opt / skybg_opt.max()
In [51]:
# load the reference sky spectrum (conviently convolved to approximatley the spectral resulution of our data)
skyref = np.genfromtxt(os.path.join(data_dir,'UVES_nightsky_lowres.dat'), names='wav, flux')

Plot the reference and extracted sky spectra on two separate plots and adjust the wavelength range of the reference spectrum until the two roughly match. This will help you see the mapping between column numbers and wavelengths.

In [52]:
# compare the reference spectrum and the extracted sky spectrum
fig, axarr = plt.subplots(2)
axarr[0].plot(skyref['wav'], skyref['flux'])
axarr[0].axes.set_ylabel('Flux Density ($10^{16} f_{\lambda}$)')
axarr[0].axes.set_xlabel('Wavelength ($\AA$)')

# adjust the wavlength range to approx. match the extracted sky spectrum range
axarr[0].axes.set_xlim(5100, 10400)

# plot the extracted sky spectrum 
axarr[1].plot(xs, bgspec)
axarr[1].axes.set_ylabel('Counts')
axarr[1].axes.set_xlabel('Column Number');

By scaling the extracted counts to roughly match the reference spectrum and shifting/scaling the column number values to match the reference spectrum wavelengths, we can determine column number to wavelength mapping.

Let's define a function that takes our column numbers and counts spectra and returns a spectrum with wavelengths and fluxes that roughly match the reference spectrum. Note you usually match a model to the data, but here we go the other way.

If we fit the extracted spectrum over a limited wavelength range, we can approximate the mapping with a quadratic function.

In [53]:
def model_sky(params, dx, counts):
    wav0, a, b, scale, c0 = params
    
    dtype = []
    dtype.append(('wav', float))
    dtype.append(('flux', float))
    model = np.zeros(counts.shape, dtype=dtype)
    model['wav'] = wav0 + a * dx + b * dx**2
    model['flux'] = c0 + scale * counts
    
    return model
In [54]:
# roughly match column values to wavelength values
col = 2000
wav = 7730
params = (wav, 1.2, 0, 4, 0)
model = model_sky(params, xs - col, bgspec)

plt.plot(model['wav'], model['flux'], label='Extracted Sky Background')
plt.plot(skyref['wav'], skyref['flux'], label='Reference Sky Spectrum')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Counts or Flux')
plt.xlim(5000, 10200)
plt.legend();

You may find that over the full wavelength range the quadratic mapping is not sufficient. One way to handle this would be to adopt a higher order mapping function. Instead, we can use the quadratic function to find the mapping over a number of sections of limited wavelength range.

The procedure is similar to what has already been done above: make an initial guess for the model parameters and then minimize the difference between the scaled model and the reference spectrum.

In [55]:
# we must interpolate the reference spectrum to the model wavelengths
skyref_interp = interp1d(skyref['wav'], skyref['flux'], bounds_error=False)

def get_sky_difference(params, dx, specdata, skyatlas):
    model = model_sky(params, dx, specdata)
    
    # residual
    res = model['flux'] - skyref_interp(model['wav'])
    
    return np.sum(res**2 / np.sqrt(np.abs(model['flux'])))

Select sections of the extracted spectrum to fit. Since there is a lower density of lines in the blue (left) portion of the spectrum, we need larger bins to ensure we have sufficient information to fit.

In [56]:
cols = [200 * (i + 1) for i in range(8)]
cols.extend([max(cols) + 100 * (i + 1) for i in range(51)])
cols = np.array(cols)

hwidths = (cols[1:] - cols[:-1]) // 2
cols = cols[:-1]
w = (cols + hwidths) < nx
cols = cols[w]
hwidths = hwidths[w]
cols, hwidths
Out[56]:
(array([ 200,  400,  600,  800, 1000, 1200, 1400, 1600, 1700, 1800, 1900,
        2000, 2100, 2200, 2300, 2400, 2500, 2600, 2700, 2800, 2900, 3000,
        3100, 3200, 3300, 3400, 3500, 3600, 3700, 3800, 3900]),
 array([100, 100, 100, 100, 100, 100, 100,  50,  50,  50,  50,  50,  50,
         50,  50,  50,  50,  50,  50,  50,  50,  50,  50,  50,  50,  50,
         50,  50,  50,  50,  50]))

Set up an array to hold the best-fit parameters from each section.

In [57]:
nparams = 5
fitparams = np.zeros((cols.size, nparams)) + np.nan

Write a function to show the fitting results.

In [58]:
def plot_sky_model(skyref, model):
    plt.plot(model['wav'], model['flux'], label='Extracted Sky Background');
    plt.plot(skyref['wav'], skyref['flux'], label='Reference Sky Spectrum');
    plt.xlim(model['wav'].min(), model['wav'].max());
    plt.ylim(model['flux'].min() - 0.25, model['flux'].max() + 0.25)

    plt.xlabel('Wavelength ($\AA$)')
    plt.ylabel('Scaled Counts or Flux')
    plt.legend();

Start by guessing a solution for one of the middle sections

In [59]:
guesses = {}
guesses[0] = (5600, 1.2, 0, 40, 0)
guesses[3] = (6295, 1.2, 0, 40, 0)
guesses[11] = (7719, 1.2, 0, 4, 0)
guesses[21] = (8920, 1.2, 0, 4, 0)
In [60]:
ind = 11
col = cols[ind]
hwidth = hwidths[ind]
#guess = (7719, 1.2, 0, 4, 0)
guess = guesses[ind]
w = (xs > (col - hwidth)) & (xs < (col + hwidth))
model = model_sky(guess, xs[w] - col, bgspec[w])
plot_sky_model(skyref, model)