Aligning HST images to an Absolute Reference Catalog#
Table of Contents#
Introduction
Import Packages
1. Download the Observations with astroquery
1.1 Check image header data
1.2 Inspect the alignment
2. Query the Reference Catalogs
2.1 Identify Coordinates
2.2 SDSS Query
2.3 Gaia Query
3. Align images with TweakReg
3.1 SDSS Alignment
3.2 Inspect the shift file and fit residuals
3.3 Gaia Alignment
3.4 Inspect the shift file and fit residuals
3.5 Overplot Matched Sources on the image
3.6 Update header with optimal parameters
4. Combine the images using AstroDrizzle
Conclusions
About this Notebook
Introduction#
The alignment of HST exposures is a critical step in image stacking or combination performed with software such as AstroDrizzle. Generally, a relative alignment is performed to align one or multiple images to another image that is designated as the reference image. The reference image is generally the deepest exposure and/or that covering the largest area of all the exposures. This process aligns the images 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 with an absolute world coordinate system (WCS). In this example, we will provide a workflow to query catalogs such as SDSS and Gaia using the astroquery package, and then align the images to that catalog using TweakReg.
The workflow in this notebook for aligning images to Gaia is based on WFC3 ISR 2017-19: Aligning HST Images to Gaia: a Faster Mosaicking Workflow and contains a subset of the information and code found in this repository. For more information, see the notebook in that repository titled Gaia_alignment.ipynb.
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 Packages#
Table of Contents
The following Python packages are required to run the Jupyter Notebook:
os - change and make directories
glob - gather lists of filenames
shutil - remove directories and files
numpy - math and array functions
matplotlib - make figures and graphics
astropy - file handling, tables, units, WCS, statistics
astroquery - download data and query databases
drizzlepac - align and combine HST images
import os
import glob
import shutil
import numpy as np
import matplotlib.pyplot as plt
from collections import defaultdict
from IPython.display import Image, clear_output
import astropy.units as u
from astropy.io import ascii, fits
from astropy.table import Table
from astropy.units import Quantity
from astropy.coordinates import SkyCoord
from astropy.visualization import ZScaleInterval
from astroquery.mast import Observations
from astroquery.sdss import SDSS
from astroquery.gaia import Gaia
from drizzlepac import tweakreg, astrodrizzle
from drizzlepac.processInput import getMdriztabPars
%config InlineBackend.figure_format = 'retina' # Improves the resolution of figures rendered in notebooks.
Gaia.MAIN_GAIA_TABLE = 'gaiadr3.gaia_source' # Change this to the desired Gaia data release.
Gaia.ROW_LIMIT = 100000 # Show all of the matched sources in the printed tables.
# Set the locations of reference files and retrieve the MDRIZTAB recommended drizzle parameters.
os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'
os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'
os.environ['CRDS_PATH'] = './crds_cache'
os.environ['iref'] = './crds_cache/references/hst/wfc3/'
1. Download the Observations with astroquery#
Table of Contents
For this notebook, we will download WFC3/UVIS images of NGC 6791 from program 12692 in the F606W filter acquired in Visit 01. The three FLC images have been processed by the HST WFC3 pipeline (calwf3), which includes bias subtraction, dark current and CTE correction, cosmic-ray rejection, and flat-fielding:
ibwb01xqq_flc.fits
ibwb01xrq_flc.fits
ibwb01xxq_flc.fits
MAST queries may be done using query_criteria, where we specify:
\(\rightarrow\) obs_id, proposal_id, and filters
MAST data products may be downloaded by using download_products, where we specify:
\(\rightarrow\) products = calibrated (FLT, FLC) or drizzled (DRZ, DRC) files
\(\rightarrow\) type = standard products (CALxxx) or advanced products (HAP-SVM)
obs_ids = ['ibwb01*']
props = ['12692']
filts = ['F606W']
obsTable = Observations.query_criteria(obs_id=obs_ids, proposal_id=props, filters=filts)
products = Observations.get_product_list(obsTable)
data_prod = ['FLC'] # ['FLC','FLT','DRC','DRZ']
data_type = ['CALWF3'] # ['CALACS','CALWF3','CALWP2','HAP-SVM']
Observations.download_products(products, productSubGroupDescription=data_prod, project=data_type, cache=True)
Move to the files to the local working directory:
for flc in glob.glob('./mastDownload/HST/*/*flc.fits'):
flc_name = os.path.basename(flc)
os.rename(flc, flc_name)
if os.path.exists('mastDownload'):
shutil.rmtree('mastDownload')
1.1 Check image header data#
The cell below retrieves values for specific keywords from the first and second extensions of the image header. We see that the 1st exposure is 30 seconds and the 2nd and 3rd exposures are 360 seconds in duration. The 3rd exposure is dithered by approximately 82 arcseconds in the y-direction which is approximately the width of one WFC3 UVIS detector chip. A full list of header keywords is available in Section 2.4 of the WFC3 Data Handbook.
paths = sorted(glob.glob('*flc.fits'))
data = []
keywords_ext0 = ["ROOTNAME", "ASN_ID", "TARGNAME", "DETECTOR", "FILTER", "exptime",
"ra_targ", "dec_targ", "postarg1", "postarg2", "DATE-OBS"]
keywords_ext1 = ["orientat"]
for path in paths:
path_data = []
for keyword in keywords_ext0:
path_data.append(fits.getval(path, keyword, ext=0))
for keyword in keywords_ext1:
path_data.append(fits.getval(path, keyword, ext=1))
data.append(path_data)
keywords = keywords_ext0 + keywords_ext1
table = Table(np.array(data), names=keywords, dtype=['str', 'str', 'str', 'str', 'str', 'f8', 'f8', 'f8', 'f8', 'f8', 'str', 'f8'])
table['exptime'].format = '7.1f'
table['ra_targ'].format = table['dec_targ'].format = '7.4f'
table['postarg1'].format = table['postarg2'].format = '7.3f'
table['orientat'].format = '7.2f'
table
1.2 Inspect the Alignment#
Check the active WCS solution in the image header. If the image is aligned to a catalog, the list the number of matches and the fit RMS in milli-arcseconds. Next, convert the fit RMS values to pixels for comparison with the alignment results performed later in this notebook using Tweakreg.
ext_0_keywords = ['DETECTOR', 'EXPTIME'] # extension 0 keywords.
ext_1_keywords = ['WCSNAME', 'NMATCHES', 'RMS_RA', 'RMS_DEC'] # extension 1 keywords.
# Define the detector plate scales in arcsec per pixel.
DETECTOR_SCALES = {
'IR': 0.1283,
'UVIS': 0.0396,
'WFC': 0.05
}
formatted_data = {}
column_data = defaultdict(list)
for fits_file in sorted(glob.glob('*flc.fits')):
column_data['filename'].append(fits_file)
header0 = fits.getheader(fits_file, 0)
header1 = fits.getheader(fits_file, 1)
for keyword in ext_0_keywords:
column_data[keyword].append(header0[keyword])
for keyword in ext_1_keywords:
if 'RMS' in keyword:
value = np.around(header1[keyword], decimals=1)
else:
value = header1[keyword]
column_data[keyword].append(value)
for keyword in ['RMS_RA', 'RMS_DEC']:
rms_value = header1[keyword] / 1000 / DETECTOR_SCALES[header0['DETECTOR']]
column_data[f'{keyword}_pix'].append(np.round(rms_value, decimals=2))
wcstable = Table(column_data)
wcstable
2. Query the Reference Catalogs#
Table of Contents
Now that we have the FITS images, we will download the reference catalogs from both SDSS and Gaia using astroquery.
2.1 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.
RA = table['ra_targ'][0]
Dec = table['dec_targ'][0]
2.2 SDSS Query#
We now give the RA and Dec values to an astropy SkyCoord object, which we will pass to the SDSS query. Additionally, we use an astropy Quantity object to create a radius for the SDSS query in physical units. We set the radius to 5 arcminutes to comfortably cover the area of our images. For reference, the UVIS detector field of view is approximately 2.7’\(\times\)2.7’ and a y-dither of 82 arcseconds covers a total area on the sky of approximately 2.7’\(\times\)4.1’.
coord = SkyCoord(ra=RA, dec=Dec, unit=(u.deg, u.deg))
radius = Quantity(3., u.arcmin)
Next, we 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, as well as a g-band magnitude if we want to remove very dim and/or bright objects from the catalog, as those are likely measured poorly. Details on selecting objects by magnitude may be found in the ‘Gaia_alignment’ notebook. Many other fields are available in the SDSS query and are documented on this webpage.
sdss_query = SDSS.query_region(coord, radius=radius, spectro=False, fields=['ra', 'dec', 'g'])
sdss_query
We then need to save the catalog to use with TweakReg. As the returned value of the query is an astropy.Table, saving the file is straightforward:
sdss_query.write('sdss.cat', format='ascii.commented_header', overwrite=True)
2.3 Gaia Query#
Similarly to SDSS, we can query Gaia catalogs for our target via astroquery.gaia. This may be preferable in many cases because the spaced-based astrometry from Gaia is generally very accurate and precise. We can use the same coord and radius search parameters from our earlier SDSS query.
gaia_query = Gaia.query_object_async(coordinate=coord, radius=radius)
gaia_query
This query has returned very large number of columns, so we select the most useful columns to make the catalog easier to use with TweakReg.
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', overwrite=True)
3. Align with TweakReg#
Table of Contents
With the catalogs downloaded and the headers populated, we now 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 case, we set the updatehdr parameter to False to avoid altering the WCS of the images until we are satisfied with the alignment by inspecting both the TweakReg shift file and the astrometric residual plots.
3.1 SDSS Alignment#
refcat = 'sdss.cat'
wcsname = 'SDSS' # Specify the WCS name for this alignment.
cw = 3.5 # 2x the FWHM of the PSF = 3.5 pixels for ACS/WFC, WFC3/UVIS or 2.5 for WFC3/IR.
tweakreg.TweakReg('*flc.fits', # Pass the input images to tweakreg.
imagefindcfg={'threshold': 500., 'conv_width': cw}, # Detection parameters that vary for different datasets.
refcat=refcat, # Use user supplied catalog (Gaia).
interactive=False, # Don't show the interactive interface.
see2dplot=False, # Suppress additional figures in this test.
shiftfile=True, # Save output shift file to examine shifts later.
outshifts='SDSS_shifts.txt', # Name of the shift file that will be saved.
wcsname=wcsname, # Give our WCS a new name defined above.
reusename=True, # Overwrite the WCS if it exists.
ylimit=1.5,
fitgeometry='rscale', # Allow for translation, rotation, and plate scaling.
updatehdr=False) # Don't update the header with new WCS until later.
clear_output()
# If the alignment is unsuccessful then stop the notebook.
with open('SDSS_shifts.txt', 'r') as shift:
for line_number, line in enumerate(shift, start=1):
if "nan" in line:
raise ValueError('nan found in line {} in shift file'.format(line_number))
else:
continue
3.2 Inspect the shift file and fit residuals#
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]
shift_table = Table.read('SDSS_shifts.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
Show the 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)
We can see in the plots above that the RMS is fairly large at about 0.5 pixels. This is generally considered a poor fit, with the desired RMS being < 0.1 pixels depending on the number of sources. This is likely because the SDSS astrometric precision is insufficient for a high-quality 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.
3.3 Gaia Alignment#
refcat = 'gaia.cat'
wcsname = 'Gaia' # Specify the WCS name for this alignment.
cw = 3.5 # 2x the FWHM of the PSF = 3.5 pixels for ACS/WFC, WFC3/UVIS or 2.5 for WFC3/IR.
tweakreg.TweakReg('*flc.fits', # Pass the input images to tweakreg.
imagefindcfg={'threshold': 500., 'conv_width': cw}, # Detection parameters that vary for different datasets.
refcat=refcat, # Use user supplied catalog (Gaia).
interactive=False, # Don't show the interactive interface.
see2dplot=False, # Suppress additional figures in this test.
shiftfile=True, # Save output shift file to examine shifts later.
outshifts='Gaia_shifts.txt', # Name of the shift file that will be saved.
wcsname=wcsname, # Give our WCS a new name defined above.
reusename=True, # Overwrite the WCS if it exists.
ylimit=0.3, # The ylimit for the residuals plot.
fitgeometry='rscale', # Allow for translation, rotation, and plate scaling.
searchrad=0.1, # With many sources use a smaller search radius for convergence.
updatehdr=False) # Don't update the header with new WCS until later.
clear_output()
# If the alignment is unsuccessful then stop the notebook.
with open('Gaia_shifts.txt', 'r') as shift:
for line_number, line in enumerate(shift, start=1):
if "nan" in line:
raise ValueError('nan found in line {} in shift file'.format(line_number))
else:
continue
3.4 Inspect the shift file and fit residuals#
We can similarly look at the shift file from alignment to the Gaia catalog. As expected, the Gaia catalog does quite a bit better, with rms residuals ~0.05 pixels.
The columns are: Filename, X Shift [pixels], Y Shift [pixels], Rotation [degrees], Scale, X RMS [pixels], Y RMS [pixels]
shift_table = Table.read('Gaia_shifts.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
Print the Number of Gaia Matches per image from the TweakReg output.
rootname1 = 'ibwb01xqq'
rootname2 = 'ibwb01xrq'
rootname3 = 'ibwb01xxq'
match_tab1 = ascii.read(rootname1+'_flc_catalog_fit.match') # load match file in astropy table
match_tab2 = ascii.read(rootname2+'_flc_catalog_fit.match') # load match file in astropy table
match_tab3 = ascii.read(rootname3+'_flc_catalog_fit.match') # load match file in astropy table
print(rootname1+': Number of Gaia Matches =', len(match_tab1))
print(rootname2+': Number of Gaia Matches =', len(match_tab2))
print(rootname3+': Number of Gaia Matches =', len(match_tab3))
Compare the RMS and number of matches from TweakReg with the MAST alignment results.
wcstable
The new Gaia solution looks better than the MAST solution, so we rerun TweakReg below and update the image headers with the new WCS solution.
For more details on absolute astrometry in HST images, see Section 4.5 in the DrizzlePac Handbook. Show the 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)
3.5 Overplot Matched Sources on the Image#
Now, let’s plot the sources that were matched between all images on the “bottom” UVIS chip. This is referred to as chip 2 or SCI, 1 or extension 1. The cell below reads in the *_fit.match file as an astropy table. Unfortunately, it doesn’t automatically name columns so we look at the header to grab the columns.
To verify that TweakReg obtained an acceptable fit between matched source catalogs, it is useful to inspect the results before updating the image header WCS. In the figure below, sources that are matched with Gaia are overplotted on one of the input FLC frames with the two chips merged together.
It can be useful to check that TweakReg locked onto stars and not hot pixels or other detector artifacts before proceeding to update the image header.
# Choose one of the files to view.
# rootname = 'ibwb01xqq'
rootname = 'ibwb01xrq'
# rootname = 'ibwb01xxq'
plt.figure(figsize=(7, 7), dpi=140)
chip1_data = fits.open(rootname+'_flc.fits')['SCI', 2].data
chip2_data = fits.open(rootname+'_flc.fits')['SCI', 1].data
fullsci = np.concatenate([chip2_data, chip1_data])
z1, z2 = ZScaleInterval().get_limits(fullsci)
plt.imshow(fullsci, cmap='Greys', origin='lower', vmin=z1, vmax=z2)
match_tab = ascii.read(rootname+'_flc_catalog_fit.match') # Load the match file in astropy table.
match_tab_chip1 = match_tab[match_tab['col15'] == 2] # Filter the table for sources on chip 1.
match_tab_chip2 = match_tab[match_tab['col15'] == 1] # Filter the table for sources on chip 2.
x_cord1, y_cord1 = match_tab_chip1['col11'], match_tab_chip1['col12']
x_cord2, y_cord2 = match_tab_chip2['col11'], match_tab_chip2['col12']
plt.scatter(x_cord1, y_cord1+2051, s=50, edgecolor='r', facecolor='None', label='Matched Sources')
plt.scatter(x_cord2, y_cord2, s=50, edgecolor='r', facecolor='None')
plt.title(f'Matched sources UVIS to Gaia: N = {len(match_tab)}', fontsize=14)
plt.show()
3.6 Update header with optimal parameters#
Once you’ve decided on the optimal set of parameters for your alignment, it is safe to update the header of the input data, (assuming that the new solution is better than the MAST alignment.) 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'
wcsname = 'Gaia' # Specify the WCS name for this alignment
cw = 3.5 # 2x the FWHM of the PSF = 3.5 for ACS/WFC, WFC3/UVIS or 2.5 for WFC3/IR
tweakreg.TweakReg('*flc.fits', # Pass the input images to tweakreg.
imagefindcfg={'threshold': 500., 'conv_width': cw}, # Detection parameters that vary for different datasets.
refcat=refcat, # Use user supplied catalog (Gaia).
interactive=False, # Don't show the interactive interface.
see2dplot=False, # Suppress additional figures in this test.
shiftfile=True, # Save output shift file to examine shifts later.
outshifts='Gaia_shifts.txt', # Name of the shift file that will be saved.
wcsname=wcsname, # Give our WCS a new name defined above.
reusename=True, # Overwrite the WCS if it exists.
ylimit=0.3, # The ylimit for the residuals plot.
fitgeometry='rscale', # Allow for translation, rotation, and plate scaling.
searchrad=0.1, # With many sources use a smaller search radius for convergence.
updatehdr=True) # Update the header with new WCS solution.
clear_output()
4. Combine the Images using AstroDrizzle#
Table of Contents
While the three sets of FLC files are now aligned, in this example 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 the signal transitions from only being present in short exposures near the core to being present at larger radii in the longer exposures. The result is a discontinuity in the PSF radial profile and a resulting flux that is too low in some regions. For photometry of saturated objects, the short exposures should be drizzled separately from the long exposures.
Next, we combine the images through drizzling, and retrieve some recommended values for this process from the MDRIZTAB reference file. The parameters in this file are different for each detector and are based on the number of input frames. These are a good starting point for drizzling and may be adjusted based on your science goals.
input_images = sorted(glob.glob('ibwb01x[rx]q_flc.fits'))
mdz = fits.getval(input_images[0], 'MDRIZTAB', ext=0).split('$')[1]
print('Searching for the MDRIZTAB file:', mdz)
!crds sync --hst --files {mdz} --output-dir {os.environ['iref']}
def get_vals_from_mdriztab(input_images,
kw_list=[
'driz_sep_bits',
'combine_type',
'driz_cr_snr',
'driz_cr_scale',
'final_bits']):
'Get only selected parameters from the MDRIZTAB.'
mdriz_dict = getMdriztabPars(input_images)
requested_params = {}
print('Outputting the following parameters:')
for k in kw_list:
requested_params[k] = mdriz_dict[k]
print(k, mdriz_dict[k])
return requested_params
selected_params = get_vals_from_mdriztab(input_images)
Note that the parameter final_bits= ‘96’ is equivalent to final_bits= ‘64, 32’ because they are added in a bit-wise manner.
# To override any of the above values:
selected_params['driz_sep_bits'] = '256, 64, 32'
selected_params['final_bits'] = '256, 64, 32'
# selected_params['combine_type'] = 'median'
# selected_params['driz_cr_snr'] = '4.0 3.5'
# selected_params['driz_cr_scale'] = '1.5 1.2'
astrodrizzle.AstroDrizzle(input_images,
output='f606w',
preserve=False,
clean=True,
build=False,
context=False,
skymethod='match',
**selected_params)
clear_output()
Display the combined science and weight images. In this example there appear to be three detector chips due to the large dither that was used in the observations. We note that the exposures overlapped in the center region, as evidenced by the increased value of the weight map corresponding to a longer effective exposure time. In addition, the resulting mosaic is better cleaned of artifacts and cosmic rays in this region due to the multiple exposures, as seen by the uniform background in the science mosaic.
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.0, vmax=1000, cmap='Greys_r', origin='lower')
ax1.title.set_text('Science Image')
ax2.title.set_text('Weight Image')
ax1.set_xlabel('X (pixels)')
ax1.set_ylabel('Y (pixels)')
ax2.set_xlabel('X (pixels)')
ax2.set_ylabel('Y (pixels)')
For more details on inspecting the drizzled products after reprocessing, see Section 7.3 in the DrizzlePac Handbook.
Conclusions#
Table of Contents
Many other services have interfaces for querying catalogs that can also be used to align HST images. In general, Gaia works very well for HST due to its 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 a way to make data comparable across many epochs/detectors/observatories, and in many cases, makes the alignment easier to complete.
About this Notebook#
Created: 14 Dec 2018; V. Bajaj
Updated: 03 Jul 2024; M. Revalski, V. Bajaj, & J. Mack
Source: GitHub spacetelescope/hst_notebooks
Additional Resources#
Below are some additional resources that may be helpful. Please send any questions through the HST Help Desk, selecting the DrizzlePac category.
Citations#
If you use Python packages such as astropy, astroquery, drizzlepac, matplotlib, or numpy for published research, please cite the authors.
Follow these links for more information about citing various packages:
Top of Page
