Sky Subtraction

If you are doing stellar photometry, you will be doing local sky subtraction using the APPHOT and DAOPHOT procedures, and you will not need to do the global subtraction that is described here. This page is for people working on extended sources.

1. Simple Sky Subtraction

The simplest form of sky subtraction is very straightforward. Find a rectangular region of blank sky far away from everything, containing no detectable stars, residual cosmic rays, etc. Determine the coordinates of the region corners, then use 'imstat' to measure the average brightness level in counts/pixel. For example:

cl> imstat zn2013[745:930,810:950]
#               IMAGE      NPIX      MEAN    STDDEV       MIN       MAX
 zn2013[745:930,810:950]     26226     81.62     10.38     42.46     173.2

If you want to be somewhat safer, pick a few sky regions around the image and average the results. And if you want to be safer still, epar imstat and ask it to return the median rather than the mean. The median is less influenced by outlier events such as cosmic rays, hot pixels, etc. Once you have determined the sky level you need only to subtract it from the image:

cl> imexpr
expression : a-b
output image : zn2013sub
operand a : zn2013
operand b : 81.62
10% 20% 30% 40% 50% 60% 70% 80% 90% 100% - done

It pays to think a little about how precisely you need to determine the sky value. In the example, the standard deviation (the root-mean-square deviation from the mean) for all pixels in the region is 10.38. The uncertainty in the mean value is the standard deviation divided by the square root of the number of pixels. Thus the sky level is 81.62 +/- 0.06. This is quite a precise measurement--- IF the errors are dominated by pure counting statistics and not by other systematic effects.

2. Getting More Sophisticated

Whether this simple procedure will work for you depends on your data and what you are trying to do with it. If you are working at low surface brightness levels it you will probably have to do better, because it is not safe to assume that the "sky" (i.e. the true sky plus whatever residual junk still lurks in your images after basic reductions) is flat.

There are at least 3 reasons why your sky may not be flat:

Your flat field correction may not quite have removed all of the variations in sensitivity and/or illumination across the detector;
Your data may be contaminated with uncorrected dark current, or light leak;
You may have a bright star in or near your field giving you scattered light or just bright PSF wings.

If you have 2 or more blank sky exposures (no bright or extended sources in the field), you may be able to deal with either of the first two possibilities. The hard part is that if all you see is a spatial variation in the sky level, it's very hard to know what the cause is.

The first thing to do is to combine your blank sky frames. (This does not necessarily mean all of your blank sky frames, e.g., if you have interleaved object and sky exposures. Which you combine depends on what you did at the telescope.) Use the generic 'combine' task to compute the average of the exposures if they all have the same exposure time, or the sum if they do not. Turn off any image scaling by setting "scale" to "none". Recommended rejection algorithms are "crreject" or "minmax"; if you use the latter, start with "nhigh" and "nlow" both set to 1 if you are combining more than 4 frames. Check that the combined image has no sources in it. Next, use 'mkskycor' to smooth the image. (See section 3.11.2 and Fig. 18 of the User's Guide to CCD Reductions with IRAF.) The defaults will probably work, though you may want to increase the clipping sigma parameters.

Now, what you are going to do with this smoothed image depends on whether you think you have an additive contamination (light leak, dark current) or a multiplicative correction (flat-fielding error). Our current state of experience with the ST8 camera suggests that the former is probably the cause of big (say, 20%) variations in the apparent sky level.

To correct an additive contamination, first figure out what the effective exposure time is for the smoothed combined sky image. If you averaged frames of the same exposure time, then the resulting image retains that exposure time; if you summed frames, then you want the sum of the exposure times. Scale this to the exposure time of the image you are correcting, and subtract. For example, suppose you computed the smoothed sum of 4 blank sky exposures of 300, 300, 600, and 600 seconds. The sum has an effective exposure of 1800 s. To use this to correct a 300 s exposure of a galaxy, you would do the following:

cl> imexpr
expression : "a-b/6."
output image : zn2013corrected
operand a : zn2013
operand b : SmoothedSky
10% 20% 30% 40% 50% 60% 70% 80% 90% 100% - done

(Note the quotation marks around the expression, necessary if the expression contains a "/".) You may find that the best subtraction is obtained with a scale factor other than that implied by the ratio of the exposure times. This could happen if the sky brightness was changing in time, or if your contamination is a nonlinear dark current. A little trial and error may be necessary to get the best results. The bad news is that it appears that this procedure works poorly for GOT/ST8 data unless the sky and target frames have the same exposure times to begin with.

If, instead, you want to use the blank sky frame to correct for residual flat-fielding errors, then you need only take another pass through 'ccdproc'. Set the input name to the image you want to correct. Turn on the illumination correction, and set the name of the illumination correction image to your smoothed sky image. Run the task. Note that you need to do this before you do any other sky subtraction.

Last updated 2008 May 29. Written and maintained by Tom Statler