Radio Beam

radio-beam provides a tools for manipulating and utilizing two-dimensional Gaussian beams within the astropy framework. It is primarily built for handling radio astronomy data and is integrated into the spectral-cube <https://spectral-cube.readthedocs.io>_ package, amongst others.

radio-beam also handles operations on sets of beams, for example from a spectral cube with varying resolution in spectral channels. Of note are the algorithms for identifying the smallest common beam in a set (i.e., the minimum enclosing ellipse area).

Getting started

Basic Examples

Handling a single beam

Beam handles operations on individual beams.

Read a beam from a FITS header:

>>> from radio_beam import Beam
>>> from astropy.io import fits
>>> header = fits.getheader('file.fits')  
>>> my_beam = Beam.from_fits_header(header)  
>>> print(my_beam)  
Beam: BMAJ=0.038652855902928 arcsec BMIN=0.032841067761183604 arcsec BPA=32.29655838013 deg

To add the beam parameters to a FITS header:

>>> header.update(my_beam.to_header_keywords())  

This will return new or add values for the BMAJ, BMIN, and BPA keywords.

Create a new circular beam:

>>> from astropy import units as u
>>> my_beam = Beam(0.5*u.arcsec)
>>> my_beam
Beam: BMAJ=0.5 arcsec BMIN=0.5 arcsec BPA=0.0 deg

Beam assumes a circular beam when a minor full-width-half-max (FWHM) is not given. To create an elliptical beam, the minor FWHM and position angle need to be given:

>>> my_beam_ellip = Beam(major=0.5*u.arcsec, minor=0.25*u.arcsec, pa=30*u.deg)
>>> my_beam_ellip
Beam: BMAJ=0.5 arcsec BMIN=0.25 arcsec BPA=30.0 deg

The beam area in steradians is:

>>> my_beam_ellip.sr 
<Quantity 3.3290795e-12 sr>

Or projected into physical units given a distance:

>>> my_beam_ellip.beam_projected_area(840*u.kpc).to(u.pc**2) 
<Quantity 2.3489985 pc2>

A common unit conversion in radio astronomy is Jy/beam to K, which depends on the beam area. A Beam object can be used for this unit conversion with units:

>>> (1*u.Jy).to(u.K, u.brightness_temperature(25*u.GHz, my_beam)) 
<Quantity 7821.572919292681 K>

Or alternatively with:

>>> (1*u.Jy).to(u.K, my_beam.jtok_equiv(25*u.GHz)) 
<Quantity 7821.572919292681 K>

To get the value of 1 Jy in K for a given beam:

>>> my_beam.jtok(25*u.GHz) 
<Quantity 7821.572919292681 K>

Two beams can be convolved:

>>> my_asymmetric_beam = Beam(0.75*u.arcsec, 0.25*u.arcsec, 0*u.deg)
>>> my_other_asymmetric_beam = Beam(0.75*u.arcsec, 0.25*u.arcsec, 90*u.deg)
>>> my_asymmetric_beam.convolve(my_other_asymmetric_beam) 
Beam: BMAJ=0.790569415042 arcsec BMIN=0.790569415042 arcsec BPA=45.0 deg

And also deconvolved:

>>> my_big_beam = Beam(1.0*u.arcsec, 1.0*u.arcsec, 0*u.deg)
>>> my_little_beam = Beam(0.5*u.arcsec, 0.5*u.arcsec, 0*u.deg)
>>> my_big_beam.deconvolve(my_little_beam) 
Beam: BMAJ=0.866025403784 arcsec BMIN=0.866025403784 arcsec BPA=0.0 deg

An error is returned if the beam area is too small to deconvolve from the other:

>>> my_little_beam.deconvolve(my_big_beam)  

To find the smallest common beam between any two beams:

>>> my_asymmetric_beam.commonbeam_with(my_other_asymmetric_beam) 
Beam: BMAJ=0.75 arcsec BMIN=0.75 arcsec BPA=90.0 deg

Handling a sets of beams

Beams handles operations on sets of beams.

To read a table of beams from a FITS table:

>>> from radio_beam import Beams
>>> from astropy.io import fits
>>> bin_hdu = fits.open('file.fits')[1]  
>>> beams = Beams.from_fits_bintable(bin_hdu)  

In the above example, the second FITS extension contains the beam tables, while the first would have the spectral cube data.

To read a table of beams from a CASA image (must be run inside a CASA environment!):

>>> beams = Beams.from_casa_image('file.image')  

Create a table of beams:

>>> my_beams = Beams([1.5, 1.3] * u.arcsec, [1., 1.2] * u.arcsec, [0, 50] * u.deg)

Beams acts like a numpy array and can be sliced in the same way:

>>> my_beams[0]
Beam: BMAJ=1.5 arcsec BMIN=1.0 arcsec BPA=0.0 deg
>>> my_beams[1]
Beam: BMAJ=1.3 arcsec BMIN=1.2 arcsec BPA=50.0 deg

Find the largest beam in the set:

>>> my_beams.largest_beam()
Beam: BMAJ=1.3 arcsec BMIN=1.2 arcsec BPA=50.0 deg

Find the smallest common beam for the set (see here for more on common beams):

>>> my_beams.common_beam() 
Beam: BMAJ=1.50671729431 arcsec BMIN=1.25695643792 arcsec BPA=6.69089813778 deg

Return the smallest and largest beams in a set (by beam area):

>>> smallest_beam, largest_beam = my_beams.extrema_beams()
>>> smallest_beam
Beam: BMAJ=1.5 arcsec BMIN=1.0 arcsec BPA=0.0 deg
>>> largest_beam
Beam: BMAJ=1.3 arcsec BMIN=1.2 arcsec BPA=50.0 deg

Optionally mask out a beam (to exclude from the calculation):

>>> import numpy as np
>>> beam_mask = np.array([True, False])
>>> smallest_beam, largest_beam = my_beams.extrema_beams(includemask=beam_mask)
>>> smallest_beam
Beam: BMAJ=1.5 arcsec BMIN=1.0 arcsec BPA=0.0 deg
>>> largest_beam
Beam: BMAJ=1.5 arcsec BMIN=1.0 arcsec BPA=0.0 deg

This masking can be applied to most operations, including common_beam to exclude large outliers in the set. One useful example is if a channel is blanked in a spectral-cube, and the beam is a NaN. To make a mask to select only finite beams:

>>> my_beams.isfinite
array([ True,  True])