Reducing Palomar 200-inch LFC Data

Reducing Palomar 200-inch LFC Data
Part 1 - Image Processing

Roy Gal
Updated 1 May 2006 - how to deal with pointing offsets, some easier file manipulation

NEW FOR FALL 2005!: DATA CAN NOW BE SAVED AT THE TELESCOPE AS SIGNED 16-BIT FITS FILES WITH BZERO=32768 (using the iraf command). THUS, LFCASSEMBLE NO LONGER MODIFIES UNSIGNED/BZERO=0 DATA. IF YOU HAVE OLDER DATA you can use lfcassemble_old.cl which is now included in lfc_scripts.tar. This page provides a relatively detailed, step-by-step description of how I've reduced data from the LFC (Large Format Camera), mounted at the prime focus of the 200-inch Hale Telescope at Palomar Observatory. These notes combine:

Detailed mosaic reduction techniques from the NOAO Deep Wide Field Survey
Matt Hunt's lfcred reduction page. Note: as of March 2005 this page is gone.
My own experience reducing LFC data taken from 2001 to the present.

In order to perform these reductions, you will need:
1. The external IRAF packages mscred and mscdb from the IRAF external package ftp site. Currently (March 2004), mscred is at V4.8, and mscdb is at V21243. These must be installed in your IRAF installation in /iraf/extern or some similar location, and task definitions added to your login.cl; see the installation instructions inside the downloaded folders.

2. The lfc-specific database files. These should be located in a palomar/lfc subdirectory of your mscdb installation, such as /iraf/extern/mscdb/palomar/lfc. You should download this tar file (lfc_mscdb.tar) and untar it in that directory.

3. A variety of iraf scripts for preparing and processing the data. You should download this tarfile (lfc_scripts.tar) of lfcred programs and put them in your local IRAF script location, with task definitions added to your loginuser.cl file. The versions of these scripts you download from here includes an updated version of lfcassemble.cl which adds read noise and saturation keywords to the image headers; if you have an old one from Matt Hunt you can replace it with this one. This tarfile also includes crplusbpmask.cl, which will be used for cosmic ray masking. Place this file in your iraf scripts directory, and add the task definition to your loginuser.cl. For instance, to add all the lfc tasks, include the following in your login.cl or loginuser.cl:

 
#LFC data reduction tasks from Roy Gal, Matt Hunt and NOAO
task lfcassemble ="scripts$lfcassemble.cl" 
task lfcbleed="scripts$lfcbleed.cl" 
task lfcmask="scripts$lfcmask.cl" 
task crplusbpmask="scripts$crplusbpmask.cl"

In order to properly reduce a night's worth of observations, you need the following data:

A set of bias frames taken that night; 7-10 frames is enough. Matt Hunt has observed that there is quite a bit of structure in the bias frames (and I concur); I also see that it changes somewhat night to night, and almost certainly between runs. This is due to 60Hz electronics noise.
A set of dome flats in each filter you have used for observing targets.
Your target and standard star observations.
An astrometry calibration field which is rich in USNO-A2 stars (preferably low galactic latitude). If astrometry is very important to you, you may want frames in all your filters. This is actually not necessary; the solutions provided in the database files are good enough to allow matching to USNO stars to derive final solutions good to ~0.3".
Note: Twilight flats are not used in this process.

1. Preparation	2. Initial Reduction of Calibrations	3. Make the Bias Frame	4. Make the Dome Flat
5. Process Data Frames	6. Make Object Masks	7. Make Super Sky Flat	8. Make Fringe Correction (i and z data only)
9. Apply Fringe Correction (i and z data only)	10. Apply Super Sky Flat	11. Improve the World Coordinate System	12. Cosmic Ray Cleaning
13. Generate Tangent-Plane Projected Images	14. Fix Satellite Trails	15. Subtract Global Sky Pattern	16. Stack the Mosaicked Images
17. Prepare Standards for Photometry

Preparation

Drink a Scotch. Have the rest of the bottle handy because this will take a while.

Okay, let's get started. Begin by firing up iraf (you must be in your iraf home directory to do this), and then go to your data directory:

% cl
cl> cd /Users/gal/Work/supercluster/lfc/may13_2001/

We need to load the ccdred and mscred packages if they are not loaded already, and set the instrument keywords to point to the LFC description:

cl> 
cl> ccdred
cl> epar ccdred
PACKAGE = imred
   TASK = ccdred

(pixelty=            real real) Output and calculation pixel datatypes
(verbose=                  yes) Print log information to the standard output?
(logfile=              logfile) Text log file
(plotfil=                     ) Log metacode plot file
(backup =                     ) Backup directory or prefix
(instrum= mscdb$palomar/lfc/lfc.dat) CCD instrument file
(ssfile =              subsets) Subset translation file
(graphic=             stdgraph) Interactive graphics output device
(cursor =                     ) Graphics cursor input
(version=      2: October 1987)
(mode   =                   ql)

Exit with ":wq" to save changes.

cl> mscred
cl> epar mscred
PACKAGE = clpackage
   TASK = mscred

(pixelty=            real real) Output and calculation pixel datatypes
(verbose=                  yes) Print log information to the standard output?
(logfile=              logfile) Text log file
(plotfil=                     ) Log metacode plot file
(backup =                 none) Backup data (none|once|all)?
(bkuproo=                 Raw/) Backup root (directory or prefix)
(instrum= mscdb$palomar/lfc/lfc.dat) CCD instrument file
(ampfile=                 amps) Amplifier translation file
(ssfile =              subsets) Subset translation file
(im_bufs=                   4.) Image I/O buffer size (in Mbytes)
(graphic=             stdgraph) Interactive graphics output device
(cursor =                     ) Graphics cursor input
(version=  V4.8: July 30, 2003)
(mode   =                   ql)
($nargs =                    0)

Exit with ":wq" to save changes. Note that I've set backup=none to save disk space. I strongly urge you to make a subdirectory where you keep your raw data, and make a copy of the raw files to the location!

The lfc data is originally stored as 6 separate FITS files, one per CCD. However, the mscred package wants these in a single multi-extension FITS, or MEF, file. As a bonus, the initial files don't have the .fits extensions so you'll need to rename them.

NEW FOR FALL 2005!: It is now possible to save data at the telescope with the .fits extension. If so, you can skip this small step.

cl> !ls -1 lfc* | awk '{print "mv",$1,$1".fits"}' > go
cl> !chmod +x go
cl> !.go

If your files already had the .fits extension, continue here.

Now, let's go to the directory where our data is, and we'll use lfcassemble to do this. lfcassemble will:

Put the 6 existing files for each readout (lfc.142.0...lfc.142.5) into an MEF file with extensions lfc.142[1]...lfc.142[6]. Note the offset in extension name relative to chip number!!!
Adds DETSEC header which allows "mscdisplay" to display the array.
Produces or changes OBSTYPE, RA, DEC, DATASEC, CCDSEC
Adds BZERO headers and changes pixel values so that IRAF correctly interperets values larger than 32k. NEW FOR FALL 2005!: DATA CAN NOW BE SAVED AT THE TELESCOPE AS SIGNED 16-BIT FITS FILES (using the iraf command). THUS, LFCASSEMBLE NOW SKIPS THIS STEP. IF YOU HAVE OLDER DATA you can use lfcassemble_old.cl which is now included in lfc_scripts.tar.
Adds BPM and XTALKFIL header keywords pointing to the bad pixel masks and crosstalk definitions in the mscdb I've provided, as well as specifying the trim and overscan regions.
Adds RDNOISE, BLEEDVAL and SATURATE header keywords (RRG version of code ONLY!). If you are using Matt Hunt's code, you'll need to use hedit to add these values from the table below.

I usually keep a copy of the post-lfcassemble data as my raw backup. Also, we want to add header keywords for the bleed trails. We'll set the bleed threshold to 5000 below saturation (see table below). My version of lfcassemble.cl does this as well.

 
cl> !mkdir raw 
cl> !mv lfc*.fits raw/ 
cl> cd raw
cl> !ls -1 lfc*.0.fits | sed 's/\.0\.fits//g' > list.all 
cl> lfcassemble @list.all

The current version of lfcassemble adds the bleed, readnoise, and saturation values to the headers, but if you need to do it manually, just use:

cl> hedit lfc*.fits[1] bleedval 46000 add+ ver-
cl> hedit lfc*.fits[2] bleedval 48400 add+ ver- 
cl> hedit lfc*.fits[3] bleedval 35800 add+ ver- 
cl> hedit lfc*.fits[4] bleedval 39000 add+ ver- 
cl> hedit lfc*.fits[5] bleedval 16700 add+ ver- 
cl> hedit lfc*.fits[6] bleedval 56400 add+ ver-

Chip	Filenames	Read Noise	Saturation	Bleed
0	lfc*.fits[1]	20	59100	46000
1	lfc*.fits[2]	13	53400	48400
2	lfc*.fits[3]	15	40800	35800
3	lfc*.fits[4]	16	44000	39000
4	lfc*.fits[5]	14	21700	16700
5	lfc*.fits[6]	14	61400	56400

Copy the files to your working directory.

cl> !cp lfc[0-9][0-9][0-9].fits ../ 
cl> cd ../

We are now in our working directory, and we have a subdirectory called raw/ with files we won't touch.

It's a good idea to check your images to make sure they are ok. The first image (even a bias) taken each night may be saturated if the CCDs were not wiped of charge buildup by an earlier readout. Start up your favorite display tool and have a look at the images:

cl> !ds9 & 
cl> mscdisplay lfc001 frame=1

Another problem to keep an eye out for is horizontal striping in target observations. I imagine this is due to some sort of readout error.

Initial Reduction of Calibrations

We begin by performing a number of initial reduction steps on the calibration (i.e., bias, dark, domeflat) frames. These steps include applying the crosstalk and overscan corrections, and trimming off the overscan regions. Make a text file called calibrations.list listing these frames (one per line), then get ready to run ccdproc:

cl> epar ccdproc
PACKAGE = mscred
   TASK = ccdproc
    
images  =   @calibrations.list  List of Mosaic CCD images to process
(output =                     ) List of output processed images
(bpmasks=                     ) List of output bad pixel masks
(ccdtype=                     ) CCD image type to process
(noproc =                   no) List processing steps only?

(xtalkco=                  yes) Apply crosstalk correction?
(fixpix =                   no) Apply bad pixel mask correction?
(oversca=                  yes) Apply overscan strip correction?
(trim   =                  yes) Trim the image?
(zerocor=                   no) Apply zero level correction?
(darkcor=                   no) Apply dark count correction?
(flatcor=                   no) Apply flat field correction?
(sflatco=                   no) Apply sky flat field correction?
(split  =                   no) Use split images during processing?
(merge  =                   no) Merge amplifiers from same CCD?

(xtalkfi=            !xtalkfil) Crosstalk file
(fixfile=                  BPM) List of bad pixel masks
(saturat=            !saturate) Saturated pixel threshold
(sgrow  =                    1) Saturated pixel grow radius
(bleed  =                INDEF) Bleed pixel threshold
(btrail =                   15) Bleed trail minimum length
(bgrow  =                    1) Bleed pixel grow radius
(biassec=             !biassec) Overscan strip image section
(trimsec=             !trimsec) Trim data section
(zero   =                     ) List of zero level calibration images
(dark   =                     ) List of dark count calibration images
(flat   =                     ) List of flat field images
(sflat  =                     ) List of secondary flat field images
(minrepl=                  0.9) Minimum flat field value

(interac=                   no) Fit overscan interactively?
(functio=             legendre) Fitting function
(order  =                    1) Number of polynomial terms or spline pieces
(sample =                    *) Sample points to fit
(naverag=                    1) Number of sample points to combine
(niterat=                    1) Number of rejection iterations
(low_rej=                   3.) Low sigma rejection factor
(high_re=                   3.) High sigma rejection factor
(grow   =                   0.) Rejection growing radius
(fd     =                     )
(fd2    =                     )
(mode   =                   ql)

Exit the epar with ":go" to set it running.

Make the Bias Frame

Now we are ready to create our first correction frame: the bias, or zero, frame. Make a list of your bias frames in a file called zero.list. Then, we use zerocombine:


cl> epar zerocombine
PACKAGE = mscred
   TASK = zerocombine
    
input   =           @zero.list  List of zero level images to combine
(output =                 Zero) Output zero level name
(combine=              average) Type of combine operation
(reject =             crreject) Type of rejection
(ccdtype=                     ) CCD image type to combine
(process=                  yes) Process images before combining?
(delete =                   no) Delete input images after combining?
(scale  =                 none) Image scaling
(statsec=                     ) Image section for computing statistics
(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 =                   3.) 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)
(pclip  =                 -0.5) pclip: Percentile clipping parameter
(blank  =                   0.) Value if there are no pixels
(mode   =                   ql)

Exit with ":go" to execute. As usual, it is a good to idea to take a look at the output:

cl> mscdisplay Zero frame=1

Your output should look something like this.

Because there is 60 MHz pattern noise in the bias frames, we will do a light medianing to reduce its effect:

cl> epar mscmedian
PACKAGE = mscred
   TASK = mscmedian
    
input   =                 Zero  Input mosaic images
output  =              medZero  Output mosaic images
xwindow =                   15  X window size of median filter
ywindow =                   15  Y window size of median filter
(outtype=               median) Output type (median|difference)
(zloreje=                INDEF) Lowside pixel value cutoff
(zhireje=                INDEF) High side pixel value cutoff
(verbose=                  yes) Print messages about actions taken by the task

                                # Fast median
(fmedian=                   no) Use fast median algorithm?
(hmin   =                   -5) Minimum histogram bin
(hmax   =                  100) Maximum histogram bin
(zmin   =                  -5.) Pixel value corresponding to hmin
(zmax   =                 100.) Pixel value corresponding to hmax
(fd     =                     )
(mode   =                   ql)

Make the Dome Flat

The next step is to combine the domeflats taken in each filter. We'll want to make sure the bad pixel masks (in the mscdb$palomar/lfc directory) are applied. NOTE: I have observed that the masks provided by Matt Hunt don't match the bad pixels I see in data taken in September 2001. I urge you to check the masks! Mask checking can be done by displaying an image with the bad pixel mask overlaid

cl> mscdisplay lfc012.fits bpdisplay=overlay bpmask=BPM

If you are not happy with the masks, you can make new ones if you have domeflats taken in the same filter but with different exposure times. If you take the ratio of these two images, the bad pixels will not be corrected, and you can use the ccdmask task in IRAF on each extension. Or, you can make a text file listing the bad pixels, and convert it to a .pl file using text2mask. Please read the help for those tasks for more details. Assuming we have good mask files, we will first set up ccdproc to use these masks. We will NOT execute ccdproc, just edit the parameters, as the next command (flatcombine) checks these parameters.

cl> epar ccdproc
PACKAGE = mscred
   TASK = ccdproc
    
images  =                       List of Mosaic CCD images to process
(output =                     ) List of output processed images
(bpmasks=                     ) List of output bad pixel masks
(ccdtype=                     ) CCD image type to process
(noproc =                   no) List processing steps only?

(xtalkco=                  yes) Apply crosstalk correction?
(fixpix =                  yes) Apply bad pixel mask correction?
(oversca=                  yes) Apply overscan strip correction?
(trim   =                  yes) Trim the image?
(zerocor=                  yes) Apply zero level correction?
(darkcor=                   no) Apply dark count correction?
(flatcor=                   no) Apply flat field correction?
(sflatco=                   no) Apply sky flat field correction?
(split  =                   no) Use split images during processing?
(merge  =                   no) Merge amplifiers from same CCD?
(xtalkfi=            !xtalkfil) Crosstalk file
(fixfile=                  BPM) List of bad pixel masks
(saturat=            !saturate) Saturated pixel threshold
(sgrow  =                    1) Saturated pixel grow radius
(bleed  =                INDEF) Bleed pixel threshold
(btrail =                   15) Bleed trail minimum length
(bgrow  =                    1) Bleed pixel grow radius
(biassec=             !biassec) Overscan strip image section
(trimsec=             !trimsec) Trim data section
(zero   =              medZero) List of zero level calibration images
(dark   =                     ) List of dark count calibration images
(flat   =                     ) List of flat field images
(sflat  =                     ) List of secondary flat field images
(minrepl=                  0.9) Minimum flat field value

(interac=                   no) Fit overscan interactively?
(functio=             legendre) Fitting function
(order  =                    1) Number of polynomial terms or spline pieces
(sample =                    *) Sample points to fit
(naverag=                    1) Number of sample points to combine
(niterat=                    1) Number of rejection iterations
(low_rej=                   3.) Low sigma rejection factor
(high_re=                   3.) High sigma rejection factor
(grow   =                   0.) Rejection growing radius
(fd     =                     )
(fd2    =                     )
(mode   =                   ql)

Use ":wq" to save the parameters and exit the epar.

Now, it is time for to make the dome flats using flatcombine. Make a file domeflats.list listing your dome flat images, then:

cl> epar flatcombine


PACKAGE = mscred
   TASK = flatcombine
    
input   =      @domeflats.list  List of flat field images to combine
(output =                DFlat) Output flat field root name
(combine=              average) Type of combine operation
(reject =              ccdclip) Type of rejection
(ccdtype=                     ) CCD image type to combine
(process=                  yes) Process images before combining?
(subsets=                  yes) Combine images by subset parameter?
(delete =                   no) Delete input images after combining?
(scale  =                 mode) Image scaling
(statsec=                     ) Image section for computing statistics
(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 =                   3.) 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)
(pclip  =                 -0.5) pclip: Percentile clipping parameter
(blank  =                   1.) Value if there are no pixels
(mode   =                   ql)

Use ":go" to exit and execute. The resulting dome flats should show plenty of structure; here is an example of an r-band dome flat.

Process Data Frames

Next, we'll apply the corrections we've generated to all of our data frames (long exposures, astrometric fields, standards) frames. We'll need a list of these frames, called objects.list. We also will be outputting bad pixel masks for each frame, so we'll need a list of mask files, which we can make easily with

cl> !sed 's/lfc/bpm/g' objects.list | sed 's/\.fits//g' > masks.list

Then, as usual, we need to set up ccdproc, and this time we'll execute it.

cl> epar ccdproc
PACKAGE = mscred
   TASK = ccdproc
    
images  =        @objects.list  List of Mosaic CCD images to process
(output =                     ) List of output processed images
(bpmasks=          @masks.list) List of output bad pixel masks
(ccdtype=                     ) CCD image type to process
(noproc =                   no) List processing steps only?

(xtalkco=                  yes) Apply crosstalk correction?
(fixpix =                  yes) Apply bad pixel mask correction?
(oversca=                  yes) Apply overscan strip correction?
(trim   =                  yes) Trim the image?
(zerocor=                  yes) Apply zero level correction?
(darkcor=                   no) Apply dark count correction?
(flatcor=                  yes) Apply flat field correction?
(sflatco=                   no) Apply sky flat field correction?
(split  =                   no) Use split images during processing?
(merge  =                   no) Merge amplifiers from same CCD?
(xtalkfi=            !xtalkfil) Crosstalk file
(fixfile=                  BPM) List of bad pixel masks
(saturat=            !saturate) Saturated pixel threshold
(sgrow  =                    1) Saturated pixel grow radius
(bleed  =            !bleedval) Bleed pixel threshold
(btrail =                   15) Bleed trail minimum length
(bgrow  =                    1) Bleed pixel grow radius
(biassec=             !biassec) Overscan strip image section
(trimsec=             !trimsec) Trim data section
(zero   =              medZero) List of zero level calibration images
(dark   =                     ) List of dark count calibration images
(flat   =          DFlat*.fits) List of flat field images
(sflat  =                     ) List of secondary flat field images
(minrepl=                  0.9) Minimum flat field value

(interac=                   no) Fit overscan interactively?
(functio=             legendre) Fitting function
(order  =                    1) Number of polynomial terms or spline pieces
(sample =                    *) Sample points to fit
(naverag=                    1) Number of sample points to combine
(niterat=                    1) Number of rejection iterations
(low_rej=                   3.) Low sigma rejection factor
(high_re=                   3.) High sigma rejection factor
(grow   =                   0.) Rejection growing radius
(fd     =                     )
(fd2    =                     )
(mode   =                   ql)

Exit with ":go" to execute. This will take a long time to run on all the frames in a single night (a few hours at least).

Make Object Masks

We will need object masks for all our data frames (already in the file objects.list. These masks are used to avoid bright objects when generating flat fields, etx. We will make rough object masks for each image using the objmasks task in the nproto package. This task has some hidden variables; we need to edit these in the objmasks1 task, and then epar objmasks. It is also important that the output masks be written as pixel list (.pl) files! You will also need to create some lists of output filenames. Below is an example of how to make lists of all the files in your directory. If you have i or z data, you will need to make separate lists!!! This is because they need different parameters for some of the tasks below! It turns out that fringing is so strong, especially in z, that parts of the fringe pattern are masked as objects if you use the same parameters in objmask for r and the redder filters.

cl> !sed 's/\.fits//g' objects.list | sed 's/lfc/om/g' > om.list
cl> !sed 's/om/sky/g' om.list > sky.list
cl> nproto
np> set masktype = pl      <=====IMPORTANT!!!!
np> epar objmasks1

PACKAGE = nproto
   TASK = objmasks1
    
(exps   =                     ) List of exposure maps
(gains  =                     ) List of gain maps
(catalog=                     ) List of catalogs
(catdefs=                     ) List of catalog definitions
(dodetec=                  yes) Detect objects?
(dosplit=                   no) Split merged objects?
(dogrow =                  yes) Grow object regions?
(doevalu=                   no) Evaluate objects?
(skytype=                block) Type of sky estimation
(fitstep=                   10) Line step for sky sampling
(fitblk1=                   10) Block average for line fitting
(fithcli=                   3.) High sky clipping during 1D sky estimation
(fitlcli=                   3.) Low sky clippling during 1D sky estimation
(fitxord=                    2) Sky fitting x order
(fityord=                    2) Sky fitting y order
(fitxter=                 half) Sky fitting cross terms
(blknsub=                    2) Number of subblocks per axis
(updates=                  yes) Update sky during detection?
(sigavg =                   4.) Sigma of mean flux cutoff
(sigmax =                   4.) Sigma of maximum pixel
(bpval  =                INDEF) Output bad pixel value
(splitma=                INDEF) Maximum sigma above sky for splitting
(splitst=                  0.4) Splitting steps in convolved sigma
(splitth=                   5.) Splitting threshold in sigma
(sminpix=                    8) Minimum number of pixels in split objects
(ssigavg=                  10.) Sigma of mean flux cutoff
(ssigmax=                   5.) Sigma of maximum pixel
(magzero=                INDEF) Magnitude zero point
(mode   =                   ql)

Exit with ":wq" to save these settings.

Note the different settings below for (i and z) versus (g and r) data in objmasks! You will have to set up separate lists of objects, masks, and sky files for g/r and i/z.

np> epar objmasks

PACKAGE = nproto
   TASK = objmasks

images  =        @objects.list  List of images or MEF files
objmasks=             @om.list  List of output object masks
(omtype =              numbers) Object mask type
(skys   =            @sky.list) List of input/output sky maps
(sigmas =                     ) List of input/output sigma maps
(masks  =                 !BPM) List of input bad pixel masks
(extname=                     ) Extension names
(logfile=               STDOUT) List of log files

(blkstep=                    1) Line step for sky sampling
(blksize=                  -10) Sky block size (+=pixels, -=blocks)
(convolv=            block 3 3) Convolution kernel
(hsigma =                    5) Sigma threshold above sky
*** HSIGMA MAY VARY BASED ON FILTER AND SKY BRIGHTNESS ***
(lsigma =                   7.) Sigma threshold below sky
(hdetect=                  yes) Detect objects above sky?
(ldetect=                   no) Detect objects below sky?
(neighbo=                    8) Neighbor type"
(minpix =                    6) Minimum number of pixels in detected objects
(ngrow  =                    3) Number of grow rings
(agrow  =                   2.) Area grow factor
(mode   =                   ql)

Exit with ":go" to execute.

If you have ONLY g and r data, AND you have existing Sky Flats, skip to the apply sky flat section.

If you have i and/or z data, AND you have existing Fringe correction frames, skip to the apply fringe correction section.

If you need to make Sky Flats, continue directly below.

Make a Super Sky Flat

Once the dome flat has been applied, check your images. The background should be much smoother, and the small-scale imperfections should be gone. However, if you adjust the contrast, you will still see large-scale background variations. These are due to uneven dome illumination for the dome flat, and because the color of the dome differs from that of the night sky, and the CCD sensitivity is color-dependent. So, we want to create a flatfield that has the same color as the dark night sky. We do this by combining our imaging data, with objects masked, to create a super-sky flat. We will use only long exposures where there are significant counts in the sky background!

Next, we need to combine the images in each filter to produce the super sky flats. Make sure you return to the mscred package. Remember, we will want to use only the long exposures for this! Make a file called targets.list listing these long exposures. This might be easily done by copyng and editing the existing objects.list.

np> mscred
ms> epar sflatcombine

PACKAGE = mscred
   TASK = sflatcombine
    
input   =        @targets.list  List of images to combine
(output =                SFlat) Output sky flat field root name
(combine=              average) Type of combine operation
(reject =            avsigclip) Type of rejection
(ccdtype=                     ) CCD image type to combine
(subsets=                  yes) Combine images by subset parameter?
(masktyp=             !objmask) Mask type
(maskval=                   0.) Mask value
(scale  =                 mode) Image scaling
(statsec=                     ) Image section for computing statistics
(nkeep  =                    1) Minimum to keep (pos) or maximum to reject (neg)
(nlow   =                    1) minmax: Number of low pixels to reject
(nhigh  =                    1) minmax: Number of high pixels to reject
(mclip  =                  yes) Use median in sigma clipping algorithms?
(lsigma =                   4.) Lower sigma clipping factor
(hsigma =                  2.5) Upper sigma clipping factor
(rdnoise=              RDNOISE) ccdclip: CCD readout noise (electrons)
(gain   =                 GAIN) ccdclip: CCD gain (electrons/DN)
(snoise =                   0.) ccdclip: Sensitivity noise (fraction)
(pclip  =                 -0.5) pclip: Percentile clipping parameter
(blank  =                   1.) Value if there are no pixels
(grow   =                   3.) Radius (pixels) for neighbor rejection
(fd     =                     )
(mode   =                   ql)

The images we just created contain information on the smooth sensitivity variations across the chip, but has pixel-to-pixel noise. So, we smooth them with a sliding boxcar that is pretty large. In the task below, I like to set the parameters zlo, zhi, hmin, hmax. zmin, zmax to the median of the image +/-1000! These values will be different for each filter!

ms> epar mscmedian

PACKAGE = mscred
   TASK = mscmedian
    
input   =               SFlati  Input mosaic images
output  =            MedSFlati  Output mosaic images
xwindow =                  328  X window size of median filter
ywindow =                  328  Y window size of median filter
(outtype=               median) Output type (median|difference)
(zloreje=               10500.) Lowside pixel value cutoff
(zhireje=               12500.) High side pixel value cutoff
(verbose=                  yes) Print messages about actions taken by the task

                                # Fast median
(fmedian=                  yes) Use fast median algorithm?
(hmin   =                10500) Minimum histogram bin
(hmax   =                12500) Maximum histogram bin
(zmin   =                10500) Pixel value corresponding to hmin
(zmax   =                12500) Pixel value corresponding to hmax
(fd     =                     )
(mode   =                   ql)

Exit with ":go" to execute. You'll need to do the same thing for your SFlatz. Examine the output image to ensure that no remnants of fringing are present. An example MedSFlati is here. VERY VERY IMPORTANT: The MedSFlats you have just created are the final super sky flats in the g and r (or other similar) filters. They are NOT the final flats for any redder filters (such as i and z).

Make the Fringe Correction (i and z ONLY)

For the i and z filters, we will use this first median-filtered SFlat to create a fringe correction frame, apply this to our data, and then make a new SFlat that contains only the residual flat field correction, without fringing.

To do this, we will subtract the median-filtered image from the initial (unfiltered) SFlat, which creates an output frame with only the fringing:

ms> epar mscarith

PACKAGE = mscred
   TASK = mscarith

operand1=               SFlati  Operand image or numerical constant
op      =                    -  Operator
operand2=            MedSFlati  Operand image or numerical constant
result  =              FringeI  Resultant image
(extname=                     ) Extension names to select
(title  =                     ) Title for resultant image
(divzero=                   0.) Replacement value for division by zero
(hparams=                     ) List of header parameters
(pixtype=                     ) Pixel type for resultant image
(calctyp=                     ) Calculation data type
(verbose=                   no) Print operations?
(noact  =                   no) Print operations without performing them?

(fd1    =                     )
(fd2    =                     )
(fd3    =                     )
(mode   =                   ql)

Exit with ":go" to execute. You'll need to do the same thing for your SFlatz. Examine the output image to check that no large structures remain, and that the fringing is well represented. An example Fringei is here. If you are happy with this output, we'll now remove the fringing from all of our data. Remember you'll need to do this for both the i and z data.

At this point, we have a fringe correction frame for i and z.

Apply Fringe Correction

Now, we apply the fringe correction to ALL of our data frames in i and z (long and short exposures), and run sflatcombine yet again. To apply the fringe correction, we use the task rmfringe. We'll want the output files to have some distinctive name. Make lists of all the data frames in the i and z filters (targets, standards, etc.) called idata.list and zdata.list. We will also need lists of object masks and sky frames that correspond to this list of frames.

ms> !awk '{print $1"f"}' idata.list > outfringei.list
ms> !sed 's/\.fits//g' idata.list | sed 's/lfc/om/g' > omi.list
ms> !sed 's/om/sky/g' omi.list > skyi.list 
ms> mscred          (just making sure we are in mscred package!)
ms> epar rmfringe

PACKAGE = mscred
   TASK = rmfringe

input   =         @idata.list  List of input images
output  =    @outfringei.list  List of output corrected images
fringe  =             Fringei  Fringe or list of fringe patterns
masks   =           @omi.list  List of object/bad data masks
(fringem=                     ) Fringe masks
(backgro=           @skyi.list) Lisk of input image backgrounds
(ncblk  =                    5) Column smoothing
(nlblk  =                    5) Line smoothing
(extfit =                     ) Extensions to use in scaling fit
(logfile=                     ) Logfile
(verbose=                  yes) Verbose?
(mode   =                   ql)

Exit with ":go" to execute. Repeat for the z data.

If you have existing Sky Flats for your i and z data, skip to the apply sky flat section.

If not, follow along to make the necessary sky flats.

Now we will use the output files to make super-sky flats that have no fringing, only the large-scale sensitivity gradients. We'll delete the old SFlati and SFlatz since they were just temporary. Remember to run sflatcombine on both the i and z data:

ms> imdel SFlati,SFlatz,MedSFlati,MedSFlatz
ms> epar sflatcombine
PACKAGE = mscred
   TASK = sflatcombine
    
input   =     @outfringez.list  List of images to combine
(output =                SFlat) Output sky flat field root name
(combine=              average) Type of combine operation
(reject =            avsigclip) Type of rejection
(ccdtype=                     ) CCD image type to combine
(subsets=                  yes) Combine images by subset parameter?
(masktyp=             !objmask) Mask type
(maskval=                   0.) Mask value
(scale  =                 mode) Image scaling
(statsec= [1000:1500,2000:3000]) Image section for computing statistics
(nkeep  =                    1) Minimum to keep (pos) or maximum to reject (neg)
(nlow   =                    1) minmax: Number of low pixels to reject
(nhigh  =                    1) minmax: Number of high pixels to reject
(mclip  =                  yes) Use median in sigma clipping algorithms?
(lsigma =                   4.) Lower sigma clipping factor
(hsigma =                  2.5) Upper sigma clipping factor
(rdnoise=              RDNOISE) ccdclip: CCD readout noise (electrons)
(gain   =                 GAIN) ccdclip: CCD gain (electrons/DN)
(snoise =                   0.) ccdclip: Sensitivity noise (fraction)
(pclip  =                 -0.5) pclip: Percentile clipping parameter
(blank  =                   1.) Value if there are no pixels
(grow   =                   3.) Radius (pixels) for neighbor rejection
(fd     =                     )
(mode   =                   ql)

Exit with ":go" to execute. Have a look at the output to make sure that it doesn't show artifacts from bright objects or other problems. If it looks ok, we will median filter these with a large smoothing kernel to produce the final Dark Sky Flats:

ms> epar mscmedian

PACKAGE = mscred
   TASK = mscmedian
    
input   =               SFlati  Input mosaic images
output  =            MedSFlati  Output mosaic images
xwindow =                  328  X window size of median filter
ywindow =                  328  Y window size of median filter
(outtype=               median) Output type (median|difference)
(zloreje=               10500.) Lowside pixel value cutoff
(zhireje=               12500.) High side pixel value cutoff
(verbose=                  yes) Print messages about actions taken by the task

                                # Fast median
(fmedian=                  yes) Use fast median algorithm?
(hmin   =                10500) Minimum histogram bin
(hmax   =                12500) Maximum histogram bin
(zmin   =                10500) Pixel value corresponding to hmin
(zmax   =                12500) Pixel value corresponding to hmax
(fd     =                     )
(mode   =                   ql)

Exit with ":go" to execute. You'll need to do the same thing for your SFlatz. You now have the final super-sky flat for your i and z data.

Applying your Sky Flats

Now we are ready to apply our super-sky flats to all of our data (targets, astrometric fields, standards, etc.). Recall that the g and r data files all have names like lfc123.fits, but the i and z data have been fringe corrected and are named lfc123f.fits. Make a copy of objects.list called tosflat.list, and edit the filenames to be the correct input files. We also want new output filenames to make sure we are applying a good correction. We create this list and run ccdproc one more time:

ms> !awk '{print $1"s"}' objects.list > sflatout.list
ms> epar ccdproc

PACKAGE = mscred
   TASK = ccdproc
    
images  =        @tosflat.list  List of Mosaic CCD images to process
(output =       @sflatout.list) List of output processed images
(bpmasks=                     ) List of output bad pixel masks
(ccdtype=                     ) CCD image type to process
(noproc =                   no) List processing steps only?

(xtalkco=                  yes) Apply crosstalk correction?
(fixpix =                  yes) Apply bad pixel mask correction?
(oversca=                  yes) Apply overscan strip correction?
(trim   =                  yes) Trim the image?
(zerocor=                  yes) Apply zero level correction?
(darkcor=                   no) Apply dark count correction?
(flatcor=                  yes) Apply flat field correction?
(sflatco=                  yes) Apply sky flat field correction?
(split  =                   no) Use split images during processing?
(merge  =                   no) Merge amplifiers from same CCD?
(xtalkfi=            !xtalkfil) Crosstalk file
(fixfile=                  BPM) List of bad pixel masks
(saturat=            !saturate) Saturated pixel threshold
(sgrow  =                    1) Saturated pixel grow radius
(bleed  =            !bleedval) Bleed pixel threshold
(btrail =                   15) Bleed trail minimum length
(bgrow  =                    1) Bleed pixel grow radius
(biassec=             !biassec) Overscan strip image section
(trimsec=             !trimsec) Trim data section
(zero   =              medZero) List of zero level calibration images
(dark   =                     ) List of dark count calibration images
(flat   =          DFlat*.fits) List of flat field images
(sflat  =       MedSFlat*.fits) List of secondary flat field images
(minrepl=                  0.9) Minimum flat field value

(interac=                   no) Fit overscan interactively?
(functio=             legendre) Fitting function
(order  =                    1) Number of polynomial terms or spline pieces
(sample =                    *) Sample points to fit
(naverag=                    1) Number of sample points to combine
(niterat=                    1) Number of rejection iterations
(low_rej=                   3.) Low sigma rejection factor
(high_re=                   3.) High sigma rejection factor
(grow   =                   0.) Rejection growing radius
(fd     =                     )
(fd2    =                     )
(mode   =                   ql)

Exit with ":go" to execute. Examine some of the output images for each filter to make sure they are ok.

Improving the World Coordinate System

Due to the significant distortions in wide field imaging cameras, it is necessary to derive and apply corrections to the astrometric solutions. We begin by using the default solutions for the LFC. There exist versions for the I, R_s, and g' filters, located in the mscdb$palomar/lfc directory. If you're data are taken in a filter for which there is no default solution, use the one from the most similar filter. This means you may have to make separate file lists for different filters, and run mscsetwcs for each set.

ms> epar mscsetwcs
PACKAGE = mscred
   TASK = mscsetwcs

images  =       @sflatout.list  Mosaic images
database= mscdb$palomar/lfc/lfc-I.db  WCS database   <=== USE APPROPRIATE ONE FOR YOUR DATA!
(ra     =                   ra) Right ascension keyword (hours)
(dec    =                  dec) Declination keyword (degrees)
(equinox=              equinox) Epoch keyword (years)
(ra_offs=                   0.) RA offset (arcsec)
(dec_off=                   0.) Dec offset (arcsec)
(extlist=                     )
(mode   =                   ql)

Exit with ":go" to execute.
Next, we need to match catalogs of USNO stars to each of our fields. We first edit the parameters for the mscgetcatalog task, which will go out on the web and get the necessary data:

ms> epar mscgetcat
PACKAGE = mscred
   TASK = mscgetcatalog

input   =                       List of Mosaic files
output  =                       Output file of sources
(magmin =                  14.) Minimum magnitude
(magmax =                 19.5) Maximum magnitude
(catalog=         CADC:USNO-A2) Catalog
(rmin   =                  21.) Minimum radius (arcmin)
(mode   =                   ql)

Exit with ":wq" to save. Now, we run msccmatch to do the matching. Make an input list of the sky-flattened long exposures called targets.list, and then:
NOTE: for standard stars or shallow exposures, you should set the magnitude limits in mscgetcat to something like 10...16, and the parameters in msccmatch as
(nsearch= 10) Maximum number of positions to use in search
(search = 10.) Translation search radius (arcsec)

ms> epar msccmatch
PACKAGE = mscred
   TASK = msccmatch
    
input   =        @targets.list  List of input mosaic exposures
coords  =     !mscgetcat $I $C  Coordinate file (ra/dec)
(outcoor=                     ) List of updated coordinate files
(usebpm =                  yes) Use bad pixel masks?
(verbose=                  yes) Verbose?

                                # Coarse Search
(nsearch=                   30) Maximum number of positions to use in search
(search =                  30.) Translation search radius (arcsec)
(rsearch=                  0.2) Rotation search radius (deg)

                                # Fine Centroiding
(cbox   =                   11) Centering box (pixels)
(maxshif=                  3.5) Maximum centering shift to accept (arcsec)
(csig   =                  0.1) Maximum centering uncertainty to accept (arcsec)
(cfrac  =                  0.5) Minimum fraction of accepted centers
(listcoo=                  yes) List centered coordinates in verbose mode?

                                # WCS Fitting
(nfit   =                   40) Min for fit (>0) or max not found (<=0)
(rms    =                  0.5) Maximum fit RMS to accept (arcsec)
(fitgeom=              general) Fitting geometry
(reject =                   3.) Fitting rejection limit (sigma)
(update =                  yes) Update coordinate systems?
(interac=                  yes) Interactive?
(fit    =                  yes) Interactive fitting?
(graphic=             stdgraph) Graphics device
(cursor =                     ) Graphics cursor

accept  =                  yes  Accept solution?
(mode   =                   ql)

Exit with ":go" to execute. Some information about the matching will spew to the screen, and then a window showing fit info will appear. Hitting "r" will show x residuals as a function of x positions. Hitting "y" shows y residuals as a function of y position, and "s" shows y residuals as a function of x position. So, something has gone wrong with part of the fit. Looking at the plot you get when you hit "s", we will delete the minority of points that are in the upper left corner. Move your cursor to that area, and press "d". This deletes the point nearest the cursor (it changes to an "X"). Keep holding the "d" key until all the points in this region are deleted. Now, hit "f" to refit. While the fit isn't perfect, it will do as a first pass. Hit "q" to quit. In your xgterm window, you will be asked

Accept solution? (yes):

Just hit return to answer yes. You should get a message

Coordinate system updated.

You'll go through all the target exposures, fixing this erroneous matching.

Now, we will rerun msccmatch again; with our first solution in place, things should go much better.

ms> msccmatch
List of input mosaic exposures (lfc262s): @targets.list 
Coordinate file (ra/dec) (!mscgetcat $I $C):

As before, stuff spews to the screen, and the graphic window pops up. When you hit the "x", "y", "r", and "s" keys, things should look much better. An example is shown here. You can adjust some of the fit parameters to improve the residuals (shown on the upper left of the graphic window). Type ":reject 2.5" to set the rejection threshold to 2.5 sigma. Type ":maxiter 3" to set the number of rejection iterations to 3. Type ":order 5" to set the fitting order to 5. Typing any of these commands without the argument will show you the current value of that parameter. Finally, hit "f" to refit. Rejected points will be marked with a circle. You can delete any points you want by placing the cursor on the point and hitting "d". Hit "u" to undelete. Remember to hit "f" to refit. Once the residuals are less than about 0.3-0.4 arcsec, you are good to go. Hit "q" to quit. Go through this procedure for all of the target (long exposure) images.

If you have astrometry problems then maybe the telescope pointing is off. Read the description of how to fix this.

Cosmic Ray Cleaning

We next need to create a mask for the cosmic rays in the target exposures. These pixels should be excluded from the image combine operations in future steps. We will use the iraf task craverage, located in the crutil package. Unfortunately, this task doesn't work on MEF images, so we'll make a little script that calls it for each CCD in the mosaic.

In a terminal window, outside of iraf:

% rm foo gocraver.cl
% sed 's/lfc//g' targets.list | sed 's/s//g'> nums.list
% echo set masktype=pl > gocraver.cl
% foreach i (`cat nums.list`)
foreach? foreach j ( 1 2 3 4 5 6)
foreach? echo $i $j > foo
foreach? awk '{print "craverage lfc"$1"s["$2"] output=\"\" crmask=\"crmask_"$1"_"$2"\""}' foo >> gocraver.cl
foreach? end
foreach? end

This will produce a file with lines like
craverage lfc261[1] output="" crmask="crmask_lfc261_1"

Next, we'll set up the parameters for the craver task as follows:

ms> crutil
ms> epar craverage

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   =                   13) Block average box size
(nrej   =                   11) Number of high pixels to reject from the average
(nbkg   =                   15) 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 =                  4.0) High cosmic ray sigma outside object
(crgrow =                   1.) 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)

Exit with ":wq" to save changes. We execute our script with

cr> cl < gocraver.cl

This task is very slow, and it may take many hours (or even days!) for the entire script to run.

Now, we need to combine our existing bad pixel masks with the newly generated cosmic ray pixel list. Make sure you have crplusbpmask.cl in your scripts directory, and the task definition

task crplusbpmask="scripts$crplusbpmask.cl"

in your loginuser.cl file. This task will take the existing bad pixel masks (bpm_###/bpm_lfc[0-5]) and combine them with the cosmic ray masks (crmask_lfc###_[1-6]). The output files will be called bpm_###/bpmcr_lfc[0-5], and the headers of the sky-flatfielded images updated. NOTE: image names, paths to bpms, etc. are all hard coded in this script, so if you are not following my naming conventions, you'll have to edit the script accordingly! To run the script, just execute

ms> crplusbpmask @nums.list

NOTE: The interactive checks of cosmic ray detection below are NOT necessary, especially if you have a large (say, >8) number of pointings in a given filter on a field. If so, you can skip to fixing the bad pixels.

Examine the masks by displaying the images with the new bad pixel masks overlaid. I like to do this on a chip-by-chip basis:

ms> displ lfc261s[1] bpm=BPM bpdisp=overlay bpcolor=red frame=1

You should see cosmic rays, bad columns, saturated star centers, and bleed trails all masked. You will see that sometimes parts of large cosmic ray strikes or other defects are not masked. We need to add these to the bad pixel file. To do so, we will use imexam to mark the bad regions.

cl> epar imexam

PACKAGE = tv
   TASK = imexamine
    
input   =                       images to be examined
(output =                     ) output root image name
(ncoutpu=                   10) Number of columns in image output
(nloutpu=                   10) Number of lines in image output
frame   =                    1  display frame
image   =                       image name
(logfile=           badpix.log) logfile
(keeplog=                  yes) log output results
(defkey =                    a) default key for cursor list input
(autored=                  yes) automatically redraw graph
(allfram=                   no) use all frames for displaying new images
(nframes=                    0) number of display frames (0 to autosense)
(ncstat =                    5) number of columns for statistics
(nlstat =                    5) number of lines for statistics
(graphcu=                     ) graphics cursor input
(imagecu=                     ) image display cursor input
(wcs    =              logical) Coordinate system
(xformat=                     ) X axis coordinate format
(yformat=                     ) Y axis coordinate format
(graphic=             stdgraph) graphics device
(display= display(image='$1',frame=$2)) display command template
(use_dis=                  yes) enable direct display interaction
(mode   =                   ql)

Exit with :wq to save parameters. When using imexam the b key allows us to mark the bottom left and top right corners of a region, which are saved in the logfile we specified in the epar (badpix.log). Note that this file is always appended to, so we need to delete it each time we go on to a new image. So, after you've displayed the image:

cl> !rm badpix.log
cl> imexam
...move around the image systematically, searching for missed bad regions
...use the b key to mark the bottom left and top right corners of bad area
...when done, hit "q" to quit

If any bad pixels are missed, we'll use imreplace to edit the pixel mask, making sure to use the correct chip number (lfc0-5). REMEMBER, if you are displaying lfc261[1], the corresponding mask is lfc0! When I am done with imexam on an image, I use awk to create an input file for imreplace, using the badpix.log we just made. An example awk command for lfc261.fits[1] is

% awk '{print "bpm_261/bpmcr_lfc0.pl["$1":"$2","$3":"$4"]"}' badpix.log > gofix

Then I can use imreplace in iraf to edit the pixel mask:

ms> imutil
im> epar imreplace

PACKAGE = imutil
   TASK = imreplace

images  =               @gofix  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)

Remember that you need to display each subimage and run imexam/imreplace for each subimage of each image!

When you are happy with your bad pixel masks, prepare to run fixpix. First, we prepare a list of individual CCD frames to run it on, using a terminal session:

% rm fixpix.list
% foreach i ( `cat nums.list` )
foreach? foreach j ( 1 2 3 4 5 6)
foreach? echo lfc$i\s\.fits\[$j\] >> fixpix.list
foreach? end
foreach? end

Your resulting list should look like

lfc261s.fits[1]
lfc261s.fits[2]
lfc261s.fits[3]
lfc261s.fits[4]
lfc261s.fits[5]
lfc261s.fits[6]
lfc262s.fits[1]
etc.

Then, in iraf, we run fixpix, which requires the loading of the crutil package:

im> crutil
cr> epar fixpix

PACKAGE = proto
   TASK = fixpix

images  =         @fixpix.list  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)

Exit with ":go" to execute.

Generating the Tangent-Plane Projected Images

Next, we want to create single images from our MEF files, that have all been tangent-plane projected, with the same pixel scale and center position. We will use one of the chips of the center pointing as our reference. Typically, the center pointing will be the first exposure in a sequence taken in the same filter; check the observing log to see which one it is. We will use the task mscimage to perform the projection. We also want to make sure that the proper mean background is used for areas in each image with no data. This is done using iterstat.

At this point, we will be working on sets of images taken of the same target AND in the same filter. Create a list of image numbers for each such data set. For instance, list all image numbers for all SC1604 superclusterpointing 1 images in sc1604_1.list. Then, in a regular shell session:

% awk '{print "iterstat lfc"$1"s.fits[1][700:1500,1500:3500]"}' sc1604.list > goiter.1604

Now, inside your IRAF session, epar iterstat. If it is not part of your user package, you can load it with:

ms> stsdas
st> hst_calib
hs> nicmos

So, now let's set it up:

ms> epar iterstat
PACKAGE = user
   TASK = iterstat

image   =                       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)

Exit with ":wq" to save changes. Now, we will run our little script goiter, and grab the median background values:

ms> cl < goiter > iterout.1604
lfc264s.fits[1][700:1500,1500:3500]: mean=2933.466  rms=40.76279  npix=1593639  median=2930.953  mode=2930.882
etc.

At this point your output file should have lines like

lfc039s.fits[1][700:1500,1500:3500]: mean=4232.546  rms=47.98962  npix=1588082  median=4229.516  mode=4228.3
38

Now, in your shell session, munge this and the numbers file to create a boundary extension list:

% sed 's/=/ /g' iterout.1604 | awk '{print $9}' > foo
% paste sc1604.list foo > boundaries.sc1604

We now have what we need to create the script to run mscimage. Make sure you set the proper parameters in the mscimage task.
IMPORTANT!! You must use the same fiducial (reference) ccd frame for ALL observations of a given target. This means using one image as the reference for all observations, over all nights, in all filters, of that target. To do this you will have to decide on the reference image, and copy it into each of the directories for the separate nights.

% awk '{print "mscimage lfc"$1"s.fits mos"$1".fits boundary=\"constant\" constant="$2}' boundaries_sc1604i.list > gomos

Back to IRAF, edit the mscimage parameters:

ms> set masktype=pl
ms> epar mscimage
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
(wcssour=                image) Output WCS source (image|parameters)
(referen=     lfc279s.fits[1]) Reference image
(ra     =		 INDEF) RA of tangent point (hours)
(dec    =		 INDEF) DEC of tangent point (degrees)
(scale  =                INDEF) Scale (arcsec/pixel)
(rotatio=                INDEF) Rotation of DEC from N to E (degrees)

                                # Resampling parmeters
(blank  =                   0.) Blank value
(interpo=               sinc17) Interpolant for data
(minterp=               linear) Interpolant for mask
(boundar=             constant) Boundary extension
(constan=                   0.) Constant boundary extension value
(fluxcon=                   no) Preserve flux per unit area?
(ntrim  =                    7) Edge trim in each extension
(nxblock=                 2100) X dimension of working block size in pixels
(nyblock=                 4200) 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_ext =                     )
(fd_coor=                     )
(mode   =                   ql)

Exit with ":wq" to save changes. Make sure you have set the values for interpolant and boundary correctly! Now, we run our little script:

ms> cl < gomos

Fix Satellite Trails

At this point, we have a single mosaicked image for each of our long-exposure targets. We need to remove any satellite trails before we stack them. We use a perl script called sat-b-gon written by Matt Hunt. Put this either in your local directory or someplace in your path.

Display each mosaicked image and examine them for satellite trails:

ms> mscdispl mos279 1

If there is a satellite trail, record its xmin ymin xmax ymax in a file called something like mos279.satin. Also, note how wide the satellite trail is; you should adjust the image stretch so you can see the faintest parts of the trail. I typically make the grow radius one or two pixels bigger than I think it should be just to be sure. Then, in a shell, execute:

% sat-b-gon -g width mos279.satin > mos279.satout

where width is the half-width (grow radius) of the trail. The output file will have a list of x,y that are in the trail. Now, we need to convert this to a pixel list (.pl file) and add these bad pixels to our mask for the mosaic. We need to get the dimension of each image with a satellite trail. In iraf:

cl> imhead mos279
mos279[8288,8237][real]: sc1604_1
cl> set masktype=pl
cl> epar text2mask
PACKAGE = proto
   TASK = text2mask

text    =        mos279.satout  Text file description of mask regions
mask    =         mos279sat.pl  Mask image name
ncols   =                 8288  Number of columns in mask
nlines  =                 8237  Number of lines in mask
(linterp=                    1) Mask code for rectangles narrower along lines
(cinterp=                    2) Mask code for rectangles narrower along columns
(square =                    3) Mask code for squares
(pixel  =                    4) Mask code for single pixels
(mode   =                   ql)

Exit with ":go" to execute.

Now we combine the bad pixel mask for the mosaic with the satellite trail mask using imexpr:

cl> imutil
cl> set masktype=pl
cl> epar imexpr
PACKAGE = imutil
   TASK = imexpr
    
expr    =             max(a,b)  expression
output  =       combmask279.pl  output image
(dims   =                 auto) output image dimensions
(intype =                 auto) minimum type for input operands
(outtype=                 auto) output image pixel datatype
(refim  =                 auto) reference image for wcs etc
(bwidth =                    0) boundary extension width
(btype  =              nearest) boundary extension type
(bpixval=                   0.) boundary pixel value
(rangech=                  yes) perform range checking
(verbose=                  yes) print informative messages
(exprdb =                 none) expression database
(lastout=                     ) last output image  (don't worry about this one)
a       =         mos279sat.pl  operand a
b       =       mos279_bpm.pl  operand b
c       =                       operand c
d       =                       operand d
....etc.

Exit with ":go" to execute. Now, we should check this combined mask to make sure we masked the entire trail properly:

cl> mscdispl mos279 bpd=overlay bpm=combmask279.pl

Zoom in on the satellite trail and adjust the contrast. If it is good, we can rename the combined mask to the mosaic mask

cl> mv combmask279.pl mos279_bpm.pl

If the mask is no good, rerun sat-b-gon with a different width or endpoints, rerun text2mask and imexpr. You'll need to delete the intermediate .pl files (mos279sat.pl, combmask279.pl) each time.

Of course, if there are no trails in an image, you don't need to do the above steps.

Subtract Global Sky Pattern

Despite our best attempts to flatfield our data, there are typically very large scale background gradients in the sky level across the entire mosaic. You can see this by displaying a mosaicked long exposure image and adjusting the contrast so the background becomes bright. To remove these gradients, we fit a 2-dimensional surface to the background level and subtract it using mscskysub.

cl> !ls -1 mos*fits > skysubin
cl> !awk '{print "mos"$1"s"}' nums.list > skysubout
cl> epar mscskysub
PACKAGE = mscred
   TASK = mscskysub
    
input   =            @skysubin   Input images to be fit
output  =           @skysubout  Output images
xorder  =                    2  Order of function in x
yorder  =                    2  Order of function in y
(type_ou=             residual) Type of output (fit,residual,response,clean)
(functio=                  leg) Function to be fit (legendre,chebyshev,spline3)
(cross_t=                  yes) Include cross-terms for polynomials?
(xmedian=                  200) X length of median box
(ymedian=                  200) Y length of median box
(median_=                  25.) Minimum fraction of pixels in median box
(lower  =                   3.) Lower limit for residuals
(upper  =                   2.) Upper limit for residuals
(ngrow  =                    0) Radius of region growing circle
(niter  =                    2) Maximum number of rejection cycles
(regions=                 mask) Good regions (all,rows,columns,border,sections,c
(rows   =                    *) Rows to be fit
(columns=                    *) Columns to be fit
(border =                   50) Width of border to be fit
(section=                     ) File name for sections list
(circle =                     ) Circle specifications
(mask   =                  BPM) Mask
(div_min=                INDEF) Division minimum for response output
(mode   =                   ql)

Exit with ":go" to execute. Check the output; you should see that the background gradients are reduced. I've not been able to eliminate them completely, but they seem to be well below the 1% level.

Stack the Mosaicked Images

We're now ready to combine all the images of a single field in a given filter into a single, deep image. To do this, gather all of the sky-subtracted mosaicked images you would like to stack into a single directory, if they aren't already. You will need the actual images (mos_###s.fits) and the corresponding bad pixel masks (mos_###_bpm.pl). For each filter, make a file listing the images you want to combine, called something like sc1604_i.list (or whatever your object/filter are). We'll use mscstack to do the combining:

ms> epar mscstack
PACKAGE = mscred
   TASK = mscstack
    
input   =       @sc1604_i.list  List of images to combine
output  =             sc1604_i  Output image
(headers=                     ) List of header files (optional)
(bpmasks=         sc1604_i_bpm) List of bad pixel masks (optional)
(rejmask=                     ) List of rejection masks (optional)
(nrejmas=                     ) List of number rejected masks (optional)
(expmask=         sc1604_i_exp) List of exposure masks (optional)
(sigmas =                     ) List of sigma images (optional)

(combine=              average) Type of combine operation (median|average)
(reject =            avsigclip) Type of rejection
(masktyp=            goodvalue) Mask type
(maskval=                   0.) Mask value
(blank  =                   0.) Value if there are no pixels

(scale  =                 none) Image scaling
(zero   =               median) Image zero point offset
(weight =                 none) Image weights
(statsec= [3000:6000,3000:6000]) Image section for computing statistics

(lthresh=                 100.) Lower threshold
(hthresh=               60000.) 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 =                   4.) Lower sigma clipping factor
(hsigma =                  3.5) 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)

Depending on your application, you may also want to coadd ALL of your data for a single field (but different filters) to make a super-deep detection image. IMPORTANT: Note that the output images in the different filters (and the super-deep image) will NOT have objects in the same x,y locations because the dithers used during observations are typically not identical! This means that if you want to run detection on the deep image, and use the same apertures for the individual filters, you will need to calculate and apply offsets to align everything in x and y! This is a big (but necessary) pain in the butt.

IMPORTANT
For some reason the EQUINOX header keyword is not set properly in the image headers. This is needed for SExtractor to output correct coordinates. Just do:

cl> hedit sc1604_i.fits EQUINOX 2000.0 add+ ver-

and similarly for all the mosaicked images.

Prepare Standards for Photometry

In order to run a photometry package such as SExtractor on any standard star observations, they also must be mosaicked, and a weight map (for blank pixels) produced. Make a list that contains the file numbers of all of the observed standards, say stdnums.list. Then we proceed as we did for the target frames:

% cp stdnums.list boundaries_stds.list
% awk '{print "iterstat lfc"$1"s.fits[1][700:1500,1500:3500]"}' stdnums.list > goiter
% emacs boundaries_stds.list &

Now, inside your IRAF session:

ms> epar iterstat
PACKAGE = user
   TASK = iterstat

image   =                       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)

Exit with ":wq" to save changes. Now, we will run our little script goiter, and put the output median values into the boundaries file that is open for editing:

ms> cl < goiter
lfc703s.fits[1][700:1500,1500:3500]: mean=49.30294  rms=8.8368  npix=1601215  median=48.38079  mode=47.52637
etc.

At this point your boundaries file should have lines like

703 48
704 61
etc.

Now, in your shell session, create the script to run mscimage. We need to do this once for EACH of the standards. For each standard, copy the appropriate lines for just that standard from the boundaries file into a distinct file (say, sa108.list). Then, use awk to make a script to run mscimage for that standard:

% awk '{print "mscimage lfc"$1"s.fits mos"$1".fits boundary=\"constant\" constant="$2}' sa108.list > gomos_sa108.cl

Next, set the proper parameters in the mscimage task. You should use some fiducial ccd frame for your target (set in the reference keyword). I typically use the first observation of the standard I am working on. We must change this reference image for each standard!

% awk '{print "mscimage lfc"$1"s.fits mos"$1".fits boundary=\"constant\" constant="$2}' boundaries_sc1604i.list > gomos

Back to IRAF, edit the mscimage parameters:

ms> set masktype=pl
ms> epar mscimage
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
(wcssour=                image) Output WCS source (image|parameters)
(referen=     lfc703s.fits[1]) Reference image  <====Make sure this is correct
(ra     =		 INDEF) RA of tangent point (hours)
(dec    =		 INDEF) DEC of tangent point (degrees)
(scale  =                INDEF) Scale (arcsec/pixel)
(rotatio=                INDEF) Rotation of DEC from N to E (degrees)

                                # Resampling parmeters
(blank  =                   0.) Blank value
(interpo=               sinc17) Interpolant for data
(minterp=               linear) Interpolant for mask
(boundar=             constant) Boundary extension
(constan=                   0.) Constant boundary extension value
(fluxcon=                   no) Preserve flux per unit area?
(ntrim  =                    7) Edge trim in each extension
(nxblock=                 2100) X dimension of working block size in pixels
(nyblock=                 4200) 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_ext =                     )
(fd_coor=                     )
(mode   =                   ql)

Exit with ":wq" to save changes. Make sure you have set the values for reference, interpolant and boundary correctly! Now, we run our little script:

ms> cl < gomos_sa108.cl

Repeat the above steps until you are done with all your standards. In addition to the mosaicked images (mos###.fits) you will also have bad pixel masks (mos_###_bpm.pl). For SExtractor to use these as weight maps, these need to be converted to fits images, wherein the good pixels have a value of 1 and the bad pixels have a value of 0. We will make an IRAF script to do this. This is easiest done in from the shell:

% echo stsdas > makeweights.cl
% echo toolbox >> makeweights.cl
% echo imgtools >> makeweights.cl
% awk '{print "imcopy mos"$1"_bpm.pl foo"$1".fits",ORS"imcalc foo"$1".fits weight_"$1".fits \"if im1 .gt. 0. then 0. else 1.\"",ORS"imdel foo"$1".fits"}' stdnums.list >> makeweights.cl

Then we can just run this script within IRAF:

cl> cl < makeweights.cl

I also like to have all of my data in counts/second. This makes calculating and applying the photometric zero points much easier, because SExtractor doesn't know about exposure times. To deal with this, we simply divide each image of a standard by its exposure time (for the standards). Make a copy of the list of the image numbers for the standards, and edit it to include the exposure times:

% cp stdnums.list stdexp.list
% emacs stdexp.list

Your file stdexp.list should have lines like

703 3.
704 2.
etc.

We make a script to normalize each standard by its exposure time, and run the script:

ms> !awk '{print "imarith mos"$1".fits / "$2" mos"$1"_norm.fits"}' stdexp.list > normstd.cl
ms> cl < normstd.cl

At this point, we have everything we need to run SExtractor on our standards.