Analyzing HST imaging data with GIM2D

Introduction

Analyzing HST images with GIM2D is really easy. The following is a detailed step-by-step procedure to go from calibrated HST images to a structural parameter catalog. The use of GIM2D is by no means restricted to HST images. I ran it without any problems on ground-based data. In fact, GIM2D runs a LOT faster on ground-based images because it does not have to deal with pixel undersampling such as WFPC2's. Almost everything that follows can be directly applied to the analysis of ground-based data. If you have never worked with HST/WFPC2 data, you are strongly encouraged to look through the WFPC2 Instrument Handbook and the WFPC2 Data Handbook.

Index

Procedure In a Nutshell

Names in italics must be specified by the user.
  1. Obtain data from archive (!!)
  2. Use the task CRREJ in stsdas.hst_calib.wfpc to remove cosmic rays and combine images in the stack to get the science image.
  3. Run SExtractor on the science image and then SExreduce the output catalog to produce a .gfxt file (X_IMAGE Y_IMAGE BACKGROUND BKGSIG ISOAREA_IMAGE)
  4. Run GSCRIPTER to produce a .xgal file and a GIMFIT2D control IRAF script reduce.cl
  5. Combine Data Quality Files of stack using DQFCOMB
  6. Apply combined DQF image to the SExtractor segmentation image using GBADPIX
  7. Extract postage-size images from the science and segmentation images using XGAL
  8. Create a PSF for every object that will be analyzed. Run tiny1 to create a TinyTim parameter file psf.par. Run tiny2 to generate the PSFs.
  9. Run GIMFIT2D to perform bulge/disk decompositions
  10. Run GRESIDUAL to create a galaxy-subtracted image

Obtaining HST imaging data

HST images can be obtained either from the Space Telescope Science Institute (STScI) Archive or from the Canadian Astronomy Data Center (CADC). I prefer the latter. The CADC was the first to provide on-the-fly recalibration of HST data, and their WWW interface makes more sense. As an added bonus, you are able to cross-correlate your HST search with ground-based archives (CFHT, NTT, UKIRT, La Palma). This is no small bonus! CADC also offers WFPC2 image associations. Associations are groups of images taken at the same positions on the sky by one or many different observers, and the CADC pipeline groups them for you. CADC people are also a very friendly and helpful bunch.

Archive data comes in two parts: the science frames (.c0 files) and the data quality files (DQF, the .c1 files). Each science frame comes with a DQF frame. This frame flags pixels in the science frame which have been corrupted for one reason or another. The flag codes are given in the table below:

Flag ValueDescription
0Good pixel
1Reed-Solomon decoding error
2Calibration file defect
4Permanent camera defect
8A/D converter saturation
16Missing data
32Bad pixel

The DQF files are used in the GIM2D analysis to exclude bad pixels from the fits. Click here if you want more information on DQF files.

To convert the archive FITS images to GEIS format (.hhh), make sure that the IRAF variable imtype is set to hhh by typing "reset imtype=hhh" in the IRAF CL. You must use the STSDAS task STRFITS to convert your images. The HST archive images are three-dimensional stacks. The first slice is the PC image and the last three slices are the WF images. It is a good idea to create a directory for each WFPC2 chip and analyze each one separately. You can use the IRAF task IMCOPY to break down the image stacks.

IMPORTANT NOTE:You should really trim your HST archive images to remove the bad strip around the edges of those images to avoid problems later.

Cosmic ray rejection

There is a very nice cosmic ray rejection task in STSDAS called CRREJ. It is tailored for HST images, and it can be found in the HST_CALIB subpackage. It combines a stack of consecutive exposures while eliminating cosmic rays and scaling the remaining pixels to the total exposure time of the stack. Here is a good default set of CRREJ parameters: 

PACKAGE = wfpc
   TASK = crrej
    
input   =         @f606_3.list  Input images
output  =               f606_3  Output image name
masks   =                       Output masks (optional)
sigmas  =                  6,4  Rejections levels in each iteration
(radius =                   0.) CR expansion radius in pixels
(pfactor=                   0.) CR expansion discriminator reduction factor
(initial=                  min) Initial value estimate scheme (min or med)
(lower  =                 -99.) Lower limit of usable data
(upper  =                4096.) Upper limit of usable data
(sky    =                 mode) How to compute the sky
(expname=              exptime) Exposure time keyword name
(readnoi=                 0.76) Read noise(s) (in DN) of the detector(s)
(atodgai=                    7) Detector gain(s) in electrons/DN
(scaleno=                    5) Multiplicative term (in PERCENT) from noise mode
(dq     =                     ) Data quality filter pset
(skyname=                     ) Name of the group parameter to be updated with t
(crdqval=                  128) Data quality value for pixels flagged as cosmic
(fillval=               -9999.) Fill value for pixels having no good values in a
(verbose=                  yes) Print out verbose messages?
(mode   =                   al)
The value of the parameter sigmas is very important. If sigmas is too high, you will miss too many cosmic rays. On the other hand, if sigmas is too low, you will start eating away at the noise in your background pixels. Unfortunately, there is no clear-cut way of determining sigmas. Here is what I do. I operate on my image stack with sigmas set to ridiculously high values such that no pixels are rejected by CRREJ. Let us call this image HIGH. I then use CRREJ again on the stack with sigmas around "10,8" and I blink the resulting image against HIGH. Now, if sigmas is too low, you will modified pixel values in background areas where no cosmic rays hit. Increase sigmas if this is the case. If not, decrease sigmas. Iterate as many times as necessary. What about low-energy cosmic rays you say? Well, once I have gone through the iterative rejection process, I use the GIM2D task GBKG to study the pixel intensity histogram of the background. If low-energy cosmic rays are important, you should see a significant deviation i.e. a tail from a pure Gaussian pixel intensity distribution.

Running SExtractor

There are at least two galaxy photometry packages out there: FOCAS and SExtractor. Although I initially experimented with FOCAS, I opted for SExtractor because it produces what is called a "segmentation" image which is extremely useful to deblend galaxy images. Every pixel that SExtractor judges belong to a given object are assigned the same pixel flag value in the segmentation image. Background pixels all have the same pixel flag value of zero. It is then trivial for GIM2D to concentrate on the pixels belonging to a given object and the background pixels surrounding it.

You can obtain SExtractor via anonymous FTP at ftp.iap.fr where you will find a tar file containing the source code, a User's Manual and the Astronomy and Astrophysics article describing SExtractor. SExtractor is basically used as a structural preprocessor to detect objects, create a segmentation image and determine the following five quantities for each object: X and Y coordinates (X_IMAGE, Y_IMAGE), local sky background (BACKGROUND), local sky sigma (BKGSIG) and isophotal area (ISOAREA_IMAGE). Note that BKGSIG is not among the standard set of parameters produced by SExtractor. You must add it in yourself following the procedure given in the User Manual or contact Luc.Simard(at)nrc-cnrc.gc.ca. IMPORTANT NOTE: This is not required if you set the GIMFIT2D/C2GIMFIT2D/MGIMFIT2D parameter "dobkg" to yes. If dobkg is set to yes, GIM2D will automatically re-calculate the local sky level and sigma, and the background sigma information in the .gfxt file (below) will be ignored.

SExtractor runs on FITS images, and the command is simply: 
sex image.fits -c image.sex
where image is the name of the image to be analyzed, and image.sex is the name of the SExtractor configuration to be used. Here is a an example of such a configuration file. The parameter DEBLEND_MINCONT is particularly important as it specifies the fraction of the total object flux a "branch" must have in order for this branch to be considered as a separate object. The lower this parameter, the more eager SExtractor will be to separate a given object into multiple sub-components. You will also need a SExtractor parameter file such as this one. Your SExtractor catalog should now list X_IMAGE, Y_IMAGE, BACKGROUND, BKGSIG, and ISOAREA_IMAGE for each object. If the catalog file has a .cat suffix, you can now change that suffix to .gfxt.

Going back to IRAF/GIM2D ... the .gfxt file is used as input for the GIM2D task GSCRIPTER. GSCRIPTER generates a file used by XGAL to extract postage-size images around each object from the science and segmentation images. GSCRIPTER also generates an IRAF script that can be used to run GIMFIT2D in batch mode. See the GSCRIPTER page for more information.

Data Quality Files

Data Quality Files (DQF) are used by GIM2D to reject bad pixels from bulge/disk decompositions. Rejecting bad pixels in GIM2D requires two steps. First, you must combine the DQF files associated with your image stack with the task DQFCOMB. This task identifies pixels in the stack that have NO good pixel values and flags them with the value -2.0. The DQFCOMB command is of the form: 
dqfcomb image imagedqf @dqf.list
where dqf.list is the name of the file containing a list of the DQF images associated with your image stack.

The second step consists of applying the combined DQF mask image to the SExtractor segmentation image with the task GBADPIX: 
gbadpix segimage imagedqf
Note that the badflag should always be set to -2 since this is the bad pixel flag recognized by GIMFIT2D.

Extracting galaxy images

Extract a postage-size image around each object from both the science image and the segmentation image to be analyzed by GIMFIT2D is done in one simple command: 
xgal image segimage objects.xgal
where image is your science image, segimage is the segmentation image associated with your science image, and objects.xgal is the name of the file produced by GSCRIPTER with objects (X,Y) and postage image size. XGAL extracts a postage-size image around each object in objects.xgal and subtracts the local sky value measured by SExtractor from it.

TINYTIM PSF modelling

The TinyTim software written by J. Krist can be used to generate WFPC2 Point-Spread-function at arbitrary positions and in any of the available filters. I usually generate a PSF at the location of every object I plan to analyze with GIM2D. Tinytim can be used in batch mode to do so, but it is important to make sure that the parameter MAX_POSITIONS in the include file tinytim.h is larger than the total number of PSFs you plan to create in one batch. Increase MAX_POSITIONS as necessary before you compile TinyTim on your system. The object x and y coordinates can be extracted from the .xgal file using something like the UNIX awk command.

IMPORTANT: I also recommend generating oversampled PSFs. Wide-Field camera TinyTim PSFs should be oversampled by a factor of 5 whereas Planetary Camera should be oversampled by a factor of 3.

Jitter can also be important, and, unfortunately, there is no information in HST image headers about the amount of jitter present in a given image. To determine the amount of jitter in your image (if it is a major concern), I suggest the following procedure although I have not tried it myself yet: One runs Tinytim in two stages. Tiny1 queries the user to make a TinyTim parameter file (click here for an example), and Tiny2 actually creates the PSFs using the tiny1 parameter file.

IMPORTANT TIP: If you use an oversampling factor of 5 for WF TinyTim PSFs, use a PSF diameter of 2.4 arcsec so that the PSF image will be 120X120. The idea is to keep the PSF image size below 128X128 so that GIMFIT2D does not have to step up to the next power of 2 (256X256) to keep various Fourier Transforms happy.

Running GIMFIT2D

This step could not be easier if you have correctly perform all the previous steps. The command is simply: 
cl < reduce.cl >& reduce_machinename.log &
where reduce.cl is the IRAF script produce by GSCRIPTER and reduce_machinename.log which tells us which images have been processed so far and when. A nice feature of GIMFIT2D is that you can issue the above command with the same reduce.cl on as many different machines as you want, and these machines will all go down the same master list of galaxy images to analyze without stepping on each other's toes. This means you can use multiple CPUs to crunch the same frame simultaneously by having each CPU work on a different galaxy.

You can use GRESIDUAL to generate a galaxy-subtracted version of your input science image.

Last Update: 29 October 2001
Comments/questions to Luc.Simard(at)nrc-cnrc.gc.ca