Canadian Astronomy Data Centre
Herzberg Institute of Astrophysics

ST-ECF
ESA-ESO

The WFPC2 Association Project Version 2 Stacking algorithm

Association Example: a tale of the Chandra Deep Field
Single Exposure	Association Product

Please note that this page is in construction. More details will appear very soon

Pipeline details

The wfpc2 associations pipeline is divided in multiple tasks. They are:

IMSTACK: Stacking the members of an association together

IMSTACK VERSION: 4.1

Introduction

The imstack process in 6 sentences:

Rescales all the input images to common gain, exposure time and common background sky level;
Shift all members to the same sky position, with the provision that imshift is not applied if the shift is lower than a given tolerance (set to be 0.15 pixels);
Runs the Skepticism Algorithm (P.Stetson, 1979) a first time to estimate a better weight map for a subsequent second run;
Runs the Skepticism Algorithm a second time to stack together the imshifted images, so to obtain:
1. the FINAL stack, normalised at the AVERAGE of the exposure times
2. the FINAL weight map
Computes, and stores in the header, all the Characterisation VO parameters, including coverage, sampling and resolution of the spatial, spectral and temporal axes;
Stores in the header all the used zeropoints (ZEROSTS,ZEROSTE,ZEROABS,ZEROABE).

The full description

The entire pystack.py is here described. It consists of the following sections:

Initialisation;
Prepare the weight maps
1. Prepare a mask using the data quality (c1 extension) and the flat field information;
2. Prepare a variance map for each of the input (member) images, taking into account:
  - photon noise
  - readout noise
  - dark contribution
  - flat fielding
3. Obtain a weight map (of type variance) by multiplying the mask and the variance map;
Scale the input science images to same gain, exposure time, and same sky level;
Shift both science and weight images to common grid;
Run the skepticism algorithm (pystack) a first time, to obtain a better guess of what each input image should have looked like;
Scale back the obtained pystack science image to each original member level (gain, exptime, and sky level);
Compute the better weight map for each member image, using the image obtained in the previous step, and reapplying the same method;
Run the skepticism algorithm (pystack) again, with the new weight maps;
Compute the Characterisation VO parameters, including coverage, sampling and resolution of the spatial, spectral and temporal axes (e.g. CVOSRES, CVOWSPA, etc.).

INITIALISATION
Some initialisation, of which worth noting are:
- setting readnoise and gain for the four chips (see Table GAIN and READNOISE below)
- setting the maximum shift tolerance to 0.15 pixels
  members shifted by less than shift_tolerance pixels are considered well aligned, and no imshift will occur.
- computing the min and max of the x and y shifts among the members (from the .asn file)
- checking whether the members have all the same gain or not (quot;mixed gain")
- determining the nominal gain of the coadded output image (stack_image),to be:
  
  nominal gain = 14.0 if all members have ATODGAIN=15
  
  nominal gain = 7.0 in all other cases(all ATODGAIN=7, or mixed cases of 7 and 15).
- computing the average exposure time

PREPARE THE WEIGHT MAPS

For each member, and for each chip:

Prepare a mask to flag artifacts

using the data quality (c1.fits extension) and the flat field information,
(1 - good, 0 - bad):

find which (inverse) flatfield was used from the member c0.fits file
build a flatfield mask (MASK_{FF_inv}) from the original flatfield:

set to val=1
all pixels which exhibit an inverse flatfield lower than 2.

set to val=0
all bad pixels which exhibit an inverse flatfield greater than 2,
to cope with the bad vignetting.

build a flag mask from the c1.fits (MASK_C1):

set to val=1: all the good pixels (where c1 pixel val is 0-good or 1024-corrected warmpixel)
set to val=0: all the bad pixels (where c1 pixel value is not 0 AND not 1024).; For an explanation of the c1 pixel values, please refer to the WFPC2 Data Quality table.

The mask that will be used will hence be:

MASK_member,chip = MASK_{FF_inv} * MASK_C1 (1)

Prepare a variance map
for each of the input (member) images;
1. Variance of a C0 member file:
  1. Given an image in electrons, its variance can be computed using the following formula:
    
    (2)
    
    whereby the three addictive terms cover the noise contribution respectively of the poissonian photon noise, the read out noise, and the noise contribution of the dark current:
    
    C0_e^-
    stands for the actual member image values in electrons
    
    g_e^-/adu
    is the gain setting of the chip in electrons per adu
    
    RN
    the readout noise of the given chip
    
    P_dark(T)
    is the number of electrons expected per pixel and per second,
    which depends on the temperature (T) of the detector (see below)
    
    t + 46
    is the total time during which the dark current accumulated,
    as per the WFPC2 Instrument Handbook Exposure Time Calculator paragraph.
  2. The science image is usually expressed in counts (adu) and not in electrons:
    
    (3)
    
    where
    
    C0_adu
    stands for the actual member image pixel values in adu
    
    FF_inv
    stands for the inversed flat field pixel values
  3. Hence the variance image transforms similarly:
    
    (4)
  4. Applying both (3) and (4) to the equation (2) we obtain the final equation for the variance map in counts (adu):
    
    (5)
Obtain a weigth map (of type variance) by multiplying the mask (1) and the variance map (5).
Hence, the weight map which will be used for each member is expressed by the formula:

WEIGHT_MAP_member,chip = MASK_member,chip * VAR_{adu_{FF_inv}} (6)

Masking the weight map (setting it to 1.E20) and the mask (to zero) for the vignetted region of the chip

MASK_member,chip = 0 if ( x < vignetted or y < vignetted )
WEIGHT_MAP_member,chip = 1.E20 if ( x < vignetted or y < vignetted )
(7)

where the vignetted limits in X and Y vary with the chip number:

From Table 1.2, p.10, WFPC2 Data Handbook v4.0, Jan 2002
Chip X limit Y limit

1 44 52

2 46 26

3 30 47

4 43 44

Scale the input science images
to same gain, exposure time, and same sky level
Each member (subscript j) science data is rescaled to the same nominal gain ( gain_stack ), and to the average exposure time ( avg_exptime ) of the association:

Scaled_Image_j = Image_j * (gain_j / gain_stack) (avg_exptime/ exptime_j)
(9)
Shift science, weight, flat and mask images
Both the science image, and accompanying weight map, of each chip are shifted onto a common grid.
At this stage we also shift the flat field and the mask images, since they will be needed for a second iteration.
We use the imshift iraf task with linear interpolation for shift values > 0.15 pixels. That is, a lower shift is considered to be negligible.
XXX PLEASE, TOMORROW CONTINUE FROM HERE XXX

Stacking the images: iteration 1

Data_js and Variance_js
where Data_js = Data_j shifted

Then we need to shift the MAP_j and the Flat_j into MAP_js and Flat_js

Then we process these 2 input images lists with pystack and the result is
a stack image with bad values masked to 10e06 called P₁

Processing Steps for WFPC2 V2: Pass two (for J images)

It is essentially the same equation than above but the have first to "descale" the P₁ image with the ration of gain and the ration of exposure time, i.e.: (gain_j / gain_stack) and (t_m/ t_j)

The variance J2

Variance_j2 = [ (Flat_js * P₁ / gain_j ) * (t_j/t_m) * (gain_stack / gain_j ) + Flat_js² * R_j²/ gain_j² ] * (gain_j / gain_stack) * (t_m/ t_j)* MAP_js

The data J2 is simply

Data_j2 = Data_js

Then pystack pass 2

The resulting stacked image and stack are composing the final product.

Astrometric correction (using USNOB)

Image caracterisation

catalogue extraction

Storage in AD

Readout Noise Table

Readout noises
used by the `imstack.py` procedure
Nominal Gain	Readout Noise (electrons)
Nominal Gain	chip 1	chip 2	chip 3	chip 4
7	5.24	5.51	5.22	5.19
15	7.02	7.84	6.99	8.32

Gain Table

Real Gain
used by the `imstack.py` procedure
Nominal Gain	Real Gain
Nominal Gain	chip 1	chip 2	chip 3	chip 4
7	7.12	7.12	6.90	7.10
15	13.99	14.50	13.95	13.95

GAIN manipulation

From Source extractor for dummies V4

Method Effective Gain Magnitude zeropoint Input image

1 gain * total_exp_time Zeropoint(1 sec) Image in count/sec

2 gain Zeropoint(1 sec) * log₁₀(total_exp_time) Sum of N frames

3 N * gain Zeropoint(1 sec) * log₁₀(average_exp_time) Average of N frames

4 2 * N * gain / 3 Zeropoint(1 sec) * log₁₀(average_exp_time) Median of N frames

WFPC2 Data Quality Flag Values

AS from the table 3.4 in WFPC2 Instrument Handbook


Flag Value	Description
0	Good pixel
1	Reed-Solomon decoding error. This pixel is part of a packet of data in which one or more pixels may have been corrupted during transmission.
2	Calibration file defect-set if pixel flagged in any calibration file. Includes charge transfer traps identified in static pixel mask file (`.r0h`).
4	Permanent camera defect. Static defects are maintained in the CDBS database and flag problems such as blocked columns and dead pixels. (Not currently used.)
8	A/D converter saturation. The actual signal is unrecoverable but known to exceed the A/D full-scale signal (4095).¹
16	Missing data. The pixel was lost during readout or transmission. (Not currently used.)
32	Bad pixel that does not fall into above categories.
128	Permanent charge trap. (Not currently used.)
256	Questionable pixel. A pixel lying above a charge trap which may be affected by the trap.
512	Unrepaired warm pixel.
1024	Repaired warm pixel.

Estimate of Dark Current contribution to the noise

From WFPC2 Instrument book c14/15, page 88, the dark count rate in electrons per second per pixel is tabularised as follows:

Temperature Dark counts
(electrons/second/pixel)

-20 10

-30 3.

-40 1.

-50 0.3

-70 0.03

-77 0.016

-83 0.008

-88 0.0045

The detector temperature is provided by the FITS keywords: UCH1CJTM, UCH2CJTM, UCH3CJTM, UCH4CJTM for each of the four chips.

Interpolation is required, and for that we use the following formula which approximates quite well the above table:

P_dark = (273 + T) * 10^{^( 0.04719*T - 0.4592 )}

Magnitude systems

STMAG

ZEROPT_second = -2.5 * log₁₀(PHOTFLAM) + PHOTZPT
ZEROPT_exptime = -2.5 * log₁₀(PHOTFLAM) + PHOTZPT + 2.5 * log₁₀(exptime)

PHOTZPT = -21.1 (always)
m_stmag = -2.5 * log₁₀(counts/exptime) + ZEROPT_exptime
m_stmag = -2.5 * log₁₀(counts/sec) + ZEROPT_seconds
m_stmag = -2.5 * log₁₀( PHOTFLAM * DN / EXPTIME) + PHOTZPT

AB mag

m_AB = -2.5 log₁₀( Flux_nu ) - 48.594
where flux_nu is in ergs cm^-2 s^-1 HZ^-1
m_AB = -2.5 log₁₀( Flux_lambda ) - 2.402 - 5.0 * log₁₀( lambda )
where flux_lambda is in ergs cm^-2 s_-1 A_-1
and lambda in A
F_nu = f_lambda * lambda² / c

flam = 3x10^18 * fnu / lam^2
fnu = (lam^2 * flam) / 3x10^18

Converting an image to flux

flux = DN * PHOTFLAM / exposure_time
where flux is in erg cm^-2s^-1A^-1

flux = DN * PHOTFLAM / exposure_time * PHOTPLAM**2 / 2.9972E18 * 1E+23
where flux is in Jy

JANSKY

1 Jy = 1.0 * 10^-26 (W/m²/Hz) = 1.0E-23 erg/s/cm²/Hz

FWHM

        seeing_fwhm = 1.22 * (photplam / 24000000000.0 ) * 206264.806247
	seeing_fwhm = math.sqrt( (seeing_fwhm * seeing_fwhm) + (pixscale[chip] * pixscale[chip])  )
	24000000000.0 = size of the mirror
	206264.806247 = 180 * 3600. / pi


        $diameter=2.4e10; #diameter of the telescope

	$arc_rad = 206264.8; #number of arcsec/rad

	$rayleigh=1.22*$photplam * $arc_rad / $diameter; #1.22 * lambda/D

	$resolution  = sqrt( ($rayleigh**2) + ($pixscale**2))

CVORES

cvosres = (( photplam / 10.0E10 ) / 2.4 ) * (2.063 * (10**5)) cvonyqr = cvosres / pixscale[chip] cvosres = cvosres / 3600.0 The Nyquist ratio (spatial_resolution/spatial_sample). Values less than 2.5 are undersampled.

equations

# [Y ABnu] = -2.5 * log([X ergs/cm^2/s/Hz]) - 48.594
# [Y ABnu] = -2.5 * log([X ergs/cm^2/s/A]) - 2.402 - 5.0 * log(lambda A)
# [Y ABnu] = -2.5 * log([X W/m^2/Hz]) - 56.094
# [Y ABnu] = -2.5 * log([X photon/cm^2/s/A]) + 16.852 - 2.5 * log(lambda A)
# [Y ABnu] = -2.5 * log([X photon/cm^2/s/um]) + 16.852 - 2.5 * log(lambda um)
[Y Jy] = 1.0E+23 * [X erg/cm^2/s/Hz]

Daniel Durand, April 2004

Alberto Micol, September 2004

Photometric systems: zeropoints and formulae


Generic: m = -2.5 * log10( DN/EXPTIME ) + ZEROPOINT

ST sys:  m_st = -2.5 * log10( F_lambda ) - 21.10

AB sys:  m_AB = -2.5 * log10( F_nu ) - 48.60

where F_lambda is in erg/cm2/sec/A
  and F_nu        in erg/cm2/sec/Hz

F_lambda = countRate * PHOTFLAM

F_nu     = F_lambda * lambda² / c
                   = countRate * PHOTFLAM * PHOTPLAM² / c

where:

 - counteRate is the number of counts diveded by the exptime

 - PHOTFLAM is the flux of a source with constant flux per unit wavelength
            (in erg/cm2/sec/A) which produces a count rate of 1DN/sec.
   Note: Obviously a change of GAIN will have repercussions on PHOTFLAM

 - c must be given in A/sec, that is: 2.99792E+18 A/sec
     since PHOTPLAM is in A, and PHOTFLAM in erg/cm2/sec/A.

Hence,

 m_st = -2.5 * log10( countRate * PHOTFLAM ) - 21.10
     = -2.5 * log10( countRate ) -2.5*log10( PHOTFLAM ) - 21.10

     = -2.5 * log10( count ) -2.5*log10( PHOTFLAM ) - 21.10 + 2.5 * log10( EXPTIME )


 m_AB = -2.5 * log10( countRate * PHOTFLAM * PHOTPLAM² / c ) - 48.60
     = -2.5 * log10( countRate ) -2.5*log10(PHOTFLAM*PHOTPLAM²/c) - 48.60

     = -2.5 * log10( count ) -2.5*log10(PHOTFLAM*PHOTPLAM²/c) - 48.60 + 2.5*log10( EXPTIME )

Effective Gain

The effective gain is to be used within PIPE only when running SExtractor,
that is, when the noise of the image must be computed.


The effective gain is g for the noise of an single image.

The effective gain is g for the noise of an coadded image (pure sum).

The effective gain is N * g for the noise of an averaged image.

   That's because the poissonian noise can be computed only from the Total number of electrons,

   not from the average number of electrons; and
   the total number of electrons is N * g * the number of counts

    of an averaged image.


All the formulae to do with gain

FLUX COMPUTATION


For a source which has 1DN/sec in an image with gain g,
the flux is PHOTFLAM, where PHOTFLAM is to be scaled to g7 or g15.

Flux(erg/cm2/sec/A) = N_DN/exptime * PHOTFLAM(g=g7 or g15)

For an averaged image, whose members are images at g=g7,
the flux is computed this way:

Flux = Total N_e / g7 / Total Exptime * PHOTFLAM(g7)

     = N * g7 * M_p / g7 / Total Exptime * PHOTFLAM(g7)

     = N/Total Exptime * M_p * PHOTFLAM(g7)

     = M_p / AVGExpTime * PHOTFLAM(g7)


For a summed image, whose members are images at g=g7,
the flux is computed this way:

Flux = Total N_e / g7 / Total Exptime * PHOTFLAM(g7)

     = g7 * T_p / g7 / Total Exptime * PHOTFLAM(g7)

     =  T_p / Total Exptime * PHOTFLAM(g7)


SExtractor outputs always FLUXes as counts;
hence,
a SExtractor run onto an averaged image, will return the average counts;
a SExtractor run on a total image will return the total number of counts.

-->

The WFPC2 Association Project Version 2Stacking algorithm