SPECTRA module

spectra is a module that implements some classes/tools for handling basic spectral data reduction: tracing and extraction of spectra, wavelength calibration, and flux calibration. The routines allow for both longslit data, multiobject data, echelle data, and a combination (multiple orders with longslit, or multiple slitlets).

Data are read and processed using pyvista Data structures, which carry along both the image data and an associated uncertainty array, and allow wavelength, response, and sky attributes to be added.

The routines are implemented with three basic classes, the Trace object, the WaveCal object, and the FluxCal. In more detail:

Trace

The Trace object is used to store functional forms and locations of spectral traces. It supports multiple traces per image, e.g., in a long-slit spectrograph, a multiobject spectrograph (fiber or slitlets), or a multiorder spectrograph.

Trace objects can be defined from one image and used on other images, allowing for small global shifts of the traces in the spatial direction. They can be saved and read from disk as FITS files.

To create a trace from scratch, instantiate a trace object with rows= (to specify window, i.e. length of slit), lags= (to specify how far search can extend for translation of trace), degree= (to specify degree of polynomial for fit to trace), and rad= (to specify a radius in pixels to use for centroiding a trace). Internally, pyvista spectral routines will assume that wavelength runs horizontally, i.e. across columns, so if your spectrograph has vertical spectra, set transpose=True in the Trace object.

Alternatively, you can load a saved Trace object from disk by instantiating using the file= keyword (use file=’?’ to see a listing of files distributed with the package).

To make the trace, the code will generally start at the center of the image in the wavelength direction (unless the sc0 attribute is set to another value), centroid the spectrum there, then work in both directions to get a centroid each pixel (or every skip pixels if the skip= keyword is given); the centroid for each wavelength is used as a starting guess for the subsequent wavelength if the peak is above a S/N threshold), Finally, a function will be fit to the derived centroids and the model (currently astropy Polynomial1D) wll be stored in the trace structure.

To make a trace manually use the trace(spectrum,row) method, where row is the starting guess for the spatial postion at the center of the spectrum (or sc0). You can use the findpeak() method to find object(s) automatically and pass the returned value(s) to trace(). Alternatively, you can use find() which will cross-correlate a saved trace crossection and identify a location from the highest peak; if you specify inter=True with find(), then the image will be displayed (display= must also be specified to give the TV object to use), and you can mark the location interactively. If row is a list (i.e. for multiple objects or orders), then multiple traces will be made.

For spectrographs where the traces are relatively stable in location (e.g. multiple orders), a saved model can be convenient to use to remake traces. To retrace, use the retrace(spectrum) method. This will do a cross-correlation of a saved cross-section with the input spectrum to determine a shift, then will use that as a starting location to retrace.

The default method for finding the trace locations is to use a centroid. If you add gaussian=True to trace(), then the center will be determined using a Gaussian fit to the profile. This will also determine the derived Gaussian widths and fit a polynomial model to these, which is required if one wants to do optimal extraction using a Gaussian profile.

If you pass a pyvista TV object with the display= keyword, the display will show the location of the calculated centers, and the polynomial fit to these: centers in red are those rejected during the fiting process. You will need to hit the space bar in the display window to continue.

Extraction is done using the extract(spectrum,rad=,[back=[[b1,b2],[b3,b4]]) method, which uses boxcar extraction with the specified radius.

Optionally, a background value as determined from one or more background windows can be subtracted, where the window locations are specified in pixels relative to the object trace position. Note that if there is non-negligible line curvature that this can lead to poor subtraction of sky emission lines. In this case, you might want to determine a 2D wavelength solution (see below), and resample the sky spectra to the wavelength scale of the object (or to some other wavelength scale) before subtraction.

If you pass a pyvista TV object with the display= keyword, then the display will show the location of the extraction and background windows (if any). You will need to hit the space bar in the display window to continue.

Attributes

  • type : type of astropy model to use for trace shape (currently, just Polynomial1D)

  • model : array of astropy models for each trace

  • degree : polynomial degree for model fit

  • rad : radius in pixels to use for calculating centroid

  • sc0 : starting column for trace, will work both directions from here

  • spectrum : reference spatial slice at sc0

  • pix0 : derived shift of current image relative to reference

  • lags : array of lags to try when finding object location

  • transpose : boolean for whether axes need transposing to put spectra along rows

  • rows : range of slit length

  • index : identification for each trace

Methods (for details, see docstrings below)

  • find()

  • findpeak()

  • retrace()

  • trace()

  • extract()

  • extract2d()

WaveCal

Wavelength calibration proceeds as follows:

If there is a previous WaveCal object with wavelength calibration for this instrument, that can be used to facilitate wavelength calibration. A previous WaveCal object can be saved with the write() method. What is saved is the type and degree of the wavelength model, a list of pixel and wavelengths of identified lines, and the spectrum from which these were identified. A WaveCal object can be instantiated using the file= keyword to read in a previous WaveCal (use file=’?’ to see a listing of files distributed with the package).

Given this WaveCal instantiation, identify() is called to identify and fit a wavelength solution. This is achieved as follows:

  1. if input wav array/image is specified, use this to identify lines

  2. if WaveCal object has an associated spectrum, use cross correlation to identify shift of input spectrum, then use previous solution to create a wavelength array. Cross correlation lags to try are specified by lags=range(dx1,dx2), default range(-300,300)

  3. if inter==True, prompt user to identify 2 lines

  4. use header cards DISPDW and DISPWC for dispersion and wave center or as specified by input disp=[dispersion] and wref=[lambda,pix]

Given the wavelength guess array, identify() will identify lines from input file of lamp/reference wavelengths, or, if no file given, the lines saved in the WaveCal structure.

Lines are identified by looking for peaks within rad pixels of initial guess that have S/N exceeded the threshold given by thresh= keyword (default 100). Pixel centers of the lines are determined by a Gaussian fit to the line; the derived line widths are saved in a fwhm attribute along with the line centers.

After line identification, fit() is automatically called, unless fit=False. During the fit, if plot=True, the user can remove discrepant lines (use ‘n’ to remove line nearest to cursor in wavelength, ‘l’ to remove all lines to left, ‘r’ to remove all lines to right); removal is done by setting the weights of these lines to zero in the WaveCal structure. After removing lines to get a better initial solution, it may be desirable to re-enable these lines to see if they can be more correctly identified given a better initial wavelength guess.

An example, given an extracted spectrum spec, might proceed as follows:

wav=spectra.WaveCal('KOSMOS/KOSMOS_blue_waves.fits')
wav.identify(spec,plot=True)
# if you removed a bunch of lines, especally at short and long wavelength ends,
#   you can try to recover them using your revised solution:
wav.weights[:] = 1.
wav.identify(spec,plot=True,lags=range(-50,50))

If there is no previous WaveCal, a new WaveCal can be instantiated, specifing the type and degree of the model. Lines are identified given some estimate of the wavelength solution, either from an input wavelength array (wav= keyword) or from a [wave,pix] pair plus a dispersion (wref= and disp=) keyword, as described above. This solution will be used to try to find lines as specified in an input reference line file (file=) : a centroid around the input position of width given by the rad= keyword is computed.

One you have an acceptable wavelength solution, it can be used to create a wavelength array and transferred to the wave attribute of an object using the add_wave() method.

The WaveCal object also includes a skyline() method, which can be used to adjust the wavelength solution for an object based on sky emission lines. The object should already have a wavelength solution, e.g. from arcs, attached. The skyline() routine will use this to look for skylines and then redo the wavelength solution, but only allowing for the 0th order term to change, i.e. a constant wavelength shift. If multiple sky lines across a broad wavelength range are available, you can also fit for a dispersion change if you set linear=True.

FluxCal

The FluxCal object allows for the creation of a spectral response curve given observations of one or more flux standard stars. It also includes a spectral extinction correction using a set of mean extinction coefficients as a function of wavelength. The derived response curve can be applied to an input spectrum to provide a flux-calibrated spectrum. The accuracy of the flux calibration may be limited by differential refraction effects in either the calibration stars or the object itself. The accuracy of the absolute calibration is further limited by light losses outside the slit.

A FluxCal object is instantiated, specifying either the polynomial degree (degree=) to be used for the response curve fit, or else you can use a median or mean of the observed flux standards (degree=-1). Observed 1d (i.e., extracted) spectra of spectrophotometric stars are added using the addstar(spectrum, waves, file=) method, where the input file gives the true fluxes as a function of wavelength; libraries of standard star spectra from ESO are included. Alternatively, you can specify the flux of the standard star using the stdflux= keyword, passing an astropy Table with (at least) columns wave, flux, and bin. The observed spectrum is corrected for extinction and saved, along with the true spectrum.

Once or or more stars have been loaded, a response curve is derived using the response() method. If degree>=0 is used, a polynomial is fit to the logarithm of the ratio of the extinction-corrected observed to true fluxes, allowing for an independent scale factor for each star. If degree<0, then a median of the flux curves is used, or a mean if mean=True is passed to response(). A median is probably only appropriate if multiple flux standards have been observed in a single observation (i.e., a multiobject spectrograph), otherwise a mean is probably more useful, since it is probably unlikely that all of the individual star response curves will be at the same level, given variations in throughput from exposure to exposure and the assumption of a mean extinction curve.

Given the response curve, input spectra can be corrected for extinction and instrumental response using the correct() method.

spectra functions

class pyvista.spectra.WaveCal(file=None, type='chebyshev', degree=2, ydegree=2, pix0=0, index=0, hdu=1, orders=None)[source]

Class for wavelength solutions

Parameters

  • file (str, optional) – filename for FITS file with WaveCal attributes

  • type (str, optional, default='chebyshev') – astropy model function

  • degree (int, optional, default=2) – polynomial degree for wavelength direction

  • ydegree (int, optional, default=2) – polynomial degree for spatial or cross-dispersed direction

  • pix0 (int, optional, default=0) – reference pixel

Variables

  • type (str) – astropy model function for wavelength solution

  • degree (int) – polynomial degree for wavelength direction

  • ydegree (int) – polynomial degree for spatial or cross-dispersed direction

  • pix0 (int) – reference pixel

  • orders (array_like) – spectral order for each row in spatially resolved or cross-dispersed

  • pix (array_like) – pixel positions of identified lines in spectrum

  • waves (array_like) – wavelength positions of identified lines in spectrum

  • weights (array_like) – weights for fitting of identified lines in spectrum

  • y (array_like) – spatial or cross-dispersed array position

  • spectrum (array_like) – reference spectrum to be used to identify lines

  • model (astropy model, or list of models) – Model(s) for the wavelength solution

write(file, append=False)[source]

Save object to file

Parameters

  • file (str, name of output file to write to, FITS format) –

  • append (bool, append to existing file (untested)) –

wave(pixels=None, image=None, domain=False)[source]

Wavelength from pixel using wavelength solution model

With pixels=[pixels,rows] keyword, return wavelengths for input set of pixels/rows With image=(nrow,ncol), returns wavelengths in an image

Parameters

  • pix (array_like, optional) – input pixel positions [x] or [y,x]

  • image (tuple, optional) – for input image size [nrows,ncols], return wavelengths at all pixels

Return type

wavelength

add_wave(hd)[source]

Add wavelength attribute to input image using current wavelength solution

Parameters

hd (Data object) – Image to add wavelength array to

getmod()[source]

Return model for current attributes

plot(hard=None)[source]

Plot current solution

dispersion()[source]

approximate dispersion from 1st order term”

fit(degree=None, reject=3, inter=True)[source]

do a wavelength fit

If a figure has been set in identify, will show fit graphically and allow for manual removal of lines in 1D case. In 2D case, outliers are detected and removed

Parameters

degree (int, optional) – degree of polynomial in wavelength, else as previously set in object

set_spectrum(spectrum)[source]

Set spectrum used to derive fit

get_spectrum()[source]

Set spectrum used to derive fit

identify(spectrum, sky=False, file=None, wav=None, wref=None, inter=False, orders=None, verbose=False, rad=5, thresh=100, fit=True, maxshift=10000000000.0, disp=None, display=None, plot=None, pixplot=False, domain=False, plotinter=True, xmin=None, xmax=None, lags=range(-300, 300), nskip=None, rows=None)[source]

Given some estimate of wavelength solution and file with lines, identify peaks and centroid, via methods:

  1. if input wav array/image is specified, use this to identify lines

  2. if WaveCal object as associated spectrum, use cross correlation to identify shift of input spectrum, then use previous solution to create a wavelength array. Cross correlation lags to try are specified by lags=range(dx1,dx2), default range(-300,300)

  3. if inter==True, prompt user to identify 2 lines

  4. use header cards DISPDW and DISPWC for dispersion and wave center or as specified by input disp=[dispersion] and wref=[lambda,pix]

Given wavelength guess array, identify lines from input file, or,

if no file given, lines saved in the WaveCal structure

Lines are identified by looking for peaks within rad pixels of

initial guess

After line identification, fit() is called, unless fit=False

With plot=True, plot of spectrum is shown, with initial wavelength

guess. With pixplot=True, plot is shown as function of pixel

skyline(hd, plot=True, thresh=50, inter=True, linear=False, file='skyline.dat', rows=None)[source]

Adjust wavelength solution based on sky lines

Parameters

  • hd (Data object) – input pyvista Data object, must contain wave attribute with initial wavelengths

  • plot (bool, default=True) – display plot results

  • thresh (float, default=50) – minimum S/N for line detection

  • rows (array-like, default=None) – if specified, only use specified rows for sky spectrum, relevant for 2D correction to ignore object rows

  • inter (bool, default=True) – allow for interactive removal of lines

  • linear (bool, default=False) – if True, allow for dispersion to be ajusted as well as wavelength zeropoint requires at least two sky lines!

  • file (str, default='skyline.dat') – file with sky lines to look for, if you want to override default

scomb(hd, wav, average=True, usemask=True)[source]

Resample onto input wavelength grid

Uses current wavelength solution, linearly interpolates to specified

wavelengths, on a row-by-row basis. Allows for order overlap.

Parameters

  • hd (array or CCDData) – input image to resample

  • wav (array_like) – new wavelengths to interpolate to

  • average (bool, optional, default=True) – if overlapping orders, average if True, otherwise sum

  • usemask (bool, optional, default=True) – if True, skip input masked pixels for interpolation

correct(hd, wav)[source]

Resample input image to desired wavelength scale

Uses current wavelength solution, linearly interpolates to specified

wavelengths, on a row-by-row basis.

Parameters

  • hd (Data, input image to resample) –

  • wav (array_like, new wavelengths to interpolate to) –

class pyvista.spectra.Trace(file=None, inst=None, type='Polynomial1D', degree=2, sigdegree=0, pix0=0, rad=5, spectrum=None, model=None, sc0=None, rows=None, transpose=False, lags=None, channel=None, hdu=1)[source]

Class for spectral traces

Variables

  • type (str) – type of astropy model to use

  • degree (int) – polynomial degree to use for trace

  • sigdegree (int) – polynomial degree to use for fitting gaussian sigma trace width

  • sc0 (int) – starting column for trace, will work in both directions from here

  • pix0 (int) – derived shift of current image relative to reference image

  • spectrum (array_like) – reference spatial slice at sc0, used to determine object location

  • rad (int) – radius in pixels to use for calculating centroid

  • lags (array_like) – range of lags to use to try to find object locations

Parameters

file (str, optional) – filename for FITS file with Trace attributes

write(file, append=False)[source]

Write trace information to FITS file

trace(im, srows, sc0=None, plot=None, display=None, rad=None, thresh=20, index=None, skip=10, gaussian=False, verbose=False)[source]

Trace a spectrum from starting position

Parameters

  • im (Data) – input image

  • srows (array-like) – location(s) at sc0 for initial trace location(s) guess(es)

  • rad (float, optional, default=self.rad) – radius of window to use to find trace locations

  • index (integer, optional, default=None) – index to label trace(s) with

  • skip (integer, optional, default=10) – measure trace center every skip pixels, using median of data from -skip/2 to skip/2

  • gaussian (bool, optional, default=False) – if True, use gaussian fit for trace location instead of centroid. with gaussian=True, will also fit trace widths into sigmodel, with polynomial of degree self.sigdegree

  • sc0 (integer, optional, default=ncol/2) –

  • plot (bool, optional, default=None) –

  • display (TV object, optional, default=None) –

plot()[source]

Plots points and fits for traces

retrace(hd, plot=None, display=None, thresh=20, gaussian=False, skip=10, rad=3)[source]

Retrace starting with existing model

findpeak(hd, sc0=None, width=100, thresh=50, plot=False, sort=False, back_percentile=10, method='linear', smooth=5, diff=10000, bundle=10000, verbose=False)[source]

Find peaks in spatial profile for subsequent tracing

Parameters

  • hd (Data object) – Input image

  • sc0 (int, default=None) – pixel location of wavelength to make spatial profile around if none, use sc0 defined in trace

  • width (int, default=100) – width of window around specfied wavelength to median to give spatial profile

  • thresh (float, default = 50) – threshold for finding objects, as a factor to be multiplied by the median uncertainty

  • smooth (float, default = 5) – smoothing FWHM (pixels) for cross-section before peak finding

  • sort (bool, default=False) – return peaks sorted with brightest first

Returns

tuple – peak locations can be passed to trace()

Return type

list of peak locations, and list of indices

find(hd, sc0=None, width=100, lags=None, plot=None, display=None, inter=False, rad=3)[source]

Determine shift from existing trace to input frame

Parameters

  • hd (Data object) – Input image

  • width (int, default=100) – width of window around central wavelength to median to give spatial profile

  • lags (array-like, default=self.lags) – range of cross-correlation lags to allow

  • rad (int, default=3) – radius around xcorr peak to do polynomial fit to

  • display (pyvista.tv object, default=None) – if not None, tv object to display in

immodel(im, ext, threads=0)[source]

Create model 2D image from input fluxes

extract(im, rad=None, back=[], fit=False, old=False, display=None, plot=None, medfilt=None, nout=None, threads=0)[source]

Extract spectrum given trace(s)

Parameters

  • hd (Data object) – Input image

  • rad (float, default=self.rad) – radius for extraction window

  • back (array-like of array-like) – list of two-element lists giving start and end of background window(s), in units of pixels relative to trace location

  • nout (integer, default=None) – used for multi-object spectra. If not None, specifies number of rows of output image; each extracted spectrum will be loaded into indices loaded into index attribute, with an index for each trace

extract2d(im, rows=None, plot=None, display=None, buffer=0)[source]

Extract 2D spectrum given trace(s)

Assumes all requested rows uses same trace, just offset, not a 2D model for traces. Linear interpolation is used.

findslits(im, smooth=3, thresh=1500, display=None, cent=None, degree=2, skip=50, sn=False)[source]

Find slits in a multi-slit flat field image

findslits() attempts to find slit locations by looking for peaks in the derivative of the flux (or, if sn=True, in the derivative of S/N) in a window of skip//2 pixels around the center of the image (or around the location given by cent=). Taking those values, it then looks for corresponding edges moving to the left and right of the center. Finally, it fits a polynomial to each set of slit edges to populate the model attribute of the Trace object. It also populates the rows attribute of the Trace object with bottom and top locations of the slit in the center of the detector.

Parameters

  • data (array or Data object, wavelength dimension horizontal) – input data to find slit edges, usually a flat field

  • smooth (float, optional, default=3) – smoothing radius for calculated derivatives

  • sn (boolean, optional, default=False) – if True, look for edges in delta(S/N)

  • thresh (float, optional, default=1500) – threshold for detecting edges in delta(signal)

  • skip (integer, optional, default=50) – spacing of where along spectra to identify edges

  • display (TV object, optional) – TV object to display derivatives and slit locations

  • cent (int, optional, default=None) – location of center of spectra, if None use chip center

  • degree (int, option, degree=2) – polynomial degree to use to fit edge locations

class pyvista.spectra.FluxCal(degree=3, median=False, extinct='flux/apo_extinct.dat')[source]

Class for flux calibration

Parameters

  • degree (int, default=3) – Order of polynomial for response curve, -1 for median/mean

  • median (bool, default=False) – if degree<0, use median response curve rather than mean

  • extinct (Table/str, default='flux/apo_extinct.dat') – mean extinction curve to use, can be input as astropy Table with columns [‘wave’,’mag’], or a filename from which these can be read

write(file, append=False)[source]

Save object to file UNTESTED, and NO READ yet!

Parameters

  • file (str, name of output file to write to, FITS format) –

  • append (bool, append to existing file (untested)) –

extinct(hd, wave)[source]

Correct input image for atmospheric extinction

Parameters

  • hd (Data object with data) – Input image

  • wave (float, array-like) – Wavelengths for input image (separate from attribute to allow slicing in hd)

  • file (str, default='flux/apo_extinct.dat') – Two column file ([‘wave’,’mag’]) with extinction curve

addstar(hd, wave, pixelmask=None, file=None, stdflux=None, exptime='EXPTIME', extinct=True, badval=None)[source]

Derive flux calibration vector from standard star

Parameters

  • hd (Data object with standard star spectrum) – Input image

  • file (str, optional) – File with calibrated fluxes, with columns [‘wave’,’flux’,’bin’], must be readable by astropy.io.ascii with format=’ascii’

  • stdflux (astropy Table, optional) – Table with calibrated fluxes in columns ‘wave’,’flux’,’bin’

  • extinct (bool, default=True) – Use mean extinction curve to correct for atmospheric extinction

  • exptime (str, default='EXPTIME') – Card name to use for exptime normalization, set to None for no norm

response(degree=None, inter=False, plot=True, legend=True, hard=None, medfilt=None)[source]

Create response curve from loaded standard star spectra and fluxes

Parameters

  • degree (integer, default=self.degree) – polynomial degree

  • medfilt (integer, default=None) – width of median filter to apply to mean/median response curve

  • plot (bool, default=False) – set to True to see plot

  • legend (bool, default=True) – label stars on plot with a legend

  • hard (str, default=None) – file name for hardcopy plot

correct(hd, waves, extinct=True, exptime='EXPTIME')[source]

Apply flux correction to input spectrum

Parameters

  • hd (Data object with spectrum to be corrected) – Input image

  • waves (array-like) – Wavelength array for hd (separate from hd to allow slicing of hd)

  • extinct (bool, default=True) – Apply atmospheric extinction correction

  • exptime (str, default='EXPTIME') – Card name to use for exptime normalization, set to None for no norm

pyvista.spectra.gfit(data, x0, rad=10, sig=3, back=False)[source]

Fit 1D gaussian

pyvista.spectra.gauss(x, *p)[source]

Gaussian function

pyvista.spectra.findpeak(x, thresh, diff=10000, bundle=0, verbose=False, sort=False)[source]

Find peaks in vector x above input threshold attempts to associate an index with each depending on spacing

Parameters

  • x (float, array-like) – input vector to find peaks in

  • thresh (float) – threshold for peak finding

  • diff (int) – maximum difference in pixels between traces before incrementing fiber index

  • bundle (int) – number of fibers after which to allow max distance to be exceeded without incrementing

  • sort (bool, optional, default=False) – return peak locations and indices sorted with brightest first

pyvista.spectra.fitpeak(data, peaks, smooth=3, width=3, desc=False)[source]

Find slit edges given input set of locations

Parameters

  • data (1D array) – input data to find slit edges, usually a flat field

  • smooth (integer, optional, default=3) – smoothing radius for derivative calculation

  • width (integer, optional, default=3) – window around peak of derivative to do polynomial fit for edge

  • desc (book, optional, default=False) – if True, find descending peaks

pyvista.spectra.getinput(prompt, fig=None, index=False)[source]

Get input from terminal or matplotlib figure

pyvista.spectra.model_col(pars)[source]

Extract a single column, using boxcar extraction for multiple traces

pyvista.spectra.extract_col_fit(pars)[source]

Extract a single column, using boxcar extraction for multiple traces

pyvista.spectra.extract_col_old(pars)[source]

Extract a single column, using boxcar extraction for multiple traces

pyvista.spectra.extract_col(pars)[source]

Extract a series of columns, using boxcar extraction for multiple traces

pyvista.spectra.chebyshev2D_shift(x, y, *coeff, xdegree=0, ydegree=0, xdomain=[-1, 1], ydomain=[-1, 1], apers=None)[source]

Evaluate chebyshev2d function