Gridding#

mpol.gridding.verify_no_hermitian_pairs(uu: ndarray[Any, dtype[floating[Any]]], vv: ndarray[Any, dtype[floating[Any]]], data: ndarray[Any, dtype[complexfloating[Any, Any]]], test_vis: int = 5, test_channel: int = 0) → bool[source]#

Check that the dataset does not contain Hermitian pairs. Because the sky brightness \(I_\nu\) is real, the visibility function \(\mathcal{V}\) is Hermitian, meaning that

\[\mathcal{V}(u, v) = \mathcal{V}^*(-u, -v).\]

Most datasets (e.g., those extracted from CASA) will only record one visibility measurement per baseline and not include the duplicate Hermitian pair (to save storage space). This routine attempts to determine if the dataset contains Hermitian pairs or not by choosing one data point at a time and then searching the dataset to see if its Hermitian pair exists. The routine will declare that a dataset contains Hermitian pairs or not after it has searched test_vis number of data points. If 0 Hermitian pairs have been found for all test_vis points, then the dataset will be declared to have no Hermitian pairs. If test_vis Hermitian pairs were found for test_vis points searched, then the dataset will be declared to have Hermitian pairs. If more than 0 but fewer than test_vis Hermitian pairs were found for test_vis points, an error will be raised.

Gridding objects like mpol.gridding.DirtyImager will naturally augment an input dataset to include the Hermitian pairs, so that images of \(I_\nu\) produced with the inverse Fourier transform turn out to be real.

Parameters:

uu (numpy array) – array of u spatial frequency coordinates. Units of [\(\lambda\)]
vv (numpy array) – array of v spatial frequency coordinates. Units of [\(\lambda\)]
data (numpy complex) – array of data values
test_vis (int) – the number of data points to search for Hermitian ‘matches’
test_channel (int) – the index of the channel to perform the check

Returns:

True if dataset does contain Hermitian pairs, False if not.

Return type:

boolean

class mpol.gridding.GridderBase(coords=<class 'mpol.coordinates.GridCoords'>, uu=numpy.ndarray[typing.Any, numpy.dtype[numpy.floating[typing.Any]]], vv=numpy.ndarray[typing.Any, numpy.dtype[numpy.floating[typing.Any]]], weight=numpy.ndarray[typing.Any, numpy.dtype[numpy.floating[typing.Any]]], data_re=numpy.ndarray[typing.Any, numpy.dtype[numpy.floating[typing.Any]]], data_im=numpy.ndarray[typing.Any, numpy.dtype[numpy.floating[typing.Any]]])[source]#

This class is not designed to be used directly, but rather to be subclassed.

Subclasses will need to implement a _grid_visibilities(self,…) method.

The GridderBase object uses desired image dimensions (via the cell_size and npix arguments) to define a corresponding Fourier plane grid as a GridCoords object. A pre-computed GridCoords can be supplied in lieu of cell_size and npix, but all three arguments should never be supplied at once. For more details on the properties of the grid that is created, see the GridCoords documentation.

Subclasses will accept “loose” ungridded visibility data and store the arrays to the object as instance attributes. The input visibility data should be the set of visibilities over the full \([-u,u]\) and \([-v,v]\) domain, and should not contain Hermitian pairs (an error will be raised, if they are encountered). The visibilities can be 1d for a single continuum channel, or 2d for an image cube. If 1d, visibilities will be converted to 2d arrays of shape (1, nvis). Like the ImageCube class, after construction, GridderBase assumes that you are operating with a multi-channel set of visibilities. These routines will still work with single-channel ‘continuum’ visibilities, they will just have nchan = 1 in the first dimension of most products.

Parameters:

coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
uu (numpy array) – (nchan, nvis) array of u spatial frequency coordinates. Units of [\(\lambda\)]
vv (numpy array) – (nchan, nvis) array of v spatial frequency coordinates. Units of [\(\lambda\)]
weight (2d numpy array) – (nchan, nvis) array of thermal weights. Units of [\(1/\mathrm{Jy}^2\)]
data_re (2d numpy array) – (nchan, nvis) array of the real part of the visibility measurements. Units of [\(\mathrm{Jy}\)]
data_im (2d numpy array) – (nchan, nvis) array of the imaginary part of the visibility measurements. Units of [\(\mathrm{Jy}\)]

property ground_cube: ndarray[Any, dtype[floating[Any]]]#

The visibility FFT cube fftshifted for plotting with imshow.

Returns:: the FFT of the image cube, in sky plane format.
Return type:: (torch.complex tensor, of shape (nchan, npix, npix))

class mpol.gridding.DataAverager(coords=<class 'mpol.coordinates.GridCoords'>, uu=numpy.ndarray[typing.Any, numpy.dtype[numpy.floating[typing.Any]]], vv=numpy.ndarray[typing.Any, numpy.dtype[numpy.floating[typing.Any]]], weight=numpy.ndarray[typing.Any, numpy.dtype[numpy.floating[typing.Any]]], data_re=numpy.ndarray[typing.Any, numpy.dtype[numpy.floating[typing.Any]]], data_im=numpy.ndarray[typing.Any, numpy.dtype[numpy.floating[typing.Any]]])[source]#

The DataAverager object uses desired image dimensions (via the cell_size and npix arguments) to define a corresponding Fourier plane grid as a GridCoords object. A pre-computed GridCoords can be supplied in lieu of cell_size and npix, but all three arguments should never be supplied at once. For more details on the properties of the grid that is created, see the GridCoords documentation.

The DataAverager object accepts “loose” ungridded visibility data and stores the arrays to the object as instance attributes. The input visibility data should be the set of visibilities over the full \([-u,u]\) and \([-v,v]\) domain, and should not contain Hermitian pairs (an error will be raised, if they are encountered). (Note that, unlike DirtyImager, this class will not augment the dataset to include Hermitian pairs. This is by design, since Hermitian pairs should not be used in likelihood calculations).

The input visibilities can be 1d for a single continuum channel, or 2d for image cube. If 1d, visibilities will be converted to 2d arrays of shape (1, nvis). Like the ImageCube class, after construction, the DataAverager assumes that you are operating with a multi-channel set of visibilities. These routines will still work with single-channel ‘continuum’ visibilities, they will just have nchan = 1 in the first dimension of most products.

Once the DataAverager object is initialized with loose visibilities, you can average them and export them for use in Regularized Maximum Likelihood imaging with the mpol.gridding.DataAverager.to_pytorch_dataset() routine.

Example:

averager = gridding.DataAverager(
    coords=coords,
    uu=uu,
    vv=vv,
    weight=weight,
    data_re=data_re,
    data_im=data_im,
    )

dset = averager.to_pytorch_dataset()

Parameters:

coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
uu (numpy array) – (nchan, nvis) array of u spatial frequency coordinates. Units of [\(\lambda\)]
vv (numpy array) – (nchan, nvis) array of v spatial frequency coordinates. Units of [\(\lambda\)]
weight (2d numpy array) – (nchan, nvis) array of thermal weights. Units of [\(1/\mathrm{Jy}^2\)]
data_re (2d numpy array) – (nchan, nvis) array of the real part of the visibility measurements. Units of [\(\mathrm{Jy}\)]
data_im (2d numpy array) – (nchan, nvis) array of the imaginary part of the visibility measurements. Units of [\(\mathrm{Jy}\)]

to_pytorch_dataset(check_visibility_scatter: bool = True, max_scatter: float = 1.2) → GriddedDataset[source]#

Export gridded visibilities to a PyTorch mpol.datasets.GriddedDataset object.

Parameters:

check_visibility_scatter (bool) – whether the routine should check the standard deviation of visibilities in each within each \(u,v\) cell (\(\mathrm{cell}_{i,j}\)) defined by self.coords. Default is True. A RuntimeError will be raised if any cell has a scatter larger than max_scatter.
max_scatter (float) – the maximum allowable standard deviation of visibility values in a given \(u,v\) cell (\(\mathrm{cell}_{i,j}\)) defined by self.coords. Defaults to a factor of 120%.

Returns:

GriddedDataset with gridded visibilities.

class mpol.gridding.DirtyImager(coords: GridCoords, uu: ndarray[Any, dtype[floating[Any]]], vv: ndarray[Any, dtype[floating[Any]]], weight: ndarray[Any, dtype[floating[Any]]], data_re: ndarray[Any, dtype[floating[Any]]], data_im: ndarray[Any, dtype[floating[Any]]])[source]#

This class is mainly used for producing diagnostic “dirty” images of the visibility data.

The DirtyImager object uses desired image dimensions (via the cell_size and npix arguments) to define a corresponding Fourier plane grid as a GridCoords object. A pre-computed GridCoords can be supplied in lieu of cell_size and npix, but all three arguments should never be supplied at once. For more details on the properties of the grid that is created, see the GridCoords documentation.

The DirtyImager object accepts “loose” ungridded visibility data and stores the arrays to the object as instance attributes. The input visibility data should be the normal set of visibilities over the full \([-u,u]\) and \([-v,v]\) domain; internally the DirtyImager will automatically augment the dataset to include the complex conjugates, i.e. the ‘Hermitian pairs.’

The input visibilities can be 1d for a single continuum channel, or 2d for image cube. If 1d, visibilities will be converted to 2d arrays of shape (1, nvis). Like the ImageCube class, after construction, the DirtyImager assumes that you are operating with a multi-channel set of visibilities. These routines will still work with single-channel ‘continuum’ visibilities, they will just have nchan = 1 in the first dimension of most products.

Example:

imager = gridding.DirtyImager(
    coords=coords,
    uu=uu,
    vv=vv,
    weight=weight,
    data_re=data_re,
    data_im=data_im,
    )

img, beam = imager.get_dirty_image(weighting="briggs", robust=0.0)

Parameters:

cell_size (float) – width of a single square pixel in [arcsec]
npix (int) – number of pixels in the width of the image
coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
uu (numpy array) – (nchan, nvis) array of u spatial frequency coordinates. Units of [\(\lambda\)]
vv (numpy array) – (nchan, nvis) array of v spatial frequency coordinates. Units of [\(\lambda\)]
weight (2d numpy array) – (nchan, nvis) array of thermal weights. Units of [\(1/\mathrm{Jy}^2\)]
data_re (2d numpy array) – (nchan, nvis) array of the real part of the visibility measurements. Units of [\(\mathrm{Jy}\)]
data_im (2d numpy array) – (nchan, nvis) array of the imaginary part of the visibility measurements. Units of [\(\mathrm{Jy}\)]

classmethod from_tensors(coords: GridCoords, uu: Tensor, vv: Tensor, weight: Tensor, data: Tensor) → DirtyImager[source]#: mpol.gridding.DirtyImager is written internally in numpy. However, one may wish to use mpol.gridding.DirtyImager to image residual visibilities, which are commonly of type torch.Tensor. This classmethod converts the tensors to numpy arrays and instantiates a mpol.gridding.DirtyImager instance.

get_dirty_beam_area(ntheta: int = 24, single_channel_estimate: bool = True) → ndarray[Any, dtype[floating[Any]]][source]#

Compute the effective area of the dirty beam for each channel. Assumes that the beam has already been generated by running get_dirty_image(). This is an approximate calculation involving a simple sum over all pixels out to the first null (zero crossing) of the dirty beam. This quantity is designed to approximate the conversion of image units from \([\mathrm{Jy}\,\mathrm{beam}^{-1}]\) to \([\mathrm{Jy}\,\mathrm{arcsec}^{-2}]\), even though units of \([\mathrm{Jy}\,\mathrm{dirty\;beam}^{-1}]\) are technically undefined.

Parameters:

ntheta (int) – number of azimuthal wedges to use for the 1st null calculation. More wedges will result in a more accurate estimate of dirty beam area, but will also take longer.
single_channel_estimate (bool) – If True (the default), use the area estimated from the first channel for all channels in the multi-channel image cube. If False, calculate the beam area for all channels.

Returns:

(1D numpy array float) beam area for each channel in units of \([\mathrm{arcsec}^{2}]\)

get_dirty_image(weighting: str = 'uniform', robust: float | None = None, taper_function: Callable[[ndarray[Any, dtype[floating[Any]]], ndarray[Any, dtype[floating[Any]]]], ndarray[Any, dtype[floating[Any]]]] | None = None, unit: str = 'Jy/beam', check_visibility_scatter: bool = True, max_scatter: float = 1.2, **beam_kwargs) → tuple[ndarray[Any, dtype[floating[Any]]], ndarray[Any, dtype[floating[Any]]]][source]#

Calculate the dirty image.

Parameters:

weighting (string) – The type of cell averaging to perform. Choices of "natural", "uniform", or "briggs", following CASA tclean. If "briggs", also specify a robust value.
robust (float) – If weighting='briggs', specify a robust value in the range [-2, 2]. robust=-2 approxmately corresponds to uniform weighting and robust=2 approximately corresponds to natural weighting.
taper_function (function reference) – a function assumed to be of the form \(f(u,v)\) which calculates a prefactor in the range \([0,1]\) and premultiplies the visibility data. The function must assume that \(u\) and \(v\) will be supplied in units of \(\mathrm{k}\lambda\). By default no taper is applied.
unit (string) – what unit should the image be in. Default is "Jy/beam". If "Jy/arcsec^2", then the effective area of the dirty beam will be used to convert from "Jy/beam" to "Jy/arcsec^2".
check_visibility_scatter (bool) – whether the routine should check the standard deviation of visibilities in each within each \(u,v\) cell (\(\mathrm{cell}_{i,j}\)) defined by self.coords. Default is True. A RuntimeWarning will be raised if any cell has a scatter larger than max_scatter.
max_scatter (float) – the maximum allowable standard deviation of visibility values in a given \(u,v\) cell (\(\mathrm{cell}_{i,j}\)) defined by self.coords. Defaults to a factor of 120%.
**beam_kwargs – all additional keyword arguments passed to get_dirty_beam_area() if unit="Jy/arcsec^2".

Returns:

2-tuple of (image, beam) where image is an (nchan, npix, npix) numpy array of the dirty image cube in units unit. beam is an numpy image cube with a dirty beam (PSF) for each channel. The units of the beam are always Jy/{dirty beam}, i.e., the peak of the beam is normalized to 1.0.

Gridding

Contents

Gridding#