[Creating Improved Bad-pixel Masks]
[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 coordinate 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.) Identify cosmic rays and add them to our bad pixel masks.
2.) Use ccdproc FIXPIX option to replace the bad pixels with interpolated values from adjacent pixels. Note that we never use these pixels in our final stacked image, but replacing their value with an interpolation from 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 under development to perform this last task, but in the next version of these notes we hope to describe the basic proceedure in more detail.
[Creating Improved Bad-pixel Masks]
[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. That is why we added such pixels to our bad pixel masks during our initial run of ccdproc. There are other pixels that will cause problems (will ring) during projection of the images if we do not get their positions into the bad pixel masks. Cosmic-ray hits are a prime example. Therefore we must add such pixels to our bad-pixels masks. We used to use the script xzap to identify cosmicray events in the images. It worked well, but was very slow. We now use craverage to find the cosmicray effects on our images. This task is still a bit slow, but faster than xzap.
Here is the parameter file for craverage (you will need to have loaded the package crutil):PACKAGE = crutil TASK = craverage input = List of input images output = List of output images (crmask = ) List of output cosmic ray and object masks (average= ) List of output block average filtered images (sigma = ) List of output sigma images (navg = 7) Block average box size (nrej = 3) Number of high pixels to reject from the average (nbkg = 5) Background annulus width (nsig = 50) Box size for sigma calculation (var0 = 0.) Variance coefficient for DN^0 term (var1 = 0.) Variance coefficient for DN^1 term (var2 = 0.) Variance coefficient for DN^2 term (crval = 1) Mask value for cosmic rays (lcrsig = 10.) Low cosmic ray sigma outside object (hcrsig = 3.75) High cosmic ray sigma outside object (crgrow = 0.) Cosmic ray grow radius (objval = 0) Mask value for objects (lobjsig= 10.) Low object detection sigma (hobjsig= 5.5) High object detection sigma (objgrow= 0.) Object grow radius (mode = ql)
You can run this in script mode by creating a script that looks like the example below. A hint for easy creation of this script follows.
craverage obj016f.fits[1] output="" crmask="crmask016_1" craverage obj016f.fits[2] output="" crmask="crmask016_2" craverage obj016f.fits[3] output="" crmask="crmask016_3" craverage obj016f.fits[4] output="" crmask="crmask016_4" craverage obj016f.fits[5] output="" crmask="crmask016_5" craverage obj016f.fits[6] output="" crmask="crmask016_6" craverage obj016f.fits[7] output="" crmask="crmask016_7" craverage obj016f.fits[8] output="" crmask="crmask016_8" # craverage obj017f.fits[1] output="" crmask="crmask017_1" craverage obj017f.fits[2] output="" crmask="crmask017_2" craverage obj017f.fits[3] output="" crmask="crmask017_3" craverage obj017f.fits[4] output="" crmask="crmask017_4" craverage obj017f.fits[5] output="" crmask="crmask017_5" craverage obj017f.fits[6] output="" crmask="crmask017_6" craverage obj017f.fits[7] output="" crmask="crmask017_7" craverage obj017f.fits[8] output="" crmask="crmask017_8" etc...
There are many ways to generate ones input lists (including ccdlist). An example of one quick way that avoids a lot of editing is to run mscstat with parameters set as follows:
PACKAGE = mscred TASK = mscstat images = *.fits Images (extname= ) Extension name selection (gmode = no) Global mode statistics? (fields = image) Fields to be printed (lower = INDEF) Lower cutoff for pixel values (upper = INDEF) Upper cutoff for pixel values (binwidt= 0.1) Bin width of histogram in sigma (format = yes) Format output and print column labels? (fd1 = ) (fd2 = ) (mode = ql)
Do this: mscstat > cravg_script You now have a file called cravg_script that looks like:
# IMAGE obj234pfs.fits[im1] obj234pfs.fits[im2] obj234pfs.fits[im3] obj234pfs.fits[im4] obj234pfs.fits[im5] obj234pfs.fits[im6] obj234pfs.fits[im7] obj234pfs.fits[im8] obj235pfs.fits[im1] obj235pfs.fits[im2] obj235pfs.fits[im3] obj235pfs.fits[im4] obj235pfs.fits[im5] obj235pfs.fits[im6] obj235pfs.fits[im7] obj235pfs.fits[im8]
So basically you have a listing of the images and all their groups. After that you can just use your favorite editor to edit this list to look like the script above. Now you can run your script: cl < cravg_script This will take a significantly long time, depending on your system.
Once craverage finishes you'll have a crmask***.fits (or .pl) mask file for each group of each image. You'll want to blink a couple with the images to make sure it's getting all the cosmic rays (or the vast majority of them anyway).
After this you'll want to combine this cosmic ray mask with your previously created bad pixel masks. We'll do this with a script called crplusbpmask.cl.
The source for this file is available here:
crplusbpmask.cl, a script modified
by Jenna Claver from scripts by Chris Greer, Lisa Storrie-Lombardi,
and Buell Jannuzi.
This script takes a list of "object numbers". For instance: 145 146 147 etc.
It assumes that the prefix is obj and that there is no suffix before the .fits.
If you have a suffix of s or pfs you should modify this script to account for your naming convention. This task also assumes that your previous bad pixel masks are in directories named bpm123,etc. and containing masks that are named bpm_im1, bpm_im2, etc.. It assumes your cosmic ray masks are named crmask123_im1, etc. It will create combined masks in the same directory, bpm123/, with the names bpmcr_im1 etc.. If you're using the correct script for your naming convention it will update the headers with the new mask name and location. This takes almost no time at all. It processes about as fast as it can print what it's doing to the screen.
The script crplusbpmask 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 crplusbpmask='/any/path/to/the/scriptfile/crplusbpmask.cl'
To run the task type:
crplusbpmask @yourlist
[Creating Improved Bad-pixel Masks]
[Using Bad-pixel masks (fixpix)]
[Tangent-Plane Projection (mscimage)]
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.
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. Run fixpix with the following parameters.
PACKAGE = proto TASK = fixpix images = @fixpix_in 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)
The input file is a list of image names, with extentions (so 8 entrees per image), as shown below. You can create such a listing in the same manner that you created your craverage script. Just use mscstt with only "image" output and pipe to a file, which you can then edit.
obj033pfs.fits[1] obj033pfs.fits[2] obj033pfs.fits[3] obj033pfs.fits[4] obj033pfs.fits[5] obj033pfs.fits[6] obj033pfs.fits[7] obj033pfs.fits[8] obj034pfs.fits[1] obj034pfs.fits[2] obj034pfs.fits[3] obj034pfs.fits[4] obj034pfs.fits[5] obj034pfs.fits[6] obj034pfs.fits[7] obj034pfs.fits[8] etc.
It is important at this point to inspect each MEF image for any bleed trails, cosmic rays, or other obvious defects 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 ]
[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 the same uniform (linear arcsecond per pixel) scale and will have been projected relative to the same position. In general, the reference position can be 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), can be from a reference image constructed in some other manner, or can be an RA and DEC postion specified as parameters in the task. For the NDWFS we constructed 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. We have been choosing to produce our output images with a fixed scale of 0.258 arcseconds per pixel, North up, East to the left. For NDWFS team members, the reference files are on /chomper/d10/moscal directory and are called: bootescenterNup.fits and cetuscenterNup.fits. For general users, however, you can either leave the "reference" parameter blank and the task will use the first input image in your input list as the reference image, or you can specify your reference preferences by using the ra,dec,scale and rotation parameters. (Also set the wcssource parameter = "parameters" rather than "image".)
We use a sinc-function (specified by interpo=sinc17) for the necessary interpolation during projection because it 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 it from the command line so that we can create a script that will run each image with a constant level that we set by running iterstat on a section of each image to find the mean somewhere close to the middle of the image. Mscimage would be run for an individual frame with a command something like:
mscimage ("obj156pfsz","mos156",reference="bootescenterNup.fits",
boundary="constant",constan=6255.)
PACKAGE = mscred
TASK = mscimage
input = List of input mosaic exposures
output = List of output images
(format = image) Output format (image|mef)
(pixmask= yes) Create pixel mask?
(verbose= yes) Verbose output
# Output WCS parameters
(wcssource = "image") Output WCS source (image|parameters)
(reference = "bootescenterNup.fits") Reference image
(ra = INDEF) RA of tangent point (hours)
(dec = INDEF) DEC of tangent point (degrees)
(scale = INDEF) Scale (arcsec/pixel)
(rotation = INDEF) Rotation of DEC from N to E (degrees)
# 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= 4200) X dimension of working block size in pixels
(nyblock= 2100) 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)