Stacking Process

After reducing, calibrating, and performing photometry on data, we are finally able to stack. Stacking, in short, consists of aligning and literally stacking images of like fields on top of one another to increase flux and depth. "Home base" is pawnee:/data/stack.

The file combine.doc contains an overview of the steps involved in stacking, for quick reference later. I will explain here, in greater detail, what you can find in combine.doc.

Because of the relative sensitivities of different filters, for any given field a different number of images must be stacked for different filters. The number of fields needed to stack are as follows:
      7598 and 7661: 4 images
      8132 and 8199: 12 images
      8615 and 8685: 20 images
      9155 and 9233: 32 images
Stacking should be done in filter pairs (i.e., 7598 with 7661, 8132 with 8199, and 8615 with 8685). Stacking cannot be performed on a field/filter until enough images have been taken in BOTH filters in a pair. Each stack will consist of images of only one filter, but you must perform the stacking on both filters in a pair one after the other.

Get Airmasses
1. Open get_airmass.cl. Copying lines of this file to iraf will return airmasses for desired images. In order to get airmasses for whatever field/filter pair we are stacking, we first must add the necessary lines into get_airmass.cl. Do so, taking care to follow the order of the file: ELAIS_N1 7598/7661, LOCKMAN 7598/7661, WYSCH 7598/7661; ELAIS_N1 8132/8199, LOCKMAN 8132/8199, WYSCH 8132/8199; ELAIS_N1 8615/8685, LOCKMAN 8615/8685, WYSCH 8615/8685. When the lines are added, copy them into iraf to retrieve the airmasses. NOTE: Simply highlighting all desired lines and pasting them into iraf all at once is the efficient way of doing this, but be careful. Depending on how you highlight the last line in your desired block of lines, you may need to hit Enter once in iraf after pasting to get it to return the last line.
2. Add these airmasses to imagesinfofiles.list. This file follows the same order as get_airmass.cl. Be sure to use proper formatting, and remember to include an "images.info" line before each field/filter block, and make sure there are no spaces after "images,info"; get_scaling will not run if there are spaces. NOTE: You must include an "images.info" line before EACH FIELD AND FILTER, not each filter pair. Make sure you have copied all of the appropriate .mag.1 files.

Scale Images and Create .txt Files
5. In unix, in /stack, run get_scaling.f (a FORTRAN code): compile by typing "f77 get_scaling.f," which creates a file called a.out. This should already have been done; unless you have made changes to the .f file, a.out has already been renamed get_scaling. Run this by typing ./get_scaling. This will create .txt files in the relevant field/filter directories, alongside the .mag.1 files. These .txt files are the coordinates of the 5 stars used in the photometry step.
6. Due to a gain change in the CCD chip used at WIRO Prime occuring around Jan 2007 to Dec 2008 you have to be carefull which images are used for scaling. This will only be a problem where there are images before, after, and perhaps during this gain change in the same field and filter. Therefore, the get_scaling program will output the ratio of every image star counts (those 5 stars which photometry was performed) to the highest image star counts if that ratio is below 95%. This is an arbitrary cut-off. Usually images with that ratio between 90-95% indicate low star counts that are within acceptable limits. However, If the ratio is below 90% then it could indicate acceptable low star counts or an image with a different gain. Check this by opening the images.info file in stack/field/filter. This file has the flux counts for the 5 stars of all images. If your field/filter has normal star counts then go to step 8; if it has gain issues go to the next step.
7. You want to scale all images to the most reliable, high count image. The scaling program has a flag parameter that will not allow that image to be the scale image. This flag is set in the imagesinfofiles.list and can either be a 0 or a 1 at the end of each image. Zero indicates nothing wrong with that particular image, and a ONE indicates that the image should not be used as the highest count image, even if it has the highest counts. The following are some possible permutations of images with different gains:
  - If the majority of images are higher in counts compared to a few low count images, then the scaling program will correctly scale the images.
  - If the majority of images are lower in counts compared to the few high count images, then you must flag the few high count images with a 1. Otherwise, all images will be scaled to the few high count images.
  - If there are close to equal number of high count images and low count images, you must look at the dispersion of star counts. The sample (high or low) with the lowest dispersion will be the most accurate and therefore the most reliable to scale to. If the most reliable is the highest count, then the scaling program will work properly. However, if the low count sample is the most reliable then you must flag the high count sample images with a 1.
8. Open correct_for_airmass_n_clouds.cl. Copy newly created (by get_scaling.f) lines to iraf. Again, make sure you've entered all lines; you may have to press Enter to enter the final line into iraf. This step creates a*'date'A.fit files in the photometry directory which are now scaled.

Align Images
The following steps will align your images around the center image you chose. Do all following steps separately for both filters in a filter pair.

9. Enter into a terminal: ls *.txt > ToCoordinates.list
10. Enter into iraf: geomap @ToCoordinates.list Coordinates.list 1 2037 1 2048. Hit "q" in the Textronic window to cycle through all the images. IMPORTANT!! Watch carefully as the Textronic window cycles through the images. The small plus signs should stay in relatively the same place, and the grid of lines should remain basically horizontal and vertical. If the grid becomes askew (you'll know), or if the plus signs suddenly move around, something's wrong. It usually means that you accidentally hit spacebar on the wrong star during phot, or you spacebar'd the stars in the wrong order. Go back and redo those steps.
11. Enter: ls /data/WYSCH/19/8132/*A.fit > ToAlign.list. Obviously, update this line for whatever field/filter with which you're working.
12. Enter: cp ToAlign.list Aligned.list. Open up Aligned.list and add "R" after each "A" in the file names. This "R" indicates that these images are aligned.
13. Enter into iraf: geotran @ToAlign.list @Aligned.list Coordinates.list @ToCoordinates.list.
14. Go through the newly aligned images by entering "prevu @Aligned.list" into iraf. Press spacebar to cycle through the images. If you place a region circle on a star, you will be able to check if the images are truly aligned. If the star moves around, the alignment failed.
Summary of Lines to Enter
      ls *.txt > ToCoordinates.list
      geomap @ToCoordinates.list Coordinates.list 1 2037 1 2048 (press "q" to cycle)
      ls /data/WYSCH/19/8132/*A.fit > ToAlign.list
      cp ToAlign.list Aligned.list (add "R" to Aligned.list)
      geotran @ToAlign.list @Aligned.list Coordinates.list @ToCoordinates.list
      prevu @Aligned.list (check alignment)

Stack Images
To stack the images, all you need to do is combine them using imcombine, using certain parameters. The most important of these is hsigma. Normally, the default hsigma works fine. However, every once in awhile you'll get a field whose stack looks very streaky or grainy. If this is the case, remake the stack with increasingly higher values of hsigma, beginning at 3 and going up in increments of 0.1. Keep going up until the streaking/graininess disappears. You can be pretty liberal with this; continue going higher until you're absolutely sure all streaking and graininess are gone. NOTE: There was one case I ran into where, instead of increasing hsigma, I had to increase lsigma. This should rarely, if ever, happen, but if higher and higher hsigmas don't seem to be working, perhaps consider trying to tweak lsigma. Also, once you reach a good hsigma value the image will not change significatly.

If hsigma (or lsigma) does indeed need to be changed, note the value you used in combine.doc. Here's the imcombine command to enter into iraf (leave out the hsigma parameter at first):
13. imcombine @Aligned.list stack.fit combine=average reject=avsigclip weight=@weighting.bin lthresh=INDEF hthresh=INDEF hsigma=3.2
14. Repeat the above alignment and stacking steps for the other filter in the filter pair.

Add WCS Coordinates to Header
This final task adds WCS coordinates to the headers of the images. Unlike alignment and stacking, you only need to do this set of steps once (i.e., do not repeat for both filters in a filter pair). However, begin in the image's lower-filter directory (i.e., wysch/19/8132 or wysch/19/8615).

15. Load the "gasp" package in iraf. This can be found in "stsdas" followed by "analysis."
16. Open all .txt files (i.e., n *.txt &). Clicking through these, you'll notice that the first two columns are always the same, while the second two columns differ for each file. These columns are coordinates of the stars you chose for photometry: the first two columns are the coordinates of the center image, and the second two columns are the coordinates of each specific image. Keep these open; we will use the first two columns presently.
17. Load stack.fit into frame 1 of ds9 (i.e., disp stack.fit 1), and load the corresponding DSS image in frame 2 (i.e., disp ../../../fields/dss/w19[0] 2). Remember to use [0] when displaying the DSS images. NOTE: Your stack.fit should be xy-flipped, but the DSS image should not be. Fiddle with the zoom on each image until you can blink the images and be able to map, by eye, a star from your stack to the corresponding star in the DSS image.
18. You now need to examine the DSS image and mark the stars that correspond to the stars used in photometry. To do this, type "imexam logfile=dss.txt keep+". This prepares a file called "dss.txt" that will record each keystroke within imexam. Look at the first two columns of a .txt file (they should still be open). For each pair of coordinates, find the star to which they point in your stack image, press "tab" to move to the DSS image, and press "a" on the corresponding DSS image star. Do this for all five stars, then quit out of imexam. NOTE: Because the DSS image is somewhat offset from your stacked image, if one of your stars is somewhat close to the edge, it will be out of the frame in the DSS image. You DO NOT need to go back and do everything. This step merely lines up the two images, and using the five stars you used before is done strictly out of convenience. So just choose another star. But REMEMBER which star it is, because in step 20 below, you will have to enter its coordinates (gotten by using phot on the star IN THE IMAGE YOU USED AS THE CENTER) in a file.
19. Enter the following into iraf, remembering to update the DSS image. This will create a file called dss.coords:
      xyeq iminfo=yes image=../../../fields/dss/w19[0] xy=dss.txt pix_center=iraf xcol=1 ycol=2 nskip=2 orig=yes new=no ra_h=yes ra_f=%h dec_f=%h > dss.coords
20. You must now create a new file, called ccmap.in. Into this file you will paste four columns: the first two columns will be the first two columns of any .txt file; and the second two columns will be the RA and DEC columns of dss.coords. No columns headers are needed; only five lines of text, one per star. After doing this, you can close dss.coords, ccmap.in, and all .txt files. NOTE: If you had to use a different star in step 18, use phot on your center image to get its coordinates, and put these in the proper place in your ccmap.in table. Note that you only need to replace the X and Y coordinates of the star; you can still use the RA and DEC in dss.coords. (For example, if your second star was out of frame and you replaced it, the first two columns in the second row will be the X and Y coordinates of your new star.)
21. Last step. Enter the first line below into iraf, CD INTO THE OTHER FILTER'S DIRECTORY (i.e. wysch/19/8199 or wysch/19/8685), and enter the second line below. IMPORTANT!! Each of those lines will display some text in the iraf window. Towards the bottom of this text, in a section entitled "Coordinate mapping parameters," there will be a line that enumerates the pixel scale ("X and Y scale"). MAKE SURE IT'S CLOSE TO 0.523"/pix. It doesn't need to be exact, but if it's off by, say, 0.004 or more, it should warrant examination of your steps to see if anything went wrong.
Enter the following (remember to update the filter in the second line):
      ccmap input=ccmap.in database=ccmap.database images=stack.fit lngunits=hours latunits=degrees results=ccmap.out update=yes interactive=no
      cd into other filter's directory
      ccmap input=../8132/ccmap.in database=ccmap.database images=stack.fit lngunits=hours latunits=degrees results=ccmap.out update=yes interactive=no Also, if you are redoing a stack you must redo the ccmap step; it adds information to the header.
22. Ok, I lied, one more step. But this one's trivial. Update the reduction status page by putting "St" in the appropriate boxes of the tables, to indicate that those fields have been stacked. These fields are now complete, as indicated by OTZFQISPSt.
Summary of Lines to Enter
      load stsdas, analysis, gasp
      n *.txt &
      disp stack.fit 1
      disp ../../../fields/dss/w19[0] 2
      imexam logfile = dss.txt keep+ (press "a" on corresponding DSS stars)
      xyeq iminfo=yes image=../../../fields/dss/w19[0] xy=dss.txt pix_center=iraf xcol=1 ycol=2 nskip=2 orig=yes new=no ra_h=yes ra_f=%h dec_f=%h > dss.coords
      n ccmap.in & (paste first two columns of .txt file, and RA and DEC columns from dss.coords)
      ccmap input=ccmap.in database=ccmap.database images=stack.fit lngunits=hours latunits=degrees results=ccmap.out update=yes interactive=no (check pixel scale, cd to other filter)
      ccmap input=../8132/ccmap.in database=ccmap.database images=stack.fit lngunits=hours latunits=degrees results=ccmap.out update=yes interactive=no (check pixel scale)