docs/notebooks/single_dish_imaging.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "single_dish_imaging.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YP24b2iukunS"
      },
      "source": [
        "# Single-Dish Imaging \n",
        "\n",
        "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.\n",
        "\n",
        "The CASA task used for imaging single-dish data is [sdimaging](../api/casatasks.rst#single-dish).  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](../api/casatasks.rst#input-output).\n",
        "\n",
        "The sections below first describe the general process of gridding single-dish data followed by the actual procedures invoked with CASA.\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CLcdEBd6kunT"
      },
      "source": [
        "## Theoretical Description \n",
        "\n",
        "A theoretical description of single-dish image generation and data gridding\n",
        "\n",
        "For single-dish observations, ALMA uses on-the-fly mapping.  The technique is described in Mangum et al. (2007) [\\[1\\]](#Bibliography).\n",
        "\n",
        "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:\n",
        "\n",
        "1.  Forming the image grid\n",
        "2.  Populating the image grid\n",
        "3.  Smoothing the image data\n",
        "\n",
        "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.\n",
        "\n",
        "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.\n",
        "\n",
        "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.\n",
        "\n",
        "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.  \n",
        "\n",
        "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.\n",
        "\n",
        " \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cXQqTAybkunU"
      },
      "source": [
        "## Procedure for SD Imaging \n",
        "\n",
        "CASA uses the [sdimaging](../api/casatasks.rst#single-dish) task to grid a single-dish image or cube. While the steps are described in detail here, an example of the full single-dish calibration and imaging processes can be found in the [M100 Band 3 Single Dish CASAguide](https://casaguides.nrao.edu/index.php/M100_Band3_SingleDish_4.3).\n",
        "\n",
        "**Default operation**\n",
        "\n",
        "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.\n",
        "\n",
        "```\n",
        "sdimaging(infiles=sd_ms+'.bl',outfile=imagename)\n",
        "```\n",
        "\n",
        "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.\n",
        "\n",
        "**Customized operation**\n",
        "\n",
        "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](https://safe.nrao.edu/wiki/bin/view/ALMA/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.\n",
        "\n",
        "These parameters can be left to be determined by **sdimaging**, or they can be determined using CASA tools.\n",
        "\n",
        "\n",
        "**Image dimensions and pixel interval**\n",
        "\n",
        "The image extent can be explicitly determined using **aU.getTPSampling**:\n",
        "\n",
        "```\n",
        "xSampling, ySampling, maxsize = aU.getTPSampling(refvis, showplot=True)\n",
        "```\n",
        "\n",
        "which returns an image output showing the scans, their sky angles, and positions in RA-Dec, as shown here:\n",
        "\n",
        "![683703caa9ae3cec3d89b02760ba73df66041b2d](https://github.com/casangi/casadocs/blob/master/docs/notebooks/media/683703caa9ae3cec3d89b02760ba73df66041b2d.png?raw=1)\n",
        "\n",
        "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.\n",
        "\n",
        "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](https://safe.nrao.edu/wiki/bin/view/ALMA/PrimaryBeamArcsec)).\n",
        "\n",
        "```\n",
        "freq=115.27e+9\n",
        "fwhmfactor=1.13\n",
        "diameter=12\n",
        "theorybeam = aU.primaryBeamArcsec(frequency=freq*1e-9, fwhmfactor=fwhmfactor, diameter=diameter)\n",
        "cell = theorybeam/9.0\n",
        "imsize = int(round(maxsize/cell)*2)\n",
        "```\n",
        "\n",
        "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.\n",
        "\n",
        "```\n",
        "phasecenter='J2000 12h22m54.9 +15d49m15'\n",
        "```\n",
        "\n",
        "**Frequency and/or velocity axis**\n",
        "\n",
        "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:\n",
        "\n",
        "```\n",
        "msmd.open(vislist[0])\n",
        "freq = msmd.meanfreq(spw)\n",
        "msmd.close()\n",
        "print \"SPW %d: %.3f GHz\" % (spw, freq*1e-9)\n",
        "```\n",
        "\n",
        "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:\n",
        "\n",
        "```\n",
        "veltype='radio',\n",
        "outframe='lsrk',\n",
        "```\n",
        "\n",
        "and the rest frequency can be specified with:\n",
        "\n",
        "```\n",
        "restfreq='115.271204GHz'\n",
        "```\n",
        "\n",
        "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:\n",
        "\n",
        "```\n",
        "nchan=70,\n",
        "mode='velocity',\n",
        "start='1400km/s',\n",
        "width='5km/s',\n",
        "```\n",
        "\n",
        "**Gridding parameters**\n",
        "\n",
        "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](../api/casatasks.rst) and Mangum et al. 2007 [\\[1\\]](#Bibliography) for more details on these parameters).\n",
        "\n",
        "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. \n",
        "\n",
        "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.\n",
        "\n",
        "**Example script**\n",
        "\n",
        "Fully specified, a call to **sdimaging** might look like the following:\n",
        "\n",
        "```\n",
        "sdimaging(infiles=sd_ms+'.bl',\n",
        "    field='M42',\n",
        "    spw='%s'%(spw),\n",
        "    nchan=70,\n",
        "    mode='velocity',\n",
        "    start='1400km/s',\n",
        "    width='5km/s',\n",
        "    veltype='radio',\n",
        "    outframe='lsrk',\n",
        "    restfreq='%sGHz'%(freq/1e+9),\n",
        "    gridfunction='SF',\n",
        "    convsupport=6,\n",
        "    stokes='I',\n",
        "    phasecenter='J2000 12h22m54.9 +15d49m15',\n",
        "    ephemsrcname='',\n",
        "    imsize=imsize,\n",
        "    cell='%sarcsec'%(cell),\n",
        "    overwrite=True,\n",
        "    outfile=imagename)\n",
        "```\n",
        "\n",
        "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.\n",
        "\n",
        " \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "37ZA-RI4kunT"
      },
      "source": [
        "## Bibliography\n",
        "\n",
        "1. Mangum, et al. 2007, A&A, 474, 679-687 (http://www.aanda.org/articles/aa/pdf/2007/41/aa7811-07.pdf)\n",
        "\n",
        "\n",
        "\n"
      ]
    }
  ]
}