First go at spectral WCS

docs/specutils/spectrum1d.rst

@@ -31,5 +31,15 @@ Once a spectrum is instantiated, one can access the `flux` and `dispersion`::
     array([ 0.34272852,  0.92834782,  0.64680224, ...,  0.03348069,
             0.10291822,  0.33614334])
 
+**Playing with WCS**:
+
+    >>> flux = np.random.random(1000)
+    >>> wave = np.sort(np.random.random(1000))
+    >>> spec1d = Spectrum1D.from_array(wave, flux)
+    >>> spec1d = Spectrum1D.interpolate(Spectrum1DLinearWCS(6000, dispersion_delta=1, pixel_index=np.arange(1000)))
+    >>> #Now it can be written to a fits file as it is linear encoded which is representable in FITS headers.
+
+
+
 .. automodapi:: specutils
     :no-inheritance-diagram:

specutils/spectrum1d.py

@@ -3,14 +3,14 @@
 
 from __future__ import print_function, division
 
-__all__ = ['Spectrum1D', 'spec_operation']
+__all__ = ['Spectrum1D']
 
 from astropy import log
 from astropy.nddata import NDData, FlagCollection
 
 from astropy.utils import misc
 
-from .wcs import Spectrum1DLookupWCS, Spectrum1DLinearWCS
+from specutils.wcs import BaseSpectrum1DWCS, Spectrum1DLookupWCS
 
 import numpy as np
 
@@ -183,7 +183,7 @@ def interpolate(self, new_dispersion, kind='linear', bounds_error=True, fill_val
 
         Parameters
         ----------
-        new_disp : `~numpy.ndarray`
+        new_dispersion : `~numpy.ndarray`
             The dispersion array to interpolate the flux on to.
 
         kind : `str` or `int`, optional
@@ -217,53 +217,22 @@ def interpolate(self, new_dispersion, kind='linear', bounds_error=True, fill_val
         """
 
         # Check for SciPy availability
-        try:
-            from scipy import interpolate
-        except ImportError as e:
-            raise ImportError("Could not import interpolate from scipy; cannot"+
-                              " interpolate to new dispersion map without this"+
-                              " (need scipy.interpolate.interp1d)")
-
-        spectrum_interp = interpolate.interp1d(self.dispersion,
-                                               self.flux,
-                                               kind=kind,
-                                               bounds_error=bounds_error,
-                                               fill_value=fill_value)
-
-        new_flux = spectrum_interp(new_dispersion)
-
-        # We need to perform error calculation for the new dispersion map
-        if self.uncertainty is None:
-            new_uncertainty = None
-        else:
-            # After having a short think about it, it seems reasonable to me only to
-            # take the nearest uncertainty for each interpolated dispersion point
-
-            new_uncertainty = interpolate.interp1d(self.dispersion,
-                                             self.flux,
-                                             kind=1, # Nearest
-                                             bounds_error=bounds_error,
-                                             fill_value=fill_value)
-
-        # The same should also apply for masks
-        if self.mask is None:
-            new_mask = None
-        else:
-            new_mask = interpolate.interp1d(self.dispersion,
-                                            self.flux,
-                                            kind=1, # Nearest
-                                            bounds_error=bounds_error,
-                                            fill_value=fill_value)
-
-        # As for flags it is not entirely clear to me what the best behaviour is
-        # In the face of uncertainty, for the time being, I am discarding flags
-
-        return self.__class__.from_array(new_dispersion,
-                                         new_flux,
-                                         uncertainty=new_uncertainty,
-                                         mask=new_mask,
-                                         meta=self.meta)
+
+
+        if kind != 'linear':
+            raise ValueError('No other kind but linear supported')
+
+        if not isinstance(new_dispersion, BaseSpectrum1DWCS):
+            new_dispersion = Spectrum1DLookupWCS(np.array(new_dispersion))
+
+
+        new_pixel = self.wcs.invert(new_dispersion.lookup_table)
+
+        new_flux = np.interp(new_pixel, self.wcs.pixel_index, self.flux, left=np.nan, right=np.nan)
 
+
+        return self.__class__(new_flux, wcs=new_dispersion, meta=self.meta)
+
 
     def slice_dispersion(self, start=None, stop=None):
         """Slice the spectrum within a given start and end dispersion value.
@@ -284,8 +253,8 @@ def slice_dispersion(self, start=None, stop=None):
 
         For example::
 
-            >>> from astropy.specutils import Spectrum1D
-            >>> from astropy.units import Units as unit
+            >>> from specutils import Spectrum1D
+            >>> from astropy.units import Unit as unit
             >>> import numpy as np
 
             >>> dispersion = np.arange(4000, 5000, 0.12)

specutils/wcs/__init__.py

@@ -1 +1 @@
-from .specwcs import Spectrum1DLinearWCS, Spectrum1DLookupWCS
+from .specwcs import Spectrum1DLinearWCS, Spectrum1DLookupWCS, BaseSpectrum1DWCS

specutils/wcs/specwcs.py

@@ -8,14 +8,13 @@
 
 from astropy.utils import misc
 from astropy.io import fits
+from astropy import modeling
 
-class NDenModelsPlaceHolder(object):
 
-    pass
 class BaseSpectrum1DWCSError(Exception):
     pass
 
-class BaseSpectrum1DWCS(NDenModelsPlaceHolder):
+class BaseSpectrum1DWCS(modeling.Model):
     """
         Base class for a Spectrum1D WCS
     """
@@ -32,9 +31,9 @@ def dispersion2pixel(self, dispersion_value):
         """
         return self.invert(dispersion_value)
 
-    def create_lookup_table(self, pixel_indices):
-        self.lookup_table = self(pixel_indices)
-
+    @misc.lazyproperty
+    def lookup_table(self):
+        return self(self.pixel_index)
 
 
 
@@ -49,24 +48,32 @@ class Spectrum1DLookupWCS(BaseSpectrum1DWCS):
         lookup table for the array
     """
 
-    def __init__(self, lookup_table, unit=None):
+    param_names = ['lookup_table_parameter']
+
+    def __init__(self, lookup_table, unit=None, lookup_table_interpolation_kind='linear'):
         self.unit = unit
-        self.lookup_table = lookup_table
+        self._lookup_table_parameter = modeling.parameters.Parameter('lookup_table_parameter', lookup_table, self, 1)
 
-        #check that array gives a bijective transformation (that forwards and backwards transformations are unique)
+        self.lookup_table_interpolation_kind = lookup_table_interpolation_kind
+        super(Spectrum1DLookupWCS, self).__init__(self.param_names, n_inputs=1, n_outputs=1, param_dim=1)
 
-        if len(self.lookup_table) != len(np.unique(self.lookup_table)):
+        #check that array gives a bijective transformation (that forwards and backwards transformations are unique)
+        if len(self.lookup_table_parameter[0]) != len(np.unique(self.lookup_table_parameter[0])):
             raise BaseSpectrum1DWCSError('The Lookup Table does not describe a unique transformation')
-
+        self.pixel_index = np.arange(len(self.lookup_table_parameter[0]))
 
     def __call__(self, pixel_indices):
-        if misc.isiterable(pixel_indices) and not isinstance(pixel_indices, basestring):
-            pixel_indices = np.array(pixel_indices)
-        return self.lookup_table[pixel_indices]
+        if self.lookup_table_interpolation_kind == 'linear':
+            return np.interp(pixel_indices, self.pixel_index, self.lookup_table_parameter[0], left=np.nan, right=np.nan)
+        else:
+            raise NotImplementedError('Interpolation type %s is not implemented' % self.lookup_table_interpolation_kind)
 
-    def invert(self, dispersion_values):
 
-        return np.searchsorted(self.lookup_table, dispersion_values)
+    def invert(self, dispersion_values):
+        if self.lookup_table_interpolation_kind == 'linear':
+            return np.interp(dispersion_values, self.lookup_table_parameter[0], self.pixel_index, left=np.nan, right=np.nan)
+        else:
+            raise NotImplementedError('Interpolation type %s is not implemented' % self.lookup_table_interpolation_kind)
 
 
 
@@ -76,28 +83,30 @@ class Spectrum1DLinearWCS(BaseSpectrum1DWCS):
 
     """
 
+    param_names = ['dispersion0', 'dispersion_delta']
 
     @classmethod
     def from_fits(cls, fname, unit=None, **kwargs):
         header = fits.getheader(fname, **kwargs)
         return cls(header['CRVAL1'], header['CDELT1'], header['CRPIX1'] - 1, unit=unit)
 
 
-    def __init__(self, dispersion0, dispersion_delta, dispersion_pixel0=0, unit=None):
+    def __init__(self, dispersion0, dispersion_delta, pixel_index, unit=None):
         self.unit = unit
+        self.pixel_index = pixel_index
         self.dispersion0 = dispersion0
         self.dispersion_delta = dispersion_delta
-        self.dispersion_pixel0 = dispersion_pixel0
+
 
     def __call__(self, pixel_indices):
         if misc.isiterable(pixel_indices) and not isinstance(pixel_indices, basestring):
             pixel_indices = np.array(pixel_indices)
-        return self.dispersion0 + self.dispersion_delta * (pixel_indices - self.dispersion_pixel0)
+        return self.dispersion0 + self.dispersion_delta * (pixel_indices - self.pixel_index[0])
 
     def invert(self, dispersion_values):
         if misc.isiterable(dispersion_values) and not isinstance(dispersion_values, basestring):
             dispersion_values = np.array(dispersion_values)
-        return (dispersion_values - self.dispersion0) / self.dispersion_delta + self.dispersion_pixel0
+        return (dispersion_values - self.dispersion0) / self.dispersion_delta + self.pixel_index[0]
 
 
 #### EXAMPLE implementation for Chebyshev