🌈⏰🌊 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
.ipynb
file) 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/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/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.14'
🧑💻 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 2 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].
🌈🤖 The 2048 input wavelengths were not monotonically increasing. <🌈(2048w, 280t)> has been sorted from lowest to highest wavelength. If you want to recover the original wavelength order, the original wavelength indices are available in `rainbow.original_wave_index`.
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()
normalized.pcolormesh(vmin=vmin, vmax=vmax, filename="unbinned-2D-flux.png");
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)
tidied.pcolormesh(vmin=vmin, vmax=vmax, filename="trimmed-2D-flux.png");
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()
binned.pcolormesh(vmin=vmin, vmax=vmax, filename="binned-2D-flux.png");
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")
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")
<IPython.core.display.Image object>
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")
/Users/zach/Dropbox/zach/code/chromatic/chromatic/rainbows/visualizations/utilities.py:169: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored plt.scatter(x, y, **scatterkw)
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)
really_binned.plot(spacing=(vmax - vmin) / 2, filename="stack-of-lightcurves.png");
🌫 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()
)
out_of_transit.pcolormesh(
vmin=vmin, vmax=vmax, filename="trimmed-2D-flux-with-transit-removed.png"
);
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")
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")
/Users/zach/Dropbox/zach/code/chromatic/chromatic/rainbows/visualizations/wavelike/noise_comparison.py:132: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored plt.scatter(plot_x, yplot, **this_scatterkw)
/Users/zach/Dropbox/zach/code/chromatic/chromatic/rainbows/visualizations/wavelike/noise_comparison.py:132: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored plt.scatter(plot_x, yplot, **this_scatterkw)
/Users/zach/Dropbox/zach/code/chromatic/chromatic/rainbows/visualizations/wavelike/noise_comparison.py:132: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored plt.scatter(plot_x, yplot, **this_scatterkw)
/Users/zach/Dropbox/zach/code/chromatic/chromatic/rainbows/visualizations/wavelike/noise_comparison.py:132: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored plt.scatter(plot_x, yplot, **this_scatterkw)
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")
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")
/Users/zach/Dropbox/zach/code/chromatic/chromatic/rainbows/visualizations/wavelike/noise_comparison.py:132: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored plt.scatter(plot_x, yplot, **this_scatterkw) 🌈🤖 The `remove_trends` function was applied to this `Rainbow`, making it very plausible that some long-timescale signals and/or noise have been suppressed. Be suspicious of binned scatters on long timescales.
/Users/zach/Dropbox/zach/code/chromatic/chromatic/rainbows/visualizations/wavelike/noise_comparison.py:132: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored plt.scatter(plot_x, yplot, **this_scatterkw) 🌈🤖 The `remove_trends` function was applied to this `Rainbow`, making it very plausible that some long-timescale signals and/or noise have been suppressed. Be suspicious of binned scatters on long timescales.
/Users/zach/Dropbox/zach/code/chromatic/chromatic/rainbows/visualizations/wavelike/noise_comparison.py:132: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored plt.scatter(plot_x, yplot, **this_scatterkw) 🌈🤖 The `remove_trends` function was applied to this `Rainbow`, making it very plausible that some long-timescale signals and/or noise have been suppressed. Be suspicious of binned scatters on long timescales.
/Users/zach/Dropbox/zach/code/chromatic/chromatic/rainbows/visualizations/wavelike/noise_comparison.py:132: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored plt.scatter(plot_x, yplot, **this_scatterkw) 🌈🤖 The `remove_trends` function was applied to this `Rainbow`, making it very plausible that some long-timescale signals and/or noise have been suppressed. Be suspicious of binned scatters on long timescales.
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()
Hooray for you! You asked for help on what you can do with this 🌈 object. Here's a quick reference of a few available options for things to try. ----------- | actions | ----------- 🌈🧮📝 | +-*/ Do basic math operations with two Rainbows. 🌈🗂🔪 | .[:,:]() Index, slice, or mask a Rainbow to get a subset. 🌈🚧🌊 | .align_wavelengths() Align spectra with different wavelength grids onto one shared axis. 🌈🧺🧱 | .bin() Bin to a new wavelength or time grid. 🌈🧑🤝🧑🌈 | .compare() Connect to other 🌈 objects for easy comparison. 🌈🐈⏰ | .concatenate_in_time() Stitch together two Rainbows with identical wavelengths. 🌈🐈🌊 | .concatenate_in_wavelength() Stitch together two Rainbows with identical times. 🌈🚩👀 | .flag_outliers() Flag outlier data points. 🌈⏲🎞 | .fold() Fold times relative to a particular period and epoch. 🌈🧺⏰ | .get_average_lightcurve_as_rainbow() Bin down to a single integrated light curve. 🌈🧺🌊 | .get_average_spectrum_as_rainbow() Bin down to a single integrated spectrum. 🌈🎧🎲 | .inject_noise() Inject (uncorrelated, simple) random noise. 🌈🎧🎹 | .inject_systematics() Inject (correlated, wobbly) systematic noise. 🌈⭐️👻 | .inject_spectrum() Inject a static stellar spectrum. 🌈🪐🚞 | .inject_transit() Inject a transit signal. 🌈🫓😑 | .normalize() Normalize by dividing through by a typical spectrum (and/or light curve). 🌈🏴☠️🛁 | .remove_trends() Remove smooth trends in time and/or wavelength. 🌈🚇🌊 | .shift() Doppler shift wavelengths. 🌈🍱💇 | .trim() Trim away wavelengths or times. ---------------- | get/timelike | ---------------- ⏰🐋🌈 | .get_average_lightcurve() Get the weighted average light curve. ⏰🔎🌈 | .get_for_time() Get a quantity associated with a time index. ⏰🖖🌈 | .get_median_lightcurve() Get the median light curve. ⏰👌🌈 | .get_ok_data_for_time() Get a quantity associated with a time index. ⏰🛰🌈 | .get_times_as_astropy() Get the times as an astropy Time object. ⏰🚀🌈 | .set_times_from_astropy() Set the times from an astropy Time object (modifies in-place). ---------------- | get/wavelike | ---------------- 🌊🐋🌈 | .get_average_spectrum() Get the weighted average spectrum. 🌊🔎🌈 | .get_for_wavelength() Get a quantity associated with a wavelength index. 🌊🎯🌈 | .get_measured_scatter() Get the measured scatter on the time series for each wavelength. 🌊🖖🌈 | .get_median_spectrum() Get the median spectrum. 🌊👌🌈 | .get_ok_data_for_wavelength() Get a quantity associated with a wavelength index. 🌊💎🌈 | .get_spectral_resolution() Get the spectral resolution (R=w/dw). ----------- | helpers | ----------- 🙋🌈📄 | .help() Get one-line help summaries of available methods. 📓🌈🪵 | .history() Get the history that went into this Rainbow. 💾🌈📼 | .save() Save this Rainbow out to a permanent file. ------------------ | visualizations | ------------------ 🎨📽⏰ | .animate_lightcurves() Animate a sequence of light curves across different wavelengths. 🎨📽🌊 | .animate_spectra() Animate a sequence of spectra across different times. 🎨🖼📺 | .imshow() Paint a map of flux across wavelength and time. 🎨🕹📺 | .imshow_interact() Show flux map and lightcurves with interactive wavelength selection. 🎨🖌📺 | .pcolormesh() Paint a map of flux across wavelength and time (with non-uniform grids). 🎨🖌🧶 | .plot() Plot a sequence of light curves with vertical offsets. ------------------------------ | visualizations/diagnostics | ------------------------------ 🎨🗂📺 | .imshow_quantities() Show multiple 2D (wavelength and time) quantities as imshow maps. 🎨🗂🧶 | .plot_quantities() Show multiple 1D (wavelength or time) quantities as scatter plots. --------------------------- | visualizations/timelike | --------------------------- 🎨⏰🐋 | .plot_average_lightcurve() Plot the weighted average flux per time. 🎨⏰🖖 | .plot_median_lightcurve() Plot the median flux per time. --------------------------- | visualizations/wavelike | --------------------------- 🎨🌊🔭 | .plot_average_spectrum() Plot the weighted average flux per wavelength. 🎨🌊💎 | .plot_spectral_resolution() Plot the spectral resolution per wavelength. 🎨🌊🎧 | .plot_noise_comparison() Plot the measured and expected scatter per wavelength. 🎨🌊🥦 | .plot_noise_comparison_in_bins() Plot measured and expected scatter in different size bins.