Current as of May 6, 2002.
We use the task mscstack to combine the images for a field into one deep exposure. Note that at the present time we do not make any effort to match the DIQ of the input images, although some groups using MOSAIC have successfully done this. We do, however, take some pains to enable the deepest possible final image, free of defects like faint residual images from satellite trails in an individual image. This can be accomplished if all such defects are masked out of the individual input images. In order to include as many defects in the input masks as possible, we use software written by Frank Valdes, which is still in active development, to compare each of the individual images to an initial average image that is created using a harsh rejection criteria. The software used to make the comparison, part of the ACE package Frank Valdes is writing, also produces a new bad-pixel mask that can be used in our "final" stacking of the images.
To make the "first pass" combined image that will be used to identify defects in the individual exposures, we use mscstack with a low threshold for rejection. Note that with the release of Iraf 2.12 you have a choice of using a fits format for your masks, or sticking to the old style of .pl masks. In either case, make sure that your masktype is set appropriately: masktype = fits or masktype = pl. If you have 2.12 available, we recommend using .fits type so that it is easier to write masks to fits tapes.
PACKAGE = mscred TASK = mscstack input = @infile.list List of images to combine output = Stackharsh.fits Output image (headers = "") List of header files (optional) (bpmasks = "harsh_bpm.fits") List of bad pixel masks (optional) (rejmasks = "") List of rejection masks (optional) (nrejmasks = "") List of number rejected masks (optional) (expmasks = "") List of exposure masks (optional) (sigmas = "") List of sigma images (optional)\n (combine= average) Type of combine operation (median|average) (reject = ccdclip) Type of rejection (masktyp= goodvalue) Mask type (maskval= 0.) Mask value (blank = 50000.) Value if there are no pixels (scale = !mscscale) Image scaling (zero = !msczero) Image zero point offset (weight = none) Image weights (statsec= ) Image section for computing statistics (lthreshold= 1.) Lower threshold (hthreshold= INDEF) Upper threshold (nlow = 1) minmax: Number of low pixels to reject (nhigh = 1) minmax: Number of high pixels to reject (nkeep = 1) Minimum to keep (pos) or maximum to reject (neg) (mclip = yes) Use median in sigma clipping algorithms? (lsigma = 5.) Lower sigma clipping factor (hsigma = 2.9) Upper sigma clipping factor (rdnoise= rdnoise) ccdclip: CCD readout noise (electrons) (gain = gain) ccdclip: CCD gain (electrons/DN) (snoise = 0.) ccdclip: Sensitivity noise (fraction) (sigscal= 0.1) Tolerance for sigma clipping scaling corrections (pclip = -0.5) pclip: Percentile clipping parameter (grow = 0.) Radius (pixels) for neighbor rejection (mode =
The task,diffdetect, is run to identify differences between one of the individual exposures and the average or stacked image. The mask that was created by mscstack is just a record of the number of input exposures that did not contribute to each pixel in the output image. The value for each pixel in the mask is an integer from 0, if that given pixel was contributed to by every input frame, to n, where n is the number of images in the stack whose data were excluded from being used in the determination of the average. The task diffdetect knows which pixels have little or no information (like the positions of saturated stars, which in general will have been saturated in all of the input images). These positions will not be flagged further.
diffdetect compares each image with the stack and notices what is different taking into account the noise properties of both the single exposure and the combined stack. You may need to load the tasks xdwred and aceold to see diffdetect. An example parameter file.
PACKAGE = ace
TASK = diffdetect
images = mos068S List of images
(masks = BPM) List of bad pixel masks
(skys = SKYFIT) List of skys
(sigmas = SKYSIG) List of sky sigmas
(scales = !mscscale) List of image intensity scale factors
# Reference Image(s)
(rimages= Stackharsh.fits) List of reference images
(rmasks = harsh_bpm.fits) List of reference bad pixel masks
(rskys = SKYFIT) List of reference skys
(rsigmas= SKYSIG) List of reference sky sigmas
(rscales= 1.) List of reference intensity scale factors
# Output
objmasks= diff068S List of output object masks
catalogs= diff068S List of output catalogs
(logfile= STDOUT) List of log files
# Sky
(newsky = yes) Determine new sky fit if one already exists?
(nskylin= 100) Number of sky sample lines
(skyblk1= 10) Sky block size for 1D sky estimation
(skyhcli= 3.) High sky clipping during 1D sky estimation
(skylcli= 3.) Low sky clippling during 1D sky estimation
(skyxord= 4) Sky fitting x order
(skyyord= 4) Sky fitting y order
(skyxter= half) Sky fitting y order
# Detection
(convolv= block 3 3) Convolution kernel
(hsigma = 3.) Sigma threshold above sky
(lsigma = 10.) Sigma threshold below sky
(minpix = 10) Minimum number of pixels in detected objects
(sigavg = 6.) Sigma of mean flux cutoff
(sigmax = 6.) Sigma of maximum pixel
(rfrac = 0.5) Minimum fraction of reference flux
(bpval = 1) Output bad pixel value
# Growing
(ngrow = 2) Number of grow rings
(mode = ql)
This will generate two files, a list of objects in an STScI
tables format, and a bad pixel mask.
We use acedisplay and reviewproto to examine the image
and make sure that 1) no good areas of the image are included in the
new mask by some error, and 2) that everything we want in the mask has
been caught. First look at the stack and make sure there are no
surviving satellite trails or scattered light. If there are, you might
have to go through the individual images and mask these out by hand in
the new ace-generated bad pixel mask using satzap for satellite
trails and blkzap for bad areas in general (and tasks to "undo"
corrections, unsatzap and unblkzap. The source for
these scripts is available at these links: satzap.cl by James Rhoads and Arjun Dey
blkzap.cl by Arjun Dey, based on script by James
Rhoads
and (unblkzap.cl and (unsatzap.cl).
Here is the parameter file for acedisplay. Remember to make sure that your stdimage is set to imt8800 in your login.cl so that you have proper resolution for this task to work. Also note that the parameter "overlay" is pointing to a header keyword called objmask. This keyword gets added when you run diffdetect. If you want to use acedisplay before you run diffdetect set overlay to BPM so that it will overlay the original bad pixel mask.
PACKAGE = ace TASK = acedisplay image = mos068S image to be displayed frame = 1 frame to be written into (bpmask = BPM) bad pixel mask (bpdispl= overlay) bad pixel display (none|overlay|interpolate) (bpcolor= red) bad pixel colors (overlay= !objmask) overlay mask (ocolors= 1=red,+203) overlay colors (erase = yes) erase frame (border_= no) erase unfilled area of window (select_= yes) display frame being loaded (repeat = no) repeat previous display parameters (fill = no) scale image to fit display window (zscale = yes) display range of greylevels near median (contras= 0.25) contrast adjustment for zscale algorithm (zrange = yes) display full image intensity range (zmask = ) sample mask (nsample= 1000) maximum number of sample pixels to use (xcenter= 0.5) display window horizontal center (ycenter= 0.5) display window vertical center (xsize = 1.) display window horizontal size (ysize = 1.) display window vertical size (xmag = 1.) display window horizontal magnification (ymag = 1.) display window vertical magnification (order = 0) spatial interpolator order (0=replicate, 1=linear) (z1 = ) minimum greylevel to be displayed (z2 = ) maximum greylevel to be displayed (ztrans = linear) greylevel transformation (linear|log|none|user) (lutfile= ) file containing user defined look up table (mode = ql)
Another useful task is tinfo. If you run tinfo on the output diff*.tab that you got from diffdetect it will tell you how many rows were written to each table. Each row represents a portion of the mask. You want that number to be somewhere less than 200 (unless you're using the engineering grade ccd chips). If it's much more than 200 you should look at that image to see if some large area which should have been blocked out by a large mask was blocked out by many little ones, ie- a satellite trail or bad image edge.
The task, reviewproto, will display three images: the stacked image you made with mscstack, an individual image, and a masked area of that image. By blinking these three back and forth, or more easily, using the tile option on your imagetool, you can compare these three images. Look at the stacked image to see if the object is real or not, then look at the individual image and the mask. For example, if the object is not on the stacked image, but is on the individual image, chances are it is a cosmic ray, so it should be masked. If you find that something was masked when it should not have been, you will need to use imreplace on the mask file to un-mask that portion of the stacked image. You will need to display the mask file, and determine the image section of the offending portion of the mask and then replace that section with 0. The parameter files for reviewproto,imreplace, blkzap, and satzap are below.
PACKAGE = ace TASK = reviewproto catalog = diff051S Catalog (nooverl= yes) Display image without overlays (overlay= yes) Display image with overlays (compari= yes) Display comparison image (compima= Stackharsh.fits) Comparison image (box = 200) Box size (pixels) (ocolors= 204) Object mask color (lcolor = 205) Label color (fd = ) (mode = ql) PACKAGE = imutil TASK = imreplace images = diff051S.pl[7782:7828,950:994] Images to be edited value = 0. Replacement pixel value (imagina= 0.) Imaginary component for complex (lower = INDEF) Lower limit of replacement window (upper = INDEF) Upper limit of replacement window (radius = 0.) Replacement radius (mode = ql) PACKAGE = user TASK = blkzap object = mos051S Input image x1 = 1247. lower left corner: location 1 - X y1 = 3705. lower left corner: location 1 - Y x2 = 1851. upper right corner: location 2 - X y2 = 3787. upper right corner: location 2 - Y (mode = al) PACKAGE = user TASK = satzap object = mos051S Input image x1 = 1151. Satellite trail: location 1 - X y1 = 3968. Satellite trail: location 1 - Y x2 = 1596. Satellite trail: location 2 - X y2 = 4011. Satellite trail: location 2 - Y (hwidth = 10.) half-width of flagged track? (mode = al)
Once you're satisfied with the quality of your images, you've discarded trailed or problem images, and you're happy with your mask files, the final stack can be made! To get to this point may take several iterations of mscstack and diffdetect. The final stacked image is obtained by using mscstack again, but with a much less harsh rejection setting, counting on our improved bad-pixel masks to reject the contribution of unwanted pixels. Remember to use hedit on your images to update the BPM keyword with the new diff*.pl mask files you've just created with diffdetect. You will also have to create a weight file before you can run mscstack. That is explained below.
PACKAGE = mscred TASK = mscstack input = @infile.list List of images to combine output = 21q1_final.fits Output image (headers = "") List of header files (optional) (bpmasks = "21q1_final_bpm.fits") List of bad pixel masks (optional) (rejmasks = "21q1_final_rej.fits") List of rejection masks (optional) (nrejmasks = "21q1_final_nrej.fits") List of number rejected masks (optional) (expmasks = "21q1_final_exp.fits") List of exposure masks (optional) (sigmas = "") List of sigma images (optional)\n (combine= average) Type of combine operation (median|average) (reject = none) Type of rejection (masktyp= goodvalue) Mask type (maskval= 0.) Mask value (blank = 50000.) Value if there are no pixels (scale = !mscscale) Image scaling (zero = !msczero) Image zero point offset (weight = @21q1_weights) Image weights (statsec= ) Image section for computing statistics (lthresh= 1.) Lower threshold (hthresh= INDEF) Upper threshold (nlow = 1) minmax: Number of low pixels to reject (nhigh = 1) minmax: Number of high pixels to reject (nkeep = 1) Minimum to keep (pos) or maximum to reject (neg) (mclip = yes) Use median in sigma clipping algorithms? (lsigma = 5.) Lower sigma clipping factor (hsigma = 3.) Upper sigma clipping factor (rdnoise= rdnoise) ccdclip: CCD readout noise (electrons) (gain = gain) ccdclip: CCD gain (electrons/DN) (snoise = 0.) ccdclip: Sensitivity noise (fraction) (sigscal= 0.1) Tolerance for sigma clipping scaling corrections (pclip = -0.5) pclip: Percentile clipping parameter (grow = 0.) Radius (pixels) for neighbor rejection (mode = ql)
The weight parameter for mscstack is simply a list of numbers by which each image will be weighted. To obtain the proper weight for each image use the equation: one over the square of the product of mscscale and the rms of the background sky, ie- 1/(mscscale * (rms of the background sky))**2.
mscscale is simply a header keyword you can read with imhead or hedit. To obtain the rms of the background sky you can use iterstat or some similar task on a fairly large image section near the center of each image that is free of masks, and bright objects.
PACKAGE = user TASK = iterstat image = mos073S.fits[3871:3990,4304:4453] Input image(s) (nsigrej= 5.) Number of sigmas for limits (maxiter= 10) Maximum number of iterations (print = yes) Print final results? (verbose= no) Show results of iterations? (lower = INDEF) Initial lower limit for data range (upper = INDEF) Initial upper limit for data range (mean = ) Returned value of mean (sigma = ) Returned value of sigma (median = ) Returned value of sigma (valmode= ) Returned value of mode (inimgli= ) (mode = ql)You might want to create a little file to compute the weights by putting your mscscale header keyword values and results of iterstat into a file such as the one below, call it spud.
= 1/(1*(13.97367))**2 = 1 / (0.977049948652306*(13.52238))**2 = 1/(0.95787457643229*(13.31077))**2 = 1/(0.989892959525147*(11.89777))**2 = 1/(0.931024572149885*(11.73618))**2 = 1/ (0.926741617301131*(11.40967))**2 = 1/ (0.913667453628092*(10.6781))**2 = 1/ (0.924982168253485*(10.38573))**2Then you can redirect the above file to the cl and redirect the output to a weights file:
cl < spud > 21q1_weights
Your final output image from mscstack should be your final reduced product. Thus ends the "basic" reduction of the images. Next it is time to do the photometric calibrations and generate catalogs from these images.