5.1 Data Reduction and Analysis Applications
Most of the software tools for operating on COS FITS files are contained in three Python packages that come with the stenv python distribution: https://stenv.readthedocs.io
The three Python packages used with COS are:
- calcos: Contains the calcos package, used for processing data through the calibration pipeline. Section 3.6 describes how to run calcos.
- costools: Contains COS-specific tasks to calibrate and interpret data. Many of these tasks are described in Chapter 3. A complete listing is presented in Section 5.1.1.
- astropy: Contains packages and modules to read and operate on FITS tables, which is the format used for many COS products including 1-D spectra and most COS reference files. We provide specific examples of their use here in Section 5.1.2.
We also present a brief summary of spectral display and analysis tasks in Section 5.1.3. To run the examples in Sections 5.1.1 to 5.4, one must activate the stenv environment before starting Python.
In addition to what is below, users are encouraged to use the Jupyter Notebooks for COS. Those available include:
| Name | Topic |
|---|---|
| Setup | Setting up your computer environment to work with COS data |
| DataDI | Downloading COS Data from the archive |
| ViewData | Beginning to work with COS data in Python: plotting, binning, calculating SNR, & evaluating a spectrum |
| AsnFile | Modifying or creating an association file |
| CalCOS | Running the COS pipeline |
| DayNight | Filtering out COS data taken during the day or night |
| LSF | Working with the COS Line Spread Function (LSF) |
| Extract | Editing the extraction boxes in a BOXCAR-method spectral extraction file (XTRACTAB) |
5.1.1 COS-Specific Tasks
In Chapter 3, we gave detailed discussions of the use of the data reduction pipeline calcos. This task is contained in the PyPl package, calcos. Other modules useful for reducing and analyzing COS data are contained in a separate package called costools. A complete listing and brief description of these modules is given in Table 5.1. All of these tasks are run in Python. Consult the online help for each module for more information. Some of these modules will be discussed in greater detail in the remainder of this chapter.
Table 5.1: COS-Specific Modules.
Module | Description | Package |
|---|---|---|
calcos | Process COS data through the calibration pipeline | calcos |
saamodel | Get vertices for SAA model by model number | costools |
splittag | Split a | costools |
timefilter | Filter a corrtag table based on the TIMELINE extension | costools |
5.1.2 FITS Table Tasks
COS 1-D spectra and most COS reference files are stored in FITS tables (see Section 2.4.1 for a description of the structure of the table extension files for 1-D spectra and TIME-TAG data). Tasks designed to handle this format can be used to examine and extract information from the tables. Here we give specific examples of the use of routines in astropy.table.Table and astropy.io.fits to help you get started. A sample output is given after each command.
Use the Table module in the astropy.table package to find out what information is given in the columns of a FITS table.
> from astropy.table import Table
> t = Table.read('ldel05k2q_x1d.fits')
> t.info # to print the column info
<Table length=2>
name dtype shape unit format
-------------------- ------- ------- ---------------------- -------
SEGMENT bytes4
EXPTIME float64 s {:8.3f}
NELEM int32 {:6d}
WAVELENGTH float64 (1274,) Angstrom
FLUX float32 (1274,) erg / (Angstrom cm2 s)
ERROR float32 (1274,) erg / (Angstrom cm2 s)
ERROR_LOWER float32 (1274,) erg / (Angstrom cm2 s)
VARIANCE_FLAT float32 (1274,)
VARIANCE_COUNTS float32 (1274,)
VARIANCE_BKG float32 (1274,)
GROSS float32 (1274,) ct / s
GCOUNTS float32 (1274,) ct
NET float32 (1274,) ct / s
BACKGROUND float32 (1274,) ct / s
DQ int16 (1274,)
DQ_WGT float32 (1274,)
DQ_OUTER int16 (1274,)
BACKGROUND_PER_PIXEL float32 (1274,) ct / (pix s)
NUM_EXTRACT_ROWS int16 (1274,)
ACTUAL_EE float64 (1274,)
Y_LOWER_OUTER float64 (1274,)
Y_UPPER_OUTER float64 (1274,)
Y_LOWER_INNER float64 (1274,)
Y_UPPER_INNER float64 (1274,)
Use print() to look at the contents of the table:
> from astropy.table import Table
> t = Table.read('ldel05k2q_x1d.fits')
> print(t)
SEGMENT EXPTIME ... Y_UPPER_INNER [16384]
s ...
------- ------- ... ---------------------
FUVA 503.744 ... 421.0 .. 423.0
Note that the number of columns displayed is limited by the width of the window that you are working in. To see more columns, you can adjust your window to be wider and re-type the command:
> print(t)
SEGMENT EXPTIME NELEM ... Y_UPPER_OUTER [16384] Y_LOWER_INNER [16384] Y_UPPER_INNER [16384]
s ...
------- ------- ----- ... --------------------- --------------------- ---------------------
FUVA 503.744 16384 ... 427.0 .. 430.0 414.0 .. 410.0 421.0 .. 423.0
The table output indicates that some of the columns contain arrays of 16384 elements rather than a single value. For those columns, printing the table element t displays the value of the first and last elements in the array.
The user may wish to look at values in a COS reference file FITS table, which generally have many rows. Each row typically characterizes a specific operating mode, location on the detector, a value of a parameter to be used in the reduction, etc.
> from astropy.io import fits
> from astropy.table import Table
> hdulist = fits.open('s7g1700hl_disp.fits')
> tbdata = hdulist[1].data
> print(Table(tbdata))
SEGMENT OPT_ELEM APERTURE CENWAVE NELEM COEFF [4]
------- -------- -------- ------- ----- ---------------------
FUVA G130M PSA 1291 2 1272.64958352 .. -0.0
FUVA G130M WCA 1291 2 1272.21660012 .. -0.0
FUVA G130M BOA 1291 2 1272.64958352 .. -0.0
FUVB G130M PSA 1291 2 1119.38291962 .. -0.0
FUVB G130M WCA 1291 2 1118.98191865 .. -0.0
FUVB G130M BOA 1291 2 1119.38291962 .. -0.0
FUVA G130M PSA 1300 2 1282.84550637 .. -0.0
FUVA G130M WCA 1300 2 1282.41252297 .. -0.0
FUVA G130M BOA 1300 2 1282.84550637 .. -0.0
FUVB G130M PSA 1300 2 1129.58656362 .. -0.0
FUVB G130M WCA 1300 2 1129.18055776 .. -0.0
FUVB G130M BOA 1300 2 1129.58656362 .. -0.0
...
> hdulist.close()
In Python, there are often many methods to accomplish the same goal. For instance, in the above example, one could have used fits.getdata() instead of fits.open() (as shown in examples below). When using fits.open(), it is important to close the fits file afterwards as shown in the last line of this example.
5.1.3 General Spectral Display and Analysis Tasks
Table 5.2 lists several useful Python packages for displaying and analyzing COS spectral data.
Table 5.2: Spectral Analysis Tasks.
Task | Input Formats | Purpose |
|---|---|---|
specviz | 1-D fits or ASCII files | Generating 1-D plots, fitting curves to spectral features, and much more with a graphical user interface |
matplotlib.pyplot.plot() | 1-D arrays | Generating 1-D plots from columns in tables |
matplotlib.pyplot.imshow() | 2-D image arrays | Generating 2-D plots from 2-D image arrays |
astropy.modeling | 1-D and 2-D arrays | Fitting curves to spectral features |
scipy.optimize.curve_fit() | 1-D arrays | Fitting curves to spectral features |
Specvis is a powerful spectral plotting and fitting package that provides the user with many options for modeling emission and absorption features, continuum, and much more. It is controlled via a graphical user interface, and has documentation and examples available: https://specviz.readthedocs.io/. Please note that this is a new package which has not been fully tested and verified as of May 2018.
Plotting COS Spectral Images (flt)
The COS FUV data consist of separate files for each of the two FUV detector segments, segment A and segment B. In general, each segment contains one target spectrum and one wavecal. On the other hand, for the NUV data there are three disjointed portions of the spectrum present on the image for both the target and the wavecal. The following examples illustrate the use of the imshow() function in matplotlib.pyplot module to plot COS FUV and NUV spectra from the 2-D flt files. The parameters in plt.ylim() refer to the Y location of the spectrum for each particular grating/segment combination. One may need to adjust the scaling factors, vmin and vmax, to suit the data range of any particular dataset. A range of different colormaps may be used, and can be found here:
https://matplotlib.org/tutorials/colors/colormaps.html.
The parameters in plt.ylim() refer to the Y location of the spectrum for each particular grating/segment combination. For the FUV detector, the Y location also depends on lifetime position (LP). The example below shows both segments of a lifetime position 4, G160M spectrum.
# target spectrum, FUV segment A:
> import matplotlib.pyplot as plt
> from astropy.io import fits
> fltdata = fits.getdata('ldel05ivq_flt_a.fits')
> plt.imshow(fltdata, cmap=plt.get_cmap('gist_yarg'), \
> aspect='auto', vmin=0.0, vmax=0.2, origin='lower')
> plt.ylim([410., 430.])
> plt.xlabel('X FULL')
> plt.ylabel('Y FULL')
> plt.colorbar()
> plt.show()
# wavecal spectrum, FUV segment A:
> import matplotlib.pyplot as plt
> from astropy.io import fits
> fltdata = fits.getdata('ldel05ivq_flt_a.fits')
> plt.imshow(fltdata, cmap=plt.get_cmap('gist_yarg'), \
> aspect='auto', vmin=0.0, vmax=0.01, origin='lower')
> plt.ylim([510., 535.])
> plt.xlabel('X FULL')
> plt.ylabel('Y FULL')
> plt.colorbar()
> plt.show()
# target spectrum, NUV stripes A, B, C:
> import matplotlib.pyplot as plt
> from astropy.io import fits
> fltdata = fits.getdata('lcwj10qtq_flt.fits')
# stripe A
> plt.imshow(fltdata, cmap=plt.get_cmap('gist_yarg'), \
> aspect='auto', vmin=0.0, vmax=0.01, origin='lower')
> plt.xlabel('XFULL')
> plt.ylabel('YFULL')
> plt.ylim([175.0, 190.0])
> plt.colorbar()
> plt.show()
# stripe B
> plt.imshow(fltdata, cmap=plt.get_cmap('gist_yarg'), \
> aspect='auto', vmin=0.0, vmax=0.01, origin='lower')
> plt.xlabel('XFULL')
> plt.ylabel('YFULL')
> plt.ylim([275.0, 290.0])
> plt.colorbar()
> plt.show()
# stripe C
> plt.imshow(fltdata, cmap=plt.get_cmap('gist_yarg'), \
> aspect='auto', vmin=0.0, vmax=0.01, origin='lower')
> plt.xlabel('XFULL')
> plt.ylabel('YFULL')
> plt.ylim([420.0, 435.0])
> plt.colorbar()
> plt.show()
Plotting COS Tabular Spectra (x1d)
The following example illustrates the use of the Python package matplotlib.pyplot to plot COS tabular data. See Section 5.3.1 for more details and examples.
> import matplotlib.pyplot as plt
> from astropy.io import fits
> x1ddata = fits.getdata('ldel05k2q_x1d.fits')
> # select only segment A
> x1d_fuva = x1ddata[x1ddata['segment'] == 'FUVA']
> wavelength = x1d_fuva['wavelength'][0]
> flux = x1d_fuva['flux'][0]
> plt.plot(wavelength, flux)
> plt.xlabel('Wavelength')
> plt.ylabel('Flux')
> plt.show()
The above example selects the portion of the data that corresponds to the segment we want to plot. The segment selector can be either 'FUVA' or 'FUVB' for FUV data and 'NUVA', 'NUVB', or 'NUVC' for NUV data.
The next example shows how one would plot both FUV segments by using the concatenate() module in numpy.
> import matplotlib.pyplot as plt
> from astropy.io import fits
> import numpy as np
> x1ddata = fits.getdata('ldel02cxq_x1d.fits')
> wavelength = np.concatenate((x1ddata["wavelength"][1],\
x1ddata["wavelength"][0]))
> flux = np.concatenate((x1ddata["flux"][1], \
x1ddata["flux"][0]))
> plt.plot(wavelength, flux)
> plt.xlabel('Wavelength')
> plt.ylabel('Flux')
> plt.show()
We note that instead of using concatenate(), one could also use the flatten() module as shown in the example in Section 5.3.1.