How to handle echelle spectra

[hamogu]

This issue is here to raise the issue of echelle spectra. They are also 1D spectra in the sense that they only have one axis (wavelength), however, they have multiple orders. They should be kept in one object, as all the echelle order typically have the same metadata (GAIN, OBS-TIME etc.).

We could subclass Spectrum1D for that purpose or just allow dispersion and flux to be 2 D arrays and add a method like flatten(some options here) to merge all echelle orders.

Clearly flatten is a major undertaking that requires many different options and methods (cross-correltate features to improve wavelength solution? How to set the merge points? Rebin to common grid?), but the first step is to decide if echelle spectra need a separate object.

[andycasey]

Could we just have the readers be written such that an echelle spectrum returns a list of Spectrum1D objects?

[andycasey]

When considering the metadata, I wonder what would be the best way to keep this in a single NDData object.. thoughts?

[keflavich]

As a first step, I see no harm in sticking to a "list of Spectrum1D objects", but I think it would be worthwhile to make a "list of spectra" wrapper, or a Spectra1D class, that deals with broadcasting some operations. I've implemented that in pyspeckit: http://pyspeckit.bitbucket.org/html/sphinx/spectrum.html?highlight=spectra#pyspeckit.spectrum.classes.Spectra and frequently make use of it for quick-look plots of a whole echelle spectrum. Then again, there's no reason this can't be stored as an NDData object as long as the x-axis class can be similarly n-dimensional. Can NDData have N-d metadata?

[hamogu]

I see three ways to handle echelle spectra. Let me make a list here, so that I can add a few pros and cons to each of them:

1. Subclass Spectrum1D to make a Spetrum1DEchelle
2. use a list of Spectrum1D objects
3. Use a Spectrum1D object with a 2D flux and wavelength array

(2) duplicates the meta information (the meta field, which in most cases will hold information ffrom the fits header and the units of the dispersion and the flux). I often want to do things that affect all spectra in the list, e.g. change the units of the flux (if they were not set properly by the reader), change the name of the object (e.g. from HD 1234 to AA Lup), recalculate the airmass. Thus, it seems natural for me to hold that info in the same object.

(1) The name of the object already says what's in there. Could add methods for merging echelle orders.

(3) Is almost like one (you get individual orders as spec.array[3,:]) but requires no new programming at this point.

[keflavich]

I think (1) is the right approach, but perhaps implemented internally as an NDArray.

[andycasey]

I''m in favour of (3) actually. (2) has the disadvantage that all of the items in the list need to have the same metadata, and if there are any changes needing to be made to that metadata then they likely need to be made for every item in the list.

If we did (1), the question arises on how you would access the individual orders. Would it be my_echelle = Spectrum1DEchelle(...) and then my_echelle[0].disp to access the dispersion array for the first order, or would it be my_echelle.disp[0, :]? If it's the former then we're just extending list, adding a method or two and calling it Spectrum1DEchelle.

In terms of (1), what about the case of multi-beam spectrographs with 2 or 4+ completely separate orders? It's not so much that they are necessary echelle orders, but simply multiple segments of spectra (overlapping or separate) from the same exposure. You could have Spectrum1DEchelle and Spectrum1DMultiBeam but I think that's adding an unnecessary inception level. We don't have to go deeper.

Since Spectrum1D inherits from NDData -- and I assume that NDData would be able to have a 2D flux and wavelength array -- shouldn't we try to leverage from that and start from (3)? We can always have a merge util somewhere where multiple overlapping orders in this format could be fed into it.

[hamogu]

Most spectra we will use in astropy will come as fits files for the foreseeable future. Thinking of it in term of fits files, I suggest that everything that can be in a single fits extension can be usually be in a single object (at least initially - echelle orders may be merged by the user or multi-object spectra may be split in Spectrum1D objects to get one per object).

On the other hand, if a spectrograph like UVES has different arms each with different headers (GAIN, read-out, obstime) etc. I would return separate Spectrum1D objects initially and leave it to the user to merge them into a single spectrum later.

[hamogu]

After thinking about this some more, I now tend towards a Spectra1D class (could also be called Spectrum1DEchelle as I suggested above), that is almost identical to Spectrum1D except that it has flux and disp as 2D arrays. The main reason for this is that then Spectrum1D can test in the initialization that input arrays are 1D - I think that is important to catch many common types of input error.

@andycasey: A single order would then be my_echelle.disp[0, :]. That's OK with me.

It seems that we have exchanged most arguments on how to handle this. @wkerzendorf: I am happy to help you on an API document - do you plan to write one? I know that came up in the discussion of another issue as well...

I am sure further discussion would happen there.

[keflavich]

Would this same approach be used for multiple spectra with the same disp axis, then? i.e., how would you handle a list of spectra all taken with the same spectrograph & configuration but of different objects?

I'm concerned there may be ambiguity between Echelle spectra (2d data, 2d dispersion) and "blocks" of 1D spectra with the same dispersion axis (2d data, 1d dispersion). I think it would be better to have a Spectrum1DEchelle and Spectra1D to distinguish the structures explicitly.

[wkerzendorf]

@hamogu The last thing that I think is still on the table is the discussion about the specwcs - I think if that is solved (in one way or another). We should start putting the API together. Thanks for offering to help!

[iancze]

Hi everyone! I am new to the specutils project, but have happily used it to load some echelle spectra (specifically from TRES, which coincidentally is also in the example spectra). I was wondering if there has been any more progress on the API document mentioned in this thread?

Related: I would like to help develop the capabilities of this package (and specifically echelle functionality), but I'm a little unsure where to start. Is there a -dev mailing list for specutils to ask general questions? Thank you!

[keflavich]

@iancze - right now, there's no mailing list as there are only a handful of us. I've been sitting in #specutils on freenode.

I think the first step in dealing with any data format is to build a reader. If you can already read in Echelle spectra, perhaps a writer is a good next step.

[hamogu]

@iancze As you might have seen when you read the TRES file, this discussion stalled and the current implementation just returns a list of Spectrum1D objects. That's not necessary the best solution, but it was easy to implement and at this point it's probably better to have something easy than not have something.

[hamogu]

@keflavich : Personally, I don't think writers are the highest priority at this point. We can read spectra, but we have little functionality to do things with them (coadd, shift, cross-match, fit) and nothing to extract them from raw data. Personally, I suggest adding more operations on spectra, for Echelle spectra e.g. combining orders, fitting and removing a blaze function etc.

But clearly my view is biased by what I would like to use ;-)

[keflavich]

@hamogu - I agree, writers are a low-ish priority, but I've had need for them recently. For the use cases you've mentioned, I use pyspeckit, so I'm less worried about implementing them, but I do plan to get a student working on incorporating pyspeckit utilities into specutils (or vice-versa).

You're right that we really need some astropy tools to extract spectra from raw data. I think we're waiting on a handful of developments from, e.g., ccdproc to help enable that.

Perhaps we should add a page on the wiki describing the tools needed (combining orders, fitting & removing blaze function, etc.) and perhaps a hint of how to implement them, e.g. links to the appropriate IRAF packages. The more we can make individual, digestible projects clear, the more likely we'll all work on implementing them. I'll get this started...

[wkerzendorf]

• Regarding the specutils-dev mailing list - I'll create one in the next few days and add you (will also announce it to astropy-dev). I'll invite the people here on this list.
• It is very trivial to implement reading of Echelle data in a two dimensional mode. if you glance over at specwcs.py you can see what we have done for 1D-spectra for the wavelength solution. For a 2D-echlle spectrum one just needs to write a WCS that uses the row to know which 1D wcs solution to pick and then evaluate the right column pixel.
• extracting spectra from raw data does not really wait on any other development. The right way to go forward is implement the Horne 86 algorithm (http://adsabs.harvard.edu/abs/1986PASP...98..609H) - which is one of the most widely used ones. The only data input that would be required is an ndarray (which can be ccdproc'ed or not). The fitting should be implemented using the Models framework. This project is part of Google summer of Code, but might not necessarily be picked up as it is relatively advanced.
• coadding, shifting, crossmatching of spectra is at the moment limited by our interpolation capacity. Interpolation is intricately linked with WCS and the main part missing to do a good general interpolation is a description of the pixel reference frame.

[keflavich]

@wkerzendorf - be cautious about trivializing echelle reading. Overlapping orders means that the data may be in an effectively useless format - a list of 1d spectra, say? I don't think writing a WCS is ever trivial.

Could you write a breakdown of the individual tasks in the Horne 86 algorithm on a wiki page? If we see the individual tasks laid out, perhaps we can start taking them on as bite-sized code snippets.

[wkerzendorf]

@keflavich reading the echelle data is trivial. Combining orders is a non-trivial task but is not part of reading. I have extensively worked with Echelle spectra and in most cases had them as a list of 1D spectra (i.e. measuring line strengths ...) . Things like removing the blaze function can be done much better when they are in a list of 1D spectra - so I would not call that a useless format.

I believe that we can not all work on an optimal extraction routine together - but someone needs to read it and program it. @mhvk has written this algorithm in Fortran - maybe he can comment on the way to move forward on that.

[keflavich]

A wrapper around the fortran code would be very useful.

If the echelle data is stored in a well-understood format and has already been extracted, then I guess it's trivial, yes.

[ladsantos]

Hey folks, shall I unearth this 3-year old issue? I will start working on an echelle spectra extraction code for a new instrument, and I thought this might be a good opportunity to contribute to specutils. Any suggestions on where to start? What about this Fortran code?

[wkerzendorf]

@RogueAstro what specifically do you need for your echelle spectra.

[ladsantos]

@wkerzendorf Sorry for not being specific, I later realized that my project probably does not fit well to the objectives of specutils, I think. I will mostly deal with the initial data reduction (calibration and so on), and then later move on to writing a echelle spectral extraction routine.

[crawfordsm]

@RogueAstro You may want to check out ccdproc and pyhrs for deailing with reduction of echelle spectra. I'd be happy to discuss your needs with your data reductions needs with you.

https://github.com/astropy/ccdproc https://github.com/saltastro/pyhrs

[jgussman]

Heyyo! Is there a way to read echelle spectra now that are in fits format now? Or is still the "best" way making a 1Dspectra list? I don't need to reduce the data. Thanks

[matiruizdiaz]

Hello! I didn't understood at all how to read a echelle spectra. Apparently it is trivial according the comments but there is no way to get a image of my spectra. There is some documentation or anything? Thank you

[hamogu]

@matiruizdiaz Thanks for asking. As you may see from the dates and the text above, this is an issue where we discussed how to handle echelle spectra in the code about a decade ago. That's been implemented, and thus the discussion is closed.

General documentation on specutils is here: https://specutils.readthedocs.io/en/stable/ Unfortunately, I can't help you in more detail, because your post does not provide any details of what you want.

Specutils is for reading in reduced data. That is, you've already run an instrument pipeline and you have a set of fits files with e.g. one order per file, or, one file that holds several orders. However, there is no "image of my spectra". Since we are dealing with reduced spectra, each spectrum is (essentially, there a slightly different formats), a table of wavelength vs. flux. That's a 1D format, not an image. That's what's referred to as "trivial" above, as in "it is easy to use the existing code to read a 1D table of wavelength vs. flux from a fits file". Of course, dealing with echelle spectra in general is not trivial at all!

Maybe you are looking for a way to reduce raw data, that's still in the format of the images captured by the CCD? You might want to a hve a look at specreduce or other similar data reduction pipelines, such as [PipeIt](https://specreduce.readthedocs.io/en/latest/_, DRAGONS, pyreduce, or any pipeline software that's commonly used for your instrument.

Either way, I'm not sure this issue is the right place to discuss echelle data reduction. If you have a specific suggestion what this code should do, or what should be documented where, please open a new issue. For a more general discussion, I recommend the community discussion forum: https://community.openastronomy.org/c/astropy/8