From 3a5ea3568892da0b6502eebc150c3c9e412feeb6 Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Sat, 22 Jun 2013 00:10:02 -0400 Subject: [PATCH 1/5] fixing __all__ in spectrum1d --- specutils/spectrum1d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/specutils/spectrum1d.py b/specutils/spectrum1d.py index 7d8c4dc2c..0035368c9 100644 --- a/specutils/spectrum1d.py +++ b/specutils/spectrum1d.py @@ -3,7 +3,7 @@ from __future__ import print_function, division -__all__ = ['Spectrum1D', 'spec_operation'] +__all__ = ['Spectrum1D'] from astropy import log from astropy.nddata import NDData, FlagCollection From 177f80cbcd0f70e7778a8f14d2da34ad17683371 Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Sat, 22 Jun 2013 14:17:10 -0400 Subject: [PATCH 2/5] added first implementation of generalized spectrum wcs --- specutils/wcs/specwcs.py | 33 ++++++++++++++++++++------------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/specutils/wcs/specwcs.py b/specutils/wcs/specwcs.py index 4f0bf4e94..3820db03e 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 """ @@ -49,24 +48,32 @@ class Spectrum1DLookupWCS(BaseSpectrum1DWCS): lookup table for the array """ - def __init__(self, lookup_table, unit=None): + param_names = ['lookup_table'] + + def __init__(self, lookup_table, unit=None, lookup_table_interpolation_kind='linear'): self.unit = unit - self.lookup_table = lookup_table + self._lookup_table = modeling.parameters.Parameter('lookup_table', 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[0]) != len(np.unique(self.lookup_table[0])): raise BaseSpectrum1DWCSError('The Lookup Table does not describe a unique transformation') - + self.pixel_index = np.arange(len(self.lookup_table[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[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[0], self.pixel_index, left=np.nan, right=np.nan) + else: + raise NotImplementedError('Interpolation type %s is not implemented' % self.lookup_table_interpolation_kind) From 3cc320eb7b2f207aec0029614797004cff5e97ce Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Sat, 22 Jun 2013 15:33:10 -0400 Subject: [PATCH 3/5] tried to implement the interpolation using WCSs --- specutils/spectrum1d.py | 65 ++++++++++----------------------------- specutils/wcs/__init__.py | 2 +- 2 files changed, 18 insertions(+), 49 deletions(-) diff --git a/specutils/spectrum1d.py b/specutils/spectrum1d.py index 0035368c9..4e5783310 100644 --- a/specutils/spectrum1d.py +++ b/specutils/spectrum1d.py @@ -10,7 +10,7 @@ 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) + + 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. 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 From 8e79c90e3c0b4e339506d1bd8dfb57162cf1b16a Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Sat, 22 Jun 2013 16:26:30 -0400 Subject: [PATCH 4/5] finished first implementation of WCS in the specutils class --- specutils/spectrum1d.py | 6 +++--- specutils/wcs/specwcs.py | 28 +++++++++++++++------------- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/specutils/spectrum1d.py b/specutils/spectrum1d.py index 4e5783310..174262916 100644 --- a/specutils/spectrum1d.py +++ b/specutils/spectrum1d.py @@ -228,7 +228,7 @@ def interpolate(self, new_dispersion, kind='linear', bounds_error=True, fill_val new_pixel = self.wcs.invert(new_dispersion.lookup_table) - new_flux = np.interp(new_pixel, self.wcs.pixel_index, self.flux) + 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) @@ -253,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/specwcs.py b/specutils/wcs/specwcs.py index 3820db03e..93c4b55ca 100644 --- a/specutils/wcs/specwcs.py +++ b/specutils/wcs/specwcs.py @@ -31,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) @@ -48,30 +48,30 @@ class Spectrum1DLookupWCS(BaseSpectrum1DWCS): lookup table for the array """ - param_names = ['lookup_table'] + param_names = ['lookup_table_parameter'] def __init__(self, lookup_table, unit=None, lookup_table_interpolation_kind='linear'): self.unit = unit - self._lookup_table = modeling.parameters.Parameter('lookup_table', lookup_table, self, 1) + self._lookup_table_parameter = modeling.parameters.Parameter('lookup_table_parameter', lookup_table, self, 1) 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) #check that array gives a bijective transformation (that forwards and backwards transformations are unique) - if len(self.lookup_table[0]) != len(np.unique(self.lookup_table[0])): + 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[0])) + self.pixel_index = np.arange(len(self.lookup_table_parameter[0])) def __call__(self, pixel_indices): if self.lookup_table_interpolation_kind == 'linear': - return np.interp(pixel_indices, self.pixel_index, self.lookup_table[0], left=np.nan, right=np.nan) + 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): if self.lookup_table_interpolation_kind == 'linear': - return np.interp(dispersion_values, self.lookup_table[0], self.pixel_index, left=np.nan, right=np.nan) + 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) @@ -83,6 +83,7 @@ class Spectrum1DLinearWCS(BaseSpectrum1DWCS): """ + param_names = ['dispersion0', 'dispersion_delta'] @classmethod def from_fits(cls, fname, unit=None, **kwargs): @@ -90,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 From f9fd2d40d5def92bc92136df18595236f662190f Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Sat, 22 Jun 2013 16:37:56 -0400 Subject: [PATCH 5/5] added a little bit of Documentation --- docs/specutils/spectrum1d.rst | 10 ++++++++++ 1 file changed, 10 insertions(+) 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