Aligning Deep Exposures of Sparse Fields

Note: The notebook in this repository 'Initialization.ipynb' goes over many of the basic concepts such as the setup of the environment/package installation and should be read first if you are new to HST images, DrizzlePac, or Astroquery.

Introduction

This notebook demonstrates aligning long exposures which have relatively few stars and a large number of cosmic rays. It is based on the example described in the ISR linked here (ACS ISR 2015-04: Basic Use of SExtractor Catalogs With TweakReg - I), but uses a much simpler methodology.

Rather than making use of external software (e.g. SExtractor) and going through the extra steps to create 'cosmic-ray cleaned' images for each visit, this notebook demonstrates new features in TweakReg designed to mitigate false detections.

TweakReg’s source finding task imagefind now includes parameters to exclude false detections and allow the software to more easily solve for the image offsets using matched sources lists. For example, dqbits is a list of DQ flag values to include as 'good' or to exclude as 'bad' before generating and matching source lists. For ACS/WFC, setting dqbits=-16 will mask hot pixels flagged by the instrument team, eliminating a common problem where TweakReg locks onto hot pixels and solves for the dither pattern. This can occur when users set the detection threshold value too low and hot pixels outnumber astronomical sources. Other new parameters allow selection for sharpness and roundness, which give users better control over source selection criteria and the mitigation of artifacts. More details on imagefindpars options may be found on the following webpage.

In [1]:
from astropy.table import Table
from astropy.io import fits
from astroquery.mast import Observations
from ccdproc import ImageFileCollection
import glob
import matplotlib.pyplot as plt
import os
import shutil
from drizzlepac import tweakreg
from drizzlepac import astrodrizzle
from IPython.display import Image
The following task in the stsci.skypac package can be run with TEAL:
                                    skymatch                                    

ModuleNotFoundErrorTraceback (most recent call last)
ModuleNotFoundError: No module named 'numpy.core._multiarray_umath'
 Coordinate transformation and image resampling library NOT found!

 Please check the installation of this package to insure C code was built successfully.

ImportErrorTraceback (most recent call last)
/opt/conda/envs/notebooks_env/lib/python3.6/site-packages/drizzlepac/ablot.py in <module>
     22 try:
---> 23     from . import cdriz
     24 except ImportError:

ImportError: numpy.core.multiarray failed to import

During handling of the above exception, another exception occurred:

ImportErrorTraceback (most recent call last)
<ipython-input-1-b02bebe9f452> in <module>
      7 import os
      8 import shutil
----> 9 from drizzlepac import tweakreg
     10 from drizzlepac import astrodrizzle
     11 from IPython.display import Image

/opt/conda/envs/notebooks_env/lib/python3.6/site-packages/drizzlepac/__init__.py in <module>
     19 from .version import *
     20 
---> 21 from . import ablot
     22 from . import adrizzle
     23 from . import astrodrizzle

/opt/conda/envs/notebooks_env/lib/python3.6/site-packages/drizzlepac/ablot.py in <module>
     26     print('\n Coordinate transformation and image resampling library NOT found!')
     27     print('\n Please check the installation of this package to insure C code was built successfully.')
---> 28     raise ImportError
     29 
     30 from .version import *

ImportError: 

1. Download the Data

The data in this example are comprised of 4 exposures in the F814W filter, all from Visit 0X of HST program 10092. Each exposure was dithered by ~60 pixels along the y-axis in order to obtain coverage in the area of the CCD chip gap. The X and Y dithers are given in arcseconds by the POSTARG1 and POSTARG2 keywords recorded in the image header.

The following commands query MAST, download the calibrated, CTE-corrected FLC files, and place them in the same 'working' directory as this notebook.

In [ ]:
# Query MAST for the F814W files.
science_list = Observations.query_criteria(proposal_id='10092', filters='F814W', obs_id='j8xi0x*')
Observations.download_products(science_list['obsid'], mrp_only=False, download_dir='./science',
                               productSubGroupDescription=['FLC'])

science_files = glob.glob('science/mastDownload/HST/*/*fits')
for im in science_files:
    root = im.split('/')[-1]
    os.rename(im, './'+root)
shutil.rmtree('science/')

2. Inspect the image headers

In [ ]:
collect = ImageFileCollection('./',
                              keywords=["asn_id", "detector", "filter2", "exptime", "postarg1", "postarg2"],
                              glob_include="*flc.fits", ext=0)
out_table = collect.summary
out_table

3. TweakReg Alignment

Use TweakReg to align the FLC frames based on sources in the image. The provided input list (input_flc.list) is used to align the frames in the specified order with j8xi0xsaq_flc.fits as the reference image. The parameter conv_width specifies the convolution kernel width in pixels, with recommended values ~2x the PSF FWHM for detecting point sources in the FLC frame. For ACS/WFC & WFC3/UVIS, this parameter is typically set to 3.5 pixels and for WFC3/IR to 2.5 pixels, but the value can be increased in order to use compact objects such as small galaxies for alignment.

3a. Use 'default' parameters (Test1)

In [ ]:
tweakreg.TweakReg('@input_flc.list',
                  imagefindcfg={'threshold': 100,'conv_width': 3.5},
                  shiftfile=True, outshifts='shift814_flc_test1.txt',
                  updatehdr=False, interactive=False, ylimit=0.4)
In [ ]:
# Give the 'fit residual plots' a unique name for comparison with other tests.
residual_pngs = glob.glob("residual*png")
for png in residual_pngs: 
    path = os.path.abspath(os.path.join(os.curdir, png))
    new_path = os.path.abspath(os.path.join(os.curdir, 'test1_{}'.format(png)))
    os.rename(path, new_path)

To verify that TweakReg has found a good fit, it is important to inspect the fit residuals. Below are the dx,dy residuals for each FLC file with respect to the reference image j8xi0xsaq_flc.fits. TweakReg finds >100 matches per frame, but the RMS of the fit residuals is quite large, ~0.7 pixels, and the points are not nicely clustered around dx,dy=0.0, as expected for a good fit.

In [ ]:
Image(filename='test1_residuals_j8xi0xs0q_flc.png')
In [ ]:
Image(filename='test1_residuals_j8xi0xs3q_flc.png')
In [ ]:
Image(filename='test1_residuals_j8xi0xs6q_flc.png')
In [ ]:
# Inspect the shift file for Test1
shift_table = Table.read('shift814_flc_test1.txt', format='ascii.no_header',
               names=['file', 'dx', 'dy', 'rot', 'scale', 'xrms', 'yrms'])
formats = ['.2f', '.2f', '.3f', '.5f', '.2f', '.2f']
for i, col in enumerate(shift_table.colnames[1: ]):
    shift_table[col].format = formats[i]
shift_table

3b. Adjust conv_width to find extended objects (Test 2)

In order for TweakReg to use compact galaxies rather than point sources for alignment, we increase the convolution kernel width parameter conv_width from 3.5 to 6.0 pixels in order to find sources with a FWHM ~3 pixels in the FLC frames.

In [ ]:
tweakreg.TweakReg('@input_flc.list',
                  imagefindcfg={'threshold': 100, 'conv_width': 6.0},
                  shiftfile=True, outshifts='shift814_flc_test2.txt',
                  updatehdr=False, interactive=False, ylimit=0.4)
In [ ]:
# Give the 'fit residual plots' a unique name for comparison with subsequent tests.
residual_pngs = glob.glob("residual*png")
for png in residual_pngs: 
    path = os.path.abspath(os.path.join(os.curdir, png))
    new_path = os.path.abspath(os.path.join(os.curdir, 'test2_{}'.format(png)))
    os.rename(path, new_path)

TweakReg now finds matches ~50 objects per frame, and the fit for the first matched image j8xi0xs0q_flc.fits looks good, with an RMS ~0.1 pixels and with the residuals dx,dy clustered around 0.0. For the other two frames j8xi0xs3q_flc.fits and j8xi0xs6q_flc.fits, the RMS is ~0.2 pixels, and the points are not centered around dx,dy=0 pixels.

In [ ]:
Image(filename='test2_residuals_j8xi0xs0q_flc.png')
In [ ]:
Image(filename='test2_residuals_j8xi0xs3q_flc.png')
In [ ]:
Image(filename='test2_residuals_j8xi0xs6q_flc.png')
In [ ]:
# Inspect the shift file for Test2
shift_table = Table.read('shift814_flc_test2.txt', format='ascii.no_header', 
                       names=['file', 'dx', 'dy', 'rot', 'scale', 'xrms', 'yrms'])
formats = ['.2f', '.2f', '.3f', '.5f', '.2f', '.2f']
for i, col in enumerate(shift_table.colnames[1:]):
    shift_table[col].format = formats[i]
shift_table

3c. Exclude flagged pixels with dqbits (Test 3)

To further improve the alignment, we make use of flags in the DQ array of the FLC files. The source finding parameters in TweakReg may be modified to exclude flagged pixels when generating lists of objects in each image.

Setting the parameter dqbits=0 will consider all non-zero pixels in the DQ mask to be “bad” pixels, and the corresponding image pixels will not be used for source finding. The default value of 'None' will turn off the use of image’s DQ array for source finding. In this case, AstroDrizzle was already run by MAST on visit 0X, and cosmic-ray flags were populated in the DQ array of the FLC frames. Since the exposures within this visit were already well aligned, the 4096 flags for cosmic rays are useful for excluding false detections.

In [ ]:
tweakreg.TweakReg('@input_flc.list',
                  imagefindcfg={'threshold': 100, 'conv_width': 6.0, 'dqbits': 0},
                  shiftfile=True, outshifts='shift814_flc_test3.txt',
                  updatehdr=False, interactive=False, ylimit=0.4)
In [ ]:
# Give the 'fit residual plots' a unique name for comparison with other tests.
residual_pngs = glob.glob("residual*png")
for png in residual_pngs: 
    path = os.path.abspath(os.path.join(os.curdir, png))
    new_path = os.path.abspath(os.path.join(os.curdir, 'test3_{}'.format(png)))
    os.rename(path, new_path)   

In this third test, TweakReg again finds ~40 matches per frame, but with spurious detections eliminated has an easier time locking onto the correct solution. The fit residuals below all have an RMS ~0.1 pixels and the points are all clustered around dx,dy=0.

In [ ]:
Image(filename='test3_residuals_j8xi0xs0q_flc.png')
In [ ]:
Image(filename='test3_residuals_j8xi0xs3q_flc.png')
In [ ]:
Image(filename='test3_residuals_j8xi0xs6q_flc.png')
In [ ]:
# Inspect the shift file for Test2
shift_table = Table.read('shift814_flc_test3.txt', format='ascii.no_header',
                         names=['file', 'dx', 'dy', 'rot', 'scale', 'xrms', 'yrms'])
formats = ['.2f', '.2f', '.3f', '.5f', '.2f', '.2f']
for i, col in enumerate(shift_table.colnames[1:]):
    shift_table[col].format = formats[i]
shift_table

3d. Rerun TweakReg to update the header WCS

Now run TweakReg with updatehdr=True to update the image headers with this solution.

In [ ]:
tweakreg.TweakReg('@input_flc.list',
                  imagefindcfg={'threshold': 100, 'conv_width': 6.0, 'dqbits': 0},
                  shiftfile=False, updatehdr=True, interactive=False)

4. Drizzle the aligned frames

Combine the aligned FLC files with AstroDrizzle. The ACS team now corrects for stable hot pixels (DQ flag=16) via the dark reference files, so these pixels can be considered 'good'. Full well saturated pixels (DQ flag=256) and warm pixels (DQ flag=64) may also be treated as good. More details on the recommended drizzle parameters for ACS may be found in ISR 2017-02.

In [ ]:
astrodrizzle.AstroDrizzle('*flc.fits',
                          output='f814w',
                          preserve=False,
                          clean=True,
                          build=False,
                          context=False,
                          driz_sep_bits='256,64,16',
                          combine_type='median',
                          final_bits='256,64,16',
                          runfile='f814w_driz.log')

5. Display the original FLC image and the combined DRC data

In [ ]:
sci2 = fits.getdata('j8xi0xsaq_flc.fits', ('sci', 2))
sci1 = fits.getdata('j8xi0xsaq_flc.fits', ('sci', 1))

fig = plt.figure(figsize=(50, 50))
ax1 = fig.add_subplot(2, 2, 1)
ax2 = fig.add_subplot(1, 2, 1)

ax1.imshow(sci2, vmin=50, vmax=100, cmap='Greys_r', origin='lower')
ax2.imshow(sci1, vmin=50, vmax=100, cmap='Greys_r', origin='lower')
In [ ]:
sci = fits.getdata('f814w_drc_sci.fits')
fig = plt.figure(figsize=(20, 20))
plt.imshow(sci, vmin=-0.01, vmax=0.1, cmap='Greys_r', origin='lower')

Summary

TweakReg may be used to align HST images based on the position of objects in the frame. Point sources are typically used for this purpose, but compact sources such as background galaxies may also be used by increasing the value of the parameter conv_width in imagefindpars. The data quality arrays of the input calibrated frames may also be used to further improve the fits by telling TweakReg to ignore pixels with specific flags inthe DQ array via the parameter dqbits in imagefindpars.

About this Notebook

Author: J. Mack, STScI WFC3 Team
Updated: January 8, 2019