Timefilter
This tutorial outlines the steps needed to filter out unwanted COS data in an exposure. To do this, we will use the TimeFilter python module to select bad data in COS corrtag files, and then re-extract the data into x1d spectra.
This notebook is written in Python, as that is the language of the calibration pipeline (CalCOS) and the Timefilter module. This task can also be run from PyRAF, and is located in the stsdas.hst_calib.hstcos package.
#-- Necessary library imports for the tutorial
import os
import glob
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from astropy.io import fits
from calcos import calcos
from costools import timefilter
#-- iPython notebook 'magic' command to make plots appear within the
#-- browser instead of opening a separate plotting window.
%matplotlib inline
#-- Increase the default fontsize for the plots
mpl.rcParams['font.size'] = 16
Understanding Airglow
HST UV spectra are subject to contamination from geocoronal airglow emission. The most prominent example of this is the strong emission from Hydrogen Lyman Alpha at 1215.6 Anstroms; present in all COS observations that cover that wavelength.
A list of COS airglow observations, along with possible UV airglow lines, can be found on the airglow webpage
Airglow in data
Three of the strongest lines seen in COS data are given below. The strengths of these lines varies from observation to observation, due to changes in the solar activity and changes in HST's orbital position.
- Lyman Alpha 1215.67
- OI 1302.2 - 1306.0
- NI 1199.5 - 1200.7
To see this, we will look at some COS observations taken of blank sky, and thus contain only airglow features. The top panel in the figure below shows the region around the OI emission, while the bottom panel shows both the Lyman Alpha and the faint NI line. Two observations have been overplotted, and the variation in airglow strength can been seen.
fig = plt.figure(figsize=(15, 15))
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
for item in glob.glob('data/lbh*x1d.fits'):
hdu = fits.open(item)
ax1.plot(hdu[1].data['WAVELENGTH'].ravel(), hdu[1].data['FLUX'].ravel(), label=os.path.split(item)[1][:9])
ax2.plot(hdu[1].data['WAVELENGTH'].ravel(), hdu[1].data['FLUX'].ravel(), label=os.path.split(item)[1][:9])
for axis in [ax1, ax2]:
axis.grid(True)
axis.set_xlabel("Wavelength")
axis.legend(shadow=True, numpoints=1, loc='upper left')
ax1.set_ylim(0, 5e-14)
ax1.set_xlim(1295, 1315)
ax2.set_ylim(0, 2e-14)
ax2.set_xlim(1190, 1230)
ax1.set_ylabel('FUVA')
ax2.set_ylabel('FUVB')
ax1.set_title('Airglow Emission')
Airglow Variation
Airglow line strengths are driven by the solar activity influencing the Earth's atmosphere. This impact is severely lessened when the sun is below the horizon, and thus airglow contamination be diminished or eliminated by utilizing only data taken during orbital night.
To facilitate filtering out specific data, a TIMELINE extension was added to corrag data files to track various parameters throughout an exposure. The different parameters are shown below, and each one is calculated on a second-by-second basis throughout the exposure.
hdu = fits.open('data/lbh408vaq_corrtag_b.fits')
print hdu['timeline'].data.names
Now let's plot the strength of two airglow lines against the Sun's altitude during the exposure. You can see that, although Lyman Alpha only diminishes during orbital night, OI fades away. Multiple datasets are over-plotted, and all show the same trend.
fig = plt.figure(figsize=(15, 15))
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
for item in glob.glob('data/lbh*corrtag_b.fits'):
hdu = fits.open(item)
ax1.plot(hdu['timeline'].data['SUN_ALT'], hdu['timeline'].data['LY_ALPHA'], 'o', label=os.path.split(item)[1][:9])
for item in glob.glob('data/lbh*corrtag_a.fits'):
hdu = fits.open(item)
ax2.plot(hdu['timeline'].data['SUN_ALT'], hdu['timeline'].data['OI_1304'],'o', label=os.path.split(item)[1][:9])
for axis in [ax1, ax2]:
axis.grid(True)
axis.axvspan(xmin=-45, xmax=0, ymin=0, ymax=1, label='Orbital Night', alpha=.5, color='grey')
axis.set_xlim(-45, 90)
axis.set_xlabel("Sun altitude above horizon (degrees)", fontsize=16)
axis.legend(shadow=True, numpoints=1, loc='upper left', fontsize=16)
ax1.set_ylabel('Lyman Alpha (cnts/s)', fontsize=16)
ax2.set_ylabel('Oxygen I (cnts/s)', fontsize=16)
ax1.set_title('Airglow Lines with Altitude')
fig.savefig('airglow_strength.pdf', bbox_inches='tight')
Dealing with airglow in data
We can use this TIMELINE extension in coordination with the Timefilter module to filter out selected times of input data. Any of the parameters listed in the TIMELINE extension can be used, depending on your interest, but for our current purpose we will focus on the task of eliminating as much airglow as possible.
First we will plot some real data contaminated by airglow. The data shown below is taken by PI Andrew Fox as part of Program 12604. We will be looking at the G130M/1291 data to view the OI airglow line. In the plotted spectra, there are very prominent and variable emission features in the region of OI superimposed on the source continuum.
for dataset in glob.glob('data/lbry*x1d.fits'):
hdu = fits.open(dataset)
print dataset, hdu[0].header['OPT_ELEM'], hdu[0].header['CENWAVE'], hdu[1].header['EXPTIME']
fig = plt.figure(figsize=(17, 8))
ax = fig.add_subplot(1, 1, 1)
for dataset in glob.glob('data/lbry*x1d.fits'):
hdu = fits.open(dataset)
ax.plot(hdu[1].data['wavelength'].ravel(), hdu[1].data['flux'].ravel(), label=hdu[0].header['rootname'])
#-- Label with arrows
ax.annotate('Geo-coronal Oxygen I',
xycoords='data',
xy=(1302, 1e-13),
xytext=(1310, 1.5e-13),
arrowprops=dict(facecolor='black', shrink=0.05),
fontsize=16)
ax.annotate('',
xycoords='data',
xy=(1305, 1e-13),
xytext=(1310, 1.5e-13),
arrowprops=dict(facecolor='black', shrink=0.05),
fontsize=16)
ax.set_ylim(0, 2e-13)
ax.set_xlim(1280, 1340)
ax.legend(shadow=True, numpoints=1, fontsize=16)
ax.set_ylabel('Flux')
ax.set_xlabel('Wavelength $(\AA)$')
ax.set_title('Spectrum Contaminated by Airglow')
fig.savefig('spectrum_before.pdf', bbox_inches='tight')
Timefilter
#-- timefilter and CalCOS need reference files to function.
#-- this will set the lref environment variable if you don't have it set.
#--
#-- This is the correct location for users internal to STScI, you may have to change it.
if not 'lref' in os.environ:
os.environ['lref'] = '/grp/hst/cdbs/lref/'
for dataset in glob.glob('data/lbry*corrtag*.fits'):
filepath, filename = os.path.split(dataset)
print "Filtering ", filename
timefilter.TimelineFilter(input=dataset, filter='SUN_ALT > 0')
The corrtags have had the unwanted data flagged as bad, but we still need to re-extract the spectrum. To do this, simply call the CalCOS pipeline on the corrtag files. The pipeline will re-run the needed steps and extract the spectrum on the filtered data.
Note:
- CalCOS will automatically search for the FUVB segment if only FUVA data is supplied (and vice-versa)
- When running CalCOS on corrtag data (instead of rawtags) a different output directory MUST be supplied
#-- This step will take a little while as we re-un CalCOS on the filtered datasets.
for dataset in glob.glob('data/lbry*corrtag_a*'):
calcos(dataset, outdir='data/out/')
fig = plt.figure(figsize=(17, 8))
ax = fig.add_subplot(1, 1, 1)
for dataset in glob.glob('data/out/lbry*x1d.fits'):
hdu = fits.open(dataset)
ax.plot(hdu[1].data['wavelength'].ravel(), hdu[1].data['flux'].ravel(), label=hdu[0].header['rootname'])
ax.set_ylim(0, 2e-13)
ax.set_xlim(1280, 1340)
ax.legend(shadow=True, numpoints=1, fontsize=16)
ax.set_title('Airglow Filtered Spectrum')
ax.set_ylabel('Flux')
ax.set_xlabel('Wavelength $(\AA)$')
fig.savefig('spectrum_after.pdf', bbox_inches='tight')
You may have noticed that only a single spectrum was visible here, though two were seen before filtering. This is not an error, but rather that only 1 of the two datasets actually contains data taken during orbital night. For a dataset containing only data with SUN_ALT > 0, no events would still be remaining after filtering, and the spectra are empty. Even for the one that remains, the signal-to-noise is considerably reduced. That is because only a small portion of the data was taken during night.
Again, the TIMELINE extension can help visualize this: below is plotted the SUN_ALT of each observation vs time. Only lbry01i6q contains data taken during orbital night.
fig = plt.figure(figsize=(13, 6))
ax = fig.add_subplot(1, 1, 1)
for dataset in glob.glob('data/out/lbry*corrtag_a.fits'):
hdu = fits.open(dataset)
file_path, file_name = os.path.split(dataset)
#-- Plot the sun altitude against time
ax.plot(hdu['timeline'].data['time'], hdu['timeline'].data['sun_alt'], lw=3)
#-- Label dataset name at the end of track
ax.annotate(file_name[:9],
xycoords='data',
xy=(hdu['timeline'].data['time'][-1], hdu['timeline'].data['sun_alt'][-1]),
xytext=(hdu['timeline'].data['time'][-1], hdu['timeline'].data['sun_alt'][-1]),
fontsize=16)
ax.axhspan(ymin=-20, ymax=0, xmin=0, xmax=1, color='grey', alpha=.5, label='Orbital Night')
ax.set_ylim(-20, 120)
ax.set_title('Sun Altitude During Exposure')
ax.set_xlabel('Time into observation (seconds)')
ax.set_ylabel('Sun Altitude (degrees)')
ax.legend(loc='upper left', shadow=True, numpoints=1)
