Open in Colab: https://colab.research.google.com/github/casangi/casadocs/blob/fe2b08a/docs/notebooks/single_dish_imaging.ipynb
Single-Dish Imaging
The data should already be T\(_{sys}\) and T\(_{sky}\) calibrated (at least into antenna temperature units, \(T_A^*\) [K]) , according to the process described in the single-dish calibration pages.
The CASA task used for imaging single-dish data is sdimaging. This task will return a four-dimensional array with two position axes, one frequency or velocity axis, and one polarization axis. The output of sdimaging is a CASA image that can be explored, analyzed, and manipulated in CASA, or exported into a versatile FITS image format via exportfits.
The sections below first describe the general process of gridding single-dish data followed by the actual procedures invoked with CASA.
Theoretical Description
A theoretical description of single-dish image generation and data gridding
For single-dish observations, ALMA uses on-the-fly mapping. The technique is described in Mangum et al. (2007) [1].
Converting single-dish observations into an image or cube is done almost entirely in the image domain. After taking and calibrating the data, the process follows three steps:
Forming the image grid
Populating the image grid
Smoothing the image data
The fundamental parameter relevant to image quality is the sampling interval. There are a number of sampling functions that need to be considered: the sky sampling function, the image grid sampling function, and the response function of the single-dish beam. These functions all convolve against each other to yield an effective image resolution somewhat poorer than the actual theoretical FWHM of the telescope primary beam.
The dimensions and extent of the image grid are determined by the mapped area on the sky. The gridding pixel size must be at least one half the size of the theoretical beam when convolved with the sky-sampling function. Since the sky sampling function is typically 1/3 to 1/5 of the primary beam, and the effective FWHM of the telescope and sky sampling function is close to that of the telescope anyway, it’s safe to use a pixel dimension that is 1/3th the width of the primary beam.
For example, a 30” telescope beam with a 6” sky sampling function has an effective FWHM of \(\sim \sqrt{(30^2+6^2)}\simeq\) 30.6”. Therefore, computing an image pixel size that is 30”/3 = 10”, is appropriately oversampling the effective beam FWHM and sampling interval.
After the coordinates of the data are transformed into sky coordinates, the image grid is formed with dimensions either consistent with the user specifications, or so that the image fully encompasses the observed sky positions.
For each pixel in the grid (e.g. in RA-Dec space), the gridding process searches through the data for measurements taken within some cutoff radius (specified by convsupport). Depending on their distance from the grid coordinate, the observation is weighted according to the kernel type and added together in the spatial domain (i.e. entire spectra are added together). If the clipminmax function is invoked, the maximum AND minimium data in the ensemble (prior to weighting) are rejected before summing. This process is repeated iteratively for each element in the grid.
Procedure for SD Imaging
CASA uses the sdimaging task to grid a single-dish image or cube. While the steps are described in detail here, examples of the full single-dish calibration and imaging processes can be found in the ALMA CASA Guides.
Default operation
The sdimaging task can determine and populate almost all the gridding parameters by default. Simply invoking sdimaging with the single-dish MeasurementSet and output file name will work. This will produce a single, potentially very large cube having as many channels as necessary to span the entire spectral range, with a spectral resolution equal to that of the observation.
sdimaging(infiles=sd_ms+'.bl',outfile=imagename)
The default parameter choices for imaging are selected as follows: the image pixel size is 1/3 the primary beam, the primary beam itself is computed based on the standard \(\frac{\lambda}{D}\), with some empirically-validated tapering applied. The image dimensions are determined by the spatial extent of the mapped area in the observation, and by default, all channels and all spectral windows are imaged, along with all antennas and all polarizations.
Customized operation
Of course, users can tune their data products by specifying the image size and dimensions, the frequency/velocity characteristics, the gridding and data filtering and smoothing parameters, and so on. The defaults for sdimaging for image resolution (i.e. cellsize in arcsec) is determined from the rest frequency of the 0th spectral window so that there are three pixel elements across the beam, the beam being calculated with \(b\times\frac{\lambda}{D}\). See information here about tapering: PrimaryBeamArcsec. The image extent is computed by default so that the sampled area is completely encompassed in a single rectangle, and the pixel dimension follows naturally from maxsize/cellsize. The default image center (the somewhat inappropriately-named phasecenter parameter) is computed simply as the center of that region.
These parameters can be left to be determined by sdimaging, or they can be determined using CASA tools.
Image dimensions and pixel interval
The image extent can be explicitly determined using aU.getTPSampling:
xSampling, ySampling, maxsize = aU.getTPSampling(refvis, showplot=True)
which returns an image output showing the scans, their sky angles, and positions in RA-Dec, as shown here:
Note that getTPSampling MUST operate on the original MeasurementSet (i.e. one that is not split or subselected). getTPSampling also yields the sampling rates in the x and y (i.e. azimuth and elevation) axes, as well as the maximum size of the image, in arcseconds.
The beam size used by sdimaging is determined using the aU.primaryBeamArcsec task, though this can also be invoked by the user and used to compute, for example, a cellsize and image size. The default for aU.primaryBeamArcsec corresponds to a 12m antenna with normal tapering. Setting the fwhmfactor modifies the beam taper (see discussion in PrimaryBeamArcsec).
freq=115.27e+9
fwhmfactor=1.13
diameter=12
theorybeam = aU.primaryBeamArcsec(frequency=freq*1e-9, fwhmfactor=fwhmfactor, diameter=diameter)
cell = theorybeam/9.0
imsize = int(round(maxsize/cell)*2)
The center of the image can be modified using the phasecenter parameter. Single-dish images actually have many phase centers, so the name is somewhat misleading. However it is preserved here for consistency with the interferometer terminology. In the context of single-dish data, phasecenter refers only to coordinates that will align with the center of the image, and this can be in J2000 or Az/El, e.g.
phasecenter='J2000 12h22m54.9 +15d49m15'
Frequency and/or velocity axis
The default rest frequency is the mean frequency of the first spectral window (i.e. that having the lowest spectral window ID). Of course it can instead be set by the user, or a different spectral window frequency can be used, extracted from the data using msmd tools:
msmd.open(vislist[0])
freq = msmd.meanfreq(spw)
msmd.close()
print "SPW %d: %.3f GHz" % (spw, freq*1e-9)
The third axis of the image cube can be specified using the veltype and outframe parameters. Many spectral-line observers will prefer to change these so the output has a velocity axis in the radio convention as follows:
veltype='radio',
outframe='lsrk',
and the rest frequency can be specified with:
restfreq='115.271204GHz'
The velocity extent of the image cube can be specified by selecting a spectral window (via the spw parameter), the channel range (via the nchan and start parameters), and the frequency/velocity resolution (via the width parameter). For example:
nchan=70,
mode='velocity',
start='1400km/s',
width='5km/s',
Gridding parameters
The gridding kernel defaults to a box shape, but it can be specified as a spherical (‘SF’), Primary beam (‘PB’), Gaussian (‘GAUSS’) or Gaussian*Jinc (GJINC). The recommended setting for ALMA data is a spherical (‘SF’) kernel. The convsupport parameter defines the cut-off radius for ‘SF’ in units of pixels, defaulting to 3 pixels. However, the recommended value for ALMA data is convsupport=6 (see sdimaging and Mangum et al. 2007 [1] for more details on these parameters).
The parameter stokes specifies the stokes product. At present, the weighting for stokes I is computed consistently with historical usage: I=XX/2+YY/2. While this is mathematically consistent with the computation of stokes I, it is an incorrect treatment since the computation necessarily must incorporate the contributions from Q and U. Ordinarily, these terms cancel out from the computation of stokes I, but their error parameters must be incorporated, and historically, this is not respected.
CASA development is seeking to make the computation of the weights consistent with a proper computation of stokes I, and this is done in sdfit, but it is not yet completed for sdimaging. However, to emphasize, while the current implementation of computation for stokes I by sdimaging is consistent with convention, the convention is formally incorrect.
Example script
Fully specified, a call to sdimaging might look like the following:
sdimaging(infiles=sd_ms+'.bl',
field='M42',
spw='%s'%(spw),
nchan=70,
mode='velocity',
start='1400km/s',
width='5km/s',
veltype='radio',
outframe='lsrk',
restfreq='%sGHz'%(freq/1e+9),
gridfunction='SF',
convsupport=6,
stokes='I',
phasecenter='J2000 12h22m54.9 +15d49m15',
ephemsrcname='',
imsize=imsize,
cell='%sarcsec'%(cell),
overwrite=True,
outfile=imagename)
The products here are the image data, returned in the variable ‘imagename’, and also a map of weights: <imagename>.weight. The weights indicate the robustness of the gridded data on a per-pixel basis, and are important when performing further computations and analysis with the image products.
Bibliography
Mangum, et al. 2007, A&A, 474, 679-687 (http://www.aanda.org/articles/aa/pdf/2007/41/aa7811-07.pdf)