Deep Full-sky Coadds from Three Years of WISE and NEOWISE Observations

There are two main types of full-sky coadded WISE data products presently available. First, there are the Atlas stacks created by the WISE/NEOWISE team. The algorithms and code (called "AWAIC") employed to construct the Atlas coadds are described in Masci & Fowler (2009). AWAIC uses a point response function kernel to interpolate from single exposures to the coadd, resulting in stacked images that are optimally matched filtered for point source detection (see e.g., Zackay & Ofek 2017a, 2017b). As a result of this matched filtering, the Atlas coadds have angular resolution, which is reduced relative to the native WISE beam. AWAIC does include a resolution enhancement option called "HiRes," but this functionality was only applied to select sky regions (e.g., Jarrett et al. 2012) and was never included in the Atlas data products. The most recent full-sky set of Atlas coadds is that of the AllWISE release and is based solely on pre-hibernation exposures.

Alternatively, the unWISE (Lang 2014) line of coadds uses Lanczos interpolation (for a review see, e.g., Getreuer 2011) to preserve the native WISE resolution. The present work is part of an ongoing effort to upgrade the W1/W2 unWISE coadds with each new release of NEOWISER L1b exposures, enhancing the depths achieved and using added redundancy to remove artifacts. The unWISE stacks are designed for forced photometry, a use case for which the coadds need not be optimized for point source detection, because the locations and morphologies of all sources are fixed based on prior knowledge from external data sets. Our coaddition is slightly suboptimal in that we do not attempt to account for small variations in the W1/W2 PSFs, either with time or across the focal plane. Over the course of the WISE instrument's lifetime, the PSF has broadened by ∼3% (∼2.5%) in W1 (W2).⁶ The rms variation of the PSF FWHM across the focal plane is 2.5% in both W1 and W2.

unWISE coaddition is meant to be lightweight/minimalist, with the plurality of compute time consumed by reading in the large number of L1b frames that overlap a given coadd footprint, and subsequently writing the output coadd images. As such, no truly robust outlier rejection methods (e.g., median filtering the resampled L1b intensities at each coadd pixel location) are employed, as these would dramatically increase the total computational cost.

Given that we have already produced full-depth coadds including all exposures through the first-year NEOWISER release (Meisner et al. 2017b), one might imagine a scheme in which we update our coadds by only processing the "new" exposures, adding them onto the existing coadd. Such an approach would not fully take advantage of improved outlier rejection enabled by the two newly added sky passes from second-year NEOWISER observations. For instance, we would inherit artifacts present in previous iterations of the unWISE coadds, rather than leveraging the newly added redundancy of the latest set of exposures to improve artifact removal (see e.g., Figure 5). Such an approach would also likely require us to implement a series of contrived shortcuts/hacks.

Therefore, with every new release of NEOWISER exposures, we opt to rebuild all of our full-depth coadds from scratch. Even when jointly processing all exposures spanning 2010 January to 2015 December, as we have chosen to do, the total full-sky computational expense is only of the order of tens of thousands of CPU hours.

In combining observations that span the full WISE lifetime, we also strive to achieve the best possible relative calibration of all exposures. To accomplish this, we adopt a custom photometric calibration using repeat measurements of calibrator sources near the ecliptic poles, as described in Section 4 of Meisner et al. (2017b). We have now performed this photometric calibration analysis for the second year of NEOWISER observations. The results are shown in Figure 1. Changes in the zero-point with time are predominantly due to variation of the WISE beam-splitter assembly temperature (Cutri et al. 2015), and the ≲1% zero-point changes we measure are small enough to be explained by this effect.⁷

Zoom In Zoom Out Reset image size

One final component of our philosophy in generating WISE coadds is to throw away as little data as possible. This motivated our work in Meisner et al. (2017b) to recover exposures corrupted by scattered moonlight. In this same vein, we retain all frames at/near the ecliptic poles, despite the added computational cost associated with the exceptionally large depth of coverage in these regions. As a result, the tiles at the North/South ecliptic poles have peak integer coverage of ∼17,000 frames per band, for a total integration time of ∼1.5 days.

unWISE coadds employ the same set of 18240 tile centers as do the AllWISE Atlas stacks. We adopt a pixel scale matching that of the W1/W2 L1b exposures themselves, 275 per pixel. Each coadd astrometric footprint has its −x (+y) axis aligned with celestial east (north). We refer to a coadd astrometric footprint as a "tile," and identify tiles by their coadd_id values, which are strings encoding their central (R.A., decl.) coordinates, e.g., "0000p000." Our coadds are 156 (2048 pixels) on a side.

Here, we briefly summarize the core steps of the unWISE coaddition methodology, with an emphasis on the outlier rejection procedures. Refer to Lang (2014) for additional details.

Each pairing of a coadd_id and a WISE band (in this work, either W1 or W2) corresponds to a set of coadd outputs that must be generated, and each such pair is processed independently (i.e., a W1 coadd does not take into account any information about its W2 counterpart or its neighboring coadd images in either band). The following list of steps is undertaken for each (coadd_id, band) pair.

• 1.  
  We start by identifying all exposures in the relevant band that are near enough to the coadd footprint of interest to contribute to the coadded image, using a look-up table containing central coordinates for all L1b images. We ignore exposures with a frame-level quality score (qual_frame) of zero.
• 2.  
  Frames that are flagged by the WISE team as potentially corrupted by scattered moonlight are rejected if their pixel value histograms indicate an unusually high robust standard deviation relative to other frames in the same sky location that are not affected by the Moon. These assessments can be made prior to reading in the L1b images, using frame-level metadata provided by the WISE team.
• 3.  
  For the remaining valid exposures, we read them in and resample them onto the coadd astrometry. In doing so, we simply make use of the L1b WCS solutions provided by the WISE team. Performing a custom recalibration of the L1b WCS or tweaking the alignment of the many L1b images relative to each other would likely improve the astrometric consistency of the contributing exposures, although we have not yet implemented any such procedures. The resampling is performed using using third order Lanczos interpolation. We employ a kernel that is normalized so as to preserve the multiplicative scaling of the input images when resampling.
• 4.  
  At this stage, we have a set of typically ∼1000 resampled L1b images held in memory. We begin with a simplistic "first round" of coaddition, which consists merely of computing the mean of the resampled L1b values and its variance at each pixel in coadd space. The resulting first-round coadds will serve as the basis for rejection of outlier pixel values downstream. Section 4.3.1 describes a newly introduced aspect of the unWISE first-round coaddition, whereby a "min/max" rejection step is applied to avoid contamination of the first-round coadd by severe defects that may be present in a small number of L1b images.
• 5.  
  Using coadd-space images of the mean pixel value and its standard deviation provided by the first-round coadd products, we re-examine all of the resampled L1b images. If a resampled L1b pixel is >5σ discrepant relative to the mean in the first-round coadd, it is flagged as an outlier. We reject exposures for which >1% of pixels are flagged as outliers.
• 6.  
  With the resampled images and their corresponding outlier masks from the previous step in hand, we perform a final "second round" of coaddition, accumulating resampled L1b pixel values at each location in coadd space. During both first and second round coadd accumulation, inverse variance weights are applied, the details of which can be found in Lang (2014). We create two versions of each coadded image, one in which we interpolate over the outlier pixels identified in each resampled L1b exposure before coadding (output files with names containing "-u-") and a second in which we simply ignore such pixels by giving them zero weight (output files with names containing "-m-").

The foremost improvement achieved in this work relative to the coadds of Meisner et al. (2017b) is the ∼50% enhancement in depth of coverage obtained by including second-year NEOWISER observations. However, because we uniformly reprocess all exposures following every NEOWISER release, each new, deeper set of coadds also represents an opportunity to build further processing improvements into the unWISE coaddition pipeline. The following subsections describe newly introduced facets of the unWISE processing and data products.

4.3.1. Min/Max Rejection

As mentioned previously, unWISE coaddition does not use robust statistics to perform outlier rejection, instead relying on naive first-round coadd mean and standard deviation images to flag problematic pixels in the L1b exposures. This computationally cheap approach is vulnerable to the appearance of extreme outliers in a small number of exposures. The most conspicuous example of this phenomenon arises due to large numbers of cosmic ray strikes in frames affected by the South Atlantic Anomaly, which mainly impacts observations at −45° < δ < −10°. This leads to first-round coadd results that are strongly contaminated by cosmic rays, and therefore are relatively ineffective at flagging outliers (see left panel of Figure 2).

Zoom In Zoom Out Reset image size

To combat such corruption of our first-round coadds, we have now implemented a "min/max" rejection step during first-round coaddition. Specifically, for each pixel in coadd space, we precompute which L1b exposures contribute the highest and lowest single-frame intensity values at that location. Then, when computing the first-round coadd mean and standard deviation, we disregard the minimum and maximum single-exposure values at each coadd pixel location. Precomputing the minimum and maximum values at each coadd pixel location requires ∼0.01 s per exposure, which is negligible relative to the overall runtime of coaddition, and typically leads to ∼10 additional seconds of CPU time per tile at low ecliptic latitude.

The right panel of Figure 2 shows an example of dramatically reduced cosmic ray contamination in a first-round coadd mean image achieved by incorporating our min/max rejection step. Obtaining a cleaner first-round estimate of the per-pixel mean and standard deviation improves our ability to reject outliers during the second round of coaddition. This is illustrated by the example region shown in Figure 3, where min/max rejection has eliminated the imprint of cosmic ray strikes from our final coadd image. We emphasize that our min/max rejection step is applicable to other types of outliers beyond cosmic rays and that there already exists an extensive literature on outlier rejection in the context of image processing (e.g., Hasler et al. 2004; Ieng & Charbonnier 2007; Zibetti & Mayer 2007).

Zoom In Zoom Out Reset image size

Note that our min/max rejection step does not necessarily result in a reduction of integer coverage for our final outputs at every coadd pixel, as it is only applied during the first-round coaddition process. In future iterations of our full-depth coadds, as the typical integer coverage continues to increase, we may investigate more aggressive min/max rejection in which the highest/lowest N single-exposure values at each coadd pixel location are discarded, with N > 1.

4.3.2. Bitmasks

One feature that has been absent from previous unWISE data releases is a pixel-level bitmask corresponding to each coadd image. Such a bitmask data product could be very helpful for selection of spectroscopic samples (e.g., for DESI) in order to improve efficiency by flagging sources associated with artifacts. It could also be quite useful in aiding rare object searches, by indicating when unusual objects might in fact be of instrumental rather than astrophysical origin. For instance, a close (W1-only) variant of the bitmask procedure presented here was used to great effect in the Planet Nine search of Meisner et al. (2017a).

In this work, we address these needs by creating a pixel-level bright star mask output for each coadd tile, with suffix -msk.fits.gz. We construct one such bitmask image per tile, containing information about bright stars in both W1 and W2, rather than generating two separate bitmasks with one corresponding to each band.

The simplest bright star mask one might imagine would flag a circular region about every sufficiently bright star, with the radius scaled based on each star's total flux according to some empirically calibrated prescription. Indeed, the eBOSS collaboration currently employs masks of bright WISE stars created in this way (Myers et al. 2015; Prakash et al. 2016). Our bitmasks are designed to go beyond this simplistic approach, taking into account both the WISE PSF and the WISE scan direction to create highly detailed masks that include features such as diffraction spikes and optical ghosts.

As the basis for our bitmasks, we must select samples of very bright WISE sources in each of W1 and W2. We do so by making use of the positions and fluxes in the AllWISE catalog. We create a full-sky list of W1-bright sources by selecting objects from the AllWISE catalog with w1mpro < 9.5. For W2, we adopt a brightness threshold of w2mpro < 8.3. Our masks are primarily designed with extragalactic/cosmology applications in mind, and therefore are optimized for high Galactic latitude sky regions. In locations far enough off the Galactic plane to be part of DESI's footprint ($| {b}_{\mathrm{gal}}| \gt 18^\circ$), our bright star sample contains on average ∼29 (∼11) sources per square degree in W1 (W2).

For each coadd_id astrometric footprint, we begin constructing its bitmask by identifying all bright stars that may have profiles at least partially falling inside of the tile boundaries. We properly account for all bright stars that overlap the tile footprint of interest with any portion of their PSF wings/diffraction spikes, even those for which the centroid falls outside of the tile boundaries. For each bright star, we then render a model of its appearance given its AllWISE centroid and flux. To render these models, we make use of the W1 and W2 PSFs of Meisner & Finkbeiner (2014). These PSFs, displayed in Figure 4 of Lang (2014), are 149 on a side and therefore extend quite far into the wings, capturing details such as diffraction spikes and the W2 optical ghost.

One subtlety involved in rendering our bright star profiles is that the Meisner & Finkbeiner (2014) PSFs provide models of the single-exposure PSF, not the effective PSF obtained after resampling onto the unWISE tile astrometric footprints and coadding. For instance, in L1b images, the diffraction spikes always emanate from bright sources at 45°, 135°, 225°, and 315° degrees from the $+{x}_{L1b}$ direction. However, this is not the case for unWISE tiles, which are oriented along the equatorial cardinal directions, whereas the L1b exposures are very nearly oriented exactly along the ecliptic cardinal directions.

A further subtlety is that the Meisner & Finkbeiner (2014) PSF models require specification of detector (${x}_{L1b}$, ${y}_{L1b}$) coordinates, as they incorporate PSF variation across the single-exposure field of view (FOV). Therefore, before rendering any bright star models during bitmask construction, we precompute mean versions of the W1 and W2 Meisner & Finkbeiner (2014) PSFs, averaging over (${x}_{L1b}$, ${y}_{L1b}$) detector location.

A final subtlety is that, except for a small area near $| \beta | =90^\circ$, there will be two discrete WISE scan directions contributing to each coadd, corresponding to the ecliptic north and south directions. The PSF is not perfectly symmetric with respect to swapping the scan direction. The most notable asymmetry in the bands of interest is the W2 ghost, which is offset by ∼5' from its parent bright source centroid, at a position angle that is fixed relative to the scan direction. In practice, this means that, in our coadds, the ghost will appear on opposite sides of its parent bright star in exposures with opposite scan directions. This explains why the doughnut-shaped ghost in Figure 4 appears twice.

Zoom In Zoom Out Reset image size

In order to create a coadd-level model rendering of a bright star, we compute the angle between the scan direction and celestial north, rotate the FOV-averaged PSF by this angle, and scale the rotated profile so that its total flux matches that quoted by the AllWISE catalog. For each bright star, we create two such renderings, one for each scan direction. Pixels in the bright star models above a threshold of 13.2 (38.3) Vega nanomaggies per square arcsecond in W1 (W2) are flagged according to the mask bits listed in Table 1. At very high latitudes, $| {b}_{\mathrm{gal}}| \gt 80^\circ$, ∼0.4% of the total sky area is masked, on average. This fraction ramps up toward lower $| {b}_{\mathrm{gal}}|$, reaching 1.6% at $| {b}_{\mathrm{gal}}| =18^\circ$, the closest DESI will observe to the Galactic plane.

Table 1.  Description of Bright Star Mask Bits

Bit	Description	Scan Direction
0	W1 bright star	south
1	W1 bright star	north
2	W2 bright star	south
3	W2 bright star	north

Download table as:  ASCIITypeset image

Our pixel-level bitmasks have now been propagated into the DESI imaging Legacy Survey DR4 catalogs. An example bitmask image near an extremely bright star is shown in Figure 4. Much room remains to extend our present bitmasks in future unWISE coadd releases, making them even more intricately detailed and elaborate. The primary avenues for doing so are as follows.

• 1.  
  Near the Galactic plane, our current bitmasks flag an excessively large fraction of pixels, rendering the masks unhelpful in these regions. In the future, we could counteract this issue by making our magnitude thresholds for bright star masking dependent on, e.g., Galactic latitude, so that larger flux would be required to trigger masking of a star at low $| {b}_{\mathrm{gal}}|$ than at high $| {b}_{\mathrm{gal}}|$.
• 2.  
  As seen in Figure 4, the diffraction spikes of extremely bright stars can sometimes extend beyond the 149 boundary of the Meisner & Finkbeiner (2014) PSF models. In the future, we could calibrate look-up tables of diffraction spike length versus magnitude in W1 and W2, and use these to create geometric masks composed of lines emanating from the centroids of extremely bright stars at angles of 45°, 135°, 225°, and 315° relative to ecliptic north. Such a procedure was employed in W3 by Meisner & Finkbeiner (2014), as described in their Section 5.2.4.
• 3.  
  "Latents" represent another class of defects associated with stars bright enough to saturate in their cores. Latents are persistence artifacts that appear as diffuse positive blobs at the detector positions of saturated pixels, but in the frame immediately after imaging of the parent bright star. The locations of all latents can be predicted exactly, given L1b metadata and a catalog of bright stars, so it would be possible to reliably flag these artifacts by adding new mask bits.
• 4.  
  Our assumption that there are exactly two discrete scan directions corresponding to ecliptic north and south is violated at the ecliptic poles, where continuous coverage yields a continuum of scan directions. As a result, very near $| \beta | =90^\circ$, bright star diffraction spikes become spread out into disks in the coadds and are correspondingly attenuated. One approach for taking these effects into account has been described in Section 5.2.4 of Meisner & Finkbeiner (2014).

Relative to the coadds of Meisner et al. (2017b), we have in general increased the depth of coverage by ∼50% while further mitigating artifacts. The mean integer coverage over the entire sky is 108 (107) frames in W1 (W2), and every pixel has integer coverage of at least 33 (30) in W1 (W2). Figure 5 illustrates that folding in a third full year of observations has reduced the impact of scattered moonlight, the dominant systematic problem with the W1/W2 unWISE coadds on large angular scales. Because the Moon affects different ranges of ecliptic longitude during different sky passes, any remaining traces of scattered moonlight contamination become further suppressed as additional NEOWISER coverage is incorporated. Figure 6 shows a typical example of the reduction in statistical noise achieved for coadds at low ecliptic latitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

DR4 of the DESI imaging Legacy Survey, hereafter referred to as simply "DR4," contains the only currently available photometry based on the full-depth W1/W2 coadds presented in this work. This photometry is forced rather than WISE-selected, adopting source locations and morphologies derived from deep optical data. DR4 covers an ∼4500 square degree extragalactic footprint limited to δ ≳ +30°, as shown in Figure 3.19 of DESI Collaboration et al. (2016a).

In validating our coadds, we first seek catalog-level confirmation of the decrease in statistical noise suggested by the images in Figure 6. We examine objects with W1/W2 forced photometry available from both Lang et al. (2016), based on the original unWISE coadds and DR4. We select a comparison sample drawn from the region 170° < α < 230°, 425 < δ < 48°, using a 1'' match radius and restricting to objects with DR4 forced photometry signal-to-noise in the range 10 ± 1. For these 127,000 (144,000) sources in W1 (W2), we find that forced photometry flux uncertainties from DR4 are smaller than those of Lang et al. (2016) by median factors of 1.68× (1.53×). These values indicate that forced photometry based on the present coadds is 0.56 (0.46) mags deeper in W1 (W2) than forced photometry based solely on pre-hibernation W1/W2 imaging. For both the W1 and W2 samples, the median increase in integer coverage is a factor of 2.94×, leading us to expect flux uncertainties reduced by 1.71× under the assumption that all WISE observations have maintained the same sensitivity regardless of mission phase.

We additionally check that forced photometry fluxes derived from our new coadds are consistent with the AllWISE catalog. We select a comparison sample drawn from the region 170° < α < 230°, 425 < δ < 48°, using a 1'' match radius. We further require a one-to-one DR4-AllWISE match, DR4 morphological type PSF and w?cc_map = 0 in the band of interest. This yields 626,000 (612,000) sources in W1 (W2). Figure 7 summarizes our DR4-AllWISE comparison for this sample. There is good overall agreement across a wide range of fluxes, and the (DR4 − AllWISE) offset asymptotes to 4.1 (−1.6) mmag in W1 (W2) toward the bright end. The faint end upturn has previously been noted in Lang et al. (2016) and Meisner et al. (2017b) and is due to Malmquist bias, as the AllWISE sources are WISE-selected while the DR4 objects are not.

Zoom In Zoom Out Reset image size

The full-depth W1/W2 coadds described in this work are publicly available via the unWISE web interface at http://unwise.me.