diff --git a/docs/specutils/spectrum1d.rst b/docs/specutils/spectrum1d.rst index b98f88d87..b97c82257 100644 --- a/docs/specutils/spectrum1d.rst +++ b/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: \ No newline at end of file diff --git a/specutils/spectrum1d.py b/specutils/spectrum1d.py index 7d8c4dc2c..174262916 100644 --- a/specutils/spectrum1d.py +++ b/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) diff --git a/specutils/wcs/__init__.py b/specutils/wcs/__init__.py index c2e09e2e1..bbc385768 100644 --- a/specutils/wcs/__init__.py +++ b/specutils/wcs/__init__.py @@ -1 +1 @@ -from .specwcs import Spectrum1DLinearWCS, Spectrum1DLookupWCS +from .specwcs import Spectrum1DLinearWCS, Spectrum1DLookupWCS, BaseSpectrum1DWCS diff --git a/specutils/wcs/specwcs.py b/specutils/wcs/specwcs.py index 4f0bf4e94..93c4b55ca 100644 --- a/specutils/wcs/specwcs.py +++ b/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,6 +83,7 @@ class Spectrum1DLinearWCS(BaseSpectrum1DWCS): """ + param_names = ['dispersion0', 'dispersion_delta'] @classmethod def from_fits(cls, fname, unit=None, **kwargs): @@ -83,21 +91,22 @@ def from_fits(cls, fname, unit=None, **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