CCDLAB: A Graphical User Interface FITS Image Data Reducer, Viewer, and Canadian UVIT Data Pipeline

Given hardware and data rate limitations, it is not generally possible to download newly generated onboard UVIT data every time an observation is made with UVIT. Thus, observations actually performed in one orbit may be downloaded in a later orbit, and if the orbit number of data-download is used for identification purposes (and it is), then this can lead to confusion as to when downloaded data actually originates. That is, there may be two or more L1 files with unique file names and unique file sizes, and with metadata indicating that the data contained therein is from different orbits; however, upon inspection it can be found that the two (or more) files contain at least partially, or totally, identical data sets.

While all of the data is eventually downloaded and the correlation of downloaded data can be made to the observation target to which it corresponds, the uniqueness of L1 data files for a particular target is not guaranteed. That is, observation data for a target may be replicated at least once in the L1 files which are subsequently generated for the target and which are intended for ingest in any down-stream pipeline. For example, data for a 20-minute stare using some filter and some channel may be replicated among one or more L1 files with different orbit numbers.

Given that this is expected to be a general condition of the L1 data files, then any pipeline must ensure that the set of L1 files for a particular target do not become reduced as if all data is unique data. That is, combining two (or more) data sets, assuming that they are unique observations, when the data sets are actually of duplicated data, is not scientifically meaningful. The Canadian UVIT Pipeline in CCDLAB ensures that this does not occur in its reduction of UVIT data by examining the centroid timestamps from different data files for repetition, and discarding such.

More recently, L1 data files are also being provided as a "merged" format in which data repeats have already been removed, but the procedure described here can still be run to ensure data fidelity.

Related to the above-mentioned problem with data download limitations, and due to the fact that the filter electronics system and target identification metadata are packed into the L1 fits headers separately from the UVIT DS data, then the filter and target identification occasionally may not match to what was actually observed. To correct this, the L1 data set also contains auxiliary "housekeeping" files, as FITS tables, which correlate the time of observation to the filter and target ID. The time of observation used is an ASTROSAT-specific modified Julian date, and this table entry is found in both the L1 detector data table and the auxiliary housekeeping table. A pipeline must ensure that the filter and target ID actually correspond to what was observed, and by using the auxiliary housekeeping fits table this may be done.

More recently, L1 data files are being provided as a "1.2" version, in which the filter metadata has already been corrected, but the procedure described here can still be run to ensure data fidelity.

Each detector has unique field distortion which should be corrected for astrometric accuracy (Girish et al. 2017). Field distortion is probably mainly due to the fiber optic taper, and runs in the range of 10–15 arcseconds.

To get the astrometric accuracy and source resolution for which UVIT is designed, subpixel-accuracy centroiding is used, but the centroiding is of a very poorly sampled pulse distribution. A photon pulse on the CMOS chip falls almost entirely within only 3 × 3 pixels, with a FWHM of ∼1.5 pixels. A weighted-mean centroid performed on such a poorly sampled profile leads to a preferential centroid computation bias tending toward the center pixel of the 3 × 3 centroiding kernel. This centroid bias must be corrected for so that the subpixel resolution final science images are not affected by the astrometric and qualitative image problems that this bias creates at the one-half pixel scale. Correction of the centroiding bias is described in Hutchings et al. (2007) and Postma et al. (2011), and the bias correction table is provided with the UVIT Calibration Database at http://www.ucalgary.ca/uvit/.

As a single photon pulse is typically a few thousand counts on the CMOS chip and is recorded as a single count for its centroid location, the flat response of the CMOS chip itself is of much less consequence as compared with traditional direct imaging with a CMOS or CCD. However, the photocathode, phosphor, microchannel plate, and fiber optic taper do still provide for some modulation of the quantum efficiency (QE) across the detector field. The flat field for each detector was measured in Calgary using a UV-reflective integrating sphere with the detector window placed at the output of the sphere where its response is flat. The FUV and NUV channels must apply the inverse of the flat value as a weighting for each centroid, given its location on the detector, when the 2D sum of all centroids is performed to produce the final science image. The same would need to occur for VIS centroids; however, VIS centroids typically will not be produced for science purposes (although it can be done if needed).

The flat field as measured by a detector must also first have the distortion correction applied to its centroids for that detector. If there is spatial pinching and stretching of the detection field across the detector due to field distortion, then pinching transforms the initial uniform distribution of the light source into a higher frequency of counts in the region where there is pinching, and stretching causes the opposite lower frequency of counts to occur in regions of the detector where there is stretching. This will make it look like regional variations in QE, but it is not a QE effect because it is actually due to field distortion bunching together or diluting the incoming flat-distribution of photon pulse centroids in a region. To get the actual flat field for QE effects, the flat centroids must be corrected for detector distortion first before the flat-field centroids are summed into a 2D image. The centroids of science data are, of course, corrected for distortion too, so their corresponding flat-field weight in the distortion-corrected location where they are detected must be from the QE effect only and not from an effect due to distortion.

The CCDLAB pipeline creates a flat-field list where each entry is the flat weighting for the centroids in the centroid list; this weighting is (optionally) applied to each centroid when they are summed into a 2D image.

The flat fields for each channel are provided with the UVIT Calibration Database at http://www.ucalgary.ca/uvit/.

To spread the charge depletion of the MCP over its surface area, ASTROSAT is commanded to oscillate its pointing with a sine drift about the center of the observation field with an amplitude of about 50'' and a rate of about 1'' per second. While this induced oscillation is very well behaved and low frequency, a generally unpredictable high-frequency field translation usually enters most observations, which is caused by movement of the scanning sky monitor (SSM) camera (Ramadevi et al. 2016). The induced field drift and the field movement caused by the SSM are, conveniently, translational only, with field rotation for a given observational pointing being negligible or undetectable. The field drift must be subtracted out of the centroid list.

Field rotation can come in for different orbits observing the same field, and if the observer desires a deep-field stack of all centroids for all observation orbits of a target, then each observation must be registered with both field translation and field rotation between them. CCDLAB provides an interactive routine to easily register all orbit observations for field translation and rotation.

The field drift implies that not all portions of the imaging field are exposed to an equal total integration time. Typically, this effect only applies around the periphery of the imaging field while the central target of observation would be fully exposed. To correct for the areas that are not exposed to the full integration time, then one can utilize the drift series from the previous section to create an "exposure array" which is a normalized 2D map of the ratio of the local total integration time within the map to the total integration time for the fully exposed part of the image. It helps to pad the imaging field to a larger array size than is actually imaged for this process, as some field-edge sources may drift out of the mean field of view. The exposure array then provides weighting values for the final image which enter in the same way that a flat-field correction would; that is, areas of the image which are exposed to a normalized value of less than unity are divided by this value to scale those areas up to the total integration time.

While the VIS channel of the UVIT telescopes is intended to be used for drift correction, due to problems with temporal alignment of data from different channels during first-light and PV (payload verification) phase an alternative algorithm was developed which allows the NUV and FUV channels to track their own drift.

During any single observation it is found that field rotation is negligible, thus drift determination and correction only requires a solution which determines field translation. Thus, reducing the x and y axes of the background-subtracted frame to 1D vectors by binning the frame along either axis produces a sort of "emission spectrum" where point sources become emission spikes in the vector. These data vectors can be cross-correlated between a reference frame and comparison frame to get the x and y lag shift in the position of the "emission sources" between the frames. A cosine-bell function (i.e., Hanning window) is applied to the images before they are summed into vectors for cross-correlation to suppress ringing and to suppress bright edge sources that may come in and out of view during drift. The correlation lag, or drift shift, is then determined by fitting a quadratic to the top three points of the cross-correlation and using the center of the quadratic as the lag. The author used this technique in undergraduate work at the University of Western Ontario under Dr. David Gray, for determining radial velocity variations in stellar absorption spectra at subpixel accuracy; here, the author simply "inverted the spectrum." This method of image registration can be used in general as long as there is no or very little field rotation; it is usually helpful to subtract out any background gradients in the source (reference and comparison) images, and the cross-correlation vectors should be zeroed-out by subtracting their mean value.

Each individual frame read from the UVIT DS can only contain a single count at each source location, plus randomly placed single counts across the frame due to background "noise." Therefore, if a number of seconds worth of frames is first summed (stacked) together, then the sources begin to come out while the background noise remains spread out and generally singular in count. The slew rate is slow enough that for a stack of several seconds, a source is more or less stationary at the CMOS pixel scale. The CCDLAB drift-determination menu gives an option for stacking 1 through 20 seconds worth of centroid data for usage by the algorithm, and 2 seconds is typically found sufficient. Large stack intervals degrade the resolution of the drift series particularly when the high-frequency translations occur due to movement of the SSM camera. The algorithm automatically sets all count values equal to 1 in the stacked frame to be equal to zero because these will be background noise; source counts will be higher than one in the stacked frame.

For example, 2 sequential seconds of comparison centroid data will be stacked together into an image (a reference image will have been created for the initial 2 seconds of the centroid list), counts equal to one in the image are set to zero (as with the reference image), the image binned to x and y vectors after having the Hanning Window applied (as with the reference image), the vectors cross-correlated to the reference ones, and their lag shift recorded. The lag shifts are then interpolated for each frame time within the 2-second stack with the center time of the stacked frame set being the time of the lag shift.

When operating upon a single centroid list at a time, the user can optionally plot the drift curve determined by the algorithm; or, using the search functionality of the Windows open file dialog, one can search the main directory for all centroid lists just created by L1 digestion and then CCDLAB will process and apply the drift determination to each centroid list. A cancellable wait bar will show the progress. The procedure will also write FITS images of the drift-corrected data with the total integration time for the image written to the header, and at that point one basically has a science-ready image aside from normalization to the integration time and filter-dependent factors applied to transform counts/second into whatever other units are desired (Tandon et al. 2017).

For the 42 unique data sets which came out of the L1 digestion of the NGC 1851 data, drift determination and correction for all sets required a total of 8 minutes to process, including the creation of flat-corrected images of 1/8th pixel resolution (4096 × 4096 pixel images at double precision).

The user can re-run the drift-determination algorithm on the drift-corrected list to tighten up the PSF of the sources a little bit more. Simply select the new drift-corrected centroid lists from the same open file dialog for drift correction. In our example of NGC 1851 data, this required another 8 minutes to process. This is called second-order drift correction and it essentially is an attempt to pick up any residuals which may remain in the field position relative to the drift series calculated for a given position; this procedure has reliably been seen to reduce the PSF in processed images. If the first drift series was calculated with a 2-second cadence, then choosing a disharmonious cadence of 3 or 5 seconds for the second order may help to find residuals in the drift series, although simply running it at 2 seconds again has also been found to help. Third- and higher-order drift corrections can be performed as needed, but typically it is not necessary unless one is dealing with a particularly difficult (faint) source. Typical final PSF results are in the range of 1.1'' for NUV and 1.3'' for FUV with this method.

For some targets, the flux in FUV and/or NUV channels is too low and the channel will not be able to track its own drift. If a larger time stack is used for tracking to gather more signal in each stack, then the resolution of the spacecraft drift pointing becomes degraded as the point sources get smeared through the running image stack, resulting in poor drift determination and thus degraded or useless final science images. Thus, VIS tracking of point sources can be used and the drift series determined from the VIS channel can be transformed and used for the NUV and FUV channels. The VIS channel is always expected to have at least one point source for a target's observation.

The CCDLAB UVIT Pipeline can determine the VIS drift series in two ways. The first is an automated routine that determines the x and y lags between the reference frame (first frame) and subsequent frames by the cross-correlation method discussed earlier. This routine requires at least one strong point source in the VIS images, which does not fall out of the field of view (due to drift) during the observation.

The second routine uses point-source tracking of one or more point sources in the image, and these point sources can be tracked at a much weaker source intensity. The routine is interactive, with the user selecting the point sources in the reference image, and then the routine continues to automatically track those point sources in subsequent images, gathering the mean differential in point-source positions relative to the reference frame to determine the drift series. Again, the point sources must not leave the field.

Once the drift series is determined for the tracking channel, then the series can be transformed and applied to another channel by the pipeline. One can repeat the drift-tracking procedure again to get a second-order correction, but this time using self-tracking at a stacking time of, say, 20 seconds if the source is quite faint. The VIS channel correction will have gotten most of the high- and low-frequency drift shifts, but performing a second-order self correction can pick up whatever residuals it may be able to detect. The VIS-only correction will give PSF results in the same range as with self correction, but when VIS correction and then second-order self correction has been able to be performed this is where PSF results down to 0.9'' for NUV and 1.1'' for FUV have been seen.

The frame rate sets a limit on the source photon rate detectable without approaching saturation, although statistical corrections can be made to photon rates which enter into saturation of the detector read rate. CCDLAB provides a saturation correction utility for point sources in reduced images. The saturation correction is made by examining the centroid list in the graphically specified local region of the saturated source, and counting the number of frames where no photons were detected in that region. Given the measured probability of detecting zero counts per frame read, the Poisson probability equation can be solved for the actual count rate as the natural log of the inverse of the probability of detecting zero counts per frame.

For example, in an observation of the white dwarf HZ4 during PV phase, the nearby star TYC 659-370-1 gave a probability of no photon arrival per frame of 0.418. This gives an actual photon rate of ln(1/0.418) = −ln(0.418) = 0.872 per frame, in comparison with a detected photon rate of 0.598 per frame. This corrected value would then be used, for example, for the photometry and/or flux calibration as in Tandon et al. (2017).

The uncertainty on desaturated photometry is as follows. For an actual rate of "X/frame" photon events, the observed number is N(1−exp(−X)) for N frames. The variance of the observed counts is then N*exp(−X)*(1−exp(−X)). The rate of change of observed counts with the actual counts is d(N(1−exp(−X)))/d(NX) = exp(−X). This relation can be used to relate estimates of rms variation of the observed counts to the rms variation of the actual counts and write: rms variation on the actual count of incoming photons (i.e., N X) = (N* exp(−X) (1−exp(−X)))o 0.5/(exp(−X)) =(observed counts/exp(−X))o 0.5.

A region of interest can be graphically selected from which to extract the corresponding centroid data only, thus one can generate the centroid list for a particular target. This centroid list can be used for timing analysis of the centroid arrival times for the object selected in the region of interest.

Given that these times should be absolute times and not just detector clock times, the L1 data provides a file from which one can correlate detector clock times to ground-based UTC time. CCDLAB uses this to then transform the ground-based time into Barycentric Julian Day (BJD) times for each centroid for the target's coordinates, using the methods from the 2005 Astronomical Almanac which are stated to be accurate to about 0.1 minutes. More accurate methods will be incorporated in a future release via the International Astronomical Union's Standards of Fundamental Astronomy (SOFA) software library. The BJD times are generated as part of the main L1 file data reduction, and are extracted with the centroid data for a particular target in this method.