Merge pull request #506 from keflavich/several_fixes

Several serious bug fixes

.travis.yml

@@ -43,9 +43,13 @@ matrix:
 
         # Do a coverage test in Python 3.
         # make sure *everything* is run for coverage tests
-        - python: 3.6
-          env: SETUP_CMD='test --coverage'
-               CONDA_DEPENDENCIES='matplotlib aplpy yt bottleneck pytest pytest-xdist astropy-helpers joblib'
+        ## COVERAGE TESTS DISABLED: they don't work any more.  I have found no combination of anything that fixes this.
+        #- python: 3.6 # (travis doesn't support py3.7 natively) use py3.7 because it may not have this bug? https://github.com/pytest-dev/pytest/issues/2854
+        #  env: SETUP_CMD='test --coverage'
+        #       PYTHON_VERSION=3.7
+        #       CONDA_DEPENDENCIES='matplotlib aplpy yt bottleneck pytest pytest-xdist astropy-helpers joblib'
+        #       # pytest 3.9.3 suggested by @e-koch, @astrofrog.  3.10.1 doesn't work, 4.x doesn't work
+        #       # overrode python version too, per bsipocz
 
         # Check for sphinx doc build warnings - we do this first because it
         # may run for a long time

docs/index.rst

@@ -69,6 +69,7 @@ Getting started
    spectral_extraction.rst
    continuum_subtraction.rst
    examples.rst
+   visualization.rst
 
 Advanced
 ^^^^^^^^

docs/visualization.rst

@@ -0,0 +1,51 @@
+Visualization
+=============
+
+Spectral-cube is not primarily a visualization package, but it has several
+tools for visualizing subsets of the data.
+
+All lower-dimensional subsets,
+`~spectral_cube.lower_dimensional_structures.OneDSpectrum`, and
+`~spectral_cube.lower_dimensional_structures.Projection`, have their own
+``quicklook`` methods
+(`~spectral_cube.lower_dimensional_structures.OneDSpectrum.quicklook` and
+`~spectral_cube.lower_dimensional_structures.Projection.quicklook`,
+respectively).  These methods will plot the data with somewhat properly labeled
+axes.
+
+The two-dimensional viewers default to using `aplpy <http://aplpy.github.io/>`_.
+Because of quirks of how aplpy sets up its plotting window, these methods will
+create their own figures.  If ``use_aplpy`` is set to ``False``, and similarly
+if you use the ``OneDSpectrum`` quicklook, the data will be overplotted in the
+latest used plot window.
+
+
+In principle, one can also simply plot the data.  For example, if you have a cube,
+you could do::
+
+    >>> plt.plot(cube[:,0,0]) # doctest: +SKIP
+
+to plot a spectrum sliced out of the cube or::
+
+    >>> plt.imshow(cube[0,:,:]) # doctest: +SKIP
+
+to plot an image slice. 
+
+.. warning:: 
+
+   There are known incompatibilities with the above plotting approach:
+   matplotlib versions ``<2.1`` will crash, and you will have to clear the plot
+   window to reset it.
+
+
+Other Visualization Tools
+=========================
+To visualize the cubes directly, you can use some of the other tools we provide
+for pushing cube data into external viewers.
+
+See :doc:`yt_example` for using yt as a visualization tool.
+
+
+The `spectral_cube.SpectralCube.to_glue` and
+`spectral_cube.SpectralCube.to_ds9` methods will send the whole cube to glue
+and ds9.  This approach generally requires loading the whole cube into memory.

spectral_cube/lower_dimensional_structures.py

@@ -157,13 +157,14 @@ def __getitem__(self, key, **kwargs):
         return new
 
     def __array_finalize__(self, obj):
-        #log.debug("Finalizing self={0}{1} obj={2}{3}"
-        #          .format(self, type(self), obj, type(obj)))
         self._wcs = getattr(obj, '_wcs', None)
         self._meta = getattr(obj, '_meta', None)
         self._mask = getattr(obj, '_mask', None)
         self._header = getattr(obj, '_header', None)
         self._spectral_unit = getattr(obj, '_spectral_unit', None)
+        self._fill_value = getattr(obj, '_fill_value', np.nan)
+        self._wcs_tolerance = getattr(obj, '_wcs_tolerance', 0.0)
+
         super(LowerDimensionalObject, self).__array_finalize__(obj)
 
     @property
@@ -178,6 +179,12 @@ def array(self):
         """
         return np.asarray(self)
 
+    @property
+    def _data(self):
+        # the _data property is required by several other mixins
+        # (which probably means defining it here is a bad design)
+        return self.array
+
     @property
     def quantity(self):
         """
@@ -278,11 +285,40 @@ def shrink_mask(self):
             self._mask = nomask
         return self
 
-class Projection(LowerDimensionalObject, SpatialCoordMixinClass):
+    def _initial_set_mask(self, mask):
+        """
+        Helper tool to validate mask when originally setting it in __new__
+
+        Note that because this is intended to be used in __new__, order
+        matters: ``self`` must have ``_wcs``, for example.
+        """
+        if mask is None:
+            mask = BooleanArrayMask(np.ones_like(self.value, dtype=bool),
+                                    self._wcs, shape=self.value.shape)
+        elif isinstance(mask, np.ndarray):
+            if mask.shape != self.value.shape:
+                raise ValueError("Mask shape must match the {0} shape."
+                                 .format(self.__class__.__name__)
+                                )
+            mask = BooleanArrayMask(mask, self._wcs, shape=self.value.shape)
+        elif isinstance(mask, MaskBase):
+            pass
+        else:
+            raise TypeError("mask of type {} is not a supported mask "
+                            "type.".format(type(mask)))
+
+        # Validate the mask before setting
+        mask._validate_wcs(new_data=self.value, new_wcs=self._wcs,
+                           wcs_tolerance=self._wcs_tolerance)
+
+        self._mask = mask
+
+class Projection(LowerDimensionalObject, SpatialCoordMixinClass,
+                 MaskableArrayMixinClass):
 
     def __new__(cls, value, unit=None, dtype=None, copy=True, wcs=None,
                 meta=None, mask=None, header=None, beam=None,
-                read_beam=False):
+                fill_value=np.nan, read_beam=False, wcs_tolerance=0.0):
 
         if np.asarray(value).ndim != 2:
             raise ValueError("value should be a 2-d array")
@@ -294,7 +330,11 @@ def __new__(cls, value, unit=None, dtype=None, copy=True, wcs=None,
                                   copy=copy).view(cls)
         self._wcs = wcs
         self._meta = {} if meta is None else meta
-        self._mask = mask
+        self._wcs_tolerance = wcs_tolerance
+
+        self._initial_set_mask(mask)
+
+        self._fill_value = fill_value
         if header is not None:
             self._header = header
         else:
@@ -338,6 +378,65 @@ def with_beam(self, beam):
 
         return self
 
+    def with_fill_value(self, fill_value):
+        """
+        Create a new :class:`Projection` or :class:`Slice` with a different
+        ``fill_value``.
+        """
+        return self._new_projection_with(fill_value=fill_value)
+
+    def _new_projection_with(self, data=None, wcs=None, mask=None, meta=None,
+                             fill_value=None, spectral_unit=None, unit=None,
+                             header=None, wcs_tolerance=None, beam=None,
+                             **kwargs):
+
+        data = self._data if data is None else data
+        if unit is None and hasattr(data, 'unit'):
+            if data.unit != self.unit:
+                raise u.UnitsError("New data unit '{0}' does not"
+                                   " match unit '{1}'.  You can"
+                                   " override this by specifying the"
+                                   " `unit` keyword."
+                                   .format(data.unit, self.unit))
+            unit = data.unit
+        elif unit is None:
+            unit = self.unit
+        elif unit is not None:
+            # convert string units to Units
+            if not isinstance(unit, u.Unit):
+                unit = u.Unit(unit)
+
+            if hasattr(data, 'unit'):
+                if u.Unit(unit) != data.unit:
+                    raise u.UnitsError("The specified new cube unit '{0}' "
+                                       "does not match the input unit '{1}'."
+                                       .format(unit, data.unit))
+            else:
+                data = u.Quantity(data, unit=unit, copy=False)
+
+        wcs = self._wcs if wcs is None else wcs
+        mask = self._mask if mask is None else mask
+        if meta is None:
+            meta = {}
+            meta.update(self._meta)
+        if unit is not None:
+            meta['BUNIT'] = unit.to_string(format='FITS')
+
+        fill_value = self._fill_value if fill_value is None else fill_value
+
+        if beam is None:
+            if hasattr(self, 'beam'):
+                beam = self.beam
+
+        newproj = self.__class__(value=data, wcs=wcs, mask=mask, meta=meta,
+                                 unit=unit, fill_value=fill_value,
+                                 header=header or self._header,
+                                 wcs_tolerance=wcs_tolerance or self._wcs_tolerance,
+                                 beams=beam,
+                                 **kwargs)
+
+        return newproj
+
     @property
     def beam(self):
         return self._beam
@@ -595,24 +694,8 @@ def __new__(cls, value, unit=None, dtype=None, copy=True, wcs=None,
         self._meta = {} if meta is None else meta
         self._wcs_tolerance = wcs_tolerance
 
-        if mask is None:
-            mask = BooleanArrayMask(np.ones_like(self.value, dtype=bool),
-                                    self._wcs, shape=self.value.shape)
-        elif isinstance(mask, np.ndarray):
-            if mask.shape != self.value.shape:
-                raise ValueError("Mask shape must match the spectrum shape.")
-            mask = BooleanArrayMask(mask, self._wcs, shape=self.value.shape)
-        elif isinstance(mask, MaskBase):
-            pass
-        else:
-            raise TypeError("mask of type {} is not a supported mask "
-                            "type.".format(type(mask)))
-
-        # Validate the mask before setting
-        mask._validate_wcs(new_data=self.value, new_wcs=self._wcs,
-                           wcs_tolerance=self._wcs_tolerance)
+        self._initial_set_mask(mask)
 
-        self._mask = mask
         self._fill_value = fill_value
         if header is not None:
             self._header = header
@@ -630,10 +713,6 @@ def __new__(cls, value, unit=None, dtype=None, copy=True, wcs=None,
         if beams is not None:
             self.beams = beams
 
-        # HACK: OneDSpectrum should eventually become not-a-quantity
-        # Maybe it should be a u.Quantity(np.ma)?
-        self._data = self.value
-
         return self
 
     def __repr__(self):
@@ -1013,14 +1092,5 @@ def __getattribute__(self, attrname):
         else:
             return super(OneDSpectrum, self).__getattribute__(attrname)
 
-    def __array_finalize__(self, obj):
-        #from astropy import log
-        #log.debug("in OneDSpectrum, Finalizing self={0}{1} obj={2}{3}"
-        #          .format(self, type(self), obj, type(obj)))
-        self._fill_value = getattr(obj, '_fill_value', np.nan)
-        self._data = self.view(np.ndarray)
-        self._wcs_tolerance = getattr(obj, '_wcs_tolerance', 0.0)
-        super(OneDSpectrum, self).__array_finalize__(obj)
-
 class VaryingResolutionOneDSpectrum(OneDSpectrum, MultiBeamMixinClass):
     pass

spectral_cube/masks.py

@@ -1,6 +1,7 @@
 from __future__ import print_function, absolute_import, division
 
 import abc
+import warnings
 
 import numpy as np
 from numpy.lib.stride_tricks import as_strided
@@ -10,6 +11,7 @@
 from astropy.extern.six.moves import zip
 
 from . import wcs_utils
+from .utils import WCSWarning
 
 
 __all__ = ['MaskBase', 'InvertedMask', 'CompositeMask', 'BooleanArrayMask',
@@ -112,6 +114,22 @@ def include(self, data=None, wcs=None, view=(), **kwargs):
         self._validate_wcs(data, wcs, **kwargs)
         return self._include(data=data, wcs=wcs, view=view)
 
+    # Commented out, but left as a possibility, because including this didn't fix any
+    # of the problems we encountered with matplotlib plotting
+    def view(self, view=()):
+        """
+        Compatibility tool: if a numpy.ma.ufunc is run on the mask, it will try
+        to grab a view of the mask, which needs to appear to numpy as a true
+        array.  This can be important for, e.g., plotting.
+
+        Numpy's convention is that masked=True means "masked out"
+
+        .. note::
+            I don't know if there are broader concerns or consequences from
+            including this 'view' tool here.
+        """
+        return self.exclude(view=view)
+
     def _validate_wcs(self, new_data=None, new_wcs=None, **kwargs):
         """
         This method can be overridden in cases where the data and WCS have to
@@ -485,11 +503,19 @@ def _validate_wcs(self, new_data=None, new_wcs=None, **kwargs):
             raise ValueError("data shape cannot be broadcast to match mask shape")
         if new_wcs is not None:
             if new_wcs not in self._wcs_whitelist:
-                if not wcs_utils.check_equality(new_wcs, self._wcs,
-                                                warn_missing=True,
-                                                **kwargs):
-                    raise ValueError("WCS does not match mask WCS")
-                else:
+                try:
+                    if not wcs_utils.check_equality(new_wcs, self._wcs,
+                                                    warn_missing=True,
+                                                    **kwargs):
+                        raise ValueError("WCS does not match mask WCS")
+                    else:
+                        self._wcs_whitelist.add(new_wcs)
+                except InconsistentAxisTypesError:
+                    warnings.warn("Inconsistent axis type encountered; WCS is "
+                                  "invalid and therefore will not be checked "
+                                  "against other WCSes.",
+                                  WCSWarning
+                                  )
                     self._wcs_whitelist.add(new_wcs)
 
     def _include(self, data=None, wcs=None, view=()):