🌈⏰🌊 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
Once chromatic
is installed, you should have access to all the tools you need!
from chromatic import *
version()
'0.4.5'
🧑💻 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).
๐๐ค 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/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/opt/anaconda3/envs/chromatic/lib/python3.10/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, /Users/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/opt/anaconda3/envs/chromatic/lib/python3.10/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, /Users/zach/opt/anaconda3/envs/chromatic/lib/python3.10/site-packages/matplotlib/axes/_base.py:2539: UserWarning: Warning: converting a masked element to nan. xys = np.asarray(xys) /Users/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/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/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/opt/anaconda3/envs/chromatic/lib/python3.10/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, /Users/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/opt/anaconda3/envs/chromatic/lib/python3.10/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, /Users/zach/opt/anaconda3/envs/chromatic/lib/python3.10/site-packages/matplotlib/axes/_base.py:2539: UserWarning: Warning: converting a masked element to nan. xys = np.asarray(xys) /Users/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/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/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. ๐๐ฉ๐ | .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 that contain no good data. ---------------- | 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() Get the history that went into this Rainbow. ------------------ | 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.