Rocky Worlds DDT : Software

Introduction

The Rocky Worlds DDT is sharing High level Science Products (HLSPs) through MAST and at the same time we provide a few Python code examples for working with the available data.

For additional pieces of code on how to work with and analyze data from Rocky Worlds DDT, access the rocky-worlds-utils Github repository form the link below.

Access the RW GitHub repository

Instructions

In order to run the code provided in this webpage, you will need to install a few Python packages in your working environment. Here is the code to execute in a Terminal :


pip install numpy xarray matplotlib
pip install matplotlib astroquery

Once those packages are available, simply copy the codes below to your Notebook and execute the codes as usual. You should be able to replicate the plots shown here.

General Python Package Imports Package Imports



# This is the set of package imports and plotting paramters 
# needed to execute the Python codes provided in this webpage.

# Import necessary packages
import xarray as xr
import numpy as np
from astropy.stats import sigma_clip
from astropy.io import fits
from astroquery.mast import Observations

# Plotting 
import matplotlib
import matplotlib.pyplot as plt

# Set global Matplotlib parameters
matplotlib.rcParams['figure.figsize'] = (14, 6)  # Figure size in inches
matplotlib.rcParams['figure.dpi'] = 100  # DPI
matplotlib.rcParams['figure.facecolor'] = 'white'  # Figure background color
matplotlib.rcParams['figure.edgecolor'] = 'black'  # Figure border color
matplotlib.rcParams['axes.linewidth'] = 1  # Axes border width
matplotlib.rcParams['xtick.major.width'] = 1  # X-axis major tick width
matplotlib.rcParams['ytick.major.width'] = 1  # Y-axis major tick width
matplotlib.rcParams['xtick.major.size'] = 8  # X-axis major tick length
matplotlib.rcParams['ytick.major.size'] = 8  # Y-axis major tick length
matplotlib.rcParams['xtick.labelsize'] = 10  # X-axis tick label size
matplotlib.rcParams['ytick.labelsize'] = 10  # Y-axis tick label size
matplotlib.rcParams['xtick.major.pad'] = 12  # Padding for x-axis tick labels
matplotlib.rcParams['ytick.major.pad'] = 12  # Padding for y-axis tick labels
matplotlib.rcParams['xtick.color'] = 'black'  # X-axis tick color
matplotlib.rcParams['ytick.color'] = 'black'  # Y-axis tick color
matplotlib.rcParams['xtick.direction'] = 'in'  # X-axis tick direction
matplotlib.rcParams['ytick.direction'] = 'in'  # Y-axis tick direction
matplotlib.rcParams['xtick.top'] = True  # X-axis ticks on top
matplotlib.rcParams['ytick.right'] = True  # Y-axis ticks on right
matplotlib.rcParams['axes.titlesize'] = 14  # Title font size (from commented code)
matplotlib.rcParams['axes.labelsize'] = 12  # Axis label size (matching tick label size)
matplotlib.rcParams['legend.fontsize'] = 12  # Legend font size (if needed)

Accessing RW Data via MAST Astroquery


# Note that this code assumes that you want to download 
# all products from the Rocky Worlds HLSP

# Search for all ROCKY-WORLDS products
all_obs = Observations.query_criteria(provenance_name="rocky-worlds")
data_products = Observations.get_product_list(all_obs)

# Print the number of data products that would be downloaded
print('Downloading ', str(len(data_products)) , 'data products')

# Download data
Observations.download_products(data_products)

Plotting JWST Light Curves JWST Light Curves



# This code reads a given JWST light curve and plots the raw 
# and cleaned flux with fitted models.
# Note that the input H5 light curve file is needed to 
# successfully run this code.

# Load the input light curve file
filename_lc = 'hlsp_rocky-worlds_jwst_miri_gj3929b-ecl001_f1500w_v1.0_lc.h5'
ds_lc = xr.load_dataset(filename_lc)
 
# Compute the x-axis offset for tidy x-ticks
t_offset = int(np.floor(np.nanmin(ds_lc.time.values)))
 
# Setup the 3 panel figure
fig, axs = plt.subplot_mosaic('A;B;C', sharex=True, gridspec_kw={'hspace':0.075})
  
# Plot the raw light curve measurements and the fitted model
axs['A'].errorbar(ds_lc.time-t_offset, sigma_clip(ds_lc.rawFlux[0], 10), \
                  ds_lc.rawFluxErr[0], fmt='.', color='blue', alpha=0.1)
axs['A'].plot(ds_lc.time-t_offset, ds_lc.fullModel[0], '-', color='red', \
                  lw=1, label='Full Fitted Model', zorder=np.inf)
 
# Plot the systematics-removed measurements and the fitted eclipse model
axs['B'].errorbar(ds_lc.time-t_offset, ds_lc.cleanedFlux[0], ds_lc.rawFluxErr[0], \
                  fmt='.', color='blue', alpha=0.1)
axs['B'].plot(ds_lc.time-t_offset, ds_lc.astroModel[0], '-', color='red', \
                  lw=1, label='Fitted Eclipse Model', zorder=np.inf)
 
# Plot the data-model residuals in units of ppm
axs['C'].errorbar(ds_lc.time-t_offset, (ds_lc.rawFlux[0]-ds_lc.fullModel[0])*1e6, \
                  ds_lc.rawFluxErr[0]*1e6, fmt='.', color='blue', alpha=0.1)
axs['C'].axhline(0, color='navy', lw=1, zorder=np.inf)
 
# Setup the figure axes and other details
axs['A'].set_title( ds_lc.attrs['PLANET'] + f'  Fiducial Light Curve Fit')
axs['A'].set_ylabel("Raw Flux")
axs['B'].set_ylabel("Cleaned Flux")
axs['C'].set_ylabel("Residuals\n(ppm)")
axs['C'].set_xlabel(f"Time (BJD_TDB - {t_offset})")
axs['A'].legend(loc=1, frameon=False)
axs['B'].legend(loc=1, frameon=False)
fig.align_ylabels([axs['A'], axs['B'], axs['C']])

plt.show()

Raw and cleaned flux with fitted models

Plotting HST Spectra HST Spectra



# This code reads the HST stellar spectrum and plots 
# both the flux as a function of wavelength and the 
# reconstructed spectrum

# Note that the input _spec.fits file is needed to 
# successfully run this code.


# Load the HST spectrum file
filename = "hlsp_rocky-worlds_hst_stis_gj3929_g140m_v1.0_spec.fits"
observed_spectrum = fits.getdata(filename, ext=1)

# Plot full observed spectrum
fig, ax = plt.subplots(1, 1, figsize=(14, 4), dpi=100, facecolor='white' ) 

ax.plot(observed_spectrum['WAVELENGTH'],
        observed_spectrum['FLUX'],
        label='Stellar flux')

ax.axvspan(xmin=min(model_spectrum['WAVELENGTH']), \
           xmax=max(model_spectrum['WAVELENGTH']), \
           alpha=0.3, color='cyan')

ax.set_xlabel(r'Wavelength [${\rm \AA}$]')
ax.set_ylabel(r'Flux density' + '\n' + r'[erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$]')
ax.set_xlim(min(observed_spectrum['WAVELENGTH']), \
            max(observed_spectrum['WAVELENGTH']))
plt.legend(loc=2, frameon=False)
plt.show()

# Plot observed Lyman-alpha spectrum and the reconstructed model
model_spectrum = fits.getdata(filename, ext=2)
fig, ax = plt.subplots(1, 1, figsize=(14, 4), dpi=100, facecolor='white' ) 

ax.errorbar(observed_spectrum['WAVELENGTH'],
            observed_spectrum['FLUX'],
            yerr=observed_spectrum['FLUXERROR'],
            color='blue',
            label='Observed spectrum')
ax.errorbar(model_spectrum['WAVELENGTH'],
            model_spectrum['FLUX'],
            yerr=model_spectrum['FLUXERROR'],
            color='red',
            label='Reconstruction model')

plt.legend(loc=2, frameon=False)
ax.set_xlabel(r'Wavelength [${\rm \AA}$]')
ax.set_ylabel(r'Flux density' + '\n' + r'[erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$]')
ax.set_xlim(min(model_spectrum['WAVELENGTH']), max(model_spectrum['WAVELENGTH']))

plt.show()

HST Stellar Spectrum

HST Reconstructed Spectrum

Raw and cleaned flux with fitted models

HST Stellar Spectrum

HST Reconstructed Spectrum

Change Log

ABOUT

NEWS

TARGETS

STATUS

DATA

TEAMS

SUPPORT