🌈⏰🌊 time-series spectra¶
This example notebook provides quick visualizations to characterize a set of time-series stellar spectra for a transiting exoplanet, using the chromatic toolkit. This notebook can be run before you have fit any models to the data, requiring only the data themselves and a few details about the planet. To access this notebook, you might want to:
- Download it from GitHub (click 'Raw' to download the
.ipynbfile) and run on your own computer. - Open in Google Colab, load your data into that interactive session, and run everything online.
You should be able to modify just code cells in the 💻 Load the Data + 🪐 Describe the Planet section, and then run the entire notebook to automatically generate visualizations.
💾 Make Sure chromatic is Installed¶
If you don't already have chromatic installed, run the following command.
!pip install chromatic-lightcurves --upgrade --quiet
DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/virga_exo-0.4-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/bokeh-3.4.3-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/photutils-1.13.0-py3.12-macosx-11.1-arm64.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/cftime-1.6.4-py3.12-macosx-11.1-arm64.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/shapely-2.0.6-py3.12-macosx-11.1-arm64.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/netCDF4-1.7.1.post2-py3.12-macosx-11.1-arm64.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/PyMieScatt-1.8.1.1-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/bibtexparser-2.0.0b7-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/pylatexenc-3.0a30-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/platon-5.3-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/numba-0.60.0-py3.12-macosx-11.1-arm64.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/picaso-3.2.2-py3.12.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330 DEPRECATION: Loading egg at /Users/zach/opt/anaconda3/envs/2024/lib/python3.12/site-packages/llvmlite-0.43.0-py3.12-macosx-11.1-arm64.egg is deprecated. pip 24.3 will enforce this behaviour change. A possible replacement is to use pip for package installation.. Discussion can be found at https://github.com/pypa/pip/issues/12330
Once chromatic is installed, you should have access to all the tools you need!
from chromatic import *
version()
'0.4.12'
🧑💻 Load the Data + 🪐 Describe the Planet¶
Let's load the data into a chromatic 🌈 object. If you encounter errors loading your dataset, try specifying the file format with read_rainbow(filename, format=...) as described in the Reading/Writing a 🌈. Let's also associated a title with this object, which will automatically appear in most plots.
filename = "example-datasets/stsci/jw02734*x1dints.fits"
rainbow = read_rainbow(filename)
rainbow.title = "WASP-96b | NIRISS | x1dints"
🌈🤖 This file contains data for 3 spectroscopic orders. Because no `order=` keyword was supplied, we're defaulting to first order. You can hide this warning by expliciting stating which order you want to load. For this file, the options include [1 2 3].
🌈🤖 No `int_times` extension was found in the first file jw02734002001_04101_00001-seg001_nis_x1dints.fits 🌈🤖 Times were set by linearly interpolating between the exposure start and end points, as estimated from the 'SCI' extension using the 'TDB-BEG', 'TDB-END', and 'EFFINTTM' keywords. Times may be off by a few seconds and possibly up to the duration of one integration (= 76.916s).
--------------------------------------------------------------------------- KeyError Traceback (most recent call last) Cell In[4], line 2 1 filename = "example-datasets/stsci/jw02734*x1dints.fits" ----> 2 rainbow = read_rainbow(filename) 3 rainbow.title = "WASP-96b | NIRISS | x1dints" File ~/Dropbox/zach/code/chromatic/chromatic/rainbows/__init__.py:29, in read_rainbow(filepath, **kw) 8 def read_rainbow(filepath, **kw): 9 """ 10 A friendly wrapper to load time-series spectra and/or 11 multiwavelength light curves into a `chromatic` Rainbow (...) 27 The loaded data! 28 """ ---> 29 r = Rainbow(filepath, **kw) 30 if "model" in r.fluxlike: 31 return RainbowWithModel(**r._get_core_dictionaries()) File ~/Dropbox/zach/code/chromatic/chromatic/rainbows/rainbow.py:227, in Rainbow.__init__(self, filepath, format, wavelength, time, flux, uncertainty, wavelike, timelike, fluxlike, metadata, name, **kw) 225 # then try to initialize from a file 226 elif isinstance(filepath, str) or isinstance(filepath, list): --> 227 self._initialize_from_file(filepath=filepath, format=format, **kw) 229 # finally, tidy up by guessing the scales 230 self._guess_wscale() File ~/Dropbox/zach/code/chromatic/chromatic/rainbows/rainbow.py:463, in Rainbow._initialize_from_file(self, filepath, format, **kw) 461 # pick the appropriate reader 462 reader = guess_reader(filepath=filepath, format=format) --> 463 reader(self, filepath, **kw) 465 # validate that something reasonable got populated 466 self._validate_core_dictionaries() File ~/Dropbox/zach/code/chromatic/chromatic/rainbows/readers/x1dints.py:360, in from_x1dints(rainbow, filepath, order, **kw) 356 current_integration = hdu[e].header["int_num"] 358 # in case of missing segments, convert integration to index 359 current_time_index = np.nonzero( --> 360 timelike["integration_number"] == current_integration 361 )[0][0] 362 # print(current_integration, current_time_index) 363 364 # loop through all the columns in the data extension 365 for column in hdu[e].columns: 366 367 # get a lower case name for the unit KeyError: 'integration_number'
Let's define a few basic parameters describing the planet. These will help make the plots easier to interpret and be used to mask out the transit for noise characterization. If you need a period, epoch, and duration for your transit, you might find them in the NASA Exoplanet Archive.
period = 3.4252577 * u.day
t0 = 2459111.30170 * u.day
duration = 2.4264 * u.hour * 1.1
Let's set some very basic plotting defaults. These will try to sync up the vertical limits of light curve plots with the color maps of 2D flux maps.
vmin = 0.98
vmax = 1.005
ylim = [vmin, vmax]
🎨 Make Basic Visualizations¶
Let's make some plots to get an overall sense for the dataset. First of all, let's normalize by dividing through by the median stellar spectrum and display a 2D map of the flux.
normalized = rainbow.normalize()
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[7], line 1 ----> 1 normalized = rainbow.normalize() AttributeError: module 'chromatic.rainbows.rainbow' has no attribute 'normalize'
normalized.pcolormesh(vmin=vmin, vmax=vmax, filename="unbinned-2D-flux.png");
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[8], line 1 ----> 1 normalized.pcolormesh(vmin=vmin, vmax=vmax, filename="unbinned-2D-flux.png"); NameError: name 'normalized' is not defined
Let's trim any bad wavelengths off the edges and phase-fold these data to the planet's known period and time of mid-transit. This latter step is just to make the times easier to interpret.
tidied = normalized.trim().fold(period=period, t0=t0)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[9], line 1 ----> 1 tidied = normalized.trim().fold(period=period, t0=t0) NameError: name 'normalized' is not defined
tidied.pcolormesh(vmin=vmin, vmax=vmax, filename="trimmed-2D-flux.png");
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[10], line 1 ----> 1 tidied.pcolormesh(vmin=vmin, vmax=vmax, filename="trimmed-2D-flux.png"); NameError: name 'tidied' is not defined
Now let's bin to a fixed wavelength resolution $R=\lambda/\Delta\lambda$ and cadence $dt$. Be averaging together multiple wavelengths and/or times, we'll decrease the noise (at the cost of worse resolution); some features may become more apparent with stronger binning, as long as the bins aren't so large to smooth features away. We'll normalize again after binning, once the noise has been averaged down a little bit.
binned = tidied.bin(R=100, dt=2 * u.minute).normalize()
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[11], line 1 ----> 1 binned = tidied.bin(R=100, dt=2 * u.minute).normalize() NameError: name 'tidied' is not defined
binned.pcolormesh(vmin=vmin, vmax=vmax, filename="binned-2D-flux.png");
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[12], line 1 ----> 1 binned.pcolormesh(vmin=vmin, vmax=vmax, filename="binned-2D-flux.png"); NameError: name 'binned' is not defined
If we want to see how these flux maps line up to light curves, we can use the interactive exploration tool to drag and select different wavelength ranges, and plot the (unweighted) average light curves for them. Being able to select particular wavelength ranges allows us to zoom in on features of interest or trends with wavelength.
binned.imshow_interact(ylim=ylim, filename="interactive-2D-flux.html")
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[13], line 1 ----> 1 binned.imshow_interact(ylim=ylim, filename="interactive-2D-flux.html") NameError: name 'binned' is not defined
Let's make an animation the flips through wavelength bins. Animations can be useful way to to let your eyes take in a lot of light curves in quick succession, to see how trends or noise sources change with wavelength.
binned.animate_lightcurves(ylim=ylim, filename="animated-lightcurves.gif")
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[14], line 1 ----> 1 binned.animate_lightcurves(ylim=ylim, filename="animated-lightcurves.gif") NameError: name 'binned' is not defined
Let's look at the weighted average light curve, integrating together all the available wavelengths. The binned uncertainties should be smaller than any individual wavelength bin, so this should hopefully be a very precise light curve!
binned.plot_average_lightcurve(ylim=ylim, filename="integrated-lightcurve.png")
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[15], line 1 ----> 1 binned.plot_average_lightcurve(ylim=ylim, filename="integrated-lightcurve.png") NameError: name 'binned' is not defined
Finally, let's bin the data to very low spectral resolution, and then plot a stack of light curves for the separate wavelengths.
really_binned = binned.bin(R=5)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[16], line 1 ----> 1 really_binned = binned.bin(R=5) NameError: name 'binned' is not defined
really_binned.plot(spacing=(vmax - vmin) / 2, filename="stack-of-lightcurves.png");
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[17], line 1 ----> 1 really_binned.plot(spacing=(vmax - vmin) / 2, filename="stack-of-lightcurves.png"); NameError: name 'really_binned' is not defined
🌫 Characterize the Noise¶
After looking closely at the transit in the data, it can be useful to filter or mask out the transit. By creating a dataset that we expect to be mostly flat, we can characterize the noise properties. First, let's simply mask out the transit and look only at the data before and after it.
out_of_transit = (
rainbow.trim().mask_transit(period=period, t0=t0, duration=duration).normalize()
)
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[18], line 2 1 out_of_transit = ( ----> 2 rainbow.trim().mask_transit(period=period, t0=t0, duration=duration).normalize() 3 ) AttributeError: module 'chromatic.rainbows.rainbow' has no attribute 'trim'
out_of_transit.pcolormesh(
vmin=vmin, vmax=vmax, filename="trimmed-2D-flux-with-transit-removed.png"
);
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[19], line 1 ----> 1 out_of_transit.pcolormesh( 2 vmin=vmin, vmax=vmax, filename="trimmed-2D-flux-with-transit-removed.png" 3 ); NameError: name 'out_of_transit' is not defined
Because unbinned data have the small wavelength and time intervals, and therefore contain the relatively few photons per bin, their expected uncertainties will be large. As we bin to larger wavelength intervals (lower resolution) and longer time intervals, the expected uncertainties will decrease, revealing increasingly subtle features.
resolutions = [None, 30, 10, 3]
if rainbow.ntime > 1e4:
dt = 1 * u.minute
else:
dt = None
binned_rainbows = {}
N = len(resolutions)
fi, ax = plt.subplots(N, 1, sharex=True, figsize=(8, N * 3))
for i, R in enumerate(resolutions):
binned_rainbows[R] = out_of_transit.bin(R=R, dt=dt).normalize()
binned_rainbows[R].pcolormesh(ax=ax[i])
plt.title(f"R={R}, dt={binned_rainbows[R].dt:.2} | {ax[i].get_title()}")
rainbow.savefig("binned-2D-flux-no-transit.png")
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[20], line 2 1 resolutions = [None, 30, 10, 3] ----> 2 if rainbow.ntime > 1e4: 3 dt = 1 * u.minute 4 else: AttributeError: module 'chromatic.rainbows.rainbow' has no attribute 'ntime'
If the data were perfect and simple, we would expect the scatter in the out-of-transit flux to match the expected uncertainties for each wavelength. Let's compare the expected and measured scatters, after binning down in wavelength and time. If the measured scatter does not decrease as expected when averaging bins together, there is probably some systematic noise that is correlated across wavelength and/or time that should be addressed.
fi, ax = plt.subplots(
N, 2, sharey=True, figsize=(8, N * 3), gridspec_kw=dict(width_ratios=[2, 1])
)
for i, R in enumerate(binned_rainbows):
binned_rainbows[R].plot_noise_comparison(ax=ax[i, 0], method="MAD")
plt.title(
f"R={R}, dt={binned_rainbows[R].dt:.2} | with transit removed\n{ax[i,0].get_title()}"
)
binned_rainbows[R].plot_noise_comparison_in_bins(ax=ax[i, 1], method="MAD")
rainbow.savefig("noise-comparison-no-transit.png")
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[21], line 2 1 fi, ax = plt.subplots( ----> 2 N, 2, sharey=True, figsize=(8, N * 3), gridspec_kw=dict(width_ratios=[2, 1]) 3 ) 4 for i, R in enumerate(binned_rainbows): 5 binned_rainbows[R].plot_noise_comparison(ax=ax[i, 0], method="MAD") NameError: name 'N' is not defined
It's possible that excess scatter might be the result of very smooth trends in time. Let's repeat the above analysis of the dataset binned to different wavelength and time intervals, but also remove a smooth quadratic trend in time, which could imagine might be either instrumental or astrophysical.
fi, ax = plt.subplots(N, 1, sharex=True, figsize=(8, N * 3))
for i, R in enumerate(binned_rainbows):
with_trend_removed = binned_rainbows[R].remove_trends(method="polyfit", deg=2)
with_trend_removed.pcolormesh(ax=ax[i])
plt.title(f"R={R}, dt={binned_rainbows[R].dt:.2} | {ax[i].get_title()}")
rainbow.savefig("binned-2D-flux-no-transit-no-trends.png")
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[22], line 1 ----> 1 fi, ax = plt.subplots(N, 1, sharex=True, figsize=(8, N * 3)) 2 for i, R in enumerate(binned_rainbows): 3 with_trend_removed = binned_rainbows[R].remove_trends(method="polyfit", deg=2) NameError: name 'N' is not defined
fi, ax = plt.subplots(
N, 2, sharey=True, figsize=(8, N * 3), gridspec_kw=dict(width_ratios=[2, 1])
)
for i, R in enumerate(binned_rainbows):
with_trend_removed = binned_rainbows[R].remove_trends(method="polyfit", deg=2)
with_trend_removed.plot_noise_comparison(ax=ax[i, 0], method="MAD")
plt.title(
f"R={R}, dt={binned_rainbows[R].dt:.2} | with transit + trends removed\n{ax[i,0].get_title()}"
)
with_trend_removed.plot_noise_comparison_in_bins(ax=ax[i, 1], method="MAD")
rainbow.savefig("noise-comparison-no-transit-no-trends.png")
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[23], line 2 1 fi, ax = plt.subplots( ----> 2 N, 2, sharey=True, figsize=(8, N * 3), gridspec_kw=dict(width_ratios=[2, 1]) 3 ) 4 for i, R in enumerate(binned_rainbows): 5 with_trend_removed = binned_rainbows[R].remove_trends(method="polyfit", deg=2) NameError: name 'N' is not defined
That's it! Hopefully these automated visualizations can serve as a useful starting place for understanding the dataset you're working with and for comparing to other analyses.
🗺 Explore Further!¶
You can explore other options to visualize or work with your 🌈 data by using the built-in .help() method. On any Rainbow object, run this to list available options.
rainbow.help()
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[24], line 1 ----> 1 rainbow.help() AttributeError: module 'chromatic.rainbows.rainbow' has no attribute 'help'