Recovering the number of chunks used by makemap

This post has been superceded - see the new post "Recovering the number of chunks used by makemap - take 2".

As of today, the number of chunks of continuous time-series data used by makemap is recorded in the output map. It is stored in the NCHUNK item within the SMURF extension of the output map. So for instance, it can be viewed by doing:

   % kappa
   % setext mymap.sdf smurf get nchunk loop=no

If NCHUNK has a value of 1, then all the input time-stream data was processed as a single chunk. If NCHUNK is larger than 1, then either the computer had insufficient memory to process the data as a single chunk, or the data was not taken as a single continuous observation.

An rsync of the Hilo stardev system is required to use this new feature.

A new Starlink release contains notable updates to the SCUBA-2 configuration files.

The latest Starlink release - 2014A - has been made public. For details please read the release notes provided at: http://starlink.jach.hawaii.edu/starlink/2014A

As part of this new release we want to highlight one significant update and a couple of new additions to the SCUBA-2 reduction arsenal of config files.

Updates to bright_extended config file

This config file 'dimmconfig_bright_extended.lis' has always been intended for reducing data containing bright extended sources. It has remained untouched for a couple of years now despite advances in our understanding of SCUBA-2 reduction of bright regions. The config file now contains the following parameters/links:

   ^$STARLINK_DIR/share/smurf/dimmconfig.lis

   numiter=-40
   flt.filt_edge_largescale=480
   ast.zero_snr = 3
   ast.zero_snrlo = 2

   ast.skip = 5
   flt.zero_snr = 5
   flt.zero_snrlo = 3

In previous Starlink releases i.e. Hikianalia, the bright_extended configuration file only contained the following:

   numiter = -40
   ast.zero_snr = 5

   flt.filt_edge_largescale = 600

New 'FIX' config file

Two new config parameter files have been added These are intended to be used with one of the existing dimmconfig files. They provide new values for selected parameters aimed at solving a particular problem ("blobs" in the final map, or very slow convergence).

• dimmconfig_fix_blobs.lis 
• These parameters attempt to prevent smooth bright blobs of emission appearing in the final map. It does this by 1) identifying and flagging samples that appear to suffer from ringing a soft-edged Butterworth filter in place of the normal hard-edged filter, and 3) rejecting samples for which the separate sub-arrays see a markedly different common-mode signal.

• dimmconfig_fix_convergence.lis
• The parameters defined by this file attempt to aid the convergence process, and should be used for maps that will not converge within a reasonable number of iterations.

Spotting chunking caused by insufficient memory

It has been shown that better maps are created if all available time stream data are processed in a single chunk to create a single map. If the available memory is too small to allow this to happen, the data will be split into multiple chunks and a separate map made from each chunk, which are all co-added at the end to make the final map.

Thus it can be important to spot whether or not your map was "chunked" or not. Until recently the only option was to trawl through the makemap screen output looking for the warning about insufficient memory. Now, however, you also have an option to tell makemap to abort immediately with an informative error message if the data would be chunked due to insufficient memory (no output map is created in this case). To do this, add "memcheck=1" to your config.

Note, this only affects what happens if the data is processed in chunks due to lack of memory - if chunking is caused by the time-series data being intrinsically discontiguous then a map is still created even if memcheck is set to 1.

It will be necessary to rsync the stardev system from Hilo to use this feature. See starlink.jach.hawaii.edu/starlink/rsyncStarlink

Calibrating SCUBA-2 data in surface brightness units

In most cases the default calibration for SCUBA-2 data processed by the ORAC-DR pipeline is mJy beam^-1. The exception is the recipe for extended sources, REDUCE_SCAN_EXTENDED_SOURCES, which calibrates data in mJy arcsec^-2.

Unfortunately there was an error in an earlier version of this recipe which meant that the FCF was applied incorrectly. The corrected method is available now with an update of ORAC-DR (either from github or via rsync from JAC). If you have data processed with this recipe (either by running it yourself, or downloading processed products from CADC) then re-calibrating the data is easy: simply divide by the pixel area using KAPPA cdiv.

There is a new PICARD recipe for easy calibration of maps produced by running makemap by hand. CALIBRATE_SCUBA2_DATA allows data to be calibrated in in per-beam and surface brightness units. With no parameters, this recipe will calibrate data in mJy beam^-1. For surface brightness calibration, set the recipe parameter USEFCF to 1 and FCF_CALTYPE to ARCSEC, and the recipe will then use the default ARCSEC FCF for the wavelength of the given data.

The recipe can also convert the calibration from one type to another. If your data are already calibrated in mJy beam^-1, they can be given to CALIBRATE_SCUBA2_DATA with the FCF_CALTYPE recipe parameter above, and the recipe will create a new file (with suffix _cal) with units of mJy arcsec^-2. The value and units of the FCF are written into the FITS header of the calibrated file.

The companion recipe, UNCALIBRATE_SCUBA2_DATA, will undo the current calibration, reverting the units to pW in the output file (which has a suffix of _uncal).

Using ORAC-DR or PICARD to perform the (un)calibration is preferred to simply multiplying your data by the FCF as they also set the units correctly for the output files(s), and write the value of the FCF used into the FITS header of the file.

However, there is one note to highlight: the recommended way to calibrate data (either from raw or when changing from per beam to per square-arcsec) is to calibrate the individual observations first, and then coadd those (re)calibrated files. Calibrating or re-calibrating coadds will fail because the coadding step was recently updated to remove FITS header entries that differ between the input files. These usually include the UTDATE which is used by the ORAC-DR calibration system. A future upgrade will provide a workaround though the recommendation to calibrate individual observations stands.

Change to FLT Ringing Filter

If you use the AST.SKIP and FLT.RING_BOX1 configuration parameters to prevent artificial blobs of emission appearing in your SCUBA-2 maps, then you may have seen cases where  the centres of real bright sources are badly affected, as in the following example spotted by Harriet:

This issue may also cause complete holes to appear around bright sources. These problems are caused by the fact that the ringing filter is applied on the initial "pre-iterations" specified by AST.SKIP, as well as on all subsequent "real" iterations. If the ringing filter is restricted to just the "real" iterations, the following map is produced (left uses the same scale as above, right is a close up of the central region):


The stardev system at Hilo has been updated to include this fix.

For those interested in the details... The problem is caused by the fact that no FLT masking can be done on the first pre-iteration, since there is no map available at that point to use as a mask. Therefore bright sources cause heavy ringing on the first iteration. This ringing is then detected and the corresponding samples flagged by the ringing filter. This means that the bright sources are blanked out in the map created at the end of the first iteration. On the second pre-iteration, this map is used to locate the bright sources that should be masked out when estimating the FLT model - but the bright source is not present in this map and so is not masked out on the second iteration either, leading to heavy ringing again. And round it goes doing the same thing on each iteration. Restricting the ringing filter to "real" iterations solves the problem because bright sources have been masked out in the FLT model by the time the ringing filter is used for the first time, and so have not caused the heavy ringing that would otherwise have been detected and flagged by the ringing filter. In fact, in may have been sufficient just to inhibit the ringing filter on the first pre-iteration, rather than on all pre-iterations, but this would have been of little benefit since blobs only start to develop once the AST model is included in the iterative scheme (i.e. once all the pre-iterations have completed).

Interpreting SCUBA-2 map Quality Arrays

SCUBA-2 maps created using "AST-masking", such as is done by  dimmconfig_bright_extended, will contain a Quality array indicating the background pixels that were masked (i.e. forced to zero). In itself this is fairly simple - background pixels have a non-zero Quality value and source pixels have a Quality value of zero. However, it is also possible to have independent masks for the FLT and COM models, in addition to the AST mask. For instance, if the "auto-FLT masking" scheme is used the Quality array contains both a FLT and an AST mask and so some care needs to be used when interpreting the Quality array.

Each value in the Quality array is restricted to taking integer values between 0 and 255, and so can be thought of as 8 separate bits. Each "bit plane" within the Quality array holds a single mask - AST, FLT or COM.   The showqual command can be used to find out which mask is held by which bit plane:

% kappa
% showqual fred.sdf
   AST (bit 1) - "Set iff AST model is zeroed at the output pixel"
   FLT (bit 2) - "Set iff FLT model is blanked at the output pixel"

This means that the AST mask is stored in bit 1 (the least significant bit), the FLT mask is stored in bit 2, and there is no COM mask. Note, if a map was produced using FLT masking but no AST masking, then the FLT mask would be stored in bit 1.

The decimal integer value of any element of the Quality array is equal to the binary value formed from the bits listed by showqual. So in the above case the maximum Quality value is 3 (the decimal equivalent of binary "11" - both bits set). Remembering that a bit is set (i.e. is 1) for background pixels and cleared (i.e. is 0) for source pixels, it follows that the four possible decimal Quality values in the above case (0-3) are:

0 - neither bit set, so the pixel is inside both the AST and the FLT mask (a source pixel).

1 - bit 1 set but not bit 2, so the pixel is outside the AST mask but inside the FLT mask (a border-line pixel).

2 - bit 2 set but not bit 1, so the pixel is inside the AST mask but outside the FLT mask (a border-line pixel).

3 - both bits set, so the pixel is inside neither mask (a background pixel).

In the following screen shot of a Quality array displayed using GAIA, the black areas have value zero and are thus inside both masks, the dark brown areas have value 2 and are thus inside the AST mask but outside the FLT mask. The light brown areas have value 3 and are inside neither mask. In this particular case, there are no  areas with a quality value of 1, so the FLT mask is contained entirely within the AST mask.

There are several commands within kappa that manipulate Quality arrays in various ways. For instance, the setbb command allows pixel data values to be set bad if the associated Quality value has a specified collection of set bits. Thus:

% setbb fred 1

will set all pixels bad except for those inside the AST mask. Likewise, 

% setbb fred 2

will set all pixels bad except for those inside the FLT mask.  Note, the change made by setbb is temporary - it can be undone by doing:

% setbb fred 0

To display the fred.sdf map and then overlay the AST mask in blue and the FLT mask in red, do:

% display fred mode=perc percentiles=\[2,98\]
% setbb fred 1
% contour fred clear=no mode=good labpos=! style='colour=blue'
% setbb fred 2
% contour fred clear=no mode=good labpos=! style='colour=red'
% setbb fred 0

Map Variances when using FLT.FILT_EDGE_LARGESCALE_LAST

The FLT.FILT_EDGE_LARGESCALE_LAST makemap configuration parameter allows you to create a map in which the background regions (i.e. the regions zeroed out by the AST mask) are filtered on a shorter scale than the source regions, thus producing a pleasingly flat background without removing structure from the source regions.

However, when using this option care needs to be taken when interpreting the variances stored with the map. Setting FLT.FILT_EDGE_LARGESCALE_LAST  causes the filter size to be changed on the last iteration. Whilst this does not affect the data values in the source regions much (since the source signal has been extracted into a separate model by the time the final iteration occurs), it does affect the variances. This is because the variances are formed from the time-stream residuals that are left after removal of the astronomical signal estimated on the penultimate iteration, and all these residuals - no matter where they are on the sky - are filtered using the smaller filter size given by   FLT.FILT_EDGE_LARGESCALE_LAST.

The upshot of all this is that the basic scheme described above results in artificially low variances in the source regions.  To avoid this, a change was introduced into SMURF on 13th February so that now, if the FLT.FILT_EDGE_LARGESCALE_LAST parameter is set, the variances stored in the map are the ones generated on the penultimate iteration (that is, the final iteration that used the basic filter size rather than the special last-iteration filter size).  This should give more realistic variances for the source regions, but at the cost of  unrealistically high variances in the background regions.