Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

Spec WCS #22

Merged
merged 11 commits into from

4 participants

@wkerzendorf
Owner

Hey guys,
I finally did my spec wcs implementation. This also solved the problem that @kbarbary and @keflavich had with the missing dispersion_unit (a I think we should switch all disp to dispersion as that was discussed and agreed upon - next PR).

Here's a little taster of what it can do:

from astropy.table import Table
from astropy import units as u
from specutils import Spectrum1D
from specutils.wcs import Spectrum1DLinearWCS

wavelength = np.arange(6000, 7000, 0.1)
flux = np.random.random(10000)
uncertainty = np.sqrt(flux)

mydata = Table(dict(disp=wavelength, flux=flux, uncertainty=uncertainty))
mydata['disp'].units = u.angstrom
mydata['flux'].units = u.Unit('W/(m^2 angstrom^2)')
#initializing from_table
myspec = Spectrum1D.from_table(mydata, uncertainty_column='uncertainty')

print myspec.disp
[ 6000.   6000.1  6000.2 ...,  6999.7  6999.8  6999.9]


print myspec.flux
[ 0.70864274  0.54175854  0.69471822 ...,  0.13276595  0.58555977
  0.88461505]


print myspec.dispersion_unit
Angstrom
print myspec.flux_unit
W / (Angstrom2 m2)

#initializing with a linear wcs
linear_wcs = Spectrum1DLinearWCS(dispersion0=6000, dispersion_delta=0.1, dispersion_pixel0=0, unit=u.angstrom)
myspec2 = Spectrum1D(data=flux, wcs=linear_wcs, unit=u.Unit('W / (Angstrom2 m2)'))


print myspec2.disp
[ 6000.   6000.1  6000.2 ...,  6999.7  6999.8  6999.9]

print myspec2.wcs.pixel2dispersion([0, 9999])
[ 6000.  7000.]

print myspec2.wcs.dispersion2pixel([6000, 6999.9])
[    0.  9999.]

What do you guys think?

@keflavich keflavich commented on the diff
specutils/wcs/specwcs.py
((12 lines not shown))
+class NDenModelsPlaceHolder(object):
+
+ pass
+class BaseSpectrum1DWCSError(Exception):
+ pass
+
+class BaseSpectrum1DWCS(NDenModelsPlaceHolder):
+ """
+ Base class for a Spectrum1D WCS
+ """
+
+ def pixel2dispersion(self, pixel_index):
+ """
+ This should be the forward transformation in normal WCS classes
+ """
+ return self(pixel_index)
@keflavich Owner

self[pixel_index] ? Will NDenModelsPlaceHolder have a __call__ method or a __getitem__ method? I think flushing out the inherited class would help a little.

@wkerzendorf Owner

It does - I can just copy it accross from her PR. I just didn't want to do that so we don't have that in the commits. To give you an overview: https://github.com/nden/astropy/blob/master/astropy/models/models.py. I think it's days away from being merged.

@keflavich Owner

Got it, that's good enough. Thanks.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@keflavich keflavich commented on the diff
specutils/wcs/specwcs.py
((44 lines not shown))
+
+ Parameters
+ ----------
+
+ lookup_table : ~np.ndarray
+ lookup table for the array
+ """
+
+ def __init__(self, lookup_table, unit=None):
+ self.unit = unit
+ self.lookup_table = lookup_table
+
+ #check that array gives a bijective transformation (that forwards and backwards transformations are unique)
+
+ if len(self.lookup_table) != len(np.unique(self.lookup_table)):
+ raise BaseSpectrum1DWCSError('The Lookup Table does not describe a unique transformation')
@keflavich Owner

I like this, but I think it means that you can't use Spectrum1DLookupWCS with Spectrum1DEchelle as @hamogu has suggested, since it will require 2D-dispersion. That's probably OK, I'm just noting it.

@wkerzendorf Owner

So imho I think I want to focus on a Spectrum1D first. If you open up two dimension it opens up a whole can of worms (IFUs, Echelle, spectral-spatial spectra, ..... ).

@keflavich Owner

Yes, that's fine. Just something to keep in mind when we extend to >1D

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@keflavich
Owner

Overall, I like this, but keeping in mind the discussion from #18, I think performance costs will probably lead to a lookup table being preferable. But I think maybe that's what you're suggesting? I'm not entirely sure the code example matches the committed code... the last few examples suggest that, upon creation of a linear WCS, a lookup_table will be created the first time it is needed. Is that right?

If so, it might be good to make these:

        self.dispersion0 = dispersion0
        self.dispersion_delta = dispersion_delta
        self.dispersion_pixel0 = dispersion_pixel0

properties that will update the lookup_table appropriately if they are modified

@hamogu
Collaborator

I am not an optical spectroscopist, but for now I am not quite sure what the use case for this construction is. myspec contains all the information I need to plot, fit, interpolate, add, etc.
myspec.dispersion is a numpy array (or an ndarray) that is plain, simple, fast and easy.

As far as I understand, all that the trouble with expressing the dispersion by a WCS system (linear or polynomial) is going to buy us is that we save a few bytes on disk when we write fits files (because we don't need to save a wavelength array, but instead we put a few coefficients in the header).
Thus, I would do all this in a file writer.
I may just be ignorant (if so, please let me know and I'll shut up), but why do we want to keep using the WCS while the spectrum is in momory?

@hamogu hamogu commented on the diff
specutils/wcs/specwcs.py
@@ -0,0 +1,112 @@
+# This module provides a basic and probably temporary WCS solution until astropy has a wcs built-in
+#This is all built upon @nden's work on the models class
+
@hamogu Collaborator
hamogu added a note

I appreciate the work you put in this. But do we really need a temporary solution before nede's work is merged into astropy? It's probably going to happen very soon.

@keflavich Owner

Will the framework change much with incorporation of nden's work? I think we still need much of this, but I really don't know much about the models framework yet.

@wkerzendorf Owner

@hamogu. Well it's not quite a temporary solution.I was involved in the API shaping of some of the wcs part of @nden 's model framework so I know that the things I have written here are not really duplicating the @nden's work. I think not much will change in this specwcs - but will be greatly enriched by actually inheriting from the astropy-model. One example is that we will get an easy way to fit these models - so we can ditch IRAF for the wavelength calibration part pretty much right away ;-)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@keflavich
Owner

I think there's no harm in tracking the WCS if it is linear or easily representable as a polynomial.

One major use case is for data cubes, in which you'd want to express the WCS axes as a few bytes in the header rather than another cube of the same size storing the WCS information for each pixel. I have plenty of ~GB size data cubes (and ALMA/EVLA will make LOTS of these) that I would very much like to keep as small as possible.

Where it's not possible to simplify the WCS axis with a few coefficients, we should use a look-up table, but if the overhead is small and the technical challenge not too great (which I think @wkerzendorf has shown adequately), we should try to track the dispersion axis as a series of coefficients.

As you say, it may generally be best to do this kind of work in a file writer rather than dealing with it in memory, but again the data cube case comes to mind: it may not even be possible to load a full cube's worth of dispersion axes into memory. But that can perhaps be dealt with in other ways.

@hamogu hamogu commented on the diff
specutils/wcs/specwcs.py
((21 lines not shown))
+ """
+
+ def pixel2dispersion(self, pixel_index):
+ """
+ This should be the forward transformation in normal WCS classes
+ """
+ return self(pixel_index)
+
+ def dispersion2pixel(self, dispersion_value):
+ """
+ This should be the inverse transformation in normal WCS classes
+ """
+ return self.invert(dispersion_value)
+
+ def create_lookup_table(self, pixel_indices):
+ self.lookup_table = self(pixel_indices)
@hamogu Collaborator
hamogu added a note

Should we default that to use all pixels?

def create_lookup_table(self, pixel_indices = None):
if pixel_indices is None:
    pixel_indixes = np.arange(self.dispersion0, self.dispersion_delta * len(self.flux))
@keflavich Owner

+1: I think that is the intent illustrated in the example at the top of this PR. Though I think it should be

def create_lookup_table(self, pixel_indices = None):
    if pixel_indices is None:
        pixel_indices = np.arange(self.npix) + 1
    self.lookup_table = self(pixel_indices)

?

@wkerzendorf Owner

So two comments to that. I thought about adding an npix or equivalent to the WCS class - but decided against it because it "felt" cleaner to keep NAXIS1 and WCS separate. Maybe an interesting compromise might be to tell wcs about its parent and thus it could access this information.

The second design decision was to give the first pixel the index 0. I know, that IRAF doesn't do it that way (that's why CRPIX1=1 for the first pixel).

@keflavich Owner

I agree with giving the first pixel index 0.... does that mean when you load in CRPIX from a FITS header, you will do:
self.crpix = header.get('CRPIX')-1 ?

@wkerzendorf Owner

@keflavich precisely

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@wkerzendorf
Owner

(@keflavich) Overall, I like this, but keeping in mind the discussion from #18, I think performance costs will probably lead to a lookup table being preferable. But I think maybe that's what you're suggesting? I'm not entirely sure the code example matches the committed code... the last few examples suggest that, upon creation of a linear WCS, a lookup_table will be created the first time it is needed. Is that right?

That is right.

(@hamogu) As far as I understand, all that the trouble with expressing the dispersion by a WCS system (linear or polynomial) is going to buy us is that we save a few bytes on disk when we write fits files (because we don't need to save a wavelength array, but instead we put a few coefficients in the header).

There are multiple reasons for keeping track of the actual WCS:
The few bytes on disk don't matter. The main problem is that I want to use my astropy with other pieces of code that read in FITS files (e.g. IRAF) they need some WCS, which we can only easily write if we know coefficients. As an optical spectroscopist ;-) I had to do that multiple times already (going back and forth between different codes).

In addition, interpolation needs to happen from the WCS side, so it can make right sub-pixel adjustments - that can't happen if all information that you have is lookup table.

(@hamogu) I may just be ignorant (if so, please let me know and I'll shut up), but why do we want to keep using the WCS while the spectrum is in memory?

That's why I create a lookup table the first time you want one. myspec.dispersion does that and keeps that one in memory.

@keflavich
Owner

OK, so I certainly get your intent (make a WCS, keep track of the coefficients, make a lookup table), but the code doesn't do that all quite yet. What is needed before it's fully functional? Will this "just work" once models is merged?

@wkerzendorf
Owner

@keflavich So it works just now. The example I gave is what the code does (I ran it and put in the printouts).

@keflavich
Owner

@wkerzendorf Ahhh, I see. I had missed some things because git hid them and others because I didn't realize how much was inherited from nddata. Cool!

@hamogu
Collaborator

@wkerzendorf that convinces me!
Cool + 1!

@wkerzendorf
Owner

@hamogu, @keflavich : Great. Well I would like to merge this if there are not further suggestions. talking about the API and the immediate future: I would love to have some of @hamogu's fits wcs reader making wcs's and opening spectra. If that works, I think we should write the API suggestion and showcase it.

When talking about WCS - I think we are at the perfect position to pilot that as we only have to deal with one dimension and have an easier way of trying things out.

Thoughts?

specutils/__init__.py
@@ -6,6 +6,7 @@
package, but for now is being developed independently.
"""
+from spectrum1d import Spectrum1D
@astrofrog Owner

This should be a relative import

@wkerzendorf Owner

I see - great I'll fix that just now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
specutils/wcs/__init__.py
@@ -0,0 +1 @@
+from specwcs import Spectrum1DLinearWCS, Spectrum1DLookupWCS
@astrofrog Owner

This should be a relative import

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@wkerzendorf
Owner

@astrofrog I fixed the relative import and another small fix. So looking through this, what do you think?

@keflavich
Owner
@wkerzendorf wkerzendorf merged commit 362ef5f into astropy:master
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
View
3  specutils/__init__.py
@@ -6,6 +6,7 @@
package, but for now is being developed independently.
"""
+from .spectrum1d import Spectrum1D
#this indicates whether or not we are in the package's setup.py
try:
@@ -128,3 +129,5 @@ def test(package=None, test_path=None, args=None, plugins=None,
del e
del os, warn, config_dir # clean up namespace
+
+
View
41 specutils/spectrum1d.py
@@ -8,6 +8,8 @@
from astropy import log
from astropy.nddata import NDData
+from .wcs import Spectrum1DLookupWCS, Spectrum1DLinearWCS
+
import numpy as np
@@ -20,8 +22,8 @@ class Spectrum1D(NDData):
@classmethod
- def from_array(cls, disp, flux, uncertainty=None, mask=None, flags=None, meta=None,
- units=None, copy=True):
+ def from_array(cls, disp, flux, dispersion_unit=None, uncertainty=None, mask=None, flags=None, meta=None,
+ unit=None, copy=True):
"""Initialize `Spectrum1D`-object from two `numpy.ndarray` objects
Parameters:
@@ -72,14 +74,20 @@ def from_array(cls, disp, flux, uncertainty=None, mask=None, flags=None, meta=No
if disp.ndim != 1 or disp.shape != flux.shape:
raise ValueError("disp and flux need to be one-dimensional Numpy arrays with the same shape")
-
- return cls(data=flux, wcs=disp, *args, **kwargs)
+ spec_wcs = Spectrum1DLookupWCS(disp, unit=dispersion_unit)
+ return cls(data=flux, wcs=spec_wcs, unit=unit)
@classmethod
- def from_table(cls, table, uncertainty=None, mask=None, disp_col='disp', flux_col='flux'):
- flux = table[flux_col]
- disp = table[disp_col]
- return cls(data=flux, wcs=disp, uncertainty=uncertainty, mask=mask)
+ def from_table(cls, table, mask=None, dispersion_column='disp', flux_column='flux', uncertainty_column=None):
+ flux = table[flux_column]
+ disp = table[dispersion_column]
+
+ if uncertainty_column is not None:
+ uncertainty = table[uncertainty_column]
+ else:
+ uncertainty = None
+
+ return cls.from_array(flux=flux.data, disp=disp.data, uncertainty=uncertainty, dispersion_unit=disp.units, unit=flux.units)
@@ -94,7 +102,7 @@ def from_ascii(cls, filename, uncertainty=None, mask=None, dtype=np.float, comme
if raw_data.shape[1] != 2:
raise ValueError('data contained in filename must have exactly two columns')
- return cls(data=raw_data[:,1], wcs=raw_data[:,0], uncertainty=uncertainty, mask=mask)
+ return cls.from_array(disp=raw_data[:,0], flux=raw_data[:,1], uncertainty=uncertainty, mask=mask)
@classmethod
def from_fits(cls, filename, uncertainty=None):
@@ -115,7 +123,18 @@ def flux_setter(self, flux):
@property
def disp(self):
#returning the disp
- return self.wcs
+ if not hasattr(self.wcs, 'lookup_table'):
+ self.wcs.create_lookup_table(np.arange(len(self.flux)))
+
+ return self.wcs.lookup_table
+
+ @property
+ def dispersion_unit(self):
+ return self.wcs.unit
+
+ @property
+ def flux_unit(self):
+ return self.unit
def interpolate(self, new_disp, kind='linear', bounds_error=True, fill_value=np.nan):
@@ -245,7 +264,7 @@ def slice_dispersion(self, start=None, stop=None):
"""
# Transform the dispersion end points to index space
- start_index, stop_index = self.disp.searchsorted([start, stop])
+ start_index, stop_index = self.wcs.dispersion2pixel([start, stop])
return self.slice_index(start_index, stop_index)
View
1  specutils/wcs/__init__.py
@@ -0,0 +1 @@
+from .specwcs import Spectrum1DLinearWCS, Spectrum1DLookupWCS
View
119 specutils/wcs/specwcs.py
@@ -0,0 +1,119 @@
+# This module provides a basic and probably temporary WCS solution until astropy has a wcs built-in
+#This is all built upon @nden's work on the models class
+
@hamogu Collaborator
hamogu added a note

I appreciate the work you put in this. But do we really need a temporary solution before nede's work is merged into astropy? It's probably going to happen very soon.

@keflavich Owner

Will the framework change much with incorporation of nden's work? I think we still need much of this, but I really don't know much about the models framework yet.

@wkerzendorf Owner

@hamogu. Well it's not quite a temporary solution.I was involved in the API shaping of some of the wcs part of @nden 's model framework so I know that the things I have written here are not really duplicating the @nden's work. I think not much will change in this specwcs - but will be greatly enriched by actually inheriting from the astropy-model. One example is that we will get an easy way to fit these models - so we can ditch IRAF for the wavelength calibration part pretty much right away ;-)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+# Add when nden's models are merged
+# from astropy import models
+
+import numpy as np
+
+from astropy.utils import misc
+from astropy.io import fits
+
+class NDenModelsPlaceHolder(object):
+
+ pass
+class BaseSpectrum1DWCSError(Exception):
+ pass
+
+class BaseSpectrum1DWCS(NDenModelsPlaceHolder):
+ """
+ Base class for a Spectrum1D WCS
+ """
+
+ def pixel2dispersion(self, pixel_index):
+ """
+ This should be the forward transformation in normal WCS classes
+ """
+ return self(pixel_index)
@keflavich Owner

self[pixel_index] ? Will NDenModelsPlaceHolder have a __call__ method or a __getitem__ method? I think flushing out the inherited class would help a little.

@wkerzendorf Owner

It does - I can just copy it accross from her PR. I just didn't want to do that so we don't have that in the commits. To give you an overview: https://github.com/nden/astropy/blob/master/astropy/models/models.py. I think it's days away from being merged.

@keflavich Owner

Got it, that's good enough. Thanks.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+
+ def dispersion2pixel(self, dispersion_value):
+ """
+ This should be the inverse transformation in normal WCS classes
+ """
+ return self.invert(dispersion_value)
+
+ def create_lookup_table(self, pixel_indices):
+ self.lookup_table = self(pixel_indices)
@hamogu Collaborator
hamogu added a note

Should we default that to use all pixels?

def create_lookup_table(self, pixel_indices = None):
if pixel_indices is None:
    pixel_indixes = np.arange(self.dispersion0, self.dispersion_delta * len(self.flux))
@keflavich Owner

+1: I think that is the intent illustrated in the example at the top of this PR. Though I think it should be

def create_lookup_table(self, pixel_indices = None):
    if pixel_indices is None:
        pixel_indices = np.arange(self.npix) + 1
    self.lookup_table = self(pixel_indices)

?

@wkerzendorf Owner

So two comments to that. I thought about adding an npix or equivalent to the WCS class - but decided against it because it "felt" cleaner to keep NAXIS1 and WCS separate. Maybe an interesting compromise might be to tell wcs about its parent and thus it could access this information.

The second design decision was to give the first pixel the index 0. I know, that IRAF doesn't do it that way (that's why CRPIX1=1 for the first pixel).

@keflavich Owner

I agree with giving the first pixel index 0.... does that mean when you load in CRPIX from a FITS header, you will do:
self.crpix = header.get('CRPIX')-1 ?

@wkerzendorf Owner

@keflavich precisely

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+
+
+
+
+class Spectrum1DLookupWCS(BaseSpectrum1DWCS):
+ """
+ A simple lookup table wcs
+
+ Parameters
+ ----------
+
+ lookup_table : ~np.ndarray
+ lookup table for the array
+ """
+
+ def __init__(self, lookup_table, unit=None):
+ self.unit = unit
+ self.lookup_table = lookup_table
+
+ #check that array gives a bijective transformation (that forwards and backwards transformations are unique)
+
+ if len(self.lookup_table) != len(np.unique(self.lookup_table)):
+ raise BaseSpectrum1DWCSError('The Lookup Table does not describe a unique transformation')
@keflavich Owner

I like this, but I think it means that you can't use Spectrum1DLookupWCS with Spectrum1DEchelle as @hamogu has suggested, since it will require 2D-dispersion. That's probably OK, I'm just noting it.

@wkerzendorf Owner

So imho I think I want to focus on a Spectrum1D first. If you open up two dimension it opens up a whole can of worms (IFUs, Echelle, spectral-spatial spectra, ..... ).

@keflavich Owner

Yes, that's fine. Just something to keep in mind when we extend to >1D

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+
+
+ 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]
+
+ def invert(self, dispersion_values):
+
+ return np.searchsorted(self.lookup_table, dispersion_values)
+
+
+
+class Spectrum1DLinearWCS(BaseSpectrum1DWCS):
+ """
+ A simple linear wcs
+
+ """
+
+
+ @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):
+ self.unit = unit
+ 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)
+
+ 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
+
+
+#### EXAMPLE implementation for Chebyshev
+#class ChebyshevSpectrum1D(models.ChebyshevModel):
+class ChebyshevSpectrum1D(BaseSpectrum1DWCS):
+
+ @classmethod
+ def from_fits_header(cls, header):
+ pass
+ ### here be @hamogu's code ###
+ #degree, parameters = hamogu_read_fits(header)
+ #return cls(degree, **parameters)
+
+
+
+
+
+
+
Something went wrong with that request. Please try again.