🌈 Quickstart¶

This page shows how to load a time-series spectroscopic dataset, do some basic calculations with it, generate some visualizations, and then save it out to another format.

In [1]:

from chromatic import read_rainbow, version
from astroquery.mast import Observations
import astropy.units as u

from chromatic import read_rainbow, version from astroquery.mast import Observations import astropy.units as u

In [2]:

version()

version()

Out[2]:

'0.4.12'

💾 Download¶

Let's download the JWST Early Release Observation of HAT-P-18b, one of first exoplanet transit datasets to be gathered by the telescope. We'll get the default pipeline x1dints (Stage 3) outputs; there are lots of reasons why we shouldn't use these particular pipeline files for science, but they're useful for a quick initial look.

In [3]:

Observations.download_file(
    f"mast:JWST/product/jw02734-o001_t001_niriss_clear-gr700xd-substrip256_x1dints.fits"
);

Observations.download_file( f"mast:JWST/product/jw02734-o001_t001_niriss_clear-gr700xd-substrip256_x1dints.fits" );

🧑‍💻 Read¶

Next, let's load that transit dataset into a Rainbow (🌈) object. These chromatic 🌈 objects keep track of how the brightness of source changes across both wavelength and time.

In [4]:

rainbow = read_rainbow(
    "jw02734-o001_t001_niriss_clear-gr700xd-substrip256_x1dints.fits"
)

rainbow = read_rainbow( "jw02734-o001_t001_niriss_clear-gr700xd-substrip256_x1dints.fits" )

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[4], line 1
----> 1 rainbow = read_rainbow(
      2     "jw02734-o001_t001_niriss_clear-gr700xd-substrip256_x1dints.fits"
      3 )

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:214, in from_x1dints(rainbow, filepath, order, **kw)
    211 filenames = expand_filenames(filepath)
    213 # extract the versions of the calibration pipeline
--> 214 cal_versions = get_versions_from_x1dints_files(filenames)
    216 # loop over file (each one a segment)
    217 for i_file, f in enumerate(filenames):
    218 
    219     # open this fits file

File ~/Dropbox/zach/code/chromatic/chromatic/rainbows/readers/x1dints.py:173, in get_versions_from_x1dints_files(filenames)
    171     versions = []
    172     for f in filenames:
--> 173         with fits.open(f) as hdu:
    174             versions.append(hdu["PRIMARY"].header["CAL_VER"])
    175 except KeyError:

File ~/opt/anaconda3/envs/2024/lib/python3.12/site-packages/astropy/io/fits/hdu/hdulist.py:222, in fitsopen(name, mode, memmap, save_backup, cache, lazy_load_hdus, ignore_missing_simple, use_fsspec, fsspec_kwargs, decompress_in_memory, **kwargs)
    219 if not name:
    220     raise ValueError(f"Empty filename: {name!r}")
--> 222 return HDUList.fromfile(
    223     name,
    224     mode,
    225     memmap,
    226     save_backup,
    227     cache,
    228     lazy_load_hdus,
    229     ignore_missing_simple,
    230     use_fsspec=use_fsspec,
    231     fsspec_kwargs=fsspec_kwargs,
    232     decompress_in_memory=decompress_in_memory,
    233     **kwargs,
    234 )

File ~/opt/anaconda3/envs/2024/lib/python3.12/site-packages/astropy/io/fits/hdu/hdulist.py:486, in HDUList.fromfile(cls, fileobj, mode, memmap, save_backup, cache, lazy_load_hdus, ignore_missing_simple, **kwargs)
    467 @classmethod
    468 def fromfile(
    469     cls,
   (...)
    477     **kwargs,
    478 ):
    479     """
    480     Creates an `HDUList` instance from a file-like object.
    481 
   (...)
    484     documentation for details of the parameters accepted by this method).
    485     """
--> 486     return cls._readfrom(
    487         fileobj=fileobj,
    488         mode=mode,
    489         memmap=memmap,
    490         save_backup=save_backup,
    491         cache=cache,
    492         ignore_missing_simple=ignore_missing_simple,
    493         lazy_load_hdus=lazy_load_hdus,
    494         **kwargs,
    495     )

File ~/opt/anaconda3/envs/2024/lib/python3.12/site-packages/astropy/io/fits/hdu/hdulist.py:1168, in HDUList._readfrom(cls, fileobj, data, mode, memmap, cache, lazy_load_hdus, ignore_missing_simple, use_fsspec, fsspec_kwargs, decompress_in_memory, **kwargs)
   1165 if fileobj is not None:
   1166     if not isinstance(fileobj, _File):
   1167         # instantiate a FITS file object (ffo)
-> 1168         fileobj = _File(
   1169             fileobj,
   1170             mode=mode,
   1171             memmap=memmap,
   1172             cache=cache,
   1173             use_fsspec=use_fsspec,
   1174             fsspec_kwargs=fsspec_kwargs,
   1175             decompress_in_memory=decompress_in_memory,
   1176         )
   1177     # The Astropy mode is determined by the _File initializer if the
   1178     # supplied mode was None
   1179     mode = fileobj.mode

File ~/opt/anaconda3/envs/2024/lib/python3.12/site-packages/astropy/io/fits/file.py:218, in _File.__init__(self, fileobj, mode, memmap, overwrite, cache, use_fsspec, fsspec_kwargs, decompress_in_memory)
    216     self._open_fileobj(fileobj, mode, overwrite)
    217 elif isinstance(fileobj, (str, bytes)):
--> 218     self._open_filename(fileobj, mode, overwrite)
    219 else:
    220     self._open_filelike(fileobj, mode, overwrite)

File ~/opt/anaconda3/envs/2024/lib/python3.12/site-packages/astropy/io/fits/file.py:651, in _File._open_filename(self, filename, mode, overwrite)
    648 ext = os.path.splitext(self.name)[1]
    650 if not self._try_read_compressed(self.name, magic, mode, ext=ext):
--> 651     self._file = open(self.name, IO_FITS_MODES[mode])
    652     self.close_on_error = True
    654 # Make certain we're back at the beginning of the file
    655 # BZ2File does not support seek when the file is open for writing, but
    656 # when opening a file for write, bz2.BZ2File always truncates anyway.

FileNotFoundError: [Errno 2] No such file or directory: 'jw02734-o001_t001_niriss_clear-gr700xd-substrip256_x1dints.fits'

The 🌈 object we just loaded provides easy access to the different dimensions we might want from the dataset, arrays like wavelength, time, flux, or uncertainty. If appropriate, quantities will have astropy Units.

In [5]:

rainbow.wavelength

rainbow.wavelength

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 1
----> 1 rainbow.wavelength

NameError: name 'rainbow' is not defined

In [6]:

rainbow.time

rainbow.time

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[6], line 1
----> 1 rainbow.time

NameError: name 'rainbow' is not defined

In [7]:

rainbow.flux

rainbow.flux

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 1
----> 1 rainbow.flux

NameError: name 'rainbow' is not defined

In [8]:

rainbow.uncertainty

rainbow.uncertainty

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 1
----> 1 rainbow.uncertainty

NameError: name 'rainbow' is not defined

🧮 Calculate¶

The absolute flux doesn't matter that much for many transit analyses, so let's use .normalize() to normalize out the median spectrum of the star, converting the data to relative brightness within each wavelength. This 🌈 action returns another 🌈 object, just with the brightness normalized.

In [9]:

normalized = rainbow.normalize()

normalized = rainbow.normalize()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[9], line 1
----> 1 normalized = rainbow.normalize()

NameError: name 'rainbow' is not defined

The dataset is really large, so for making some simple visualization it might help to average over bins of wavelength and/or time. Let's use .bin() to bin onto a (logarithmically) uniform wavelength grid, returning the binned 🌈.

In [10]:

binned = normalized.bin(R=200, dt=4 * u.minute)

binned = normalized.bin(R=200, dt=4 * u.minute)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[10], line 1
----> 1 binned = normalized.bin(R=200, dt=4 * u.minute)

NameError: name 'normalized' is not defined

The times in this dataset are measured relative to some arbitrary time in the distant past. To make them easier to interpret we can phase-fold the times so they're measured relative to the mid-transit time, when the planet is directly between the star and us, according to the planet's orbital properties from the NASA Exoplanet Archive. Let's use .fold() to change the times, returning a phase-folded 🌈.

In [11]:

folded = binned.fold(period=5.5080232 * u.day, t0=2459033.31735 * u.day)

folded = binned.fold(period=5.5080232 * u.day, t0=2459033.31735 * u.day)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[11], line 1
----> 1 folded = binned.fold(period=5.5080232 * u.day, t0=2459033.31735 * u.day)

NameError: name 'binned' is not defined

🎨 Visualize¶

We can visualize the dataset by making a map of the star's brightness across both wavelength and time, an image in which each the brightness along row corresponds to a transit light curve at that wavelength. Let's use .imshow() to create this map for the normalized, binned, folded 🌈.

In [12]:

folded.imshow();

folded.imshow();

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[12], line 1
----> 1 folded.imshow();

NameError: name 'folded' is not defined

It might be nice to look closely at the light curves within a particular wavelength range. Let's use .imshow_interact() to interactively explore the 🌈. Click and drag on the panel on the left to select the wavelegnth range to display as a light curve on the right.

In [13]:

folded.imshow_interact()

folded.imshow_interact()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[13], line 1
----> 1 folded.imshow_interact()

NameError: name 'folded' is not defined

What a dataset! It looks like there's something a little odd happening at about 2 microns (probably contamination from another star or spectrograph order) and a starspot crossing just after mid-transit, but otherwise it's a remarkably beautiful transit from a very impressive telescope!

🧶 Build¶

Because many of the actions possible with 🌈 objects return other 🌈 objects, it's possible to connect multiple steps into a single command, building up complicated analysis stories with relatively succinct code.

In [14]:

(
    rainbow.normalize()
    .flag_outliers()
    .bin(R=10)
    .fold(period=5.5080232 * u.day, t0=2459033.31735 * u.day)
    .plot()
);

( rainbow.normalize() .flag_outliers() .bin(R=10) .fold(period=5.5080232 * u.day, t0=2459033.31735 * u.day) .plot() );

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[14], line 2
      1 (
----> 2     rainbow.normalize()
      3     .flag_outliers()
      4     .bin(R=10)
      5     .fold(period=5.5080232 * u.day, t0=2459033.31735 * u.day)
      6     .plot()
      7 );

NameError: name 'rainbow' is not defined

💾 Save¶

Let's convert these data into a different format by saving it as a new file, which we might send around to share with our colleagues or publish along with a paper. chromatic can read and save 🌈 datasets with a variety of formats, to try to ease collaboration across different pipelines and toolkits.

In [15]:

rainbow.save("jwst-hatp18b.rainbow.npy")

rainbow.save("jwst-hatp18b.rainbow.npy")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[15], line 1
----> 1 rainbow.save("jwst-hatp18b.rainbow.npy")

NameError: name 'rainbow' is not defined

📚 Learn¶

That's it! This quick tutorial highlighted chromatic's abilities to...

• load in time-series spectra or multiwavelength light curves from formats like x1dints
• access core data variables like wavelength, time, flux, uncertainty
• perform calculations like .normalize, .bin, .fold
• visualize the data with .imshow, .imshow_interact, .plot

Hopefully, you're now curious to read through the User Guide to learn more about options for reading 🌈s, doing actions with 🌈s, visualizing 🌈s in different ways, and more! You can also run the .help() method associated with any 🌈 object to get a quick summary of what other methods are available for it>

In [16]:

rainbow.help()

rainbow.help()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[16], line 1
----> 1 rainbow.help()

NameError: name 'rainbow' is not defined