Python Example of Longslit Spectroscopy

3. Python Example of Longslit Spectroscopy#

In this note, I show how a longslit spectrum data can be reduced.

Original data: 2-D FITS files (bias, dark, flat, lamp, science image)
Input data for this note: Output (preprocessed data) of 45P Spectroscopic Preprocessing
Final data from this note: 1-D spectrum (flux at each λ)

Note

Let me repeat again: This lecture note is ignoring anything about preprocessing. If preprocessing-related lecture is required, I recommend Matthew Craig’s CCD Data Reduction Guide.

See 45P Spectroscopic Preprocessing for preprocessing the 45P data that will be used in this lecture note.

3.1. About the Data#

The data we will use is the longslit spectrum data of a periodic comet, 45P/Honda–Mrkos–Pajdušáková.

Details of the Observational

Nishi-Harima Astronomical Observatory (NHAO).
Nayuta Telescope (D=2.0m, F/12, Nasmyth)
In 2017 February
FLI PL230-BI CCD (pixel 15µm), 0.34 arcsec/pixel.
Used MALLS (official website in Japanese).
- Grating: 150 lines/mm, so-called “Low1” (or just “Low”)
- Order cut filter: WG320
- λ = 3700–9500Å, total wavelength range = 5700Å
- Slit width = 1.6 arcsec
- resolution R~600 (@550nm, under slit width 1.2”)
- Basic info is written in file headers.

3.2. Postprocessing#

After the preprocessing, there’s still a long way to go. There is no single best way, so we need a lot of interaction with computers. That includes

Wavelength identification (ID): First, you cut a small fraction of the comparison lamp image and let the computer calculate the wavelength \(\lambda = \lambda(x)\) for given pixel value x. This is a “crude version” of the full 2-D map of wavelength (\(\lambda(x, y)\)).
Reidentification (REID): Then you must check that this function makes sense throughout the whole image. REID checks it. Finally, it will generate the full 2-D map of wavelength (\(\lambda(x, y)\)).
Aperture trace (aptrace): The object spectrum cannot be perfectly parallel to the dispersion direction. Aptrace (in IRAF, APALL) is to find how the source moves along the dispersion direction. If, e.g., x-axis is the dispersion direction, aptrace is the central location of the object’s peak, \(y_\mathrm{obj}(x)\).
Sky subtraction: The sky should be subtracted at each wavelength. This is based on the pixel values at the same λ at different spatial location.
Aperture sum: We then sum the pixel values within the selected aperture along the spatial direction.

3.3. Preparation#

When you use the code for your data, you must tune each parameter for your own instrument.

# %matplotlib notebook
%config InlineBackend.figure_format = 'retina'
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = 'last_expr'
from matplotlib import pyplot as plt

from pathlib import Path
import numpy as np

from numpy.polynomial.chebyshev import chebfit, chebval

# %matplotlib notebook
# comment above out when you copy this code to your spyder, etc.

# We need to do it in a separate cell. See:
# https://github.com/jupyter/notebook/issues/3385
plt.style.use('default')
plt.rcParams.update({
    'font.family': 'latin modern math', 'font.size':12, 'mathtext.fontset':'stix',
    'axes.formatter.use_mathtext': True, 'axes.formatter.limits': (-4, 4),
    'axes.grid': True, 'grid.color': 'gray', 'grid.linewidth': 0.5,
    'xtick.top': True, 'ytick.right': True,
    'xtick.direction': 'inout', 'ytick.direction': 'inout',
    'xtick.minor.size': 2.0, 'ytick.minor.size': 2.0,  # default 2.0
    'xtick.major.size': 4.0, 'ytick.major.size': 4.0,  # default 3.5
    'xtick.minor.visible': True, 'ytick.minor.visible': True
})
from matplotlib import gridspec, rcParams, rc
from matplotlib.widgets import Cursor

from astropy.table import Table, Column
from astropy.io import fits
from astropy.stats import sigma_clip, gaussian_fwhm_to_sigma
from astropy.modeling.models import Gaussian1D, Chebyshev2D
from astropy.modeling.fitting import LevMarLSQFitter
from astropy.nddata import CCDData

from photutils.aperture import RectangularAperture, aperture_photometry

from skimage.feature import peak_local_max
import _tool_visualization as vis



def disable_mplkeymaps():
    rc('keymap',
       fullscreen='',
       home='',
       back='',
       forward='',
       pan='',
       zoom='',
       save='',
       quit='',
       grid='',
       yscale='',
       xscale='',
       all_axes=''
       )
FONTSIZE = 12 # Change it on your computer if you wish.
rcParams.update({'font.size': FONTSIZE})

#%%
CALPATH = Path('../../Tutorial_Data/45P/calib')
REDPATH = Path('../../Tutorial_Data/45P/reduced')
DISPAXIS = 1  # 1 = line = python_axis_1 // 2 = column = python_axis_0

COMPIMAGE = CALPATH/"comp.fits" # Change directory if needed!
OBJIMAGE  = DATAPATH/'mls170213_0025.fits'
LINE_FITTER = LevMarLSQFitter()

# Parameters for IDENTIFY
FITTING_MODEL_ID = 'Chebyshev'
ORDER_ID = 4
NSUM_ID = 10
FWHM_ID = 4 # rough guess of FWHM of lines in IDENTIFY (pixels)

# Parameters for REIDENTIFY
FITTING_MODEL_REID = 'Chebyshev' # 2-D fitting function
ORDER_SPATIAL_REID = 6
ORDER_WAVELEN_REID = 6
STEP_REID = 15  # Reidentification step size in pixels (spatial direction)
NSUM_REID = 10
TOL_REID = 5 # tolerence to lose a line in pixels

# Parameters for APALL (sky fitting and aperture extract after sky subtraction)
## parameters for finding aperture
NSUM_AP = 10
FWHM_AP = 10
STEP_AP = 10  # Recentering step size in pixels (dispersion direction)
## parameters for sky fitting
FITTING_MODEL_APSKY = 'Chebyshev'
ORDER_APSKY = 3
SIGMA_APSKY = 3
ITERS_APSKY = 5
## parameters for aperture tracing
FITTING_MODEL_APTRACE = 'Chebyshev'
ORDER_APTRACE = 3
SIGMA_APTRACE = 3
ITERS_APTRACE = 5
# The fitting is done by SIGMA_APTRACE-sigma ITERS_APTRACE-iters clipped on the
# residual of data.

#%%
lamphdu = fits.open(COMPIMAGE)
objhdu = fits.open(OBJIMAGE)
lampimage = lamphdu[0].data
objimage  = objhdu[0].data

if lampimage.shape != objimage.shape:
    raise ValueError('lamp and obj images should have same sizes!')

if DISPAXIS == 2:
    lampimage = lampimage.T
    objimage = objimage.T
elif DISPAXIS != 1:
    raise ValueError('DISPAXIS must be 1 or 2 (it is now {:d})'.format(DISPAXIS))

EXPTIME = objhdu[0].header['EXPTIME']
OBJNAME = objhdu[0].header['OBJECT']
# Now python axis 0 (Y-direction) is the spatial axis
# and 1 (X-direciton) is the wavelength (dispersion) axis.
N_SPATIAL, N_WAVELEN = np.shape(lampimage)
N_REID = N_SPATIAL//STEP_REID # No. of reidentification
N_AP = N_WAVELEN//STEP_AP # No. of aperture finding

# ``peak_local_max`` calculates the peak location using maximum filter:
#   med1d_max = scipy.ndimage.maximum_filter(med1d, size=10, mode='constant')
# I will use this to show peaks in a primitive manner.
MINSEP_PK = 5   # minimum separation of peaks
MINAMP_PK = 0.01 # fraction of minimum amplitude (wrt maximum) to regard as peak
NMAX_PK = 50
print("setting done!")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 69
     66 DISPAXIS = 1  # 1 = line = python_axis_1 // 2 = column = python_axis_0
     68 COMPIMAGE = CALPATH/"comp.fits" # Change directory if needed!
---> 69 OBJIMAGE  = DATAPATH/'mls170213_0025.fits'
     70 LINE_FITTER = LevMarLSQFitter()
     72 # Parameters for IDENTIFY

NameError: name 'DATAPATH' is not defined

Python Example of Longslit Spectroscopy

Contents

3. Python Example of Longslit Spectroscopy#

3.1. About the Data#

3.2. Postprocessing#

3.3. Preparation#