The Dynamic Iterative Map-Maker, hereafter just referred to as the map-maker is the tool you will use to produce SCUBA-2 maps, and is implemented by the Smurf makemap command. It performs all pre-processing steps to clean the data, followed by solving for multiple signal components using an iterative algorithm, and binning the resulting time-series data to produce a final science map.
The makemap command can be invoked in two ways: 1) directly by typing “
makemap” in response to a unix
shell prompt (see Chapter 5), or 2) indirectly as part of an ORAC-DR recipe (see Chapter 4).
This chapter describes how the map-maker produces a science image from raw SCUBA-2 data. It should be considered essential reading as it provides an understanding of how your reduced image was produced. This is particularly true if you wish to modify the default map-maker parameters. If you prefer to jump straight in to the data reduction go to Chapter 4.
The time-varying signal recorded by a bolometer, , contains contributions from a variety of sources:
where is the atmospheric extinction, is the astronomical signal (time varying as the telescope scans across the sky), is the uncorrelated white noise, is an optional scalefactor which is unique to each bolometer, is the correlated (or, common mode) signal, and is the (predominantly low-frequency) noise which does not have a simple correlation relationship that would be included in the noise term.
The map-maker works by producing individual models of the various components that make up the signal recorded by each bolometer. It models and removes the extinction and noise components in order of decreasing magnitude, ultimately leaving just the astronomical signal plus residual noise. The modelled components are listed in Table 3.1 and described more fully in Section 3.3.
A configuration file should usually be supplied when running makemap directly. This file holds a set of parameter values that control all aspects of the map-maker, including details of the pre-processing steps, which model components to include, parameters that control the determination of each model, and the stopping (or convergence) criteria. If no configuration file is supplied when running makemap1, a set of default parameter values will be used. There are a great number of these parameters, but fortunately not all of them are of interest to the typical user. Appendix H documents the parameters you are more likely to find of interest, including the default value that is assigned to each parameter if you do not give it a value in your configuration file. For a full list of all available parameters, see the appendix “Configuration Parameters” within SUN/258.
It is possible to create a map without supplying a configuration file to makemap (i.e. leaving all configuration
parameters set to their default values—“
makemap) but it is not recommeneded
since it will not give optimal results for your particular observation. For this reason, specialised configuration
files have been developed which are tailored to different science goals, be they detecting faint galaxies or
mapping large molecular clouds. A description of these specialised configuration files can be found in
Note, when using the pipeline to create a map, rather than running makemap directly, the pipeline will always use a configuration file—one of the standard configuration files will be selected if none is specified by you.
This section describes the basic map-making process used by the default configuration. It may be modified in many ways by supplying a configuration file containing alternative parameter values.
Figure 3.1 shows the flow chart of the basic map-making process. It is divided into two sections: the pre-processing stage where the data are cleaned, then the iterative stage where the different models are subtracted, a sky map created and the convergence checked.
These time-series are then re-sampled at a rate that matches the requested pixel size—the equivalent to applying a low-pass filter. Down-sampling saves time and memory usage when running the map-maker, without losing any significant information. The degree of down-sampling applied depends on how fast the telescope was moving during the observation, the requested pixel size. In general, slower scans will be more heavily down-sampled than faster scans, and fast scans may not be down-sampled at all.
A number of cleaning steps are then run: bolometers that have unusually high noise levels are identified and excluded from further processing, any sudden steps in the base-line of each bolometer are identified and corrected, and a polynomial estimate of the base-line is removed from each bolometer.
Next comes the iterative stage. In each iteration, estimates are produced for each model
component. These are removed from the cleaned time-series data and the remaining data values
are binned into a map. However, in general, the presence of astronomical signal within the original
time series will have upset the estimates of the
FLT models, causing the final map
to be inaccurate. But now that we have a map (albeit an inaccurate map), we can sample the
map at the position of each bolometer value to get an estimate of the astronomical component
within the original time-series data. We then remove this astronomical signal from the original
time-series data and re-estimate the models. These new model estimates should be more accurate
since they are not so heavily influenced by the astronomical signal. In turn, this allows us to create
a more accurate map. We repeat this process until the map does not change significantly with
further iterations. Whilst not mathematically rigorous, we expect the process to converge because
the astronomical signal is in general much smaller than any of the other components and so each
iteration introduces a very small fractional change in the model estimates.
Within each iteration, the first components to be modelled and removed are
work together to calculate and remove the average signal template of all bolometers, allowing each
bolometer to have an arbitrary gain and offset. In addition, the
GAI model will flag as unusable
any bolometer in which the signal looks very different to the common-mode.
The next model (
EXT) applies a multiplicative extinction correction to the data. Following this,
FLT model filters each bolometer time-stream independently to remove low frequencies that
correspond to angular scales larger than 600 arc-sec and 300 arc-sec at 450m
After these models have been removed, the AST model (i.e. the estimate of the astronomical signal)
produced by the previous iteration, is added back onto the remaining time-series data2,
and the new values are binned up on the sky to produce a new estimate of the final science map.
Since many samples typically contribute to the estimate of the signal in a given pixel, the noise
is greatly reduced compared with the time-series data. The variance of each map pixel value is
determined from the spread of sample values that contribute to the pixel. The value of this map
at the position of each bolometer sample is then found. This forms the new
AST model which is
removed from the time-streams, leaving just the residual noise.
NOI model then measures the noise in the residual signals for each bolometer to establish
weights for the data as they are placed into the map in subsequent iterations. This is only done on
the first iteration. Subsequent iterations re-use the weights established on the first iteration.
Convergence is checked against the parameters detailed in the configuration file. Convergence
is achieved either when the requested number of iterations has been completed, or when the
mean change in the map pixel values is less than a specified fraction of a standard deviation (see
maptol). If more iterations are required or (in the latter case) the map is still changing
significantly, all the model—except for AST—are added back onto the residuals, thus reconstructing
the original time-series but without the astronomical signal. This is the signal upon which the
model estimates will be based on the next iteration.
For full details of the map-maker see Chapin et al. (2013) .
| ||Common-mode signal|
| ||Gains that scale each bolometer to the common-mode|
| ||Extinction correction|
| ||Filter that removes low frequencies|
| ||Astronomical signal|
| ||Residual noise|
The more data the map-maker processes at once, the better chance it has of determining the difference between sky signal and background noise.
Information about chunking can be obtained in several ways:
The final number indicates the number of chunks being used (one in this case).
indicates that the map in file
mymap.sdf did not suffer from chunking as a result of low memory.
indicates that the time-series data from which
mymap.sdf was created was broken into two contiguous
groups of samples. This indicates a potential problem with the data—the NCONTIG header should
normally be one.
In daisy mode chunking is less of a concern as the entire map area is covered many times in the space of a single observation.
For pong maps chunking is a bigger concern, with the maximum number of chunks that can be tolerated dependent on the number of map rotations. For example, a 40-minute pong map with eight rotations may get divided into three or four chunks. Although not ideal, this will mean that each point is still covered by two or three passes. Fewer passes than this however and the map-maker become less effective.
Be aware that focus observations are intentionally dis-contiguous. If you attempt to make a map from a focus observation, it will be processed in several chunks no matter how much memory your computer has.
The list of models to be evaluated (and removed) by the map-maker, and the order in which they are used
during the iterative stage, is given by the
modelorder parameter in the configuration file. These
models are modular however, so their order may be changed. The default model order is given
The only configuration file not following this model order is the one tailored for blank field maps which does
not include a
FLT model in the iterative stage but instead includes an equivalent high-pass filter as part of the
cleaning stage (see Section 3.7).
Below is an introduction to each model. More complete descriptions of the models and all the associated caveats
can be found in Chapin et al. (2013) . Details of each model are controlled by a set of configuration parameters
that begin with the name of the model. For instance, all parameters related to the
COM model will be of the form
Extended emission on a scale larger than the array footprint on the sky will contribute a signal indistinguishable from a common-mode signal. This puts an upper limit on the spatial scale of astronomical emission that can be recovered.
Examples of the time traces for a single bolometer from these models is shown in Figure 3.2. These traces cover
a subset of an observation of the secondary calibrator CRL 2688. You can clearly see the dominance of the
model which is removed first. The
FLT model stores the data removed by the high-pass filter. In the
CRL 2688 is clearly seen as positive spikes which appear when the bolometer passes over the source. Finally,
the residual signal stored in
RES is flat, indicating that most of the signal has been successfully accounted for by
the other model components.
The map-maker will stop processing either when the requested number of iterations has been completed OR when the convergence criterion specified in the configuration file is reached.
Option 1: Fixed number of iterations
Specifying the number of iterations in the configuration file is done via the
numiter parameter. This is shown
dimmconfig_jsa_generic.lis (JSA stands for “JCMT Science Archive”).
A positive value for
numiter means that the requested number of iterations will be executed. A negative
value, as in the example above, indicates that no more than this number of iterations should be
performed, but that it may stop at fewer if convergence (according to the noise criterion below) has been
Option 2: Convergence parameter
The convergence criterion is set by the
When running the map-maker,
maptol gives the average normalised change in the value of
map pixels between subsequent iterations. Convergence is reached when the mean change
across all pixels is less than the parameter
maptol. It has units of the noise in the map, thus
maptol0.05 means a mean change in pixel value of
This option has the advantage of directly assessing the noise in the resulting map.
The map-maker displays the normalised change in pixel values between iterations at the end of
each iteration. This value typically drops rapidly to begin with and then flattens out, decreasing
increasingly slowly. The printed values can be inspected to check convergence is proceeding as
expected. Be aware that setting
maptol to much lower than 0.05 will dramatically increase the
length of time needed to produce the final map, and may possibly result in convergence never being
If you wish to perform further checks on the progress of convergence while running the map-maker, you can
itermap option. This causes the map created by each iteration to be dumped to disk for inspection (see
The term masking refers to the inclusion or exclusion of selected bolometer samples from the estimate of specific
models, based on the positions of those samples within the output map. Masking can be applied
independently to the
COM models, as described in the later sub-sections within this
A mask is a selection of pixels within the output map. Normally these pixels form one or more contiguous regions that enclose the areas where significant sub-mm emission is expected (the “source” regions). Pixels outside this mask are referred to as “background pixels”. The use to which such masks are put differs from model to model as described later in model-specific sub-sections.
Masks fall into two broad categories:
Static masks are appropriate in cases where the areas containing significant sub-mm emission are
known in advance. If this is not the case, then the mask needs to be based on evidence within the
maps created at the end of each iteration, resulting in the masks evolving from iteration to
iteration as the maps converge towards the final solution. Such dynamic masks have a potential
problem in that they can upset the convergence process by introducing extra flexibility into the
algorithm. For this reason, an option is provided to freeze these masks after a fixed number of
ast.zero_freeze) and this is used within the specialist configuration file
in Section 3.8.
Each of the three models—
COM —have an equivalent set of parameters to specify the masking (if
any) that should be used. Each set starts with the model name, followed by a dot, followed by the parameter
name. To enable masking, one or more of these parameters needs to be set to a suitable value to define a mask,
as described below.
Static masks can be defined either as circles of specified radius, centred on the source position
(suitable for compact sources—see for instance
ast.zero_circle), or by an external image that uses
special pixel values to indicate the source regions and the background regions (suitable for extended
sources—see for instance
ast.zero_mask). Within such “external masks” all background pixels
should be set to the Starlink bad value (see “Bad-pixel Masking” in SUN/95). Pixels with any other
value are treated as source pixels. Section 6.6 describes some ways in which external masks can be
Dynamic masks can be defined in several different ways:
ast.zero_snrlo). This allows a low S/N threshold to be used (which can sometimes help to avoid negative bowling around the edges of extended sources) without introducing lots of isolated pixels into the mask because of noise spikes in the map.
ast.zero_lowhits). Because of the scanning strategy used with SCUBA-2 (see Section 2.2), pixels near the edge of the map will receive fewer bolometer samples than those near the centre, and will therefore be noisier. Thus this scheme allows the noisy edges of the map to be masked out. Whilst this is technically a dynamic mask—because the number of samples falling in each pixel could potentially change from iteration to iteration—in fact a mask based on hits per pixel will usually only change very slightly from iteration to iteration.
ast.zero_snr_ffclean). This can be a useful alternative to simple S/N thresholding in cases where bright artificial large-scale structures appear in the map. Such structures will result in artificially large source regions in a simple S/N mask, which may in turn cause the artificial structures to become even brighter and larger in subsequent iterations. Using an ffclean-based mask instead will restrict the mask to the brighter astronomical cores on top of such artificial large-scale structures.
Multiple masks can be used by selecting the appropriate configuration parameters as described above. Whether
the combined mask represents the union or intersection of the individual masks is selected using
com.zero_union. By default, the union is
The masks used on the last iteration are stored in the
Quality component of the NDF containing the output
map. A separate bit of the quality component is used for each model (
COM) that is being masked.
These bits can be viewed in various ways—see Section 8.6 and Section 5.3.2.
The iterative map-maker algorithm described earlier in this chapter divides each cleaned bolometer value into four additive parts:
Whilst the value of each model is constrained to some extent by the manner in which the model is calculated,
there remains a large degree of degeneracy in the algorithm. In other words, there are many different ways in
which each bolometer value can be divided up into the four values listed above. For instance, any arbitrary
change in the
AST model can be accommodated provided that equal and opposite changes are
introduced into the other model values so that the total value of all four components remains the
The degeneracy between
AST is particularly problematic, since it allows large scale features
to appear in the map (on the scale of the sub-array footprint) that are corrected for by equal and
opposite features in the
COM model. Such features can appear as gradients or “saddles” across the final
The purpose of masking the
AST model is to apply an extra constraint to the system to help break this
degeneracy, and so avoid these large-scale features. It functions by forcing the
AST model to zero at the end of
each iteration for all samples that fall outside the masked area on the sky (i.e. samples that fall in the
“background” areas of the map). Forcing the
AST model to zero in this way means that the
COM model is more
Whilst this constraint is applied only in the background regions, its effects are also felt part way into the source
region because of the way in which the
COM model averages over a sub-array footprint. Thus source regions that
are much larger than a sub-array will still suffer from the degeneracy between
AST, allowing arbitrary
large-scale background features to develop within the region, but these features will become weaker as the
edges of the mask are approached.
In order to avoid the final map containing zero values in the background region, it is necessary to avoid using
AST-masking on the final iteration. For this reason, by default one further iteration is performed—without
AST-masking—once the basic iterative algorithm has converged (see
All the configuration files supplied with Smurf use
AST masking, with the exception of
FLT model estimates and removes a low frequency background from each bolometer time-stream. To do this
it uses an FFT-based high-pass filter. Ringing around sudden transients—such as astronomical sources — is a
natural consequence of all such filters, resulting in alternating bright and dark rings around sources in the map.
In most cases, the iterative nature of the map-making algorithm eventually results in the rings reducing to a
much lower level, as more of the astronomical signal is removed from the time stream on each iteration, causing
the ring-inducing transients to become weaker. However, this is not always completely effective, and also slows
For these reasons, it is often helpful to mask the
FLT model. This causes each bolometer time stream to be
blanked out as the bolometer passes over the source regions in the mask. Such blanked sections are replaced by
linear interpolation between the data on either side of the blanked section. The effect of this is to remove the
transients that cause the ringing in the
FLT model, and thus get a better estimate of the low frequency
background in each bolometer.
Such masking is most useful at the start of the iterative process, before the bulk of the astronomical signal has
been removed from the time-streams. Indeed, if it is allowed to continue for too many iterations it can cause its
own problems that tend to slow down convergence. For this reason, if
FLT masking is used, it is (by default)
limited to two iterations—see
A difficulty with
FLT masking is that the
FLT model is evaluated before the map is made on each iteration. This
means that when the
FLT model is evaluated on the first iteration, there is no map available. This is not a
problem if a static mask is used, but if a dynamic mask is used, then it is not possible to calculate the
on the first iteration (subsequent iterations will create the mask on the basis of the map created at the end of the
previous iteration). This is a shame since
FLT masking is most beneficial on the first few iterations, whilst the
bulk of the astronomical signal is still present in the time streams. A scheme to get around this problem is
described in Section 3.6.
All the configuration files supplied with Smurf use
FLT masking, with the exception of
The value of the
COM model at a specific time-slice is the mean of all remaining bolometer values at that time
slice. If the
COM model is masked, bolometer samples that fall within the source regions defined by the mask are
excluded from this mean, thus removing the bias introduced into the mean by bright astronomical sources. This
results in the
COM model more accurately representing the varying common-mode background value at each
time slice. This can help to reduce negative bowling around sources in the map, and can help to
reduce the amount of data that is rejected as being too poorly correlated with the common-mode
COM masking has the same difficulty as
FLT masking in that the
COM model is evaluated before the map is made on
each iteration, and so it is difficult to use dynamic masks with
COM masking. However, the solution described in
Section 3.6 can be used for
COM masking as for
Currently, none of the configuration files supplied with Smurf use
As discussed in Section 3.5,
FLT masking is most effective on the early iterations (particularly the first iteration).
This is because the bulk of the astronomical signal is still present in the time streams and is likely to cause bad
FLT masking is not used. However, since the
FLT model is evaluated before the map is made
within each iteration, there is no map available to use as the basis for a
FLT mask on the very first
One way to get round this is simply to use a static mask that does not require a map to be available. However, this may not be possible, particularly for extended sources where no a-priori knowledge of the source areas in the sub-mm may be available.
Another scheme is provided by the
ast.skip configuration parameter. If
ast.skip is set to a positive integer
, the subtraction
AST model from the time-series that normally occurs at the end of each iteration is not performed on the first
iterations. This means
that for the first
iterations, the time-series data remain unchanged, and consequently the models and the resulting
map are unchanged. Of itself, this does nothing useful—it merely delays the first useful iteration.
However, when combined with
FLT masking, it can be very beneficial because it allows a mask to
be created before the first “real” iteration occurs (i.e. the first iteration that subtracts off the
model). This means that the first “real” map will have far less ringing around bright sources than
would otherwise be the case. This in turn means that far fewer iterations are required to correct this
It is also possible to use
COM masking, but the benefits are not so noticeable as for
The default value for
ast.skip is zero, but the
dimmconfig_jsa_generic configuration files supplied with Smurf both set
= 5 (and also use
Maps made by the pipeline using the REDUCE_SCAN recipe will use
some alternative configuration has been specified. This configuration file is intended to provide a reasonably
good map for all types of observations. However, compromises have been made to reach that balance, the main
one being that some real extended structure is sacrificed in order to avoid introducing artificial large scale
dimmconfig_jsa_generic.lis is always a good place to start, you will want to follow this up with a
specialised configuration that will suit your observation. A few specialised configuration files are supplied with
Smurf and can be found in
Note, any parameter for which no value is specified within the supplied configuration will take on the default value listed in Appendix H.
Below is a description of each of the specialised configuration files. The tables following each description list the
parameters set by each file. The files themselves are stored in
$STARLINK_DIR/share/smurf/—they are simple
text files and so can be inspected easily.
dimmconfig_jsa_generic configuration file is designed to handle all science data. This is a good
configuration file to start with if you are unsure how to proceed with your data reduction. You
may notice that elements of this configuration are common to the other specialised configuration
This configuration is optimised to suppress the creation of artificial large scale structure within the final map. Consequently, some real extended emission may be lost if this configuration is used. This goal is achieved by:
flt.filt_edge_largescale) to 200 arc-seconds so that the
FLTmodel uses a heavier than normal filter.
com.perarrayto 1, indicating that four separate
COMmodels should be used, one for each sub-array4.
ASTmodel to zero in the background regions after each iteration (see parameters
ast.zero_snrlo—also see Section 3.5.1).
In addition the S/N threshold for DC steps (
dcthresh) is relaxed from 25 in the default file to 100 to avoid
problems bright sources triggering the step detection algorithm.
The other major feature of this configuration is that the subtraction of the
AST model is skipped on the first five
iterations. In conjunction with the use of
FLT masking, this helps to speed up convergence and minimise dark
rings around bright source (see Section 3.6.
dimmconfig_blank_field configuration is tuned for blank field surveys for which the goal is to detect
extremely low signal-to-noise point sources within otherwise blank fields.
Applying a high-pass filter (
FLT) on each iteration can result in convergence problems when there is little or no signal in
the map. Instead, a single, harsher high-pass filter is applied as a pre-processing step (corresponding to 200-arc-sec scales
at both 450 m
and 850 m).
There are also more conservative cuts to remove noisy/problematic bolometers. Only 4 (positive) iterations are
requested as there is no signal to confuse to models.
com.perarray = 1 requires the
COM model to be fit to each sub-array independently. This improves
the overall fit but with the loss of any structure on scales larger than a single sub-array—not an issue for blank
Figure 3.3 shows the sharp contrast in the output map between reducing data with the default configuration
file and using
Blank-field maps commonly have a matched filter applied as part of the post-processing to aid source detection (see Section 8.7), however this is not applied by the map-maker.
dimmconfig_bright_compact configuration should be used for producing maps of bright, isolated
compact sources that are positioned at the centre of the map (e.g. calibrators).
The addition of the
ast.zero_circle parameter is used to constrain the map to zero beyond a radius of
1 arc-min from the source centre (see Section 3.5.1). This strategy helps with map convergence significantly, and
can provide good maps of bright sources, even in cases where scan patterns failed to complete in
com.perarray is set to 1 indicating that a
COM model should be fit separately for each sub-array. This is not
advised for extended sources as signal on scales larger than a single sub-array is lost, but is fine
for a compact central source. Likewise, the filtering is tighter. The S/N threshold for DC steps
dcthresh) is relaxed from 25 in the default file to 100 to avoid problems associated with bright
The values of the parameters controlling the
COM model are modified to relax the test that compares each
bolometer to the common mode. This is because the algorithm that compares each bolometer to the
common-mode can be fooled by a bright compact source passing across individual bolometers. Likewise the test
for DC steps within each bolometer baseline is also relaxed (parameter
dimmconfig_bright_extended configuration is for reducing maps that contain emission that is extended
to some degree and contains some bright regions. Here, the combination of
ast.zero_snrlo is used to identify components of the
AST model wherever the S/N is higher than
3 and then to expand this basic
mask down to an SNR equal to 2
without introducing any new isolated source areas (see Section 3.5.1).
numiter has been raised to
as more iterations are required to maximise the sensitivity to large dynamic signal ranges in the
The other major feature of this configuration is that the subtraction of the
AST model is skipped on the first five
iterations. In conjunction with the use of
FLT masking, this helps to speed up convergence and minimise dark
bowls around bright source (see Section 3.6.
Figure 3.4 shows a comparison between maps reduced with the default configuration file and using
dimmconfig_bright_extended.lis; the most noticeable difference is the improvement in the bowling around
Two additional configuration files are available that can be used to supplement your configuration file in order to solve specific problems:
dimmconfig_fix_convergenceto help your map converge. It does this by preventing masks from changing, and by freezing the flags that identify aberrant bolometers. Both of these restrictions are imposed only after 10 iterations have been performed.
dimmconfig_fix_blobsto remove large “blooms” or smooth blobs of spurious emission in your map. These can be caused by ringing induced by unidentified steps or spikes in the bolometer time-streams. The fix is to look for oscillations within the time-streams following the
FLTmodel, and flag such time-slices as unusable. In addition, time slices at which the four sub-arrays see significantly different common-mode values are also flagged as unusable.
These files provide a few additional settings that should be added to some other configuration in order to
achieve the desired fix. For instance, a configuration file containing the following two lines would indicate that
values from both
dimmconfig_fix_blobs.lis should be
This reads a set of configuration parameter from the
dimmconfig_bright_extended.lis file, and then assigns
new values for just those parameters defined in
dimmconfig_fix_blobs.lis. These values override the values
The parameters in these files are discussed further Chapter 6.
2This step is omitted on the first iteration since no estimate of AST has yet been created.
3Although by default masks are never frozen.
4By default, a single
COM model is used that represents the mean data value across all four sub-arrays.