Toward Long-Slit Data Reduction in Python

I recently collected a bunch of long-slit spectroscopy data at the Apache Point Observatory 3.5m telescope. Now, it's been about 20 yrs since I've reduced a long-slit spectrum. I wasn't keen to re-learn IRAF, and so I thought I'd give it a try in Python. I've managed to get as far as correcting the distortion of the wavelength grid based on observations of neon arc lines. The correction is subtle: essentially, in the raw image, an emission line gently curves so that an image column does not map uniquely to a specific wavelength. The correction removes this curvature.

This wavelength grid correction makes sky subtraction a lot easier, and my first results are shown in the figure below. The top panel is a raw spectrum of an active galaxy; actually, if you look closely, you'll see spectra of two galaxies whose nuclei fell within the slit. The spectrum is dominated by sky glow (bright vertical lines). The bottom panel shows the image after the wavelength grid correction, and the sky glow has been subtracted by interpolating over the source. The quality of the sky subtraction has been affected by cosmic ray hits, but it's a good start. To be clear: all work was done in Python with no help from IRAF!