Spec WCS #22

Merged
merged 11 commits into from Apr 7, 2013

Conversation

Projects
None yet
4 participants
Member

wkerzendorf commented Apr 3, 2013

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 Apr 3, 2013

specutils/wcs/specwcs.py
+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

keflavich Apr 3, 2013

Contributor

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

wkerzendorf Apr 3, 2013

Member

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

keflavich Apr 3, 2013

Contributor

Got it, that's good enough. Thanks.

@keflavich keflavich commented on the diff Apr 3, 2013

specutils/wcs/specwcs.py
+
+ 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

keflavich Apr 3, 2013

Contributor

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

wkerzendorf Apr 3, 2013

Member

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

keflavich Apr 3, 2013

Contributor

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

Contributor

keflavich commented Apr 3, 2013

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

Member

hamogu commented Apr 3, 2013

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 Apr 3, 2013

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

hamogu Apr 3, 2013

Member

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

keflavich Apr 3, 2013

Contributor

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

wkerzendorf Apr 3, 2013

Member

@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 ;-)

Contributor

keflavich commented Apr 3, 2013

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 Apr 3, 2013

specutils/wcs/specwcs.py
+ """
+
+ 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

hamogu Apr 3, 2013

Member

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

keflavich Apr 3, 2013

Contributor

+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

wkerzendorf Apr 3, 2013

Member

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

keflavich Apr 3, 2013

Contributor

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 ?

Member

wkerzendorf commented Apr 3, 2013

(@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.

Contributor

keflavich commented Apr 3, 2013

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?

Member

wkerzendorf commented Apr 3, 2013

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

Contributor

keflavich commented Apr 4, 2013

@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!

Member

hamogu commented Apr 4, 2013

@wkerzendorf that convinces me!
Cool + 1!

Member

wkerzendorf commented Apr 5, 2013

@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?

@astrofrog astrofrog and 1 other commented on an outdated diff Apr 5, 2013

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

astrofrog Apr 5, 2013

Owner

This should be a relative import

@wkerzendorf

wkerzendorf Apr 5, 2013

Member

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

@astrofrog astrofrog commented on an outdated diff Apr 5, 2013

specutils/wcs/__init__.py
@@ -0,0 +1 @@
+from specwcs import Spectrum1DLinearWCS, Spectrum1DLookupWCS
@astrofrog

astrofrog Apr 5, 2013

Owner

This should be a relative import

Member

wkerzendorf commented Apr 5, 2013

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

Contributor

keflavich commented Apr 5, 2013

I say go ahead, lets get some working examples and tests with readers, then
we can see whether any major changes are needed.

@wkerzendorf wkerzendorf added a commit that referenced this pull request Apr 7, 2013

@wkerzendorf wkerzendorf Merge pull request #22 from wkerzendorf/specwcs
Spec WCS - first implementation
362ef5f

@wkerzendorf wkerzendorf merged commit 362ef5f into astropy:master Apr 7, 2013

@astrofrog astrofrog added a commit to astrofrog/specutils that referenced this pull request Oct 12, 2013

@astrofrog astrofrog Merge pull request #22 from astrofrog/improve-travis-config
Test with different Numpy versions, only test Sphinx with Python 2.7.
9376c25
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment