Drizzling new WFPC2 FLT data products with chip-normalization

---

This notebook requires creating and activating a virtual environment using the requirements file in this notebook's repository. Please also review the README file before using the notebook.

Table of Contents

  1. Introduction
  2. Download the Data
  3. Update the WCS
  4. Align the Images
  5. Drizzle-combine the images
  Conclusions

1. Introduction

Table of Contents

This notebook shows how to work with a new type of WFPC2 calibrated data product in MAST. The new data products combine the previously separate files containing the science array c0m.fitsand data quality array c1m.fits into a new file with suffix flt.fits, similar to calibrated images for ACS and WFC3.

These flt files have now been corrected for differences in the inverse sensitivity the science arrays of each chip using the software ‘photeq’ so that a single PHOTFLAM (or PHOTFNU) value may be used for photometry.

The new flt files include four key differences from the older style c0m, c1m WFPC2 files:

1.) The SCI & DQ arrays are merged into a single file with multiple extensions.
2.) The SCI arrays are corrected for different 'inverse sensitivities' across chips.
3.) The SCI arrays from different chips may now be drizzled together in a mosaic.
4.) The SCI header WCS includes improved absolute astromentry (i.e. WCSNAME='*GAIA*').

This example uses WFPC2 observations of Omega Centauri (NGC 5139) in the F555W filter taken from CAL proposal 8447. The two exposures have been dithered to place the same stars on the WF2 and WF4 chips to measure relative photometry across the two chips. After using TweakReg to realign the images, we show how to combine the images with AstroDrizzle. We inspect the resulting science and weight images and compare the results with those derived from the older MAST data products.

Acronyms:

• Hubble Space Telescope (HST)

• Wide Field Camera 3 (WFC3)

• Advanced Camera for Surveys (ACS)

• Wide Field and Planetary Camera 2 (WFPC2)

• Mikulski Archive for Space Telescopes (MAST)

• Flexible Image Transport System (FITS)

• World Coordinate System (WCS)

• Hubble Advanced Products (HAP)

• Single Visit Mosaic (SVM)

• Planetary Chip (PC)

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 - copy files between folders

• numpy - math and array functions

• matplotlib - create graphics

  • pyplot - make figures and graphics

• astroquery - download astronomical data

  • mast.Observations - query MAST database

• astropy - model fitting and file handling

  • io.fits - import FITS files

  • table.Table - tables without physical units

  • visualization.LogStretch - display image normalization

  • visualization.ImageNormalize - display image normalization

• ipython - cell formatting and interactives

  • display.Image - display image files

  • display.clear_output - clear text from cells

• photutils - aperture photometry tools

  • photutils.aperture - aperture photometry tools

    • aperture_photometry - measure aperture photometry

    • CircularAperture - measure photometry in circle

    • CircularAnnulus - measure photometry in annulus

    • ApertureStats - calculate aperture statistics

• drizzlepac - combine HST images into mosaics

  • astrodrizzle - combine exposures by drizzling

  • tweakreg - align exposures to a common WCS

• stwcs - WCS based distortion models and transformations

  • updatewcs - modify WCS in HST file headers

import os
import glob
import shutil
import numpy as np
import matplotlib.pyplot as plt
from astroquery.mast import Observations
from astropy.table import Table
from astropy.io import fits
from astropy.visualization import LogStretch, ImageNormalize
from IPython.display import Image, clear_output
from photutils.aperture import CircularAnnulus, CircularAperture, ApertureStats, aperture_photometry
from drizzlepac import tweakreg, astrodrizzle
from stwcs.updatewcs import updatewcs

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

The following task in the stsci.skypac package can be run with TEAL:
                                    skymatch                                    

The following tasks in the drizzlepac package can be run with TEAL:
    astrodrizzle       config_testbed      imagefindpars           mapreg       
       photeq            pixreplace           pixtopix            pixtosky      
  refimagefindpars       resetbits          runastrodriz          skytopix      
     tweakback            tweakreg           updatenpol

2. Download the Data

Table of Contents

The WFPC2 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.

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\) data_prod = original c0m, c1m files and new flt calibrated images and their corresponding drw drizzled files

    \(\rightarrow\) data_type = standard products CALWFPC2 or Hubble Advanced Products: ‘Single Visit Mosaics’ HAP-SVM

Depending on your connection speed this cell may take a few minutes to execute.

obs_ids = ['U5JX0108R', 'U5JX010HR']
filts = ['F555W']

obsTable = Observations.query_criteria(obs_id=obs_ids, filters=filts)
products = Observations.get_product_list(obsTable)

data_prod = ['C0M', 'C1M', 'FLT', 'DRW']                        
data_type = ['CALWFPC2'] # Options: ['CALWFPC2', 'HAP-SVM']
                           
Observations.download_products(products, productSubGroupDescription=data_prod, project=data_type)

clear_output()

The MAST pipeline uses the new flt files to create two new types of WFPC2 drizzled products. Unlike the original u*_drz.fits files, which were drizzled to an output scale of the WF chips (0.0996”/pix), the new drizzled files u*_drw.fits are drizzled to an output scale of the PC chip (0.0455”/pix), where the new drw suffix reflects the improved WCS.

A second drizzled product hst_*drz.fits combines all images in the same filter and may include improved relative alignment across filters in the same visit. These are referred to as Hubble Advanced Products and are listed below, but are beyond the scope of this notebook. For more details, see the HAP Webpage and the Drizzlepac Handbook.

Drizzled Image                             Dimensions    Scale     Input      Updated
                                                        ("/pix)    File       WCS?
______________________________________________________________________________________
u5jx010hr_drz.fits                         [1515, 1495]  0.0996   c?m.fits    N
u5jx010hr_drw.fits                         [4171, 4143]  0.0455   flt.fits    Y
hst_8447_01_wfpc2_pc_f555w_u5jx01_drz.fits [6258, 6215]  0.0455   flt.fits    Y    

Move files to the current working directory

fits_files = glob.glob('mastDownload/HST/*/*.fits')
for f in fits_files:
    base_name = os.path.basename(f)
    if os.path.isfile(base_name):
        os.remove(base_name)
    shutil.move(f, '.')   
shutil.rmtree('mastDownload')

Inspect the File Structure

While the flt.fits files use the same notation for different chips, EXT=('SCI',1) the order of the FITS extensions ('SCI', 'ERR', 'DQ') is different for WFC3/UVIS and WFPC2. For example:

    uvis_flt.fits[1]   =  uvis_flt.fits['SCI',1]   --> UVIS2, SCI
    uvis_flt.fits[2]   =  uvis_flt.fits['ERR',1]   --> UVIS2, ERR
    uvis_flt.fits[3]   =  uvis_flt.fits['DQ',1]    --> UVIS2, DQ
    uvis_flt.fits[4]   =  uvis_flt.fits['SCI',2]   --> UVIS1, SCI
    uvis_flt.fits[5]   =  uvis_flt.fits['ERR',2]   --> UVIS1, ERR
    uvis_flt.fits[6]   =  uvis_flt.fits['DQ',2]    --> UVIS1, DQ

    wfpc2_flt.fits[1]  = wfpc2_flt.fits['SCI',1)   -->  PC1, SCI
    wfpc2_flt.fits[2]  = wfpc2_flt.fits['SCI',2]   -->  WF2, SCI
    wfpc2_flt.fits[3]  = wfpc2_flt.fits['SCI',3]   -->  WF3, SCI
    wfpc2_flt.fits[4]  = wfpc2_flt.fits['SCI',4]   -->  WF4, SCI
    wfpc2_flt.fits[5]  = wfpc2_flt.fits['DQ',1]    -->  PC1, DQ
    wfpc2_flt.fits[6]  = wfpc2_flt.fits['ERR',1]   -->  PC1, ERR        
    wfpc2_flt.fits[7]  = wfpc2_flt.fits['DQ',2]    -->  WF2, DQ
    wfpc2_flt.fits[8]  = wfpc2_flt.fits['ERR',2]   -->  WF2, ERR
    wfpc2_flt.fits[9]  = wfpc2_flt.fits['DQ',3]    -->  WF3, DQ
    wfpc2_flt.fits[10] = wfpc2_flt.fits['ERR',3]   -->  WF3, ERR
    wfpc2_flt.fits[11] = wfpc2_flt.fits['DQ',4]    -->  WF4, DQ
    wfpc2_flt.fits[12] = wfpc2_flt.fits['ERR',4]   -->  WF4, ERR

fits.info('u5jx010hr_c0m.fits')

Filename: u5jx010hr_c0m.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     251   ()      
  1  SCI           1 ImageHDU       143   (800, 800)   float32   
  2  SCI           2 ImageHDU       132   (800, 800)   float32   
  3  SCI           3 ImageHDU       132   (800, 800)   float32   
  4  SCI           4 ImageHDU       132   (800, 800)   float32   

fits.info('u5jx010hr_flt.fits')

Filename: u5jx010hr_flt.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     241   ()      
  1  SCI           1 ImageHDU       206   (800, 800)   float32   
  2  SCI           2 ImageHDU       199   (800, 800)   float32   
  3  SCI           3 ImageHDU       199   (800, 800)   float32   
  4  SCI           4 ImageHDU       199   (800, 800)   float32   
  5  DQ            1 ImageHDU       107   (800, 800)   int16   
  6  ERR           1 ImageHDU       141   (800, 800)   float32   
  7  DQ            2 ImageHDU       107   (800, 800)   int16   
  8  ERR           2 ImageHDU       132   (800, 800)   float32   
  9  DQ            3 ImageHDU       107   (800, 800)   int16   
 10  ERR           3 ImageHDU       132   (800, 800)   float32   
 11  DQ            4 ImageHDU       107   (800, 800)   int16   
 12  ERR           4 ImageHDU       137   (800, 800)   float32   
 13  D2IMARR       1 ImageHDU        16   (800, 1)   float32   
 14  HDRLET        1 HeaderletHDU     21   ()      
 15  WCSCORR       1 BinTableHDU     59   28R x 24C   [40A, I, A, 24A, 24A, 24A, 24A, D, D, D, D, D, D, D, D, 24A, 24A, D, D, D, D, J, 40A, 128A]   
 16  HDRLET       16 HeaderletHDU     21   ()      
 17  HDRLET        3 HeaderletHDU     21   ()      

The drizzled products, on the other hand, contain 3 data extensions: SCI, WHT, and CTX.

fits.info('u5jx010hr_drw.fits')

Filename: u5jx010hr_drw.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     797   ()      
  1  SCI           1 ImageHDU       100   (4171, 4143)   float32   
  2  WHT           1 ImageHDU        89   (4171, 4143)   float32   
  3  CTX           1 ImageHDU        92   (4171, 4143)   int32   
  4  HDRTAB        1 BinTableHDU    557   4R x 274C   [9A, 3A, K, D, D, D, D, D, D, D, D, D, D, D, D, K, D, D, D, D, 8A, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, K, D, 8A, 18A, D, D, D, D, D, D, 8A, D, 6A, 19A, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, 23A, D, D, D, D, D, D, D, D, D, D, D, D, D, 12A, 12A, 8A, 18A, D, D, 19A, 10A, D, D, D, 2A, D, 3A, D, D, 7A, D, K, D, 6A, 9A, D, D, D, 4A, 18A, 3A, K, D, D, D, D, 8A, D, D, D, D, D, 23A, 1A, D, 23A, D, D, D, 3A, L, D, 5A, 21A, D, 6A, D, D, D, D, D, D, D, D, D, D, D, L, K, K, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, D, 8A, 8A, D, D, D, D, K, K, D, 12A, D, 37A, D, D, 1A, 1A, D, K, D, D, 1A, 1A, D, D, D, D, D, D, D, D, D, D, D, 4A, D, K, D, D, D, D, D, 23A, D, D, K, D, D, D, D, D, 2A, D, D, K, 20A, 1A, 1A, 2A, D, D, D, D, D, D, 4A, D, D, D, D, D, D, D, D, D, D, D, D, D, D, 4A, 18A, D, L, D, D, D, D, D, D, D, K, K, 6A, D, D, D, D, 1A, 12A, D, 3A, 8A, D, D, K, 21A, D, D]   

Compare photometric keywords from the handbook with the image headers

PHOTFLAM ‘inverse sensitivity’ values are listed in Table 5.1 of the WFPC2 DATA HANDBOOK Section 5.1 Photometric Zeropoint for the gain 7 setting (which is the default setting for most science programs). The PHOTFLAM values for gain 14 can be obtained by multiplying by the gain ratios: 1.987 (PC1), 2.003 (WF2), 2.006 (WF3), and 1.955 (WF4), where the values are from Holtzman et al. (1995).

Filter   PC1           WF2           WF3           WF4
-----    ---------     ---------     ---------     ---------
f555w 	 3.483e-18     3.396e-18     3.439e-18     3.507e-18 

The PHOTFLAM values include differences in the gain (e- DN-1) between chips, i.e. Table 4.2 in Section 4.4: “Read Noise and Gain Settings” in the 2008 WFPC2 Instrument Handbook, reproduced below:

Gain     PC1           WF2           WF3           WF4
-----    ---------     ---------     ---------     ---------
 "7"      7.12 ± 0.41   7.12 ± 0.41   6.90 ± 0.32   7.10 ± 0.39
 "15"    13.99 ± 0.63  14.50 ± 0.77  13.95 ± 0.63  13.95 ± 0.70

While WFC3 and ACS flt.fits data products are in units of electrons (WFC3/UVIS, ACS/WFC) or electrons/sec (WFC3/IR), WFPC2 images are in units of ‘Data Number’ (DN) or ‘COUNTS’ as reflected in the BUNIT keyword.

Here we inspect a subset of photometry keywords in the c0m and flt files:

with fits.open('u5jx0108r_c0m.fits') as h:
    phot1 = h[1].header['PHOTFLAM']
    phot2 = h[2].header['PHOTFLAM']
    phot3 = h[3].header['PHOTFLAM']
    phot4 = h[4].header['PHOTFLAM']
    gain = h[0].header['ATODGAIN']
    bunit = h[1].header['BUNIT']

print('PHOTFLAM for PC1: ', phot1, ' Factor Relative to WF3: ', round(phot1/phot3, 4))
print('PHOTFLAM for WF2: ', phot2, ' Factor Relative to WF3: ', round(phot2/phot3, 4))
print('PHOTFLAM for WF3: ', phot3, ' Factor Relative to WF3: ', round(phot3/phot3, 4))
print('PHOTFLAM for WF4: ', phot4, ' Factor Relative to WF3: ', round(phot4/phot3, 4))
print('Gain Value is: ', gain)    
print('Units are: ', bunit)    

PHOTFLAM for PC1:  3.482944e-18  Factor Relative to WF3:  1.0128
PHOTFLAM for WF2:  3.395806e-18  Factor Relative to WF3:  0.9874
PHOTFLAM for WF3:  3.439046e-18  Factor Relative to WF3:  1.0
PHOTFLAM for WF4:  3.506729e-18  Factor Relative to WF3:  1.0197
Gain Value is:  7.0
Units are:  COUNTS

Note the unique PHOTFLAM values for each chip in the original c0m data products (above).

In the new flt data products, the PHOTFLAM values are all set to the value for WF3 (below).

with fits.open('u5jx0108r_flt.fits') as h:
    phot1 = h[1].header['PHOTFLAM']
    phot2 = h[2].header['PHOTFLAM']
    phot3 = h[3].header['PHOTFLAM']
    phot4 = h[4].header['PHOTFLAM']
    gain = h[0].header['ATODGAIN']
    bunit = h[1].header['BUNIT']

print('PHOTFLAM for PC1: ', phot1, ' Factor Relative to WF3: ', round(phot1/phot3, 4))
print('PHOTFLAM for WF2: ', phot2, ' Factor Relative to WF3: ', round(phot2/phot3, 4))
print('PHOTFLAM for WF3: ', phot3, ' Factor Relative to WF3: ', round(phot3/phot3, 4))
print('PHOTFLAM for WF4: ', phot4, ' Factor Relative to WF3: ', round(phot4/phot3, 4))
print('Gain Value is: ', gain)    
print('Units are: ', bunit)  

PHOTFLAM for PC1:  3.439046e-18  Factor Relative to WF3:  1.0
PHOTFLAM for WF2:  3.439046e-18  Factor Relative to WF3:  1.0
PHOTFLAM for WF3:  3.439046e-18  Factor Relative to WF3:  1.0
PHOTFLAM for WF4:  3.439046e-18  Factor Relative to WF3:  1.0
Gain Value is:  7.0
Units are:  COUNTS

The new drw data products carry the WF3 PHOTFLAM values. While the BUNIT keyword says ‘counts’, it should say ‘counts/sec’

with fits.open('u5jx010hr_drw.fits') as h:
    phot1 = h[1].header['PHOTFLAM']
    gain = h[0].header['ATODGAIN']
    bunit = h[1].header['BUNIT']

print('PHOTFLAM for Combined image', phot1, ' Factor Relative to WF3: ', round(phot1/phot3, 4))
print('Gain Value is: ', gain)    
print('Units are: ', bunit)  

PHOTFLAM for Combined image 3.439046e-18  Factor Relative to WF3:  1.0
Gain Value is:  7.0
Units are:  counts

Because WF3 is the chip with the most accurate flux calibration, this was selected as the ‘reference’ chip and the inverse sensitivities in the new flt files have been equalized to that value. For example:

photeq.photeq(files='*_flt.fits', ref_phot_ext=3) # or ('sci',3)

Note that photeq adjusts the pixels values in the SCI array so that photometry can be obtained using a single PHOTFLAM. The ratio of the count rates in the chips will therefore be equal to the ratio of the PHOTFLAM values.

Compare the equation for computing STMAG using the new and the old workflow:

NEW: Use flt files

Drizzle the ‘flt.fits’ files and do photometry for all chips at once;

        STMAG = -21.1 -2.5*log (countrate * PHOTFLAM)  

OLD: Use c0m files

option A.) Drizzle the ‘c0m.fits’ files and do photometry for each chip separately;

        STMAG = -21.1 -2.5*log (countrate * PHOTFLAM_1,2,3,4)  

option B.) Run phot_eq on the ‘c0m.fits’ files prior to drizzling, normalizing to the WF3 chip. Use the first equation with a single PHOTFLAM value which equal to PHOTFLAM for WF3 in the second equation.

With the older c0m data products, caution must be taken when using AstroDrizzle to combine chips of different 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 value, 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 can affect the cosmic-ray rejection algorithm used by AstroDrizzle, and this may result in the misidentification of some good pixels as cosmic rays. Using the new FLT style files solves these issues by normalizing the chips to a single PHOTFLAM value.

Check image header data

In this section, we inspect keywords in the image headers which are useful for understanding the data. Note the PHOTMODE keyword in column 6 below, which lists the CCD gain. Additionally, we inspect the keywords WCSNAME and NMATCHES to determine whether the images have improved astrometry.

paths = sorted(glob.glob('*flt.fits'))
data = []
keywords_ext0 = ["ROOTNAME", "TARGNAME", "FILTNAM1", "EXPTIME"]
keywords_ext1 = ["ORIENTAT",  "PHOTMODE", "PHOTFLAM", "WCSNAME", "NMATCHES"]

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', 'f8', 'f8', 'str', 'f8', 'str', 'i8'])
table['EXPTIME'].format = '7.1f'  
table['ORIENTAT'].format = '7.2f'
table['PHOTFLAM'].format = '7.6g' 
table

Table length=2

ROOTNAME	TARGNAME	FILTNAM1	EXPTIME	ORIENTAT	PHOTMODE	PHOTFLAM	WCSNAME	NMATCHES
str32	str32	str32	float64	float64	str32	float64	str32	int64
u5jx0108r	NGC5139-POS2	F555W	16.0	162.26	WFPC2,1,A2D7,F555W,,CAL	3.43905e-18	IDC_ta81040lu-FIT_IMG_GAIAeDR3	140
u5jx010hr	NGC5139-POS2	F555W	16.0	162.28	WFPC2,1,A2D7,F555W,,CAL	3.43905e-18	IDC_ta81040lu-FIT_IMG_GAIAeDR3	112

Here we see that the images have been aligned to Gaia eDR3 with ~ 100 matched objects in the HST data.

3. Update the WCS

Table of Contents

Unlike the new flt products, the older WFPC2 c0m files 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 in the DrizzlePac Handbook, chapter 4.4.5 under section: ‘Making WFPC2 Images Compatible with AstroDrizzle’. NOTE: Running updatewcs is not necessary when working with flt files.

First we download the necessary reference files from the CRDS website.

os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'
os.environ['CRDS_PATH'] = os.path.abspath(os.path.join('.', 'reference_files'))
os.system('crds bestrefs --files u5jx01*_c0m.fits --sync-references=1 --update-bestrefs')
os.environ['uref'] = os.path.abspath(os.path.join('.', 'reference_files', 'references', 'hst', 'wfpc2')) + os.path.sep

clear_output()

NOTE: This next step may raise warnings which may be ignored.

updatewcs('u*c0m.fits', use_db=False)
clear_output()

While the new flt images have already been aligned to Gaia, we want to directly compare (ratio) the two sets of drizzled products. Thus we reset the WCS in the flt files to use the ‘distortion-only’ solution. When there are a large number of GAIA matches (\(\geq\) 20), users will want to use the GAIA-based WCS solutions for optimal alignment, especially when there is a large dither across chips, which can introduce a significant rotational offset between exposures (as seen in the TweakReg shift files below). In this notebook we focus on photometric differences between the WF2 and WF4 chips, so simply align the exposures using bright stars in common between the two flt images.

updatewcs('u*flt.fits', use_db=False)
clear_output()

4. Align the Images

Table of Contents

After resetting the WCS, we check for any residual pointing errors between the two dithered exposures. Precise aligment is needed to achieve the best drizzle-combined products. The expected pointing accuracy for various observing scenerios is summarized in the DrizzlePac Handbook Chapter 4.4. 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 in this repository. Note we use a farily large searchrad of 1 arcsec, because we have removed the Gaia WCS and are using the original ‘distortion only’ solution which can have large errors in the case of large dithers.

First we align the flt files and inspect the results before aligning the c0m files and verifying that output shiftfile values are identical. While testing the alignment, we recommend setting updatehdr = False to confirm that the alignment results look accurate and then rerunning with the keyword set to True.

input_images = sorted(glob.glob('u*flt.fits'))

tweakreg.TweakReg(input_images, 
                  updatehdr=True, 
                  clean=True,
                  reusename=True, 
                  interactive=False,
                  conv_width=3.0, 
                  threshold=200.0,
                  ylimit=1,
                  shiftfile=True, 
                  outshifts='shift_flt.txt',
                  searchrad=1,
                  tolerance=3,
                  minobj=7)
clear_output()

If the alignment is unsuccessful, stop the notebook

with open('shift_flt.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

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]

Note that there is a small shift in x=1.8 pix and y=-0.4 pix and large rotation of 0.014 degrees between the two images, likely due to the large chip-sized dither. These residuals are larger than the RMS scatter of 0.15 pixels and are therefore likely to be real.

shift_flt_table = Table.read('shift_flt.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_flt_table.colnames[1:]):
    shift_flt_table[col].format = formats[i]
shift_flt_table

Table length=2

file	dx	dy	rot	scale	xrms	yrms
str18	float64	float64	float64	float64	float64	float64
u5jx0108r_flt.fits	0.00	0.00	0.000	1.00000	0.00	0.00
u5jx010hr_flt.fits	1.81	-0.40	0.014	1.00001	0.15	0.15

Inspect the Astrometric residual scatter plots

Good practice includes inspecting not only the shift file residuals, but also the accompanying diagnostic plots. The residuals in the plot below are clustered around 0.0 with no obvious systematics.

Image(filename='residuals_u5jx010hr_flt.png', width=500, height=300)

Here we align the c0m files with the same parameter settings:

input_images = sorted(glob.glob('u*c0m.fits'))

tweakreg.TweakReg(input_images, 
                  updatehdr=True, 
                  clean=True,
                  reusename=True, 
                  interactive=False,
                  conv_width=3.0, 
                  threshold=200.0,
                  ylimit=1,
                  shiftfile=True, 
                  outshifts='shift_c0m.txt',
                  searchrad=1,
                  tolerance=3,
                  minobj=7)

clear_output()

Finally, we compare the shift file computed from the flt and c0m files to confirm that the results are identical.

shift_c0m_table = Table.read('shift_c0m.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_c0m_table.colnames[1:]):
    shift_c0m_table[col].format = formats[i]
    
shift_c0m_table

Table length=2

file	dx	dy	rot	scale	xrms	yrms
str18	float64	float64	float64	float64	float64	float64
u5jx0108r_c0m.fits	0.00	0.00	0.000	1.00000	0.00	0.00
u5jx010hr_c0m.fits	1.81	-0.40	0.014	1.00001	0.15	0.15

shift_flt_table

Table length=2

file	dx	dy	rot	scale	xrms	yrms
str18	float64	float64	float64	float64	float64	float64
u5jx0108r_flt.fits	0.00	0.00	0.000	1.00000	0.00	0.00
u5jx010hr_flt.fits	1.81	-0.40	0.014	1.00001	0.15	0.15

5. Drizzle-combine the images

Table of Contents

All four chips are now drizzled together. For speed in the notebook, we set the output pixel scale to that of the WF chips, but this can be reset to the PC scale if desired, to match the scale of the new drw drizzled products.

While we have turned sky subtraction off in the notebook, recommend using skymethod = ‘localmin’ and not ‘match’ for WFPC2 to avoid issues with overlapping vignetted regions along chip boundaries. See Figure 3.1.2 in the WFPC2 Instrument Handbook. For more information drizzling HST data, see the Drizzlepac Handbook and the list of commands in the readthedocs.

Users may wish turn on the sky subtraction and cosmic-ray rejection steps by setting skysub, median, blot, and driz_cr parameters to True. In that case, the recommended parameter settings for WFPC2 are driz_cr_snr='5.5 3.5', and driz_cr_scale='2.0 1.5'. For more discussion on drizzling parameters for WFPC2 when using the output scale of the WF chips, see the Prior Drizzling Example.

astrodrizzle.AstroDrizzle('u*flt.fits', 
                          output='wfpc2_flt', 
                          preserve=False, 
                          build=False, 
                          context=False, 
                          skysub=False, 
                          driz_sep_wcs=True, 
                          driz_sep_scale=0.0996, 
                          driz_sep_bits='8,1024', 
                          driz_sep_fillval=-1, 
                          median=False, 
                          blot=False,
                          driz_cr=False, 
                          final_fillval=None, 
                          final_bits='8,1024', 
                          final_wcs=True, 
                          final_scale=0.0996)
clear_output()

astrodrizzle.AstroDrizzle('u*c0m.fits', 
                          output='wfpc2_c0m', 
                          preserve=False, 
                          build=False, 
                          context=False, 
                          static=True, 
                          skysub=False, 
                          driz_sep_wcs=True, 
                          driz_sep_refimage='u5jx010hr_single_sci.fits', 
                          driz_sep_bits='8,1024', 
                          driz_sep_fillval=-1, 
                          median=False, 
                          blot=False, 
                          driz_cr=False, 
                          final_fillval=None, 
                          final_bits='8,1024', 
                          final_wcs=True, 
                          final_scale=0.0996, 
                          final_refimage='wfpc2_flt_drw_sci.fits')
clear_output()

Clean up unnecessary mask files.

fits_files = glob.glob('*ask.fits')
for f in fits_files:
    base_name = os.path.basename(f)
    if os.path.isfile(base_name):
        os.remove(base_name)

Inspect the single drizzled output frames derived from the dithered flt inputs

Here we compare the single drizzled frames. Note that WF4 and WF2 overlap in the footprint on the sky.

drz1_single = fits.getdata('u5jx010hr_single_sci.fits')
drz2_single = fits.getdata('u5jx0108r_single_sci.fits')

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(20, 10))
im1 = ax1.imshow(drz1_single, cmap='gray', vmin=0, vmax=4, origin='lower')
ax1.set_title('Target on WF4')
im2 = ax2.imshow(drz2_single, cmap='gray', vmin=0, vmax=4, origin='lower')
ax2.set_title('Target on WF2')

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')
plt.show()

Inspect the single drizzled output frames derived from the dithered c0m inputs

drz1_c0m_single = fits.getdata('u5jx010hr_c0m_single_sci.fits')
drz2_c0m_single = fits.getdata('u5jx0108r_c0m_single_sci.fits')

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(20, 10))
im1 = ax1.imshow(drz1_single, cmap='gray', vmin=0, vmax=4, origin='lower')
ax1.set_title('Target on WF4')
im2 = ax2.imshow(drz2_single, cmap='gray', vmin=0, vmax=4, origin='lower')
ax2.set_title('Target on WF2')

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')
plt.show()

Inspect the ratio of the two drizzled images to see changes due to photometric normalization

ratio_flt_c0m = drz1_single/drz1_c0m_single
fig = plt.figure(figsize=(7, 7))
img_ratio = plt.imshow(ratio_flt_c0m, vmin=0.97, vmax=1.03, cmap='Greys_r', origin='lower')
plt.colorbar(img_ratio, orientation='vertical')
plt.show()

Compare the combined drizzled image and the weight image.

flt_drw_sci = fits.getdata('wfpc2_flt_drw_sci.fits')
flt_drw_wht = fits.getdata('wfpc2_flt_drw_wht.fits')
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(20, 10))
im1 = ax1.imshow(flt_drw_sci, cmap='gray', vmin=0, vmax=0.5, origin='lower')
ax1.set_title('SCI')
im2 = ax2.imshow(flt_drw_wht, cmap='gray', vmin=0, vmax=50, origin='lower')
ax2.set_title('WHT')
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')
plt.show()

Note the Planetary Chip (PC) has ~4x WHT because the native pixel scale is ~2x smaller.

scale_ratio = 0.0996/0.0455
area_ratio = scale_ratio**2 
print('')
print('PIXEL SCALE RATIO: WF/PC = (0.0996/0.0455) =    ', scale_ratio)
print('PIXEL AREA  RATIO: WF/PC = (0.0996/0.0455)**2 = ', area_ratio)

PIXEL SCALE RATIO: WF/PC = (0.0996/0.0455) =     2.189010989010989
PIXEL AREA  RATIO: WF/PC = (0.0996/0.0455)**2 =  4.7917691100108675

Next, we compare the two drizzled images, zooming in to highlight where the two chips overlap. Note that the one on the right has a slightly more uniform background across the field of view. Finally, we then identify three isolated stars for photometry.

c0m_drw_sci = fits.getdata('wfpc2_c0m_drw_sci.fits')
flt_drw_sci = fits.getdata('wfpc2_flt_drw_sci.fits')

im1 = c0m_drw_sci[700:1600, 700:1600]
im2 = flt_drw_sci[700:1600, 700:1600]

fig, ax = plt.subplots(1, 2, figsize=(15, 10))
norm1 = ImageNormalize(im1, vmin=0.1, vmax=8, stretch=LogStretch())
ax[0].imshow(im1, norm=norm1, cmap='gray', origin='lower')
ax[1].imshow(im2, norm=norm1, cmap='gray', origin='lower')
ax[0].set_title('C0M drz_sci (zoom)', fontsize=20)
ax[1].set_title('FLT drz_sci (zoom)', fontsize=20)
plt.show()

Plot the ratio of the combined drizzled images from FLT and C0M

ratio_drw = flt_drw_sci/c0m_drw_sci
fig = plt.figure(figsize=(7, 7))
img_ratio = plt.imshow(ratio_drw, vmin=0.97, vmax=1.03, cmap='Greys_r', origin='lower')
plt.colorbar(img_ratio, orientation='vertical')
plt.show()

Perform aperture photometry on three stars in the WF2, WF4 overlap region

Here we perform aperture photometry in a radius of 5 pixels to show the change in count rates in each drizzle-combined image.

positions = [(1180, 1381), (1388, 1190), (1232, 1210)]    # Coords are (Y+1, X+1)
aperture = CircularAperture(positions, r=5)
annulus_aperture = CircularAnnulus(positions, r_in=10, r_out=15)

print('Aperture radius:', aperture)
print('Annulus:', annulus_aperture)

Aperture radius: Aperture: CircularAperture
positions: [[1180., 1381.],
            [1388., 1190.],
            [1232., 1210.]]
r: 5.0
Annulus: Aperture: CircularAnnulus
positions: [[1180., 1381.],
            [1388., 1190.],
            [1232., 1210.]]
r_in: 10.0
r_out: 15.0

data = fits.getdata('wfpc2_flt_drw_sci.fits')
aperstats = ApertureStats(data, annulus_aperture)
bkg_mean = aperstats.mean

phot_table = aperture_photometry(data, aperture)

for col in phot_table.colnames:
    phot_table[col].info.format = '%.8g'  # for consistent table output

print('drz_flt_combined: "wfpc2_flt_drw_sci.fits"')
print('     SKY values :', bkg_mean)
print('---------------------------------')
print(phot_table)

drz_flt_combined: "wfpc2_flt_drw_sci.fits"
     SKY values : [0.31061331 0.25361143 0.25064782]
---------------------------------
 id xcenter ycenter aperture_sum
--- ------- ------- ------------
  1    1180    1381    1139.9086
  2    1388    1190    484.00455
  3    1232    1210    210.30805

data = fits.getdata('wfpc2_c0m_drw_sci.fits')
aperstats = ApertureStats(data, annulus_aperture)
bkg_mean = aperstats.mean

phot_table = aperture_photometry(data, aperture)

for col in phot_table.colnames:
    phot_table[col].info.format = '%.8g'  # for consistent table output

print('drz_c0m_combined: "wfpc2_c0m_drw_sci.fits"')
print('     SKY values :', bkg_mean)
print('---------------------------------')
print(phot_table)

drz_c0m_combined: "wfpc2_c0m_drw_sci.fits"
     SKY values : [0.31046403 0.25329725 0.25064313]
---------------------------------
 id xcenter ycenter aperture_sum
--- ------- ------- ------------
  1    1180    1381    1136.0271
  2    1388    1190    482.57234
  3    1232    1210    209.71238

While the total change in photometry is only a few DN/sec, properly accounting for the relative gain across chips is necessary to achieve the best possible photometric precision. As we’ve demonstrated in this notebook, this is accomplished by using the new flt style files, or running photeq on the older c0m style files.

Conclusions

Table of Contents

MAST now include new WFPC2 flt.fits data products to replace the old c0m.fits, c1m.fits files. The new products include a new flux normalization between chips so that only a single zeropoint (inverse sensitivity, PHOTFLAM value) is needed for all chips. This allows for different chips to be drizzled together, for example, when combining images at multiple orientations.

About this Notebook

Table of Contents

---

Created: 14 Dec 2018;     M. Cara
Updated: 16 Nov 2023;     K. Huynh & J. Mack
Updated: 10 Dec 2024;     J. Mack & M. Revalski

Source: spacetelescope/hst_notebooks

Additional Resources

Below are some additional resources that may be helpful. Please send any questions to the HST Help Desk.

• WFPC2 Website

• WFPC2 2002 Data Handbook

• WFPC2 2008 Instrument Handbook

• Drizzlepac Handbook

• WFPC2 Prior Drizzling Example

Software Citations

If you use Python packages such as astropy, astroquery, drizzlepac, matplotlib, or numpy for published research, please cite the authors. Follow the links below for more information about citing various packages:

• Citing astropy

• Citing astroquery

• Citing drizzlepac

• Citing matplotlib

• Citing numpy

• Citing photutils

---

Top of Page