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
In this example notebook, archival WFPC2 images are used to demonstrate advanced reprocessing using
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
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-651af469802b> in <module> 9 from astroquery.mast import Observations 10 from stwcs.updatewcs import updatewcs ---> 11 from drizzlepac import tweakreg, astrodrizzle, photeq 12 from stsci.skypac import skymatch 13 /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:
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.
# 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.
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') 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.
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-18
Upon 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.header['PHOTFLAM'] phot2 = h.header['PHOTFLAM'] phot3 = h.header['PHOTFLAM'] phot4 = h.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.header['ATODGAIN'] gain1 = WFPC2_GAINS[cmdgain] gain2 = WFPC2_GAINS[cmdgain] gain3 = WFPC2_GAINS[cmdgain] gain4 = WFPC2_GAINS[cmdgain] # 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.data[:, :] = np.random.poisson(conv1a * sky, h.data.shape) * conv1b h.data[:, :] = np.random.poisson(conv2a * sky, h.data.shape) * conv2b h.data[:, :] = np.random.poisson(conv3a * sky, h.data.shape) * conv3b h.data[:, :] = np.random.poisson(conv4a * sky, h.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() x2 = ax2.get_position().get_points().flatten() - x1 ax_cbar = fig.add_axes([x1, 0, x2, 0.03]) plt.colorbar(im1, cax=ax_cbar, orientation='horizontal')