Coadding the Exposures

Note: the following assumes you have created a file called swarp.zeropoints, by following the instructions given in the preceding page.

Some preparation is necessary before the single-extension FITS files can be properly coadded into a final mosaic. My practice has been to set all unexposed pixels to a large negative value (-10000, or some such value) so that they can be excluded during the coadding process by thresholding. Although this is a bit awkward, I've found it's necessary to use this approach because Megacam images are so large that coadding a large number of them, even at 2x2 binning, routinely crashes IRAF/imcombine. The crashing problem is significantly worse on the RedHat Linux machines than under Sun/Solaris. Frank Valdez of NOAO has invested a considerable amount of effort investigating the problem (including downloading Megacam data to reduce natively in Tucson) and reproduced the symptoms but could not determine the cause. It appears to be triggered only when outlier pixel rejection is being used, but that's not something you will want to turn off.

Given this situation, the first order of business is to inspect every single-extension frame so that contaminated regions can be clobbered with the large negative values. For this purpose, one can use imreplace to mask out rectangular regions, and satbgon to mask out trails from meteors, satellites, and planes. The usage is as follows:

me> imreplace EGS7.0195.swarp.fits[x1:x2,y1:y2] -10000 
me> satbgon EGS7.0195.swarp.fits x1 y1 x2 y2 width 

satbgon is a simple script that can be obtained here; if you use it you'll also need to have the parameter file. I like to accumulate all these imreplace commands in a single IRAF macro called imreplace.macro so that it can be executed via cl redirection:

cl < imreplace.macro

...because I find myself often iterating Megacam reductions and it's handy to have a macro such as this ready to clean things up.

Next, hedit all the *.swarp.fits files just created, so that their headers contain a keyword BPM that lists the corresponding *.weight.fits files. Then one can use the *weight.fits files to insert the 'magic' values (my convention: -10000.0) in the blank pixels of the input files, e.g., by using the script convert.cl which is can be obtained here.

ls *.swarp.fits > swarplist	
ls *weight.fits > weightlist	
!paste swarplist weightlist | awk '{printf "hedit %s BPM %s add+ delete- update+ verify-\n",$1,$2}' > hedit.macro	
!paste swarplist | awk '{printf "hedit %s BLANK -10000.0 add+ delete- update+ verify-\n",$1}' >> hedit.macro	
cl < hedit.macro                # puts in BPM fields, and BLANK = -10000.0	
list = "swarplist"	
convert	
del *weight.fits 

Next one must generate the input files used by imcombine. There are many ways to do this. I did it as follows. NOTE: For the imstat command, one should choose upper and lower thresholds large enough to measure variations in the RMS. They should be several times the pixel-pixel standard devations; if they're not then they must be increased. The example IRAF commands listed below are cribbed from a reduction in which the input images had RMS values that varied from about 10-20; that's the reason for the choice of plus/minus 70 counts for the imstat command; it gives at least 3-sigma sensitivity to variations in the RMS.

imstat.fields="image,npix,mean,midpt,mode,stddev,min,max"
imstat *.swarp.fits format- lower=-70.0 upper = 70.0 > swarp.stats 
!cat swarp.stats | awk '{printf "%s\n", $1}' > swarp.inlist 
!cat swarp.stats | awk '{printf "%6.2f\n",(-1.0 * $5)}' > swarp.modes 
!cat swarp.zeropoints | awk '{printf "%9.7f\n",(exp(-0.9211 * $2))}' > swarp.scales 
!paste swarp.scales swarp.stats | awk '{printf "%6.4f\n", (100.0 / $1 / $1 / $7 / $7 )}' > swarp.weights 

Note: you ought to check the swarp.weights file to make sure that the values have enough significant figures to be meaningful. If you're coadding a lot of exposures, you might need to use a number larger than 100.0 in the !paste command to ensure that the weights are accurate.

After the files are made, configure and run imcombine as follows:

del tmp*
imcombine.hsigma=3.0
imcombine.lsigma=3.0
imcombine.offsets = "none"
imcombine.masktype = "none"
imcombine.lthresh = -300.0
imcombine.zero = ""
imcombine.weight = "@swarp.weights"
imcombine.scale = "@swarp.scales"
imcombine.nrejmasks="f1.g.coadd.nrej"
imcombine.expmasks="f1.g.coadd.cov"
imcombine.outlimits=""
imcombine.combine= "average"
imcombine.reject ="avsigclip"
imcombine @swarp.inlist output="f1.g.coadd.fits" bpmask="f1.g.coadd.bpm"
imcopy f1.g.coadd.cov.pl f1.g.coadd.cov.fits

Last step: The above procedure will create an optimally-weighted mosaic under conditions of stable seeing. If that mosaic is coextensive with SDSS, the final reduction step is to use sdssmatch to compute the absolute calibration from the Linux prompt:

sdssmatch final_coadd_name.fits 

sdssmatch will use SExtractor to find, photometer, and match SDSS sources. It will then generate a file called, in this generic example, final_coadd_name.sdss.matched.tab containing the SDSS and Megacam photometry matched row-by-row that can be used to determine the flux calibration with whatever method you find most convenient. Brian McLeod has kindly provided this prescription, which I used iteratively to determine a Megacam-SDSS color correction simultaneously with the absolute calibration:

< final_coadd_name.sdss.matched.tab column -a resid color | 
compute 'color=r-i;resid=Mag-i+0.13*color+32.79' | 
row 'Mag<-13&&i_err<0.1' | plottable color,resid marker circle

The command above will plot up the data together with a correction based on the r-i color; if the zeropoint of 32.79 mag (this is SDSS, so it's on the AB system) is correct, then the residuals will form a horizontal line at zero on the vertical axis. You can iterate the command until you've chosen a coefficient for the color term that makes the residuals' distribution horizontal in this magnitude-color space, and iterate the zeropoint until that distribution lies along the residual==0 axis. This sounds more complicated than it really is.

Last updated August 15, 2008.