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.
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
For this example, we will use WFC3/UVIS images of NGC 6791 from Visit 01 of proposal 12692.
# 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
# Download all the images above
Observations.download_products(filtered_products, mrp_only=False)
# 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)
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.
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
Now that we have the images, we will download the reference catalogs from both SDSS and Gaia using astroquery.
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.
RA = table['ra_targ'][0]
Dec = table['dec_targ'][0]
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'.
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.
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:
sdss_query.write('sdss.cat', format='ascii.commented_header')
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.
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:
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.
reduced_query.write('gaia.cat', format='ascii.commented_header')
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.
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:
for line in open('SDSS_shifts.txt').readlines():
print(line)
# Astrometric residual plots
Image(filename='residuals_ibwb01xqq_flc.png',width=500, height=300)
Image(filename='residuals_ibwb01xrq_flc.png',width=500, height=300)
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.
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:
for line in open('Gaia_shifts.txt').readlines():
print(line)
# Astrometric residual plots
Image(filename='residuals_ibwb01xqq_flc.png',width=500, height=300)
Image(filename='residuals_ibwb01xrq_flc.png',width=500, height=300)
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:
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
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.
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')
# 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')
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.
Author: V. Bajaj, STScI WFC3 Team
Updated: December 14, 2018