Masking and flattening a spectrum¶
Often in high resolution spectroscopy we wish to designate the continuum—regions in the spectrum devoid of spectral lines. The continuum is conspicuous in hot stars, but becomes ill-defined in cool M-dwarfs. In this tutorial we show how to define the continuum using muler and specutils.
[1]:
from muler.hpf import HPFSpectrum
import numpy as np
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 1
----> 1 from muler.hpf import HPFSpectrum
2 import numpy as np
File ~/checkouts/readthedocs.org/user_builds/muler/envs/latest/lib/python3.8/site-packages/muler/hpf.py:14
12 import warnings
13 import logging
---> 14 from muler.echelle import EchelleSpectrum, EchelleSpectrumList
15 from muler.utilities import resample_list
16 import numpy as np
File ~/checkouts/readthedocs.org/user_builds/muler/envs/latest/lib/python3.8/site-packages/muler/echelle.py:31
29 from scipy.ndimage import median_filter, gaussian_filter1d
30 import specutils
---> 31 from muler.utilities import apply_numpy_mask, is_list, resample_list
34 # from barycorrpy import get_BC_vel
35 from astropy.coordinates import SkyCoord, EarthLocation
File ~/checkouts/readthedocs.org/user_builds/muler/envs/latest/lib/python3.8/site-packages/muler/utilities.py:15
13 from astropy.convolution import convolve, Gaussian1DKernel
14 from scipy.ndimage import binary_dilation
---> 15 from astroquery.simbad import Simbad
16 Simbad.add_votable_fields('V', 'B', 'J', 'H', 'K', 'parallax')
17 LinInterpResampler = LinearInterpolatedResampler()
ModuleNotFoundError: No module named 'astroquery'
[2]:
import matplotlib.pyplot as plt
%config InlineBackend.figure_format='retina'
plt.rcParams['axes.facecolor']='white'
plt.rcParams['savefig.facecolor']='white'
[3]:
path = 'https://github.com/OttoStruve/muler_example_data/raw/main/HPF/01_A0V_standards/'
filename = 'Goldilocks_20210801T083618_v1.0_0036.spectra.fits'
raw_spectrum = HPFSpectrum(file = path+filename, order=6)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 3
1 path = 'https://github.com/OttoStruve/muler_example_data/raw/main/HPF/01_A0V_standards/'
2 filename = 'Goldilocks_20210801T083618_v1.0_0036.spectra.fits'
----> 3 raw_spectrum = HPFSpectrum(file = path+filename, order=6)
NameError: name 'HPFSpectrum' is not defined
[4]:
spectrum = raw_spectrum.sky_subtract(method='vector').remove_nans().deblaze().normalize()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 1
----> 1 spectrum = raw_spectrum.sky_subtract(method='vector').remove_nans().deblaze().normalize()
NameError: name 'raw_spectrum' is not defined
[5]:
spectrum.plot(ylo=0);
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 1
----> 1 spectrum.plot(ylo=0);
NameError: name 'spectrum' is not defined
Look at those conspicuous sharp spectral lines and the smooth dip near 10940 Angstroms. Let’s divide the spectrum by a tuned model.
[6]:
from gollum.phoenix import PHOENIXSpectrum
[7]:
native_template = PHOENIXSpectrum(teff=9600, logg=5.0, download=True).normalize()
[8]:
doctored_template = native_template\
.rotationally_broaden(105.0)\
.rv_shift(-40.0)\
.resample(spectrum)\
.normalize()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 4
1 doctored_template = native_template\
2 .rotationally_broaden(105.0)\
3 .rv_shift(-40.0)\
----> 4 .resample(spectrum)\
5 .normalize()
NameError: name 'spectrum' is not defined
[9]:
ax = spectrum.plot(ylo=0, label='A0V Star');
doctored_template.plot(ax=ax, label='Doctored A0V Template'); ax.legend();
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 1
----> 1 ax = spectrum.plot(ylo=0, label='A0V Star');
2 doctored_template.plot(ax=ax, label='Doctored A0V Template'); ax.legend();
NameError: name 'spectrum' is not defined
[10]:
telluric_response = spectrum / doctored_template
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 1
----> 1 telluric_response = spectrum / doctored_template
NameError: name 'spectrum' is not defined
Normalize by the top 10th percentile, rather than the median…
[11]:
spectrum_peak = np.percentile(telluric_response.flux, 65)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[11], line 1
----> 1 spectrum_peak = np.percentile(telluric_response.flux, 65)
NameError: name 'np' is not defined
[12]:
telluric_response = telluric_response.divide(spectrum_peak)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 1
----> 1 telluric_response = telluric_response.divide(spectrum_peak)
NameError: name 'telluric_response' is not defined
[13]:
ax = telluric_response.plot(ylo=0);
ax.axhline(1.0, color='k', linestyle='dotted');
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 1
----> 1 ax = telluric_response.plot(ylo=0);
2 ax.axhline(1.0, color='k', linestyle='dotted');
NameError: name 'telluric_response' is not defined
The continuum should all hover near \(1.0\). You can see that the spectrum has a slight downward curvature to it, with the blue-most values about \(5\%\) lower than \(1.0\), and some red-end values exceeding \(1.0\) by 2 percentage points.
Let’s mask all the lines, fit a smooth function to the continuum, and then flatten by that function, evaluated at all wavelength points.
Masking the lines and fitting the continuum¶
[14]:
import astropy.units as u
from specutils import SpectralRegion
from specutils.manipulation import extract_region
from specutils.fitting import find_lines_derivative, fit_continuum
[15]:
continuum_subtracted = telluric_response-1.0
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 1
----> 1 continuum_subtracted = telluric_response-1.0
NameError: name 'telluric_response' is not defined
Here we find all the absorption lines above some threshold.
[16]:
lines = find_lines_derivative(continuum_subtracted, flux_threshold=0.05)
lines = lines[lines['line_type'] == 'absorption']
len(lines)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 1
----> 1 lines = find_lines_derivative(continuum_subtracted, flux_threshold=0.05)
2 lines = lines[lines['line_type'] == 'absorption']
3 len(lines)
NameError: name 'continuum_subtracted' is not defined
[17]:
ax = telluric_response.plot(ylo=0);
ax.axhline(1.0, color='k', linestyle='dashed');
for line in lines:
ax.axvline(line['line_center'].value, linestyle='dotted', color='k', alpha=0.5)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 1
----> 1 ax = telluric_response.plot(ylo=0);
2 ax.axhline(1.0, color='k', linestyle='dashed');
4 for line in lines:
NameError: name 'telluric_response' is not defined
Neat, we see that we have found the most conspicuous lines. Let’s mask these lines with a \(\pm1\;\)A width.
[18]:
line_regions = SpectralRegion.from_line_list(lines, width=2.0*u.Angstrom)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 line_regions = SpectralRegion.from_line_list(lines, width=2.0*u.Angstrom)
NameError: name 'lines' is not defined
Fit the continuum while excluding the spectral lines.
[19]:
%%capture
continuum_model = fit_continuum(telluric_response, exclude_regions=line_regions)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 1
----> 1 continuum_model = fit_continuum(telluric_response, exclude_regions=line_regions)
NameError: name 'telluric_response' is not defined
[20]:
continuum_fit = continuum_model(telluric_response.wavelength)
continuum_spectrum = HPFSpectrum(flux=continuum_fit,
spectral_axis = telluric_response.wavelength)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 1
----> 1 continuum_fit = continuum_model(telluric_response.wavelength)
2 continuum_spectrum = HPFSpectrum(flux=continuum_fit,
3 spectral_axis = telluric_response.wavelength)
NameError: name 'continuum_model' is not defined
[21]:
ax=telluric_response.plot(ylo=0, label='Telluric response')
ax.axhline(1.0, color='k', linestyle='dashed');
continuum_spectrum.plot(ax=ax, label='Continuum fit');
ax.legend();
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[21], line 1
----> 1 ax=telluric_response.plot(ylo=0, label='Telluric response')
2 ax.axhline(1.0, color='k', linestyle='dashed');
3 continuum_spectrum.plot(ax=ax, label='Continuum fit');
NameError: name 'telluric_response' is not defined
Nice! That looks OK. You can see that we are fitting to weak telluric lines at 10910 A, which pulls down the fit from 10890. So it’s not perfect.
[22]:
flattened_telluric = telluric_response/continuum_spectrum
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[22], line 1
----> 1 flattened_telluric = telluric_response/continuum_spectrum
NameError: name 'telluric_response' is not defined
[23]:
ax=flattened_telluric.plot(ylo=0)
ax.axhline(1.0, color='k', linestyle='dashed');
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[23], line 1
----> 1 ax=flattened_telluric.plot(ylo=0)
2 ax.axhline(1.0, color='k', linestyle='dashed');
NameError: name 'flattened_telluric' is not defined
Hooray! We have flattened the continuum! :) In practice this process takes some more fine-tuning of the model complexity, etc, and one may prefer a strategy of iteratively fitting and masking. Ideally we would use the GP smoothing technique, which is possible.