Extra care must be taken when using AstroDrizzle to combine observations from detectors comprised of multiple chips of varying sensitivity. AstroDrizzle works with calibrated images in units of counts (electrons or Data Numbers) or count rates and not in units of flux. It assumes that all input frames can be converted to physical flux units using a single inverse-sensitivity factor, recorded in the FITS image headers as PHOTFLAM, and the output drizzled product simply copies the PHOTFLAM keyword value from the first input image. When this occurs, the inverse-sensitivity will vary across the final drizzled product, and users will need to keep track of which sources fell on which chip when doing photometry. Moreover, varying detector sensitivities will affect the cosmic-ray rejection algorithm used by AstroDrizzle, and this may result in the misidentification of some good pixels as cosmic rays.
This is a typical situation when drizzle-combining images from HST instruments with different chip sensitivities, e.g. Wide Field and Planetary Camera 2 (WFPC2). For more detail, see the section on Gain Variation under 'Position-Dependent Photometric Corrections' in the WFPC2 Data Handbook. As a result, each of the four chips requires a unique PHOTFLAM header keyword value. A similar situation may occur when drizzle-combining observations taken over a span of several years as detector's sensitivity declines over time, see e.g. ACS ISR 2016-03.
One approach is to rescale the input data so that AstroDrizzle can properly assume the images/chips have the same sensitivity; that is, a single PHOTFLAM value can be used to convert re-scaled image counts (or count-rates) to physical integrated flux units. The photeq task in Drizzlepac automates this image intensity rescaling to a single inverse-sensitivity factor PHOTFLAM.
In this example notebook, archival WFPC2 images are used to demonstrate advanced reprocessing using TweakReg and AstroDrizzle for alignment and image combination. The notebook is based on a prior WFPC2 example but includes additional information about equalizing the chip sensitivities prior to combining.
NOTE: It is important to note that photeq only adjusts image counts so that integrated physical fluxes can be obtained using a single PHOTFLAM. It does nothing to account for different throughtputs at different wavelengths.
import shutil
import glob
import os
import subprocess
import matplotlib.pyplot as plt
from astropy.io import fits
from astroquery.mast import Observations
from stwcs.updatewcs import updatewcs
from drizzlepac import tweakreg, astrodrizzle, photeq
from stsci.skypac import skymatch
# ONLY needed for the simulation section:
import numpy as np
from stwcs.wcsutil import HSTWCS
from drizzlepac.wfpc2Data import WFPC2_GAINS
%matplotlib inline
This example uses WFPC2 observations of Messier 2 in the F814W filter. The data come from GO proposal 11100 "Two new 'bullets' for MOND: revealing the properties of dark matter in massive merging clusters". Four images were acquired using a 4-pt dither box pattern, followed by two images offset with a dither-line pattern.
The data are downloaded using the astroquery API to access the MAST archive. The astroquery.mast documentation has more examples for how to find and download data from MAST.
# Retrieve the observation information.
if os.path.isdir('mastDownload'):
shutil.rmtree('mastDownload')
obs_table = Observations.query_criteria(obs_id='ua0605*', filters='F814W', obstype='ALL')
products = Observations.get_product_list(obs_table)
# Download only the ua0605*_c0m.fits and ua0605*_c1m.fits (DQ) images:
Observations.download_products(products, mrp_only=False, productSubGroupDescription=['C0M', 'C1M'], extension='fits')
# Move the files from the mastDownload directory to the current working
# directory and make a backup of the files.
fits_files = glob.glob('mastDownload/HST/ua*/ua*c?m.fits')
for fil in fits_files:
base_name = os.path.basename(fil)
if os.path.isfile(base_name):
os.remove(base_name)
shutil.move(fil, '.')
# Delete the mastDownload directory and all subdirectories it contains.
shutil.rmtree('mastDownload')
WFPC2 images downloaded from the archive contain World Coordinate System (WCS) information based on an older-style description of image distortions. Before these images can be processed with drizzlepac, their WCS must be converted to a new format. This can be achieved using updatewcs() function from the stwcs package. More details may be found here: 'Making WFPC2 Images Compatible with AstroDrizzle'. Note that updatewcs is no longer a parameter in AstroDrizzle or TweakReg and must be run separately before processing the data.
First we download the reference files from the CRDS website. See the initialization notebook in this repository for more information.
os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'
os.environ['CRDS_PATH'] = os.path.abspath(os.path.join('.', 'reference_files'))
subprocess.check_output('crds bestrefs --files ua0605*_c0m.fits --sync-references=1 --update-bestrefs', shell=True, stderr=subprocess.DEVNULL)
os.environ['uref'] = os.path.abspath(os.path.join('.', 'reference_files', 'references', 'hst', 'wfpc2')) + os.path.sep
NOTE: This next step may raise warnings because the Astrometry database is in progress and currently does not cover the entire sky. Please ignore these warnings. The WCS will still be updated.
updatewcs('ua*c0m.fits', use_db=True)
In a later section we will generate simulated data to illustrate the effects of drizzling WFPC2 images without sensitivity equalization. For that purpose we will need a copy of an original image that has the original inverse-sensitivity values (PHOTFLAM) in their headers. Here we create a backup copy of the first image.
NOTE: This step is needed for illustration purpose in this notebook only. It is not needed when processing data.
orig_image = glob.glob('ua*c0m.fits')[0]
backup_image = 'simulation.fits'
if os.path.isfile(backup_image):
os.remove(backup_image)
shutil.copy2(orig_image, backup_image)
Due to small pointing errors, the image header WCS typically needs to be updated in order to achieve the best drizzle-combined products. The expected pointing accuracy for various observing scenerios is summarized in the DrizzlePac Handbook Appendix B. Input images must first be aligned so that when the coordinates of a given object (in detector space) are converted to sky coordinates (using the WCS), that object's sky coordinates must be approximately equal in each frame.
The DrizzlePac task TweakReg may be used to correct for any errors in the image header WCS. First, TweakReg finds sources in each image, matches sources in common across images, and finds a separate linear transformation to align each image. TweakReg then computes a new WCS for each image based on this linear transformation.
Here we show a basic image alignment procedure. For a more detailed illustration of image alignment, please refer to other example notebooks such as the mosaic example in this repository.
tweakreg.TweakReg('ua*c0m.fits', updatehdr=True, reusename=True, interactive=False,
conv_width=3.0, threshold=300.0, peakmin=100, peakmax=10000)
This step adjusts image data values so that all images and chips appear (to AstroDrizzle) to have a single inverse sensitivity (PHOTFLAM). This can be achieved using the photeq task in Drizzlepac. This task adjusts image data so that when these data are multiplied by the same single PHOTFLAM value, the correct flux is obtained.
photeq.photeq(files='ua*_c0m.fits', ref_phot_ext=3, readonly=False)
In the above command, we instruct photeq to "equalize" all chips of all input images using the PHOTFLAM for the WF3 chip (ref_phot_ext=3), using the first image as a reference. This reference PHOTFLAM value is reported in the log file:
REFERENCE VALUE FROM FILE: 'ua060502m_c0m.fits['SCI',1]'
REFERENCE 'PHOTFLAM' VALUE IS: 2.507987e-18Upon the completion, photeq will not only adjust image data but also update PHOTFLAM values for all chips to this specific reference value.
All four chips are now drizzled together with an output pixel scale set to that of the WF chips:
astrodrizzle.AstroDrizzle('ua*c0m.fits',
preserve=False,
driz_sep_bits='8,1024',
driz_sep_wcs=True,
driz_sep_scale=0.0996,
combine_type='median',
driz_cr_snr='5.5 3.5',
driz_cr_scale='2.0 1.5',
final_fillval=None,
final_bits='8,1024',
final_wcs=True,
final_scale=0.0996)
drz = fits.getdata('final_drz_sci.fits')
plt.figure(figsize=(15, 15))
plt.imshow(drz, cmap='gray', vmin=-0.1, vmax=0.5, origin='lower')
The effect of drizzling images with different detector sensitivies, while tangible, is sometimes difficult to visualize in noisy data, especially when drizzling multiple dithered images that can blur the borders between chips.
In this section we produce a simple simulation of observing a constant intensity "blank sky". We then make a copy of this image and apply sensitivity equalization to it. Finally we drizzle both images and compare them side-by-side.
In this simple simulation we assume only Poisson noise.
with fits.open(backup_image, mode='update') as h:
# get chip inverse-sensitivity:
phot1 = h[1].header['PHOTFLAM']
phot2 = h[2].header['PHOTFLAM']
phot3 = h[3].header['PHOTFLAM']
phot4 = h[4].header['PHOTFLAM']
# get chip WCS:
w1 = HSTWCS(h, ext=1)
w2 = HSTWCS(h, ext=2)
w3 = HSTWCS(h, ext=3)
w4 = HSTWCS(h, ext=4)
ref_pscale = w4.idcscale
# get chip gain:
cmdgain = h[0].header['ATODGAIN']
gain1 = WFPC2_GAINS[1][cmdgain][0]
gain2 = WFPC2_GAINS[2][cmdgain][0]
gain3 = WFPC2_GAINS[3][cmdgain][0]
gain4 = WFPC2_GAINS[4][cmdgain][0]
# final drizzle scale:
scale = 0.0996
# simulated sky level ("true" sky is constant):
sky = 10 * phot3
# simulate observed counts assuming only Poisson noise:
conv1a = gain1 * (w1.idcscale / ref_pscale)**2 / phot1
conv1b = (gain4 / gain1**2) * (ref_pscale / scale)**2
conv2a = gain2 * (w2.idcscale / ref_pscale)**2 / phot2
conv2b = (gain4 / gain2**2) * (ref_pscale / scale)**2
conv3a = gain3 * (w3.idcscale / ref_pscale)**2 / phot3
conv3b = (gain4 / gain3**2) * (ref_pscale / scale)**2
conv4a = gain4 * (w4.idcscale / ref_pscale)**2 / phot4
conv4b = (1.0 / gain4) * (ref_pscale / scale)**2
h[1].data[:, :] = np.random.poisson(conv1a * sky, h[1].data.shape) * conv1b
h[2].data[:, :] = np.random.poisson(conv2a * sky, h[2].data.shape) * conv2b
h[3].data[:, :] = np.random.poisson(conv3a * sky, h[3].data.shape) * conv3b
h[4].data[:, :] = np.random.poisson(conv4a * sky, h[4].data.shape) * conv4b
# make a copy of this file:
photeq_image = 'simulation_eq.fits'
if os.path.isfile(photeq_image):
os.remove(photeq_image)
shutil.copy2(backup_image, photeq_image)
# apply equalization to the image copy:
photeq.photeq(files=photeq_image, ref_phot_ext=3, readonly=False)
astrodrizzle.AstroDrizzle(
backup_image,
output='nonequalized.fits',
stepsize=1,
preserve=False,
restore=False,
in_memory=True,
context=False,
build=False,
static=False,
skysub=False,
median=False,
blot=False,
driz_cr=False,
final_fillval=None,
final_bits='',
final_wcs=True,
final_scale=scale)
astrodrizzle.AstroDrizzle(
photeq_image,
output='equalized.fits',
stepsize=1,
preserve=False,
restore=False,
in_memory=True,
context=False,
build=False,
static=False,
skysub=False,
median=False,
blot=False,
driz_cr=False,
final_fillval=None,
final_bits='',
final_wcs=True,
final_scale=scale)
drz_noneq = fits.getdata('nonequalized_sci.fits')
drz_eq = fits.getdata('equalized_sci.fits')
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(24, 10))
im1 = ax1.imshow(drz_noneq, cmap='gray', vmin=0.88, vmax=1.08, origin='lower')
ax1.set_title('Original Simulation (Not equalized)')
im2 = ax2.imshow(drz_eq, cmap='gray', vmin=0.88, vmax=1.08, origin='lower')
ax2.set_title('Equalized Simulated Image')
x1 = ax1.get_position().get_points().flatten()[0]
x2 = ax2.get_position().get_points().flatten()[2] - x1
ax_cbar = fig.add_axes([x1, 0, x2, 0.03])
plt.colorbar(im1, cax=ax_cbar, orientation='horizontal')
Author: M. Cara, STScI Data Analysis Tools Branch
Updated: December 14, 2018