Learn.Astropy / Tutorials

Download

Outline

1. Authors
2. Learning Goals
3. Keywords
4. Summary

Spectroscopic Data Reduction Part 1: Tracing

Authors

Adam Ginsburg, Kelle Cruz, Lia Corrales, Jonathan Sick, Adrian Price-Whelan

Learning Goals

• Open a two-dimensional spectrum from an image file (bitmap)
• Fit a spectroscopic trace

Keywords

Spectroscopy

Summary

This tutorial will walk through the derivation of a spectroscopic trace model and extraction using astropy tools.

A spectroscopic trace is the path of a point source (star) spectrum through a two-dimensional dispersed spectrum. The trace is needed because, in general, spectra are not perfectly aligned with the axes of a detector.

with open('requirements.txt') as f:
    print(f"Required packages for this notebook:\n{f.read()}")

Required packages for this notebook:
astropy
astroquery>=0.4.8.dev9474  # 2024-09-24 pinned for Gaia column capitalization issue
IPython
numpy
pillow
matplotlib

Step 1: Examine the spectrum

We'll work with a 2D spectrum that contains no attached metadata, so we have to infer many of the features ourselves.

All we know is that this is a spectrum of a star, Aldebaran.

Our data are in the form of .bmp (bitmap) files, so we need PIL (Python Imaging Library) to open them.

While .bmp files are not astronomical standard FITS files, as are commonly delivered from professional observatories, image formats like .bmp, .jpg, .raw, .png, etc. produced by consumer cameras may also be used for spectroscopy.

In this case, our images are monochromatic, which is similar to standard FITS images.

from PIL import Image
import numpy as np

import matplotlib.pyplot as plt
plt.rcParams['image.origin'] = 'lower'
plt.style.use('dark_background') # Optional!

spectrum_filename = "aldebaran_3s_1.bmp"

image_data = Image.open(spectrum_filename)

# because this is an image, simply entering on the command line will show it
image_data

# our data are unsigned 8-bit integers (0-255) representing a monochromatic image
# we can see this by printing the array version of the image
# we can also see its shape, verifying that it is indeed 2-dimensional
image_array = np.array(image_data)
image_array, image_array.shape

(array([[41, 43, 50, ..., 52, 28, 27],
        [33, 22, 48, ..., 37, 35, 26],
        [41, 64, 30, ..., 21, 30, 33],
        ...,
        [33, 28, 27, ..., 26, 33, 28],
        [35, 25, 23, ..., 23, 46, 46],
        [32, 35, 22, ..., 23, 40, 29]], dtype=uint8),
 (1200, 1600))

# but we'd like to see it with axes labeled
plt.imshow(image_data, cmap='gray')
plt.colorbar(); # the semicolon at the end of the last line prevents ipython from printing out the object

The main goal of the trace is to obtain a model f(x) defining the vertical position of the light (the signal) along the detector.

We're going to start by assuming that wavelength dispersion is in the X-direction and the Y-direction is entirely spatial.

This is an approximation made by inspecting the image by eye.

Step 2a. Try to find the spine to trace using argmax

To obtain the trace, we first measure the Y-value at each X-value. we'll start with the trivial approach of using argmax:

yvals = np.argmax(image_data, axis=0)
xvals = np.arange(image_data.width)

plt.plot(xvals, yvals, 'x')
plt.ylabel("Argmax trace data")
plt.xlabel("X position");

There's a pretty clear line going through the center, which represents our signal, but there are also a lot of erroneous data points.

We can get rid of most of the bad data just by filtering it out using a pixel mask

bad_pixels = (yvals < 400) | (yvals > 500)

plt.plot(xvals, yvals, 'x')
plt.plot(xvals[bad_pixels], yvals[bad_pixels], 'rx')
plt.ylabel("Argmax trace data")
plt.xlabel("X position");

plt.plot(xvals[~bad_pixels], yvals[~bad_pixels], 'x')
plt.ylabel("Argmax trace data")
plt.xlabel("X position");

We can be a little more precise by 'zooming in' along the y-axis, so we refine the mask again to be over a narrower range:

bad_pixels = (yvals < 425) | (yvals > 460)

plt.plot(xvals[~bad_pixels], yvals[~bad_pixels], 'x')
plt.ylabel("Argmax trace data")
plt.xlabel("X position");

The stuff at x>1100 looks bad, but there's still signal out there.

We can see there is clear signal out to nearly pixel ~1400:

plt.imshow(image_array[425:475,:])
plt.gca().set_aspect(10)

Step 2b: Use moment analysis to extract a spine to trace

We can use moments) to provide a different, possibly better, estimate of where the trace's center is. The advantage of moment analysis is that we're using all of the data to estimate the vertical position, not just the single brightest value, which is what we used above.

Note that we need to subtract off the background to avoid a bias toward the center, so we use the median of the whole image as our background estimate.

(the first-order moment is the intensity-weighted mean position, $$m_1 = \frac{\Sigma_i x_i f(x_i)}{\Sigma_i f(x_i)}$$ where $x_i$ is the position and $f({x_i})$ is the intensity at that position. $f(x_i)$ must be zero in the signal-free region for $m_1$ to return an accurate estimate of the location of the peak)

# we use a cutout around the traced line, so the Y-values are from that cutout
# the `repeat` command here is used to extend our Y-axis position values, which are 425, 426, ... 475
# along the X-direction.  The indexing with [:, None] adds a "dummy" axis along the second (x) dimension,
# then `repeat` copies our Y-axis values.  The resulting array has the same shape as our weight array,
# which is image_array[425:475, :] minus the median
yaxis = np.repeat(np.arange(425, 475)[:,None],
                  image_array.shape[1], axis=1)
background = np.median(image_array)
# moment 1 is the data-weighted average of the Y-axis coordinates
weighted_yaxis_values = np.average(yaxis, axis=0,
                                   weights=image_array[425:475,:] - background)

plt.plot(xvals, weighted_yaxis_values, 'x')
plt.xlabel("X Coordinate")
plt.ylabel("Moment-1 estimated Y-value trace");

Overplot the "weighted", centroid locations on the data to verify they look reasonable.

# we need to use the 'extent' keyword to have the axes correctly labeled
plt.imshow(image_array[425:475,:],
          extent=[0,image_array.shape[1],425,475])
plt.gca().set_aspect(10) # we stretch the image out by 10x in the y-direction
plt.plot(xvals, weighted_yaxis_values, 'wx', alpha=0.5)
plt.axis((0,1600,425,475));

We can also compare the argmax and weighted approaches. They agree well at x<1200, but there are simply more points from the weighted approach at x>1200.

plt.plot(xvals, weighted_yaxis_values, 'x', label="Weighted", alpha=0.5)
plt.plot(xvals[~bad_pixels], yvals[~bad_pixels], '+', label="Argmax", alpha=0.5)
plt.legend(loc='best');

That's a decent set of data, we'll use the moments instead of the argmax. There's still some data to flag out, though:

bad_moments = (weighted_yaxis_values > 460) | (weighted_yaxis_values < 430)

plt.plot(xvals[~bad_moments], weighted_yaxis_values[~bad_moments], 'x', label="Weighted", alpha=0.5)
plt.plot(xvals[~bad_pixels], yvals[~bad_pixels], '+', label="Argmax", alpha=0.5)
plt.legend(loc='best');

Step 3. Fit the trace profile

We want a model f(x) that gives the y-value of the centroid as a function of x.

from astropy.modeling.polynomial import Polynomial1D
from astropy.modeling.fitting import LinearLSQFitter

# We fit a 2nd-order polynomial
polymodel = Polynomial1D(degree=2)
linfitter = LinearLSQFitter()
fitted_polymodel = linfitter(polymodel, xvals[~bad_moments], weighted_yaxis_values[~bad_moments])

fitted_polymodel

<Polynomial1D(2, c0=454.80138106, c1=-0.02381282, c2=0.00001121)>

plt.plot(xvals[~bad_moments], weighted_yaxis_values[~bad_moments], 'x', alpha=0.5)
plt.plot(xvals, fitted_polymodel(xvals), color='r');

We plot and examine the residuals to visually inspect whether the fit is good:

plt.plot(xvals[~bad_moments],
        weighted_yaxis_values[~bad_moments] - fitted_polymodel(xvals[~bad_moments]), 'x')
plt.ylabel("Residual (data-model)");

The curvature seen at the left is a sign of a suboptimal fit. Specifically, curvature in the residual indicates that we need to use a higher order model - i.e., we need more terms in the polynomial. We change degree=2 to degree=3.

polymodel = Polynomial1D(degree=3)
fitted_polymodel = linfitter(polymodel, xvals[~bad_moments], weighted_yaxis_values[~bad_moments])
fitted_polymodel

<Polynomial1D(3, c0=453.31622273, c1=-0.01256666, c2=-0.00000651, c3=0.00000001)>

plt.plot(xvals[~bad_moments], weighted_yaxis_values[~bad_moments], 'x', alpha=0.5)
plt.plot(xvals, fitted_polymodel(xvals), color='r');

plt.plot(xvals[~bad_moments],
      weighted_yaxis_values[~bad_moments] - fitted_polymodel(xvals[~bad_moments]), 'x')
plt.ylabel("Residual (data-model)");

Arguably, we should toss out the data at >1400 pixels since there's no clear signal there. We'll come back to this...

plt.imshow(image_array[425:475,:],
          extent=[0,image_array.shape[1],425,475])
plt.gca().set_aspect(10);

Again, we should verify the trace by overplotting on the original data:

plt.imshow(image_array[425:475,:], extent=[0,image_array.shape[1],425,475])
plt.gca().set_aspect(10)
plt.plot(xvals, fitted_polymodel(xvals), 'w')
plt.axis((0,1600,425,475));

Seeing the curve up in the model to the right - which we do not observe in the data - suggests we should re-fit without including the x>1200 data at all:

polymodel = Polynomial1D(degree=3)
fitted_polymodel = linfitter(polymodel, xvals[(~bad_moments) & (xvals < 1200)],
                             weighted_yaxis_values[(~bad_moments) & (xvals < 1200)])
fitted_polymodel

<Polynomial1D(3, c0=453.63074008, c1=-0.01596396, c2=0.00000083, c3=0.)>

We now have a satisfactory fit:

plt.imshow(image_array[425:475,:], extent=[0,image_array.shape[1],425,475])
plt.gca().set_aspect(10)
plt.plot(xvals, fitted_polymodel(xvals), 'w')
plt.axis((0,1600,425,475));

plt.plot(xvals[~bad_moments], weighted_yaxis_values[~bad_moments], 'x', alpha=0.5)
plt.plot(xvals, fitted_polymodel(xvals), color='r');

plt.plot(xvals[~bad_moments & (xvals < 1200)],
        weighted_yaxis_values[~bad_moments & (xvals < 1200)] - fitted_polymodel(xvals[~bad_moments & (xvals < 1200)]), 'x')
plt.plot(xvals[~bad_moments & (xvals > 1200)],
        weighted_yaxis_values[~bad_moments & (xvals > 1200)] - fitted_polymodel(xvals[~bad_moments & (xvals > 1200)]), 'r+', alpha=0.5)
plt.ylim(-5, 5)
plt.ylabel("Residual (data-model)");

Step 4. Obtain a trace profile

Now we can extract the data along that trace.

We want to take a "profile" of the trace to see how many pixels on either side of the line we should include.

plt.imshow(image_array[425:475,:], extent=[0,image_array.shape[1],425,475])
plt.gca().set_aspect(10)
plt.fill_between(xvals, fitted_polymodel(xvals)-15,
                fitted_polymodel(xvals)+15,
                color='orange', alpha=0.5)
plt.axis((0,1600,425,475));

# start by taking +/- 15 pixels
npixels_to_cut = 15
trace_center = fitted_polymodel(xvals)
cutouts = np.array([image_array[int(yval)-npixels_to_cut:int(yval)+npixels_to_cut, ii]
                    for yval, ii in zip(trace_center, xvals)])
cutouts.shape

(1600, 30)

That last step deserves some explanation:

cutouts = np.array([image_array[int(yval)-npixels_to_cut:int(yval)+npixels_to_cut, ii]
                    for yval, ii in zip(trace_center, xvals)])

• [... for yval, ii in zip(trace_center, xvals)] takes each trace y-value and each x-value and 'zips' them together, so each iteration of the for loop has one x, y pair
• image_array[int(yval)-npixels_to_cut:int(yval)+npixels_to_cut, ii] is taking a single pixel along the x-direction (the second dimension, ii) and a range of pixels along the y-direction, i.e., y+/-n
• these are put together in a loop, so we have a y+/-n pixel region for each x-pixel
• then we make them all into an array

We can see the result visually:

ax1 = plt.subplot(2,1,1)
ax1.imshow(image_array[425:475,:], extent=[0,image_array.shape[1],425,475])
ax1.set_aspect(10)
ax1.set_title("We go from this...")
ax2 = plt.subplot(2,1,2)
ax2.imshow(cutouts.T)
ax2.set_title("...to this")
ax2.set_aspect(10)

Then we average along the X-direction to get the trace profile:

mean_trace_profile = (cutouts - background).mean(axis=0)
trace_profile_xaxis = np.arange(-npixels_to_cut, npixels_to_cut)
plt.plot(trace_profile_xaxis, mean_trace_profile)
plt.xlabel("Distance from center")
plt.ylabel("Average source profile");

We want to fit that profile with a Gaussian for future use, so we import the Gaussian model profile and non-linear fitter and run a fit:

from astropy.modeling.models import Gaussian1D
from astropy.modeling.fitting import LevMarLSQFitter
lmfitter = LevMarLSQFitter()
guess = Gaussian1D(amplitude=mean_trace_profile.max(), mean=0, stddev=5)
fitted_trace_profile = lmfitter(model=guess, x=trace_profile_xaxis, y=mean_trace_profile)
model_trace_profile = fitted_trace_profile(trace_profile_xaxis)
fitted_trace_profile

<Gaussian1D(amplitude=123.84846797, mean=0.17719819, stddev=5.10872134)>

plt.plot(trace_profile_xaxis, mean_trace_profile, label='data')
plt.plot(trace_profile_xaxis, model_trace_profile, label='model')
plt.legend(loc='best')
plt.xlabel("Distance from center")
plt.ylabel("Average source profile");

Both the empirical trace profile mean_trace_profile and the modeled model_trace_profile can reasonably be used; the latter is more convenient to serialize (i.e., write to disk or on paper)

Step 5. Extract the traced spectrum

We can obtain our spectrum by directly averaging the pixels along the trace:

average_spectrum = (cutouts - background).mean(axis=1)
plt.plot(average_spectrum);

Or, we can obtain our spectrum by taking the trace-weighted average:

trace_avg_spectrum = np.array([np.average(
    image_array[int(yval)-npixels_to_cut:int(yval)+npixels_to_cut, ii] - background,
    weights=mean_trace_profile)
                               for yval, ii in zip(trace_center, xvals)])

We can also do this with the Gaussian weights:

gaussian_trace_avg_spectrum = np.array([np.average(
    image_array[int(yval)-npixels_to_cut:int(yval)+npixels_to_cut, ii] - background,
    weights=model_trace_profile)
                               for yval, ii in zip(trace_center, xvals)])

plt.plot(average_spectrum, label="Direct Average")
plt.plot(trace_avg_spectrum, label='Trace-weighted average')
plt.plot(gaussian_trace_avg_spectrum, label='Gaussian-model-Trace-weighted average', alpha=0.5, linewidth=0.5, color='r')
plt.legend(loc='best');

In general, the trace-weighted average will have higher signal-to-noise, as seen here (while we haven't measured the noise, it is approximately constant across the image).

Note that the Gaussian model and the direct trace yield nearly identical results

Step 6: Repeat for another star

In this last step, we go through all the above steps again for another star (Deneb), but with less explanation.

image_array_2 = np.array(Image.open('deneb_3s_13.63g_1.bmp'))
plt.imshow(image_array_2, cmap='gray')
plt.colorbar();

plt.imshow(image_array_2[470:520,:], extent=[0,1600,470,520])
plt.gca().set_aspect(10)

yaxis2 = np.repeat(np.arange(470, 520)[:,None], image_array_2.shape[1], axis=1)
weighted_yaxis_values2 = np.average(yaxis2, axis=0, weights=image_array_2[470:520,:] - np.median(image_array_2))
polymodel2 = Polynomial1D(degree=3)
fitted_polymodel2 = linfitter(polymodel2, xvals, weighted_yaxis_values2)
trace_center2 = fitted_polymodel2(xvals)

plt.plot(xvals, weighted_yaxis_values2, 'x', alpha=0.5)
plt.plot(xvals, trace_center2, color='r');

plt.imshow(image_array_2[470:520,:], extent=[0,1600,470,520])
plt.plot(xvals, weighted_yaxis_values2, 'w+', alpha=0.25)
plt.plot(xvals, trace_center2, color='r')
plt.gca().set_aspect(10)

spectrum2 = np.array([np.average(image_array_2[int(yval)-npixels_to_cut:int(yval)+npixels_to_cut, ii] - np.median(image_array_2),
                                weights=mean_trace_profile)
                     for yval, ii in zip(trace_center2, xvals)])

plt.plot(spectrum2);

In the next tutorial, Spectroscopic Data Reduction 2, we'll work on the wavelength calibration.

Astropy footer