New method brainstorming: Ideas for more actions you can do to a spectrum, ranked

[gully]

Here is a brainstorm of things we can do to a spectrum. For each item we will put the estimated Implementation Difficulty (:keycap_ten: being the harder) and estimated Benefit to the community (:keycap_ten: being very valuable to the community).

1) .mask_telluric_regions( ) Telluric regions are sometimes either uncorrected, or corrected so poorly that their residuals would ruin a science application. Simply masking these regions can be an adequate solution for many applications. This method would require having a list of spectral regions affected by telluric lines, presumably a line list or experimentally determined mask. Which approach we choose will have differences in complexity to implement, with line lists being possibly outside the scope of muler. Implementation difficulty: :seven: Benefit to community: :eight:

2) .refine_wavelength_calibration(method='telluric') Sometimes the wavelength calibration is poor and can be improved/standardized based on either information in the stellar spectrum itself or on the conspicuous telluric absorption lines. In practice this method is delicate and may be instrument-specific, and really should be the domain of the instrument pipelines that have access to more diagnostic information in the first place. Implementation difficulty: :nine: Benefit to community: :six:

3) .to_pandas() and/or .to_excel(). Often you want to post-process a spectrum and then simply save it as a csv file to send to a collaborator or post on a website. Or you may have some existing analysis code that assume dataframes. So being able to convert to pandas would help. We have this method in astropy Tables, so it should be straighforward to implement. Implementation difficulty: :two: Benefit to community: :five:

4) .combine_spectra() Sometimes you have two or more spectra of the same object, wavelength range, and instrumental setup, and you want to combine them into a single high signal-to-noise ratio composite spectral template. In principle we can currently just add spectra together, one after another net_spec = spec1.add(spec2) or net_spec = spec1 + spec2 + spec3. The problem with this naive approach is what-to-do-with-the-metadata: should we be a bit smart about how to combine the dates? Should we warn the user to conduct RV registration first? The implementation difficult depends on how "smart" we want to be. To some extent the implementation is made easier by specutils' built-in methods. Should this be a method on a SpectrumList instead? Implementation difficulty: :six: Benefit to community: :five:

5) .cross_correlate(template) It's common to want to cross correlate an observed, noisy spectrum, with a high SNR or noise-free template. RV applications and other detection-based techniques apply cross-correlation. specutils already has this method as a built-in, so we would only have to provide a lightweight wrapper. The main trick is making sure that the pre-processing and flattening steps have occurred so the cross-correlation is meaningful and effective. So we would want a high quality tutorial to go along with the implementation.
Implementation difficulty: :four: Benefit to community: :six:

[gully]

Open question to current or prospective users: which of these methods appear most useful to your science? Which other methods would you like to see? I suspect a key method folks would want to see is "find best fit template" or "estimate teff" or something along those lines. I think those methods are better suited to https://gollum-astro.readthedocs.io/en/latest/ which has all the knowledge of the models. Let me know your thoughts--- which methods belong in muler and which in gollum.

[kfkaplan]

Going off functions I have in my old IGRINS plotspec code, here are some suggestions if they don't already exist:

• Continuum fitting (e.g. a running median) for either subtracting the continuum or normalizing by it. This is especially useful for dealing with things like planetary nebula where you are interested in the emission lines and want to subtract out any continuum.
  • Implementation difficulty: 8
  • Benefit to community: 7
• Bad pixel (nan) filling with the median of nearby pixels. This is limited only to small gaps, so large gaps (e.g. masked telluric) don't get filled.
  • Implementation difficulty: 4
  • Benefit to community: 5
• Output to a .fits file
  • Implementation difficulty: 2
  • Benefit to community: 9
• Output to a .txt file
  • Implementation difficulty: 2
  • Benefit to community: 9
• Built in plotting function(s)
  • Implementation difficulty: 8
  • Benefit to community: 8
• Dereddening to correct for Galactic extinction
  • Implementation difficulty: 6
  • Benefit to community: 5

[kfkaplan]

I think we should add "band math" to the EchelleSpectrum and EchelleSpectrumListclasses including addition, subtraction, division, and multiplication. This will help simplify scripting things later.

• Implementation difficulty: 3
• Benefit to community: 7

I would be happy to implement this.

[gully]

Regarding "band math"--- EchelleSpectrum should already have support for adding, subtracting multiplying, and dividing for two vectors, since it inherits from specutils, which already supports two argument band math. By band math you are implying support for more than a mere two spectra, with 3 or more groups of spectra combined instantaneously. EchelleSpectrumList does not currently support band math for any number of spectra, but could and probably should, so that would be a welcomed addition! Thank you! 🙏

There's an Open Pull Request attempting to add a "combine" function: #36

One hiccup arises with spectra possessing different (enough) wavelength solutions. I documented this RV-based "gotcha" in this tutorial where I come up with this observation:

Whoa! Specutils pretends like the signals are aligned to Spectrum 1. That is probably not the desired behavior for such an extreme offset as this one, but may be “good enough” for spectra that are either exactly aligned or within a pixel. It depends on your science application. PRV applications should not just round-to-the-nearest pixel, since they are trying to infer changes at sub-pixel levels.

This "gotcha" bit me yesterday. I thought one way around it is to extend specutils' band math with an optional keyword such as RV_policy='resample' or RV_policy='error' or RV_policy='warn', which would institute different policies depending on the keyword. resample would simply resample the second spectrum to the first. error would refuse to add the two-or-more spectra unless the spectral axes were exactly the same, etc.

Based on these considerations I would say the implementation difficulty is a bit higher in order to consider the effect of mis-aligned axes and how to propagate the error under resampling.

[kfkaplan]

Oh I didn't realize the EchelleSpectrum class was inherited from Specutils and already had bandmath built in! I basically meant just arithmetic between two spectra with the same wavelength solution. I'll look into seeing if EchelleSpectrumList can have the same sort of bandmath added to it.

As for dealing with different wavelength solutions, it's not something I need at the moment but probably will be useful to people, especially if working on RVs and precise wavelength solutions are required.

[kfkaplan]

Turns out it was pretty simple to propagate the band math from EchelleSpectrum to EchelleSpectrumList. I will include it in the pull request I will send later today. Here is what the definitions look like for EchelleSpectrumList:

    def __add__(self, other):
        """Bandmath addition
        """
        spec_out = copy.deepcopy(self)
        for i in range(len(self)):
            spec_out[i] = self[i] + other[i]
        return spec_out

    def __sub__(self, other):
        """Bandmath subtraction
        """
        spec_out = copy.deepcopy(self)
        for i in range(len(self)):
            spec_out[i] = self[i] - other[i]
        return spec_out

    def __mul__(self, other):
        """Bandmath multpilication
        """
        spec_out = copy.deepcopy(self)
        for i in range(len(self)):
            spec_out[i] = self[i] * other[i]
        return spec_out

    def __truediv__(self, other):
        """Bandmath division
        """
        spec_out = copy.deepcopy(self)
        for i in range(len(self)):
            spec_out[i] = self[i] / other[i]
        return spec_out