How Should I Mosaic My SCUBA-2 Data (Redux)?

Mosaicking your SCUBA-2 with wcsmosaic or makemos is fine if all you're interested in is the data. But what if you want to know the exposure time per pixel in your map to determine if the noise is reasonable? That's easy for each map written by makemap: the exposure time image is stored in the .MORE.SMURF.EXP_TIME NDF extension (which you can view in Gaia). But wcsmosaic doesn't know about this extension, so the EXP_TIME data in a mosaic is only that of the first file.

There are two ways to get properly mosaicked data:

• Use the ORAC-DR pipeline and process from raw;
• Use the PICARD recipe MOSAIC_JCMT_IMAGES on processed data.

I'll write a post about using the pipeline in the near future but for now I'll highlight the second method.

PICARD comes for free with ORAC-DR (which comes for free with Starlink). It's basically a tool for processing and analyzing reduced data which takes advantage of the same infrastructure used by the pipeline.

At the command line type:

picard -log s MOSAIC_JCMT_IMAGES myfiles*.sdf

It'll tell you what it's doing and exit. This will create a file called "mylastfile_mos.sdf" which has the correct EXP_TIME extension (where "mylastfile" is the name of the last file in the list given to MOSAIC_JCMT_IMAGES). Try it and see.

Note that all the files must be of the same source.

At the moment, it only supports wcsmosaic as the mosaicking task but will support makemos in the future. Keep rsync'ing and one day it'll just be there...

[This post was updated on 20100510 to reflect the change in the recipe name and the file-naming behaviour]

How Should I Mosaic My SCUBA-2 Data?

I've been asked this a number of times so I thought I'd write something down. Currently the SMURF iterative map-maker calculates models for each chunk of data independently (a chunk is a continuous sequence so many times that will be a single observation but it might be smaller than an observation if a flatfield is inserted mid-way through). It then combines them using a weighted coadd to give you the final map. This means that there is nothing to be gained in throwing all your observations at MAKEMAP apart from it taking much longer to make a map before you can see the result.

In terms of flexibility it is better to make a map from each observation separately. You can then look at each map to see how it looks before combining it. You may also want to remove any low frequency structure at this point as well. To combine these maps into a single map you have two choices using Starlink software:

1. Use KAPPA WCSMOSAIC with an interpolation or rebinning scheme of your choice. Make sure that you set the VARIANCE parameter to true to enable variance weighting of your mosaic.
2. If you have already ensured that the maps are made on the same grid (you can use the REF argument in MAKEMAP to ensure this) then you can try CCDPACK MAKEMOS. This can be used to do median or clipped mean stacking and does no interpolation or smoothing. If your images are not aligned on the same grid it is probably better to remake them so that they are, but if that is difficult you can use the KAPPA WCSALIGN command to match up the grids first. Use the GENVAR and USEVAR parameters in MAKEMOS.

In the future it may be that the map-maker will be able to do a better job handling all the data at once itself instead of external mosaicking but at the moment this is not the case.

If you want to keep a track of exposure times you need to track that separately. We will cover that in a follow up post.

Large scale structure divergence

First, my apologies for the lousy image layout! Next post will be better...

As many of you have noticed, large-scale noise structures in the maps appear to grow with the number of iterations (clearly not the desired behaviour!). This blog entry describes the problem in more detail with an example, and discusses a partial fix.

The map on the left was produced from about 6 min. of scans across a blank-field at 450um using the dimmonfig_faint.lis configuration file, but turning up the number of iterations to 100. As you can see there is a strong gradient across the map -- the more iterations I added, the more the gradient grew.

In the next plot I show the common-mode signal from this map solution (the average signal in all the detectors as a function of time). Mostly we see the fridge temperature oscillations (about every 30s), plus a more slowly varying component (maybe sky?)

Next, I set ast.zero_lowhits = 0.5 in the config file. What this option does is identify all of the pixels in the output map with less than 50% of the mean number of hits, and then sets those pixels to zero after each iteration (i.e. a zero boundary condition). In the next two plots I show the map, and the difference between the common-mode estimated from the two solutions:








What we can can see is that the gradient has been drastically reduced. So what's happening? The common-mode difference plot shows that the map without edge constraints has a large additional periodic signal (about 10% of the main common mode signal) with a frequency that is a factor of a few larger than the 30s fridge oscillations. In fact, this signal is nicely correlated with the higher-frequency scan pattern.

Basically this shows that, left to its own devices, the solution is perfectly happy to pump large-scale structure into the map. It then still manages a flat residual because it simply compensates for this gradient by putting negative signal of the correct amplitude into the common-mode. In other words, the largest scales and the common-mode are degenerate model parameters.

However, I have noticed that even with this edge constraint turned on, there is other strange structure in the map, like this convex curvature. This fact, combined with the substantially worse noise at 850um, shows us that the model is insufficient for accurately describing the data. The edge constraint helps, but that curvature indicates some kind of tension (like the model trying to fit some sort of gradient across the array and failing). I suspect that, lacking any other conditions on the model, the large-scale noise is just growing without bound.

We are presently looking at ways to improve the model for the data and/or introduce other reasonable constraints to help convergence to something sensible! For now, it is worth turning ast.zero_lowhits on if you don't believe you have any strong emission features around the edge of the map (tricky for maps of structure in the galactic plane).

Detection and Correction of DC steps

Sudden changes in the DC level are seen to occur within bolometer time streams. The following algorithm attempts to identify and correct these steps. Each bolometer time stream is processed independently using the following two part process

Part 1 - Detecting Candidate Steps:

In this part we search the bolometer time stream for contiguous blocks of samples that have unusually large gradients compared to the neighbouring samples. Each such block of samples is considered to be a candidate DC jump. The starting and ending time of each such block is noted. Once all candidate jumps have been found for a bolometer, the algorithm proceeds to part 2. Candidate jumps are found as follows:

1. Smooth the bolometer time stream using a median filter of width DCSMOOTH (defaults to 50) samples. This reduces the effects of the noise in the data. Using a median filter rather than a mean filter helps to preserve the sharp edges at each DC jump, and is robust against spikes, etc. The following  figure shows a typical step, with  the original bolometer data values in red and the median-smoothed data values in blue.
   
2. At each sample, find the difference between the two adjacent values in the median filtered bolometer data. This is proportional to the gradient of the median filtered bolometer data. These gradients include both the gradient of the underlying noise free median-filtered data stream ,and the extra gradient caused by noise, spikes, DC jumps, etc. The following figure shows these differences for the above step. The high gradients at the step are clearly visible.
   
3. We now need to split these differences up into differences caused by variation of the underlying  noise-free data, and differences caused by noise, steps, spikes, etc.To do this we smooth the above differences, and use the smoothed differences as the gradient of the underlying noise free data. In order to get an accurate result, we need to do this smoothing with a mean filter rather than a median filter. But a mean filter is badly affected by aberrant values, so we first take a copy of the differences found above, and remove extremely large differences by flagging them as bad. Specifically, differences more than three times the RMS difference (taken over the whole time stream) are set bad. The RMS is then re-calculated and a further two such rejection iterations are performed. This gets rid of  high gradient values caused by spikes, jumps, point sources, etc, leaving the following:
   
4. The differences close to a step are often not typical of those in the wider neighbourhood, so each block of contiguous values removed above is doubled in width, leaving the following:
   
5. Smooth the remaining differences using a mean filter of width DCSMOOTH.  The smoothing process fills in small holes in the array (up to 0.8*DCSMOOTH), but will not fill in larger holes. Such larger holes are filled in by linear interpolation. The resulting smoothed differences - shown in blue below overlayed on the original differences in red -  estimate the local gradient of the underlying noise-free median-filtered data stream.
   
6. Subtract the smoothed differences from the total differences found in step 2 above. This gives the residual differences caused by noise, spikes, DC jumps, etc., shown below:
   
7. To decide if a residual difference value is "unusually high" we need an estimate of the local RMS of these residual differences at each point in the data stream. To get this we smooth the squared residual gradients using a median filter of width DCSMOOTH2 (fixed at 200), and then take the square root of the resulting values.
   
8. Divide the residual differences by the local RMS of the residual differences to get a "signal-to-noise" ratio for each residual difference value:
    
9. Find all blocks of contiguous absolute SNR values greater than DCTHRESH (defaults to 25). Such blocks are allowed to contain small numbers of values that do not meet this criterion. The maximum length for such a section of low SNR is DCFILL (fixed at 40) samples. For the example step shown above, the block starts at 13696 and ends at 13703. Blocks are ignored if they fail any of the following tests:

• The width of the block must be less than DCMAXWIDTH (fixed at 60)
• The total of all the SNR values in the block must be at least DCTHRESH. This test avoids spikes being treated as two very close steps, since such steps will  have SNR values of opposite signs and will thus cancel out.
• The total of all SNR values in the block must exceed half the value of the largest single SNR value in the block.
•  (this tends to rejects spikes since they will have large positive and negative SNR values within the block, which will cancel out).

Part 2 - Measuring and Correcting Steps

In this part, each of the candidate steps found above is processed to determine the height of the step.  This is done by doing a least squares linear fit to the median-smoothed bolometer data just before the step, and another fit to the data just after the step. These two fits are used to predict the value at the centre of the step, and the difference between these two predicted values is used as the step height. Various other tests are performed to identify and reject spurious step detections, and finally the original bolometer data stream is corrected to take out the step. For a given bolometer, each candidate step is processed as follows:

1. Correct the step start and end to include adjacent SNR values that are below DCTHRESH but are still significantly high. Specifically, the step start is moved to the end of the last block of DCNLOW (fixed at 5) SNR values that are all less than DCSIGLOW (fixed at 8.0), and which occurs before the original step start. Likewise, the step end is moved to the start of the earliest block of DCNLOW SNR values that are all less than DCSIGLOW, and which occurs after the original step end.For the example step shown above, the start moves from 13696 to 13693, and the end moves from 13703 to 13711. 
2. Perform a linear fit to the median-smoothed bolometer data before the step. The fit is over a box containing DCFITBOX (defaults to30) samples. A gap of DCFITBOX samples is left between the end of the box and the start of the step.
3. Perform a linear fit to the median-smoothed bolometer data after the step. The fit is over a box containing DCFITBOX samples, and a gap of DCFITBOX samples is left between the end of the step and the start of the box.
4. These two fits are evaluated at the centre of the step. The difference between these values gives an estimate of the step height. However, moving the fitting box slightly can often result in big changes in the line gradient, and thus big changes to the estimate of the jump height. Steps where this is the case are ignored (i.e. left unfixed). To check this, more linear fits are performed at increasing distances from the step, and the consistency of the estimated jump heights from these extra fits is checked. Specifically, the fitting box used in step 2) above is moved one sample earlier, a new  fit is performed, and the new fit is used to estimate the value at the centre of the step.  The fitting box is then moved again by one sample and a third fit performed giving a third estimate of the central data value. This process is repeated until we have performed 2*DCFITBOX fits (an efficient algorithm is used to perform these multiple fits). The variance of all these estimates of the central value is found. The same thing is then done to the fit performed in step 3) above, to give us the variance of the central data value determined from the data following the step. The uncertainty in the jump height is then taken to be the square root of the mean of the two variances values, and the jump is ignored if the jump height - determined from steps 2) and 3) - is less than DCTHRESH2 (fixed at 1.5) times the uncertainty in the jump height.
5. The jump is also ignored if it is less than DCTHRESH3 (fixed at 4.0) times the noise in the original bolometer data. The bolometer noise value used is the RMS of the differences between adjacent bolometer values, excluding aberrant points, and corrected by a factor of 1/sqrt(2).
6. If the step is not rejected,  all the original bolometer values following the centre of the step are adjusted upwards or downwards by the step height.

Once all steps have been corrected, a constant value is added back on to all bolometer samples so that the original mean value in the bolometer is retained.