Content-type: text/html Manpage of mosaic


Section: User Commands (1)
Index Return to Main Contents


mosaic - Deconvolution using mosaicing method  






Mosaic does a simultaneous deconvolution of a set of adjacent pointings, combining them into a single map. This should yield a better result than individually deconvolving the pointings. It replaces the dirty point spread function (i.e. one with many and/or deep sidelobes) with one that is well-behaved (i.e. a gaussian). [In the near future, it will be possible to include single-dish data.]

A short outline of the algorithm: start with an initial model. For each pointing, multiply the model with the primary beam pattern; convolve with the synthesized beam of that pointing. Subtract the observation and calculate the chi^2 of the residual. Then apply the formulae of the Maximum Entropy Method, in which unobserved data are constrained by a "reference map". To find a step, use Lagrange Multipliers and a Newton-Raphson iterative method. Drive the sum of the chi^2 values of all maps to a value such that on average the residual deviates by one sigma from the prediction and such that the total flux converges toward a wanted value.

The many keywords below come in five groups:

  - in, beam, pbpar, maxcorr, rms -- specify interferometer inputs
  - region#, planes, initial, reference, flux -- set up constraints
  - out, mode, beamsize, save, saveiter -- specify outputs
  - maxiters, tol, measure, lm0, meslev -- control parameters
Obligatory keywords are: in, beam, rms and out; the rest is for fine tuning.

Inputs for interferometer data are: a list of miriad datasets with input observations, a list of the corresponding synthesized beams, an optional list of primary beam descriptions and a list of values for the rms noise in each input map.

For each of the input maps it is possible to select a region outside of which the model is set to zero (see the region and reference keywords). A list of common planes can be given to simplify the selection of channels. If datasets are not on the same frequency grid, use miriad task regrid to make them fit.

The output is a miriad dataset with the model. It is also possible (and advisable) to write output datasets containing the residuals of the input map, a combined residual (without primary beam corrections), the model convolved with a "clean beam", and a final map (the sum of the primary beam corrected residual and the convolved model). The reference coordinate (crval) of the output maps is that of the first in the list of input maps. Note: the mask of the model (outside which it is zero) is only written to the mask item of the 'convolved' output map.

Calculating the secondary outputs (residuals, convolved and final map) is equivalent to using 'restor' after a 'clean' operation. They can also be calculated from a read-in model map (see the mode keyword).

Most intermediate results (like the evolving model and residuals) can be written by using the save= keyword.

The program has three other independent uses: 1) to write a linear mosaic of the primary beam corrected, variance-weighted input maps; 2) to write a noise map, i.e. the harmonic mean of the variances; 3) to write primary-beam-corrected input maps. A choice is made using the mode= keyword.

An example of standard simple usage is:

 mosaic in=map* beam=beam* rms=1,... flux=-10 out=out mode=m,f,r,c


Beware: the value given for the rms is of paramount importance in determining convergence. In fact, it is the parameter that has the most influence on the reliability of the result. If it is too low or too high convergence may not occur.

A current limitation is that the crvals of the input maps must be an integer number of pixels offset from each other. The program checks this and tells you what input to use for the invert keyword offset= to get it right. The reason for this limitation is that mosaic does not regrid the maps.

For large fields (larger than about one degree) geometrical projection effects may become important. Maps made with invert will not be on the same tangential plane. Usually this can be ignored. Strictly speaking, however, they should be reprojected (using the miriad task regrid). Unfortunately, miriad does not differentiate between pointing center and projection center, and this causes trouble. So, after using regrid, you need to fix the crval and crpix of the reprojected datasets so that they indicate the pointing center of the observation.

As the model is forced to be non-zero everywhere inside the selected region, it can also be non-zero where there is no actual emission. Thus, it may be even more important than in CLEAN to define a region where the signal is expected. Use the region#= or the reference= keyword (see there for details). One way to create a region is to first run mosaic with mode=m,c, without region selection; then take the resulting convolved model, and fiddle a bit with it; next use it as input to the reference keyword, making use of the clipping option. A choice for the region that is too small will result in missed structure in the model. Or, the bias flux that is spread around in the map may become too high. A choice for the region that is too large can result in "structure" being found in the noise.

It is furthermore important to put in a good value for the total flux (if only as initial estimate). If you don't have this, don't be surprised if the answer looks bad or convergence does not occur. You might try to iterate by changing the initial input value until it does not vary much from the value at the last iteration.

A good reference map improves convergence. Without it, in principle the final result should be the same, but you may have to wait longer if the source structure is complex. The reference map gives the algorithm reasonable data on spatial frequencies that were not observed. There are several possible ways to create a reference map (or an initial estimate):

  a) a flat map (default, usually good enough and preferable)
  b) the sum of the primary beam patterns (for sources whose diameter
     is comparable to the primary beam) (often unreliable)
  c) a single-dish map gridded to the interferometer gridspacing
  d) a map of the same source at a different frequency
  e) the result of a previous or lower-resolution observation
The latter three cases require you to specify the reference map as input. If any a-priori information is present, use this to create an initial estimate and reference map.

Officially, miriad has no header variable that represents the pointing center. The uv-variables "obsra" and "obsdec" come closest. They are defined as the phase center at the epoch of the observation. To get the phase center at the epoch of the map requires a precessional correction, which is not done (yet?). Thus, "crval1, crval2" are used instead as the coordinates of the pointing center. For miriad datasets this should almost alway be OK, but it is a non-guaranteed use of the variables. For maps imported from aips, check whether crval/crpix indeed give the pointing center. If "crval1, crval2" does not give the pointing center, the primary beam is offset and the algorithm will fail.

Be very cautious using this method. If you don't understand it, it is easy to obtain wrong output maps. Don't blame the program if you fail at the first try, fiddle with the parameters or ask for help (remember the need for an accurate rms).  


List of input maps with observations.
List of input beam datasets (one for each map). The lengths of the x and y axes must be a power of 2.
List of values giving the formula used to describe the primary beam. There should be a value for each input map, unless "..." is used. [The default is "pbpar=tel,...", i.e. use the header for all input maps, as described below.] Values can be:
  - ""         => tel
  - none       => no primary beam correction
  - tel        => use formula for telescope as described below
                  if telescope not known, make a gaussian with fwhm
                  from header item "pbfwhm"
  - gauss,fwhm => create a gaussian primary beam (cut off at 0.05)
  - airy,fwhm  => create an airy disk
  - airyc,fwhm => create an airy disk, but only out to first zero
  - ...        => repeat the last given value for remaining inputs
  If fwhm > 0  => make primary beam of given type with this fwhm
                 (if restfreq>10 GHz units are arcsec, else arcmin)
  If fwhm = 0  => read fwhm from header item "pbfwhm"; if not
                  present, get fwhm from "restfreq" and "telescop"
If "pbpar=tel", header items "telescop" and "restfreq" are used to make the primary beam: a polynomial for the WSRT, VLA and ATCA, a gaussian for FST and BIMA; these functions are cut off at a level of about 0.05 (see keyword maxcorr). If "restfreq" is not present, but the type of the "freq" axis is "FREQ", then crval, crpix and cdelt are used to get the frequency of plane 1. Example 1: with two input maps, pbpar=gauss,120,gauss,120 makes gaussian primary beams with fwhm 120 arcsec. Example 2: pbpar=airy,120,tel,gauss,0,gauss,118 with four input maps makes an airy disk with fwhm 120 arcsec for the first, uses the telescope formula for the second, makes a gaussian with fwhm found from "restfreq" and "telescop" for the third and makes a gaussian with fwhm 118 arcsec for the fourth.
Influences to how far out the primary beam correction is done: everywhere where the correction factor is less than maxcorr. Not valid for airy shapes. Mainly useful with mode=pbc. [Default is telescope dependent (20 for gauss, 16 for WSRT, 43 for VLA, 33 for ATCA).]
List of values of the pixel-to-pixel rms (in Jy/beam, as given by imstat), one for each map. Determine this from signal-free areas or channels. You can repeat the value of the rms by using "...". E.g.: "rms=0.1,..." will make the rms 0.1 for each input map. Entering "rms=0.1" with more than one input map will result in an error message, i.e., the value is not automatically repeated; this is a safety feature. Beware: getting the rms right is the most important thing to get proper convergence.
See region#.
See region#.
See region#.
region# stands for region1, region2, region3, etc. You never actually type the '#'. There can be one region keyword for each input map. For the first map use 'region1', for the second 'region2', etc. For the format of the specification, see 'doc region' or the users manual. If you are using the miriad 'shell', region can only be given for the first three input maps. To use more, use a c-shell command, making sure that all '(' and ')' are put within single quotations marks. If none of the region# keywords is given, the region is determined either by the reference keyword (if a cutoff is given), or by merging the inner quarters of all maps. For each individual region keyword, the list of boxes and polygons is farther out than the primary beam cutoff are masked out. If the input file has a mask item or if the keyword specification contains 'mask', these masks are 'anded' in. The set of input regions, one for each map, is 'anded' together to create a region in the mosaiced map. Everything outside this region will be set to zero in the final model. Example1: region1=quart(1) makes the mosaic model non-zero only in the inner quarter of the first map. Example2: region1=relpix,box(-30,-30,30,30),box(-50,-50,0,0) combines the two boxes into a single, larger region, and makes the mosaic model non-zero in this region. Example3: region1=quart(1) region2=quart(1) makes the model non-zero in the overlap region between the inner quarters of the first and second input map. It is easy to check the region that is produced, by using the save=ref keyword. This write out the reference map, which is non-zero only inside the selected region.
To make it a little easier to specify a list of input planes, use this keyword to work on same planes for all input datasets. It has the same syntax as the standard region keyword. It supersedes any image specifications in individual region keywords. Subcommands other than