Merge pull request #433 from e-koch/attach_custom_beam

Attach custom beam

CHANGES.rst

@@ -1,6 +1,8 @@
 0.4.1 (unreleased)
 ------------------
-
+ - Add SpectralCube.with_beam and Projection.with_beam for attaching
+   beam objects. Raise error for position-spectral slices of VRSCs
+   (https://github.com/radio-astro-tools/spectral-cube/pull/433)
  - Raise a nicer error if no data is present in the default or
    selected HDU
    (https://github.com/radio-astro-tools/spectral-cube/pull/424)

docs/beam_handling.rst

@@ -22,6 +22,33 @@ loading the entire cube into memory!::
    >>> kcube.unit  # doctest: +SKIP
    Unit("K")
 
+
+Adding a Beam
+-------------
+
+If your cube does not have a beam, a custom beam can be attached given::
+
+    >>> new_beam = Beam(1. * u.deg)  # doctest: +SKIP
+    >>> new_cube = cube.with_beam(new_beam)  # doctest: +SKIP
+    >>> new_cube.beam  # doctest: +SKIP
+    Beam: BMAJ=3600.0 arcsec BMIN=3600.0 arcsec BPA=0.0 deg
+
+This is handy for synthetic observations, which initially have a point-like beam::
+
+    >>> point_beam = Beam(0 * u.deg)  # doctest: +SKIP
+    >>> new_cube = synth_cube.with_beam(point_beam)  # doctest: +SKIP
+    Beam: BMAJ=0.0 arcsec BMIN=0.0 arcsec BPA=0.0 deg
+
+The cube can then be convolved to a new resolution::
+
+    >>> new_beam = Beam(60 * u.arcsec)  # doctest: +SKIP
+    >>> conv_synth_cube = synth_cube.convolve_to(new_beam)  # doctest: +SKIP
+    >>> conv_synth_cube.beam  # doctest: +SKIP
+    Beam: BMAJ=60.0 arcsec BMIN=60.0 arcsec BPA=0.0 deg
+
+Beam can also be attached in the same way for `~spectral_cube.Projection` and
+`~spectral_cube.Slice` objects.
+
 Multi-beam cubes
 ----------------
 

docs/nitpick-exceptions

@@ -1,5 +1,6 @@
 py:class spectral_cube.spectral_cube.BaseSpectralCube
 py:obj radio_beam.Beam
+py:obj Beam
 py:obj astroquery.splatalogue.Splatalogue
 py:class spectral_cube.base_class.BaseNDClass
 py:class spectral_cube.base_class.SpectralAxisMixinClass

docs/rtd-pip-requirements

@@ -2,3 +2,4 @@
 numpy
 Cython
 -e git+http://github.com/astropy/astropy.git#egg=astropy
+-e git+http://github.com/radio-astro-tools/radio-beam.git#egg=radio_beam

spectral_cube/cube_utils.py

@@ -17,6 +17,7 @@
 from astropy import units as u
 import itertools
 import re
+from radio_beam import Beam
 
 
 def _fix_spectral(wcs):
@@ -267,8 +268,6 @@ def try_load_beam(header):
     Try loading a beam from a FITS header.
     '''
 
-    from radio_beam import Beam
-
     try:
         beam = Beam.from_fits_header(header)
         return beam
@@ -304,7 +303,7 @@ def try_load_beams(data):
                 if 'BPA' in hdu_item.data.names:
                     beam_table = hdu_item.data
                     return beam_table
-        
+
         try:
             # if there was a beam in a header, but not a beam table
             return beam

spectral_cube/lower_dimensional_structures.py

@@ -10,6 +10,7 @@
 from astropy import wcs
 #from astropy import log
 from astropy.io.fits import Header, Card, HDUList, PrimaryHDU
+from radio_beam import Beam
 
 from .io.core import determine_format
 from . import spectral_axis
@@ -223,26 +224,54 @@ def __new__(cls, value, unit=None, dtype=None, copy=True, wcs=None,
         else:
             self._header = Header()
 
-        if beam is not None:
-            self._beam = beam
-            self.meta['beam'] = beam
-        else:
+        if beam is None:
             if "beam" in self.meta:
-                self._beam = self.meta['beam']
+                beam = self.meta['beam']
             elif read_beam:
                 beam = cube_utils.try_load_beam(header)
-                if beam is not None:
-                    self._beam = beam
-                    self.meta['beam'] = beam
-                else:
+                if beam is None:
                     warnings.warn("Cannot load beam from header.")
 
+        if beam is not None:
+            self.beam = beam
+            self.meta['beam'] = beam
+            # TODO: Enable header updating when non-celestial slices are
+            # properly handled in the WCS object.
+            # self._header.update(beam.to_header_keywords())
+
+        return self
+
+    def with_beam(self, beam):
+        '''
+        Attach a new beam object to the Projection.
+
+        Parameters
+        ----------
+        beam : `~radio_beam.Beam`
+            A new beam object.
+        '''
+
+        meta = self.meta.copy()
+        meta['beam'] = beam
+
+        self = Projection(self.value, unit=self.unit, wcs=self.wcs,
+                          meta=meta, header=self.header,
+                          beam=beam)
+
         return self
 
     @property
     def beam(self):
         return self._beam
 
+    @beam.setter
+    def beam(self, obj):
+
+        if not isinstance(obj, Beam):
+            raise TypeError("beam must be a radio_beam.Beam object.")
+
+        self._beam = obj
+
     @staticmethod
     def from_hdu(hdu):
         '''

spectral_cube/spectral_cube.py

@@ -24,6 +24,8 @@
 
 import numpy as np
 
+from radio_beam import Beam
+
 from . import cube_utils
 from . import wcs_utils
 from . import spectral_axis
@@ -187,14 +189,6 @@ def _new_cube_with(self, data=None, wcs=None, mask=None, meta=None,
 
         return cube
 
-    def _attach_beam(self):
-
-        beam = cube_utils.try_load_beam(self.header)
-
-        if beam is not None:
-            self.beam = beam
-            self._meta['beam'] = beam
-
     @property
     def unit(self):
         """ The flux unit """
@@ -2275,7 +2269,7 @@ class SpectralCube(BaseSpectralCube):
 
     def __init__(self, data, wcs, mask=None, meta=None, fill_value=np.nan,
                  header=None, allow_huge_operations=False, beam=None,
-                 read_beam=True, wcs_tolerance=0.0, **kwargs):
+                 wcs_tolerance=0.0, **kwargs):
 
         super(SpectralCube, self).__init__(data=data, wcs=wcs, mask=mask,
                                            meta=meta, fill_value=fill_value,
@@ -2285,36 +2279,18 @@ def __init__(self, data, wcs, mask=None, meta=None, fill_value=np.nan,
                                            **kwargs)
 
         # Beam loading must happen *after* WCS is read
-        if beam is None and read_beam:
-            self._attach_beam()
-
-        if beam is None and not read_beam and 'BUNIT' in self._meta:
-            bunit = re.sub("\s", "", self._meta['BUNIT'].lower())
-            if bunit == 'jy/beam':
-                warnings.warn("Units are in Jy/beam. Attempting to parse "
-                              "header for beam information.")
-
-                self._attach_beam()
 
-                if hasattr(self, 'beam') or hasattr(self, 'beams'):
-                    warnings.warn("Units were Jy/beam.  The 'beam' is now "
-                                  "stored in the .beam attribute, and the "
-                                  "units are set to Jy")
-                else:
-                    warnings.warn("Could not parse Jy/beam unit.  Either "
-                                  "you should install the radio_beam "
-                                  "package or manually replace the units."
-                                  " For now, the units are being interpreted "
-                                  "as Jy.")
+        if beam is None:
+            beam = cube_utils.try_load_beam(self.header)
+        else:
+            if not isinstance(beam, Beam):
+                raise TypeError("beam must be a radio_beam.Beam object.")
 
         if beam is not None:
             self.beam = beam
             self._meta['beam'] = beam
-
-            # Ensure that the beam is properly defined in the header
             self._header.update(beam.to_header_keywords())
 
-        if 'beam' in self._meta:
             self.pixels_per_beam = (self.beam.sr /
                                     (astropy.wcs.utils.proj_plane_pixel_area(self.wcs) *
                                      u.deg**2)).to(u.dimensionless_unscaled).value
@@ -2328,6 +2304,30 @@ def _new_cube_with(self, **kwargs):
 
     _new_cube_with.__doc__ = BaseSpectralCube._new_cube_with.__doc__
 
+    def with_beam(self, beam):
+        '''
+        Attach a beam object to the `~SpectralCube`.
+
+        Parameters
+        ----------
+        beam : `~radio_beam.Beam`
+            `Beam` object defining the resolution element of the
+            `~SpectralCube`.
+        '''
+
+        if not isinstance(beam, Beam):
+            raise TypeError("beam must be a radio_beam.Beam object.")
+
+        meta = self._meta.copy()
+        meta['beam'] = beam
+
+        header = self._header.copy()
+        header.update(beam.to_header_keywords())
+
+        newcube = self._new_cube_with(meta=self.meta, beam=beam)
+
+        return newcube
+
     def spatial_smooth_median(self, ksize, **kwargs):
         """
         Smooth the image in each spatial-spatial plane of the cube using a median filter.
@@ -2699,8 +2699,7 @@ def convolve_to(self, beam, convolve=convolution.convolve_fft):
                                        normalize_kernel=True)
             pb.update()
 
-        newcube = self._new_cube_with(data=newdata, read_beam=False,
-                                      beam=beam)
+        newcube = self._new_cube_with(data=newdata, beam=beam)
 
         return newcube
 
@@ -2740,7 +2739,6 @@ def __init__(self, *args, **kwargs):
             'pa':5.0}``
         """
         # these types of cube are undefined without the radio_beam package
-        from radio_beam import Beam
 
         beam_table = kwargs.pop('beam_table', None)
         beams = kwargs.pop('beams', None)
@@ -2853,6 +2851,16 @@ def __getitem__(self, view):
 
             # Slice objects know how to parse Beam objects stored in the
             # metadata
+            # A 2D slice with a VRSC should not be allowed along a
+            # position-spectral axis
+            if not isinstance(self.beams[specslice], Beam):
+                raise AttributeError("2D slices along a spectral axis are not "
+                                     "allowed for "
+                                     "VaryingResolutionSpectralCubes. Convolve"
+                                     " to a common resolution with "
+                                     "`convolve_to` before attempting "
+                                     "position-spectral slicing.")
+
             meta['beam'] = self.beams[specslice]
             return Slice(value=self.filled_data[view],
                          wcs=newwcs,
@@ -2995,7 +3003,6 @@ def identify_bad_beams(self, threshold, reference_beam=None,
         includemask : np.array
             A boolean array where ``True`` indicates the good beams
         """
-        from radio_beam import Beam
 
         includemask = np.ones(len(self.beams), dtype='bool')
 
@@ -3233,7 +3240,6 @@ def convolve_to(self, beam, allow_smaller=False,
                                meta=self.meta, fill_value=self.fill_value,
                                header=self.header,
                                allow_huge_operations=self.allow_huge_operations,
-                               read_beam=False,
                                beam=beam,
                                wcs_tolerance=self._wcs_tolerance)
 

spectral_cube/tests/test_projection.py

@@ -271,6 +271,32 @@ def test_projection_with_beam():
     assert new_proj.beam == exp_beam
     assert new_proj.meta['beam'] == exp_beam
 
+    # load beam from beam object
+    exp_beam = Beam(3.0 * u.arcsec)
+    header = hdu.header.copy()
+    del header["BMAJ"], header["BMIN"], header["BPA"]
+    new_proj = Projection(hdu.data, wcs=proj.wcs, header=header,
+                          beam=exp_beam)
+
+    assert new_proj.beam == exp_beam
+    assert new_proj.meta['beam'] == exp_beam
+
+
+def test_projection_attach_beam():
+
+    exp_beam = Beam(1.0 * u.arcsec)
+    newbeam = Beam(2.0 * u.arcsec)
+
+    proj, hdu = load_projection("55.fits")
+
+    new_proj = proj.with_beam(newbeam)
+
+    assert proj.beam == exp_beam
+    assert proj.meta['beam'] == exp_beam
+
+    assert new_proj.beam == newbeam
+    assert new_proj.meta['beam'] == newbeam
+
 
 @pytest.mark.parametrize(('LDO', 'data'),
                          zip(LDOs_2d, data_two_2d))