Rocky Worlds DDT

code-tutorials

• Reading the JWST Eclipse Catalog

  This example shows how to open the JWST eclipse catalog file in Python and how to read metadata from netCDF4/HDF5 files.

  import xarray as xr
  
  # Open Eclipse Catalog
  filename = "hlsp_rocky-worlds_jwst_miri_gj3929b-ecl001_f1500w_v1.0.1_eclipse-cat.h5"
  ds = xr.load_dataset(filename)
  
  # Prints out summary information about the dataset
  print(ds)
  
  # An example of how to get access to the metadata stored in the file
  print(ds.attrs['HLSPTARG'])

• Plotting the JWST Light Curve

  This example opens the JWST light curve file and creates a plot of the flux over time, highlighting the best fit model results for this eclipse:

  import xarray as xr
  import numpy as np
  import matplotlib.pyplot as plt
  from astropy.stats import sigma_clip
   
  filename_lc = 'hlsp_rocky-worlds_jwst_miri_gj3929b-ecl001_f1500w_v1.0.1_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', figsize=(10*0.8,7.5*0.8), 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='.', c='k', alpha=0.1)
  axs['A'].plot(ds_lc.time-t_offset, ds_lc.fullModel[0], '-', c='r', 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='.', c='k', alpha=0.1)
  axs['B'].plot(ds_lc.time-t_offset, ds_lc.astroModel[0], '-', c='r', 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='.', c='k', alpha=0.1)
  axs['C'].axhline(0, c='r', lw=1, zorder=np.inf)
   
  # Setup the figure axes and other details
  axs['A'].set_title(f"{ds_lc.attrs['PLANET']} Fiducial Light Curve Fit")
  axs['A'].set_ylabel("Raw Flux")
  axs['B'].set_ylabel("Cleaned Flux")
  axs['C'].set_ylabel("Residuals (ppm)")
  axs['C'].set_xlabel(f"Time (BJD_TDB - {t_offset})")
  axs['A'].legend(loc=1)
  axs['B'].legend(loc=1)
  fig.align_ylabels([axs['A'], axs['B'], axs['C']])
   
  plt.show()

• Plotting the HST Spectrum

  This example opens the HST spectrum file and creates a plot for GJ 3929:

  import matplotlib.pyplot as plt
  from astropy.io import fits
  
  filename = "hlsp_rocky-worlds_hst_cos-stis_gj3929_multi_v1.1_spec.fits"
  
  observed_spectrum = fits.getdata(filename, ext=1)
  
  # Plot full observed spectrum
  ax = plt.subplot()
  ax.plot(observed_spectrum['WAVELENGTH'],
          observed_spectrum['FLUX'])
  ax.set_xlabel(r'Wavelength [${\rm \AA}$]')
  ax.set_ylabel(r'Flux density [erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$]')
  ax.set_xlim(min(observed_spectrum['WAVELENGTH']), max(observed_spectrum['WAVELENGTH']))
  plt.show()
  
  # Plot observed Lyman-alpha spectrum and the reconstruction model
  model_spectrum = fits.getdata(filename, ext=2)
  ax = plt.subplot()
  ax.errorbar(observed_spectrum['WAVELENGTH'],
              observed_spectrum['FLUX'],
              yerr=observed_spectrum['FLUXERROR'],
              label='Observed spectrum')
  ax.errorbar(model_spectrum['WAVELENGTH'],
              model_spectrum['FLUX'],
              yerr=model_spectrum['FLUXERROR'],
              label='Reconstruction model')
  
  plt.legend(loc=2)
  ax.set_xlabel(r'Wavelength [${\rm \AA}$]')
  ax.set_ylabel(r'Flux density [erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$]')
  ax.set_xlim(min(model_spectrum['WAVELENGTH']), max(model_spectrum['WAVELENGTH']))
  plt.show()

• Using the JWST Checkpoint Files

  The JWST checkpoint files re-use the reductions from the individual eclipses and contain the results from a simultaneous fit to all four eclipses. This example shows how to open the checkpoint eclipse catalog and  plot the combined light curve for GJ 3929 b.

  import matplotlib.pyplot as plt
  import numpy as np
  from scipy.stats import binned_statistic
  import xarray as xr
  
  # Load the checkpoint eclipse catalog and grouped light-curve product.
  catalog = xr.load_dataset(
      "hlsp_rocky-worlds_jwst_miri_gj3929b-checkpoint02_f1500w_v1.0_eclipse-cat.h5",
      engine="netcdf4",
  )
  lightcurves = xr.open_datatree(
      "hlsp_rocky-worlds_jwst_miri_gj3929b-checkpoint02_f1500w_v1.0_lc.h5",
      engine="netcdf4",
  )
  
  # Choose the number of time bins.
  n_bins = 28
  
  # Collect all visit light curves after shifting each one to its fitted
  # mid-eclipse time. Each visit lives in its own DataTree group.
  rows = []
  for i, name in enumerate(sorted(lightcurves.children)):
      ds = lightcurves[name].ds.isel(visit=0)
      rows.append(np.column_stack((
          (ds.time.values - float(catalog.eclipseTime.isel(visit=i))) * 24.0,
          ds.cleanedFlux.values,
          ds.rawFluxErr.values,
          ds.astroModel.values,
      )))
  
  # Stack all visits, then keep only points that were included in the fit.
  # The model is NaN for integrations that were dropped before fitting.
  time_hr, flux, err, model = np.concatenate(rows).T
  fit = np.isfinite(model) & np.isfinite(time_hr + flux + err) & (err > 0)
  order = np.argsort(time_hr[fit])
  time_hr, flux, err, model = (x[fit][order] for x in (time_hr, flux, err, model))
  
  # Bin the light curve in time using inverse-variance weights from rawFluxErr.
  edges = np.linspace(time_hr.min(), time_hr.max(), n_bins + 1)
  centers = 0.5 * (edges[:-1] + edges[1:])
  weights = 1.0 / err**2
  wsum, _, _ = binned_statistic(time_hr, weights, statistic="sum", bins=edges)
  fwsum, _, _ = binned_statistic(
      time_hr, flux * weights, statistic="sum", bins=edges
  )
  bflux = fwsum / wsum
  berr = np.sqrt(1.0 / wsum)
  
  # Plot unbinned points faintly, then overlay the model and binned light curve.
  fig, ax = plt.subplots(figsize=(7.2, 3.0))
  ax.errorbar(centers, (bflux - 1) * 1e6, yerr=berr * 1e6, fmt=".", color="k", label="Binned data")
  ylim = ax.get_ylim()
  ax.plot(time_hr, (flux - 1) * 1e6, ".", ms=3, color="0.7", alpha=0.08)
  ax.plot(time_hr, (model - 1) * 1e6, color="r", label="Eclipse model")
  ax.set(
      xlabel="Time since mid-eclipse (hours)",
      ylabel="Planetary flux (ppm)",
      ylim=ylim,
  )
  ax.legend()
  fig.tight_layout()
  fig.savefig('hlsp_rocky-worlds_jwst_miri_gj3929b-checkpoint02_f1500w_v1.0_preview.png', dpi=300, bbox_inches='tight')
  plt.show(fig)