3 Visualising data

In this section several procedures are described for looking at SCUBA-2 data, as well as basic data reduction steps that can be run separately. Working through these examples will illustrate some of the features of SCUBA-2 data, but will not result in a science-grade image at the end. If you are interested only in making the best possible map with minimal effort proceed to Section 4.

3.1 Concatenating data

Since SCUBA-2 data for a given subarray are broken into many pieces by the DA system, it is useful for visualisation to first concatenate the data into single files. The Smurf task sc2concat can be used for this operation. For example, assuming the following files are in the current working directory,

combines all of the files associated with observation 15 for the s4a array into a single file called s4a20091214_00015_0002_con.sdf. sc2concat will automatically filter out any dark observations, so that the concatenated file contains only the data taken during the scan across Uranus with the shutter open. It also applies the flatfield by default, although it can be disabled using the ‘noflat’ option on the command-line. Be careful when concatenating a very long observation since the output file may be too large to handle sanely. Fifteen minute chunks (30 files) should be more than big enough.

3.2 Displaying scan patterns

The pointing of the telescope throughout a scan (as well as other state information) is stored in the MORE.SMURF.JCMTSTATE extension of a data file. The Smurf task jcmtstate2cat will convert this information into a simple ASCII tab separated table:

Figure 1: The telescope positions during observation number 7 on 20090107. The plot is created by Topcat from the output from jcmtstate2cat plotting the DRA column (right ascension offset from the map centre in arcsec) against the DDEC column (declination offset from the map centre in arcsec). The PONG scan pattern is clearly visible (see Scott & Engelen 2008[5]).

The ‘-h’ option to jcmtstate2cat can be used to find more information on the command. In particular, multiple files can be supplied to the command using standard shell wild cards (not escaped) and for SCUBA-2 data the ‘–with-mce’ option can be used to dump the low-level MCE header information.

This catalogue can be loaded into Topcat for plotting, making sure that Topcat is told that the TST format is to be used for reading.

An example plot of the scan pattern for this observation generated by Topcat can be seen in Fig. 1. All the time varying header values are available for plotting. In particular DRA and DDEC will show the RA/Dec offset of the telescope, DAZ and DEL will show the Az/El offset and the 225 GHz opacity values are also calculated from the raw WVM measurements.

3.3 Displaying data cubes

 

Figure 2: Initial Gaia windows displayed upon loading a data cube (with slight modifications to the colour table using options under the ‘File’ and then ‘Startup options’). Left: The main window, after clicking the ‘Z’ button a number of times to zoom-in, shows a map of bolometer values at a fixed sample in time. Note that these data came from an array with a number of broken columns and broken isolated bolometers, all indicated in grey. Right: The ‘Display image sections of a cube’ dialogue enables the user to navigate the time dimension. The ‘Index of plane’ slider near the top can be used to select different time slices, and the main window will automatically update.

The easiest way to visualise the bolometer time series data is to use Gaia. Loading in the concatenated file above (combining the two example files included with this Starlink release) produces two windows (Fig. 2). The main window shows a map of bolometer values at a given instant in time. The second window can be used to navigate the time axis; by moving the ‘Index of plane’ slider in the ‘Display image sections of a cube’ dialogue, different time slices may be selected, with the main Gaia window updating automatically.

For this concatenated data file (1 minute in total, as it is the combination of two 30 s files), each bolometer has a large offset relative to its neighbours as well as relative to any smaller time-varying signals, which means that little difference can be seen by moving the slider. However, Gaia can produce an automatically scaled plot of the time series for an individual bolometer by simply clicking on it in the main window. For example, clicking on the bolometer at (13,23) spawns the ‘Spectral plot’³ window shown in Fig. 3 for a single bolometer. Clicking on other bolometers over-writes the plot of the original bolometer in the same window. Looking at the vertical axis range on the left, the mean levels clearly vary significantly from bolometer to bolometer, although the time-varying component of the signals are quite similar.

Figure 3: The ‘Spectral Plot’ window is spawned automatically once a bolometer is clicked in the main window, such as (13,23) in this example, displaying its time-varying signal. The vertical red line indicates the time slice that is currently selected in the ‘Display image sections of a cube’ dialogue. The regular pattern, with a period of about 30 seconds in this particular case, is slow baseline drift correlated with variations in the SCUBA-2 fridge temperature. The green-circled narrow spike is produced by the bolometer crossing Uranus. Usually astronomical signals are not readily visible in raw time series plots such as these, since they are relatively much fainter. Note that on some systems ‘gaps’ may appear, such as those at the start of the time series, and at $\sim 34$ s in this example. These are simply rendering artefacts due to re-sampling the data to the resolution of the display (zooming-in to these regions using the ‘Lower index’ and ‘Upper index’ sliders in the ‘Display image sections of a cube’ dialogue followed by a click on ‘Re-extract’ to update the plot demonstrates this).

3.4 Regridding data into a map

A simple and quick map can be made from a data cube using the Smurf makemap task. The following will produce a map directly from the raw concatenated data by re-gridding it into a pixelated map:

The makemap task automatically scales the bounds of the image to encompass all of the data. The output map here is called ‘uranus.sdf’, and the pixel scale is 2 arcsec on a side by default at 450 $\mu$m and 4 arcsec at 850 $\mu$m (although this can be changed using the ‘pixsize=$x$’ option on the command-line, where $x$ is in arcsec)⁴. Since we already know that relative bolometer offsets are large in these data, it is unsurprising that no astronomical source can be seen with Gaia in the resulting image.

Figure 4: Map produced from raw data of a scan across Uranus. The image is completely dominated by noise in the relative signal offsets of each bolometer, and no astronomical signal can be seen. The scan (a Curvy PONG) can clearly be seen as a repetitive ‘waffle’ pattern in the image.

3.5 Cleaning data

The previous examples illustrate the need for some kind of data cleaning before there is any hope of seeing astronomical sources. A useful Smurf task for time-domain data processing is sc2clean, which can perform several different steps controlled by a range of parameters. Note that all of the algorithms available to sc2clean are also accessible within the DIMM (Section 4), so in practice the user does not issue this command directly before making the final map.

3.5.1 Removing bolometer averages and DC Steps

The following will remove a 0th-order polynomial (the mean) from each bolometer time stream, storing the cleaned data in a file called ‘clean.sdf’:

Figure 5: Top: Gaia plot of the array signal at time slice 1 after using sc2clean to remove the mean of each bolometer signal (clean.sdf). The LOW and HIGH data values for this intensity plot are $-1.7171$ pW and $+5.50352$ pW, respectively. Most of the working bolometers (not grey) now have approximately the same signal range (approximately $\pm 0.1$ pW), and therefore have nearly the same colour (brown). However, there are several outliers, for example, the orange pixel with an ‘x’ at (10,17). Bottom: The time series for bolometer (10,17) showing the presence of a large level change, or ‘DC Step’ near 30 s.

This operation results in bolometer data that lie primarily in the range $\pm 0.1$ pW, except for a handful of outliers as shown in Fig. 5. In these cases there have been large level changes, or ‘DC Steps’, during data acquisition. These can be identified and repaired using some extra parameters to sc2clean:

In addition to removing the mean as in the previous example, the extra parameters starting with ‘dc’ instruct sc2clean to smooth the data with a median filter of width 50 samples, and then look for steps in the smoothed data in excess of 25-$\sigma$. The heights of steps found in this way are measured by fitting straight lines to the smoothed data on either side of the jump using 30 samples in each fit. The corrections thus determined are applied to the original bolometer data. Note that an additional flag has been set by default, fillgaps=1, indicating that the data range around the identified steps ($\pm 400$ samples in this case) should be replaced with a constrained (smooth) realisation of noise to avoid introducing spikes into the data stream (these ranges are ultimately ignored when producing final maps). A map produced from ‘clean2.sdf’ is shown in Fig. 6.

Figure 6: Map produced from ‘clean2.sdf’ in which DC steps have been repaired and the bolometer means have been subtracted. There are still many artefacts in the data, but now Uranus is visible.

The command line can get long, so it is also possible to write the configuration parameters into a text file and specify the text file for CONFIG using the standard group notation (‘myconfig.lis’ in the following example; the leading ‘^’ is necessary, but it is not part of the file name):

An example config file containing the default values for each item can be found in
$SMURF_DIR/smurf_sc2clean.def.

3.5.2 Watching a movie

Gaia has the ability to animate the display of a data cube. We will use this feature to make a ‘movie’ of the array data. Load ‘clean2.sdf’ from the previous step into Gaia. In the ‘Display image sections of a cube’ dialogue, switch from the ‘Spectrum’ to the ‘Animation’ tab approximately half-way down. Set ‘Delay’ to 10 milliseconds (the smallest value), ‘Step’ to 5 (such that it only shows 1 in every 5 frames), and click the ‘On’ button next to ‘Looping’. Finally, click ‘Start’, and an animation of the data cube will be shown in the main Gaia window. The dominant signal is a gradual variation in the average value of all of the bolometers in unison. This common-mode signal is produced through a combination of SCUBA-2 fridge temperature variations, sky noise, telescope motion, and other drifts in the individual bolometers. At times it is also (just barely) possible to see Uranus itself as the array scans across is.

3.5.3 Frequency-domain filtering

sc2clean can also perform frequency domain filtering on the bolometers. In the following example data are concatenated, and filtered using a single call to sc2clean.

Figure 7: High-pass filtered data (clean3.sdf). Top: The time series for a single bolometer (13,23). Comparing with the un-filtered data in Fig. 3, the positive spikes produced by Uranus are now much more apparent. Bottom: a map produced from these cleaned and filtered data. The brightness of Uranus compared to the noise is now much more significant than in Fig. 6, but the filtering has caused ringing around the source-crossing in the time series, which appears as negative dips in the map along the scan directions (i.e. a negative cross around the source).

Here the text file, sc19_clean3.lis, contains the following configuration options:

As in the early examples, the first four parameters cause sc2clean to remove the mean bolometer values, and repair DC steps. Next, the apod=<undef> (combined with an internal default zeropad=0) is a slightly cryptic way to tell sc2clean to enforce continuity between the starts and ends of each bolometer time-series by temporarily adding padding that is a function of the filter frequency, and interpolating the ends using a cubic function. This step is required to avoid ringing once the FFT is taken. The alternative, and now deprecated, method for reduing ringing is to pad and apodize. Finally, the last parameter tells sc2clean to apply a high-pass filter with a hard lower cutoff frequency of 0.5 Hz.

Fig. 7 shows the filtered bolometer signal for (13,23), as well as a map produced from the data. The filtering has significantly improved the noise properties of the map compared with Fig. 6, but the filtering has now introduced ringing around the source which results in a negative cross pattern along the scan directions.

3.6 Bolometer power-spectra

The frequency-domain power spectra of SCUBA-2 bolometers can be produced with the Smurf task sc2fft. Similar to the previous example, we first produce cleaned, concatenated data. However, we omit the high-pass filtering by overriding the ‘filt_edgehigh’ parameter explicitly, since the goal now is to see what the full bolometer noise power spectrum looks like.

While there is a Starlink standard format for storing complex values as described in [2], sc2fft uses its own local format that allows both for Cartesian and polar forms: a 4-dimensional array, with the first axis indicating frequency, the second and third axes bolometer row and column, and the final axis containing the real and imaginary parts of each Fourier coefficient. By setting ‘power=true’ it switches to polar power form, such that the first element of the 4th array axis stores the square of the amplitudes, and the second element stores the arguments (phases) of each Fourier coefficient. As with the time series data, Gaia may then be used to view the data cube. The difference is that we now specify a subset of the data so that the 3rd axis of the cube is the array of squared Fourier amplitudes for each bolometer, and ignore the phases (i.e. we observe only the 1st elements along the 4th axis):

It is then possible to click on each bolometer to display its power spectrum. Once the ‘Spectral plot’ window has been spawned, it will also be necessary to modify the axis displays. Select ‘Options’$\to$‘Positive Y Only’, and ‘Options’$\to$‘Log Y Axis’. Similarly, select the corresponding settings for the ‘X’ axis. Fig. 8 shows the power spectrum of bolometer (13,23), with a $1/f$ knee apparent close to 1 Hz.

Figure 8: The power spectrum of bolometer (13,23) produced by the Smurf task sc2fft. In this particular case the noise looks flat above a $1/f$ knee of about 1 Hz. The red line simply indicates the frequency slice currently being displayed in the main Gaia window (not shown).

3.7 Data quality flagging

NDF files manipulated by Smurf use the standard Starlink named Quality mechanism (see discussion of Masking, Bad Values, and Quality in SUN/95). Quality itself is stored as an 8-bit integer for each sample in a data file, and each bit can reflect a different condition. For example, the following will indicate the number of samples flagged by sc2clean when producing ‘clean4.sdf’ in the previous example:

The ‘count’ supplied to showqual indicates that the total number of occurrences of each Quality flag in the data should be counted and displayed. This example shows that 4416000 samples in total were flagged ‘BADDA’ by the flatfielding algorithm or by the data acquisition (DA) system itself based on the array setup. Since each bolometer produced 12000 samples, this corresponds to 4416000/12000 $=$ 368 bolometers that were not operational (out of 1280). The ‘BADBOL’ flag is used by Smurf to flag every sample of a bolometer that is not being used. In this case the number of flags is 4476000 which corresponds to 373 bolometers. In other words, cleaning turned off an additional 5 bolometers. Finally, ‘DCJUMP’ indicates the number of samples that were affected by the DC step correction procedure (which include a number of samples on either side of the precise locations where the steps occurred).

During map-making flagged samples will not be used and by default they will be hidden from view for all Starlink tools. However, if you wish to see the data that were flagged you can use the Kappa setbb command to enable a particular flag or turn them all off (so that they can be seen in Gaia).

will disable quality and make visible all the underlying data. Using a value of 255 will turn all bits back on and so any non-zero quality will be treated as bad data.

Additionally, the Kappa task qualtobad may be used to permanently convert samples with given quality bits to a magic or invalid value. For example:

will set all of the samples flagged with quality DCJUMP to the magic value. When ‘badmasked.sdf’ is subsequently viewed with Gaia, none of those portions of the data will be visible.