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 AstroConda python distribution, which is available for download from STScI:
http://astroconda.readthedocs.io/en/latest/

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.

5.1.1 COS-Specific AstroConda Tasks

In Chapter 3, we gave detailed discussions of the use of the data reduction pipeline calcos. This task is contained in the AstroConda 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

AstroConda Package

calcos

Process COS data through the calibration pipeline

calcos

saamodel

Get vertices for SAA model by model number

costools

splittag

Split a corrtag file into multiple sub-files by time interval

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=1>
        name          dtype   shape            unit         
-------------------- ------- ------- ----------------
             SEGMENT    str4                                
             EXPTIME float64                               s
               NELEM   int32                                
          WAVELENGTH float64 (16384,)               Angstrom
                FLUX float32 (16384,) erg / (Angstrom cm2 s)
               ERROR float32 (16384,) erg / (Angstrom cm2 s)
               GROSS float32 (16384,)                 ct / s
             GCOUNTS float32 (16384,)                     ct
                 NET float32 (16384,)                 ct / s
          BACKGROUND float32 (16384,)                 ct / s
                  DQ   int16 (16384,)                       
              DQ_WGT float32 (16384,)                       
            DQ_OUTER   int16 (16384,)                       
BACKGROUND_PER_PIXEL float32 (16384,)           ct / (pix s)
    NUM_EXTRACT_ROWS   int16 (16384,)                       
           ACTUAL_EE float64 (16384,)                       
       Y_LOWER_OUTER float64 (16384,)                       
       Y_UPPER_OUTER float64 (16384,)                       
       Y_LOWER_INNER float64 (16384,)                       
       Y_UPPER_INNER float64 (16384,)  
                      

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.

# 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.