Use common beam instead of beam averaging

docs/beam_handling.rst

@@ -53,19 +53,122 @@ Multi-beam cubes
 ----------------
 
 Varying resolution (multi-beam) cubes are somewhat trickier to work with in
-general, though unit conversion is easy.  You can perform the same sort of unit
-conversion with `~spectral_cube.spectral_cube.VaryingResolutionSpectralCube` s
-as with regular `~spectral_cube.spectral_cube.SpectralCube` s; ``spectral-cube``
-will use a different beam and frequency for each plane.
+general, though some operations are similar to a single resolution cube.
+Unit conversion is one case where the operation is called the same for
+`~spectral_cube.spectral_cube.VaryingResolutionSpectralCube` s
+as with regular `~spectral_cube.spectral_cube.SpectralCube` s. For example, to
+convert from Jy / beam to K::
 
-You can identify channels with bad beams (i.e., beams that differ from a reference beam,
-which by default is the median beam) using
+    >>> vrsc_cube_K = vrsc_cube.to(u.K)  # doctest: +SKIP
+
+``VaryingResolutionSpectralCube`` will use a different beam and frequency for each plane.
+
+Handling variations in the beams is often a source of difficulty. Some spectral
+operations (e.g., moment maps) require a common resolution. However, a few channels
+may have large discprepancies in the beam compared to most of the others (e.g., if
+more data in a particular channel was bad and had to be flagged).
+`~spectral_cube.spectral_cube.VaryingResolutionSpectralCube` deals with these problems in two ways.
+
+Finding and convolving to a common resolution
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+When applying spectral operations on a
+`~spectral_cube.spectral_cube.VaryingResolutionSpectralCube`, it is often
+necessary to first convolve the data to have a common beam. To retain the
+maximum amount of spatial information, the data can be convolved to the
+smallest common beam.
+
+Tools for finding the smallest common beam and convolution operations are described
+in detail in :doc:`smoothing`.
+
+
+Small beam variations: Do I need to convolve first?
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Formally, operations on multiple spectral channels require the data to have a common
+resolution. In practice, data cubes with a high spectral resolution will have the beam
+vary by a very small factor between adjacent channels. It is then a safe approximation
+to ignore small changes in the beam to avoid having to convolve the entire data cube to
+a common beam (which might be a slow and expensive operation for large cubes).
+
+`~spectral_cube.spectral_cube.VaryingResolutionSpectralCube` allows for small variations
+in the beam by setting
+:meth:`~spectral_cube.spectral_cube.VaryingResolutionSpectralCube.beam_threshold`. The
+`beam_threshold` is the largest allowed fractional change relative to the smallest common beam that will be allowed in spectral operations. For example, the default 0.01 allows up
+to a 1% variation from the smallest common beam. If the beams vary more than this, a
+`ValueError` will be raised prior to the spectral operation. Otherwise, variations between
+the beams will be ignored, a warning will be printed describing this approximation, and
+the spectral operation will be computed.
+
+The beam threshold can be changed from the default 0.01 with::
+
+    >>> vrsc_cube.beam_threshold = 0.02   # doctest: +SKIP
+
+.. note::
+
+    Setting `vrsc_cube.beam_threshold = 1.0`, or greater, will allow for large
+    beam variations without raising an error. Large values for the beam threshold
+    should not be used for scientific results.
+
+For most spectral-line data cubes covering a single spectral line, the fine resolution
+and small frequency range will often allow the above approximation to work.
+However, if you are working with wideband data, this approximation will not hold and you
+will need to convolve to a common beam or extract spectral slabs from the data cube.
+
+Strict mode
+***********
+
+To disable allowing small variations in the beam,
+`~spectral_cube.spectral_cube.VaryingResolutionSpectralCube` has a strict checking mode
+than can be enabled with::
+
+    >>> vrsc_cube.strict_beam_mode = True  # doctest: +SKIP
+
+When strict mode is enabled, all beam must be identical for spectral operations to work.
+Note that this is equivalent to setting the beam threshold to 0.
+
+
+Identifying channels with bad beams
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Some cubes may have channels with beams that vary drastically over small ranges in spectral
+channels. This is often the case where a range of channels has poor data or is affected
+by radio frequency interference, leading to most of the data in that channel being flagged.
+If these channels are kept, the smallest common beam (see :doc:`smoothing`) may be much
+larger due to these channels.
+
+You can identify and mask channels with bad beams using
+`~spectral_cube.spectral_cube.VaryingResolutionSpectralCube.identify_bad_beams`::
+
+    >>> goodbeams = vrsc_cube.identify_bad_beams()  # doctest: +SKIP
+
+This will return a 1D boolean mask where `True` means the channel beam is good.
+By default,
 `~spectral_cube.spectral_cube.VaryingResolutionSpectralCube.identify_bad_beams`
-(the returned value is a mask array where ``True`` means the channel is good),
-mask channels with undesirable beams using
-`~spectral_cube.spectral_cube.VaryingResolutionSpectralCube.mask_out_bad_beams`,
-and in general mask out individual channels using
+will use
+`~spectral_cube.spectral_cube.VaryingResolutionSpectralCube.beam_threshold` (default
+of 0.01; see above).
+However, the comparison here is to the **median** major and minor axis rather than
+the smallest common beam used above.
+This is because bad beams are identified as outliers in the set of beams.
+
+To mask the channels with bad beams, use
+`~spectral_cube.spectral_cube.VaryingResolutionSpectralCube.with_bad_beams_masked`.
+
+
+    >>> masked_vrsc_cube = vrsc_cube.with_bad_beams_masked()  # doctest: +SKIP
+
+The masked cube without the bad beams will now exlude channels with bad beams and
+can be used, for example, to convolve to a better representative common beam
+resolution (see above).
+
+We also note that, in general, you can mask out individual channels using
 `~spectral_cube.spectral_cube.VaryingResolutionSpectralCube.mask_channels`.
 
-For other sorts of operations, discussion of how to deal with these cubes via
-smoothing to a common resolution is in the :doc:`smoothing` document.
+.. note::
+
+    The common beam is not used to find discrepant and bad beams since they are
+    identified as outliers from the set. We note that this is an approximate method of
+    finding channels with outlier beams and may fail in some cases. Please
+    `raise an issue <https://github.com/radio-astro-tools/spectral-cube/issues>`_ if
+    this method does not work well for your data cube.

docs/errors.rst

@@ -20,7 +20,7 @@ operations spanning the spectral axis, for example ``cube.moment0(axis=0)`` or
 This occurs if the beam sizes are different by more than the specified
 threshold factor.  A default threshold of 1% is set because for most
 interferometric images, beam differences on this scale are negligible (they
-correspond to flux measurement errors of 10^-4).  
+correspond to flux measurement errors of 10^-4).
 
 To inspect the beam properties, look at the ``beams`` attribute, for example:
 
@@ -54,9 +54,12 @@ There are several options to manage this problem:
 .. code::
 
    good_beams = cube.identify_bad_beams(threshold=0.1)
-   mcube = cube.mask_out_bad_beams(threshold=0.1)
+   mcube = cube.with_bad_beams_masked(threshold=0.1)
 
 
+Please see the more detailed explanation of the beam threshold treatment in
+:doc:`beam_handling` for more information about this warning.
+
 
 
 Moment-2 or FWHM calculations give unexpected NaNs

docs/smoothing.rst

@@ -29,12 +29,31 @@ will be different for each slice.
 
 Common Beam selection
 ^^^^^^^^^^^^^^^^^^^^^
-You may want to convolve your cube to the smallest beam that is still larger
-than all contained beams.  To do this, you can use the
-`~radio_beam.Beams.common_beam` tool.  For example::
 
-    common_beam = cube.beams.common_beam()
-    new_cube = cube.convolve_to(common_beam)
+Many spectral operations on a cube with a varying resolution will require first
+convolving to a common resolution. You may want to retain the finest spatial
+resolution possible by finding the smallest common beam that encloses all beams
+in the cube. To do this, we can use
+:meth:`~spectral_cube.VaryingResolutionSpectralCube.set_common_beam`::
+
+    cube.set_common_beam()
+
+This method wraps the `~radio_beam.Beams.common_beam`, and keyword arguments for
+the common beam algorithm can be passed directly in
+:meth:`~spectral_cube.VaryingResolutionSpectralCube.set_common_beam`.
+The common beam can then be accessed with::
+
+    cube.common_beam
+
+A new common beam will be computed each time
+:meth:`~spectral_cube.VaryingResolutionSpectralCube.set_common_beam` is called.
+This can be used to make changes to the settings for the common beam algorithm
+or if some beams are masked (see :doc:`beam_handling`).
+
+The cube can then be convolved to the smallest common beam using the normal
+call for convolution::
+
+    new_cube = cube.convolve_to(cube.common_beam)
 
 Sometimes, you'll encounter the error "Could not find common beam to deconvolve
 all beams." This is a real issue, as the algorithms we have in hand so far do
@@ -45,7 +64,7 @@ algorithm to converge to a valid common beam:
 `~radio_beam.commonbeam.getMinVolEllipse` code by
 passing ``tolerance=1e-5`` to the common beam function::
 
-    cube.beams.common_beam(tolerance=1e-5)
+    cube.set_common_beam(tolerance=1e-5)
 
 Convergence may be met by either increasing or decreasing the tolerance; it
 depends on having the algorithm not step within the minimum enclosing ellipse,
@@ -57,7 +76,7 @@ and will take longer to run.
 to overestimate the beam size, ensuring that solutions that are marginally
 smaller than the common beam will not be found by the algorithm::
 
-    cube.beams.common_beam(epsilon=1e-3)
+    cube.set_common_beam(epsilon=1e-3)
 
 The default value of ``epsilon=1e-3`` will sample points 0.1% larger than the
 edge of each beam in the set. Increasing ``epsilon`` ensures that a valid common

setup.cfg

@@ -20,6 +20,7 @@ install_requires =
     dask[array]
     joblib
     casa-formats-io
+    scipy
 
 [options.extras_require]
 test =
@@ -34,7 +35,6 @@ novis =
     pvextractor
     regions ; python_version<'3.8'
     reproject
-    scipy
 all =
     zarr
     fsspec