{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"
Guide to Long-Slit Spectral Extractions
\n",
" by Robert Quimby (SDSU)
\n",
"\n",
"\n",
"This Jupyter notebook supplements the spectroscopy lesson from the GROWTH Summer School 2020 and serves as a complete guide to long-slit spectral extractions. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-school-2020.html"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" NOTE: this version of the notebook is incomplete. With this version the reader can test their understanding as they progress through the notebook by completing the missing code statements. The intent here is not to supply a push-button solution but to teach the art of spectroscopic reductions so that the reader can adapt the techniques as necessary to their own data. However, a separate solutions guide is also available. \n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Process\n",
"\n",
"- [Step 0](#step0) - Image Processing\n",
"- [Step 1](#step1) - Find the (Initial) 1-D Wavelength Solution\n",
"- [Step 2](#step2) - Find the 2-D Wavelength Solution\n",
"- [Step 3](#step3) - Flat Field the 2D Spectra (Optional)\n",
"- [Step 4](#step4) - Model the 2D Sky\n",
"- [Step 5](#step5) - Extract the 1-D Target Spectrum\n",
"- [Step 6](#step6) - Tweak the Wavelength Solution Using Sky Lines\n",
"- [Step 7](#step7) - Flux Calibration (Optional)\n",
"- [Step 8](#step8) - Removal of Telluric Features (Optional) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Required Python modules\n",
"* python 3\n",
"* `astropy`\n",
"* `numpy`\n",
"* `scipy`\n",
"* `matplotlib`"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# python modules that we will use\n",
"import os\n",
"from glob import glob\n",
"import numpy as np\n",
"from astropy.io import fits\n",
"from astropy.stats import sigma_clip\n",
"from scipy.interpolate import interp1d, LSQUnivariateSpline, LSQBivariateSpline\n",
"from scipy.optimize import fmin, curve_fit\n",
"from scipy.signal import find_peaks\n",
"\n",
"%matplotlib inline\n",
"import matplotlib.pylab as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from platform import python_version\n",
"assert python_version() >= '3.6', \"Python version 3.6 or greater is required for f-strings\"\n",
"\n",
"import scipy\n",
"assert scipy.__version__ >= '1.4', \"scipy version 1.4 or higher is required for this module\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# change plotting defaults\n",
"plt.rc('axes', labelsize=14)\n",
"plt.rc('axes', labelweight='bold')\n",
"plt.rc('axes', titlesize=16)\n",
"plt.rc('axes', titleweight='bold')\n",
"plt.rc('font', family='sans-serif')\n",
"plt.rcParams['figure.figsize'] = (17, 7)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# where are the data located?\n",
"cwd = os.getcwd()\n",
"data_dir = os.path.join(cwd, 'data')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" \n",
" Step 0: Image Processing\n",
"
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The 2-D spectral images must be processed to remove bias. Just as with other CCD imaging data, overscan should be removed and (for best results) the residual 2-D bias should be subtracted.\n",
"\n",
"For display purposes it is also convienent to transpose the images (if necessary) so that wavelength increases to the right, but this is not a requirement. \n",
"\n",
"The example data provided have already been overscan subtracted and transposed, so if you are starting with these, no actions are required for this step."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's begin by looking at the processed 2D image of the science spectrum, `spec_sci.fits`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# load the 2D science image data\n",
"image = fits.getdata(os.path.join(data_dir,'spec_sci.fits'))\n",
"\n",
"# get the a pixel coordinate near the image center (for later use)\n",
"ny, nx = image.shape\n",
"cy, cx = ny//2, nx//2\n",
"\n",
"# create 1d arays of the possible x and y values (for later use)\n",
"xs = np.arange(nx)\n",
"ys = np.arange(ny)\n",
"\n",
"# pixel coordinates for each pixel (for later use)\n",
"yvals, xvals = np.indices(image.shape)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# this is a helper function to display images in this notebook\n",
"def show_image(image, lower=-1, upper=3, extent=None):\n",
" sample = sigma_clip(image)\n",
" vmin = sample.mean() + lower * sample.std()\n",
" vmax = sample.mean() + upper * sample.std()\n",
" plt.figure(figsize=(15, 3))\n",
" plt.imshow(image, origin='lower', cmap='gray', aspect='auto', vmin=vmin, vmax=vmax, extent=extent)\n",
" plt.xlabel('Column Number')\n",
" plt.ylabel('Row Number');"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# display the image\n",
"show_image(image)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" \n",
" Step 1: Find the (Initial) 1-D Wavelength Solution\n",
"
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Recall that when you extract photometry from images you use different flat-field images to process your science frames depending on the filter used (e.g., you use the r-band flat-field to process r-band science frames). The same idea holds for spectroscopy except each 2-D spectral image contains a range of wavelengths. Thus it helps to first (roughly) identify the effective wavelength of each pixel before you start the extraction process. Because you may have thousands of unique wavelength bands on each 2-D spectral image, this can be the most time consuming step in the whole process. But once you set the pixel wavelengths, you will have reduced the task of spectral extraction to simple photometric extraction (but on an industrial scale). \n",
"\n",
"To begin, it helps to use arc-lamp spectra, which produce a set of emission lines of known wavelengths. By matching up the features you observed to an atlas of the expected features (with known wavelengths), you can create a mapping between pixel column and row numbers to wavelengths in angstroms."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# load the 2D lamp image data\n",
"lamp_image = fits.getdata(os.path.join(data_dir,'spec_lamp.fits'))\n",
"show_image(lamp_image)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You should see that the lamp spectra consist of emission lines from the excited gas in the arc lamps (in this case mostly a combination of Neon and Argon). The wavelengths of these lines have been precisely measured in laboratory experiments, and we can use these known wavelengths and the columns where they are found to construct our wavelength solution. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load the comparison line list"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# load the \"linelist\" with precise laboratory wavelengths for these lamp lines\n",
"dtype = [('wav', float), ('id', 'U2')]\n",
"linelist = np.genfromtxt(os.path.join(data_dir, 'line_list.dat'), dtype=dtype)\n",
"linelist.sort(order='wav')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Extract the lamp spectra from a few rows near the middle"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# simple average of each column over a few rows\n",
"row1 = cy - 5\n",
"row2 = cy + 5\n",
"lamp_spec = lamp_image[row1:row2, :].mean(axis=???)\n",
"plt.plot(lamp_spec, c='g')\n",
"plt.xlabel('Column Number')\n",
"plt.ylabel('Counts');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that many of the lines peak at about 65000 counts? This indicates that these lines are saturated (too bright) and will not be useful for determining the wavelength solution."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Identify each column where the flux peaks"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We need the column value where each line peaks for our column to wavelength mapping. To start we can get approximate column values using `scipy.signal.find_peaks`. The keyword `prominence` is the desired minimum contrast (in counts) between the value at a peak and the value at the nearest troughs. Set this so that significant lines are kept but noise peaks are rejected. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"# locate approximate line centers\n",
"prominence = ???\n",
"peak_cols, properties = find_peaks(lamp_spec, prominence=prominence)\n",
"\n",
"# plot the spectra\n",
"plt.plot(lamp_spec, c='g')\n",
"\n",
"# mark the peaks\n",
"plt.scatter(peak_cols, lamp_spec[peak_cols], marker='x', c='k')\n",
"plt.xlabel('Column Number')\n",
"plt.ylabel('Counts')\n",
"plt.yscale('log')\n",
"plt.grid();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `find_peaks` function merely returns integer indices into the `lamp_spec` array. We can fit Gaussian models to each one of these peaks to find the line centers to much better than 1 column. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# use a Gaussian function to model spectral lines\n",
"def gaussian(x, *params):\n",
" amp, x0, sigma = params\n",
" return amp * np.exp(-(x - x0)**2 / 2 / sigma**2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# define a function to determine the precise column centers for each lamp line\n",
"def get_lamp_lines(lamp_spec, prominence=100):\n",
" peak_cols, properties = find_peaks(lamp_spec, prominence=prominence)\n",
"\n",
" # record in a structured array for convenience \n",
" dtype = []\n",
" dtype.append( ('col', float) )\n",
" dtype.append( ('counts', float) )\n",
" dtype.append( ('x', float) )\n",
" dtype.append( ('y', float) )\n",
" dtype.append( ('sigma', float) )\n",
" dtype.append( ('id', 'U2') )\n",
" dtype.append( ('wav', float) )\n",
" dtype.append( ('wavres', float) )\n",
" dtype.append( ('used', bool) )\n",
" lamp_lines = np.zeros(peak_cols.size, dtype=dtype)\n",
" lamp_lines['col'] = peak_cols\n",
" lamp_lines['counts'] = lamp_spec[peak_cols]\n",
" lamp_lines['x'] = np.nan\n",
" lamp_lines['y'] = np.nan\n",
" lamp_lines['sigma'] = np.nan\n",
" lamp_lines['wav'] = np.nan\n",
" lamp_lines['wavres'] = np.nan\n",
"\n",
" # fit each peak to determine precise center\n",
" cols = np.arange(lamp_spec.size)\n",
" sigma_guess = 2.5\n",
" for line in lamp_lines:\n",
" if line['counts'] > 60000:\n",
" # line is saturated...skip\n",
" continue\n",
"\n",
" i0 = max([0, int(line['col'] - 5)])\n",
" i1 = min([lamp_spec.size - 1, int(line['col'] + 5)])\n",
" guess = (line['counts'], line['col'], sigma_guess)\n",
" bounds = ((0, line['col'] - 3, 0), (np.inf, line['col'] + 3, np.inf))\n",
" try:\n",
" popt, pcov = curve_fit(gaussian, cols[i0:i1], lamp_spec[i0:i1], p0=guess, bounds=bounds)\n",
" except RuntimeError:\n",
" # curve_fit failed to converge...skip\n",
" continue\n",
"\n",
" line['x'] = popt[1]\n",
" line['y'] = gaussian(popt[1], *popt)\n",
" line['sigma'] = popt[2]\n",
"\n",
" # filter lamp_lines to keep only lines that were fit\n",
" wasfit = np.isfinite(lamp_lines['x'])\n",
" lamp_lines = lamp_lines[wasfit]\n",
" print('found center pixel values for', lamp_lines.size, 'lines')\n",
" return lamp_lines"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# helper function to mark lamp lines with and 'x'\n",
"def mark_peaks(plt, lamp_lines, xtype='x', ytype='y', c='k'):\n",
" w = np.isfinite(lamp_lines['wav'])\n",
" if w.sum() > 0:\n",
" plt.scatter(lamp_lines[xtype][w], lamp_lines[ytype][w], c=np.abs(lamp_lines['wavres'][w]), zorder=10)\n",
" if (~w).sum() > 0:\n",
" plt.scatter(lamp_lines[xtype][~w], lamp_lines[ytype][~w], c=c, marker='x')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# find the precise column centers for all lamp lines\n",
"lamp_lines = get_lamp_lines(lamp_spec, prominence=prominence)\n",
"\n",
"# plot the lamp spectra\n",
"plt.plot(lamp_spec, c='g')\n",
"plt.xlabel('Column Number')\n",
"plt.ylabel('Counts')\n",
"plt.grid()\n",
"\n",
"# mark initial values from find_peaks in red\n",
"mark_peaks(plt, lamp_lines, 'col', 'counts', c='r')\n",
"\n",
"# mark best-fit values\n",
"mark_peaks(plt, lamp_lines)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load reference spectrum"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have precise column centers for each line, we need to identify each line to set its wavelength. We can begin by looking for similar line features in a wavelength calibrated lamp spectrum atlas. \n",
"\n",
"For the example Keck/LRIS spectra provided, a reference lamp spectrum is available in the `lamp_reference.dat` file. This may also work for other low-resolution spectrographs using similar arc lamps (primarily Ne, Ar). "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# load and plot the provided spectral atlas\n",
"lamp_ref = np.genfromtxt(os.path.join(data_dir, 'lamp_reference.dat'), names='wav, counts')\n",
"plt.plot(lamp_ref['wav'], lamp_ref['counts'], c='r')\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('Counts');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define a rough mapping between column number and wavelength"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# this function will be used below to match a given line to a given list\n",
"def match_to_list(listing, values, plt=None, tol=None, revcoeff=None, c='k'):\n",
" matched = []\n",
" cols = []\n",
" for value in values:\n",
" absdiff = np.abs(value - listing)\n",
" ind = np.argmin(absdiff)\n",
" if tol is None:\n",
" bestmatch = listing[ind]\n",
" elif absdiff[ind] < tol:\n",
" bestmatch = listing[ind]\n",
" else:\n",
" bestmatch = np.nan\n",
" matched.append(bestmatch)\n",
"\n",
" if plt is not None:\n",
" plt.axvline(bestmatch, ls='dotted', c=c)\n",
" \n",
" if revcoeff is not None:\n",
" col = np.polyval(revcoeff, bestmatch)\n",
" cols.append(col)\n",
" print(f\"{bestmatch:.1f} is expected near column {col:.0f}\")\n",
"\n",
" if revcoeff is not None:\n",
" return np.array(matched), cols\n",
" \n",
" return np.array(matched)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Focus on a small wavelength range of the reference spectrum (between `wav1` and `wav2`). Pick a few good (non-blended) lines to start with."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"wav1 = 5700\n",
"wav2 = 7500\n",
"\n",
"# plot the reference spectrum in red\n",
"plt.plot(lamp_ref['wav'], lamp_ref['counts'], label='reference', c='r')\n",
"plt.xlim(wav1, wav2)\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('Counts'); plt.grid();\n",
"\n",
"# mark all wavelengths available in linelist\n",
"for row in linelist:\n",
" plt.axvline(row['wav'], ls='dashed', c='0.8')\n",
"\n",
"# pick a few lines in this plot\n",
"rough_waves = [] # <---- add a few values here\n",
"refwavs = match_to_list(linelist['wav'], rough_waves, plt=plt)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now try to find a section of your `lamp_spec` (between columns `col1` and `col2`) that looks similar."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# find a section of lamp_spec that looks similar\n",
"col1 = 300\n",
"col2 = 1800\n",
"plt.plot(lamp_spec, c='g')\n",
"plt.xlim(col1, col2)\n",
"plt.xlabel('Column Number'); plt.ylabel('Counts'); plt.grid()\n",
"\n",
"# mark lines with Gaussian-fit centers\n",
"mark_peaks(plt, lamp_lines)\n",
"\n",
"# record the rough column numbers of the same lines as above in the same order\n",
"rough_cols = [] # <---- add rough column values here\n",
"refcols = match_to_list(lamp_lines['x'], rough_cols, plt=plt)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When you have marked 3-4 matching lines, use the cells below to record the line wavelengths in `lamp_lines`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def set_line_identity(lamp_lines, linelist, x, wav):\n",
" # find the closest matching lamp_line\n",
" ilamp = np.argmin(np.abs(lamp_lines['x'] - x))\n",
" \n",
" # find the closest matching row in the linelist\n",
" ilist = np.argmin(np.abs(linelist['wav'] - wav))\n",
" \n",
" # reset values in lamp_lines\n",
" lamp_lines[ilamp]['id'] = linelist[ilist]['id']\n",
" lamp_lines[ilamp]['wav'] = linelist[ilist]['wav']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# record ids and wavelengths of matched lines\n",
"for col, wav in zip(refcols, refwavs):\n",
" set_line_identity(lamp_lines, linelist, col, wav)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Set the initial wavelength solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the matched column/wavelength pairs we can define an initial wavelength solution by assuming a polynomial (linear to start) relation."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# this routine finds the polynomial coefficients needed to transform between \n",
"# column numbers and wavelengths (and vice versa). Outlier rejection is included.\n",
"def get_wavelength_solution(lamp_lines, order=4):\n",
" wfit = np.isfinite(lamp_lines['wav'])\n",
" \n",
" # define the reverse mapping (wavelength to column)\n",
" revcoeff = np.polyfit(lamp_lines['wav'][wfit], lamp_lines['x'][wfit], order)\n",
"\n",
" # define the forward mapping (column to wavelength)\n",
" coeff = np.polyfit(lamp_lines['x'][wfit], lamp_lines['wav'][wfit], order)\n",
" \n",
" # check the fit for outliers\n",
" fit_wav = np.polyval(coeff, lamp_lines['x'])\n",
" wavres = fit_wav - lamp_lines['wav']\n",
" lamp_lines['wavres'] = wavres\n",
" sample = wavres[wfit]\n",
" sample.sort()\n",
" sample = sample[int(0.1 * sample.size) : int(0.9 * sample.size + 0.5)] \n",
" std = np.std(sample, ddof=1)\n",
" w = wfit \n",
" w[wfit] = (np.abs(lamp_lines['wavres'][wfit]) < (5 * std))\n",
" if w.sum() != lamp_lines.size:\n",
" # re-fit with outliers rejected\n",
" coeff, revcoeff = get_wavelength_solution(lamp_lines[w], order=order)\n",
" \n",
" # reset wavelength residuals using new coefficients\n",
" fit_wav = np.polyval(coeff, lamp_lines['x'])\n",
" wavres = fit_wav - lamp_lines['wav']\n",
" lamp_lines['wavres'] = wavres\n",
" \n",
" lamp_lines['used'] = w\n",
" return coeff, revcoeff\n",
" \n",
"def check_wavelength_solution(lamp_spec, lamp_lines, coeff): \n",
" wavs = col_to_wav(coeff, np.arange(lamp_spec.size))\n",
" plt.plot(wavs, lamp_spec, c='g', lw=2)\n",
" mark_peaks(plt, lamp_lines, 'wav')\n",
" plt.colorbar(label='Residual ($\\AA$)')\n",
" plt.xlabel('Wavelength ($\\AA$)')\n",
" plt.ylabel('Counts')\n",
" plt.grid()\n",
" \n",
"def col_to_wav(coeff, cols):\n",
" return np.polyval(coeff, cols)\n",
"\n",
"def wav_to_col(revcoeff, wavs):\n",
" return np.polyval(revcoeff, wavs)\n",
"\n",
"def mark_matched(lamp_lines):\n",
" for line in lamp_lines:\n",
" plt.axvline(line['wav'], ls='dotted', c='k')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# estimate a linear relation between column number and wavlength\n",
"coeff, revcoeff = get_wavelength_solution(lamp_lines, order=???)\n",
"\n",
"# apply the wavelength solution to the column numbers\n",
"wavs = col_to_wav(???)\n",
"\n",
"# plot the initial wavelength calibrated spectrum\n",
"plt.plot(wavs, lamp_spec, c='g', lw=2, label='lamps')\n",
"plt.xlim(wav1, wav2)\n",
"\n",
"# plot the reference spectrum in red\n",
"plt.plot(lamp_ref['wav'], lamp_ref['counts'], label='reference', c='r')\n",
"plt.xlim(wav1, wav2)\n",
"plt.xlabel('Wavelength ($\\AA$)'); plt.ylabel('Counts'); plt.grid();\n",
"mark_matched(lamp_lines)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Use the initial wavelength solution to match more lines"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# check for more matches in the range already fit\n",
"def match_more(lamp_lines, linelist, order=4, tol=2):\n",
" coeff, revcoeff = get_wavelength_solution(lamp_lines, order=order)\n",
" wfit = np.isfinite(lamp_lines['wav'])\n",
" minwav = lamp_lines['wav'][wfit].min()\n",
" maxwav = lamp_lines['wav'][wfit].max()\n",
" \n",
" xmin = lamp_lines['x'][wfit].min()\n",
" xmax = lamp_lines['x'][wfit].max()\n",
" \n",
" w = (lamp_lines['x'] > xmin) & (lamp_lines['x'] < xmax)\n",
" for line in lamp_lines[w]:\n",
" rough_wav = col_to_wav(coeff, line['x'])\n",
" refwav = match_to_list(linelist['wav'], [rough_wav], tol=tol)\n",
" if np.isfinite(refwav):\n",
" #print(f'matched column {line[\"x\"]:.1f} to wavelength {refwav[0]}')\n",
" set_line_identity(lamp_lines, linelist, line['x'], refwav)\n",
" \n",
"match_more(lamp_lines, linelist, order=???)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# re-fit with a higher order\n",
"coeff, revcoeff = get_wavelength_solution(lamp_lines, order=???)\n",
"check_wavelength_solution(lamp_spec, lamp_lines, coeff)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Use the revised wavelength solution to extend the range of matches"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Pick lines near the ends of the spectra to extend the solution across the full wavelength range."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"wav1 = 5300\n",
"wav2 = 10000\n",
"\n",
"# plot line list over lamp spectra\n",
"wavs = col_to_wav(coeff, np.arange(lamp_spec.size))\n",
"plt.plot(wavs, lamp_spec, c='g')\n",
"plt.xlim(wav1, wav2)\n",
"\n",
"# show available lines in the line list with dashed grey lines\n",
"for row in linelist:\n",
" plt.axvline(row['wav'], ls='dashed', c='0.8')\n",
"\n",
"# pick new lines to add\n",
"rough_waves = [] # <--- add new lines here\n",
"refwavs, rough_cols = match_to_list(linelist['wav'], rough_waves, plt=plt, revcoeff=revcoeff, c='r')\n",
"refcols = match_to_list(lamp_lines['x'], rough_cols, plt=plt)\n",
"\n",
"# mark existing matches\n",
"mark_peaks(plt, lamp_lines, 'wav')\n",
"plt.xlabel('Wavelength ($\\AA$)'); plt.ylabel('Counts');"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# check the supposed matches\n",
"col1, col2 = np.polyval(revcoeff, [wav1, wav2])\n",
"refcols = match_to_list(lamp_lines['x'], rough_cols)\n",
"\n",
"plt.plot(lamp_spec, c='g')\n",
"mark_peaks(plt, lamp_lines, 'x', 'y')\n",
"\n",
"plt.xlim(col1, col2)\n",
"refcols = match_to_list(lamp_lines['x'], rough_cols, plt=plt)\n",
"plt.xlabel('Column Number'); plt.ylabel('Counts');"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# add new lines to the matches\n",
"for col, wav in zip(refcols, refwavs):\n",
" set_line_identity(lamp_lines, linelist, col, wav)\n",
" \n",
"# auto-match more lines\n",
"match_more(lamp_lines, linelist)\n",
"\n",
"# re-fit\n",
"coeff, revcoeff = get_wavelength_solution(lamp_lines, order=???)\n",
"check_wavelength_solution(lamp_spec, lamp_lines, coeff)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot wavelength residuals"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Make sure you have a good wavelength solution by checking the residuals. For the supplied Keck/LRIS data, the wavelength solution should be accurate to a fraction of an Angstrom."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"w = lamp_lines['used']\n",
"std = np.std(lamp_lines['wavres'][w], ddof=1)\n",
"print(f'STD of wavelength residual is {std:0.2} Angstrom')\n",
"\n",
"plt.scatter(lamp_lines['wav'][w], lamp_lines['wavres'][w]) #, c=np.abs(lamp_lines['wavres'][w]))\n",
"plt.scatter(lamp_lines['wav'][~w], lamp_lines['wavres'][~w], marker='x') #c=np.abs(lamp_lines['wavres'][~w]))\n",
"plt.axhspan(-std, std, color='k', alpha=0.1)\n",
"plt.xlabel('Wavelength ($\\AA$)'); plt.ylabel('Counts'); plt.grid();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" \n",
" Step 2: Find the 2-D Wavelength Solution\n",
"
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the 1-D wavelength solution now solved for the middle rows using the arc-lamp spectra, we can turn our attention to the science image. A typical deep exposure will be dominated by emission lines from the night sky. These can be much brighter than features in our science target, and, as with simple photometry, we must subtract the sky contribution from each spectral band we wish to measure. \n",
"\n",
"We can use spectra spatially offset from the target to define the pure background (just like photometry). The quick way to do this is to define a background aperture (a set of rows) offset from the target aperture to define the sky spectra, determine separate wavelength solutions for the background and target apertures, and then subtract the value of the sky background at every target wavelength. However, this is often imperfect because the background and target spectra will have different spectral bins so we must interpolate the background spectrum to the wavelengths of the target spectrum. There are trivial ways to do this, but these can perform poorly near bright sky lines.\n",
"\n",
"Instead, we can use all the information available on the 2-D image to model the sky lines (see [Kelson 2003](https://ui.adsabs.harvard.edu/abs/2003PASP..115..688K/abstract)). Because lines of constant wavelength bend with row, each sky line is sampled in hundreds of different spectral bands. Using all of these (and excluding rows contaminated by target light, etc.) we can model the sky background as a function of wavelength on a fine scale and use this wealth of information to subtract the skylines from our target spectra. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Find the rows dominated by sky light (exclude rows with object light, bad rows)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# display the science image\n",
"show_image(image)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# compute the row averages and normalize so that the background is near 0 and the peaks are near 1\n",
"rowaverage = image.mean(axis=1)\n",
"rowaverage -= np.median(rowaverage)\n",
"rowaverage /= rowaverage.max()\n",
"plt.plot(ys, rowaverage)\n",
"plt.xlabel('Row Number (y-coordinate)'), plt.ylabel('Normalized Row Average')\n",
"plt.grid();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Notes\n",
"* the plot above should show peaks in the rows containing object light\n",
"* you may also notice a few bad rows (actually bad columns!) with counts significantly below the median."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"Then create a 2D boolean mask that is `True` for pixels that are dominated by sky light and `False` for other pixels."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# find the rows with object light\n",
"objrows = ys[???]\n",
"\n",
"# add some margin to object rows\n",
"ngrow = ??? # number of rows to add above and below object rows\n",
"newobjrows = []\n",
"for row in objrows:\n",
" newobjrows.extend([row + i for i in np.arange(-ngrow, ngrow + 1)])\n",
"objrows = np.unique(newobjrows)\n",
"\n",
"# mask to mark sky rows\n",
"skymask = np.ones(image.shape, dtype=bool)\n",
"skymask[objrows, :] = False\n",
"\n",
"# also exclude bad rows\n",
"badrows = ys[rowaverage < ???]\n",
"skymask[badrows, :] = False\n",
"\n",
"# rows with mostly sky background light\n",
"skyrows = ys[skymask.mean(axis=1) == 1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# median (unmasked) sky spectrum and standard deviation\n",
"medspec = np.median(???)\n",
"stdspec = np.std(???)\n",
"\n",
"# exclude deviant pixels from the skymask\n",
"pull = (image - medspec) / stdspec\n",
"w = pull > 5\n",
"skymask[w] = False"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# show the mask\n",
"plt.figure(figsize=(15, 3))\n",
"plt.imshow(skymask, origin='lower', aspect='auto');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Look at a small section of image up close"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# cut out a small image \"stamp\" near the center of the frame\n",
"row = cy\n",
"col = cx\n",
"hwidth = 50\n",
"stamp = image[:, col - hwidth : col + hwidth]\n",
"ys_stamp = yvals[:, col - hwidth : col + hwidth]\n",
"xs_stamp = xvals[:, col - hwidth : col + hwidth]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# show the image stamp\n",
"extent = (xs_stamp.min(), xs_stamp.max(), ys_stamp.min(), ys_stamp.max())\n",
"show_image(stamp, extent=extent);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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. \n",
"\n",
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Note\n",
"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:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot stamp values against column numbers\n",
"plt.plot(xs_stamp.ravel(), stamp.ravel(), 'r.');\n",
"plt.xlabel('Column Number'), plt.ylabel('Counts');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Map out lines of constant wavelength (determine the wavelength for each pixel)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"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`)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def get_dl_model(params, dxs, dys):\n",
" return dxs + params[0] * dys + params[1] * dys ** 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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`)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# pixel offsets from the refernece pixel\n",
"dxs = xs_stamp - col\n",
"dys = ys_stamp - row\n",
"\n",
"# parameter guess\n",
"guess = ???\n",
"\n",
"# get the wavelength offsets and plot vs. counts\n",
"dls = get_dl_model(guess, dxs, dys)\n",
"plt.plot(dls.ravel(), stamp.ravel(), 'r.')\n",
"plt.xlabel('Wavelength Offset')\n",
"plt.ylabel('Counts');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit a spline to the counts spectrum above"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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).\n",
"\n",
"We can define a function that will return the best fit spline object (in the least-squares sense) for our data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def get_profile_spl(dls, stamp):\n",
" # need to sort the data (and weights) so that the x values are increasing\n",
" x, y = dls.ravel(), stamp.ravel()\n",
" weights = np.sqrt(np.abs(y)) # not technically optimal for coadded data, but ok\n",
" wsort = x.argsort()\n",
" x, y, weights = x[wsort], y[wsort], weights[wsort]\n",
"\n",
" # set locations for spline knots\n",
" t = np.linspace(x.min() + 1, x.max() - 1, np.int(x.max() - x.min()))\n",
" spl = LSQUnivariateSpline(x, y, t, weights)\n",
" return x, y, spl"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Lets see how good this looks."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fit a spline to the data and plot\n",
"x, y, spl = get_profile_spl(dls, stamp)\n",
"\n",
"fig, axarr = plt.subplots(2, sharex=True)\n",
"axarr[0].plot(x, y, 'r.')\n",
"axarr[0].plot(x, spl(x))\n",
"axarr[1].plot(x, y - spl(x), 'r.')\n",
"plt.ylim(-200, 200)\n",
"plt.xlabel('Wavelength Offset');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Notice\n",
" * 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\n",
" * the parameter guess above may be imperfect resulting in systematic residuals (we can minimize these below)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Minimize residuals from the spline model to determine the best tilt/curvature model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def check_dl_model(params, dxs, dys, stamp):\n",
" dls = get_dl_model(params, dxs, dys)\n",
" x, y, spl = get_profile_spl(dls, stamp)\n",
" chisq = np.sum((stamp - spl(dls)) ** 2 / np.sqrt(np.abs(stamp)))\n",
" return chisq / (stamp.size - len(params))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# see how good our guess is\n",
"check_dl_model(guess, dxs, dys, stamp)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we just need to change the model parameters to minimize the residuals. We can use `fmin` for this."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# get the best model parameters for this stamp\n",
"params = fmin(check_dl_model, guess, args=(dxs, dys, stamp))\n",
"print(\"best model parameters are\", params)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can plug in these parameters for your `guess` above and re-run the cells to check the residuals if you want."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Find the lines of constant wavelength at other parts of the 2-D image"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"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)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# define the column centers for the stamps\n",
"hwidth = 200\n",
"cols = np.arange(hwidth, nx, 2 * hwidth)\n",
"cols"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# reference wavelength offsets to the center row\n",
"row = cy\n",
"\n",
"# define a 2D array to hold the wavelength offsets for each pixel\n",
"lambdas = np.zeros(image.shape) \n",
"\n",
"# loop over each central column\n",
"for col in cols:\n",
" print('col = ', col)\n",
" \n",
" # slice the data\n",
" inds = np.s_[:, col - hwidth : col + hwidth]\n",
" stamp = image[inds]\n",
" mask = skymask[inds]\n",
" dys = yvals[inds] - row\n",
" dxs = xvals[inds] - col\n",
"\n",
" # initial fit\n",
" params = fmin(check_dl_model, guess, args=(dxs[mask], dys[mask], stamp[mask]))\n",
" \n",
" # check for outliers\n",
" dls = get_dl_model(guess, dxs, dys)\n",
" x, y, spl = get_profile_spl(dls, stamp)\n",
" model = spl(dls)\n",
" pull = (stamp - model) / np.sqrt(np.abs(model))\n",
" w = (pull < 5) & mask\n",
" params2 = fmin(check_dl_model, params, args=(dxs[w], dys[w], stamp[w]))\n",
"\n",
" # record\n",
" lambdas[inds] = get_dl_model(params2, dxs, dys) + col\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Look at a few rows and see how the wavelength offsets vary with column number. We can fit a low-order polynomial to these."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# just plot offsets for a few of the rows across the image\n",
"order = ???\n",
"for y in range(10, ny, 50):\n",
" p = plt.plot(cols, lambdas[y, cols] - xs[cols], 'o')\n",
" c = np.polyfit(cols, lambdas[y, cols] - xs[cols], order)\n",
" plt.plot(xs, np.polyval(c, xs), c=p[0].get_color(), label='row {}'.format(y))\n",
"plt.legend()\n",
"plt.xlabel('Column Number')\n",
"plt.ylabel('Wavelength Offset from Middle Row');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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).\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# get the lambda values for the entire image (fit)\n",
"lambdafit = np.zeros(image.shape)\n",
"for y in range(ny):\n",
" c = np.polyfit(cols, lambdas[y, cols] - xs[cols], order)\n",
" lambdafit[y, :] = np.polyval(c, xs) + xs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" \n",
" Step 3: Flat Field the 2D Spectra (Optional)\n",
"
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To this point we have not corrected the data for pixel-to-pixel efficiencies (in other words, we have not flat-fielded the data). For many tasks this is actually not required, but if you have *good* flat field exposures available, flat fielding can improve the output data quality.\n",
"\n",
"The key for spectroscopy is to note that the CCD response is wavelength dependent (you wouldn't use a U-band flat to process an I-band image, right?), and with spectra there are many wavelength bins to account for. Just as with photometry, we do not care so much how many counts we get in a given wavelength bin, we can just determine the average number in each wavelength bin and then divide by this average to get numbers close to 1."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# load in the flat\n",
"flat = fits.getdata(os.path.join(data_dir, 'spec_flat.fits'))\n",
"show_image(flat)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that the average counts vary considerably with wavelength and there are even a few line features present in the data (there is air between the light source and the camera, and air can emit and absorb light!).\n",
"\n",
"Since we know the wavelengths for each pixel, this is no problem. We can fit a spline to this profile (which is marginalized over slit position) to determine the typical number of counts in each wavelength bin and then divide each pixel by the appropriate value from the spline fit to normalize the counts."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fit a spline\n",
"x = lambdafit.ravel()\n",
"y = flat.ravel()\n",
"weights = np.sqrt(np.abs(y)) \n",
"wsort = x.argsort()\n",
"x, y, weights = x[wsort], y[wsort], weights[wsort]\n",
"t = np.linspace(x.min() + 1, x.max() - 1, nx)\n",
"spl = LSQUnivariateSpline(x, y, t, weights)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot the spline\n",
"plt.plot(x, spl(x))\n",
"plt.xlabel('Column Number')\n",
"plt.ylabel('Counts');"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# normalize the flat field \n",
"normflat = flat / ????"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show_image(flat)\n",
"show_image(normflat)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can now divide the science images by this normalized flat to correct for pixel-to-pixel differences."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# flat field the 2D science image using the normalized flat\n",
"image /= normflat"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" \n",
" Step 4: Model the 2D Sky\n",
"
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# function to fit a 2D spline\n",
"def fit_sky(xvals, yvals, image, ky=1, dx=0.5):\n",
" # select knot points in the x (wavelength) direction\n",
" tx = np.arange(xvals.min() + 2, xvals.max() - 2, dx)\n",
" \n",
" # select knot points in the y (spatial) direction\n",
" ty = [] # if there are no knots, the fit will be a poly nomial of degree ky\n",
" \n",
" # fit the 2D spline\n",
" return LSQBivariateSpline(xvals.ravel(), yvals.ravel(), image.ravel(), tx, ty, ky=ky)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# use the (unmasked) sky background pixels and fit the 2D spline\n",
"skyfit = fit_sky(lambdafit[skymask], yvals[skymask], image[skymask])\n",
"\n",
"# evaluate the 2D sky at every pixel\n",
"sky = skyfit.ev(lambdafit, yvals)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show_image(image)\n",
"show_image(sky)\n",
"show_image(image-sky)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" \n",
" Step 5: Extract the 1-D Target Spectrum\n",
"
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can sum the counts in the sky subtracted image above to create a counts spectrum of our target. As with photometry we need to define an aperture (or a construct a probability map) to mark which pixels contain target light so that we can reject neighboring objects. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Locate object trace"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are actually two object traces in the supplied Keck/LRIS 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. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def get_profile_model(params, ys):\n",
" a1, cy1, sigma1, a2, cy2, sigma2 = params\n",
" \n",
" p1 = np.exp(-(ys - cy1)**2 / 2 / sigma1**2) \n",
" p1 /= p1.max()\n",
" \n",
" p2 = np.exp(-(ys - cy2)**2 / 2 / sigma2**2) \n",
" p2 /= p2.max()\n",
"\n",
" return a1 * p1 + a2 * p2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# get the median for each row\n",
"profile = np.median(image - sky, axis=1)\n",
"\n",
"# starting guess for the profile model\n",
"guess = ???\n",
"model = get_profile_model(guess, ys)\n",
"\n",
"plt.plot(ys, profile)\n",
"plt.plot(ys, model);\n",
"plt.xlabel('Row Number')\n",
"plt.ylabel('Median Row Counts');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course, we can improve the model by minimizing the model residuals (or $\\chi^2$/dof)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def get_profile_chisq(params, ys, profile):\n",
" model = get_profile_model(params, ys)\n",
" return np.sum( (profile - model)**2 / np.sqrt(np.abs(profile)) ) / (profile.size - len(params))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fit for the best model\n",
"params = fmin(get_profile_chisq, guess, args=(ys, profile))\n",
"print(\"best fit parameters are\", params)\n",
"\n",
"model = get_profile_model(params, ys)\n",
"plt.plot(ys, profile)\n",
"plt.plot(ys, model)\n",
"plt.xlabel('Row Number')\n",
"plt.ylabel('Median Row Counts')\n",
"plt.grid();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# fit the profile centered at these columns\n",
"hwidth = 50\n",
"cols = np.arange(hwidth, nx + 1, 2 * hwidth)\n",
"\n",
"ycenter = np.zeros( (len(cols), 2) )\n",
"for icol, col in enumerate(cols):\n",
" stamp = (image - sky)[:, col - hwidth : col + hwidth]\n",
" profile = np.mean(stamp, axis=1)\n",
" params = fmin(get_profile_chisq, guess, args=(ys, profile))\n",
" ycenter[icol, :] = params[[1, 4]]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fit the relation with a polynomial\n",
"ind = 0 # which trace 0 or 1?\n",
"t_order = ???\n",
"trace_c = np.polyfit(cols, ycenter[:, ind], t_order)\n",
"fig, axarr = plt.subplots(2, sharex=True)\n",
"axarr[0].plot(cols, ycenter[:, ind], 'ro')\n",
"axarr[0].plot(xs, np.polyval(trace_c, xs), 'r')\n",
"axarr[0].axes.set_ylabel('y-coordinate'); axarr[0].grid();\n",
"axarr[1].plot(cols, ycenter[:, ind] - np.polyval(trace_c, cols), 'ro')\n",
"axarr[1].axes.set_ylim(-0.5, 0.5)\n",
"axarr[1].axes.set_ylabel('Fit Residual (pixels)')\n",
"plt.xlabel('Column Number'); axarr[1].grid();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The typical residuals should be a small fraction of a pixel. \n",
"\n",
"**Note** we are fitting a Gaussian model to a non-Gaussian (galaxy) profile, so the center value may be systematically offset.\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# position offsets from the object trace (defined to be at slitpos = 0)\n",
"slitpos = yvals - np.polyval(trace_c, yvals)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Spatial profile of the target"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, determine the spatial profile of the target and stretch this over the whole 2D image to make a weight map for later use.\n",
"\n",
"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)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# subtract the sky\n",
"nosky = image - sky\n",
"\n",
"# normalize to the pixel brightness at the trace center\n",
"yinds = (np.round(np.polyval(trace_c, xs))).astype(int)\n",
"normed = nosky / ???"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Fit a spline to the normalized profile"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# get 1D arrays with the positions along the slit and the normalized counts\n",
"pos = slitpos.flatten()\n",
"counts = normed.flatten()\n",
"\n",
"# sort by slit position\n",
"sort_inds = pos.argsort()\n",
"pos, counts = pos[sort_inds], counts[sort_inds]\n",
"\n",
"# fit a spline to model the spatial profile\n",
"t = np.linspace(pos.min() + 2, pos.max() - 2, ny // 2) # spline knot points\n",
"profile_spl = LSQUnivariateSpline(pos, counts, t)\n",
"\n",
"# remove outliers and re-fit\n",
"diff = counts - profile_spl(pos)\n",
"sample = sigma_clip(diff)\n",
"w = ((np.abs(diff) / sample.std()) < 5) & np.isfinite(diff)\n",
"profile_spl = LSQUnivariateSpline(pos[w], counts[w], t)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot the target profile\n",
"plt.plot(pos, profile_spl(pos) )\n",
"plt.xlim(-40, 50)\n",
"plt.grid();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculate the profile for every column to make a 2D profile image."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# create the profile image\n",
"profile_image = profile_spl(???)\n",
"\n",
"# de-weight negative values in provile_image\n",
"profile_image[profile_image < 0] = 0"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show_image(profile_image, upper=50)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Extract the target spectrum (basic version)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# select which rows to sum\n",
"w = (slitpos > ???) & (slitpos < ???)\n",
"ymin, ymax = yvals[w].min(), yvals[w].max()\n",
"\n",
"# calculate the sum\n",
"spec_basic = ???\n",
"\n",
"# sky background\n",
"skybg_basic = ???\n",
"\n",
"# plot the extracted spectrum\n",
"plt.plot(xs, spec_basic)\n",
"plt.xlabel('Column Number')\n",
"plt.ylabel('Counts');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Optimal extraction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# calculate the weighted average (for each column)\n",
"spec_opt = (nosky * profile_image)[ymin:ymax, :].sum(axis=0) / profile_image.sum(axis=0)\n",
"\n",
"# calculate the bias factor needed to scale the average to a sum\n",
"bias_factor = np.median(spec_basic / spec_opt)\n",
"spec_opt *= bias_factor\n",
"\n",
"# same for the sky background\n",
"skybg_opt = ???\n",
"bias_factor_sky = ???\n",
"skybg_opt *= bias_factor_sky\n",
"\n",
"# plot the extracted spectrum\n",
"plt.plot(xs, spec_basic, label='basic extraction')\n",
"plt.plot(xs, spec_opt, label='optimal extraction')\n",
"plt.xlabel('Column Number')\n",
"plt.ylabel('Counts');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you look closely you should see that the basic extraction is noisier than the optimal."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" \n",
" Step 6: Tweak the Wavelength Solution Using Sky Lines\n",
"
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Spectrographs, including Keck/LRIS, may warp slightly as the telescope they are mounted to moves across the sky, and this can result in wavelength drift (the wavelength of a given pixel may change over time). Thus the wavelength solution we derived using the arc-lamp spectra may not apply exactly to our science data. One remedy for this is to use the wavelengths of sky lines to correct for spectral drift."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Start with the wavelength solution from the arc-lamps using the target aperture"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# re-extract lamp spectra using the profile_image\n",
"lamp_basic = lamp_image[ymin:ymax, :].mean(axis=0)\n",
"lamp_spec_opt = (lamp_image * profile_image)[ymin:ymax, :].sum(axis=0) / profile_image.sum(axis=0)\n",
"bias_factor = np.median(lamp_basic / lamp_spec_opt)\n",
"lamp_spec_opt *= bias_factor\n",
"\n",
"# re-compute lamp_lines\n",
"lamp_lines_opt = get_lamp_lines(lamp_spec_opt)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# use the old wavelength solution to identify lines\n",
"rough_waves = col_to_wav(coeff, lamp_lines_opt['x'])\n",
"lamp_lines_opt['wav'] = match_to_list(linelist['wav'], rough_waves)\n",
"\n",
"# re-fit\n",
"coeff2, revcoeff2 = get_wavelength_solution(lamp_lines_opt, order=4)\n",
"check_wavelength_solution(lamp_spec_opt, lamp_lines_opt, coeff2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"w = lamp_lines_opt['used']\n",
"plt.scatter(lamp_lines_opt['wav'][w], lamp_lines_opt['wavres'][w]) #, c=np.abs(lamp_lines['wavres'][w]))\n",
"plt.scatter(lamp_lines_opt['wav'][~w], lamp_lines_opt['wavres'][~w], marker='x') #c=np.abs(lamp_lines['wavres'][~w]))\n",
"plt.grid()\n",
"\n",
"std = np.std(lamp_lines_opt['wavres'][w], ddof=1)\n",
"plt.axhspan(-std, std, color='k', alpha=0.1)\n",
"print(f'STD of wavelength residual is {std:0.2} Angstrom')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# you may need to exclude strong outliers and re-fit\n",
"coeff2, revcoeff2 = get_wavelength_solution(lamp_lines_opt[w], order=4)\n",
"check_wavelength_solution(lamp_spec_opt, lamp_lines_opt[w], coeff2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Use night sky spectrum to tweak wavelength solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We need the precise wavelengths of all sky lines. These can be taken from the [UVES sky emission spectrum](https://www.eso.org/observing/dfo/quality/UVES/pipeline/sky_spectrum.html), which was constructed through careful analysis of high-resolution spectroscopy. Copies of the UVES nigh sky line lists are provided in the `UVES/` subdirectory of `data/`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# load skyline list\n",
"dtype = []\n",
"dtype.append( ('seq', int) )\n",
"dtype.append( ('wav', float) )\n",
"dtype.append( ('intensity', float) )\n",
"dtype.append( ('fwhm', float) )\n",
"dtype.append( ('flux', float) )\n",
"\n",
"skylines = np.array([], dtype=dtype)\n",
"wildcard = os.path.join(data_dir, 'UVES', '*dat')\n",
"fnames = glob(wildcard)\n",
"for fname in fnames:\n",
" lines = np.genfromtxt(fname, dtype=dtype, skip_header=3, skip_footer=1)\n",
" skylines = np.hstack([skylines, lines])\n",
" \n",
"SKYLINES = skylines"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can convert the line list into a 1-D spectrum by modeling each line as a Gaussian, scaling the Gaussian to the measured line flux, and then adding each line to a continuum. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# function to convert night sky linelist into spectral template\n",
"def model_sky(wav, *params):\n",
" wavshift, scale, fwhm, conta, contb = params\n",
" \n",
" # start with the continuum\n",
" model = np.polyval([conta, contb], wav)\n",
" \n",
" # add in each line\n",
" sigma = fwhm / 2.355\n",
" minwav = np.min(wav) - 3 * fwhm\n",
" maxwav = np.max(wav) + 3 * fwhm\n",
" w = ((SKYLINES['wav'] + wavshift) >= minwav) & ((SKYLINES['wav'] + wavshift) <= maxwav)\n",
" for line in skylines[w]:\n",
" # assume lines are Gaussian\n",
" model += scale * line['flux'] * np.exp( -(wav - line['wav'] - wavshift) ** 2 / 2 / sigma ** 2)\n",
" \n",
" return model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can compare this model to the data to derive the wavelength drift (and also the instrumental resolution). Start by setting the approximate `scale` and `fwhm` for the data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"# get starting guess for the fwhm and flux scaling\n",
"wavs = col_to_wav(coeff2, np.arange(skybg_opt.size))\n",
"params = ???\n",
"skymodel = model_sky(wavs, *params)\n",
"w = skymodel > 1\n",
"scale = np.mean(skybg_opt[w] / skymodel[w])\n",
"\n",
"plt.plot(wavs, skybg_opt, label='Observed')\n",
"plt.plot(wavs, scale * skymodel, label='Model', alpha=0.5)\n",
"plt.legend(); plt.grid()\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('Counts');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `fwhm`, `scale`, etc., will vary with wavelength (note the `fwhm` is related to the instrumental resolution, which is **not** constant with wavelength for most spectrographs!). The wavelength drift may also vary along the spectra, so we should fit for the model parameters in a number of wavelength ranges. Some parts of the spectra may lack significant sky lines, so we should divide up the spectra in a way that ensures each range contains at least a few lines."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# locate observed sky line peaks (roughly)\n",
"skyline_cols, properties = find_peaks(skybg_opt, prominence=???)\n",
"skyline_wavs = col_to_wav(coeff2, skyline_cols)\n",
"\n",
"plt.plot(wavs, skybg_opt, label='Observed')\n",
"plt.scatter(skyline_wavs, skybg_opt[skyline_cols]);\n",
"\n",
"# divide into sections with at least a few lines\n",
"nmin = ???\n",
"minwav = wavs.min()\n",
"wavranges = []\n",
"while (skyline_wavs > minwav).sum() > nmin:\n",
" w = skyline_wavs > minwav\n",
" maxwav = skyline_wavs[w][nmin - 1]\n",
" dwav = (maxwav - minwav) / 2\n",
" wavranges.append([minwav - dwav, maxwav + dwav]) \n",
" plt.axvline(maxwav, color='k', ls='dotted')\n",
" minwav = maxwav\n",
"wavranges.append([maxwav, wavs.max()])\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('Counts')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In each of these `wavranges` find the best set of model parameters to match the observed sky spectrum. As a bonus, you will find the instrumental resolution for each wavelength band as will. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fit for the wavelength shift section by section\n",
"dtype = []\n",
"dtype.append( ('cwav', float) )\n",
"dtype.append( ('wavshift', float) )\n",
"dtype.append( ('wavshift_err', float) )\n",
"dtype.append( ('scale', float) )\n",
"dtype.append( ('scale_err', float) )\n",
"dtype.append( ('fwhm', float) )\n",
"dtype.append( ('fwhm_err', float) )\n",
"\n",
"guess = 0, scale, 5, 0, 0\n",
"results = []\n",
"for wavrange in wavranges:\n",
" minwav, maxwav = wavrange\n",
" w = (wavs >= minwav) & (wavs < maxwav)\n",
" bounds = ((-50, 0, 0, -np.inf, -np.inf), (50, np.inf, np.inf, np.inf, np.inf))\n",
" popt, pcov = curve_fit(model_sky, wavs[w], skybg_opt[w], p0=guess, bounds=bounds)\n",
" cols = [(minwav + maxwav) / 2]\n",
" for i in range(3):\n",
" cols.append(popt[i])\n",
" cols.append(np.sqrt(np.diag(pcov)[i]))\n",
" \n",
" newrow = np.array(tuple(cols), dtype=dtype)\n",
" results.append(newrow)\n",
"results = np.array(results)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# check the last wavelength range fit\n",
"model = model_sky(wavs, *popt)\n",
"plt.plot(wavs[w], skybg_opt[w], label='Observed')\n",
"plt.plot(wavs[w], model[w], label='Model', alpha=0.5)\n",
"plt.legend(); plt.grid()\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('Counts');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can fit a polynomial to the measured wavelength shifts so that we can interpolate the results to any wavelength. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fit for the wavelength offset\n",
"cwavshift = np.polyfit(results['cwav'], results['wavshift'], ???)\n",
"plt.plot(results['cwav'], results['wavshift'])\n",
"plt.plot(wavs, np.polyval(cwavshift, wavs))\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('Wavelength Shift ($\\AA$)')\n",
"plt.grid()\n",
"\n",
"# correct wavelengths to match skyline values\n",
"wavs_corr = wavs - np.polyval(cwavshift, wavs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As a bonus, we can also fit for the instrumental resolving power as a function of wavelength, $R(\\lambda) = \\lambda / \\delta\\lambda$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fit for the resolving power (lambda / fwhm)\n",
"cres = np.polyfit(results['cwav'], results['cwav'] / results['fwhm'], ???)\n",
"plt.plot(results['cwav'], results['cwav'] / results['fwhm'])\n",
"plt.plot(wavs, np.polyval(cres, wavs))\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('Spectral Resolution')\n",
"plt.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have the drift-corrected mapping from counts to wavelengths, we can apply this to the full spectrum."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# apply the wavelength solution to the extracted spectrum\n",
"plt.plot(wavs_corr, spec_opt);\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('Counts');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Step 7: Flux Calibration (Optional)
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Sometimes you may also want to scale the counts (by a wavelength dependent function) to determine the flux densities in each wavelength bin. To do this, you need to observe a target with a known flux density spectrum. A number of such calibrated standard stars are available and should be observed. Typically you cannot observe your target and a standard star simultaneously, so you need a photometric night in order to get the absolute scaling right. However, even with clouds you can often get the *relative* flux densities of your target."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# extract standard star spectrum as above\n",
"stdimage = fits.getdata(os.path.join(data_dir, 'spec_std.fits'))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# compute the row averages and normalize so that the background is near 0 and the peaks are near 1\n",
"rowaverage = stdimage.mean(axis=1)\n",
"rowaverage -= np.median(rowaverage)\n",
"rowaverage /= rowaverage.max()\n",
"\n",
"# find the rows with object light\n",
"objrows = ys[rowaverage > ???]\n",
"\n",
"# add some margin to object rows\n",
"ngrow = ??? # number of rows to add above and below object rows\n",
"newobjrows = []\n",
"for row in objrows:\n",
" newobjrows.extend([row + i for i in np.arange(-ngrow, ngrow + 1)])\n",
"objrows = np.unique(newobjrows)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# mask to mark sky rows\n",
"stdskymask = np.ones(stdimage.shape, dtype=bool)\n",
"stdskymask[objrows, :] = False\n",
"\n",
"# also exclude bad rows\n",
"stdskymask[badrows, :] = False\n",
"\n",
"# median (unmasked) sky spectrum and standard deviation\n",
"medspec = np.median(stdimage[skyrows, :], axis=0)\n",
"stdspec = np.std(stdimage[skyrows, :], axis=0, ddof=1)\n",
"\n",
"# exclude deviant pixels from the skymask\n",
"pull = (stdimage - medspec) / stdspec\n",
"w = pull > 3\n",
"stdskymask[w] = False"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"stdskyfit = fit_sky(lambdafit[stdskymask], yvals[stdskymask], stdimage[stdskymask])\n",
"stdsky = stdskyfit.ev(lambdafit, yvals)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# compare the 2D spectra, background model, and sky subtracted spectra\n",
"show_image(stdimage)\n",
"show_image(stdsky)\n",
"show_image(stdimage - stdsky)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# select which rows to sum\n",
"w = (slitpos > ???) & (slitpos < ???)\n",
"ymin, ymax = yvals[w].min(), yvals[w].max()\n",
"\n",
"# calculate the sum\n",
"std_spec = (stdimage - stdsky)[ymin:ymax, :].sum(axis=0)\n",
"\n",
"# plot the extracted spectrum\n",
"plt.plot(wavs_corr, std_spec)\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('Counts');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The supplied standard star is called Feige 67. You can get a table of its flux densities from [ESO](https://www.eso.org/sci/observing/tools/standards/spectra.html) (and other sources)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Download the data from ESO\n",
"dtype = []\n",
"dtype.append( ('wav', float) )\n",
"dtype.append( ('flux', float) ) # units are ergs/cm/cm/s/A * 10**16\n",
"dtype.append( ('eflux', float) )\n",
"dtype.append( ('dlam', float) )\n",
"calspec = np.genfromtxt('ftp://ftp.eso.org/pub/stecf/standards/okestan/ffeige67.dat', dtype=dtype)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot the tabulated flux densities\n",
"plt.plot(calspec['wav'], calspec['flux']);\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel(r'$10^{16} \\times F_{\\lambda}$');\n",
"plt.yscale('log');"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fit a spline to the tabulated spectrum\n",
"t = np.arange(calspec['wav'][1], calspec['wav'][-2], np.int(np.median(calspec['dlam'])))\n",
"stdflux = LSQUnivariateSpline(calspec['wav'], calspec['flux'], t, calspec['eflux'])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# get the counts to flux density ratios in each wavelength bin\n",
"# exclude significant line features (fluxes near these depend on spectral resolution)\n",
"ratios = std_spec / stdflux(wavs)\n",
"w = (wavs_corr > calspec['wav'].min()) \\\n",
" & (wavs_corr < calspec['wav'].max()) \\\n",
" & (np.abs(wavs_corr - 7650) > 70) \\\n",
" & (np.abs(wavs - ???) > ???) \\\n",
" & (np.abs(wavs - ???) > ???) \\\n",
"\n",
"# fit a spline to the ratios to determine the response function\n",
"t = wavs_corr[w][1:-2:???]\n",
"respfn = LSQUnivariateSpline(wavs_corr[w], ratios[w], t)\n",
"\n",
"plt.plot(wavs_corr[w], ratios[w], 'ro')\n",
"xwav = np.linspace(wavs_corr[w][1], wavs_corr[w][-1], 1000)\n",
"plt.plot(xwav, respfn(xwav));\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('Response Function (counts/$F_{\\lambda}$)');"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# compare the tabulated and extracted flux densities (applying the response function)\n",
"plt.plot(calspec['wav'], calspec['flux'], label='Tabulated (published) Spectrum');\n",
"plt.plot(wavs_corr, std_spec / respfn(wavs), label='Extracted Spectrum')\n",
"plt.xlim(5300, 9200) # Feige 67 values are only tabulated out to ~9200 A\n",
"plt.ylim(0, 1000)\n",
"\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('$F_{\\lambda}$')\n",
"plt.ylim(50, 1000)\n",
"plt.yscale('log')\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"Now apply the response function to the target spectrum to get the flux calibrated spectrum."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.plot(wavs_corr, spec_opt / respfn(wavs))\n",
"plt.xlim(5300, 9200)\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('$F_{\\lambda}$');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" \n",
" Step 8: Removal of Telluric Features (Optional)\n",
"
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Earth's atmosphere strongly absorbs light in bands near 7600 angstroms (A-band) and 6900 angstroms (B-band), and also in a number of other bands particularly in the near-IR. These bands can be seen in the spectra of our standard star. The target spectrum is also absorbed in these bands. We can correct for this by measuring the fraction of the standard star's light absorbed as a function of wavelength compared to an extrapolation of the standard star's continuum across these bands. We can then scale this absorption strength to match the fractional absorption in the target spectrum and thus divide out the telluric absorption. This works well as long as the atmospheric conditions remain similar."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# compute the telluric spectrum (fraction of atmospheric absorption per wavelength bin)\n",
"tell_spec = (std_spec / respfn(wavs)) / stdflux(wavs)\n",
"w = (np.abs(wavs_corr - 7650) < 70) | (np.abs(wavs_corr - 6900) < 40)\n",
"tell_spec[~w] = 1\n",
"\n",
"plt.plot(wavs_corr, tell_spec)\n",
"plt.xlim(5300, 9200)\n",
"plt.ylim(0, 1.1)\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('$F_{\\lambda}$');"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# correct the target spectrum for telluric absorption\n",
"plt.plot(wavs_corr, spec_opt / respfn(wavs), label='uncorrected')\n",
"plt.plot(wavs_corr, spec_opt / respfn(wavs) / tell_spec, label='corrected')\n",
"plt.xlim(5300, 9200)\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('$F_{\\lambda}$')\n",
"plt.yscale('log')\n",
"plt.ylim(1, 2000)\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Comparison of final results to published spectra"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The provided Keck/LRIS spectra have previously been [published](https://ui.adsabs.harvard.edu/#abs/2016ApJ...830...13P/abstract). You can compare your extraction to the \n",
"published spectra, which are available from the [Weizmann Interactive Supernova Data Repository (WISeREP)](https://wiserep.weizmann.ac.il)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# download the published spectrum\n",
"dtype = []\n",
"dtype.append( ('wav', float) )\n",
"dtype.append( ('flux', float) )\n",
"dtype.append( ('sky', float) )\n",
"dtype.append( ('eflux', float) )\n",
"dtype.append( ('xpix', float) )\n",
"dtype.append( ('ypix', float) )\n",
"dtype.append( ('resp', float) )\n",
"url = 'https://wiserep.weizmann.ac.il/system/files/uploaded/general/PTF12dam_2014-04-30_10-33-34_Keck1_LRIS_None.spec'\n",
"pub = np.genfromtxt(url, dtype=dtype)\n",
"\n",
"plt.plot(wavs_corr, spec_opt / respfn(wavs) / tell_spec, label='my extraction')\n",
"plt.plot(pub['wav'], pub['flux'] * 1e16, label='published')\n",
"plt.xlim(5300, 9200)\n",
"plt.xlabel('Wavelength ($\\AA$)')\n",
"plt.ylabel('$F_{\\lambda}$')\n",
"plt.yscale('log')\n",
"plt.ylim(0.01, 2000)\n",
"plt.legend();"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.10"
}
},
"nbformat": 4,
"nbformat_minor": 2
}