HDU to Projection

CHANGES.rst

@@ -1,5 +1,7 @@
 0.4.1 (unreleased)
 ------------------
+ - Add creating a Projection from a FITS HDU
+   (https://github.com/radio-astro-tools/spectral-cube/pull/376)
  - Deprecate numpy <=1.8 because nanmedian is needed
    (https://github.com/radio-astro-tools/spectral-cube/pull/373)
  - Add tools for masking bad beams in VaryingResolutionSpectralCubes

spectral_cube/cube_utils.py

@@ -15,6 +15,8 @@
 from astropy.io.fits import BinTableHDU, Column
 from astropy import units as u
 import itertools
+import re
+
 
 def _fix_spectral(wcs):
     """
@@ -361,3 +363,32 @@ def _map_context(numcores):
     finally:
         if parallel:
             p.close()
+
+
+def convert_bunit(bunit):
+    '''
+    Convert a BUNIT string to a quantity
+
+    Parameters
+    ----------
+    bunit : str
+        String to convert to an `~astropy.units.Unit`
+
+    Returns
+    -------
+    unit : `~astropy.unit.Unit`
+        Corresponding unit.
+    '''
+
+    # special case: CASA (sometimes) makes non-FITS-compliant jy/beam headers
+    bunit_lower = re.sub("\s", "", bunit.lower())
+    if bunit_lower == 'jy/beam':
+        unit = u.Jy
+    else:
+        try:
+            unit = u.Unit(bunit)
+        except ValueError:
+            warnings.warn("Could not parse unit {0}".format(bunit))
+            unit = None
+
+    return unit

spectral_cube/lower_dimensional_structures.py

@@ -7,6 +7,7 @@
 from .io.core import determine_format
 from . import spectral_axis
 from .utils import SliceWarning
+from .cube_utils import convert_bunit
 
 import numpy as np
 
@@ -185,6 +186,29 @@ def __new__(cls, value, unit=None, dtype=None, copy=True, wcs=None,
 
         return self
 
+    @staticmethod
+    def from_hdu(hdu):
+        '''
+        Return a projection from a FITS HDU.
+        '''
+
+        if not len(hdu.shape) == 2:
+            raise ValueError("HDU must contain two-dimensional data.")
+
+        meta = {}
+
+        mywcs = wcs.WCS(hdu.header)
+
+        if "BUNIT" in hdu.header:
+            unit = convert_bunit(hdu.header["BUNIT"])
+            meta["BUNIT"] = hdu.header["BUNIT"]
+        else:
+            unit = None
+
+        self = Projection(hdu.data, unit=unit, wcs=mywcs, meta=meta,
+                          header=hdu.header)
+
+        return self
 
     def quicklook(self, filename=None, use_aplpy=True, aplpy_kwargs={}):
         """

spectral_cube/spectral_cube.py

@@ -109,19 +109,7 @@ def __init__(self, data, wcs, mask=None, meta=None, fill_value=np.nan,
             raise TypeError("If a header is given, it must be a fits.Header")
 
         if 'BUNIT' in self._meta:
-
-            # special case: CASA (sometimes) makes non-FITS-compliant jy/beam headers
-            bunit = re.sub("\s", "", self._meta['BUNIT'].lower())
-            if bunit == 'jy/beam':
-                self._unit = u.Jy
-
-
-            else:
-                try:
-                    self._unit = u.Unit(self._meta['BUNIT'])
-                except ValueError:
-                    warnings.warn("Could not parse unit {0}".format(self._meta['BUNIT']))
-                    self._unit = None
+            self._unit = cube_utils.convert_bunit(self._meta["BUNIT"])
         elif hasattr(data, 'unit'):
             self._unit = data.unit
         else:
@@ -2927,13 +2915,13 @@ def identify_bad_beams(self, threshold, reference_beam=None,
         for prop in criteria:
             val = props[prop]
             mid = getattr(reference_beam, prop)
-            
+
             diff = np.abs((val-mid)/mid)
 
             assert diff.shape == includemask.shape
-            
+
             includemask[diff > threshold] = False
-        
+
         return includemask
 
     def mask_out_bad_beams(self, threshold, reference_beam=None,

spectral_cube/tests/test_projection.py

@@ -193,3 +193,16 @@ def test_quantity_property():
 
     assert isinstance(arr, u.Quantity)
     assert not isinstance(arr, OneDSpectrum)
+
+
+@pytest.mark.parametrize(('LDO', 'data'),
+                         zip(LDOs_2d, data_two_2d))
+def test_projection_from_hdu(LDO, data):
+
+    p = LDO(data, copy=False)
+
+    hdu = p.hdu
+
+    p_new = LDO.from_hdu(hdu)
+
+    assert (p == p_new).all()