Aligning HST images to an absolute reference catalog


Note: The notebook in this repository 'Initializtion.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'.
Note: This notebook is based on WFC3 ISR 2017-19: Aligning HST Images to Gaia: a Faster Mosaicking Workflow and contains a subset of the information/code found in the repository here. For more information, see the notebook in that repository titled 'Gaia_alignment.ipynb'.

Introduction

The alignment of HST exposures is a critical step in image stacking/combination performed by software such as AstroDrizzle. Generally, a relative alignment is performed that aligns one image (or multiple images) to another image which is designated as the reference image. This makes it so the images are aligned to each other, but the pointing error of the observatory can still cause the images to have incorrect absolute astrometry.

When absolute astrometry is desired, the images can be aligned to an external catalog that is known to be on an absolute frame. In this example, we will provide a workflow to query catalogs such as SDSS and Gaia via the astroquery package, and then align the images to that catalog via TweakReg.

For more information about TweakReg, see the other notebooks in this repository or the TweakReg Documentation.

For more information on Astroquery, see the other notebooks in this repository or the Astroquery Documentation.

In [1]:
import astropy.units as u
import glob
import numpy as np
import matplotlib.pyplot as plt
import os

from astropy.io import fits
from astropy.table import Table
from astropy.units import Quantity
from astropy.coordinates import SkyCoord
from astroquery.gaia import Gaia
from astroquery.mast import Observations
from astroquery.sdss import SDSS

from ccdproc import ImageFileCollection
from IPython.display import Image

from drizzlepac import tweakreg
from drizzlepac import astrodrizzle
Created TAP+ (v1.2.1) - Connection:
	Host: gea.esac.esa.int
	Use HTTPS: True
	Port: 443
	SSL Port: 443
Created TAP+ (v1.2.1) - Connection:
	Host: geadata.esac.esa.int
	Use HTTPS: True
	Port: 443
	SSL Port: 443
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-8f19f03d4e2c> in <module>
     16 from IPython.display import Image
     17 
---> 18 from drizzlepac import tweakreg
     19 from drizzlepac import astrodrizzle

/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

For this example, we will use WFC3/UVIS images of NGC 6791 from Visit 01 of proposal 12692.

In [ ]:
# Get the observation records
obsTable = Observations.query_criteria(obs_id='ibwb01*', proposal_id=12692, obstype='all', filters='F606W')

# Get the listing of data products
products = Observations.get_product_list(obsTable)

# Filter the products for exposures
filtered_products = Observations.filter_products(products, productSubGroupDescription='FLC')

# Show the table
filtered_products
In [ ]:
# Download all the images above
Observations.download_products(filtered_products, mrp_only=False)
In [ ]:
# For convenience, move the products into the current directory.
for flc in glob.glob('./mastDownload/HST/*/*flc.fits'):
    flc_name = os.path.split(flc)[-1]
    os.rename(flc, flc_name)

Inspect the image header

The cell below shows how to query information from the image header using ImageFileCollection in ccdproc. We see that the 1st exposure is 30 seconds and the 2nd and 3rd exposures are 360 seconds. The 3rd exposure is dithered by ~82" in the Y-direction which is approximately the width of one UVIS chip.

In [ ]:
collec = ImageFileCollection('./', glob_include="*flc.fits", ext=0,
                             keywords=["targname", "ra_targ", "dec_targ", "filter", "exptime", "postarg1", "postarg2"])

table = collec.summary
table['exptime'].format = '7.1f'
table['ra_targ'].format = '7.7f'
table['dec_targ'].format = '7.7f'
table['postarg1'].format = '7.2f'
table['postarg2'].format = '7.2f'
table

2. Querying catalogs

Now that we have the images, we will download the reference catalogs from both SDSS and Gaia using astroquery.

2a. Identify Coordinates

We will first create a SkyCoord Object to point astroquery to where we are looking on the sky. Since our example uses data from NGC 6791, we will use the ra_targ and dec_targ keywords from the first image to get the coordinates of the object.

In [ ]:
RA = table['ra_targ'][0]
Dec = table['dec_targ'][0]

2b. SDSS Query

We now give those values to an astropy SkyCoord object, which we will pass to the SDSS. Additionally, we use an astropy Quantity object to create a radius for the SDSS query. We set the radius to 6 arcminutes to comfortably cover the area of our images. For reference UVIS detector field of view is ~2.7'x2.7' and a y-dither of 82" covers a total area on the sky of ~2.7'x4.1'.

In [ ]:
coord = SkyCoord(ra=RA, dec=Dec, unit=(u.deg, u.deg))
radius = Quantity(6., u.arcmin)

Then we only need to perform the query via the SDSS.query_region method of astroquery.sdss. The spectro=False keyword argument means we want to exclude spectroscopic objects, as we are looking for objects to match with an image.

In the fields parameter, we specify a list of fields we want returned by the query. In this case we only need the position, and maybe a magnitude 'g' if we want to cut very dim and/or bright objects out of the catalog, as those are likely measured poorly. Details on selecting objects by magnitude may be found in the original 'Gaia_alignment' notebook. Many other fields are available in the SDSS query and are documented here.

In [ ]:
sdss_query = SDSS.query_region(coord, radius=radius, spectro=False, fields=['ra', 'dec', 'g'])
sdss_query

We then just need to save the catalog to use with TweakReg. Since the returned value of the query is an astropy.Table, saving the file is very straightforward:

In [ ]:
sdss_query.write('sdss.cat', format='ascii.commented_header')

2c. Gaia Query

Similarly to SDSS, we can query Gaia catalogs for our target via astroquery.gaia. We can use the same coord and radius from the SDSS query.

In [ ]:
gaia_query = Gaia.query_object_async(coordinate=coord, radius=radius)
gaia_query

This query has returned very large number of columns. We want to pare down the catalog to make it easier to use with TweakReg.
We can select only the useful columns via:

In [ ]:
reduced_query = gaia_query['ra', 'dec', 'phot_g_mean_mag']
reduced_query

Then we write this catalog to an ascii file for use with TweakReg.

In [ ]:
reduced_query.write('gaia.cat', format='ascii.commented_header')

3. Running TweakReg

With the catalogs downloaded and the headers populated, we simply need to run TweakReg with each catalog passed into the refcat parameter. The steps below compare the astrometric residuals obtained from aligning to each refcat. In each test case, we set updatehdr to False until we are satisfied with the alignment by inspecting both the shift file and the astrometric residual plots.

3a. SDSS Alignment

In [ ]:
refcat = 'sdss.cat'
cw = 3.5  # Set to two times the FWHM of the PSF of the UVIS detector
wcsname = 'SDSS'  # Specify the WCS name for this alignment

tweakreg.TweakReg('*flc.fits',  # Pass input images
                  updatehdr=False,  # update header with new WCS solution
                  imagefindcfg={'threshold': 500., 'conv_width': cw},  # Detection parameters, threshold varies for different data
                  refcat=refcat,  # Use user supplied catalog (Gaia)
                  interactive=False,
                  see2dplot=False,
                  shiftfile=True,  # Save out shift file (so we can look at shifts later)
                  outshifts='SDSS_shifts.txt',  # name of the shift file
                  wcsname=wcsname,  # Give our WCS a new name
                  reusename=True,
                  ylimit=1.5,
                  fitgeometry='general')  # Use the 6 parameter fit

We can look at the shift file to see how well the fit did (or we could open the output png images for more information).

The columns are:

  • Filename
  • X Shift [pixels]
  • Y Shift [pixels]
  • Rotation [degrees]
  • Scale
  • X RMS [pixels]
  • Y RMS [pixels]
In [ ]:
for line in open('SDSS_shifts.txt').readlines():
    print(line)
In [ ]:
# Astrometric residual plots
Image(filename='residuals_ibwb01xqq_flc.png',width=500, height=300)
In [ ]:
Image(filename='residuals_ibwb01xrq_flc.png',width=500, height=300)
In [ ]:
Image(filename='residuals_ibwb01xxq_flc.png',width=500, height=300)

As we can see, the RMS is fairly large at about 0.5 pixels, which is not a great fit. This is likely because the SDSS astrometric precision is not high enough to get good HST alignment. One approach would be to align the first image to SDSS and then align the remaining HST images to one another. This would improve both the absolute and relative alignment of the individual frames.

3b. Gaia Alignment

In [ ]:
refcat = 'gaia.cat'
cw = 3.5  # Set to two times the FWHM of the PSF.
wcsname = 'Gaia'  # Specify the WCS name for this alignment

tweakreg.TweakReg('*flc.fits',  # Pass input images
                  updatehdr=False,  # update header with new WCS solution
                  imagefindcfg={'threshold':500.,'conv_width':cw},  # Detection parameters, threshold varies for different data
                  refcat=refcat,  # Use user supplied catalog (Gaia)
                  interactive=False,
                  see2dplot=False,
                  shiftfile=True,  # Save out shift file (so we can look at shifts later)
                  outshifts='Gaia_shifts.txt',  # name of the shift file
                  wcsname=wcsname,  # Give our WCS a new name
                  reusename=True,
                  sigma=2.3,
                  ylimit=0.2,
                  fitgeometry='general')  # Use the 6 parameter fit

We can similarly look at the shift file from alignment to the Gaia catalog:

In [ ]:
for line in open('Gaia_shifts.txt').readlines():
    print(line)
In [ ]:
# Astrometric residual plots
Image(filename='residuals_ibwb01xqq_flc.png',width=500, height=300)
In [ ]:
Image(filename='residuals_ibwb01xrq_flc.png',width=500, height=300)
In [ ]:
Image(filename='residuals_ibwb01xxq_flc.png',width=500, height=300)

As expected, the Gaia catalog does quite a bit better, with rms residuals less tha 0.05 pixels.

To apply these transformations to the image, we simply need to run TweakReg the same as before, but set the parameter updatehdr equal to True:

In [ ]:
refcat = 'gaia.cat'
cw = 3.5  # Set to two times the FWHM of the PSF.
wcsname = 'Gaia'  # Specify the WCS name for this alignment

tweakreg.TweakReg('*flc.fits',  # Pass input images
                  updatehdr=True,  # update header with new WCS solution
                  imagefindcfg={'threshold': 500., 'conv_width': cw},  # Detection parameters, threshold varies for different data
                  refcat=refcat,  # Use user supplied catalog (Gaia)
                  interactive=False,
                  see2dplot=False,
                  shiftfile=True,  # Save out shift file (so we can look at shifts later)
                  outshifts='Gaia_shifts.txt',  # name of the shift file
                  wcsname=wcsname,  # Give our WCS a new name
                  reusename=True,
                  sigma=2.3,
                  fitgeometry='general')  # Use the 6 parameter fit

4. Drizzle the Data

While the three sets of FLC files are now aligned, we drizzle together only the two long exposures.

When exposures are very different lengths, drizzling them together doesn't work well when 'EXP' weighting is used. For objects that saturate in the long exposures, the problem occurs at the boundary where you transition from only short exposure to short plus long. Here the pixels getting power from long exposure pixels are only getting power from pixels whose centers are outside the ring, and thus they are weighted lower than they would be if they were getting values from both inside and outside the ring. The result is a discontinuity in the PSF radial profile and a resulting flux which is too low in those boundary pixels. For photometry of saturated objects, the short exposures should be drizzled separately from the long exposures.

In [ ]:
astrodrizzle.AstroDrizzle('ibwb01x[rx]q_flc.fits', 
                          output='f606w',
                          preserve=False,
                          clean=True, 
                          build=False,
                          context=False,
                          skymethod='match',
                          driz_sep_bits='64, 32',
                          combine_type='minmed',
                          final_bits='64, 32')
In [ ]:
# Display the combined science and weight images 

sci = fits.getdata('f606w_drc_sci.fits')
wht = fits.getdata('f606w_drc_wht.fits')

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

ax1.imshow(sci, vmin=-0.05, vmax=0.4, cmap='Greys_r', origin='lower')
ax2.imshow(wht, vmin=0, vmax=1000, cmap='Greys_r', origin='lower')

Conclusions

Many other services have interfaces for querying catalogs which could also be used to align HST images. In general, Gaia works very well for HST due to it's high precision, but can have a low number of sources in some regions, especially at high galactic latitudes. Aligning images to an absolute frame provides an easy way to make data comparable across many epochs/detectors/observatories, and in many cases, makes the alignment much easier.

About this Notebook

Author: V. Bajaj, STScI WFC3 Team
Updated: December 14, 2018