Closed zjzhang42 closed 3 years ago
Specutils has a radial velocity attribute. You can set the RV to some value with units:
spectrum.radial_velocity = 18.0 * u.km/u.s
If that is the final operation you apply to a spectrum, you're done. However, if you apply additional operations specutils does not know whether you want the operation to be conducted in the rest frame or in the new frame, and for many operations it appears to just forget about the radial_velocity (which seems like a bug or at least unexpected behavior). I have found a workaround that I have to apply the following to force specutils to just adopt the RV-shifted wavelength:
spectrum.radial_velocity = 123.4 * u.km/u.s # Your RV here, requires units!
spectrum._copy(
spectral_axis= spectrum.wavelength.value * spectrum.wavelength.unit
)
which makes it simply put the RV into wavelength units and forget about keeping track of the reference frame.
Here is the line of code showing that approach in .barycentric_correct()
:
https://github.com/OttoStruve/muler/blob/4bc7cab4f695feb0f441348150251f593a436cdd/src/muler/echelle.py#L192
By the way, you could imagine adding a method called .rv_shift(val)
that shifts by some value val
, and implements this copying on the backend so the user never has to think about this workaround. Would you be open to implementing that method or adding it as a feature request?
I actually did exactly this a few weeks ago with a function called velocity_shift(val)
where val
is the velocity shift in km/s but haven't pushed the changes until now. I just pushed it to the add_resample_list branch and you can see under the add_resample_list pull request.
Okay I've modified my old velocity_shift
definition into rv_shift
. I split off most the code from barycentric_correct
which does the shift in the correct way, and moved it into rv_shift
which barycentric_correct
now calls. I've committed these changes to the add_resample_list branch.
For reference here is the function for me right now:
def rv_shift(self, velocity):
"""
Shift velocity of spectrum in astropy units (or km/s if input velocity is just a float)
"""
if type(velocity) == float: #If supplied velocity is not using astropy units, default to km/s
velocity = velocity * (u.km/u.s)
try:
self.radial_velocity = velocity
return self._copy(
spectral_axis=self.wavelength.value * self.wavelength.unit
)
except:
log.error(
"rv shift requires specutils version >= 1.2, you have: {}".format(
specutils.__version__
)
)
raise
I do find that it does update the spectral_axis
where I can see the wavelengths shift, but then when I go try to do band math, specutils seems to forget the two spectra I am dividing by each other are now on different wavelength grids. For example if I take a spectrum, shift a copy of it in velocity with the function above, and then divide the shifted by unshifted spectrum, the result is always 1. This sounds a lot like what you said Gully about specutils ignoring theradial_velocity
, but now it seems to be similarly ignoring the wavelengths. My impression was that specutils should be doing some interpolating of the flux when the wavelengths are slightly different between spectra. Am I missing something here?
Awesome work, thanks for this contribution @kfkaplan!
The behavior you are describing is summarized at the bottom of our error propagation tutorial
That tutorial also shows how to overcome this limitation. I would like to make a combine function that automatically deals with these alignment issues. @joelburke02 was working on this, and I have some new ideas for how to go about it with scipy's binned_statistic methods.
woohoo, with Kyle's new PR #45 merged, this Issue is now closed!
Does
muler
has a built-in function to RV-shift a given spectrum back to its rest frame based on a user-defined RV value? If yes, how do we use such a function? Thanks.