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.enter image description here

And here is a plot with the "removed" continuum, which does not seem quite right. enter image description here

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

  • 1
    $\begingroup$ 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$. $\endgroup$
    – ProfRob
    Commented May 13, 2021 at 12:28
  • 1
    $\begingroup$ 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). $\endgroup$ Commented May 14, 2021 at 9:18
  • 1
    $\begingroup$ "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. $\endgroup$ Commented May 14, 2021 at 9:22
  • 1
    $\begingroup$ @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.) $\endgroup$ Commented May 14, 2021 at 22:00
  • 1
    $\begingroup$ Related astronomy.stackexchange.com/questions/14518/… $\endgroup$
    – ProfRob
    Commented Feb 19, 2023 at 7:52


You must log in to answer this question.