diff --git a/CHANGES.rst b/CHANGES.rst index 71dd8324c..dbc88d039 100644 --- a/CHANGES.rst +++ b/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) diff --git a/docs/beam_handling.rst b/docs/beam_handling.rst index 8f49e15ac..e83aa8727 100644 --- a/docs/beam_handling.rst +++ b/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 ---------------- diff --git a/docs/nitpick-exceptions b/docs/nitpick-exceptions index 113ee32c8..320eba712 100644 --- a/docs/nitpick-exceptions +++ b/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 diff --git a/docs/rtd-pip-requirements b/docs/rtd-pip-requirements index e97db8007..b4f932a7e 100644 --- a/docs/rtd-pip-requirements +++ b/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 \ No newline at end of file diff --git a/spectral_cube/cube_utils.py b/spectral_cube/cube_utils.py index 177d544dd..423b019a9 100644 --- a/spectral_cube/cube_utils.py +++ b/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 diff --git a/spectral_cube/lower_dimensional_structures.py b/spectral_cube/lower_dimensional_structures.py index 9b755c4e1..a0060b7a0 100644 --- a/spectral_cube/lower_dimensional_structures.py +++ b/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): ''' diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py index 452431656..85cf8e45e 100644 --- a/spectral_cube/spectral_cube.py +++ b/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) diff --git a/spectral_cube/tests/test_projection.py b/spectral_cube/tests/test_projection.py index c244645db..001f50cf5 100644 --- a/spectral_cube/tests/test_projection.py +++ b/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)) diff --git a/spectral_cube/tests/test_spectral_cube.py b/spectral_cube/tests/test_spectral_cube.py index 0ef7676eb..abd081db5 100644 --- a/spectral_cube/tests/test_spectral_cube.py +++ b/spectral_cube/tests/test_spectral_cube.py @@ -127,12 +127,13 @@ def setup_method(self, method): ('vad', [2, 0, 1]), ('vda', [0, 2, 1]), ('adv', [0, 1, 2]), - ('vda_beams', [0, 2, 1]) ] +translist_vrsc = [('vda_beams', [0, 2, 1])] + class TestSpectralCube(object): - @pytest.mark.parametrize(('name', 'trans'), translist) + @pytest.mark.parametrize(('name', 'trans'), translist + translist_vrsc) def test_consistent_transposition(self, name, trans): """data() should return velocity axis first, then world 1, then world 0""" c, d = cube_and_raw(name + '.fits') @@ -291,6 +292,54 @@ def test_getitem(self, name, trans): #assert_allclose(c2[1,:,1].value, expected[1,:,1]) assert_allclose(c2[:,1,1].value, expected[:,1,1]) + @pytest.mark.parametrize(('name', 'trans'), translist_vrsc) + def test_getitem_vrsc(self, name, trans): + c, d = cube_and_raw(name + '.fits') + + expected = np.squeeze(d.transpose(trans)) + + # No pv slices for VRSC. + + assert_allclose(c[0,:,:].value, expected[0,:,:]) + + # Not implemented: + #assert_allclose(c[0,0,:].value, expected[0,0,:]) + #assert_allclose(c[0,:,0].value, expected[0,:,0]) + assert_allclose(c[:,0,0].value, expected[:,0,0]) + + assert_allclose(c[1,:,:].value, expected[1,:,:]) + + # Not implemented: + #assert_allclose(c[1,1,:].value, expected[1,1,:]) + #assert_allclose(c[1,:,1].value, expected[1,:,1]) + assert_allclose(c[:,1,1].value, expected[:,1,1]) + + c2 = c.with_spectral_unit(u.km/u.s, velocity_convention='radio') + + assert_allclose(c2[0,:,:].value, expected[0,:,:]) + + # Not implemented: + #assert_allclose(c2[0,0,:].value, expected[0,0,:]) + #assert_allclose(c2[0,:,0].value, expected[0,:,0]) + assert_allclose(c2[:,0,0].value, expected[:,0,0]) + + assert_allclose(c2[1,:,:].value, expected[1,:,:]) + + # Not implemented: + #assert_allclose(c2[1,1,:].value, expected[1,1,:]) + #assert_allclose(c2[1,:,1].value, expected[1,:,1]) + assert_allclose(c2[:,1,1].value, expected[:,1,1]) + + # @pytest.mark.xfail(raises=AttributeError) + @pytest.mark.parametrize(('name', 'trans'), translist_vrsc) + def test_getitem_vrsc(self, name, trans): + c, d = cube_and_raw(name + '.fits') + + expected = np.squeeze(d.transpose(trans)) + + assert_allclose(c[:,:,0].value, expected[:,:,0]) + + class TestArithmetic(object): def setup_method(self, method): @@ -1114,6 +1163,48 @@ def test_beam_attach_to_header(): assert newcube.meta['beam'] == cube.beam +def test_beam_custom(): + + cube, data = cube_and_raw('adv.fits') + + header = cube._header.copy() + beam = Beam.from_fits_header(header) + del header["BMAJ"], header["BMIN"], header["BPA"] + + newcube = SpectralCube(data=data, wcs=cube.wcs, header=header) + + # newcube should now not have a beam + assert not hasattr(newcube, "beam") + + # Attach the beam + newcube = newcube.with_beam(beam=beam) + + assert newcube.beam == cube.beam + + # Header should be updated + assert cube.header["BMAJ"] == newcube.header["BMAJ"] + assert cube.header["BMIN"] == newcube.header["BMIN"] + assert cube.header["BPA"] == newcube.header["BPA"] + + # Should be in meta too + assert newcube.meta['beam'] == cube.beam + + # Try changing the beam properties + newbeam = Beam(beam.major * 2) + + newcube2 = newcube.with_beam(beam=newbeam) + + assert newcube2.beam == newbeam + + # Header should be updated + assert newcube2.header["BMAJ"] == newbeam.major.value + assert newcube2.header["BMIN"] == newbeam.minor.value + assert newcube2.header["BPA"] == newbeam.pa.value + + # Should be in meta too + assert newcube2.meta['beam'] == newbeam + + def test_multibeam_slice(): cube, data = cube_and_raw('vda_beams.fits')