Guide to Making Custom Stacks from CP Reprojected Data Using IRAF

F. Valdes, January 2019

The NOAO DECam Community Pipeline (CP) data products in the NOAO Science Archive include reprojected and stacked versions. This guide provides a do it yourself recipe for making custom stacks from these reprojected data using IRAF. To use non-IRAF tools this guide can provide some hints as to what is involved. To go from the unreprojected data products, called InstCal and Skysub in the archive, to stacks involves an additional step to reproject or "warp" the data to a common pixel binning. How to do this is indicated in the section Projection Matching.

As as an aside, an alternative to this DIY is to to request custom stacks be produced by the CP. To explore this send email to valdes@noao.edu.

IRAF

This guide is not intended to teach you how to use IRAF. Some guides that can help are:

These are fairly old and don't cover some of the newer features of IRAF such as the extended string manipulation functions used here and documented with:

ecl> phelp strings
Due to the size of the DECam data a 64-bit version of IRAF is required (V2.15 or V2.16). The instructions and examples assume the default parameters. Because there are typically long lists and various string manipulations needed the examples use the scripting capability of the IRAF command language. They show steps done directly on the command line but it is generally easier to put the steps in a file and execute them with:
ecl> cl < myfile.cl

Remember that most IRAF tasks have on-line help. Since the primary IRAF task used here is imcombine you may want to read the help for this complex task with many options using:

ecl> phelp imcombine
Note you are also free to do the various manipulations for making lists or simple scripts that are a list of the commands with your favorite tool such as awk, python, sh, etc.

The Reprojected Data

The reprojected data to be stacked is first downloaded from the NOAO archive. In the archive these are called Resampled data products. They have filenames with _op? in the names. Note that in most cases these versions have had a sky fitting and normalization with a scale length of a fraction of a CCD. This is the same as the unresampled Skysub data product (oki).

The CP reprojects DECam exposures to a standard orientation and pixel scale. The important factor for stacking is that they have the same tangent point. The CP chooses a tangent point from a grid on the sky based on the set of data being processed. This will often be what is needed. However, if the data to be stacked is not on a common tangent projection then one must reproject the InstCal/Skysub/Resampled data which is described in the Projection Matchig section.

In this guide we will demonstrate the steps using the following public files from the NOAO archive (if you want to follow along):

c4d_160704_052117_opd_i_v1.fits.fz  c4d_160704_053537_opd_i_v1.fits.fz
c4d_160704_052117_opi_i_v1.fits.fz  c4d_160704_053537_opi_i_v1.fits.fz
c4d_160704_052828_opd_i_v1.fits.fz  c4d_160704_054247_opd_i_v1.fits.fz
c4d_160704_052828_opi_i_v1.fits.fz  c4d_160704_054247_opi_i_v1.fits.fz

The weight maps (opw) files are not used.

The reprojection geometry is given by the FITS WCS in each CCD extension. The tangent point is given by the CRVAL1/CRVAL2 keywords. There are a variety of coordinates in the header. Below we show the basic RA and DEC, which is the refinement of the telescope pointing after the astrometric solution. In the example below note that the CRVAL values are the same while the pointings move due to dithering or pointing errors. The values are the same in all extensions so we show just looking at the first image extension. In this example these exposures can be stacked without any further reprojection.

ecl> hselect *opi*[1] $I,CRVAL*,RA,DEC yes
c4d_160704_052117_opi_i_v1.fits[1] 296.016 -14.7859 19:44:56.16 -14:47:02.2
c4d_160704_052828_opi_i_v1.fits[1] 296.016 -14.7859 19:44:50.36 -14:46:59.2
c4d_160704_053537_opi_i_v1.fits[1] 296.016 -14.7859 19:44:50.32 -14:45:43.0
c4d_160704_054247_opi_i_v1.fits[1] 296.016 -14.7859 19:44:55.77 -14:45:44.4

Unpacking the Data

For the opi files from the archive they need to be uncompressed using fpack. To avoid problems use V1.7.0 or later. You can either keep the compressed version around or uncompress "in place" with the -D options:

funpack -D c4d_160704_052117_opi_i_v1.fit.fz

The data quality maps (in IRAF also know as bad pixel masks) are compressed in a format that IRAF understands directly. So rather than using fpack simple rename with the ".fz" removed:

mv c4d_160704_052117_opi_i_v1.fit.fz c4d_160704_052117_opd_i_v1.fits

Stacking

Stacking consists of taking a list of FITS images, all with the same projection though not necessarily the same origin, and combining overlapping pixels in some desired way: average, median, or sum. Because the data have already been astrometrically calibrated and the stacking software handles different, integer pixel offsets and padding the result to the full data there is no need to worry about registration. There are two ways to create stacked images from the reprojected multi-extension FITS (MEF) files. One could go directly from a list of the individual CCD extensions to a stacked image. However, because DECam has so many CCDs and, if the stack consists of many exposures, the list would be very long and strain the operating system file descriptors. The other way, which is how the CP does it, is to first "stack" the CCDs in each exposure to make a single image of the exposure and then to stack those exposure images. This guide shows this second method though the first method only differs in the input list.

Converting an MEF to a Simple Image

First we need a list of the extensions referenced as single images. The IRAF task to use is imextensions. The loop example could be put into a file to execute.

ecl> files *fits > in.list
ecl> list = "in.list"
ecl> while (fscan (list, s1) != EOF) {
>>>    s2 = substr (s1, 1, strstr(".fits",s1)) // "list"
>>>    imexten (s1, > s2)
>>>  }
ecl> list = ""
ecl> head c4d_160704_052117_opd_i_v1.list nl=4
c4d_160704_052117_opi_i_v1.fits[1]
c4d_160704_052117_opi_i_v1.fits[2]
c4d_160704_052117_opi_i_v1.fits[3]
c4d_160704_052117_opi_i_v1.fits[4]
ecl> count *opi*list
     60      60    2151 c4d_160704_052117_opi_i_v1.list
     60      60    2151 c4d_160704_052828_opi_i_v1.list
     60      60    2151 c4d_160704_053537_opi_i_v1.list
     60      60    2151 c4d_160704_054247_opi_i_v1.list
    240     240    8604 Total

Unfortunately, the bad pixel masks referenced in the headers of the image data use names from the pipeline and not the archive. So the headers have to be updated.

ecl> concat c4d*list > ccds.list
ecl> list = "ccds.list"
ecl> while (fscan (list, s1) != EOF) {
>>>    i = strstr ("opi", s1)
>>>    s2 = substr (s1, 1, i-1) // "opd" // substr (s1, i+3, 999)
>>>    hedit (s1, "BPM", s2, show-, verify-)
>>>  }
ecl> list = ""

The task for combining is imcombine. One might be aware of msccombine but that is different in that it combines MEF files by matching extensions to produce another MEF file which is not what is needed. Because the pixels cannot overlap in a single exposure there is no need to worry about combining algorithm. The pipeline also generally has matched the CCDs photometrically (i.e. flat fielded) and to a common background (as noted earler) so there is no need to worry about any scaling factors. To convert each MEF to a simple image:

ecl> files c4d*list > listoflists
ecl> list = "listoflists"
ecl> while (fscan (list, s1) != EOF) {
>>>    s2 = substr (s1, 1, strstr("_opi",s1)-1)
>>>    imcombine ("@"//s1, s2//"_im", bpm=s2//"_bpm.pl", offsets="world",
>>>        masktype="novalue", maskvalue="2", imc="")
ecl> }

Jan 10 14:26: IMCOMBINE
  combine = average, scale = none, zero = none, weight = none
  blank = 0.
  masktype = noval, maskval = 2
                Images                Offsets    Maskfile
  c4d_160704_052117_opi_i_v1.fits[1] 16516 24143 c4d_160704_052117_opd_i_v1.fits[1]
  c4d_160704_052117_opi_i_v1.fits[2] 8236 24108  c4d_160704_052117_opd_i_v1.fits[2]
  c4d_160704_052117_opi_i_v1.fits[3] 18596 21962 c4d_160704_052117_opd_i_v1.fits[3]
[...]
  c4d_160704_052117_opi_i_v1.fits[59] 16611   30 c4d_160704_052117_opd_i_v1.fits[59]
  c4d_160704_052117_opi_i_v1.fits[60] 8328    0  c4d_160704_052117_opd_i_v1.fits[60]

  Output image = c4d_160704_052117_im, ncombine = 60
  Bad pixel mask = c4d_160704_052117_bpm.pl

[...]
ecl> list = ""

This command produces a bad pixel mask as well using all the input bad pixel masks. The individual values from the input masks are not propagated. The output has value of 0 for good, 1 for no pixels, and 2 for a flagged pixel in the input masks.

Hint: These steps are also what is used to produce a simple image from the archived MEF stacks which have been sliced into contiguous tiles to limit the size of individual extensions.

Stack of Exposures

Stacking the exposures from the last step is similar and uses imcombine again. But since now there are overlapping pixels the exposures must be matched photometrically and by sky. The steps below use the magnitude zero points and sky levels in the header that are provided by the pipeline. Note that you are free to define these scalings as you find appropriate. The values can be taken from keywords (this example), files, or simple statistics from the data.

ecl> real    magref, skyref
ecl> files *im.fits > exp.list
ecl> list = "exp.list"
ecl> for (i = 1; fscan (list, s1) != EOF; i += 1) {
>>>    hselect (s1, "$MAGZERO,$SKYSUB", yes) | scan (x,y)
>>>    if (i == 1) {
>>>      magref = x
>>>      skyref = y
>>>    }
>>>    ;
>>>    z = 10. ** (0.4 * (magref - x))
>>>    hedit (s1, "STKSCALE", z, add+) 
>>>    z = skyref / z
>>>    hedit (s1, "STKZERO", z, add+) 
>>>  }
ecl> list = ""

Now the images can be stacked with flux and offset scalings as well as any of the other task parameters and the bad pixel masks.

ecl> imcombine ("@exp.list", "stack", bpmasks="stack_bpm.pl", combine="average",
>>>    offsets="world", masktype="novalue", maskvalue="2", imcmb="",
>>>    scale="!STKSCALE", zero="!STKZERO")

Jan 10 10:11: IMCOMBINE
  combine = average, scale = !STKSCALE, zero = !STKZERO, weight = none
  blank = 0.
  masktype = noval, maskval = 2
                Images       Scale  Zero    Offsets  Maskfile
  c4d_160704_052117_im.fits  1.000  2621.8    0    0 c4d_160704_052117_bpm.pl
  c4d_160704_052828_im.fits  0.872  2621.8  311   13 c4d_160704_052828_bpm.pl
  c4d_160704_053537_im.fits  0.639  2621.8  314  295 c4d_160704_053537_bpm.pl
  c4d_160704_054247_im.fits  0.775  2621.8   21  289 c4d_160704_054247_bpm.pl

  Output image = stack, ncombine = 240
  Bad pixel mask = stack_bpm.pl

As with the individual exposures, this command produces a bad pixel mask as well using the input bad pixel masks. The indivual values from the input masks are not propagated. The output has value of 0 for good, 1 for no pixels, and 2 for a flagged bad pixel in the input masks. In this case overlapping pixels allow for filling in gaps with combining only good pixels and so the no data flag will occur only where all input exposures have no data. Also where all the input pixels that do overlap have a bad value, such as the cores of saturated stars, the input pixel values are combined to provide a representative value but the pixel will be marked as bad.

As a warning, the time it takes to combine into a stack is dependent on the size and number of images. On a capabile machine a handful of exposures is fairly quick. One hundred DECam exposures will take 6 hours or more. A median is also slower than an average.

Projection Matching

If the pixel binnings do not match there is a command to reproject to a desired target projection. The command is mscimage. Defining the target projection can be done in a few ways. In the context of stacking the CP data products the simplest way is to specify one image as the target projection. This would usually be a resampled exposure representing the most common tangent point. There is an example below, which is for data before the current archive file naming convention. The BPM keywords have been set appropriately already.

ecl> hselect tu1738350[1],tu1739649[1] $I,CRVAL* yes
tu1738350[1]    1.608938600000E+02      -6.000000000000E+01
tu1739649[1]    1.629050290000E+02      -6.000000000000E+01
ecl> mscimage tu1738350.fits tu1738350_new.fits ref=tu1739649[1] \
>>>  wcssource=image format=mef interp=lsinc17 xxo=7 xyo=7 yyo=7 yxo=7
Resampling tu1738350.fits[S29] ...
Resampling tu1738350.fits[S30] ...
Resampling tu1738350.fits[S31] ...
Resampling tu1738350.fits[S25] ...
Resampling tu1738350.fits[S26] ...
[...]
ecl> hselect tu1738350_new[1],tu1739649[1] $I,CRVAL* yes
tu1738350_new[1]        162.905029      -6.000000000000E+01
tu1739649[1]    1.629050290000E+02      -6.000000000000E+01
ecl> dir *new*
tu1738350_new.fits        tu1738350_new_bpmN28.pl   tu1738350_new_bpmS2.pl
tu1738350_new_bpmN1.pl    tu1738350_new_bpmN29.pl   tu1738350_new_bpmS20.pl
tu1738350_new_bpmN10.pl   tu1738350_new_bpmN3.pl    tu1738350_new_bpmS21.pl
[...]

The masks, which are also reprojected, are output as separate files. The assignment to extensions in the reprojected image file is automatically updated in the header BPM keywords. It is possible to pack them into a single MEF as with the CP data products. How to do this may be added later to this document.

After this mscimage step the stacking can proceed as described earlier. Note that this mechanism can also be used with the unprojected data products as well.

Figure: The central part of one image and the stack.