Precision sky subtraction

I am preparing for a trip to Chile to observe extremely faint asteroids using DECam, and to visit a collaborator.  In parallel with this, I am working on processing data from my March 6-13 observing run on the Kitt Peak 2.1m telescope.  My main objective in processing this data now is to get astrometry on the faint asteroids we observed, because this is highly time-sensitive.  At the same time, it is most efficient to do some processing on the additional data while I have the files open and everything is fresh in my mind.  The additional data, which was actually the main objective of the Kitt Peak observations, consists of images for photometric monitoring of brown dwarfs in the red optical (the f814w filter, which is similar to the I-band but wider).

Optimal photometry of faint objects in the red optical requires accurate sky subtraction.  This is because the strong interference fringes produced in silicon detectors by night-sky emission lines at these wavelengths produce a bright and spatially varying sky background.  I have some tricks for accurate background subtraction that will become apparent as I describe the sequence of operations.

Observing: Every image should be taken in a different position from any other.  The telescope pointing should not be altered very much, just a small fraction of the field of view, but the pointings should be dithered in a non-redundant manner, as far as is possible.  The data set I worked with today had 63 images, each a 200 second exposure in the f814w filter.  The dithering was good: it was almost entirely non-redundant.

Processing: I begin by trimming off the overscan region (which is useless on the CCD I was using) and also trimming off edge regions extending 10-20 pixels into the valid data region, where the sky background and flatframes for this CCD behave strangely.  Then I follow the usual steps of zero subtraction and flatfielding.  After that it gets interesting.

The first step is to do an approximate sky subtraction, which is required to prepare the final, much more accurate sky subtraction.  To do this, I normalize all of the processed science images to have a sky background of 1.0, and then do a median-stack of them.  Due to the non-redundant dithering, the stars mostly vanish in this stack, leaving me with a fairly clean image representing the average sky background through the night.  Call this the master sky image.

I then perform sky subtraction on each science image by creating a version of the master sky image that is scaled to have the same sky background, and then subtracting that from the science image.  This sky subtraction is fairly effective, producing images that are substantially cleaned of fringing, etc.

I then pad, shift, rotate, and median-stack all of the science images to create a master image that has much lower noise and shows much fainter objects than each individual science image.  Here is an example of my master science image at this point:

Very faint stars and galaxies can be seen, but there are dark halos around bright stars.  This is because weak residuals of the bright stars were still present in the master sky image: despite the non-redundant dithering, these stars did not completely vanish from the median stack used to create the sky image. The master science image above also has dark and noisy borders, and the edges of some individual frames can still be seen.  This shows that the master sky image, besides retaining weak residuals of bright stars, was also imperfect in other ways.  Such imperfection is not surprising, because my processing thus far has assumed that the sky background remained exactly the same, except for its average brightness, throughout the entire night.  In fact the sky background is composed of different components: fringing from night-sky emission lines; moonlight in the sky; scattered moonlight in the dome, etc.  The relative contributions of these different components change over time, and therefore using a single sky frame for the whole night is a rather crude approximation.

I can solve this problem by dividing up the science images into several sets, and producing different sky frames for different parts of the night.  However, if I do this, the problem of the bright stars not entirely vanishing from the median stacks used to make the sky frames will get much worse.  For example, a median stack of ten images will have much stronger residuals from bright stars, even if every image was taken at a distinct pointing, than a median stack of 63 images.

What is needed is to mask the stars, especially the bright ones.  But it is not sufficient simply to mask the bright cores of the stars: those vanish in the median combine anyway.  It is the faint halos around the bright stars that cause the most trouble.  Thus, it is necessary to mask the stellar halos -- extending out even to where they are so faint that they may not be detectable on individual frames.

I deal with this be creating a star mask based on the master image, in a stepwise fashion.  For this data set, I masked all pixels within:

-A 25-pixel radius of any pixel with brightness above 50000 ADU, or
-A 20-pixel radius of any pixel with a brightness above 10000 ADU, or
-A 12-pixel radius of any pixel with a brightness above 1000 ADU, or
-A 6-pixel radius of any pixel with a brightness above 100 ADU, or
-an individual brightness above 20 ADU

Where all of these cuts are performed on the low-noise master image.  I should say that the noisy edges of the master image, where only a few science images contribute data, are trimmed off before the mask is constructed.

It remains to apply this mask to the individual science images, each of which is taken at a different pointing.  To handle this, I made my shifting and rotation routine save a log file describing exactly what shifts and rotations were applied to each image.  Thus, I can mask each science image by simply creating a version of the master mask-image with the reverse shift and rotation appropriate to that particular science image applied. This produces a set of images that are unaltered from the original, individual science frames, except that all of the stars are masked -- including even the objects too faint to be detected on individual science frames, and the faint, extended halos of bright stars.

I normalize these masked images to have a mean sky background of 1.0, and stack them in subsets to make multiple master sky images that follow the changes in the sky background throughout the night.  For the example of this data set, I created five master sky images from masked, normalized versions of science images 1-13, 14-26, 27-39, 40-52, and 53-63, respectively.  Because of the masking, these master sky images had no residual star images at all.  The non-redundant dithering prevented blank holes from appearing in these master sky images, by ensuring that there was no pixel region that was masked on all the frames used to produce a given master sky image. Thus, the master sky images had valid data everywhere, an outcome that would have been unlikely if the observers had not taken pains to make the dithers non-redundant.   The fringing in these sky images looks cool, so I show them here. Note the subtle changes in the contrast of the fringes and in the overall distribution of the background; in particular the contrast of the fringes is higher in the later images because the Moon had set:

I subtracted scaled versions of these master sky frames from the corresponding science images.  The resulting final sky subtraction was considerably better than the preliminary version, as evidenced by the new master science image made from the data with improved sky subtraction.  I show both the old and new master images together here, for comparison:

Old:

New:

While some edge effects remain, they are greatly reduced, and the dark halos around bright stars have completely vanished.  Neither fringing nor any other sky-related effect remains to introduce a spurious signal into our time-series photometry of the faint T dwarf that was the science target of these observations.