[Creating Improved Bad-pixel Masks (mkbleed,xzap)]
[Using Bad-pixel masks (fixpix)]
[Tangent-Plane Projection (mscimage)]
If you obtain more than one exposure of a field in a given filter (say as part of a mosdither command with MOSAIC), you are likely to desire to optimally combine the multiple exposures into a single, deep, image free of some of the defects of the individual images (cosmic-ray hits, gaps between chips, etc.). Before you can combine the individual exposures, they must be projected on to a standard coordiante system and image scale. The next two sections of the notes are still less complete than the earlier sections, but here are some details of the steps we go through leading up to combining the images for the NDWFS MOSAIC images. The initial steps are done to remove defective pixels from the individual image that would be propogated into other pixels by the sinc function interpolation we use (to preserve the noise properties of the images) when we tangent-plane project our images prior to combining them.
So, in brief, in this and the next two
items under III., we will cover following steps:
1.) Create bad-pixel masks that identify cosmic rays, saturated pixels (this could now be done early in the reduction process, but we have not yet taken advantage of this capability of ccdproc), and other defects. This is done using the mkbleed4 script, the xzap task, and the combinemask script. If necessary use imreplace to add additional regions to the bad-pixel masks.
2.) Use ccdproc FIXPIX option to replace the bad pixels with interpolated value from adjacent pixels. Note that we never use these pixels in our final stacked image, but replacing their value with adjacent pixels helps avoid the ringing that would otherwise result during the tangent-plane projection.
3.) Generate the tangent-plane projection using mscimage.
4.) Make final preparations for combining of the indivdual tangent-plane projected images. This includes removing any large-scale gradients in the sky background (using mscskysub) and determining the proper relative scaling between the images (using mscimatch).
5.) Generate an initial stacked image with mscstack. In some cases this might be your final stacked image. To obtain an optimally combined image it is advantageous to improve the bad-pixel masks to remove satellite trails, remaining cosmic rays, asteroids, etc., from the individual images. We are currently still using software that is still in development to perform this last task, but in the next version of these notes I will try to describe the basic proceedure in more detail.
[Creating Improved Bad-pixel Masks (mkbleed,xzap)]
[Using Bad-pixel masks (fixpix)]
[Tangent-Plane Projection (mscimage)]
Saturated pixels in general, and the charge bleed trails from saturated stars in particular, contain no useful information and also are discontinuities in the measured flux that produce ringing noise or errors in the image when we generate our tangent-plane projected images using a sinc interpolator. Therefore we must add such pixels to our bad-pixels masks. The same is true for cosmic rays. We use the scripts mkbleed4 and xzap to do accomplish this goal. Each of these scripts produces new bad-pixels masks that indicate the locations of respectively charge bleed trails and cosmic rays. We use the script combinemasks to combine these new masks along with the default masks (indicated in the header, or generated earlier from dome-flat exposures) to produce one improved set of masks for each image and update the header of the images.
After the new bad-pixel mask has been generated, we
run fixpix to interpolate over the bad
regions. This step requires the IRAF task crutil to be have been
loaded. mkbleed4, xzap, and combinemasks are
scripts and must be defined by the user, either through a
loginuser.cl file that is processed when you start
IRAF or through the command:
msc> task mkbleed4='/any/path/to/the/scriptfile/mkbleed4.cl'
The source for these files are available here:
mkbleed4.cl, written by Chris Greer
mkbleed4 calls fileroot.cl by Marc Dickinson
xzap.cl, written by Marc Dickinson
combinemasks.cl, written by Chris Greer
combinemasks_old.cl, for use on Mosaic 1 data taken with the Engineering grade CCDs.
The scripts mkbleed4 and/or xzap require the scripts iterstat.cl, fileroot.cl, minv.cl, and makemask.cl. The script mkbleed4 also calls the iraf task crgrow, which can be found under the crutil package in IRAF.
For each of these the tasks above, other IRAF tasks, or your own scripts might work as well or better. For example, although those of us involved in the NDWFS have not tested it yet, the task crreject should perform the same task as xzap, and is expected to be a bit faster. Similarly, the current version of combinemasks was written by REU summer student Chris Greer to work on files with exactly our default nomenclature. In otherwords, these scripts are not general tasks yet, but we hope to produce more general versions in the months ahead.
Before mkbleed4 is run, you need to determine the saturation
levels for each CCD and the level of the charge bleed trails. If
saturated pixels were flagged initially (which is now possible with
ccdproc, then the saturation values included in the images
headers could be used, but since we are flagging these pixels after
flatfields have been applied, these values have been rescaled and if
you are doing things according to this cook-book you should use
implot to investigate the images and determine the saturation
values for each chip. These values will be different from chip to
chip. We have also found that for some of the CCDs, the ADU level
recorded for the charge that bleeds from the saturated stars, while at
discrete, quantized, values, are not at the saturation level. We still
want to flag these pixels as bad. To do this requires a second
threshold value to be set and to look for contiguous pixels along
columns that are above this value. In mkbleed4 there are two
threshold values per chip. Threshold 1 is the value above which
pixels are in bleed trails, while threshold 2 is the value above which
pixels are saturated. Note that the parameter shift
determines how many pixels along a column need to be above the first
threshold in order for that group of pixels to be flagged. You might
need to adjust this parameter to match the properties of your images.
Note that the current version of the task insists on an input file
containing the list of rootnames for the files. For example, the
first three lines of the input file might be:
obj111pf
obj112pf
obj113pf
Here is an example parameter file:
PACKAGE = user TASK = mkbleed4 inlist = @file.list MEF image(s) for cosmic ray cleaning thrc1t1 = 37500. Value above which pixels are in bleeds for chip 1. thrc2t1 = 20500. Value above which pixels are in bleeds for chip 2. thrc3t1 = 33000. Value above which pixels are in bleeds for chip 3. thrc4t1 = 19000. Value above which pixels are in bleeds for chip 4. thrc5t1 = 6000. Value above which pixels are in bleeds for chip 5. thrc6t1 = 23000. Value above which pixels are in bleeds for chip 6. thrc7t1 = 27000. Value above which pixels are in bleeds for chip 7. thrc8t1 = 21500. Value above which pixels are in bleeds for chip 8. (thrc1t2= 39000.) Value above which pixels are saturated for chip 1. (thrc2t2= 22500.) Value above which pixels are saturated for chip 2. (thrc3t2= 36000.) Value above which pixels are saturated for chip 3. (thrc4t2= 22000.) Value above which pixels are saturated for chip 4. (thrc5t2= 27000.) Value above which pixels are saturated for chip 5. (thrc6t2= 27000.) Value above which pixels are saturated for chip 6. (thrc7t2= 30000.) Value above which pixels are saturated for chip 7. (thrc8t2= 23500.) Value above which pixels are saturated for chip 8. (radius = 1.) Grow radius (uses crgrow). (shift = 10.) Shift to be applied. (verbose= yes) Verbose. (inimgli= ) (mode = ql)
The next step is to split the object files with mscsplit in preparation for using xzap to identify the positions of pixels that have been hit by cosmic rays. Note that the filenames used in the input list for mscsplit and mscjoin should be the root names, and not include the ".fits" suffix. Example parameter file:
PACKAGE = mscred TASK = mscsplit input = @objectfiles.list List of input MEF files (output = @objectfiles.list) List of output root names (mefext = .fits) MEF filename extension (delete = yes) Delete MEF file after splitting? (verbose= yes) Verbose? (fd1 = ) (fd2 = ) (mode = ql)
The task xzap is a script that through several iterations on an individual frame identifies the location of cosmic rays. It needs to be run on each individual extension and can not yet deal with multi-extension fits files. The output here must be very specific if combinemasks is to work. The xzap output file names must have a 'z' appended to it. Here is an example. obj039pfs_1.fits ---xzap---> obj039pfsz_1.fits Here is a parameter listing.
PACKAGE = user TASK = xzap inlist = @inxzap.list Image(s) for cosmic ray cleaning outlist = @outxzap.list Output image(s) (zboxsz = 5) Box size for zapping (nsigma = 5.) Number of sky sigma for zapping threshold (nnegsig= 10.) Number of sky sigma for negative zapping (nrings = 1) Number of pixels to flag as buffer around CRs (nobjsig= 2.) Number of sky sigma for object identification (skyfilt= 15) Median filter size for local sky evaluation (skysubs= 1) Block averaging factor before median filtering (ngrowob= 0) Number of pixels to flag as buffer around objects (statsec= [550:1500,600:1500]) Image section to use for computing sky sigma (deletem= no) Delete CR mask after execution? (cleanpl= yes) Delete other working .pl masks after execution? (cleanim= yes) Delete working .imh images after execution? (verbose= yes) Verbose output? (checkli= yes) Check min and max pix values before filtering? (zmin = -32768.) Minimum data value for fmedian (zmax = 32767.) Minimum data value for fmedian (inimgli= ) (outimgl= ) (statlis= ) (mode = ql)
Where the input list would contain a list of images with names like
obj025pf_1.fits
obj025pf_2.fits
...
And the output list would contain a list of images with names like
obj025pfz_1.fits
obj025pfz_2.fits
...
After xzap has completed, you can rejoin the individual files into MEF format images. This is done with mscjoin using a parameter file like this:
PACKAGE = mscred TASK = mscjoin input = obj039pfsz List of input root names (output = obj039pfsz) List of output MEF names (delete = yes) Delete input images after joining? (verbose= yes) Verbose? (fd1 = ) (fd2 = ) (mode = ql)
You can now run combinemasks to join all the bad-pixel mask
files together and update the headers of the MEF files. The input is
very specific. Use the root names of the
same input extensions that went into xzap.
That is, without the 'z' appended. Here is an example:
obj039pfs
obj040pfs
Here is an example parameter listing. Note that there will be eight
bad-pixel mask files (.pl format) for each MEF file. They will have
names that start with crmask3_. This is a naming convention that we
use, and you are free to change it by modifying the script.
PACKAGE = user TASK = combinemasks inlist = @inxzap.list Images for wich to generate new BPMs (Assumes CCDMosaThin1) (imglist= ) (mode = ql)
[Creating Improved Bad-pixel Masks (mkbleed,xzap)]
[Using Bad-pixel masks (fixpix)]
[Tangent-Plane Projection (mscimage)]
The next step is to replace the flux values for the bad pixels in each of the individaul images with a value interpolated from surrounding pixels. This is easily done by using fixpix on each individual extension of the image. It is convienent to have all the obj???pfsz.fits files in a file and then run fixpix with the following parameters.
PACKAGE = proto TASK = fixpix images = @infixpix.list//[1] List of images to be fixed masks = BPM List of bad pixel masks (linterp= INDEF) Mask values for line interpolation (cinterp= INDEF) Mask values for column interpolation (verbose= yes) Verbose output? (pixels = no) List pixels? (mode = ql)This example will run fixpix on the first CCD of each image. Repeat for chips 2 through 8. Again, we anticipate making this step more streamlined, but we have not modernized it yet.
It is important at this point to inspect each MEF image for any bleed trails that may have been missed. If there are any, replace the area in the bad-pixel mask associated with the image using imreplace. Here is an example of parameter file necessary to replace a section of the bad-pixel mask for CCD6 of image obj056pfs.fits with a value indicating "bad" data.
PACKAGE = imutil TASK = imreplace images = crmask3_obj056pfs_6.pl[1800:1820,2748:3112] Images to be edited value = 1. 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)If you have to add regions to your masks, remember to rerun fixpix. After each image is checked and is free from bleed trails and saturated pixels it is time to generate your tangent-plane projected images.
[Creating Improved Bad-pixel Masks (mkbleed,xzap)]
[Using Bad-pixel masks (fixpix)]
[Tangent-Plane Projection (mscimage)]
We use the task mscimage to generate a tangent-plane projected image from each of our object frames. Each of the resulting images (which we name with a "mos" prefix), ultimately destined to be combined into a single, stacked, image; will have a the same uniform (linear arcsecond per pixel scale) scale and will have been projected relative to the same position. The reference position is taken from the header of either one of the images that will be combined (and should therefore be specified as the reference image for all of the the images processed with mscimage), or can be from a reference image constructed in some other manner. We construct a reference image with a specific scale, orientation, and tangent-plane reference position to be used for all images in a given survey sub-field using a development task wcsref in the xdwred package. I will detail this further in the next version of the notes, but in brief we have been choosing to produce our output images with a fixed scale of 0.258 arcseconds per pixel, North to the left, East down.
With version 4.1 of mscred, mscimage is now more sophisticated at handling the gaps and boundaries of the "obj" images and the impact that these can have when sinc17 is used for the interpolation. This is such a new (and welcome) upgrade that we have no experience yet in using the new version. We handled the problem of the edges of the individual chips by specifying a "blank" value to be used by the software when the input image had no value. The value for the "blank" parameter was chosen using iterstat on the central four CCDs of each image. So, first setup the basic parameters for mscimage. Note that you could use linear for the value of interpo, but we feel that this compromises the data quality for our applications. Using a sinc-function (specified by interpo=sinc17) preserves essentially all spatial frequencies and is the only mathematically correct interpolation method for well-sampled data. Bilinear interpolation, by contrast, smooths the data with a variable kernal, depending on where the iterpolated pixel falls with respect to the original image. The result is degraded data, and atrifact-like patterns of coherent noise. Note that while a sinc function is used for the data, a linear interpolation must be used for the projection of the bad-pixel masks.
Below is an example parameter file for mscimage.
We run iterstat on a section of each image and find the mean somewhere close to the middle of the image. This will be used to set the constant value for each input object frame. p Mscimage would be run for an individual frame with a command something like:
mscimage ("obj156pfsz","mos156",reference="pointingscalerefimage.fits",
boundary="constant",constan=6255.)
PACKAGE = mscred
TASK = mscimage
input = List of input mosaic exposures
output = List of output images
(referen= ) Reference image
(pixmask= yes) Create pixel mask?
(verbose= yes) Verbose output?
# Resampling parmeters
(blank = 0.) Blank value
(interpolant= sinc17) Interpolant for data
(minterpolant= linear) Interpolant for mask
(boundary= constant) Boundary extension
(constant= 0.) Constant boundary extension value
(fluxcon= no) Preserve flux per unit area?
(ntrim = 7) Edge trim in each extension
(nxblock= 2048) X dimension of working block size in pixels
(nyblock= 1024) Y dimension of working block size in pixels
# Geometric mapping parameters
(interac= no) Fit mapping interactively?
(nx = 10) Number of x grid points
(ny = 20) Number of y grid points
(fitgeom= general) Fitting geometry
(xxorder= 4) Order of x fit in x
(xyorder= 4) Order of x fit in y
(xxterms= half) X fit cross terms type
(yxorder= 4) Order of y fit in x
(yyorder= 4) Order of y fit in y
(yxterms= half) Y fit cross terms type
(fd_in = )
(fd_out = )
(fd_ext = )
(fd_coor= )
(mode = ql)