convolve_rs._convolve_rs¶

Classes¶

Beam

A 2-D Gaussian representation of a radio telescope's PSF (beam).

Functions¶

`common_beam`(→ Beam)	Find the smallest beam that every beam in `beams` can be convolved to.
`gauss_factor`(→ tuple[float, float, float, float, float])	Compute the flux-scaling factor for a Jy/beam convolution.
`smooth`(→ Any)	Smooth an image from `old_beam` to `new_beam`.

Package Contents¶

class convolve_rs._convolve_rs.Beam[source]¶

A 2-D Gaussian representation of a radio telescope’s PSF (beam).

All axes use FITS conventions: FWHM major and minor axes in degrees, position angle in degrees East of North.

Parameters:

major_deg – FWHM major axis in degrees (FITS BMAJ).
minor_deg – FWHM minor axis in degrees (FITS BMIN). Must be <= major_deg.
pa_deg – Position angle in degrees East of North (FITS BPA).

Raises:

ValueError – If minor_deg > major_deg or any value is non-finite.

See also

Beam.from_arcsec: Construct from arcsecond axes. Beam.from_fits_header: Construct from an astropy FITS header. Beam.from_radio_beam: Construct from a radio_beam.Beam object.

Examples

>>> from convolve_rs import Beam
>>> beam = Beam(0.005, 0.004, 30.0)
>>> beam.major_deg
0.005
>>> round(beam.major_arcsec, 6)
18.0
>>> Beam(0.004, 0.005, 0.0)
Traceback (most recent call last):
...
ValueError: invalid beam: minor axis (0.005) > major axis (0.004)

__repr__() → str[source]¶

__str__() → str[source]¶

area_sr() → float[source]¶

Solid angle of the beam in steradians.

Computed as (pi / (4 ln 2)) * major_rad * minor_rad.

Returns:: Beam solid angle in steradians.
Return type:: float

convolve(other: Beam) → Beam[source]¶

Convolve self with other.

Adds the Gaussian covariance matrices and reads off the resulting ellipse (Wild 1970).

Parameters:: other (Beam) – The beam to convolve with.
Returns:: The convolved beam.
Return type:: Beam

Examples

Convolving two circular beams adds their axes in quadrature:

>>> from convolve_rs import Beam
>>> a = Beam.from_arcsec(3.0, 3.0, 0.0)
>>> b = Beam.from_arcsec(4.0, 4.0, 0.0)
>>> round(a.convolve(b).major_arcsec, 6)
5.0

deconvolve(other: Beam) → Beam[source]¶

Deconvolve other from self (i.e. self = result ⊛ other).

Subtracts the Gaussian covariance matrices and reads off the residual ellipse (Wild 1970).

Parameters:: other (Beam) – The PSF to deconvolve from this beam.
Returns:: The deconvolved beam.
Return type:: Beam
Raises:: ValueError – If other is larger than self and deconvolution is impossible.

Examples

Deconvolution inverts convolution:

>>> from convolve_rs import Beam
>>> a = Beam.from_arcsec(3.0, 3.0, 0.0)
>>> b = Beam.from_arcsec(4.0, 4.0, 0.0)
>>> c = a.convolve(b)
>>> round(c.deconvolve(a).major_arcsec, 6)
4.0
>>> b.deconvolve(c)
Traceback (most recent call last):
...
ValueError: beam could not be deconvolved: source beam is smaller than the PSF

classmethod from_arcsec(major_arcsec: float, minor_arcsec: float, pa_deg: float) → Beam[source]¶

Construct a Beam from arcsecond axes.

Parameters:

major_arcsec – FWHM major axis in arcseconds.
minor_arcsec – FWHM minor axis in arcseconds. Must be <= major_arcsec.
pa_deg – Position angle in degrees East of North.

Returns:

The constructed beam.

Return type:

Beam

Raises:

ValueError – If minor_arcsec > major_arcsec or any value is non-finite.

Examples

>>> from convolve_rs import Beam
>>> beam = Beam.from_arcsec(18.0, 14.4, 30.0)
>>> beam.major_deg
0.005

property major_arcsec: float[source]¶: FWHM major axis in arcseconds.

property major_deg: float[source]¶: FWHM major axis in degrees (FITS BMAJ).

property minor_arcsec: float[source]¶: FWHM minor axis in arcseconds.

property minor_deg: float[source]¶: FWHM minor axis in degrees (FITS BMIN).

property pa_deg: float[source]¶: Position angle in degrees East of North (FITS BPA).

convolve_rs._convolve_rs.common_beam(beams: Sequence[Beam], tolerance: float = 0.0001, nsamps: int = 200, epsilon: float = 0.0005) → Beam[source]¶

Find the smallest beam that every beam in beams can be convolved to.

Uses the 2-beam analytic CASA algorithm when len(beams) == 2, otherwise the Khachiyan minimum-volume-enclosing-ellipse algorithm — the same as radio_beam.Beams.common_beam(method='pts').

Parameters:

beams (list[Beam]) – Input beams. Must contain at least one element.
tolerance (float) – Convergence tolerance for the Khachiyan algorithm. Default 1e-4.
nsamps (int) – Number of points sampled from each beam ellipse boundary. Default 200.
epsilon (float) – Fractional padding added to each beam before the MVE fit, to ensure the common beam can be marginally deconvolved from all inputs. Default 5e-4.

Returns:

The smallest common beam.

Return type:

Beam

Raises:

ValueError – If beams is empty or no valid common beam is found.

Examples

>>> from convolve_rs import Beam, common_beam
>>> b1 = Beam.from_arcsec(10.0, 8.0, 30.0)
>>> b2 = Beam.from_arcsec(12.0, 6.0, 60.0)
>>> cb = common_beam([b1, b2])
>>> cb.major_arcsec >= 12.0
True
>>> cb.area_sr() >= max(b1.area_sr(), b2.area_sr())
True

convolve_rs._convolve_rs.gauss_factor(conv_beam: Beam, orig_beam: Beam, dx_arcsec: float, dy_arcsec: float) → tuple[float, float, float, float, float][source]¶

Compute the flux-scaling factor for a Jy/beam convolution.

Returns the factor by which pixel values must be multiplied after convolving a Jy/beam image from orig_beam to conv_beam.

Parameters:

conv_beam (Beam) – The convolving beam (the kernel applied on top of orig_beam).
orig_beam (Beam) – The original restoring beam of the image.
dx_arcsec (float) – Pixel size along the x axis in arcseconds.
dy_arcsec (float) – Pixel size along the y axis in arcseconds.

Returns:

(fac, amp, bmaj_out, bmin_out, bpa_out_deg) where: fac is the pixel scaling factor, amp is the Gaussian kernel integral, and the remaining three are the output beam parameters (major/minor FWHM in arcseconds, PA in degrees).

Return type:

tuple

Examples

>>> from convolve_rs import Beam, gauss_factor
>>> conv = Beam.from_arcsec(5.0, 5.0, 0.0)
>>> orig = Beam.from_arcsec(10.0, 10.0, 0.0)
>>> fac, amp, bmaj, bmin, bpa = gauss_factor(conv, orig, 2.5, 2.5)
>>> fac > 0.0
True
>>> round(bmaj, 5)  # √(10² + 5²)
11.18034

convolve_rs._convolve_rs.smooth(image: Any, old_beam: Beam, new_beam: Beam, dx_deg: float, dy_deg: float, cutoff_arcsec: float | None = None, bunit: str | None = None) → Any[source]¶

Smooth an image from old_beam to new_beam.

Convolves image in the UV plane and applies the flux scaling appropriate for bunit so that the output is in the same units as the input: Jy/beam images get the Gaussian beam-area factor, Kelvin (brightness temperature) images conserve surface brightness and are left unscaled.

Parameters:

image (numpy.ndarray) – Input image, shape (ny, nx), dtype float32 or float64. The convolution runs in the input’s precision and the output keeps the same dtype.
old_beam (Beam) – Current (input) restoring beam.
new_beam (Beam) – Target (output) restoring beam. Must be larger than old_beam.
dx_deg (float) – Pixel size along the x (RA) axis in degrees (FITS CDELT1, may be negative).
dy_deg (float) – Pixel size along the y (Dec) axis in degrees (FITS CDELT2).
cutoff_arcsec (float, optional) – If given, raise ValueError if the deconvolved kernel FWHM exceeds this value in arcseconds.
bunit (str, optional) – FITS BUNIT brightness unit. If it denotes Kelvin (e.g. "K"), surface brightness is conserved and no flux scaling is applied; if it denotes Jy/beam, the Gaussian flux-scaling factor is applied. An unrecognised string emits a UserWarning and is treated as Jy/beam. Defaults to Jy/beam.

Returns:

Smoothed image, shape (ny, nx), same dtype as the: input (float32 or float64).

Return type:

numpy.ndarray

Raises:

ValueError – If new_beam is smaller than old_beam, all pixels are NaN, or the kernel exceeds cutoff_arcsec.

Warns:

UserWarning – If bunit is given but not recognised as either a Kelvin or Jy/beam unit (Jy/beam is then assumed).

Examples

Smoothing a flat Jy/beam image from a 10″ to a 20″ circular beam scales pixel values by the beam-area ratio (4); in Kelvin, surface brightness is conserved:

>>> import numpy as np
>>> from convolve_rs import Beam, smooth
>>> image = np.ones((64, 64), dtype=np.float32)
>>> old = Beam.from_arcsec(10.0, 10.0, 0.0)
>>> new = Beam.from_arcsec(20.0, 20.0, 0.0)
>>> dx = 2.5 / 3600.0
>>> jy = smooth(image, old, new, dx, dx)
>>> round(float(jy[32, 32]), 3)
4.0
>>> k = smooth(image, old, new, dx, dx, bunit="K")
>>> round(float(k[32, 32]), 3)
1.0