diff --git a/.gitignore b/.gitignore index 5fb8486d..3f7c7834 100644 --- a/.gitignore +++ b/.gitignore @@ -17,6 +17,7 @@ downloaded_hdf5_files/ *.ecsv *.fits *.gzip +*.h5 *.jpg *.jpeg *.log @@ -53,3 +54,8 @@ _build/ #pycache __pycache__ + +# Local working notes, not part of the tutorials +CLAUDE.md +OU_PR_REVIEW.md +PR_DESCRIPTION.md diff --git a/deploy_to_fornax_manifest.in b/deploy_to_fornax_manifest.in index 547bfad9..57db50d2 100644 --- a/deploy_to_fornax_manifest.in +++ b/deploy_to_fornax_manifest.in @@ -8,11 +8,9 @@ tutorials/euclid/Euclid_ERO.md tutorials/euclid/euclid-cloud-access.md tutorials/euclid/merged-objects-hats-catalog/1-euclid-q1-hats-intro.md tutorials/euclid/merged-objects-hats-catalog/4-euclid-q1-hats-magnitudes.md -tutorials/simulated-data/OpenUniverse2024/OpenUniverse2024Preview_Firefly.md +tutorials/simulated-data/OpenUniverse2024/openuniverse2024_firefly_visualization.md tutorials/simulated-data/cosmoDC2_TAP_access.md -tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_demodata_20000808_instances.csv -tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_timedomainsurvey.md -tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_wideareasurvey.md +tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_tds_supernova.md tutorials/simulated-data/OpenUniverse2024/openuniverse2024_quickstart.md tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md diff --git a/ignore_tutorials/ignore_circleci_testing b/ignore_tutorials/ignore_circleci_testing index 54a81f21..036d747f 100644 --- a/ignore_tutorials/ignore_circleci_testing +++ b/ignore_tutorials/ignore_circleci_testing @@ -1,8 +1,7 @@ # ignore these as they run out of resources on CircleCI. # Note that we don't duplicate the listing of notebooks that are ignored for all types of testing. # -tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_timedomainsurvey -tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_wideareasurvey +tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_tds_supernova tutorials/simulated-data/roman_hlss_number_density tutorials/euclid/euclid-cloud-access tutorials/techniques-and-tools/cloud-access-intro diff --git a/notebook_metadata.yml b/notebook_metadata.yml index bbf963af..12a8f26f 100644 --- a/notebook_metadata.yml +++ b/notebook_metadata.yml @@ -42,34 +42,30 @@ file: tutorials/euclid/merged-objects-hats-catalog/4-euclid-q1-hats-magnitudes.md section: Euclid description: Review the types of flux measurements available, load template-fit and aperture magnitudes, and plot distributions and comparisons for different object types. -- title: OpenUniverse2024 Visualization - file: tutorials/simulated-data/OpenUniverse2024/OpenUniverse2024Preview_Firefly.md +- title: Firefly Visualization + file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_firefly_visualization.md section: Simulated Data - description: Use Firefly to get an overview of survey structure and visualize content. + description: Use Firefly to explore the structure of the OpenUniverse2024 Roman and Rubin simulations and visualize their content. - title: CosmoDC2 Data Access file: tutorials/simulated-data/cosmoDC2_TAP_access.md section: Simulated Data description: Access, query, and visualize the CosmoDC2 catalog. -- title: OpenUniverse2024 Quickstart +- title: Quickstart file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_quickstart.md section: Simulated Data description: Access OpenUniverse2024 Roman and Rubin images and catalogs. -- title: OpenUniverse2024 TDE Light Curve +- title: TDE Light Curve file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md section: Simulated Data - description: Locate a simulated TDE, retrieve Roman images, and build a multi-epoch light curve. -- title: OpenUniverse2024 SED Fitting + description: Locate a simulated TDE in OpenUniverse2024, retrieve Roman images, and build a multi-epoch light curve. +- title: SED Fitting file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md section: Simulated Data - description: Fit spectral energy distributions for supernova host galaxies using Prospector. -- title: OpenUniverse2024 Time Domain - file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_timedomainsurvey.md + description: Fit spectral energy distributions for OpenUniverse2024 supernova host galaxies using Prospector. +- title: Roman TDS Supernova + file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_tds_supernova.md section: Simulated Data - description: Access and analyze the simulated time-domain OpenUniverse2024 survey. -- title: OpenUniverse2024 Roman Coadds - file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_wideareasurvey.md - section: Simulated Data - description: Access OpenUniverse2024 wide-area simulated survey data. + description: Explore the OpenUniverse2024 Roman Time Domain Survey observation sequence and animate a simulated SN Ia across epochs. - title: Roman HLSS Number Density file: tutorials/simulated-data/roman_hlss_number_density.md section: Simulated Data diff --git a/toc.yml b/toc.yml index c214df18..969e3887 100644 --- a/toc.yml +++ b/toc.yml @@ -67,22 +67,22 @@ project: - title: Simulated Data file: tutorials/simulated-data/simulated.md children: + - title: OpenUniverse2024 + children: + - title: Quickstart + file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_quickstart.md + - title: Firefly Visualization + file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_firefly_visualization.md + - title: Roman TDS Supernova + file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_tds_supernova.md + - title: TDE Light Curve + file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md + - title: SED Fitting + file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md - title: Roman HLSS Number Density file: tutorials/simulated-data/roman_hlss_number_density.md - - title: OpenUniverse2024 Roman Coadds - file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_wideareasurvey.md - - title: OpenUniverse2024 Visualization - file: tutorials/simulated-data/OpenUniverse2024/OpenUniverse2024Preview_Firefly.md - - title: OpenUniverse2024 Time Domain - file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_timedomainsurvey.md - title: CosmoDC2 Data Access file: tutorials/simulated-data/cosmoDC2_TAP_access.md - - title: OpenUniverse2024 Quickstart - file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_quickstart.md - - title: OpenUniverse2024 TDE Light Curve - file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md - - title: OpenUniverse2024 SED Fitting - file: tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md - title: Techniques & Tools file: tutorials/techniques-and-tools/techniques.md children: diff --git a/tutorials/simulated-data/OpenUniverse2024/OpenUniverse2024Preview_Firefly.md b/tutorials/simulated-data/OpenUniverse2024/OpenUniverse2024Preview_Firefly.md deleted file mode 100644 index 52742549..00000000 --- a/tutorials/simulated-data/OpenUniverse2024/OpenUniverse2024Preview_Firefly.md +++ /dev/null @@ -1,732 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.3 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 -authors: - - name: Jaladh Singhal - - name: Vandana Desai - - name: IRSA Team ---- - -# Using Firefly to Explore OpenUniverse2024 Data Preview Simulated Roman and Rubin Images - -+++ - -## Learning Goals - -By the end of this tutorial, you will: - -- Learn how to access cloud-hosted Roman and Rubin simulated images. - -- Learn how to launch an interactive Firefly instance inside JupyterLab. - -- Learn how to use the Firefly Jupyterlab extension to visualize cloud-hosted simulated images, overplot ds9 regions, overplot catalogs in Parquet format, and create 3 color images. - -+++ - -## Introduction - -The purpose of this tutorial is to become familiar with simulated Roman and Rubin simulated data published through the OpenUniverse 2024 data preview, and to become familiar with the Firefly JupyterLab Extension for visualizing astronomical data products. - -OpenUniverse2024 is a project to simulate spatially overlapping imaging surveys to be carried out by the Nancy Grace Roman Telescope and the Vera C. Rubin Observatory. The simulations were carried out on Argonne's Theta cluster and consist of: - -- The LSST ELAIS-S1 Deep Drilling Field (DDF) -- The Roman Time-Domain Survey (TDS) shifted to overlap the ELAIS region and LSST DDF -- Overlapping LSST Wide-Fast-Deep (WFD) survey (with rolling cadence) -- Overlapping Roman Wide-Area Survey (WAS) in the same region -- A deep-field calibration region of the Roman WAS in the same region - -This data preview release consists of a subset of data from each of the five categories above. More information about data preview can be found at [IRSA holding of this dataset](https://irsa.ipac.caltech.edu/data/theory/openuniverse2024/overview.html). - -Firefly is an open-source web-based UI library for astronomical data archive access and visualization developed at Caltech and used by multiple space- and ground-based astrophysics archives. More information on Firefly can be found [here](https://github.com/Caltech-IPAC/firefly/blob/dev/README.md). - -In addition to being used to make web applications, Firefly can be used from Python. More information on Firefly Python client can be found [here](https://caltech-ipac.github.io/firefly_client/usage/index.html). - -The Firefly JupyterLab Extension makes it particularly easy to use Firefly to efficiently visualize cloud-hosted astronomical data using JupyterLab instances running locally or on cloud. More information on Firefly JupyterLab Extension can be found [here](https://github.com/Caltech-IPAC/jupyter_firefly_extensions/blob/master/README.md). - -+++ - -## Imports - -- astropy.io.fits for accessing FITS files -- numpy for numerical computing -- s3fs for browsing cloud buckets -- matplotlib.pyplot for creating static visualizations of FITS images -- matplotlib.patches for annotating visualizations of FITS images -- astropy.wcs for dealing with astronomical world coordinate systems -- astropy.units for dealing with astronomical units -- astropy.coordinates.SkyCoord for dealing with astronomical coordinates -- firefly_client.FireflyClient for using the Firefly python client -- astropy.nddata.Cutout2D for making image cutouts -- itertools.product to support looping over Roman blocks -- reproject.reproject_interp to convert Roman coadds from STG to TAN projection -- io.BytesIO for writing a fits file to an in-memory stream - -```{code-cell} ipython3 -# Uncomment the next line to install dependencies if needed. -# !pip install numpy astropy s3fs matplotlib firefly_client reproject -``` - -```{code-cell} ipython3 -from astropy.io import fits -import numpy as np -import s3fs -from matplotlib import pyplot as plt -from matplotlib import patches -from astropy import wcs -from astropy import units as u -from astropy.coordinates import SkyCoord -from firefly_client import FireflyClient -from astropy.nddata import Cutout2D -from itertools import product -from reproject import reproject_interp -from io import BytesIO -``` - -## 1. Learn where the OpenUniverse2024 data are hosted in the cloud - -The OpenUniverse2024 data preview is hosted in the cloud via Amazon Web Services (AWS). To access these data, you need to create a client to read data from Amazon's Simple Storage Service (s3) buckets, and you need to know some information about those buckets. The OpenUniverse2024 data preview contains simulations of the Roman Wide-Area Survey (WAS) and the Roman Time Domain Survey (TDS). In this tutorial, we will focus on the WAS. - -```{code-cell} ipython3 -BUCKET_NAME = "nasa-irsa-simulations" -ROMAN_PREFIX = "openuniverse2024/roman/preview" -ROMAN_COADD_PATH = f"{ROMAN_PREFIX}/RomanWAS/images/coadds" -TRUTH_FILES_PATH = f"{ROMAN_PREFIX}/roman_rubin_cats_v1.1.2_faint" - -RUBIN_PREFIX = "openuniverse2024/rubin/preview" -RUBIN_COADD_PATH = f"{RUBIN_PREFIX}/u/descdm/preview_data_step3_2877_19_w_2024_12/20240403T150003Z/deepCoadd_calexp/2877/19" -``` - -## 2. Roman Coadds - -The Nancy Grace Roman Space Telescope will carry out a wide-area survey (WAS) in the near infrared. The OpenUniverse2024 data preview includes coadded mosaics of simulated WAS data, created with the IMCOM algorithm (Rowe et al. 2011). Bands include F184, H158, J129, K213, Y106. In this section, we define some functions that make it convenient to retrieve a given cloud-hosted simulated Roman coadd based on position and filter. - -+++ - -### Define Roman Simulated "Blocks" - -The simulated Roman coadds are arranged in 100 arcsecond blocks, as described in Hirata et al. 2024. Below we define the RA and Dec of the block centers. This cell should not be altered. - -```{code-cell} ipython3 -#Centers of roman data preview blocks. Do not alter. -ra_block_centers = np.array([9.76330352298415, 9.724522605135252, 9.68574158906671, - 9.646960496603766, 9.608179349571955, 9.56939816979703, - 9.530616979104877, 9.491835799321422, 9.453054652272561, - 9.414273559784032, 9.375492543681393, 9.336711625789874]) * u.deg -dec_block_centers = np.array([-44.252584927082495, -44.22480733304182, -44.197029724175756, - -44.16925210374898, -44.14147447502621, -44.11369684127218, - -44.08591920575162, -44.05814157172923, -44.03036394246976, - -44.0025863212379, -43.974808711298394, -43.94703111591591]) * u.deg -block_size = 100 * u.arcsec # each block is 100 arcsec across -``` - -### Define a function that returns the simulated Roman block column/row for a given RA/Dec. -Blocks are arranged in a grid with axes parallel to RA and Dec. A given RA represents a column of blocks and a given Dec represents a row of blocks. This function returns the block column or row given the block centers from the cell above and a given RA or Dec, respectively. - -```{code-cell} ipython3 -def get_block_axis(block_centers, coord, ra_or_dec): - ra_or_dec_coord = getattr(coord, ra_or_dec) - block_dist_array = np.absolute(block_centers - ra_or_dec_coord) - closest_block_idx = block_dist_array.argmin() - if (ra_or_dec_coord < block_centers.min()-block_size/2 - or ra_or_dec_coord > block_centers.max()+block_size/2): - raise ValueError(f"Chosen {ra_or_dec}: {ra_or_dec_coord} not covered by OpenUniverse 2024 data preview simulated Roman coadds") - else: - return closest_block_idx + 12 # preview covers central 12 rows 12 columns, in a grid of 36x36 blocks -``` - -### Define a function that retrieves a Roman simulated coadd given a sky position and filter. - -+++ - -Each of the cloud-hosted simulated Roman coadds can be accessed via a S3 filepath. This function returns the access path for the simulated Roman coadd that includes a specified position on the sky and desired filter. - -```{code-cell} ipython3 -def get_roman_coadd_fpath(coord, filter): - col = get_block_axis(ra_block_centers, coord, 'ra') - row = get_block_axis(dec_block_centers, coord, 'dec') - - # Construct the coadd filename from the chosen filter, row, and column. - coadd_fname_root = f"prod_{filter[0]}_{col}_{row}_map.fits" - coadd_fpath = f"{BUCKET_NAME}/{ROMAN_COADD_PATH}/{filter}/Row{row}/{coadd_fname_root}" - return coadd_fpath -``` - -Now we use this access path and prefix it with `s3://` and use astropy.fits to extract a subset of it (for more info see [this section](https://docs.astropy.org/en/stable/io/fits/usage/cloud.html#subsetting-fits-files-hosted-in-amazon-s3-cloud-storage) of astropy docs). - -We use `.section` to retrieve just the science image data from fits HDU as a 2D `numpy.array`, and extract WCS information from fits header as `astropy.wcs.WCS` object. The following function returns a dictionary of both. - -```{code-cell} ipython3 -def get_roman_coadd(coord, filter): - # retrive fits file of block/tile from the coadd mosiac - coadd_s3_fpath = get_roman_coadd_fpath(coord, filter) - coadd_s3_uri = f"s3://{coadd_s3_fpath}" - - with fits.open(coadd_s3_uri, fsspec_kwargs={"anon": True}) as hdul: - # retrieve science data from coadd fits - coadd_data = hdul[0].section[0,0, :, :] # has (2688, 2688, 15, 1) shape, with 0th layer in the cube as science image - - # make wcs using header - coadd_wcs = wcs.WCS(hdul[0].header, naxis=2) - - return {'data': coadd_data, 'wcs': coadd_wcs} -``` - -### Inspect a simulated Roman Coadd - -+++ - -Choose a filter and position that lies within the data preview region - -```{code-cell} ipython3 -coord = SkyCoord(ra=9.6055383, dec=-44.1895542, unit="deg") -filter_roman = 'H158' #F184, H158, J129, K213, and Y106 are available in the data preview -``` - -Retrieve the data and header information from the simulated Roman coadd corresponding to the chosen position and filter. - -```{code-cell} ipython3 -coadd_roman = get_roman_coadd(coord, filter_roman) -``` - -### Understand the size of a simulated Roman coadd. - -```{code-cell} ipython3 -# Number of pixels (Y, X) -coadd_roman['data'].shape -``` - -```{code-cell} ipython3 -# Pixel size (scale Y, scale X) [degrees/pixel] -coadd_roman['wcs'].proj_plane_pixel_scales() -``` - -```{code-cell} ipython3 -# Coadd size (FOV Y, FOV X) -[(num * size).to('arcsec') for num, size in zip( - coadd_roman['data'].shape, coadd_roman['wcs'].proj_plane_pixel_scales())] -``` - -The field of view of Roman coadd is ~100 arcsec. - -+++ - -### Use the WCS from the Roman simulated header to convert the specified coordinate into a pixel position. - -```{code-cell} ipython3 -def coord_to_xy(w, coord): - return w.world_to_array_index(coord)[::-1] #reverse since 0th axis is y, 1st axis is x - -coord_arr_idx = coord_to_xy(coadd_roman['wcs'], coord) -coord_arr_idx -``` - -### Use matplotlib imshow to create a static visualization of the Roman simulated coadd and overplot the selected position. - -```{code-cell} ipython3 -def stretch_color(data, clipPercent): - return np.percentile(data, (0 + clipPercent, 100 - clipPercent)) - -plt.imshow(coadd_roman['data'], origin='lower', - clim=stretch_color(coadd_roman['data'], 1) - ) - -plt.plot(*coord_arr_idx, 'r+', markersize=15) -``` - -## 3. Rubin Coadds - -The OpenUniverse2024 data preview includes coadded mosaics in the following filters: u, g, r, i, z, y. In this section, we define some functions that make it convenient to retrieve a given cloud-hosted simulated Roman coadd based on position and filter. - -+++ - -### Retrieve Rubin Coadds - -The OpenUniverse2024 data preview includes only one simulated Rubin coadd per filter. Below we define functions that make it convenient to retrieve the simulated Rubin coadd corresponding to the desired filter. We return data in same structure as the functions we defined above for Roman. - -```{code-cell} ipython3 -def get_rubin_coadd_fpath(filter): - coadd_fname_root = f"deepCoadd_calexp_2877_19_{filter}_DC2_u_descdm_preview_data_step3_2877_19_w_2024_12_20240403T150003Z.fits" - coadd_fpath = f"{BUCKET_NAME}/{RUBIN_COADD_PATH}/{filter}/{coadd_fname_root}" - return coadd_fpath -``` - -```{code-cell} ipython3 -def get_rubin_coadd(filter): - coadd_s3_fpath = get_rubin_coadd_fpath(filter) - - with fits.open(f"s3://{coadd_s3_fpath}", fsspec_kwargs={"anon": True}) as hdul: - # retrieve science data from coadd fits - coadd_data = hdul[1].section[:,:] - - # make wcs using header - coadd_wcs = wcs.WCS(hdul[1].header) - - return {'data': coadd_data, 'wcs': coadd_wcs} -``` - -### Inspect a simulated Rubin Coadd - -+++ - -Choose a filter and retrieve the data and header information from the simulated Rubin coadd corresponding to that filter. - -```{code-cell} ipython3 -filter_rubin = 'r' -coadd_rubin = get_rubin_coadd(filter_rubin) -``` - -### Understand the size of a simulated Rubin coadd. - -```{code-cell} ipython3 -# Number of pixels (Y, X) -coadd_rubin['data'].shape -``` - -```{code-cell} ipython3 -# Pixel size (scale Y, scale X) [degrees/pixel] -coadd_rubin['wcs'].proj_plane_pixel_scales() -``` - -```{code-cell} ipython3 -# Coadd size (FOV Y, FOV X) -[(num * size).to('arcsec') for num, size in zip( - coadd_rubin['data'].shape, coadd_rubin['wcs'].proj_plane_pixel_scales())] -``` - -The field of view of Rubin coadd is 840 arcsec. - -+++ - -### Use matplotlib imshow to create a static visualization of the Rubin simulated coadd and overplot the selected position. - -```{code-cell} ipython3 -plt.imshow(coadd_rubin['data'], origin='lower', - clim=stretch_color(coadd_rubin['data'], 1) - ) - -plt.plot(*coord_to_xy(coadd_rubin['wcs'], coord), 'r+', markersize=15) -``` - -### Define a function that returns the URL for a given S3 filepath -Since the OpenUniverse2024 data is available through a public S3 bucket, we can access a given S3 file using HTTPS URL as follows: - -```{code-cell} ipython3 -def https_url(s3_fpath): - s3_fpath_without_bucket = s3_fpath.split('/', 1)[1] - return f"https://{BUCKET_NAME}.s3.amazonaws.com/{s3_fpath_without_bucket}" -``` - -Let's generate URL for the Rubin coadd we plotted above. Clicking on the returned URL will allow you to download this image locally. - -```{code-cell} ipython3 -coadd_s3_fpath_rubin = get_rubin_coadd_fpath(filter_rubin) -https_url(coadd_s3_fpath_rubin) -``` - -## 4. Compare simulated Roman and Rubin cutouts for a selected position - -+++ - -### Choose cutout size - -```{code-cell} ipython3 -cutout_size = 50*u.arcsec -``` - -### Create the cutouts - -```{code-cell} ipython3 -cutout_roman = Cutout2D(coadd_roman['data'], coord, size=cutout_size, wcs=coadd_roman['wcs']) -cutout_rubin = Cutout2D(coadd_rubin['data'], coord, size=cutout_size, wcs=coadd_rubin['wcs']) -``` - -### Use matplotlib imshow to plot static side-by-side comparisons of the cutouts - -```{code-cell} ipython3 -fig, axs = plt.subplots(1, 2, figsize=(12, 6)) - - -axs[0].imshow(cutout_roman.data, origin='lower', - clim=stretch_color(cutout_roman.data, .5) - ) -axs[0].plot(*coord_to_xy(cutout_roman.wcs, coord), 'r+', markersize=15) -axs[0].set_title(f"ROMAN in filter {filter_roman}") - -axs[1].imshow(cutout_rubin.data, origin='lower', - clim=stretch_color(cutout_rubin.data, .5) - ) -axs[1].plot(*coord_to_xy(cutout_rubin.wcs, coord), 'r+', markersize=15) -axs[1].set_title(f"RUBIN in filter {filter_rubin}") - -fig.suptitle(f"Cutouts at ({coord.ra}, {coord.dec}) with {cutout_size} size", fontsize=14) -plt.tight_layout(rect=[0, 0, 1, 0.97]) -``` - -## 5. Use Firefly to interactively identify a blended source - -Clearly, the simulated Roman coadd has higher spatial resolution than the Rubin simulated coadd. Let's try to locate blended objects to compare in the simulated Rubin and Roman images. We will use Firefly's interactive visualization to make this task easier. - -+++ - -### Launch and initialize Firefly -There are two ways to initialize a Firefly client from Python, depending on whether you're running the notebook in JupyterLab or not. Assuming you have `jupyter-firefly-extensions` set up in your environment as explained [here](https://github.com/Caltech-IPAC/jupyter_firefly_extensions/blob/master/README.md), you can use `make_lab_client()` in JupyterLab, which will open the Firefly viewer in a new tab within the Lab. Otherwise, you can use `make_client()` in a Jupyter Notebook (or even a Python shell), which will open the Firefly viewer in a new web browser tab. - -You also need a Firefly server to communicate with your Firefly Python client. In this notebook, we use a public Firefly server: the IRSA Viewer (https://irsa.ipac.caltech.edu/irsaviewer). However, you can also run a local Firefly server via a [Firefly Docker image](https://hub.docker.com/r/ipac/firefly) and access it at `http://localhost:8080/firefly`. The URL of the Firefly server is read by both `make_client()` and `make_lab_client()` through the environment variable `FIREFLY_URL`. However, `make_client()` also allows you to pass the URL directly as the `url` parameter. - -```{code-cell} ipython3 -# Uncomment when using within Jupyter Lab with jupyter_firefly_extensions installed -# fc = FireflyClient.make_lab_client() - -# Uncomment for contexts other than above -fc = FireflyClient.make_client(url="https://irsa.ipac.caltech.edu/irsaviewer") - -fc.reinit_viewer() # to clean the state, if this cell ran earlier -``` - -### Send the simulated Rubin coadd to Firefly using show_fits. - -For displaying the FITS image of Rubin coadd in Firefly, we use [`show_fits`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.show_fits): - -```{code-cell} ipython3 -coadd_ff_id_rubin = 'rubin-coadd-filter-r' -fc.show_fits(url=https_url(coadd_s3_fpath_rubin), - plot_id=coadd_ff_id_rubin, - Title="Rubin Coadd" - ) -``` - -### Use ds9 region syntax to overplot the simulated Roman image blocks on the interactive display - -The Firefly client includes several methods related to controlling ds9 region overlays. To -overlay a region layer on the loaded FITS images, we can use [`overlay_region_layer`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.overlay_region_layer). - -Region data is defined in ds9 region syntax that can be found [here](https://ds9.si.edu/doc/ref/region.html). - -```{code-cell} ipython3 -# mark the roman coadd blocks as boxes -roman_regions = [ - f'icrs;box {ra_block_center.value}d {dec_block_center.value}d {block_size.value}" {block_size.value}" 0d' - for (ra_block_center, dec_block_center) in product(ra_block_centers, dec_block_centers) -] - -roman_regions_id = 'roman_regions' -fc.overlay_region_layer(region_data=roman_regions, - title='Roman Mosiac', - region_layer_id=roman_regions_id) -``` - -### Use Firefly's pan and zoom capabilities to locate a region of interest (a blended source) - -You can view the coordinates of your mouse pointer at the bottom left of the display window. To copy the coordinates for a specific coordinate: - -- Toggle the "Click Lock" to "on" in the bottom right of the image display. -- Click on the position of interest, and notice that the coordinate display is now frozen. -- Click on "EQ-J2000", the coordinate label in the bottom left of the image display. In the dialog that opens, change copy options to "[Python] Astropy SkyCoord" so that we can directly work with them in python. -- Close the dialog and click on the copy icon next to the coordinate values display. - -+++ - -### Copy the coordinates from the coordinate display to the Python notebook - -We have provided an example. You can change this based on your interests. - -```{code-cell} ipython3 -coords_of_interest = SkyCoord('0h38m25.35s -44d00m10.1s', frame='icrs') # located and copied through UI -coords_of_interest -``` - -We can now use this [astropy `SkyCoord` object](https://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html#astropy.coordinates.SkyCoord) to compare our coadds. - -+++ - -### Use ds9 region syntax to overplot the selected position - -For this we use the id of the region layer we defined above, and add more region data using [`add_region_data`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.add_region_data). - -```{code-cell} ipython3 -point_region = f'icrs;point {coords_of_interest.ra.value}d {coords_of_interest.dec.value}d # point=cross 15 text={{Blended source}}' -fc.add_region_data(region_data=point_region, region_layer_id=roman_regions_id) -``` - -## 6. Plot cutouts of the identified blended source - -```{code-cell} ipython3 -coadd_roman = get_roman_coadd(coords_of_interest, filter_roman) -``` - -```{code-cell} ipython3 -cutout_size = 20*u.arcsec -``` - -```{code-cell} ipython3 -cutout_roman = Cutout2D(coadd_roman['data'], coords_of_interest, size=cutout_size, wcs=coadd_roman['wcs']) -cutout_rubin = Cutout2D(coadd_rubin['data'], coords_of_interest, size=cutout_size, wcs=coadd_rubin['wcs']) -``` - -```{code-cell} ipython3 -fig, axs = plt.subplots(1, 2, figsize=(12, 6)) - -axs[0].imshow(cutout_roman.data, origin='lower', - clim=stretch_color(cutout_roman.data, .5) - ) -axs[0].set_title(f"ROMAN in filter {filter_roman}") - -# Let's also encircle the blended source we identified -axs[0].add_patch(patches.Circle(coord_to_xy(cutout_roman.wcs, coords_of_interest), - radius=50, color='r', fill=False, linewidth=2)) -# and a bonus blended source that is close to it -other_coords = SkyCoord(coords_of_interest.ra-9.2*u.arcsec, coords_of_interest.dec+5*u.arcsec) -axs[0].add_patch(patches.Circle(coord_to_xy(cutout_roman.wcs, other_coords), - radius=36, color='cyan', fill=False, linewidth=2)) - - -axs[1].imshow(cutout_rubin.data, origin='lower', - clim=stretch_color(cutout_rubin.data, .5) - ) -axs[1].set_title(f"RUBIN in filter {filter_rubin}") - -# Let's also encircle the source we identified -axs[1].add_patch(patches.Circle(coord_to_xy(cutout_rubin.wcs, coords_of_interest), - radius=10, color='r', fill=False, linewidth=2)) -axs[1].add_patch(patches.Circle(coord_to_xy(cutout_rubin.wcs, other_coords), - radius=8, color='cyan', fill=False, linewidth=2)) - - -fig.suptitle(f"Cutouts at ({coords_of_interest.ra:6f}, {coords_of_interest.dec:6f}) with {cutout_size} size", fontsize=14); -plt.tight_layout(rect=[0, 0, 1, 0.97]) -# plt.savefig("plot.pdf", bbox_inches='tight', pad_inches=0.2) -``` - -## 7. Use Firefly to visualize the OpenUniverse2024 data preview catalogs -Let's inspect the properties of sources in the Rubin coadd image. For this we will use the input truth files present in S3 bucket. - -The OpenUniverse2024 data preview includes the input truth files that were used to create the simulated images. These files are in Parquet and HDF5 format, and include information about the properties of galaxies, stars, and transients. - -To list the available files, we use [s3fs](https://s3fs.readthedocs.io/en/latest/index.html). We open a file-system like connection -to AWS S3 and then simply list (`ls`) truth files directory: - -```{code-cell} ipython3 -s3 = s3fs.S3FileSystem(anon=True) # to browse s3 bucket -s3.ls(f"{BUCKET_NAME}/{TRUTH_FILES_PATH}") -``` - -```{code-cell} ipython3 -# Catalog table of star properties (in parquet format) -pointsource_cat_path = f"{BUCKET_NAME}/{TRUTH_FILES_PATH}/pointsource_10307.parquet" -pointsource_cat_path -``` - -```{code-cell} ipython3 -# Catalog table of galaxy properties (in parquet format) -galaxy_cat_path = f"{BUCKET_NAME}/{TRUTH_FILES_PATH}/galaxy_10307.parquet" -galaxy_cat_path -``` - -### Use Firefly's show_table to overlay the catalogs on interactive image of coadd - -The input truth files cover a region much larger than the data preview, so we define filters on table to constraint the sources in catalog within the ra and dec bounds of preview data. (Note: you can remove filters through the table UI if you wish to see the entire data) - -+++ - -You can visualize catalogs interactively with Firefly using [`show_table`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.show_table). This capability can take many parameters. Here we will simply send our catalog to Firefly so that we can (a) see an interactive table; (b) see this table plotted over the image that we've already sent; and (c) use the GUI to quickly create exploratory plots. See if you can use the GUI to quickly determine approximately how many galaxies cover the Rubin image and what the redshift distribution of these galaxies is. - -```{code-cell} ipython3 -cat_filters = [ - f'("ra" >= {ra_block_centers.min().value} AND "ra" <= {ra_block_centers.max().value})', - f'("dec" >= {dec_block_centers.min().value} AND "dec" <= {dec_block_centers.max().value})' -] -cat_filters -``` - -```{code-cell} ipython3 -fc.show_table(url=https_url(pointsource_cat_path), - title='Stars Catalog', - tbl_id='stars_cat', - filters=" AND ".join(cat_filters)) -``` - -```{code-cell} ipython3 -gal_cat_tbl_id = 'galaxy_cat' - -# may take ~1.25min, because galaxy catalog is a big file -fc.show_table(url=https_url(galaxy_cat_path), - title='Galaxy Catalog', - tbl_id=gal_cat_tbl_id, - filters=" AND ".join(cat_filters)) -``` - -For each row in the table you can notice a marker in the image. Selecting a row or marker changes the corresponding marker or row, respectively. You can click on "Details" tab in the UI to show properties of each source selected in image/table. - -### Use Firefly's apply_table_filters to show only high-redshift galaxies - -High redshift galaxies are the most interesting, so let's filter the table we sent to Firefly to only include z>3 galaxies. Notice how the table display and image overlay change. Notice how the chart becomes a scatterplot from a heatmap because the sources reduce. You can remove this filter or add new ones through the GUI. - -For filtering, we will use [`apply_table_filters`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.apply_table_filters) method on the galaxy table we loaded above. - -```{code-cell} ipython3 -fc.apply_table_filters(tbl_id=gal_cat_tbl_id, - filters=" AND ".join(cat_filters+['"redshift" > 3'])) -``` - -You can play with the filters directly from the UI as well. Try removing adding more filters in the tables and see how markers change. - -+++ - -### Use Firefly's show_fits_3color to create a 3 color image of the simulated Rubin images - -```{code-cell} ipython3 -# [R, G, B] -ROMAN_RGB_FILTERS = ['H158', 'J129', 'Y106'] -RUBIN_RGB_FILTERS = ['r', 'g', 'u'] -``` - -We already have Rubin coadd with catalog overlaid, let's make a 3 color image to see colors of marked objects more clearly. For this we will use [`show_fits_3color`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.show_fits_3color) method: - -```{code-cell} ipython3 -coadd_ff_id_rubin_3color = 'rubin-coadd-3color' -threeC = [ - dict(url=https_url(get_rubin_coadd_fpath(filter_name)), - Title="Rubin Coadd 3 color") - for filter_name in RUBIN_RGB_FILTERS -] - -fc.show_fits_3color(three_color_params=threeC, - plot_id=coadd_ff_id_rubin_3color) -``` - -### Use Firefly's interactivity to identify a region of interest - -For example, we found a region of the sky that seems to have a high number of high redshift sources and we copy it from the image display: - -```{code-cell} ipython3 -# located and copied through UI -high_z_gal_coords = SkyCoord('0h38m00.77s -44d12m10.2s', frame='icrs') -high_z_gal_coords -``` - -```{code-cell} ipython3 -# let's also mark it in our region layer, so that it's easy to pinpoint later -point_region = f'icrs;point {high_z_gal_coords.ra.value}d {high_z_gal_coords.dec.value}d # point=cross 15 text={{z>3 mock galaxies}}' -fc.add_region_data(region_data=point_region, region_layer_id=roman_regions_id) -``` - -## 8. Plot 3-color Roman coadd containing your region of interest -Let's inspect WCS of Roman coadd first - -```{code-cell} ipython3 -coadd_roman['wcs'] -``` - -### Prepare Roman coadds for displaying in Firefly - -Roman coadds have STG projection which cannot be read by Firefly yet. Firefly can display any FITS image but it needs to read the WCS for overlaying catalogs and other interactive features. So unlike Rubin 3 color image where we directly passed URL of coadd files to Firefly, we will read Roman coadd files in Python, reproject them from STG to TAN, and write them back to FITS to pass them to Firefly. - -Let's first define functions to do so: - -```{code-cell} ipython3 -def reproject_to_TAN(coadd_roman): - # Define a new WCS with TAN projection (in CTYPE key) - output_wcs = coadd_roman['wcs'].deepcopy() - output_wcs.wcs.ctype = [ctype.replace('STG', 'TAN') for ctype in coadd_roman['wcs'].wcs.ctype] - - # Use reproject to convert a given data and wcs, to a desired wcs and shape - reprojected_data, _ = reproject_interp( - (coadd_roman['data'], coadd_roman['wcs']), - output_projection=output_wcs, - shape_out=coadd_roman['data'].shape - ) - - return {'data': reprojected_data, 'wcs': output_wcs} -``` - -```{code-cell} ipython3 -def get_fits_stream(coadd_roman): - # Create a FITS PrimaryHDU object with the coadd data - hdu = fits.PrimaryHDU(data=coadd_roman['data'], header=coadd_roman['wcs'].to_header()) - - # Write the HDU to the in-memory stream (to save I/O time) - fits_stream = BytesIO() - hdu.writeto(fits_stream, overwrite=True) - fits_stream.seek(0) # to bring reading pointer to the beginning of file - - return fits_stream -``` - -Then we perform all 3 operations we mentioned above for the RGB filters of Roman: - -```{code-cell} ipython3 -coadds_rgb = [] -coadds_rgb_reprojected = [] -coadds_rgb_fits_stream = [] - -for filter_name in ROMAN_RGB_FILTERS: - print(f'\nFILTER: {filter_name}') - print('Retrieving Roman coadd...') - coadd_roman = get_roman_coadd(high_z_gal_coords, filter_name) - coadds_rgb.append(coadd_roman) - - print('Reprojecting to TAN...') - coadd_roman_reprojected = reproject_to_TAN(coadd_roman) - coadds_rgb_reprojected.append(coadd_roman_reprojected) - - print('Writing back to fits stream...') - coadd_roman_fits_stream = get_fits_stream(coadd_roman_reprojected) - coadds_rgb_fits_stream.append(coadd_roman_fits_stream) -``` - -### Use Firefly's show_fits_3color to create a 3 color image of Roman coadds - -Now we upload each fits stream (in-memory fits file) to firefly using [`upload_fits_data()`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.upload_fits_data) and prepare color params to pass to the [`show_fits_3color()`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.show_fits_3color). - -```{code-cell} ipython3 -three_color_params = [ - { - 'file': fc.upload_fits_data(fits_stream), - 'Title': "Roman Coadd 3 color" - } for fits_stream in coadds_rgb_fits_stream] -``` - -```{code-cell} ipython3 -coadd_ff_id_roman_3color = 'roman-coadd-3color-high_z_gal' -fc.show_fits_3color(three_color_params=three_color_params, - plot_id=coadd_ff_id_roman_3color) -``` - -We can see 3 color image of Roman coadd containing the high-redshift galaxy sources. Try panning and zomming out, you can notice it spans over one block compared to the Rubin coadd which is much larger. - -### Use Firefly's pan/zoom/align methods to locate high redshift sources -Now, let's pan & zoom to the region where we located high-redshift galaxy sources. Also align & lock all images being displayed by WCS. For these operations we use these 3 methods (click on them to see their documentation): -- [`set_pan`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.set_pan) -- [`set_zoom`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.set_zoom) -- [`align_images`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.align_images) - -```{code-cell} ipython3 -fc.set_pan(plot_id=coadd_ff_id_roman_3color, x=high_z_gal_coords.ra.deg, y=high_z_gal_coords.dec.deg, coord='j2000') -fc.set_zoom(plot_id=coadd_ff_id_roman_3color, factor=1) -fc.align_images(lock_match=True) -``` - -### Use Firefly's set_stretch method to change the stretch of the image display via Python - -The image has a lot of noise that obscures our high redshift sources of interest. You can use the Firefly GUI to change the stretch of the image display. We identify that squared stretch from -2 to 10 sigma highlights the colors of our sources better. You can also use the Firefly client's [`set_stretch`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.set_stretch) to do this via Python. This is helpful for reproducibility and for scaling up to many images. - -```{code-cell} ipython3 -fc.set_stretch(plot_id=coadd_ff_id_roman_3color, stype='sigma', algorithm='squared', - band='ALL', lower_value=-2, upper_value=10) -``` - -*** - -## About This Notebook - -**Updated:** 2024-12-19 - -**Contact:** [the IRSA Helpdesk](https://irsa.ipac.caltech.edu/docs/help_desk.html) with questions or reporting problems. diff --git a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md index 9bf9704e..07a369fb 100644 --- a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md +++ b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md @@ -209,9 +209,16 @@ def assemble_SN_data(sn_flux_file, galaxy_flux_file, galaxy_info_file, *, region sn_df = pq.read_table(sn_flux_file, filesystem=fs).to_pandas() print(f" Loaded {len(sn_df)} SN entries") - # Load galaxy flux and info tables + # Load galaxy flux and info tables. The flux table is all fluxes and we want + # every band, but the info table carries a lot of shape and shear columns we + # never touch, so we ask Parquet for just the three we do use. That leaves + # about a sixth of the data to pull across and hold in memory. gal_flux = pq.read_table(galaxy_flux_file, filesystem=fs).to_pandas() - gal_info = pq.read_table(galaxy_info_file, filesystem=fs).to_pandas() + gal_info = pq.read_table( + galaxy_info_file, + filesystem=fs, + columns=["galaxy_id", "redshift", "um_source_galaxy_obs_sm"], + ).to_pandas() # Join host-galaxy flux and info tables on galaxy_id gal_joined = gal_flux.merge(gal_info, on="galaxy_id", how="inner") @@ -706,7 +713,7 @@ def plot_many_seds(df, rubin_bands, roman_bands, n_galaxies=10, loglog=False): plot_many_seds(df_scaled, rubin_bands, roman_bands, n_galaxies=50, loglog=True) ``` -+++ {"jupyter": {"source_hidden": true}} ++++ ## 3. SED fitting In this section we will use [Prospector](https://prospect.readthedocs.io/en/stable/index.html), a Bayesian SED fitting code built on [FSPS](https://dfm.io/python-fsps/current/), to infer stellar population parameters for our galaxies by fitting their SEDs. We choose Prospector because it is an open-source, powerful, yet flexible package to infer stellar population properties that is widely used in the community. diff --git a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md index 8a37e265..1b0e4919 100644 --- a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md +++ b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md @@ -82,11 +82,11 @@ import itertools ## 1. Explore the OpenUniverse2024 data directories -This section of the tutorial demonstrates how to explore the OpenUniverse2024 data directories directly on S3 and inspect simulated Roman and Rubin images without downloading large datasets locally. It establishes a connection to the public NASA IRSA simulations bucket using s3fs, defines key directory paths for the full Roman and Rubin simulations (not the preview subsets), and illustrates how to browse image files for a selected band and pointing. The accompanying functions — summarize_fits_files() and show_gallery() — provide tools for quickly summarizing FITS file metadata (e.g., number of extensions, pointing information, pixel scale) and for visualizing a small gallery of example images from the chosen directory. +This section of the tutorial demonstrates how to explore the OpenUniverse2024 data directories directly on S3 and inspect simulated Roman and Rubin images without downloading large datasets locally. It establishes a connection to the public NASA IRSA simulations bucket using s3fs, defines key directory paths for the full Roman and Rubin simulations, and illustrates how to browse image files for a selected band and pointing. The accompanying functions — summarize_fits_files() and show_gallery() — provide tools for quickly summarizing FITS file metadata (e.g., number of extensions, pointing information, pixel scale) and for visualizing a small gallery of example images from the chosen directory. In the prefix you will see that we choose "simple_model" simulations and not "truth" simulations because the simple_model images are the ones with noise and real effects, while "Truth" are noise free, perfect images. -Also in the prefix you will see that we choose the full simulation, not the preview simulation for both Roman and Rubin. Differences between the "full" and "preview" simulations are clarified in the [this](https://arxiv.org/abs/2501.05632) publication +The prefix also selects the full simulation, which covers the entire survey footprint, for both Roman and Rubin. The [OpenUniverse2024 paper](https://arxiv.org/abs/2501.05632) describes how the simulation was produced. ```{code-cell} ipython3 # Setup diff --git a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_firefly_visualization.md b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_firefly_visualization.md new file mode 100644 index 00000000..71400cf2 --- /dev/null +++ b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_firefly_visualization.md @@ -0,0 +1,987 @@ +--- +authors: +- name: Jaladh Singhal +- name: Vandana Desai +- name: IRSA Team +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.16.3 +kernelspec: + name: python3 + display_name: python3 + language: python +--- + +# Using Firefly to Explore OpenUniverse2024 Simulated Roman and Rubin Images + ++++ + +## Learning Goals + +By the end of this tutorial, you will: + +- Learn how to access cloud-hosted Roman and Rubin simulated images. + +- Learn how to find the Roman and Rubin images covering a sky position with the IRSA Simple Image Access (SIA) service. + +- Learn how to launch an interactive Firefly instance inside JupyterLab. + +- Learn how to use the Firefly Jupyterlab extension to visualize cloud-hosted simulated images, overplot ds9 regions, overplot catalogs in Parquet format, and create 3 color images. + ++++ + +## Introduction + +The purpose of this tutorial is to become familiar with the simulated Roman and Rubin data published through OpenUniverse2024, and to become familiar with the Firefly JupyterLab Extension for visualizing astronomical data products. + +OpenUniverse2024 is a project to simulate spatially overlapping imaging surveys to be carried out by the Nancy Grace Roman Telescope and the Vera C. Rubin Observatory. The simulations were carried out on Argonne's Theta cluster and consist of: + +- The LSST ELAIS-S1 Deep Drilling Field (DDF) +- The Roman Time-Domain Survey (TDS) shifted to overlap the ELAIS region and LSST DDF +- Overlapping LSST Wide-Fast-Deep (WFD) survey (with rolling cadence) +- Overlapping Roman Wide-Area Survey (WAS) in the same region +- A deep-field calibration region of the Roman WAS in the same region + +This tutorial works with the full simulation, which covers the entire survey footprint. More information about the dataset can be found at [IRSA's holding of this dataset](https://irsa.ipac.caltech.edu/data/theory/openuniverse2024/overview.html), and the [OpenUniverse2024 paper](https://arxiv.org/abs/2501.05632) describes how the simulation was produced. + +Firefly is an open-source web-based UI library for astronomical data archive access and visualization developed at Caltech and used by multiple space- and ground-based astrophysics archives. More information on Firefly can be found [here](https://github.com/Caltech-IPAC/firefly/blob/dev/README.md). + +In addition to being used to make web applications, Firefly can be used from Python. More information on Firefly Python client can be found [here](https://caltech-ipac.github.io/firefly_client/usage/index.html). + +The Firefly JupyterLab Extension makes it particularly easy to use Firefly to efficiently visualize cloud-hosted astronomical data using JupyterLab instances running locally or on cloud. More information on Firefly JupyterLab Extension can be found [here](https://github.com/Caltech-IPAC/jupyter_firefly_extensions/blob/master/README.md). + +If you are new to OpenUniverse2024, the [Quickstart](openuniverse2024_quickstart) tutorial introduces the directory layout, the parquet catalogs, and the SIA image search that this notebook builds on. + ++++ + +## Imports + +- astropy.io.fits for accessing FITS files +- numpy for numerical computing +- matplotlib.pyplot for creating static visualizations of FITS images +- matplotlib.patches for annotating visualizations of FITS images +- astropy.wcs for dealing with astronomical world coordinate systems +- astropy.units for dealing with astronomical units +- astropy.coordinates.SkyCoord for dealing with astronomical coordinates +- astroquery.ipac.irsa.Irsa for asking IRSA which images cover a position +- hpgeom for converting a sky position into a HEALPix region index +- json for reading the cloud location returned by the image search +- firefly_client.FireflyClient for using the Firefly python client +- astropy.nddata.Cutout2D for making image cutouts +- reproject.reproject_interp for rotating cutouts so both telescopes' panels point the same way +- io.BytesIO for writing a fits file to an in-memory stream + +```{code-cell} ipython3 +# Uncomment the next line to install dependencies if needed. +# !pip install numpy astropy s3fs matplotlib firefly_client astroquery hpgeom reproject +``` + +```{code-cell} ipython3 +from astropy.io import fits +import numpy as np +from matplotlib import pyplot as plt +from matplotlib import patches +from astropy import wcs +from astropy import units as u +from astropy.coordinates import SkyCoord +from astroquery.ipac.irsa import Irsa +import hpgeom +import json +from firefly_client import FireflyClient +from astropy.nddata import Cutout2D +from reproject import reproject_interp +from io import BytesIO +``` + +## 1. Learn where the OpenUniverse2024 data are hosted in the cloud + +The OpenUniverse2024 data are hosted in the cloud via Amazon Web Services (AWS). To access these data, you need to create a client to read data from Amazon's Simple Storage Service (s3) buckets, and you need to know some information about those buckets. OpenUniverse2024 contains simulations of the Roman Wide-Area Survey (WAS) and the Roman Time Domain Survey (TDS). In this tutorial, we will focus on the WAS. + +We ask IRSA's Simple Image Access (SIA) service to locate both the Roman and Rubin images by coordinate. + +```{code-cell} ipython3 +BUCKET_NAME = "nasa-irsa-simulations" +OU_PREFIX = "openuniverse2024" +TRUTH_FILES_PATH = f"{OU_PREFIX}/roman/full/roman_rubin_cats_v1.1.2_faint" +``` + +```{code-cell} ipython3 +# Point the astroquery IRSA client at the simulated-data services, which are +# separate from the ones serving IRSA's observed data. +Irsa.sia_url = "https://irsa.ipac.caltech.edu/simulated/SIA" +Irsa.tap_url = "https://irsa.ipac.caltech.edu/simulated/TAP" + +OU_ROMAN_SIA_COLLECTION = "simulated_roman_openuniverse2024" +OU_RUBIN_SIA_COLLECTION = "simulated_rubin_openuniverse2024" + +# A small radius is all we need, since we only want images containing a given point. +SEARCH_RADIUS = 1 * u.arcsec +``` + +## 2. Roman Images + +The Nancy Grace Roman Space Telescope will carry out a wide-area survey (WAS) in the near infrared. OpenUniverse2024 simulates that survey as individual exposures, each covering one detector's worth of sky. Bands include F184, H158, J129, K213, W146, Y106. + +The survey stores those exposures by pointing and detector rather than by sky position, so in this section we define a function that asks the image search which exposure covers a position, and then reads it. + ++++ + +### Define a function that retrieves a Roman simulated image given a sky position and filter. + +The Roman collection holds both the Wide-Area Survey and the Time Domain Survey, so we keep only the WAS exposures and then only the band we want. Many exposures cover any given position, so we sort by observation time and take the earliest. + +Having found the file, we prefix its path with `s3://` and open it with astropy. The science image sits in the first extension, and we take the WCS from the same header as an `astropy.wcs.WCS` object. The function returns a dictionary of both, along with the header itself so we can look at it later, and the path so that later cells can reuse it instead of searching a second time. + +```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- +def get_roman_image(coord, filter): + """Read the Roman WAS exposure covering a position in a given filter. + + Parameters + ---------- + coord : `~astropy.coordinates.SkyCoord` + Sky position to search around. + filter : str + Roman bandpass name, e.g. ``'H158'``. + + Returns + ------- + dict + Keys ``'data'`` (`~numpy.ndarray`), ``'wcs'`` (`~astropy.wcs.WCS`), + ``'header'`` (`~astropy.io.fits.Header`) and ``'fpath'`` (str, the S3 + path of the exposure). + + Raises + ------ + ValueError + If no WAS exposure in that filter covers ``coord``. + """ + results = Irsa.query_sia(pos=(coord, SEARCH_RADIUS), + collection=OU_ROMAN_SIA_COLLECTION) + + in_was = results[['WAS_simple_model' in str(r['obs_id']) for r in results]] + in_band = in_was[[str(r['energy_bandpassname']).strip() == filter + for r in in_was]] + if len(in_band) == 0: + raise ValueError(f"No Roman WAS {filter}-band images cover " + f"{coord.to_string('hmsdms')}") + + in_band.sort('t_min') + cloud_info = json.loads(in_band['cloud_access'][0])['aws'] + image_s3_fpath = f"{cloud_info['bucket_name']}/{cloud_info['key']}" + + with fits.open(f"s3://{image_s3_fpath}", fsspec_kwargs={"anon": True}) as hdul: + # retrieve science data, which sits in the first extension + image_data = hdul[1].data + + # make wcs using header + image_wcs = wcs.WCS(hdul[1].header) + + return {'data': image_data, 'wcs': image_wcs, 'header': hdul[1].header, + 'fpath': image_s3_fpath} +``` + +### Inspect a simulated Roman image + ++++ + +Choose a filter and a position within the survey footprint + +```{code-cell} ipython3 +coord = SkyCoord(ra=9.6205000, dec=-44.0641694, unit="deg") +filter_roman = 'H158' #F184, H158, J129, K213, W146, and Y106 are available +``` + +Retrieve the data and header information from the simulated Roman image corresponding to the chosen position and filter. + +```{code-cell} ipython3 +image_roman = get_roman_image(coord, filter_roman) +``` + +### Understand the size of a simulated Roman image. + +```{code-cell} ipython3 +# Number of pixels (Y, X) +image_roman['data'].shape +``` + +```{code-cell} ipython3 +# Pixel size (scale Y, scale X) [degrees/pixel] +image_roman['wcs'].proj_plane_pixel_scales() +``` + +```{code-cell} ipython3 +# Image size (FOV Y, FOV X) +[(num * size).to('arcsec') for num, size in zip( + image_roman['data'].shape, image_roman['wcs'].proj_plane_pixel_scales())] +``` + +The field of view of a Roman exposure is ~450 arcsec, at a pixel scale of ~0.11 arcsec. + ++++ + +### Look at everything else the image file contains + +So far we have read the science image, but each exposure carries more than that. Listing the extensions shows the full picture. + +```{code-cell} ipython3 +with fits.open(f"s3://{image_roman['fpath']}", + fsspec_kwargs={"anon": True}) as hdul: + hdul.info() +``` + +Each file holds a primary header carrying the observation metadata but no pixels, followed by three 4088x4088 planes: + +SCI = the simulated science image, in electrons per second + +ERR = the estimated uncertainty on each science pixel + +DQ = the data quality mask, flagging pixels that should not be trusted + ++++ + +### Use the WCS from the Roman simulated header to convert the specified coordinate into a pixel position. + +```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- +def coord_to_xy(w, coord): + """Convert a sky position into (x, y) pixel indices. + + Parameters + ---------- + w : `~astropy.wcs.WCS` + WCS of the image the indices refer to. + coord : `~astropy.coordinates.SkyCoord` + Sky position to convert. + + Returns + ------- + tuple of int + Pixel indices as ``(x, y)``. + """ + return w.world_to_array_index(coord)[::-1] #reverse since 0th axis is y, 1st axis is x +``` + +```{code-cell} ipython3 +coord_arr_idx = coord_to_xy(image_roman['wcs'], coord) +coord_arr_idx +``` + +### Use matplotlib imshow to create a static visualization of the Roman simulated image and overplot the selected position. + +Each telescope points independently, so neither survey's images arrive with north pointing up: this Roman exposure is turned about 170 degrees from north, and the Rubin visit in the next section about 160 degrees. Every figure in this notebook is therefore resampled onto a north-up grid before plotting, so that all of them, and both telescopes, can be read the same way. The rotation leaves the frame tilted inside a slightly larger box, with empty wedges in the corners. + +```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- +def stretch_color(data, clipPercent): + """Percentile clip limits for displaying an image. + + Parameters + ---------- + data : `~numpy.ndarray` + Image values; NaNs are ignored. + clipPercent : float + Percentage clipped off each end of the distribution. + + Returns + ------- + tuple of float + Lower and upper limits, for passing to ``clim``. + """ + # nan-aware: rotating to north up leaves empty corners + return np.nanpercentile(data, (0 + clipPercent, 100 - clipPercent)) + + +def north_up(data, image_wcs, size, binning=1): + """Resample onto a north-up grid, keeping the original pixel scale. + + Parameters + ---------- + data : `~numpy.ndarray` + Image to resample. + image_wcs : `~astropy.wcs.WCS` + WCS of ``data``. + size : `~astropy.units.Quantity` + Angular width of the output grid. + binning : int, optional + Sample this many input pixels per output pixel. Values greater than 1 + are enough for displaying a whole exposure and much quicker than + resampling every pixel. + + Returns + ------- + `~numpy.ndarray` + The resampled image. + `~astropy.wcs.WCS` + WCS of the resampled image. + """ + scale = abs(image_wcs.proj_plane_pixel_scales()[0].to(u.deg).value) * binning + npix = int(round(size.to(u.deg).value / scale)) + + ny, nx = data.shape + center = image_wcs.pixel_to_world((nx - 1) / 2, (ny - 1) / 2) + + target = wcs.WCS(naxis=2) + target.wcs.ctype = ['RA---TAN', 'DEC--TAN'] + target.wcs.crpix = [(npix + 1) / 2, (npix + 1) / 2] + target.wcs.crval = [center.ra.deg, center.dec.deg] + target.wcs.cdelt = [-scale, scale] + + return reproject_interp((data, image_wcs), target, shape_out=(npix, npix))[0], target + + +def full_frame_north_up(image, binning=4): + """North-up view of a whole exposure, for display only. + + Parameters + ---------- + image : dict + Image dictionary with ``'data'`` and ``'wcs'`` keys, as returned by + `get_roman_image` or `get_rubin_image`. + binning : int, optional + Sample this many input pixels per output pixel. + + Returns + ------- + `~numpy.ndarray` + The resampled image. + `~astropy.wcs.WCS` + WCS of the resampled image. + """ + ny, nx = image['data'].shape + fov = nx * image['wcs'].proj_plane_pixel_scales()[0].to(u.arcsec) + return north_up(image['data'], image['wcs'], 1.2 * fov, binning=binning) +``` + +```{code-cell} ipython3 +roman_display, roman_display_wcs = full_frame_north_up(image_roman) + +plt.imshow(roman_display, origin='lower', + clim=stretch_color(roman_display, 1) + ) + +plt.plot(*coord_to_xy(roman_display_wcs, coord), 'o', mfc='none', mec='r', + markersize=15, markeredgewidth=2) +plt.axis('off') +``` + +## 3. Rubin Images + +OpenUniverse2024 includes simulated Rubin images in the following filters: u, g, r, i, z, y. As with Roman, these are stored by visit and detector rather than by sky position, so we again let the SIA service tell us which images cover our position. In this section, we define functions that retrieve a Rubin image for a chosen position and filter, returning data in the same structure as the function we defined above for Roman. + ++++ + +### Retrieve Rubin images + +The image search returns one row per image, with the cloud location of each tucked inside a `cloud_access` JSON string. Many visits cover any given position, so we sort by observation time and take the earliest. + +```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- +def get_rubin_image_fpaths(coord, filters): + """Find the Rubin images covering a position in each of several filters. + + Parameters + ---------- + coord : `~astropy.coordinates.SkyCoord` + Sky position to search around. + filters : list of str + Rubin bandpass names, e.g. ``['g', 'r', 'i']``. + + Returns + ------- + dict + Filter name to S3 path of the image covering ``coord``. + + Raises + ------ + ValueError + If any requested filter has no image covering ``coord``. + """ + # One search covers every band, so ask once and sort the results out by band + # afterwards rather than searching again for each one. + results = Irsa.query_sia(pos=(coord, SEARCH_RADIUS), + collection=OU_RUBIN_SIA_COLLECTION) + + fpaths = {} + for filter in filters: + # Band names come back like "r_57", so compare only the part before the underscore. + in_band = results[[str(r['energy_bandpassname']).split('_')[0] == filter + for r in results]] + if len(in_band) == 0: + raise ValueError(f"No Rubin {filter}-band images cover " + f"{coord.to_string('hmsdms')}") + + in_band.sort('t_min') + cloud_info = json.loads(in_band['cloud_access'][0])['aws'] + fpaths[filter] = f"{cloud_info['bucket_name']}/{cloud_info['key']}" + + return fpaths +``` + +```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- +def get_rubin_image(coord, filter): + """Read the Rubin image covering a position in a given filter. + + Parameters + ---------- + coord : `~astropy.coordinates.SkyCoord` + Sky position to search around. + filter : str + Rubin bandpass name, e.g. ``'r'``. + + Returns + ------- + dict + Keys ``'data'`` (`~numpy.ndarray`), ``'wcs'`` (`~astropy.wcs.WCS`) and + ``'fpath'`` (str, the S3 path of the image). + """ + image_s3_fpath = get_rubin_image_fpaths(coord, [filter])[filter] + + with fits.open(f"s3://{image_s3_fpath}", fsspec_kwargs={"anon": True}) as hdul: + # retrieve science data, which sits in the first extension + image_data = hdul[1].data + + # make wcs using header + image_wcs = wcs.WCS(hdul[1].header) + + # hand back the path too, so later cells can reuse it instead of + # asking the image search for it a second time + return {'data': image_data, 'wcs': image_wcs, 'fpath': image_s3_fpath} +``` + +### Inspect a simulated Rubin image + ++++ + +Choose a filter and retrieve the data and header information from the simulated Rubin image covering our position. + +```{code-cell} ipython3 +filter_rubin = 'r' +image_rubin = get_rubin_image(coord, filter_rubin) +``` + +### Understand the size of a simulated Rubin image. + +```{code-cell} ipython3 +# Number of pixels (Y, X) +image_rubin['data'].shape +``` + +```{code-cell} ipython3 +# Pixel size (scale Y, scale X) [degrees/pixel] +image_rubin['wcs'].proj_plane_pixel_scales() +``` + +```{code-cell} ipython3 +# Image size (FOV Y, FOV X) +[(num * size).to('arcsec') for num, size in zip( + image_rubin['data'].shape, image_rubin['wcs'].proj_plane_pixel_scales())] +``` + +A single Rubin detector covers a somewhat larger patch of sky than one Roman exposure, at roughly twice the pixel scale. + ++++ + +### Use matplotlib imshow to create a static visualization of the Rubin simulated image and overplot the selected position. + +The two full-frame figures do not cover the same patch of sky. Each is centered on its own exposure and sized to its own field of view, so the Rubin panel is nearly twice as wide on the sky as the Roman one. What they share is the marked position, which falls inside both. + +```{code-cell} ipython3 +rubin_display, rubin_display_wcs = full_frame_north_up(image_rubin) + +plt.imshow(rubin_display, origin='lower', + clim=stretch_color(rubin_display, 1) + ) + +plt.plot(*coord_to_xy(rubin_display_wcs, coord), 'o', mfc='none', mec='r', + markersize=15, markeredgewidth=2) +plt.axis('off') +``` + +## 4. Compare simulated Roman and Rubin cutouts for a selected position + ++++ + +### Choose cutout size + +```{code-cell} ipython3 +cutout_size = 50*u.arcsec +``` + +### Create the cutouts + +```{code-cell} ipython3 +# Cut oversized: the rotation below needs data in the corners to draw from. +cutout_roman = Cutout2D(image_roman['data'], coord, size=1.5*cutout_size, wcs=image_roman['wcs']) +cutout_rubin = Cutout2D(image_rubin['data'], coord, size=1.5*cutout_size, wcs=image_rubin['wcs']) +``` + +### Rotate both cutouts so that north points up + +The full frames above were rotated for display; here we do the same to the cutouts, using the same helper. The two telescopes are turned about 10 degrees apart from each other, so without this the panels would not line up. Each keeps its own pixel scale, so Roman still shows the finer sampling. + +```{code-cell} ipython3 +# Rotating a same-size grid would leave empty corners, so each cutout above was cut +# oversized and is resampled down to the size we actually want. +roman_data, roman_wcs = north_up(cutout_roman.data, cutout_roman.wcs, cutout_size) +rubin_data, rubin_wcs = north_up(cutout_rubin.data, cutout_rubin.wcs, cutout_size) +``` + +### Use matplotlib imshow to plot static side-by-side comparisons of the cutouts + +```{code-cell} ipython3 +fig, axs = plt.subplots(1, 2, figsize=(12, 6)) + + +axs[0].imshow(roman_data, origin='lower', + clim=stretch_color(roman_data, .5) + ) +axs[0].set_title(f"ROMAN in filter {filter_roman}") + +axs[1].imshow(rubin_data, origin='lower', + clim=stretch_color(rubin_data, .5) + ) +axs[1].set_title(f"RUBIN in filter {filter_rubin}") + +fig.suptitle(f"Cutouts at ({coord.ra}, {coord.dec}) with {cutout_size} size", fontsize=14) +plt.tight_layout(rect=[0, 0, 1, 0.97]) +``` + +Both panels cover the same patch of sky, at the same size and with north up in each. The difference is stark: the simulated Roman exposure resolves a field full of galaxies, while the same patch in the simulated Rubin image is close to blank. A single Rubin visit is short, so only the brightest sources here rise above its noise, and from the ground its point spread function is broad enough to blend close neighbors together. Roman's finer pixels and sharper point spread function keep those sources separate. Rubin reaches its full depth by stacking many visits, which is not what we are looking at. + +## 5. Use Firefly to interactively explore the images + +Static cutouts only help once you know where to look. Choosing a position worth comparing means exploring the image interactively, and for that we hand it to Firefly. In this section we display the simulated Rubin visit, outline the area the Roman exposure covers, and pan and zoom around it to pick out a target. + ++++ + +### Launch and initialize Firefly +There are two ways to initialize a Firefly client from Python, depending on whether you're running the notebook in JupyterLab or not. Assuming you have `jupyter-firefly-extensions` set up in your environment as explained [here](https://github.com/Caltech-IPAC/jupyter_firefly_extensions/blob/master/README.md), you can use `make_lab_client()` in JupyterLab, which will open the Firefly viewer in a new tab within the Lab. Otherwise, you can use `make_client()` in a Jupyter Notebook (or even a Python shell), which will open the Firefly viewer in a new web browser tab. + +You also need a Firefly server to communicate with your Firefly Python client. In this notebook, we use a public Firefly server: the IRSA Viewer (https://irsa.ipac.caltech.edu/irsaviewer). However, you can also run a local Firefly server via a [Firefly Docker image](https://hub.docker.com/r/ipac/firefly) and access it at `http://localhost:8080/firefly`. The URL of the Firefly server is read by both `make_client()` and `make_lab_client()` through the environment variable `FIREFLY_URL`. However, `make_client()` also allows you to pass the URL directly as the `url` parameter. + +```{code-cell} ipython3 +# Uncomment when using within Jupyter Lab with jupyter_firefly_extensions installed +# fc = FireflyClient.make_lab_client() + +# Uncomment for contexts other than above +fc = FireflyClient.make_client(url="https://irsa.ipac.caltech.edu/irsaviewer") + +fc.reinit_viewer() # to clean the state, if this cell ran earlier +``` + +### Define a function that returns the URL for a given S3 filepath +Since the OpenUniverse2024 data are in a public S3 bucket, any file can also be reached over HTTPS. Firefly loads images and catalogs by URL rather than by S3 path, so we use this helper throughout the rest of the notebook. + +```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- +def https_url(s3_fpath): + """Convert an S3 path into a public HTTPS URL. + + Parameters + ---------- + s3_fpath : str + Path of the form ``'bucket/key'``. + + Returns + ------- + str + HTTPS URL for the same file. + """ + s3_fpath_without_bucket = s3_fpath.split('/', 1)[1] + return f"https://{BUCKET_NAME}.s3.amazonaws.com/{s3_fpath_without_bucket}" +``` + +Let's generate URL for the Rubin image we plotted above. Clicking on the returned URL will allow you to download this image locally. + +```{code-cell} ipython3 +image_s3_fpath_rubin = image_rubin['fpath'] +https_url(image_s3_fpath_rubin) +``` + +### Send the simulated Rubin image to Firefly using show_fits. + +For displaying the FITS image of the Rubin visit in Firefly, we use [`show_fits`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.show_fits): + +```{code-cell} ipython3 +image_ff_id_rubin = 'rubin-image-filter-r' +fc.show_fits(url=https_url(image_s3_fpath_rubin), + plot_id=image_ff_id_rubin, + Title="Rubin Image" + ) +``` + +### Use ds9 region syntax to outline the Roman exposure on the interactive display + +The Firefly client includes several methods related to controlling ds9 region overlays. To +overlay a region layer on the loaded FITS images, we can use [`overlay_region_layer`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.overlay_region_layer). + +Region data is defined in ds9 region syntax that can be found [here](https://ds9.si.edu/doc/ref/region.html). + +The Rubin visit on display covers more sky than the Roman exposure does, so it is worth seeing where the Roman data actually falls. The corners of the Roman exposure come straight from its WCS, and we draw them as a polygon. + +```{code-cell} ipython3 +roman_corners = image_roman['wcs'].calc_footprint() +roman_corners +``` + +```{code-cell} ipython3 +# outline the Roman exposure as a polygon +corner_coords = ' '.join(f'{ra}d {dec}d' for ra, dec in roman_corners) +roman_regions = [f'icrs;polygon {corner_coords} # text={{Roman {filter_roman}}}'] + +roman_regions_id = 'roman_regions' +# Name the plot explicitly. Without plot_id the layer goes to whichever plot is +# active, which is not necessarily this one if cells are run out of order. +fc.overlay_region_layer(region_data=roman_regions, + title='Roman exposure', + region_layer_id=roman_regions_id, + plot_id=image_ff_id_rubin) +``` + +### Use Firefly's pan and zoom capabilities to locate a region of interest + +You can view the coordinates of your mouse pointer at the bottom left of the display window. To copy the coordinates for a specific coordinate: + +- Toggle the "Click Lock" to "on" in the bottom right of the image display. +- Click on the position of interest, and notice that the coordinate display is now frozen. +- Click on "EQ-J2000", the coordinate label in the bottom left of the image display. The pop-up that opens lets you choose the coordinate system, and the format used when copying — pick the Python format so the value can be pasted straight into a `SkyCoord`. +- Close the pop-up and click the small clipboard icon next to the coordinate readout. + ++++ + +### Copy the coordinates from the coordinate display to the Python notebook + +The position we have been using since Section 2 is the one marked below, so by default we carry it forward. Paste a position of your own here instead to compare a different target. + +```{code-cell} ipython3 +coords_of_interest = coord # or e.g. SkyCoord('0h38m28.92s -44d03m51.01s', frame='icrs') +coords_of_interest +``` + +We can now use this [astropy `SkyCoord` object](https://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html#astropy.coordinates.SkyCoord) to compare our images. + ++++ + +### Use ds9 region syntax to overplot the selected position + +For this we use the id of the region layer we defined above, and add more region data using [`add_region_data`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.add_region_data). + +```{code-cell} ipython3 +point_region = f'icrs;point {coords_of_interest.ra.value}d {coords_of_interest.dec.value}d # point=cross 15 text={{Source of interest}}' +fc.add_region_data(region_data=point_region, region_layer_id=roman_regions_id, + plot_id=image_ff_id_rubin) +``` + +## 6. Plot cutouts of the identified source + +This is a small group of galaxies, with a close neighbour 2.7 arcsec away. Roman shows them as separate objects; in Rubin they sit at the level of the noise. The Rubin image is a single visit rather than a stack, so faint sources like these are simply not detected in it, which is worth keeping in mind when comparing the two panels below. + +```{code-cell} ipython3 +# The images from Sections 2 and 3 already cover the default position, so only +# re-read if you pasted a different one above. +if coords_of_interest.separation(coord) > 0*u.arcsec: + image_roman = get_roman_image(coords_of_interest, filter_roman) + image_rubin = get_rubin_image(coords_of_interest, filter_rubin) +``` + +```{code-cell} ipython3 +cutout_size = 8*u.arcsec +``` + +```{code-cell} ipython3 +# As in section 4, cut oversized and then rotate both to north up. +cutout_roman = Cutout2D(image_roman['data'], coords_of_interest, size=1.5*cutout_size, wcs=image_roman['wcs']) +cutout_rubin = Cutout2D(image_rubin['data'], coords_of_interest, size=1.5*cutout_size, wcs=image_rubin['wcs']) + +roman_data, roman_wcs = north_up(cutout_roman.data, cutout_roman.wcs, cutout_size) +rubin_data, rubin_wcs = north_up(cutout_rubin.data, cutout_rubin.wcs, cutout_size) +``` + +```{code-cell} ipython3 +fig, axs = plt.subplots(1, 2, figsize=(12, 6)) + +axs[0].imshow(roman_data, origin='lower', + clim=stretch_color(roman_data, .5) + ) +axs[0].set_title(f"ROMAN in filter {filter_roman}") + +# Let's also encircle the source we identified. The circle is 2 arcsec across +# on both panels, so the two cutouts can be compared directly despite their different +# pixel scales. +roman_scale = roman_wcs.proj_plane_pixel_scales()[0].to(u.arcsec).value +axs[0].add_patch(patches.Circle(coord_to_xy(roman_wcs, coords_of_interest), + radius=2/roman_scale, color='r', fill=False, linewidth=2)) + + +axs[1].imshow(rubin_data, origin='lower', + clim=stretch_color(rubin_data, .5) + ) +axs[1].set_title(f"RUBIN in filter {filter_rubin}") + +# Let's also encircle the source we identified +rubin_scale = rubin_wcs.proj_plane_pixel_scales()[0].to(u.arcsec).value +axs[1].add_patch(patches.Circle(coord_to_xy(rubin_wcs, coords_of_interest), + radius=2/rubin_scale, color='r', fill=False, linewidth=2)) + + +fig.suptitle(f"Cutouts at ({coords_of_interest.ra:6f}, {coords_of_interest.dec:6f}) with {cutout_size} size", fontsize=14); +plt.tight_layout(rect=[0, 0, 1, 0.97]) +# plt.savefig("plot.pdf", bbox_inches='tight', pad_inches=0.2) +``` + +## 7. Use Firefly to visualize the OpenUniverse2024 catalogs +Let's inspect the properties of sources in the Rubin image. For this we will use the input truth files present in S3 bucket. + +OpenUniverse2024 includes the input truth files that were used to create the simulated images. These files are in Parquet and HDF5 format, and include information about the properties of galaxies, stars, and transients. + +The catalogs are split by HEALPix sky region (nside=32, RING ordering) and the region index appears in each filename, so we can convert our sky position into that index with [`hpgeom`](https://hpgeom.readthedocs.io/en/latest/) rather than hunting through the directory listing. This is the same approach the [Quickstart](openuniverse2024_quickstart) tutorial takes. + +Each region comes as two separate catalogs: a point source catalog holding the stars, and a galaxy catalog. We load both below and overlay each on the image. + +```{code-cell} ipython3 +region = hpgeom.angle_to_pixel(32, coord.ra.deg, coord.dec.deg, lonlat=True, nest=False) +print(f"HEALPix region for our position: {region}") +``` + +```{code-cell} ipython3 +# Catalog table of star properties (in parquet format) +pointsource_cat_path = f"{BUCKET_NAME}/{TRUTH_FILES_PATH}/pointsource_{region}.parquet" +pointsource_cat_path +``` + +```{code-cell} ipython3 +# Catalog table of galaxy properties (in parquet format) +galaxy_cat_path = f"{BUCKET_NAME}/{TRUTH_FILES_PATH}/galaxy_{region}.parquet" +galaxy_cat_path +``` + +### Use Firefly's show_table to overlay the catalogs on the interactive image + +Each catalog spans a whole HEALPix region, which is much larger than the image we are looking at, so we filter the table down to the sky area the Rubin image actually covers. (Note: you can remove filters through the table UI if you wish to see the entire data) + +```{code-cell} ipython3 +# Bound the catalogs by the footprint of the Rubin image we sent to Firefly. +footprint = image_rubin['wcs'].calc_footprint() +ra_bounds = (footprint[:, 0].min(), footprint[:, 0].max()) +dec_bounds = (footprint[:, 1].min(), footprint[:, 1].max()) +``` + +You can visualize catalogs interactively with Firefly using [`show_table`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.show_table). This capability can take many parameters. Here we will simply send our catalog to Firefly so that we can (a) see an interactive table; (b) see this table plotted over the image that we've already sent; and (c) use the GUI to quickly create exploratory plots. See if you can use the GUI to quickly determine approximately how many galaxies cover the Rubin image and what the redshift distribution of these galaxies is. + +```{code-cell} ipython3 +cat_filters = [ + f'("ra" >= {ra_bounds[0]} AND "ra" <= {ra_bounds[1]})', + f'("dec" >= {dec_bounds[0]} AND "dec" <= {dec_bounds[1]})' +] +cat_filters +``` + +```{code-cell} ipython3 +fc.show_table(url=https_url(pointsource_cat_path), + title='Stars Catalog', + tbl_id='stars_cat', + filters=" AND ".join(cat_filters)) +``` + +```{code-cell} ipython3 +gal_cat_tbl_id = 'galaxy_cat' + +# may take ~1.25min, because galaxy catalog is a big file +fc.show_table(url=https_url(galaxy_cat_path), + title='Galaxy Catalog', + tbl_id=gal_cat_tbl_id, + filters=" AND ".join(cat_filters)) +``` + +For each row in the table you can notice a marker in the image. Selecting a row or marker changes the corresponding marker or row, respectively. You can click on "Details" tab in the UI to show properties of each source selected in image/table. + +### Use Firefly's apply_table_filters to show only high-redshift galaxies + +High redshift galaxies are the most interesting, so let's filter the table we sent to Firefly to only include z>2.5 galaxies. Notice how the table display and image overlay change. Notice how the chart becomes a scatterplot from a heatmap because the sources reduce. You can remove this filter or add new ones through the GUI. + +For filtering, we will use [`apply_table_filters`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.apply_table_filters) method on the galaxy table we loaded above. + +```{code-cell} ipython3 +fc.apply_table_filters(tbl_id=gal_cat_tbl_id, + filters=" AND ".join(cat_filters+['"redshift" > 2.5'])) +``` + +You can play with the filters directly from the UI as well. Try removing adding more filters in the tables and see how markers change. + ++++ + +### Use Firefly's show_fits_3color to create a 3 color image of the simulated Rubin images + +```{code-cell} ipython3 +# [R, G, B] +ROMAN_RGB_FILTERS = ['H158', 'J129', 'Y106'] +RUBIN_RGB_FILTERS = ['r', 'g', 'u'] +``` + +We already have a Rubin image with the catalog overlaid, let's make a 3 color image to see colors of marked objects more clearly. For this we will use [`show_fits_3color`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.show_fits_3color) method: + +```{code-cell} ipython3 +image_ff_id_rubin_3color = 'rubin-image-3color' +rubin_rgb_fpaths = get_rubin_image_fpaths(coord, RUBIN_RGB_FILTERS) +threeC = [ + dict(url=https_url(rubin_rgb_fpaths[filter_name]), + Title="Rubin 3 color") + for filter_name in RUBIN_RGB_FILTERS +] + +fc.show_fits_3color(three_color_params=threeC, + plot_id=image_ff_id_rubin_3color) +``` + +### Use Firefly's interactivity to identify a region of interest + +For example, we found a region of the sky that seems to have a high number of high redshift sources and we copy it from the image display: + +```{code-cell} ipython3 +# located and copied through UI +high_z_gal_coords = SkyCoord('0h38m28.62s -44d00m50.52s', frame='icrs') +high_z_gal_coords +``` + +```{code-cell} ipython3 +# let's also mark it in our region layer, so that it's easy to pinpoint later +point_region = f'icrs;point {high_z_gal_coords.ra.value}d {high_z_gal_coords.dec.value}d # point=cross 15 text={{z>2.5 mock galaxies}}' +fc.add_region_data(region_data=point_region, region_layer_id=roman_regions_id, + plot_id=image_ff_id_rubin) +``` + +## 8. Plot 3-color Roman image containing your region of interest +Let's inspect WCS of the Roman image first + +```{code-cell} ipython3 +image_roman['wcs'] +``` + +### Prepare Roman images for displaying in Firefly + +The Roman exposures are gzipped on S3, so rather than hand Firefly a URL as we did for Rubin, we read each one in Python and pass it up as an in-memory FITS file. + +```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- +def get_fits_stream(image_roman): + """Write an image to an in-memory FITS file for handing to Firefly. + + Parameters + ---------- + image_roman : dict + Image dictionary with ``'data'`` and ``'wcs'`` keys, as returned by + `get_roman_image`. + + Returns + ------- + `~io.BytesIO` + FITS file positioned at the start, ready to be read. + """ + # Create a FITS PrimaryHDU object with the image data + hdu = fits.PrimaryHDU(data=image_roman['data'], header=image_roman['wcs'].to_header()) + + # Write the HDU to the in-memory stream (to save I/O time) + fits_stream = BytesIO() + hdu.writeto(fits_stream, overwrite=True) + fits_stream.seek(0) # to bring reading pointer to the beginning of file + + return fits_stream +``` + +Then we read one exposure per RGB filter and turn each into a stream: + +```{code-cell} ipython3 +images_rgb = [] +images_rgb_fits_stream = [] + +for filter_name in ROMAN_RGB_FILTERS: + print(f'\nFILTER: {filter_name}') + print('Retrieving Roman image...') + image_roman = get_roman_image(high_z_gal_coords, filter_name) + images_rgb.append(image_roman) + + print('Writing to fits stream...') + images_rgb_fits_stream.append(get_fits_stream(image_roman)) +``` + +### Use Firefly's show_fits_3color to create a 3 color image of Roman images + +Now we upload each fits stream (in-memory fits file) to firefly using [`upload_fits_data()`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.upload_fits_data) and prepare color params to pass to the [`show_fits_3color()`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.show_fits_3color). + +```{code-cell} ipython3 +three_color_params = [ + { + 'file': fc.upload_fits_data(fits_stream), + 'Title': "Roman 3 color" + } for fits_stream in images_rgb_fits_stream] +``` + +```{code-cell} ipython3 +image_ff_id_roman_3color = 'roman-3color-high_z_gal' +fc.show_fits_3color(three_color_params=three_color_params, + plot_id=image_ff_id_roman_3color) +``` + +We can see a 3 color image of the Roman exposures containing the high-redshift galaxy sources. Try panning and zooming out; you can notice it covers a smaller patch of sky than the Rubin image. + +### Use Firefly's pan/zoom/align methods to locate high redshift sources +Now, let's pan & zoom to the region where we located high-redshift galaxy sources. Also align & lock all images being displayed by WCS. For these operations we use these 3 methods (click on them to see their documentation): +- [`set_pan`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.set_pan) +- [`set_zoom`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.set_zoom) +- [`align_images`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.align_images) + +```{code-cell} ipython3 +fc.set_pan(plot_id=image_ff_id_roman_3color, x=high_z_gal_coords.ra.deg, y=high_z_gal_coords.dec.deg, coord='j2000') +fc.set_zoom(plot_id=image_ff_id_roman_3color, factor=1) +fc.align_images(lock_match=True) +``` + +### Use Firefly's set_stretch method to change the stretch of the image display via Python + +The image has a lot of noise that obscures our high redshift sources of interest. You can use the Firefly GUI to change the stretch of the image display. We identify that squared stretch from -2 to 10 sigma highlights the colors of our sources better. You can also use the Firefly client's [`set_stretch`](https://caltech-ipac.github.io/firefly_client/api/firefly_client.FireflyClient.html#firefly_client.FireflyClient.set_stretch) to do this via Python. This is helpful for reproducibility and for scaling up to many images. + +```{code-cell} ipython3 +fc.set_stretch(plot_id=image_ff_id_roman_3color, stype='sigma', algorithm='squared', + band='ALL', lower_value=-2, upper_value=10) +``` + +*** + +## About This Notebook + +**Updated:** 2026-08-14 + +**Contact:** [the IRSA Helpdesk](https://irsa.ipac.caltech.edu/docs/help_desk.html) with questions or reporting problems. diff --git a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_quickstart.md b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_quickstart.md index 36c30c72..51f2877d 100644 --- a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_quickstart.md +++ b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_quickstart.md @@ -36,7 +36,7 @@ This tutorial is a focused introduction to data access only. It covers the three 2. **Parquet catalogs** — transient (SNANA), galaxy, and galaxy-flux tables, indexed by HEALPix sky region 3. **Image search via SIA** — querying which images cover a given sky position using astroquery and the IRSA Simple Image Access service -No astrophysical analysis is performed here. For science workflows that build on these access patterns, see the [TDE Light Curve](TDE_light_curve) and [SED Fitting](SED_fit) tutorials in this repository. +No astrophysical analysis is performed here. For science workflows that build on these access patterns, see the [TDE Light Curve](openuniverse2024_TDE_light_curve) and [SED Fitting](openuniverse2024_SED_fit) tutorials in this repository. ### Instructions @@ -77,7 +77,7 @@ from astropy.io import fits The OpenUniverse2024 data live on the cloud in a public AWS S3 bucket and can be accessed anonymously using `s3fs`. This section shows how to establish that connection, navigate the directory tree, and inspect the contents of a FITS image file. -In the path below, `simple_model` refers to the simulated images with noise and realistic instrument effects, as opposed to `truth` images which are noise-free. The `full` simulation covers the complete survey footprint; a smaller `preview` subset is also available. See the [OpenUniverse2024 paper](https://arxiv.org/abs/2501.05632) for details on the differences. A `pointing` is a unique Roman observation visit — each pointing corresponds to one placement of the 18-detector focal plane on the sky, producing up to 18 individual FITS files (one per detector). +In the path below, `simple_model` refers to the simulated images with noise and realistic instrument effects, as opposed to `truth` images which are noise-free. The `full` simulation covers the complete survey footprint. See the [OpenUniverse2024 paper](https://arxiv.org/abs/2501.05632) for details on how it was produced. A `pointing` is a unique Roman observation visit — each pointing corresponds to one placement of the 18-detector focal plane on the sky, producing up to 18 individual FITS files (one per detector). ```{code-cell} ipython3 # Create an anonymous (public read-only) connection to the NASA IRSA S3 bucket. @@ -130,6 +130,11 @@ Let's display a gallery of example images to get a sense of the data. Note this gallery can take about a minute to build. ```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- def show_gallery(files, max_images=9): """ Display a gallery of FITS images. @@ -205,10 +210,15 @@ print("Galaxy flux file: ", gal_flux_path) ### 2.1 Inspect the SNANA Transient Catalog -`inspect_parquet_columns()` reads only the Parquet metadata footer to print the row count and column names — no data is loaded into memory. +`inspect_parquet_files()` reads only the Parquet metadata footer to print the row count and column names — no data is loaded into memory. We use it here for the SNANA catalog and repeat it for the galaxy info and flux catalogs below. ```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- def inspect_parquet_files(s3_path, *, region='us-east-1'): """ Print the structure of a Parquet file on S3 without reading its data. @@ -305,6 +315,11 @@ OU_RUBIN_SIA_COLLECTION = 'simulated_rubin_openuniverse2024' ``` ```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- def get_s3_fpath(cloud_access): """Extract the S3 URI from the cloud_access JSON string in an SIA result.""" cloud_info = json.loads(cloud_access) diff --git a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_demodata_20000808_instances.csv b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_demodata_20000808_instances.csv deleted file mode 100644 index a6a3610b..00000000 --- a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_demodata_20000808_instances.csv +++ /dev/null @@ -1,40 +0,0 @@ -,filter,pointing,sca -14,R062,407,7 -19,R062,4257,10 -8,R062,4647,6 -6,R062,6182,5 -27,R062,8112,15 -31,R062,9652,17 -0,R062,11197,1 -1,R062,11582,1 -17,R062,15037,8 -23,R062,16962,11 -20,R062,18502,10 -9,R062,18507,6 -10,R062,20812,6 -37,R062,21592,18 -32,R062,21977,17 -28,R062,22357,15 -33,R062,23512,17 -2,R062,25442,1 -25,R062,27367,13 -15,R062,28512,7 -21,R062,32362,10 -11,R062,32752,6 -7,R062,34287,5 -29,R062,36217,15 -34,R062,37757,17 -3,R062,39302,1 -4,R062,39687,1 -18,R062,43142,8 -24,R062,45067,11 -22,R062,46607,10 -12,R062,46612,6 -13,R062,48917,6 -38,R062,49697,18 -35,R062,50082,17 -30,R062,50462,15 -36,R062,51617,17 -5,R062,53547,1 -26,R062,55472,13 -16,R062,56617,7 diff --git a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_timedomainsurvey.md b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_timedomainsurvey.md deleted file mode 100644 index 9c104976..00000000 --- a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_timedomainsurvey.md +++ /dev/null @@ -1,304 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.7 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 -authors: - - name: Lauren Aldoroty ---- - -# Analyzing cloud-hosted simulated Roman Time Domain Survey images - -+++ - -## Learning Goals: - -By the end of this tutorial, you will: - -1. learn more about the "observations" that make up the simulated Roman TDS preview. -2. learn how to find the locations of simulated supernovae in the preview data. -3. learn how to create aligned cutouts of simulated Roman images. -4. learn how to make an animated gif from these cutouts. - -+++ - -## Install and Import required modules - -```{code-cell} ipython3 -# Uncomment the next line to install dependencies if needed. -# !pip install astropy matplotlib numpy pandas pyarrow s3fs scipy -``` - -```{code-cell} ipython3 -# Import modules -import warnings - -import astropy.units as u -import matplotlib.pyplot as plt -import numpy as np -import pandas as pd -from astropy.coordinates import SkyCoord -from astropy.io import fits -from astropy.nddata import Cutout2D -from astropy.nddata.utils import NoOverlapError -from astropy.table import Table -from astropy.wcs import WCS, FITSFixedWarning -from matplotlib import animation -from scipy.ndimage import rotate - -# Needed to access data in the cloud -import s3fs -s3 = s3fs.S3FileSystem(anon=True) # create an S3 client - -# Filter out the FITSFixedWarning, which is consequenceless and gets thrown every time you deal with a WCS -# in a Roman openuniverse simulated image using astropy. -warnings.simplefilter('ignore', category=FITSFixedWarning) -``` - -## Define a module to get the date (mjd) of a particular pointing. - -```{code-cell} ipython3 -def get_mjd(pointing, - obseq_path=f's3://nasa-irsa-simulations/openuniverse2024/roman/preview/RomanTDS/Roman_TDS_obseq_11_6_23.fits'): - - """ - Retrieve MJD of a given pointing. - - :param pointing: Pointing ID. - :type pointing: int - :param obseq_path: Path to obseq file Roman_TDS_obseq_11_6_23.fits. - :type obseq_path: str, optional - :return: MJD of specified pointing. - :rtype: float - """ - - with fits.open(obseq_path, fsspec_kwargs={"anon": True}) as obs: - obseq = Table(obs[1].data) - mjd = float(obseq['date'][int(pointing)]) - - return mjd -``` - -## Define a module to create an animated gif from a collection of cutouts. - -```{code-cell} ipython3 -def animate_stamps(stamps, savepath, no_whitespace=True, - labels=[],labelxy=(0.05,0.95), - **kwargs): - """ - Make an animation of a sequence of image stamps. - - :param stamps: Must be in chronological order. - :type stamps: List of stamps from get_stamps or get_object_instances. - :param savepath: Path to save gif. - :type savepath: str - """ - - if no_whitespace: - with_whitespace = np.invert(np.any((np.isnan(np.array(stamps))), axis=(1,2))) # NOTE: Your first axis (first indexing value) should return one stamp. e.g. stamps[0] is the first stamp. - idx_whitespace = np.where(with_whitespace)[0] - stamps = np.array(stamps)[idx_whitespace] - if len(labels) != 0: - labels = np.array(labels)[idx_whitespace] - - fig, ax = plt.subplots(figsize=(5,5)) - fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=None, hspace=None) - plt.xticks([]) - plt.yticks([]) - - im = ax.imshow(stamps[0], animated=True) - - if len(labels) != 0: - txt = ax.text(labelxy[0],labelxy[1],labels[0],animated=True,color='white',transform=ax.transAxes,va='top',ha='left',**kwargs) - - def animate(i): - im.set_array(stamps[i]) - if len(labels) != 0: - txt.set_text(labels[i]) - - return [im] + [txt] - else: - return [im] - - writer = animation.PillowWriter() - anim = animation.FuncAnimation(fig, animate, interval=600, frames=len(stamps)) - anim.save(savepath, writer=writer) -``` - -## Read in the Observation Sequence File to learn more about the "observations" that make up the simulated Roman Time Domain Survey. - -```{code-cell} ipython3 -# Read in the (simulated) Observation Sequence File. - -BUCKET_NAME = 'nasa-irsa-simulations' -ROMAN_PREFIX = 'openuniverse2024/roman/preview' - -ROMAN_TDS_PATH = f'{ROMAN_PREFIX}/RomanTDS' -FILENAME = 'Roman_TDS_obseq_11_6_23.fits' -OBSEQ_PATH = f's3://{BUCKET_NAME}/{ROMAN_TDS_PATH}/{FILENAME}' - -obseq_hdu = fits.open(OBSEQ_PATH, fsspec_kwargs={"anon": True}) -obseq = pd.DataFrame(obseq_hdu[1].data) - -print(obseq) -``` - -## What is the spatial and temporal coverage of the openuniverse2024 Roman TDS data preview? - -```{code-cell} ipython3 -# Find the ranges of RA, Dec, and date listed in the observation sequence file. - -ra_min, dec_min = obseq[['ra','dec']].min() -ra_max, dec_max = obseq[['ra','dec']].max() -mjd_min = obseq['date'].min() -mjd_max = obseq['date'].max() - -print("ra_min, ra_max:", ra_min, ra_max) -print("mjd_min, mjd_max:", mjd_min, mjd_max) -``` - -## Read in the Supernova Analysis (SNANA) file. - -```{code-cell} ipython3 -parquet_file = f's3://{BUCKET_NAME}/{ROMAN_PREFIX}/roman_rubin_cats_v1.1.2_faint/snana_10307.parquet' -transients = pd.read_parquet(parquet_file, filesystem=s3) -``` - -## Let's find a relatively nearby SN Ia that lies within the region of the data preview. - -```{code-cell} ipython3 -#List the unique models in the SNANA file. -unique_models = pd.Series(transients['model_name']).drop_duplicates().tolist() -unique_models -``` - -```{code-cell} ipython3 -# Most of the models are non SNIa (NON1ASED). -# Choose only the SNIa -sn1a = transients[transients['model_name'] == 'SALT3.NIR_WAVEEXT'] # SNe Ia only. -print('Number of SN1a in SNANA file: ', len(sn1a)) -``` - -```{code-cell} ipython3 -# Choose the SNIa that overlap with the spatial extent of the OpenUniverse2024 Roman TDS data preview. -ra_mask = np.logical_and(sn1a['ra'] > ra_min, sn1a['ra'] < ra_max) -dec_mask = np.logical_and(sn1a['dec'] > dec_min, sn1a['dec'] < dec_max) -mjd_mask = np.logical_and(sn1a['start_mjd'] > mjd_min, sn1a['end_mjd'] < mjd_max) -all_mask = np.logical_and.reduce((ra_mask,dec_mask,mjd_mask)) -preview_sn1a = sn1a[all_mask] -print('Number of SNIa in OpenUniverse2024 data preview:', len(preview_sn1a)) -``` - -```{code-cell} ipython3 -# Choose the SNIa in the data preview that are nearby, at redshifts less than 0.7. -nearby_preview_sn1a = preview_sn1a[preview_sn1a['z_CMB'] < 0.7] -print('Number of nearby SNIa in OpenUniverse2024 data preview:', len(nearby_preview_sn1a)) -``` - -```{code-cell} ipython3 -# Let's choose SN 20000808. -oid = 20000808 -chosen_object = nearby_preview_sn1a[nearby_preview_sn1a['id'] == oid] -ra = chosen_object.get('ra') -dec = chosen_object.get('dec') -ra, dec = 9.619282, -44.313894 -coord = SkyCoord(ra*u.deg, dec*u.deg) -``` - -## Read in the auxiliary file that lists the simulated Roman TDS images covering the chosen SNIa. - -```{code-cell} ipython3 -# The auxiliary file contains all the images this thing is in. -# If you need to download this file, see https://irsa.ipac.caltech.edu/docs/notebooks/. -csvfile = './openuniverse2024_roman_demodata_20000808_instances.csv' -instances = pd.read_csv(csvfile, usecols=['filter','pointing','sca']) -instances -``` - -## Create cutouts of the chosen SNIa in the band of your choice. - -```{code-cell} ipython3 -band = 'R062' -instances = instances[instances['filter'] == band] -``` - -```{code-cell} ipython3 -#Make the cutouts; this will take a couple of minutes. -stamps = [] -mjd = [] -for i, row in enumerate(instances.itertuples()): - band, pointing, sca = row[1], row[2], row[3] - imgpath = f's3://{BUCKET_NAME}/{ROMAN_TDS_PATH}/images/simple_model/{band}/{pointing}/Roman_TDS_simple_model_{band}_{pointing}_{sca}.fits.gz' - print(imgpath) - with fits.open(imgpath, fsspec_kwargs={"anon": True}) as hdu: - img = hdu[1].data - header = hdu[0].header - wcs = WCS(header) - x, y = wcs.world_to_pixel(coord) - - # Manually rotate the images so they are all aligned. - CDmat = np.array([header['CD1_1'], header['CD1_2'], - header['CD2_1'], header['CD2_2']]).reshape(2,2) - - orientation = hdu[0].header['ORIENTAT'] - - # These chips are "flipped". - if sca % 3 == 0: - orientation += 180 - - # Build rotation matrix. - CD1_1_rot = np.cos(-orientation*np.pi/180) - CD1_2_rot = -np.sin(-orientation*np.pi/180) - CD2_1_rot = np.sin(-orientation*np.pi/180) - CD2_2_rot = np.cos(-orientation*np.pi/180) - - RotMat = np.array([CD1_1_rot, CD1_2_rot, - CD2_1_rot, CD2_2_rot]).reshape(2,2) - - RotMat_inv = np.array([CD1_1_rot, -CD1_2_rot, - -CD2_1_rot, CD2_2_rot]).reshape(2,2) - - # Apply rotation to the CDi_j header keywords. - CDmat_rot = np.dot(CDmat,RotMat_inv) - - # Update header. - header['CD1_1'], header['CD1_2'] = CDmat_rot[0] - header['CD2_1'], header['CD2_2'] = CDmat_rot[1] - header['ORIENTAT'] -= orientation - - # Rotate the image. - rot_img = rotate(img,angle=orientation,reshape=False,cval=np.nan) - hdu[1].data = rot_img - - rot_wcs = WCS(header) - - try: - # Make cutout around SN Ia location. - cutout = Cutout2D(rot_img,coord,100,wcs=rot_wcs,mode='partial') - stamps.append(cutout.data) - mjd.append(get_mjd(pointing)) - except NoOverlapError: - pass -``` - -## Make an animated gif out of the cutouts. - -```{code-cell} ipython3 -savepath = f'SN{oid}.gif' -savepath -animate_stamps(stamps, savepath, labels=mjd) -``` - -*** - -## About this notebook - -**Updated:** 2024-06-10 - -**Contact:** [the IRSA Helpdesk](https://irsa.ipac.caltech.edu/docs/help_desk.html) with questions or reporting problems. diff --git a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_wideareasurvey.md b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_wideareasurvey.md deleted file mode 100644 index 6e0544e2..00000000 --- a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_wideareasurvey.md +++ /dev/null @@ -1,174 +0,0 @@ ---- -jupytext: - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.16.2 -kernelspec: - display_name: science_demo - language: python - name: python3 -authors: - - name: Vandana Desai - - name: IPAC Science Platform team ---- - -# Analyzing cloud-hosted simulated Roman coadded images - -+++ - -## Learning Goals - -By the end of this tutorial, you will learn how to do the following: - -- Identify the row and column that contains a particular ra, dec coordinate. -- Browse the bucket containing the simulated Roman coadds. -- Identify a simulated Roman coadd by filter, row, and column. -- Take a closer look at a simulated coadd file. - -## Install and Import required modules - -```{code-cell} ipython3 -# Uncomment the next line to install dependencies if needed. -# !pip install astropy numpy s3fs -``` - -```{code-cell} ipython3 -#Import modules -from astropy.io import fits -import numpy as np -import s3fs # browse buckets -``` - -## Identify the row and column that contains a particular ra, dec coordinate. - -The full simulation covers a 1 deg x 1 deg area centered around RA = 9.5 deg and -Dec = -44.1 deg. This region is divided into 1296 blocks (36 rows and 36 columns), -each 100 arcsec across. - -The data preview presented here covers the central 20x20 arcmin, corresponding to 144 blocks (12 rows and 12 columns). - -```{code-cell} ipython3 -#Choose an RA, DEC of interest. -ra = 9.5981595 -dec = -44.2026950 - -#Centers of data preview blocks. Do not alter. -ra_block_centers = np.array([9.76330352298415, 9.724522605135252, 9.68574158906671, - 9.646960496603766, 9.608179349571955, 9.56939816979703, - 9.530616979104877, 9.491835799321422, 9.453054652272561, - 9.414273559784032, 9.375492543681393, 9.336711625789874]) -dec_block_centers = np.array([-44.252584927082495, -44.22480733304182, -44.197029724175756, - -44.16925210374898, -44.14147447502621, -44.11369684127218, - 44.08591920575162, -44.05814157172923, -44.03036394246976, - -44.0025863212379, -43.974808711298394, -43.94703111591591]) - -ra_difference_array = np.absolute(ra_block_centers-ra) -ra_block_centers_index = ra_difference_array.argmin() -closest_ra_center = ra_block_centers[ra_block_centers_index] -ra_dist = 3600. * ra_difference_array[ra_block_centers_index] -if ra_dist > 50: - print("Chosen ra not covered by OpenUniverse 2024 data preview simulated Roman coadds") -else: - COLUMN = ra_block_centers_index + 12 - print("COLUMN:", COLUMN) - print("") - -dec_difference_array = np.absolute(dec_block_centers-dec) -dec_block_centers_index = dec_difference_array.argmin() -closest_dec_center = dec_block_centers[dec_block_centers_index] -dec_dist = 3600. * dec_difference_array[dec_block_centers_index] - -if dec_dist > 50: - print("Chosen dec not covered by OpenUniverse 2024 data preview simulated Roman coadds") -else: - ROW = dec_block_centers_index + 12 - print("ROW:", ROW) -``` - -## Browse the bucket containing the simulated Roman coadds. - -```{code-cell} ipython3 -s3 = s3fs.S3FileSystem(anon=True) # create an S3 client - -BUCKET_NAME = "nasa-irsa-simulations" -ROMAN_PREFIX = "openuniverse2024/roman/preview" -COADD_PATH = f"{ROMAN_PREFIX}/RomanWAS/images/coadds" - -s3.ls(f"{BUCKET_NAME}/{COADD_PATH}") -``` - -## Identify a Roman simulated coadd by filter, row, and column. - -A simulated coadd can be uniquely identified by filter, row, and column. - -```{code-cell} ipython3 -#Choose a filter, row, and column -FILTER = 'H158' #Filters F184, H158, J129, K213, and Y106 are available in the data preview. -#ROW = 12 #Rows 12-23 are available in the data preview. -#COLUMN = 12 #Columns 12-23 are available in the data preview. - -#Construct the coadd filename from the chosen filter, row, and column. -filename_root = f"prod_{FILTER[0]}_{COLUMN}_{ROW}_map.fits" - -#Construct the full coadd path from the chosen filter, row, and column. -s3_uri = f"s3://{BUCKET_NAME}/{COADD_PATH}/{FILTER}/Row{ROW}/{filename_root}" - -#List this filename to make sure it is found. -s3.ls(s3_uri) -``` - -## Take a closer look at the simulated coadd file you identified. - -```{code-cell} ipython3 -#Show a summary of extensions for this file. - -with fits.open(s3_uri, fsspec_kwargs={"anon": True}) as hdul: - hdul.info() -``` - -The Primary HDU for the coadded image is a cube with 15 layers, i.e., its shape is 1x15x2688x2688. The layers are as follows: - -0 = simulated "Science" image (Roman+Rubin simulation, units of e/(0.11 arcsec)^2/exposure) - -1 = lab noise: based on dark frames from the April 2023 test, masked at 3 e/p/s. Units: e/(0.11 arcsec)^2/s - -2 = GalSim stars, on HEALPix resolution 14 grid, normalized to total flux of 1 - -3 = noisy stars, on HEALPix resolution 14 grid, normalized to total flux of 2.4e5 e with self-Poisson noise, including 86 e^2/input pixel background variance - -4 = stars, on HEALPix resolution 14 grid, total flux 1, but on in only one of the passes (to test transient response) - -5 = stars, on HEALPix resolution 14 grid, with total flux that varies by 5% from center to edge of the focal plane (to test what happens when the filter bandpass varies; 5% is highly exaggerated) - -6 = GalSim extended objects, on HEALPix resolution 14 grid, right now exponential profiles. The scale radius is log-distributed between 0.125 and 0.500 arcsec, and the ellipticity (g1,g2) is uniformly distributed in the disc of radius 0.5, i.e., g1^2+g2^2<0.5^2. - -7,8,9 = same objects as layer 6, but with applied shear of 0.02. The shear orientations are spaced by 60° in tangent vector space, so that in the (g1,g2)-space they are spaced by 120° and can be used for finite differences. Specifically, the directions are: layer 7 -> in East-West direction (shear PA = 270°). (g1,g2) = (0.02,0) layer 8 -> in NNW-SSE direction (shear PA = 330°). (g1,g2) = (-0.02/2,0.02√3/2) layer 9 -> in NNE-SSW direction (shear PA = 30°). (g1,g2) = (-0.02/2,-0.02√3/2) - -10 = coadded 1/f noise map, normalized to variance per ln f of 1 - -11,12,13,14 = coadded white noise maps, different seeds, normalized to variance of 1 in each input pixel - -The following HDUs contain additional information: - -CONFIG = the configuration file - -INDATA = the input images used, as a binary table. The columns are: obsid (int32) -> observation ID sca (int16) -> SCA (1 through 18, inclusive) ra (float64) -> right ascension of pointing center in degrees dec (float64) -> declination of pointing center in degrees pa (float64) -> position angle of pointing in degrees valid (logical) -> input science data file is valid (should be True) - -INWEIGHT = the mean input weights for how much each 1.25x1.25 arcsec postage stamp depends on each input exposure. The shape is 1 x Nin x 84 x 84, where Nin is the number of input images listed in INDATA. Note that each postage stamp is 32 output pixels, so 84x32=2688. - -If summed on axis 1, this will normally be something close to 1. Deviations of ~10% are common, due to plate scale variations and the normalization issues introduced by diffraction spikes. - -INWEIGHTFLAT = a reshape of INWEIGHT suitable for display as a single image in a viewer such as DS9. The contributions from the Nin input exposures are rearranged into a 1 x 84 x (N_in*84) array. - -FIDELITY, SIGMA, INWTSUM, EFFCOVER = maps of U_alpha/C, S_alpha, sum_i T_{alpha i}, and the effective coverage as rescaled int16 maps. See Rowe et al. (2011) for details on the definitions of these quantities. The comment in the 'UNIT' keyword indicates how to rescale these. - -+++ -*** - -## About this notebook - -**Updated:** 2024-06-10 - -**Contact:** [the IRSA Helpdesk](https://irsa.ipac.caltech.edu/docs/help_desk.html) with questions or reporting problems. diff --git a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_tds_supernova.md b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_tds_supernova.md new file mode 100644 index 00000000..deadbb43 --- /dev/null +++ b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_tds_supernova.md @@ -0,0 +1,376 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.16.7 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +authors: + - name: Lauren Aldoroty +--- + +# Following a simulated supernova through the Roman Time Domain Survey + ++++ + +## Learning Goals: + +By the end of this tutorial, you will: + +1. learn more about the "observations" that make up the simulated Roman Time Domain Survey (TDS). +2. learn how to find the locations of simulated supernovae in the transient catalog. +3. learn how to ask IRSA which simulated Roman images cover a given position and time. +4. learn how to create aligned cutouts of simulated Roman images. +5. learn how to make an animated gif from these cutouts. + ++++ + +## Introduction + +The Roman Time Domain Survey revisits the same patch of sky over and over, which is what makes it possible to watch a transient appear and fade. This notebook picks one simulated Type Ia supernova out of the OpenUniverse2024 transient catalog, collects every Roman image that covers it while it is bright, and stacks those images into a short movie. + +The survey stores its images by pointing and detector rather than by sky position, so which files contain a particular supernova is not something you can work out from the file paths. IRSA's Simple Image Access (SIA) service answers exactly that question, and we use it here to assemble the list of images to stack. + +If you are new to OpenUniverse2024, the [Quickstart](openuniverse2024_quickstart) tutorial introduces the directory layout, the parquet catalogs, and the image search used below. + ++++ + +## Install and Import required modules + +```{code-cell} ipython3 +# Uncomment the next line to install dependencies if needed. +# !pip install astropy matplotlib numpy pandas pyarrow s3fs scipy astroquery hpgeom +``` + +```{code-cell} ipython3 +# Import modules +import warnings +import json + +import astropy.units as u +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import hpgeom +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.nddata import Cutout2D +from astropy.nddata.utils import NoOverlapError +from astropy.wcs import WCS, FITSFixedWarning +from astroquery.ipac.irsa import Irsa +from matplotlib import animation +from scipy.ndimage import rotate +from IPython.display import HTML + +# Needed to access data in the cloud +import s3fs +s3 = s3fs.S3FileSystem(anon=True) # create an S3 client + +# Filter out the FITSFixedWarning, which is consequenceless and gets thrown every time you deal with a WCS +# in a Roman openuniverse simulated image using astropy. +warnings.simplefilter('ignore', category=FITSFixedWarning) +``` + +```{code-cell} ipython3 +# Point the astroquery IRSA client at the simulated-data services, which are +# separate from the ones serving IRSA's observed data. +Irsa.sia_url = "https://irsa.ipac.caltech.edu/simulated/SIA" +Irsa.tap_url = "https://irsa.ipac.caltech.edu/simulated/TAP" + +OU_ROMAN_SIA_COLLECTION = 'simulated_roman_openuniverse2024' +``` + +## Define a module to create an animated gif from a collection of cutouts. + +```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- +def animate_stamps(stamps, savepath, no_whitespace=True, + labels=[],labelxy=(0.05,0.95), + **kwargs): + """ + Make an animation of a sequence of image stamps. + + :param stamps: Must be in chronological order. + :type stamps: List of stamps from get_stamps or get_object_instances. + :param savepath: Path to save gif. + :type savepath: str + """ + + if no_whitespace: + with_whitespace = np.invert(np.any((np.isnan(np.array(stamps))), axis=(1,2))) # NOTE: Your first axis (first indexing value) should return one stamp. e.g. stamps[0] is the first stamp. + idx_whitespace = np.where(with_whitespace)[0] + stamps = np.array(stamps)[idx_whitespace] + if len(labels) != 0: + labels = np.array(labels)[idx_whitespace] + + fig, ax = plt.subplots(figsize=(5,5)) + fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=None, hspace=None) + plt.xticks([]) + plt.yticks([]) + + im = ax.imshow(stamps[0], animated=True) + + if len(labels) != 0: + txt = ax.text(labelxy[0],labelxy[1],labels[0],animated=True,color='white',transform=ax.transAxes,va='top',ha='left',**kwargs) + + def animate(i): + im.set_array(stamps[i]) + if len(labels) != 0: + txt.set_text(labels[i]) + + return [im] + [txt] + else: + return [im] + + writer = animation.PillowWriter() + anim = animation.FuncAnimation(fig, animate, interval=600, frames=len(stamps)) + anim.save(savepath, writer=writer) + + # Close the figure so the notebook does not also print a static copy of the first frame, + # and hand the animation back so it can be displayed with playback controls. + plt.close(fig) + return anim +``` + +## Read in the Observation Sequence File to learn more about the "observations" that make up the simulated Roman Time Domain Survey. + +```{code-cell} ipython3 +# Read in the (simulated) Observation Sequence File. + +BUCKET_NAME = 'nasa-irsa-simulations' +ROMAN_PREFIX = 'openuniverse2024/roman/full' + +ROMAN_TDS_PATH = f'{ROMAN_PREFIX}/RomanTDS' +FILENAME = 'Roman_TDS_obseq_11_6_23.fits' +OBSEQ_PATH = f's3://{BUCKET_NAME}/{ROMAN_TDS_PATH}/{FILENAME}' + +obseq_hdu = fits.open(OBSEQ_PATH, fsspec_kwargs={"anon": True}) +obseq = pd.DataFrame(obseq_hdu[1].data) + +print(obseq) +``` + +## What is the spatial and temporal coverage of the openuniverse2024 Roman TDS? + +```{code-cell} ipython3 +# Find the ranges of RA, Dec, and date listed in the observation sequence file. + +ra_min, dec_min = obseq[['ra','dec']].min() +ra_max, dec_max = obseq[['ra','dec']].max() +mjd_min = obseq['date'].min() +mjd_max = obseq['date'].max() + +print("ra_min, ra_max:", ra_min, ra_max) +print("mjd_min, mjd_max:", mjd_min, mjd_max) +``` + +## Read in the Supernova Analysis (SNANA) file. + +The transient catalogs are split by HEALPix sky region (nside=32, RING ordering), with the region index in the filename. We convert the center of the Roman TDS into that index rather than guessing at the file name. + +```{code-cell} ipython3 +# The Roman Time-Domain Survey is centered near the LSST ELAIS-S1 Deep Drilling Field. +region = hpgeom.angle_to_pixel(32, 9.45, -44.02, lonlat=True, nest=False) + +parquet_file = f's3://{BUCKET_NAME}/{ROMAN_PREFIX}/roman_rubin_cats_v1.1.2_faint/snana_{region}.parquet' +transients = pd.read_parquet(parquet_file, filesystem=s3) +``` + +## Let's find a relatively nearby SN Ia that the survey actually watched go off. + +```{code-cell} ipython3 +#List the unique models in the SNANA file. +unique_models = pd.Series(transients['model_name']).drop_duplicates().tolist() +unique_models +``` + +```{code-cell} ipython3 +# Most of the models are non SNIa (NON1ASED). +# Choose only the SNIa +sn1a = transients[transients['model_name'] == 'SALT3.NIR_WAVEEXT'] # SNe Ia only. +print('Number of SN1a in SNANA file: ', len(sn1a)) +``` + +```{code-cell} ipython3 +# Choose the SNIa that overlap with the spatial and temporal extent of the survey. +ra_mask = np.logical_and(sn1a['ra'] > ra_min, sn1a['ra'] < ra_max) +dec_mask = np.logical_and(sn1a['dec'] > dec_min, sn1a['dec'] < dec_max) +mjd_mask = np.logical_and(sn1a['start_mjd'] > mjd_min, sn1a['end_mjd'] < mjd_max) +all_mask = np.logical_and.reduce((ra_mask,dec_mask,mjd_mask)) +covered_sn1a = sn1a[all_mask] +print('Number of SNIa within the survey:', len(covered_sn1a)) +``` + +```{code-cell} ipython3 +# Choose the SNIa that are nearby, at redshifts less than 0.5, so they are bright enough +# to stand out clearly against their host galaxy. +nearby_sn1a = covered_sn1a[covered_sn1a['z_CMB'] < 0.5] +print('Number of nearby SNIa:', len(nearby_sn1a)) +``` + +A supernova is only worth animating if Roman happened to be looking at that patch of sky while it was bright. The catalog records the date each one peaks in `peak_mjd`, and the survey visits any given field in bursts rather than continuously, so we pick an object whose peak falls inside a well-visited stretch of the survey. + +The cuts above cannot check this for us: they confirm a supernova went off somewhere inside the survey's footprint and date range, not that the telescope was pointed at it at the time. If you change `oid` to explore a different object, compare the dates it was bright, `start_mjd` through `end_mjd`, against the dates Roman actually visited its position, which arrive as the `t_min` column of the image search a few cells below. Objects that pass every cut here but were only observed long before or after they erupted are common; for those, the epoch window a few cells below comes back empty and the notebook cannot build an animation at all. + +```{code-cell} ipython3 +# Let's choose SN 20131477, which peaks while its field is being visited regularly. +oid = 20131477 +chosen_object = nearby_sn1a[nearby_sn1a['id'] == oid].iloc[0] + +ra, dec = chosen_object['ra'], chosen_object['dec'] +peak_mjd = chosen_object['peak_mjd'] +coord = SkyCoord(ra*u.deg, dec*u.deg) + +print(f"SN {oid}: RA={ra:.6f}, Dec={dec:.6f}, z={chosen_object['z_CMB']:.3f}, peaks at MJD {peak_mjd:.1f}") +``` + +## Ask IRSA which simulated Roman images cover the chosen SNIa. + +We hand the position to IRSA's image search, which returns one row per image along with the time it was taken and where the file lives in the cloud. + +```{code-cell} ipython3 +sia_results = Irsa.query_sia(pos=(coord, 1 * u.arcsec), + collection=OU_ROMAN_SIA_COLLECTION) + +print(f"Images covering this position: {len(sia_results)}") +``` + +That covers every band and both Roman surveys, so we narrow it down to the Time Domain Survey images in a single band. + +```{code-cell} ipython3 +band = 'R062' + +is_tds = np.array(['TDS_simple_model' in str(obs_id) for obs_id in sia_results['obs_id']]) +is_band = np.char.strip(np.array(sia_results['energy_bandpassname'], dtype=str)) == band + +instances = sia_results[is_tds & is_band] +instances.sort('t_min') + +print(f"{band} images covering SN {oid}: {len(instances)}") +``` + +Finally we keep only the epochs around the peak. Frames from years before the explosion would add nothing to the movie beyond download time, so we take a window that starts shortly before the supernova appears and runs until it has faded. + +```{code-cell} ipython3 +epoch_mjd = np.asarray(instances['t_min'], dtype=float) +in_window = (epoch_mjd >= peak_mjd - 40) & (epoch_mjd <= peak_mjd + 80) + +instances = instances[in_window] +epoch_mjd = epoch_mjd[in_window] + +print(f"Epochs to animate: {len(instances)}, " + f"spanning MJD {epoch_mjd.min():.1f} to {epoch_mjd.max():.1f}") +``` + +The cloud location of each image arrives as a JSON string, which we unpack into an S3 path. + +```{code-cell} ipython3 +--- +jupyter: + source_hidden: true +tags: [hide-cell] +--- +def get_s3_fpath(cloud_access): + """Extract the S3 URI from the cloud_access JSON string in an image search result.""" + cloud_info = json.loads(cloud_access)['aws'] + return f"s3://{cloud_info['bucket_name']}/{cloud_info['key']}" +``` + +```{code-cell} ipython3 +image_paths = [get_s3_fpath(row['cloud_access']) for row in instances] +image_paths[:3] +``` + +## Create cutouts of the chosen SNIa. + +Roman's focal plane sits at a different angle on the sky at each visit, so the same patch of sky arrives rotated differently in every image. Before the frames can be stacked into a movie we rotate each one so that north points the same way throughout. The rotation angle is recorded in the header, and detectors in every third slot are mounted flipped, which we correct for as well. + +```{code-cell} ipython3 +#Make the cutouts; this will take a few minutes. +stamps = [] +mjd = [] +for imgpath, epoch in zip(image_paths, epoch_mjd): + print(imgpath) + with fits.open(imgpath, fsspec_kwargs={"anon": True}) as hdu: + img = hdu[1].data + header = hdu[0].header + wcs = WCS(header) + x, y = wcs.world_to_pixel(coord) + + # Manually rotate the images so they are all aligned. + CDmat = np.array([header['CD1_1'], header['CD1_2'], + header['CD2_1'], header['CD2_2']]).reshape(2,2) + + orientation = header['ORIENTAT'] + + # These chips are "flipped". + if header['SCA_NUM'] % 3 == 0: + orientation += 180 + + # Build rotation matrix. + CD1_1_rot = np.cos(-orientation*np.pi/180) + CD1_2_rot = -np.sin(-orientation*np.pi/180) + CD2_1_rot = np.sin(-orientation*np.pi/180) + CD2_2_rot = np.cos(-orientation*np.pi/180) + + RotMat = np.array([CD1_1_rot, CD1_2_rot, + CD2_1_rot, CD2_2_rot]).reshape(2,2) + + RotMat_inv = np.array([CD1_1_rot, -CD1_2_rot, + -CD2_1_rot, CD2_2_rot]).reshape(2,2) + + # Apply rotation to the CDi_j header keywords. + CDmat_rot = np.dot(CDmat,RotMat_inv) + + # Update header. + header['CD1_1'], header['CD1_2'] = CDmat_rot[0] + header['CD2_1'], header['CD2_2'] = CDmat_rot[1] + header['ORIENTAT'] -= orientation + + # Rotate the image. + rot_img = rotate(img,angle=orientation,reshape=False,cval=np.nan) + + rot_wcs = WCS(header) + + try: + # Make cutout around SN Ia location. + cutout = Cutout2D(rot_img,coord,100,wcs=rot_wcs,mode='partial') + stamps.append(cutout.data) + mjd.append(epoch) + except NoOverlapError: + pass + +print(f"Collected {len(stamps)} cutouts") +``` + +## Make an animated gif out of the cutouts. + +Watch the supernova brighten near the middle of the sequence and fade away again. + +```{code-cell} ipython3 +savepath = f'SN{oid}.gif' +savepath +anim = animate_stamps(stamps, savepath, labels=[f'MJD {m:.1f}' for m in mjd]) +``` + +The saved gif loops without stopping, which makes a brief brightening hard to follow. Displaying the animation instead gives playback controls: pause it, step one epoch at a time with the arrows, drag the slider to any frame, and read the date in the corner as you go. + +```{code-cell} ipython3 +HTML(anim.to_jshtml(default_mode='once')) +``` + +*** + +## About this notebook + +**Updated:** 2026-08-05 + +**Contact:** [the IRSA Helpdesk](https://irsa.ipac.caltech.edu/docs/help_desk.html) with questions or reporting problems. diff --git a/tutorials/simulated-data/simulated.md b/tutorials/simulated-data/simulated.md index c5ac243f..c3acc46b 100644 --- a/tutorials/simulated-data/simulated.md +++ b/tutorials/simulated-data/simulated.md @@ -5,13 +5,22 @@ Because this collection is heterogeneous in coverage, structure, and intended us Access methods are tailored to the structure and scale of each product. These tutorials are designed to help users get started with accessing, visualizing, and analyzing simulated datasets hosted at IRSA. +## OpenUniverse2024 + +[OpenUniverse2024](https://arxiv.org/abs/2501.05632) simulates the overlapping imaging surveys to be carried out by the Nancy Grace Roman Space Telescope and the Vera C. Rubin Observatory, covering roughly 70 square degrees of matched optical and infrared sky. +Start with the Quickstart tutorial, which introduces the directory layout, the parquet catalogs, and the image search that the other tutorials build on. + ```{notebook-gallery} notebook_metadata.yml -tutorials/simulated-data/roman_hlss_number_density.md -tutorials/simulated-data/cosmoDC2_TAP_access.md tutorials/simulated-data/OpenUniverse2024/openuniverse2024_quickstart.md -tutorials/simulated-data/OpenUniverse2024/OpenUniverse2024Preview_Firefly.md -tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_wideareasurvey.md +tutorials/simulated-data/OpenUniverse2024/openuniverse2024_firefly_visualization.md +tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_tds_supernova.md tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md -tutorials/simulated-data/OpenUniverse2024/openuniverse2024_roman_simulated_timedomainsurvey.md +``` + +## Other Simulated Datasets + +```{notebook-gallery} notebook_metadata.yml +tutorials/simulated-data/roman_hlss_number_density.md +tutorials/simulated-data/cosmoDC2_TAP_access.md ```