<p>I am trying to learn how to derive the line of sight velocity distributions from galaxy spectrum to extract kinematical parameters like rotational velocity and velocity dispersion. Therefore I need to do some prepare my spectra. For this example I took a stellar spectra from the Miles-catalogue <a href="http://miles.iac.es/pages/stellar-libraries/the-catalogue.php" rel="nofollow noreferrer">http://miles.iac.es/pages/stellar-libraries/the-catalogue.php</a>. I tried to rebin the spectra onto a logarithmic wavelength scale and then I tried to remove the continuum as described by Bender et. al 1994, <a href="https://ui.adsabs.harvard.edu/abs/1994MNRAS.269..785B/abstract" rel="nofollow noreferrer">Line-of-sight velocity distributions of elliptical galaxies</a> (also available <a href="https://academic.oup.com/mnras/article/269/3/785/1002169" rel="nofollow noreferrer">here</a> and in <a href="https://www.researchgate.net/publication/228075911_Line_of_sight_velocity_distributions_of_elliptical_galaxies" rel="nofollow noreferrer">researchgate</a>).</p>
<p>Here is the relevant part of the paper:</p>
<blockquote>
<p>3.1 Reduction and derivation of line-of-sight velocity
distributions</p>
<p>The CCD spectra were de-biased, dark-subtracted and flatfielded in the usual way. Hot pixels and cosmic ray events
were removed with a <span class="math-container">$\kappa$</span>—<span class="math-container">$\sigma$</span> clipping procedure. Then, the spectra were rebinned to a logarithmic wavelength scale and the
mean sky spectrum during each exposure was derived by averaging several lines from the edges of the CCD spectra. After
subtraction of the sky spectra, spectra of the same galaxy
taken at identical slit positions were centred and added. Then,
the spectra were rebinned along the slit in order to guarantee
a signal-to-noise ratio that allowed the derivation of kinematic
parameters. As the final step of the pre-processing, the galaxy
continuum was removed in two steps. First, a fourth- to sixthorder polynomial was fitted to each line-by-line spectrum.
Then, the rms variation <span class="math-container">$\sigma_s$</span> of the spectrum around the fitted
polynomial was calculated and the fit was repeated including
only those pixels of which the value fell in a range between 0<span class="math-container">$\sigma_s$</span>
and 2<span class="math-container">$\sigma_s$</span> around the first fit. Finally, the spectrum was divided
by this second fit polynomial. This procedure avoids the fact
that the stronger absorption lines disturb the fit and leads to
a satisfactory approximation of the galaxy ‘continuum’. This
method of continuum removal was found to be superior to the
usually applied procedures which include filtering in Fourier
space. The reason is that filtering in Fourier space may change
the profiles of the absorption lines in such a way that the
derived <span class="math-container">$H_4$</span> values (see below) may be systematically offset
with respect to the true values. The appropriate order of the...</p>
</blockquote>
<p>(original <a href="https://i.sstatic.net/LaSae.jpg" rel="nofollow noreferrer">screenshot</a>)</p>
<p>I am not certain if I performed this correctly. Here is my code:</p>
<pre><code># read in 1D spectral data, which is conventionally stored in FITS format
# read in stellar template spectra from miles library http://miles.iac.es/
# by Vazdekis (2010, MNRAS, 404, 1639)
miles_dir = r&quot;C:\Users\reich\OneDrive\Dokumente\Uni\Bachelorarbeit\Python\ppxf\MILES_library_v9.1_FITS&quot;
file = miles_dir + '/s0427.fits'    # s0427 K2III-star
# header data unit (HDUs) list are the highest level component of the FITS file structure
# hdu.info() summarizes content of opened fits file
hdu = fits.open(file)
# HDU objects consist of a header and a data unit/attributes
# to read the entire header one can use print(repr(hdr))
hdr = hdu[0].header
data = hdu[0].data
# extract numpy arrays for flux and wavelength from data
flux_raw = data[0]
w = WCS(hdr, naxis=1, relax=False, fix=False)
lam = w.wcs_pix2world(np.arange(len(flux_raw)), 0)[0]
# logarithmic wavelength scale to obtain constant velocity scale
# for N pixels in wavelength-array, need to rebin spectra
# to atleast N pixels in log-wavelength
log_lam_raw = np.log(lam)
# logarithmic rebinning
log_lam = np.linspace(min(log_lam_raw), max(log_lam_raw), num=len(lam))

# remove continuum of spectra
def fit_polynom(log_lam, a0, a1, a2, a3, a4, a5):
    &quot;&quot;&quot;fit func = 2 nd order polynom&quot;&quot;&quot;
    return a0+a1*log_lam+a2*(log_lam)**2+a3*(log_lam)**3+a4*(log_lam)**4+a5*(log_lam)**5
# intitialize fit paramters
init_fit_param = [0, 0, 0, 0 ,0, 0]
# fit data using scipy.curve_fit()
fit_params, fit_covariances = curve_fit(fit_polynom, log_lam, flux_raw, p0=init_fit_param)
# calculate root mean square deviation
rms_array = np.sqrt((flux_raw - fit_polynom(log_lam, *fit_params))**2)
# only use flux values for 2nd fit which have a rms value &lt; 2
fit_params2, fit_covariances2 = curve_fit(fit_polynom, log_lam[np.argwhere(rms_array&lt;=2)][:,0], flux_raw[np.argwhere(rms_array &lt;= 2)][:,0], p0=init_fit_param)
# remove continuum by dividing flux by fit polynomial
flux = (flux_raw / (fit_polynom(log_lam, *fit_params2))) -1
</code></pre>
<p>I chose a fipolynom of order 5 this can obviously be varied, but shouldnt be the main concern. This is a plot of the flux over wavelength and flux over logarithmic rebinned wavelength with the fitpolynom.<a href="https://i.sstatic.net/2bKY5.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/2bKY5.jpg" alt="enter image description here" /></a></p>
<p>And here is a plot with the &quot;removed&quot; continuum, which does not seem quite right.
<a href="https://i.sstatic.net/DUEEL.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/DUEEL.jpg" alt="enter image description here" /></a></p>
<p>So can someone please tell me if a) I did the logarithmic rebinning correct and b) if my approach to the continuum removal makes sense or should be varied. I want to further work with the continuum removed spectra and extract line of sight velocity distributions (for galaxies ), but need to be sure that I understand the preprep correctly. Thank you very much</p>

CC BY-SA 4.0
2021-05-13T12:28:07.217Z
Why is you final spectrum $\sim 0$? The paper you refers to divides by the continuum fit, which should leave the spectrum close to 1.  Obviously if you divide (or subtract)  a small number by a (uncertain) small number at the blue wavelength edge then you could get unpredictable behaviour. The x-axis values don't make any sense. Logarithmic rebinning doesn't change the x-axis values/range it just adaptively changes the size of the bins so that each one is an equal step in $\log \lambda$.
ProfRob
CC BY-SA 4.0
2021-05-14T03:44:19.510Z
I found a copy of the paper with selectable text and replaced your screenshot. Some folks use screen readers that convert text to audio and *images of text* don't process usually, so that is generally discouraged in Stack Exchange.
uhoh
CC BY-SA 4.0
2021-05-14T07:42:27.553Z
I was biased by a plot in the paper I referred to because there the flux was arround 0, that is why I subtracted 1 from the whole flux array for this plot which I just fixed. So should I only work further without the blue wave length end in this case to avoid this behaviour ? For the logarithmic rebinning I think I am lacking some understanding of the phrases used, and I am sorry if this should be obvious. What exactly is meant by bin, and size of the bin in this context ? I got an array for the flux and one for the wavelength, I thought I had to reconstruct the wavelength array ?
trynerror
CC BY-SA 4.0
2021-05-14T07:56:08.897Z
My guess it that a bin in this context means a particular range of wavelengths?
trynerror
CC BY-SA 4.0
2021-05-14T09:18:40.520Z
Your input data (the MILES template spectrum) is binned by linear wavelength: each element of the wavelength vector specifies a "bin"; for example, if lam[i] = 4800 (Angstroms), then it defines a bin containing photons with wavelengths ranging from 4800 - deltaLam/2 to 4800 + deltaLam/2, where deltaLam is the wavelength spacing (the size of the bins).
Peter Erwin
CC BY-SA 4.0
2021-05-14T09:22:27.310Z
"Logarithmically rebinning" means converting your spectrum so the wavelength vector defines bins with a constant size in log(wavelength) rather than a constant size in linear wavelength, which is what you're starting with. Your log_lam is what you want, but you also need to actually rebin the *flux* (e.g., flux_raw) so that it agrees with the rebinned log wavelength bins.
Peter Erwin
CC BY-SA 4.0
2021-05-14T09:30:39.207Z
You might look at the "log_rebin" function in Michele Cappellari's pPXF package (in the ppxf_util.py module).
Peter Erwin
CC BY-SA 4.0
2021-05-14T10:35:12.337Z
Thank you very much for this helpful explanation, I have already started to work a bit with the ppxf package to test my work. Is there a command in python where I can see the exact definition of the log_rebin() function in the ppxf_util module ?
trynerror
CC BY-SA 4.0
2021-05-14T22:00:41.703Z
@trynerror It's defined in the ppxf_util.py module, which is pure Python; you can open the Python file in a text editor. (You may have to figure out where your system installs external Python packages in order to locate the ppxf package.)
Peter Erwin
CC BY-SA 4.0
2023-02-19T07:52:17.873Z
Related https://astronomy.stackexchange.com/questions/14518/why-do-linear-velocity-redshifts-correspond-to-linear-pixel-shifts-when-the-spec
ProfRob
<observational-astronomy><galaxy><data-analysis><spectra><python>
Logarithmic rebinning and continuum removal of stellar spectra
trynerror
# Logarithmic rebinning and continuum removal of stellar spectra

(Asked by: username_0 on 2021-05-13)

I am trying to learn how to derive the line of sight velocity distributions from galaxy spectrum to extract kinematical parameters like rotational velocity and velocity dispersion. Therefore I need to do some prepare my spectra. For this example I took a stellar spectra from the Miles-catalogue http://miles.iac.es/pages/stellar-libraries/the-catalogue.php. I tried to rebin the spectra onto a logarithmic wavelength scale and then I tried to remove the continuum as described by Bender et. al 1994, Line-of-sight velocity distributions of elliptical galaxies (also available here and in researchgate).
Here is the relevant part of the paper:

3.1 Reduction and derivation of line-of-sight velocity
distributions
The CCD spectra were de-biased, dark-subtracted and flatfielded in the usual way. Hot pixels and cosmic ray events
were removed with a $\kappa$—$\sigma$ clipping procedure. Then, the spectra were rebinned to a logarithmic wavelength scale and the
mean sky spectrum during each exposure was derived by averaging several lines from the edges of the CCD spectra. After
subtraction of the sky spectra, spectra of the same galaxy
taken at identical slit positions were centred and added. Then,
the spectra were rebinned along the slit in order to guarantee
a signal-to-noise ratio that allowed the derivation of kinematic
parameters. As the final step of the pre-processing, the galaxy
continuum was removed in two steps. First, a fourth- to sixthorder polynomial was fitted to each line-by-line spectrum.
Then, the rms variation $\sigma_s$ of the spectrum around the fitted
polynomial was calculated and the fit was repeated including
only those pixels of which the value fell in a range between 0$\sigma_s$
and 2$\sigma_s$ around the first fit. Finally, the spectrum was divided
by this second fit polynomial. This procedure avoids the fact
that the stronger absorption lines disturb the fit and leads to
a satisfactory approximation of the galaxy ‘continuum’. This
method of continuum removal was found to be superior to the
usually applied procedures which include filtering in Fourier
space. The reason is that filtering in Fourier space may change
the profiles of the absorption lines in such a way that the
derived $H_4$ values (see below) may be systematically offset
with respect to the true values. The appropriate order of the...

(original screenshot)
I am not certain if I performed this correctly. Here is my code:
# read in 1D spectral data, which is conventionally stored in FITS format
# read in stellar template spectra from miles library http://miles.iac.es/
# by Vazdekis (2010, MNRAS, 404, 1639)
miles_dir = r"C:\Users\reich\OneDrive\Dokumente\Uni\Bachelorarbeit\Python\ppxf\MILES_library_v9.1_FITS"
file = miles_dir + '/s0427.fits'    # s0427 K2III-star
# header data unit (HDUs) list are the highest level component of the FITS file structure
# hdu.info() summarizes content of opened fits file
hdu = fits.open(file)
# HDU objects consist of a header and a data unit/attributes
# to read the entire header one can use print(repr(hdr))
hdr = hdu[0].header
data = hdu[0].data
# extract numpy arrays for flux and wavelength from data
flux_raw = data[0]
w = WCS(hdr, naxis=1, relax=False, fix=False)
lam = w.wcs_pix2world(np.arange(len(flux_raw)), 0)[0]
# logarithmic wavelength scale to obtain constant velocity scale
# for N pixels in wavelength-array, need to rebin spectra
# to atleast N pixels in log-wavelength
log_lam_raw = np.log(lam)
# logarithmic rebinning
log_lam = np.linspace(min(log_lam_raw), max(log_lam_raw), num=len(lam))

# remove continuum of spectra
def fit_polynom(log_lam, a0, a1, a2, a3, a4, a5):
    """fit func = 2 nd order polynom"""
    return a0+a1*log_lam+a2*(log_lam)**2+a3*(log_lam)**3+a4*(log_lam)**4+a5*(log_lam)**5
# intitialize fit paramters
init_fit_param = [0, 0, 0, 0 ,0, 0]
# fit data using scipy.curve_fit()
fit_params, fit_covariances = curve_fit(fit_polynom, log_lam, flux_raw, p0=init_fit_param)
# calculate root mean square deviation
rms_array = np.sqrt((flux_raw - fit_polynom(log_lam, *fit_params))**2)
# only use flux values for 2nd fit which have a rms value < 2
fit_params2, fit_covariances2 = curve_fit(fit_polynom, log_lam[np.argwhere(rms_array<=2)][:,0], flux_raw[np.argwhere(rms_array <= 2)][:,0], p0=init_fit_param)
# remove continuum by dividing flux by fit polynomial
flux = (flux_raw / (fit_polynom(log_lam, *fit_params2))) -1

I chose a fipolynom of order 5 this can obviously be varied, but shouldnt be the main concern. This is a plot of the flux over wavelength and flux over logarithmic rebinned wavelength with the fitpolynom.
And here is a plot with the "removed" continuum, which does not seem quite right.

So can someone please tell me if a) I did the logarithmic rebinning correct and b) if my approach to the continuum removal makes sense or should be varied. I want to further work with the continuum removed spectra and extract line of sight velocity distributions (for galaxies ), but need to be sure that I understand the preprep correctly. Thank you very much

**Tags:** observational-astronomy, galaxy, data-analysis, spectra, python
> **Comments:**
> **username_1**: Why is you final spectrum $\sim 0$? The paper you refers to divides by the continuum fit, which should leave the spectrum close to 1.  Obviously if you divide (or subtract)  a small number by a (uncertain) small number at the blue wavelength edge then you could get unpredictable behaviour. The x-axis values don't make any sense. Logarithmic rebinning doesn't change the x-axis values/range it just adaptively changes the size of the bins so that each one is an equal step in $\log \lambda$.  
> **username_2**: I found a copy of the paper with selectable text and replaced your screenshot. Some folks use screen readers that convert text to audio and *images of text* don't process usually, so that is generally discouraged in Stack Exchange.  
> **username_0**: I was biased by a plot in the paper I referred to because there the flux was arround 0, that is why I subtracted 1 from the whole flux array for this plot which I just fixed. So should I only work further without the blue wave length end in this case to avoid this behaviour ? For the logarithmic rebinning I think I am lacking some understanding of the phrases used, and I am sorry if this should be obvious. What exactly is meant by bin, and size of the bin in this context ? I got an array for the flux and one for the wavelength, I thought I had to reconstruct the wavelength array ?  
> **username_0**: My guess it that a bin in this context means a particular range of wavelengths?  
> **username_3**: Your input data (the MILES template spectrum) is binned by linear wavelength: each element of the wavelength vector specifies a "bin"; for example, if lam[i] = 4800 (Angstroms), then it defines a bin containing photons with wavelengths ranging from 4800 - deltaLam/2 to 4800 + deltaLam/2, where deltaLam is the wavelength spacing (the size of the bins).  
> **username_3**: "Logarithmically rebinning" means converting your spectrum so the wavelength vector defines bins with a constant size in log(wavelength) rather than a constant size in linear wavelength, which is what you're starting with. Your log_lam is what you want, but you also need to actually rebin the *flux* (e.g., flux_raw) so that it agrees with the rebinned log wavelength bins.  
> **username_3**: You might look at the "log_rebin" function in Michele Cappellari's pPXF package (in the ppxf_util.py module).  
> **username_0**: Thank you very much for this helpful explanation, I have already started to work a bit with the ppxf package to test my work. Is there a command in python where I can see the exact definition of the log_rebin() function in the ppxf_util module ?  
> **username_3**: @trynerror It's defined in the ppxf_util.py module, which is pure Python; you can open the Python file in a text editor. (You may have to figure out where your system installs external Python packages in order to locate the ppxf package.)  
> **username_1**: Related https://astronomy.stackexchange.com/questions/14518/why-do-linear-velocity-redshifts-correspond-to-linear-pixel-shifts-when-the-spec  

---

## 0 answers

