Spectra Extraction

The Extraction Table

The ExtractionTable is the final product of the Inverse Dispersion Solution data tuning and it is the key point to perform the spectra extraction. Given the mask and the instrument configuration, that table contains the all information to achieve the 2D spectra extraction. This table is stored as a FITS table.

As described in the previous sections, during IDS data tuning, each slit is sliced and the solution of each slice is computed independently by the others.

If we focus on a single line, the IDS solution for this line can be jagged and not as smooth as in the frame appears. It can also happen that the computation of the solution for some row fails (cosmic, bad pixels or other effects can affect the results on that row). Here we can see the effect previously described, in blue the expected positions.


To avoid these problems SpectraPy provides a methods to rectify the solutions contained in the Extraction Table and fill any gaps in case of row failures.

Extraction Table rectification

This rectification process is performed slit by slit, since the Extraction Table handles each slice separately. For each slit there are sets of coefficients which are slightly different slice by slice.

Focusing on a single slit, and fixed the order of the polynomial coefficient, SpectraPy gathers together the coefficients of all the slices. For each coefficient group, SpectraPy fits a polynomial in order to find a smoother variation of the current coefficient. The final result, is to have a smoother overall solution for the entire slit.

In the figure below for a given slit, we can see in detail the values of the 0th coefficient slice by slice; the coefficient of each slice is slightly different from the value of the next slice (the orange points). The red curve is the fit along the slits (slices position are in arbitrary units). The evaluation of the fit on each slice positions gives us the new 0th coefficients (the blue points).


Here the lines of code to rectify the solutions

>>> from spectrapy.datacalib.extractiontable import ExtTable
>>> exr = ExtTable.load("examples/tmp/ID532016.exr")
>>> exr.rectify(deg=2, margin=1)

In this case we choose a polynomial of 2nd order (parameter deg=2) and we excluded the 1 pixel (both sides) at the edges of the slit (margin=1), this is done just to avoid possible mask cutting issues.

Or, if we prefer, we can just refit the single slit

>>> exr['121'].rectify(deg=2, margin=2)

If we display the result of the rectified Extraction Table we can see the solution on a single line is smoother than before (the blue crosses before rectification, the red crosses after rectification).


Check the Wavelength Solution

Even if the wavelegenth solution contained in the Extraction Table is very precise, this solution could be not the best once applied on some science images. For example, this can append because distortions could be different night by night, or occasionally mask can slightly drift away along the night. In order to achieve these kind of problems, or at least to aware of them, the current SpectraPy release provides us one class to check and inspect the quality of the wavelength solution applied on a given image.

These checks are performed by the WavelengthCalibCheck class.

We must initialize the class with the Extraction Table to check, load the desiderd image and run the check_ext_table method. This method can be run, applying multiple constraints:

  • the input catalog: to select the lines to check

  • the wavelength range: to limit the lambda range to check

  • the slits IDs: to select only some slits instead of the whole mask

  • the rows range: to check just a region of the slit

>>> from spectrapy.check.wavelengthcheck import WavelengthCalibCheck
>>> exttable = "examples/tmp/ID532016.exr"
>>> arc_frame = "examples/data/mods1r/mods1r.20180121.0073.fits.bz2"
>>> catalog = "conf/catalogs/NeHg_hr.dat"
>>> check = WavelengthCalibCheck(exttable)
>>> check.load_image(arc_frame)
>>> results = check.check_ext_table(catalog, wrange=(6000., 9000.), slits=('49', '54'), row_start=10, row_end=20)
In this example we decide to:
  • check only slits 49 and 54

  • use catalog lines in the range between 6000 and 9000 Angstroms

  • check just the slits rows between 10 and 20

For each selected slit and row, the check_ext_table method starts from the expected positions, given by the Extraction Table. then it computes the real line position around these expect positions. This is done for each available line the catalog.

This method collects returns in a class, which allow us to display and browse these results.

Here below, how to display results

>>> results.plot(legend=True)

The class shows us a matplotlib figure containing 3 axes. In the main frame, the upper one, we can see the displacements between expected line positions and local measurement of the line position. In case the legend parameter is True, a legend with the slits ID is showed on the left side of the frame. This legend allow us to browse along slits, selecting the desired slit, clicking on the blue bullets. We can also browse along slits pressing n key for the next slit and p for the previous.

We can also inspect more details of the solution. Clicking with the left mouse button on a single marker of the main frame, the solutions related the a single row will be displayed in the bottom left panel. In case we want focus on the quality of the solution for a single line we can click with the right mouse button on a marker of that line, and the line solution will be displayed on the bottom right panel.


Spectra Extraction

The ExtractionTable table contains all the information to extract spectra for a given mask and instrument configuration.

Applying the ExtractionTable on data, we obtain a multi-extension FITS file, which contains the 2D spectra wavelength calibrated and corrected by optical distortions. In case variance (or errors) on raw data are provided, the re-binning procedure provides also the variance on the re-binned spectra. For each rebinned pixel a quality rate is computed; this is a real number from 0 to 1 which gives us the percentage of pixels used during the rebinning procedure. These rates are stored in the multi extension as well the variance and the data.

To perform the extraction we must initialize a re-binning engine

>>> from spectrapy.extraction.exponentialfilter import ExponentialFilter
>>> engine = ExponentialFilter(2, 1000)
The current implementation of SpectraPy provides the Exponential filter as re-binning engine, which requires:
  • the radius of the resampling kernel (2 pixels in this case)

  • the sub-pixels accuracy in the re-sampling for each pixel (1000 by default)

The resampling procedure requires a lot of CPU time and the extraction can take a very long time. For this reason beside the pure Python class we develop also a Cython class which is faster. In the Cython class the radius and the sub-pixels parameters can not be tuned, the class uses 2 pixels as radius and 1000 as subpixels. This choice is quite standard.

In case, for any reason, we want change these parameters, we must set the cython flag to False during the class initialization.

The re-sampling engine must also know the extraction range we are interested in, and the re-binning step of the re-sampled spectra. In this example we decided to extract from 6000 up to 9000 Angstrom, creating spectra linearly re-sampled every 0.8 Angstrom.

>>> engine.set_extraction_range(6000., 9000., 0.8)


A good choice for resampling value (0.8 Angstrom in this example), could be the nominal grism dispersion value.

Now, we can do the real job starting the extraction. By default the extraction will extract all the spectra described in the ExtractionTable, but we can also decide to extract a sub-sample of them providing a list of ID, as defined in the Mask description file.

>>> from spectrapy.datacalib.extractiontable import ExtTable
>>> #Load the extraction table and rectify it
>>> exr=ExtTable.load('examples/tmp/ID532016.exr')

>>> #Extract the spectra
>>> objlist = ('121', '119', 'gnz_lae1_1')
>>> mods_arc = "examples/data/mods1r/mods1r.20180121.0073.fits.bz2"
>>> spectra2d = engine.extract(exr, mods_arc, objlist)

and save the spectra.

>>> spectra2d.writeto('examples/tmp/mods1r_spectra.fits')

This produces a multi extension FITS file the file in figure below, where extension name is combination of: slit id, data content (in the current version only 2DCOUNTS is available) and type of data (data, variances of quality flags).


The selection can be performed also on a sub region along the cross dispersion direction. For example in the longslit LUCI case, we can decide to extract just the region around the object itself.

>>> from spectrapy.extraction.exponentialfilter import ExponentialFilter
>>> from spectrapy.datacalib.extractiontable import ExtTable

>>> exr=ExtTable.load('examples/tmp/LS075.exr')
>>> engine = ExponentialFilter()

>>> # Set the extraction range
>>> engine.set_extraction_range(15000., 22000., 4.3)

>>> # Select a sub regione along the spatial direction
>>> sc_file = "examples/data/luci1LoRes/luci1.20180202.0181.fits.bz2"
>>> spectra2d = engine.extract(exr, sc_file, row_start=250, row_end=310)

>>> #And save it
>>> spectra2d.writeto('examples/tmp/luci1_spectra.fits', overwrite=True)

The data content of the extracted spectra appears like this


Spectra Extraction Adjustment

Models are computed on a set of frames and then applied on the science frames (usually another set of frames) to perform the spectra extraction. In some cases the match between the 2 set of data can be perfect.

Data can show tiny differences along the night, due to different instrument distortions at different positions of the telescope. Or due to slight shifts of the mask along the night.

In order to compensate these changes, SpectraPy allow us to adjust the solution contained in the ExtractionTable on the data we are extracting. The extract method of the ExponentialFilter class has 2 optional parameters to perform this adjustment: xadjust and lines.

If the xadjust parameter is True the recipe recomputes both left and right edges, at the lambda reference position, of the spectra selected by the objlist. Then it computes the difference between these new edge positions and the edge positions expected by the ExtractionTable. The shift along the cross dispersion direction is the sigma clipped mean of these differences.

The parameter lines is used to compute the shift along the dispersion direction. If this parameter is a valid line catalog containing just few very bright isolated lines, SpectraPy uses these lines to compute the differences between theirs expected positions (according with the ExtractionTable solution) and theirs real position on the frame. The average the differences between expected and real positions is the shift amount along the dispersion direction.

Here just an example of the use of these parameters

>>> spectra2d_adj = engine.extract(exr, sc_file, row_start=250, row_end=310, xadjust=True, lines=[16235.376, 17880.298, 19350.119, 21802.312])