From f259c82ebfb5a029bb0141be64b57f833c5f33d9 Mon Sep 17 00:00:00 2001
From: Andreas Faisst <afaisst@ipac.caltech.edu>
Date: Fri, 8 May 2026 13:50:00 -0700
Subject: [PATCH 1/7] first version mosaic tutorial

---
 .../spherex/spherex_mosaic/spherex_mosaic.md  | 461 ++++++++++++++++++
 tutorials/spherex/spherex_psf.md              |   7 +-
 2 files changed, 466 insertions(+), 2 deletions(-)
 create mode 100644 tutorials/spherex/spherex_mosaic/spherex_mosaic.md

diff --git a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
new file mode 100644
index 00000000..2e6c0fad
--- /dev/null
+++ b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
@@ -0,0 +1,461 @@
+---
+jupyter:
+  jupytext:
+    cell_metadata_filter: -all
+    notebook_metadata_filter: all,-language_info
+    text_representation:
+      extension: .md
+      format_name: markdown
+      format_version: '1.3'
+      jupytext_version: 1.17.2
+  kernelspec:
+    display_name: Python 3 (ipykernel)
+    language: python
+    name: python3
+---
+
+---
+authors:
+- name: Andreas Faisst
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.18.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
+# Understanding and Analyzing SPHEREx Mosaic Cubes 
+
+
+## 1. Learning Goals
+
+* Open a SPHEREx mosaic cube from the *IRSA Mosaic Tool* and assign wavelengths to its planes
+* Measure a quick-view spectrum throught the cube using `PhotUtils` including background subtraction.
+* Create a PAH $3.3\,{\rm \mu m}$ emission line map from a plane in the cube.
+* Obtain the corresponding GALEX FUV image from IRSA and overlay the PAH $3.3\,{\rm \mu m}$ map.
+
+
+## 2. SPHEREx Overview
+
+SPHEREx is a NASA Astrophysics Medium Explorer mission that launched in March 2025.
+During its planned two-year mission, SPHEREx will obtain 0.75-5 micron spectroscopy over the entire sky, with deeper data in the SPHEREx Deep Fields.
+SPHEREx data will be used to:
+
+* **constrain the physics of inflation** by measuring its imprints on the three-dimensional large-scale distribution of matter,
+* **trace the history of galactic light production** through a deep multi-band measurement of large-scale clustering,
+* **investigate the abundance and composition of water and biogenic ices** in the early phases of star and planetary disk formation.
+
+The community will also mine SPHEREx data and combine it with synergistic data sets to address a variety of additional topics in astrophysics.
+
+More information is available in the [SPHEREx Explanatory Supplement](https://irsa.ipac.caltech.edu/data/SPHEREx/docs/SPHEREx_Expsupp_QR.pdf).
+
+
+## 3. Imports
+The following packages must be installed to run this notebook.
+
+```python
+# Uncomment the next line to install dependencies if needed.
+# !pip install astropy matplotlib numpy photutils reproject astroquery
+```
+
+```python
+import copy
+
+import numpy as np
+
+from astropy.io import fits
+from astropy.table import Table, vstack, hstack
+from astropy.stats import sigma_clipped_stats, SigmaClip
+from astropy.wcs import WCS
+from astropy.coordinates import SkyCoord
+from astropy.nddata import Cutout2D
+import astropy.units as u
+
+from photutils.aperture import CircularAperture, CircularAnnulus, ApertureStats, aperture_photometry
+
+from astroquery.ipac.irsa import Irsa
+
+from reproject import reproject_exact
+
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+
+
+## Plotting stuff
+def load_plotting_defaults():
+    mpl.rcParams['font.size'] = 14
+    mpl.rcParams['axes.labelpad'] = 7
+    mpl.rcParams['xtick.major.pad'] = 7
+    mpl.rcParams['ytick.major.pad'] = 7
+    mpl.rcParams['xtick.minor.visible'] = True
+    mpl.rcParams['ytick.minor.visible'] = True
+    mpl.rcParams['xtick.minor.top'] = True
+    mpl.rcParams['xtick.minor.bottom'] = True
+    mpl.rcParams['ytick.minor.left'] = True
+    mpl.rcParams['ytick.minor.right'] = True
+    mpl.rcParams['xtick.major.size'] = 5
+    mpl.rcParams['ytick.major.size'] = 5
+    mpl.rcParams['xtick.minor.size'] = 3
+    mpl.rcParams['ytick.minor.size'] = 3
+    mpl.rcParams['xtick.direction'] = 'in'
+    mpl.rcParams['ytick.direction'] = 'in'
+    #mpl.rc('text', usetex=True)
+    mpl.rc('font', family='serif')
+    mpl.rcParams['xtick.top'] = True
+    mpl.rcParams['ytick.right'] = True
+    mpl.rcParams['hatch.linewidth'] = 1.5
+    
+
+load_plotting_defaults()
+def_cols = plt.rcParams['axes.prop_cycle'].by_key()['color']
+```
+
+## Open SPHEREx Mosaic Cube
+We first show how to open the SPHEREx mosaic cube and how to assign wavelengths to the different planes in the cube.
+
+Because the [IRSA mosaic tool](https://irsa.ipac.caltech.edu/applications/spherex/tool-mosaic) does not have an API, we have created and downloaded the mosaic already and added it to the `./data/` directory. 
+
+```{tip}
+To retrieve the same mosaic as provided here, go to the [IRSA mosaic tool](https://irsa.ipac.caltech.edu/applications/spherex/tool-mosaic) and type in *M101* in the `Output Mosaic Center` text field. For the size of the mosaic, choose 30 arcminutes for both axis. For the output mosaic pixel scale choose 9 arcseconds. 
+```
+
+We first define the path to the mosaic cube.
+
+```python
+fn_spherex = "./data/M101_irsa.fits"
+```
+
+Next, we open the image. Note that we convert here the image from MJy/sr to mJy using the relation 
+
+```
+[mJy] = [MJy/sr] x 23.5045 x (pixelscale)^2 / 1E3
+```
+
+In addition, we also adjust the `BUNIT` keyword in the header to reflect the changes.
+
+```{warning}
+The header includes multiple WCS axes including a wavelength axis. In order to extract the spatial WCS (which is used, for example, for plotting or reprojection), we have to pass `astropy.wcs.WCS()` the header *and* the HDU list as well as specify the axes (the first two): `wcs = WCS(hdr, fobj=hdul, naxis=2)`. For more information, see the [astropy WCS documentation](https://docs.astropy.org/en/stable/api/astropy.wcs.WCS.html).
+In addition, some metadata on `NAXIS` is preserved in the header which gives an error later when we run the `reproject` algorithm. We therefore remove the wavelength axis here by using `hdr.pop('NAXIS3', None)`.
+```
+
+
+We can now open the cube taking into account the changes and issues we discussed above:
+
+```python
+with fits.open(fn_spherex) as hdul:
+    hdul.info()
+
+    # load wavelength table
+    wave_tab_tmp = Table( hdul['WCS-WAVE'].data )
+
+    # Load header
+    cube_hdr = hdul["IMAGE"].header
+
+    # calculate pixel scale
+    pixscale_mosaic = np.abs(cube_hdr["CDELT1"]*3600) # pixel scale in arcsec/px
+
+    # adjust image units and update header BUNIT
+    cube_img = hdul["IMAGE"].data * 23.5045 * (pixscale_mosaic)**2 / 1e3 # Mjy/sr -> mjy
+    cube_hdr['BUNIT'] = "mjy"
+    
+    # get spatial WCS
+    hdr2 = cube_hdr.copy()
+    hdr2["NAXIS"] = 2
+    hdr2.pop('NAXIS3', None) # remove wavelength axis
+    cube_wcs = WCS(hdr2 , fobj=hdul, naxis=2)
+    print(f"Loaded SPHEREx cube with pixel scale {pixscale_mosaic} arcsec/px")
+
+```
+
+Note that the FITS has different layers.
+
+* The `IMAGE` layer is the actual cube.
+* The `NHITS` layer includes the number of images that were combined in the mosaic for each pixel.
+* The `FLAGS` layer includes some useful flags to identify outlier and overflow pixels.
+* The `SPECTRAL_CHANNELS` layer summarizes some useful additional information on the spectral channels used to create the cube planes in a binary table.
+* The `WCS-WAVE` layer contains a binary table summarizing the wavelength information for each plane in the cube.
+
+Next we construct a handy table from the `WCS-WAVE` layer summarizing the wavelength information. This will also allow us to later add wavelengths to the cube planes.
+
+```python
+wave_tab = Table( [wave_tab_tmp["PLANE"][0],
+                  np.asarray([ww[0] for ww in wave_tab_tmp["WAVELENGTHS"][0] ]),
+                  np.asarray([ww[0] for ww in wave_tab_tmp["BANDWIDTHS"][0] ])],
+                  names=["plane","wavelengths","bandwidths"],
+                dtype=[int, float, float])
+```
+
+Let's have a look at this table. The `bandwidths` correspond to the width of each wavelength channel.
+
+```python
+wave_tab
+```
+
+<!-- #region -->
+## Create a Quick Look Spectrum
+
+Once we have the cube loaded, let's generate a quick-look spectrum. This has two purposes: First, it shows how to extract a spectrum from a cube using the `PhotUtils` python package. Second, it allows us to visualize the wavelength channels (i.e., planes) which we will need to chose the planes to make the final PAH $3.3\,{\rm \mu m}$ map.
+
+To obtain the spectrum, we create a function that computes the sum of the fluxes in apertures and does a background subtraction using an annulus. Both aperture and annulus sizes are user-defined. We use the `PhotUtils` package for this (see [here](https://photutils.readthedocs.io/en/latest/user_guide/aperture.html) for more information).
+The measurements are done for each cube plane and the combined.
+Furthermore, the function creates a plot showing the collapsed cube and the apertures used (this output can be turned off).
+
+
+```{tip}
+This function provides a very simple aperture photometry with background subtraction. The method can be extended. For example the function does not calculate photometric errors.  
+```
+<!-- #endregion -->
+
+```python
+def get_aperture_photometry(cube , r_aperture_px , r_inner_px , r_outer_px, MAKEPLOT):
+    '''
+    Create quick look spectrum from simple aperture photometry with background subtraction
+    from annulus.
+
+    PARAMETERS:
+    ===========
+    cube: input SPHEREx mosaic cube
+    r_aperture_px: Aperture for flux extraction in pixels
+    r_inner_px: Inner annulus bound for background calculation in pixels
+    r_outer_px: Outer annulus bound for background calculation in pixels
+    MAKEPLOT: Set to `True` if function should create a figure showing the collapsed
+              cube with overlay of aperture and annulli.
+
+    RETURNS:
+    ========
+    A table including the photometry results (sum of aperture flux, plane ID, etc).
+    
+    '''
+
+    def _helper(img , position, aperture, annulus_aperture): 
+        
+        ## Calculate background
+        sigclip = SigmaClip(sigma=3.0, maxiters=10)
+        aperstats = ApertureStats(img, annulus_aperture, sigma_clip=sigclip)
+        bkg_mean = aperstats.mean
+        aperture_area = aperture.area_overlap(img)
+        total_bkg = bkg_mean * aperture_area
+        
+        ## Get photometry and subtract background
+        #error = img * 0.3
+        phot_table = aperture_photometry(img, aperture)
+        phot_table["aperture_sum_bkgsub"] = phot_table["aperture_sum"] - total_bkg
+    
+        return(phot_table)
+
+    ## Define position and apertures
+    position = (cube.shape[1]//2,cube.shape[2]//2)
+    aperture = CircularAperture(position, r=r_aperture_px)
+    annulus_aperture = CircularAnnulus(position, r_in=r_inner_px, r_out=r_outer_px)
+
+    
+    ## Compute photometry (use helper function to iterate over planes)
+    phot_all = vstack( [_helper(cube[ii,:,:] , position, aperture, annulus_aperture) for ii in range(cube.shape[0])] ) # run all planes
+    phot_all["id"] = np.arange(cube.shape[0])+1 ## add plane numbers back
+    phot_all.rename_column("id","plane")
+
+    ## Plot if asked for
+    if MAKEPLOT:
+
+        fig = plt.figure(figsize=(5,5))
+        ax1 = fig.add_subplot(1,1,1)
+
+        # Plot the median image:
+        img_plot = np.nanmedian(cube, axis=0)
+        lims = np.nanpercentile(img_plot.flatten() , q=(5,99.9))
+        ax1.imshow(img_plot, origin="lower", vmin=lims[0], vmax=lims[1] , norm="log", cmap="inferno")
+
+        # Plot the center
+        ax1.plot(cube.shape[1]//2 , cube.shape[2]//2 , "+" , color="white")
+
+        # plot the apertures
+        aperture.plot(ax=ax1 , linestyle="-", color="white", label="Aperture")
+        annulus_aperture.plot(ax=ax1 , linestyle=":", color="white", label="Background")
+
+        ax1.legend(loc="lower right", fontsize=9, framealpha = 1, facecolor="black", labelcolor="white")
+
+        ax1.tick_params(which="both", axis="both", color="white", labelsize=11)
+
+        plt.show()
+
+    return(phot_all)
+```
+
+```python
+phot_table = get_aperture_photometry(cube = cube_img , r_aperture_px = 20, r_inner_px = 60, r_outer_px= 100, MAKEPLOT = True)
+phot_table
+```
+
+Finally, we also add the wavelength from the table we created above to the photometry table. This makes it easy to plot the spectrum afterwards.
+
+```python
+phot_table["wavelengths"] = wave_tab["wavelengths"].copy()
+```
+
+Now that we have all the information, we can finally plot the spectrum. In addition, we also add some prominent emission lines (sorted by wavelength in the legend) and indicate the cube plane numbers around the $3.3\,{\rm \mu m}$ PAH feature. These numbers will be needed later to create the emission line map.
+
+```python
+## Plot Spectrum
+fig = plt.figure(figsize=(7,4))
+ax1 = fig.add_subplot(1,1,1)
+
+# spectrum
+ax1.plot(phot_table["wavelengths"] , phot_table["aperture_sum_bkgsub"], "o", markersize=3, markerfacecolor="black", markeredgecolor="black")
+
+# some emission lines
+ax1.axvline(1.28, linestyle=":", color="gray", linewidth=0.5, label=r"Pa-$\beta$")
+ax1.axvline(1.6, linestyle=":", color="gray", linewidth=0.5, label=r"$1.6\,{\rm \mu m}$ bump")
+ax1.axvline(1.87, linestyle=":", color="gray", linewidth=0.5, label=r"Pa-$\alpha$")
+ax1.axvline(2.36, linestyle=":", color="gray", linewidth=0.5, label=r"Br-$\beta$")
+ax1.axvline(3.3, linestyle=":", color="gray", linewidth=0.5, label=r"PAH $3.3\,{\rm \mu m}$")
+ax1.axvline(4.05, linestyle=":", color="gray", linewidth=0.5, label=r"Br-$\alpha$")
+
+# add plane numbers around PAH 3.3
+sel = np.where( (phot_table["wavelengths"] > 2.6) & (phot_table["wavelengths"] < 3.8) & (~np.isnan(phot_table["aperture_sum_bkgsub"])) )[0]
+[ax1.text(phot_table["wavelengths"][ss] , phot_table["aperture_sum_bkgsub"][ss]*1.1, phot_table["plane"][ss] , fontsize=7, va="bottom", ha="center", rotation=90) for ss in sel]
+
+ax1.legend(fontsize=8)
+ax1.tick_params(which="both", axis="both", labelsize=12)
+ax1.set_xlabel(r"Wavelength [$\mu$m]", fontsize=12)
+ax1.set_ylabel(r"Flux [mJy]", fontsize=12)
+
+plt.show()
+```
+
+## Create a PAH $3.3\,{\rm \mu m}$ Map
+
+We now create the PAH $3.3\,{\rm \mu m}$ map from the cube. Using the spectrum plot above, we can identify the cube wavelength planes that include the PAH feature (plane 63) and the planes that include the continuum (planes 61/62 on the blue and 64/65 on the red). The emission line map is then simply created by subtracting the continuum planes from the plane containing the emission line.
+
+```{warning}
+Note that the plane IDs are 1-indexed (i.e. start from 1). However, the cube is an *ndarray*, which is 0-indexed. If we want to access wavelength plane 63, we will have to chose `cube[63-1,:,:]`.
+```
+
+For the cube creation, we set up a handy function:
+
+```python
+def make_map(cube , planes_feature , planes_continuum):
+    '''
+    Create SPHEREX map from planes. Note that the plane IDs are 1-indexed.
+
+    PARAMETERS:
+    ===========
+    cube: SPHEREx Mosaic cube
+    planes_feature: The cube planes including the emission line feature
+    planes_continuum: The cube planes including the continuum which will be subtracted
+    
+    '''
+    img_map = np.nansum( cube_img[np.asarray(planes_feature)-1 , :,:], axis=0 ) - np.nanmedian( cube_img[np.asarray(planes_continuum)-1 , :,:], axis=0 )
+    return(img_map)
+```
+
+And now we can easily extract the cube with the `planes_feature` and `planes_continuum` defined above.
+
+```python
+img_map = make_map(cube_img , planes_feature = [63] , planes_continuum = [61,62,64,65])
+```
+
+Finally, we can plot our PAH $3.3\,{\rm \mu m}$ SPHEREx map!
+
+```python
+fig = plt.figure(figsize=(5,5))
+ax1 = fig.add_subplot(1,1,1)
+lims = np.nanpercentile(img_map.flatten() , q=(5,99.5))
+ax1.imshow(img_map , origin="lower", vmin=lims[0], vmax=lims[1], norm="linear", cmap="inferno")
+ax1.tick_params(which="both", axis="both", color="white", labelsize=11)
+plt.show()
+```
+
+## Get Corresponding GALEX Image
+
+Next we retrieve the GALEX image corresponding to the SPHEREx mosaic. This is scientifically interesting as we can observe how the PAH emission correlates with UV emission of hot young stars. It is expected that there is an anticorrelation as the small PAH grain (traced by the $3.3\,{\rm \mu m}$ feature) are being dissociated in strong UV radiation fields.
+
+We first obtain the position in coordinates at the center of the SPHEREx cube. Note that this is a 3d cube, so we have to take axes 1 and 2 (0-indexed) to obtain the coordinates!
+
+```python
+radec = cube_wcs.all_pix2world(cube_img.shape[2]//2 , cube_img.shape[1]//2 , 0)
+print(f"Sky position at R.A. = {radec[0]} and Decl. = {radec[1]}")
+pos = SkyCoord(ra=radec[0], dec=radec[1], unit='deg')
+
+```
+
+We then query IRSA at this position to obtain the GALEX FUV image. This data product is part of the [Spitzer Local Volume Legacy Survey](https://irsa.ipac.caltech.edu/data/SPITZER/LVL/overview.html), which observes many local galaxies. Accordingly, we query the collection `spitzer_lvl` and search for `science` images from GALEX in the `FUV` energy passband.
+
+```python
+im_table = Irsa.query_sia(pos=(pos, 1 * u.arcsec), collection='spitzer_lvl')
+sel = np.where((im_table["dataproduct_subtype"] == "science") & (im_table["instrument_name"] == "GALEX") & (im_table["energy_bandpassname"] == "FUV") )[0]
+im_table[sel]
+```
+
+Once we have the table, we can access the URL of the image in the `access_url` column. With the URL we can load the FITS image and create a cutout of similar size as the SPHEREx mosaic.
+
+```{tip}
+Note that we also update the header of the image once we created the cutout. This will be necessary when we reproject the GALEX image onto the SPHEREx mosaic pixel scale.
+```
+
+```python
+image_url = im_table['access_url'][sel][0]
+size = u.Quantity([40,40], u.arcmin)
+
+with fits.open(image_url, memmap=False) as hdul:
+    hdul.info()
+    galex_cutout = Cutout2D(hdul[0].data , position=pos , size=size, wcs = WCS(hdul[0].header), mode="partial")
+    img_galex = galex_cutout.data
+    wcs_galex = galex_cutout.wcs
+    hdr_galex = hdul[0].header
+    hdr_galex.update(wcs_galex.to_header()) # update header to cutout.
+    hdr_galex["NAXIS1"] = img_galex.shape[1] # update header to cutout.
+    hdr_galex["NAXIS2"] = img_galex.shape[0] # update header to cutout.
+```
+
+Additionally, we do some post-processing on the image to convert the pixel units to milli-Jansky, similar as the SPHEREx mosaic. Note that the zero points for the GALEX FUV and NUV images are different as defined above. We refer to the [Spitzer LVL documentation](https://irsa.ipac.caltech.edu/data/SPITZER/LVL/doc/LVL_DR5_v5.pdf) for more details on the GALEX pixel units.
+
+```python
+zps = {"fuv":18.82, "nuv":20.08}
+img_galex_ab = -2.5*np.log10(img_galex) + zps["fuv"]
+img_galex = 10**(-0.4*(img_galex_ab - 23.9)) / 1e3 # AB -> ujy -> mjy
+```
+
+In a last step, we will have to reproject the GALEX image to the SPHEREx mosaic. This is necesary so we can overlap the GALEX image with the PAH $3.3\,{\rm \mu m}$ map that we have created above. For this, we use the `reproject` package.
+
+```python
+img_galex_rep, fp = reproject_exact(input_data = (img_galex , hdr_galex), output_projection = cube_wcs )
+img_galex_rep[img_galex_rep == 0] = np.nan
+img_galex_rep = img_galex_rep / np.nansum(img_galex_rep) * np.nansum(img_galex)
+```
+
+Finally, we caon overlay the PAH emission on the GALEX map, which is our final result of this tutorial notebook.
+
+```python
+fig = plt.figure(figsize=(7,7))
+ax1 = fig.add_subplot(1,1,1)
+
+# plot PAH map
+lims = np.nanpercentile(img_galex_rep.flatten() , q=(5,99.5))
+ax1.imshow(img_galex_rep , origin="lower", vmin=lims[0], vmax=lims[1], norm="linear", cmap="Blues")
+
+# Overplot contours
+ax1.contour(img_map, levels=[0.2,0.4,0.5], colors="red", linewidths=0.5)
+
+plt.show()
+```
+
+## Acknowledgements
+
+- [Caltech/IPAC-IRSA](https://irsa.ipac.caltech.edu/)
+
+## About this notebook
+
+**Updated:** 8 May 2026
+
+**Contact:** Contact [IRSA Helpdesk](https://irsa.ipac.caltech.edu/docs/help_desk.html) with questions or problems.
+
+**Runtime:** Approximately 30 seconds.
+
+```python
+
+```
diff --git a/tutorials/spherex/spherex_psf.md b/tutorials/spherex/spherex_psf.md
index 71984da3..fe2427da 100644
--- a/tutorials/spherex/spherex_psf.md
+++ b/tutorials/spherex/spherex_psf.md
@@ -3,14 +3,14 @@ authors:
 - name: Vandana Desai
 - name: Jessica Krick
 - name: Andreas Faisst
-- name: "Brigitta Sipőcz"
+- name: "Brigitta Sip\u0151cz"
 - name: Troy Raen
 jupytext:
   text_representation:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.18.1
+    jupytext_version: 1.17.2
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -637,3 +637,6 @@ To use this PSF for forward modeling or fitting, you must:
 
 **Runtime:** Approximately 30 seconds.
 
+```{code-cell} ipython3
+
+```

From 8754f081fbc2f51f1e9dcd39a91d2961f8714148 Mon Sep 17 00:00:00 2001
From: Andreas Faisst <afaisst@ipac.caltech.edu>
Date: Fri, 8 May 2026 13:53:12 -0700
Subject: [PATCH 2/7] updates on file name

---
 tutorials/spherex/spherex_mosaic/spherex_mosaic.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
index 2e6c0fad..a079080e 100644
--- a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
+++ b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
@@ -127,7 +127,7 @@ To retrieve the same mosaic as provided here, go to the [IRSA mosaic tool](https
 We first define the path to the mosaic cube.
 
 ```python
-fn_spherex = "./data/M101_irsa.fits"
+fn_spherex = "M101_irsa.fits"
 ```
 
 Next, we open the image. Note that we convert here the image from MJy/sr to mJy using the relation 

From 136885cb9214665b07e6a12f0785ffb1ee56441a Mon Sep 17 00:00:00 2001
From: Andreas Faisst <afaisst@ipac.caltech.edu>
Date: Wed, 13 May 2026 13:29:28 -0700
Subject: [PATCH 3/7] resolving comments part 1

---
 .../spherex/spherex_mosaic/spherex_mosaic.md  | 192 ++++++++++++------
 1 file changed, 130 insertions(+), 62 deletions(-)

diff --git a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
index a079080e..81970165 100644
--- a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
+++ b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
@@ -14,21 +14,6 @@ jupyter:
     name: python3
 ---
 
----
-authors:
-- name: Andreas Faisst
-jupytext:
-  text_representation:
-    extension: .md
-    format_name: myst
-    format_version: 0.13
-    jupytext_version: 1.18.1
-kernelspec:
-  display_name: Python 3 (ipykernel)
-  language: python
-  name: python3
----
-
 # Understanding and Analyzing SPHEREx Mosaic Cubes 
 
 
@@ -84,14 +69,11 @@ from reproject import reproject_exact
 
 import matplotlib.pyplot as plt
 import matplotlib as mpl
+```
 
-
-## Plotting stuff
+```python
+## Define some plotting formats
 def load_plotting_defaults():
-    mpl.rcParams['font.size'] = 14
-    mpl.rcParams['axes.labelpad'] = 7
-    mpl.rcParams['xtick.major.pad'] = 7
-    mpl.rcParams['ytick.major.pad'] = 7
     mpl.rcParams['xtick.minor.visible'] = True
     mpl.rcParams['ytick.minor.visible'] = True
     mpl.rcParams['xtick.minor.top'] = True
@@ -104,21 +86,19 @@ def load_plotting_defaults():
     mpl.rcParams['ytick.minor.size'] = 3
     mpl.rcParams['xtick.direction'] = 'in'
     mpl.rcParams['ytick.direction'] = 'in'
-    #mpl.rc('text', usetex=True)
-    mpl.rc('font', family='serif')
     mpl.rcParams['xtick.top'] = True
     mpl.rcParams['ytick.right'] = True
-    mpl.rcParams['hatch.linewidth'] = 1.5
-    
 
 load_plotting_defaults()
 def_cols = plt.rcParams['axes.prop_cycle'].by_key()['color']
 ```
 
-## Open SPHEREx Mosaic Cube
+## 4. Open SPHEREx Mosaic Cube
 We first show how to open the SPHEREx mosaic cube and how to assign wavelengths to the different planes in the cube.
 
-Because the [IRSA mosaic tool](https://irsa.ipac.caltech.edu/applications/spherex/tool-mosaic) does not have an API, we have created and downloaded the mosaic already and added it to the `./data/` directory. 
+The mosaic cube was obtained from the [IRSA SPHEREx mosaic tool](https://irsa.ipac.caltech.edu/applications/spherex/tool-mosaic) GUI. The mosaic tool computes a cube (x,y,$\lambda$) from SPHEREx data at a given sky position provided by the user. In brief, the tool gathers all the SPHEREx spectral images at that position, extracts the pixels at similar wavelenghts, and reprojects them into cube layers at the respective wavelengths. The final cube has 102 wavelength layers from $0.75$ to $5\,{\rm \mu m}$. The user can also specify the spatial size of the cube.
+
+Because the IRSA SPHEREx mosaic tool does not have an API, we have created and downloaded the mosaic already and added it to the `./data/` directory. 
 
 ```{tip}
 To retrieve the same mosaic as provided here, go to the [IRSA mosaic tool](https://irsa.ipac.caltech.edu/applications/spherex/tool-mosaic) and type in *M101* in the `Output Mosaic Center` text field. For the size of the mosaic, choose 30 arcminutes for both axis. For the output mosaic pixel scale choose 9 arcseconds. 
@@ -174,13 +154,23 @@ with fits.open(fn_spherex) as hdul:
 
 Note that the FITS has different layers.
 
-* The `IMAGE` layer is the actual cube.
+* The `IMAGE` layer is the actual cube where the 3 dimensions are (x, y, wavelength).
 * The `NHITS` layer includes the number of images that were combined in the mosaic for each pixel.
 * The `FLAGS` layer includes some useful flags to identify outlier and overflow pixels.
 * The `SPECTRAL_CHANNELS` layer summarizes some useful additional information on the spectral channels used to create the cube planes in a binary table.
 * The `WCS-WAVE` layer contains a binary table summarizing the wavelength information for each plane in the cube.
 
-Next we construct a handy table from the `WCS-WAVE` layer summarizing the wavelength information. This will also allow us to later add wavelengths to the cube planes.
+
+## 5. Create a Quick Look Spectrum
+
+Once we have the cube loaded, let's generate a quick-look spectrum. This has two purposes: First, it shows how to extract a spectrum from a cube using the `PhotUtils` python package. Second, it allows us to visualize the wavelength channels (i.e., planes) which we will need to chose the planes to make the final PAH $3.3\,{\rm \mu m}$ map.
+
+
+### 5.1 Create a Wavelength Look-Up Table
+
+Later we will measure an aperture flux on each plane (representing a different wavelength) of the mosaic cube. In order to turn this aperture flux table into a spectrum, we need to track the relevant wavelength range of each plane. Furthermore, the wavelength look-up table will be used to construct the emission line map.
+
+We construct the wavelenght look-up table from the `WCS-WAVE` layer, which summarizes the wavelength information.
 
 ```python
 wave_tab = Table( [wave_tab_tmp["PLANE"][0],
@@ -197,9 +187,7 @@ wave_tab
 ```
 
 <!-- #region -->
-## Create a Quick Look Spectrum
-
-Once we have the cube loaded, let's generate a quick-look spectrum. This has two purposes: First, it shows how to extract a spectrum from a cube using the `PhotUtils` python package. Second, it allows us to visualize the wavelength channels (i.e., planes) which we will need to chose the planes to make the final PAH $3.3\,{\rm \mu m}$ map.
+### 5.2 Extract Spectrum
 
 To obtain the spectrum, we create a function that computes the sum of the fluxes in apertures and does a background subtraction using an annulus. Both aperture and annulus sizes are user-defined. We use the `PhotUtils` package for this (see [here](https://photutils.readthedocs.io/en/latest/user_guide/aperture.html) for more information).
 The measurements are done for each cube plane and the combined.
@@ -212,27 +200,60 @@ This function provides a very simple aperture photometry with background subtrac
 <!-- #endregion -->
 
 ```python
-def get_aperture_photometry(cube , r_aperture_px , r_inner_px , r_outer_px, MAKEPLOT):
+def measure_aperture_photometry_cube(cube,
+                            r_aperture_px = 20,
+                            r_inner_px = 60,
+                            r_outer_px = 100,
+                            MAKEPLOT = True
+                           ):
     '''
     Create quick look spectrum from simple aperture photometry with background subtraction
     from annulus.
 
-    PARAMETERS:
-    ===========
-    cube: input SPHEREx mosaic cube
-    r_aperture_px: Aperture for flux extraction in pixels
-    r_inner_px: Inner annulus bound for background calculation in pixels
-    r_outer_px: Outer annulus bound for background calculation in pixels
-    MAKEPLOT: Set to `True` if function should create a figure showing the collapsed
-              cube with overlay of aperture and annulli.
-
-    RETURNS:
-    ========
-    A table including the photometry results (sum of aperture flux, plane ID, etc).
+    Parameters
+    ----------
+    cube : numpy.ndarray
+        Input SPHEREx mosaic cube.
+    r_aperture_px : float
+        Aperture for flux extraction in pixels.
+    r_inner_px : float
+        Inner annulus bound for background calculation in pixels.
+    r_outer_px : float
+        Outer annulus bound for background calculation in pixels.
+    MAKEPLOT : book
+        Set to `True` if function should create a figure showing the collapsed
+        cube with overlay of aperture and annulli.
+
+    Return
+    -------
+    astropy.QTable
+        A table including the photometry results (sum of aperture flux, plane ID, etc).
     
     '''
 
-    def _helper(img , position, aperture, annulus_aperture): 
+    ## Define helper function to compute the photometry efficiently
+    def apphot_helper(img , position, aperture, annulus_aperture):
+        '''
+        Helper function which computes the photometry for a single image.
+        Helper function will be called in series to obtain photometry for all images/planes.
+
+        Parameters
+        ----------
+        img : numpy.ndarray
+            Two-dimensional image on which aperture photometry is performed.
+        position : tuple of float
+            Center position of the aperture.
+        aperture : photutils.CircularAperture
+            Aperture for photometry extraction.
+        annulus_aperture : photutils.CircularAnnulus
+            Annulus aperture for background estimation.
+        
+        Returns
+        -------
+        astropy.QTable
+            A table including the photometry results (sum of aperture flux, plane ID, etc).
+        
+        '''
         
         ## Calculate background
         sigclip = SigmaClip(sigma=3.0, maxiters=10)
@@ -242,12 +263,12 @@ def get_aperture_photometry(cube , r_aperture_px , r_inner_px , r_outer_px, MAKE
         total_bkg = bkg_mean * aperture_area
         
         ## Get photometry and subtract background
-        #error = img * 0.3
         phot_table = aperture_photometry(img, aperture)
         phot_table["aperture_sum_bkgsub"] = phot_table["aperture_sum"] - total_bkg
     
         return(phot_table)
 
+    
     ## Define position and apertures
     position = (cube.shape[1]//2,cube.shape[2]//2)
     aperture = CircularAperture(position, r=r_aperture_px)
@@ -255,7 +276,7 @@ def get_aperture_photometry(cube , r_aperture_px , r_inner_px , r_outer_px, MAKE
 
     
     ## Compute photometry (use helper function to iterate over planes)
-    phot_all = vstack( [_helper(cube[ii,:,:] , position, aperture, annulus_aperture) for ii in range(cube.shape[0])] ) # run all planes
+    phot_all = vstack( [apphot_helper(cube[ii,:,:] , position, aperture, annulus_aperture) for ii in range(cube.shape[0])] ) # run all planes
     phot_all["id"] = np.arange(cube.shape[0])+1 ## add plane numbers back
     phot_all.rename_column("id","plane")
 
@@ -281,17 +302,24 @@ def get_aperture_photometry(cube , r_aperture_px , r_inner_px , r_outer_px, MAKE
 
         ax1.tick_params(which="both", axis="both", color="white", labelsize=11)
 
+        ax1.set_xlabel(r"$x$ [px]")
+        ax1.set_ylabel(r"$y$ [px]")
+
         plt.show()
 
     return(phot_all)
 ```
 
 ```python
-phot_table = get_aperture_photometry(cube = cube_img , r_aperture_px = 20, r_inner_px = 60, r_outer_px= 100, MAKEPLOT = True)
+phot_table = measure_aperture_photometry_cube(cube = cube_img , r_aperture_px = 20, r_inner_px = 60, r_outer_px= 100, MAKEPLOT = True)
 phot_table
 ```
 
-Finally, we also add the wavelength from the table we created above to the photometry table. This makes it easy to plot the spectrum afterwards.
+```{note}
+Note that are some NaN values in that table. This is because of missing data for these specific planes. As the SPHEREx mission progresses, these missing data will be filled in.
+```
+
+Finally, we also add the wavelength from the wavelength look-up table that created above to the photometry table. This makes it easy to plot the spectrum afterwards.
 
 ```python
 phot_table["wavelengths"] = wave_tab["wavelengths"].copy()
@@ -327,7 +355,7 @@ ax1.set_ylabel(r"Flux [mJy]", fontsize=12)
 plt.show()
 ```
 
-## Create a PAH $3.3\,{\rm \mu m}$ Map
+## 6. Create a PAH $3.3\,{\rm \mu m}$ Map
 
 We now create the PAH $3.3\,{\rm \mu m}$ map from the cube. Using the spectrum plot above, we can identify the cube wavelength planes that include the PAH feature (plane 63) and the planes that include the continuum (planes 61/62 on the blue and 64/65 on the red). The emission line map is then simply created by subtracting the continuum planes from the plane containing the emission line.
 
@@ -338,18 +366,29 @@ Note that the plane IDs are 1-indexed (i.e. start from 1). However, the cube is
 For the cube creation, we set up a handy function:
 
 ```python
-def make_map(cube , planes_feature , planes_continuum):
+def make_map(cube,
+             planes_feature = [63],
+             planes_continuum = [61,62,64,65]
+            ):
     '''
     Create SPHEREX map from planes. Note that the plane IDs are 1-indexed.
 
-    PARAMETERS:
-    ===========
-    cube: SPHEREx Mosaic cube
-    planes_feature: The cube planes including the emission line feature
-    planes_continuum: The cube planes including the continuum which will be subtracted
+    Parameters
+    ----------
+    cube : numpy.ndarray
+        SPHEREx Mosaic cube.
+    planes_feature : list of int
+        The cube planes including the emission line feature.
+    planes_continuum : list of int
+        The cube planes including the continuum which will be subtracted.
+
+    Returns
+    -------
+    numpy.ndarray
+        Two-dimensional map.
     
     '''
-    img_map = np.nansum( cube_img[np.asarray(planes_feature)-1 , :,:], axis=0 ) - np.nanmedian( cube_img[np.asarray(planes_continuum)-1 , :,:], axis=0 )
+    img_map = np.nansum( cube[np.asarray(planes_feature)-1 , :,:], axis=0 ) - np.nanmedian( cube[np.asarray(planes_continuum)-1 , :,:], axis=0 )
     return(img_map)
 ```
 
@@ -367,10 +406,15 @@ ax1 = fig.add_subplot(1,1,1)
 lims = np.nanpercentile(img_map.flatten() , q=(5,99.5))
 ax1.imshow(img_map , origin="lower", vmin=lims[0], vmax=lims[1], norm="linear", cmap="inferno")
 ax1.tick_params(which="both", axis="both", color="white", labelsize=11)
+ax1.set_xlabel(r"$x$ [px]")
+ax1.set_ylabel(r"$y$ [px]")
 plt.show()
 ```
 
-## Get Corresponding GALEX Image
+The map shows the location of the PAH $3.3\,{\rm \mu m}$ emission in the local galaxy M101. Specifically, the $3.3\,{\rm \mu m}$ feature is an observational indicator of very small carbonageous dust grains. You can see that the PAH emission is not continuous as it is tied to regions of star formation. To explore this further, we can correlate the PAH map with the GALEX far-UV (FUV) continuum map, which maps the UV emission of hot, young stars. This is done in the next section.
+
+
+## 7. Get Corresponding GALEX Image
 
 Next we retrieve the GALEX image corresponding to the SPHEREx mosaic. This is scientifically interesting as we can observe how the PAH emission correlates with UV emission of hot young stars. It is expected that there is an anticorrelation as the small PAH grain (traced by the $3.3\,{\rm \mu m}$ feature) are being dissociated in strong UV radiation fields.
 
@@ -431,19 +475,43 @@ img_galex_rep = img_galex_rep / np.nansum(img_galex_rep) * np.nansum(img_galex)
 Finally, we caon overlay the PAH emission on the GALEX map, which is our final result of this tutorial notebook.
 
 ```python
-fig = plt.figure(figsize=(7,7))
+fig = plt.figure(figsize=(9,9))
 ax1 = fig.add_subplot(1,1,1)
 
+## Main figure -------
+
 # plot PAH map
 lims = np.nanpercentile(img_galex_rep.flatten() , q=(5,99.5))
 ax1.imshow(img_galex_rep , origin="lower", vmin=lims[0], vmax=lims[1], norm="linear", cmap="Blues")
 
 # Overplot contours
-ax1.contour(img_map, levels=[0.2,0.4,0.5], colors="red", linewidths=0.5)
+ax1.contour(img_map, levels=[0.3,0.4,0.5], colors="red", linewidths=0.5)
+
+ax1.set_xlabel(r"$x$ [px]")
+ax1.set_ylabel(r"$y$ [px]")
+
+
+## Inset (zoom-in) -----
+axins = ax1.inset_axes([0.5,0.05,0.5,0.3])
+
+center = [125, 120] # center in pixels
+width = [70,50] # with in pixels
+img_galex_rep_zoom = img_galex_rep[int(center[1]-width[1]/2):int(center[1]+width[1]/2) , int(center[0]-width[0]/2):int(center[0]+width[0]/2)]
+img_map_zoom = img_map[int(center[1]-width[1]/2):int(center[1]+width[1]/2) , int(center[0]-width[0]/2):int(center[0]+width[0]/2)]
+
+# plot PAH map
+axins.imshow(img_galex_rep_zoom , origin="lower", vmin=lims[0], vmax=lims[1], norm="linear", cmap="Blues")
+
+# Overplot contours
+axins.contour(img_map_zoom, levels=[0.3,0.4,0.5], colors="red", linewidths=1)
 
 plt.show()
 ```
 
+The figure shows the PAH $3.3\,{\rm \mu m}$ emission (red contours) on top of the GALEX FUV continuum emission (blue). The inset shows a zoom in on the central regions.
+This comparison confirms that the PAH $3.3\,{\rm \mu m}$ emission generally follows the location of young hot stars. However, there are also differences (see zoom-in) where there are bright UV regions devoid of PAH emission. This may be the case where the strong ionization fields of young stars dissolve the small dust grains.
+
+
 ## Acknowledgements
 
 - [Caltech/IPAC-IRSA](https://irsa.ipac.caltech.edu/)
@@ -454,7 +522,7 @@ plt.show()
 
 **Contact:** Contact [IRSA Helpdesk](https://irsa.ipac.caltech.edu/docs/help_desk.html) with questions or problems.
 
-**Runtime:** Approximately 30 seconds.
+**Runtime:** This notebook takes about 20 seconds to run to completion on a machine with 32GB RAM and 8 CPU.
 
 ```python
 

From dbc3f89f5a66edbe0f4e1785d97e79f6a893eff6 Mon Sep 17 00:00:00 2001
From: Andreas Faisst <afaisst@ipac.caltech.edu>
Date: Wed, 13 May 2026 13:46:32 -0700
Subject: [PATCH 4/7] updated arguments in measure_aperture_photometry_cube()

---
 tutorials/spherex/spherex_mosaic/spherex_mosaic.md | 11 +++--------
 1 file changed, 3 insertions(+), 8 deletions(-)

diff --git a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
index 81970165..d19cddce 100644
--- a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
+++ b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
@@ -200,12 +200,7 @@ This function provides a very simple aperture photometry with background subtrac
 <!-- #endregion -->
 
 ```python
-def measure_aperture_photometry_cube(cube,
-                            r_aperture_px = 20,
-                            r_inner_px = 60,
-                            r_outer_px = 100,
-                            MAKEPLOT = True
-                           ):
+def measure_aperture_photometry_cube(cube, *, r_aperture_px=20, r_inner_px=60, r_outer_px=100, makeplot=True):
     '''
     Create quick look spectrum from simple aperture photometry with background subtraction
     from annulus.
@@ -281,7 +276,7 @@ def measure_aperture_photometry_cube(cube,
     phot_all.rename_column("id","plane")
 
     ## Plot if asked for
-    if MAKEPLOT:
+    if makeplot:
 
         fig = plt.figure(figsize=(5,5))
         ax1 = fig.add_subplot(1,1,1)
@@ -311,7 +306,7 @@ def measure_aperture_photometry_cube(cube,
 ```
 
 ```python
-phot_table = measure_aperture_photometry_cube(cube = cube_img , r_aperture_px = 20, r_inner_px = 60, r_outer_px= 100, MAKEPLOT = True)
+phot_table = measure_aperture_photometry_cube(cube = cube_img , r_aperture_px = 20, r_inner_px = 60, r_outer_px= 100, makeplot = True)
 phot_table
 ```
 

From bc92557acfbc96e12ee9780a943ce5150b78514f Mon Sep 17 00:00:00 2001
From: Andreas Faisst <afaisst@ipac.caltech.edu>
Date: Wed, 13 May 2026 14:37:26 -0700
Subject: [PATCH 5/7] proper jupytext conversion

---
 .../spherex/spherex_mosaic/spherex_mosaic.md  | 81 +++++++++----------
 1 file changed, 40 insertions(+), 41 deletions(-)

diff --git a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
index d19cddce..0f7dae77 100644
--- a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
+++ b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
@@ -1,21 +1,21 @@
 ---
-jupyter:
-  jupytext:
-    cell_metadata_filter: -all
-    notebook_metadata_filter: all,-language_info
-    text_representation:
-      extension: .md
-      format_name: markdown
-      format_version: '1.3'
-      jupytext_version: 1.17.2
-  kernelspec:
-    display_name: Python 3 (ipykernel)
-    language: python
-    name: python3
+jupytext:
+  cell_metadata_filter: -all
+  notebook_metadata_filter: all,-language_info
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.17.2
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
 ---
 
 # Understanding and Analyzing SPHEREx Mosaic Cubes 
 
++++
 
 ## 1. Learning Goals
 
@@ -24,6 +24,7 @@ jupyter:
 * Create a PAH $3.3\,{\rm \mu m}$ emission line map from a plane in the cube.
 * Obtain the corresponding GALEX FUV image from IRSA and overlay the PAH $3.3\,{\rm \mu m}$ map.
 
++++
 
 ## 2. SPHEREx Overview
 
@@ -39,16 +40,17 @@ The community will also mine SPHEREx data and combine it with synergistic data s
 
 More information is available in the [SPHEREx Explanatory Supplement](https://irsa.ipac.caltech.edu/data/SPHEREx/docs/SPHEREx_Expsupp_QR.pdf).
 
++++
 
 ## 3. Imports
 The following packages must be installed to run this notebook.
 
-```python
+```{code-cell} ipython3
 # Uncomment the next line to install dependencies if needed.
 # !pip install astropy matplotlib numpy photutils reproject astroquery
 ```
 
-```python
+```{code-cell} ipython3
 import copy
 
 import numpy as np
@@ -71,7 +73,7 @@ import matplotlib.pyplot as plt
 import matplotlib as mpl
 ```
 
-```python
+```{code-cell} ipython3
 ## Define some plotting formats
 def load_plotting_defaults():
     mpl.rcParams['xtick.minor.visible'] = True
@@ -106,7 +108,7 @@ To retrieve the same mosaic as provided here, go to the [IRSA mosaic tool](https
 
 We first define the path to the mosaic cube.
 
-```python
+```{code-cell} ipython3
 fn_spherex = "M101_irsa.fits"
 ```
 
@@ -123,10 +125,11 @@ The header includes multiple WCS axes including a wavelength axis. In order to e
 In addition, some metadata on `NAXIS` is preserved in the header which gives an error later when we run the `reproject` algorithm. We therefore remove the wavelength axis here by using `hdr.pop('NAXIS3', None)`.
 ```
 
++++
 
 We can now open the cube taking into account the changes and issues we discussed above:
 
-```python
+```{code-cell} ipython3
 with fits.open(fn_spherex) as hdul:
     hdul.info()
 
@@ -149,7 +152,6 @@ with fits.open(fn_spherex) as hdul:
     hdr2.pop('NAXIS3', None) # remove wavelength axis
     cube_wcs = WCS(hdr2 , fobj=hdul, naxis=2)
     print(f"Loaded SPHEREx cube with pixel scale {pixscale_mosaic} arcsec/px")
-
 ```
 
 Note that the FITS has different layers.
@@ -160,11 +162,13 @@ Note that the FITS has different layers.
 * The `SPECTRAL_CHANNELS` layer summarizes some useful additional information on the spectral channels used to create the cube planes in a binary table.
 * The `WCS-WAVE` layer contains a binary table summarizing the wavelength information for each plane in the cube.
 
++++
 
 ## 5. Create a Quick Look Spectrum
 
 Once we have the cube loaded, let's generate a quick-look spectrum. This has two purposes: First, it shows how to extract a spectrum from a cube using the `PhotUtils` python package. Second, it allows us to visualize the wavelength channels (i.e., planes) which we will need to chose the planes to make the final PAH $3.3\,{\rm \mu m}$ map.
 
++++
 
 ### 5.1 Create a Wavelength Look-Up Table
 
@@ -172,7 +176,7 @@ Later we will measure an aperture flux on each plane (representing a different w
 
 We construct the wavelenght look-up table from the `WCS-WAVE` layer, which summarizes the wavelength information.
 
-```python
+```{code-cell} ipython3
 wave_tab = Table( [wave_tab_tmp["PLANE"][0],
                   np.asarray([ww[0] for ww in wave_tab_tmp["WAVELENGTHS"][0] ]),
                   np.asarray([ww[0] for ww in wave_tab_tmp["BANDWIDTHS"][0] ])],
@@ -182,11 +186,10 @@ wave_tab = Table( [wave_tab_tmp["PLANE"][0],
 
 Let's have a look at this table. The `bandwidths` correspond to the width of each wavelength channel.
 
-```python
+```{code-cell} ipython3
 wave_tab
 ```
 
-<!-- #region -->
 ### 5.2 Extract Spectrum
 
 To obtain the spectrum, we create a function that computes the sum of the fluxes in apertures and does a background subtraction using an annulus. Both aperture and annulus sizes are user-defined. We use the `PhotUtils` package for this (see [here](https://photutils.readthedocs.io/en/latest/user_guide/aperture.html) for more information).
@@ -197,9 +200,8 @@ Furthermore, the function creates a plot showing the collapsed cube and the aper
 ```{tip}
 This function provides a very simple aperture photometry with background subtraction. The method can be extended. For example the function does not calculate photometric errors.  
 ```
-<!-- #endregion -->
 
-```python
+```{code-cell} ipython3
 def measure_aperture_photometry_cube(cube, *, r_aperture_px=20, r_inner_px=60, r_outer_px=100, makeplot=True):
     '''
     Create quick look spectrum from simple aperture photometry with background subtraction
@@ -305,7 +307,7 @@ def measure_aperture_photometry_cube(cube, *, r_aperture_px=20, r_inner_px=60, r
     return(phot_all)
 ```
 
-```python
+```{code-cell} ipython3
 phot_table = measure_aperture_photometry_cube(cube = cube_img , r_aperture_px = 20, r_inner_px = 60, r_outer_px= 100, makeplot = True)
 phot_table
 ```
@@ -316,13 +318,13 @@ Note that are some NaN values in that table. This is because of missing data for
 
 Finally, we also add the wavelength from the wavelength look-up table that created above to the photometry table. This makes it easy to plot the spectrum afterwards.
 
-```python
+```{code-cell} ipython3
 phot_table["wavelengths"] = wave_tab["wavelengths"].copy()
 ```
 
 Now that we have all the information, we can finally plot the spectrum. In addition, we also add some prominent emission lines (sorted by wavelength in the legend) and indicate the cube plane numbers around the $3.3\,{\rm \mu m}$ PAH feature. These numbers will be needed later to create the emission line map.
 
-```python
+```{code-cell} ipython3
 ## Plot Spectrum
 fig = plt.figure(figsize=(7,4))
 ax1 = fig.add_subplot(1,1,1)
@@ -360,7 +362,7 @@ Note that the plane IDs are 1-indexed (i.e. start from 1). However, the cube is
 
 For the cube creation, we set up a handy function:
 
-```python
+```{code-cell} ipython3
 def make_map(cube,
              planes_feature = [63],
              planes_continuum = [61,62,64,65]
@@ -389,13 +391,13 @@ def make_map(cube,
 
 And now we can easily extract the cube with the `planes_feature` and `planes_continuum` defined above.
 
-```python
+```{code-cell} ipython3
 img_map = make_map(cube_img , planes_feature = [63] , planes_continuum = [61,62,64,65])
 ```
 
 Finally, we can plot our PAH $3.3\,{\rm \mu m}$ SPHEREx map!
 
-```python
+```{code-cell} ipython3
 fig = plt.figure(figsize=(5,5))
 ax1 = fig.add_subplot(1,1,1)
 lims = np.nanpercentile(img_map.flatten() , q=(5,99.5))
@@ -408,6 +410,7 @@ plt.show()
 
 The map shows the location of the PAH $3.3\,{\rm \mu m}$ emission in the local galaxy M101. Specifically, the $3.3\,{\rm \mu m}$ feature is an observational indicator of very small carbonageous dust grains. You can see that the PAH emission is not continuous as it is tied to regions of star formation. To explore this further, we can correlate the PAH map with the GALEX far-UV (FUV) continuum map, which maps the UV emission of hot, young stars. This is done in the next section.
 
++++
 
 ## 7. Get Corresponding GALEX Image
 
@@ -415,16 +418,15 @@ Next we retrieve the GALEX image corresponding to the SPHEREx mosaic. This is sc
 
 We first obtain the position in coordinates at the center of the SPHEREx cube. Note that this is a 3d cube, so we have to take axes 1 and 2 (0-indexed) to obtain the coordinates!
 
-```python
+```{code-cell} ipython3
 radec = cube_wcs.all_pix2world(cube_img.shape[2]//2 , cube_img.shape[1]//2 , 0)
 print(f"Sky position at R.A. = {radec[0]} and Decl. = {radec[1]}")
 pos = SkyCoord(ra=radec[0], dec=radec[1], unit='deg')
-
 ```
 
 We then query IRSA at this position to obtain the GALEX FUV image. This data product is part of the [Spitzer Local Volume Legacy Survey](https://irsa.ipac.caltech.edu/data/SPITZER/LVL/overview.html), which observes many local galaxies. Accordingly, we query the collection `spitzer_lvl` and search for `science` images from GALEX in the `FUV` energy passband.
 
-```python
+```{code-cell} ipython3
 im_table = Irsa.query_sia(pos=(pos, 1 * u.arcsec), collection='spitzer_lvl')
 sel = np.where((im_table["dataproduct_subtype"] == "science") & (im_table["instrument_name"] == "GALEX") & (im_table["energy_bandpassname"] == "FUV") )[0]
 im_table[sel]
@@ -436,7 +438,7 @@ Once we have the table, we can access the URL of the image in the `access_url` c
 Note that we also update the header of the image once we created the cutout. This will be necessary when we reproject the GALEX image onto the SPHEREx mosaic pixel scale.
 ```
 
-```python
+```{code-cell} ipython3
 image_url = im_table['access_url'][sel][0]
 size = u.Quantity([40,40], u.arcmin)
 
@@ -453,7 +455,7 @@ with fits.open(image_url, memmap=False) as hdul:
 
 Additionally, we do some post-processing on the image to convert the pixel units to milli-Jansky, similar as the SPHEREx mosaic. Note that the zero points for the GALEX FUV and NUV images are different as defined above. We refer to the [Spitzer LVL documentation](https://irsa.ipac.caltech.edu/data/SPITZER/LVL/doc/LVL_DR5_v5.pdf) for more details on the GALEX pixel units.
 
-```python
+```{code-cell} ipython3
 zps = {"fuv":18.82, "nuv":20.08}
 img_galex_ab = -2.5*np.log10(img_galex) + zps["fuv"]
 img_galex = 10**(-0.4*(img_galex_ab - 23.9)) / 1e3 # AB -> ujy -> mjy
@@ -461,7 +463,7 @@ img_galex = 10**(-0.4*(img_galex_ab - 23.9)) / 1e3 # AB -> ujy -> mjy
 
 In a last step, we will have to reproject the GALEX image to the SPHEREx mosaic. This is necesary so we can overlap the GALEX image with the PAH $3.3\,{\rm \mu m}$ map that we have created above. For this, we use the `reproject` package.
 
-```python
+```{code-cell} ipython3
 img_galex_rep, fp = reproject_exact(input_data = (img_galex , hdr_galex), output_projection = cube_wcs )
 img_galex_rep[img_galex_rep == 0] = np.nan
 img_galex_rep = img_galex_rep / np.nansum(img_galex_rep) * np.nansum(img_galex)
@@ -469,7 +471,7 @@ img_galex_rep = img_galex_rep / np.nansum(img_galex_rep) * np.nansum(img_galex)
 
 Finally, we caon overlay the PAH emission on the GALEX map, which is our final result of this tutorial notebook.
 
-```python
+```{code-cell} ipython3
 fig = plt.figure(figsize=(9,9))
 ax1 = fig.add_subplot(1,1,1)
 
@@ -506,6 +508,7 @@ plt.show()
 The figure shows the PAH $3.3\,{\rm \mu m}$ emission (red contours) on top of the GALEX FUV continuum emission (blue). The inset shows a zoom in on the central regions.
 This comparison confirms that the PAH $3.3\,{\rm \mu m}$ emission generally follows the location of young hot stars. However, there are also differences (see zoom-in) where there are bright UV regions devoid of PAH emission. This may be the case where the strong ionization fields of young stars dissolve the small dust grains.
 
++++
 
 ## Acknowledgements
 
@@ -518,7 +521,3 @@ This comparison confirms that the PAH $3.3\,{\rm \mu m}$ emission generally foll
 **Contact:** Contact [IRSA Helpdesk](https://irsa.ipac.caltech.edu/docs/help_desk.html) with questions or problems.
 
 **Runtime:** This notebook takes about 20 seconds to run to completion on a machine with 32GB RAM and 8 CPU.
-
-```python
-
-```

From 2fbea4df32204f74b19a7ea40b7ccc8dba44454b Mon Sep 17 00:00:00 2001
From: Jessica Krick <jkrick@caltech.edu>
Date: Wed, 13 May 2026 18:24:43 -0400
Subject: [PATCH 6/7] simplify WCS extraction, improve spectrum plot labels,
 and select PAH planes by wavelength

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 .../spherex/spherex_mosaic/spherex_mosaic.md  | 86 +++++++++++--------
 1 file changed, 48 insertions(+), 38 deletions(-)

diff --git a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
index 0f7dae77..813b58c3 100644
--- a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
+++ b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
@@ -6,14 +6,14 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.17.2
+    jupytext_version: 1.19.1
 kernelspec:
-  display_name: Python 3 (ipykernel)
-  language: python
   name: python3
+  display_name: python3
+  language: python
 ---
 
-# Understanding and Analyzing SPHEREx Mosaic Cubes 
+# Understanding and Analyzing SPHEREx Mosaic Cubes
 
 +++
 
@@ -112,22 +112,21 @@ We first define the path to the mosaic cube.
 fn_spherex = "M101_irsa.fits"
 ```
 
-Next, we open the image. Note that we convert here the image from MJy/sr to mJy using the relation 
+Next, we open the image and convert the flux units from MJy/sr to mJy:
 
-```
-[mJy] = [MJy/sr] x 23.5045 x (pixelscale)^2 / 1E3
-```
+$$[\rm{mJy}] = [\rm{MJy/sr}] \times 23.5045 \times (\rm{pixel\ scale})^2 / 1000$$
 
-In addition, we also adjust the `BUNIT` keyword in the header to reflect the changes.
+We also update the `BUNIT` keyword in the header to reflect this change.
 
-```{warning}
-The header includes multiple WCS axes including a wavelength axis. In order to extract the spatial WCS (which is used, for example, for plotting or reprojection), we have to pass `astropy.wcs.WCS()` the header *and* the HDU list as well as specify the axes (the first two): `wcs = WCS(hdr, fobj=hdul, naxis=2)`. For more information, see the [astropy WCS documentation](https://docs.astropy.org/en/stable/api/astropy.wcs.WCS.html).
-In addition, some metadata on `NAXIS` is preserved in the header which gives an error later when we run the `reproject` algorithm. We therefore remove the wavelength axis here by using `hdr.pop('NAXIS3', None)`.
+The cube's WCS describes three axes — right ascension, declination, and wavelength. For spatial operations like plotting and reprojection we only need the two sky axes, which we extract with `WCS.celestial`.
+
+```{note}
+We pass `fobj=hdul` when constructing the WCS so that astropy has access to the full FITS file, not just the header. This is necessary when the WCS solution references auxiliary data — such as distortion corrections or lookup tables — stored in other extensions of the file.
 ```
 
 +++
 
-We can now open the cube taking into account the changes and issues we discussed above:
+We can now open the cube:
 
 ```{code-cell} ipython3
 with fits.open(fn_spherex) as hdul:
@@ -146,11 +145,8 @@ with fits.open(fn_spherex) as hdul:
     cube_img = hdul["IMAGE"].data * 23.5045 * (pixscale_mosaic)**2 / 1e3 # Mjy/sr -> mjy
     cube_hdr['BUNIT'] = "mjy"
     
-    # get spatial WCS
-    hdr2 = cube_hdr.copy()
-    hdr2["NAXIS"] = 2
-    hdr2.pop('NAXIS3', None) # remove wavelength axis
-    cube_wcs = WCS(hdr2 , fobj=hdul, naxis=2)
+    # extract the 2D spatial WCS (dropping the spectral axis)
+    cube_wcs = WCS(cube_hdr, fobj=hdul).celestial
     print(f"Loaded SPHEREx cube with pixel scale {pixscale_mosaic} arcsec/px")
 ```
 
@@ -322,7 +318,7 @@ Finally, we also add the wavelength from the wavelength look-up table that creat
 phot_table["wavelengths"] = wave_tab["wavelengths"].copy()
 ```
 
-Now that we have all the information, we can finally plot the spectrum. In addition, we also add some prominent emission lines (sorted by wavelength in the legend) and indicate the cube plane numbers around the $3.3\,{\rm \mu m}$ PAH feature. These numbers will be needed later to create the emission line map.
+Now that we have all the information, we can plot the spectrum. We also mark some prominent near-infrared emission features for reference.
 
 ```{code-cell} ipython3
 ## Plot Spectrum
@@ -332,19 +328,20 @@ ax1 = fig.add_subplot(1,1,1)
 # spectrum
 ax1.plot(phot_table["wavelengths"] , phot_table["aperture_sum_bkgsub"], "o", markersize=3, markerfacecolor="black", markeredgecolor="black")
 
-# some emission lines
-ax1.axvline(1.28, linestyle=":", color="gray", linewidth=0.5, label=r"Pa-$\beta$")
-ax1.axvline(1.6, linestyle=":", color="gray", linewidth=0.5, label=r"$1.6\,{\rm \mu m}$ bump")
-ax1.axvline(1.87, linestyle=":", color="gray", linewidth=0.5, label=r"Pa-$\alpha$")
-ax1.axvline(2.36, linestyle=":", color="gray", linewidth=0.5, label=r"Br-$\beta$")
-ax1.axvline(3.3, linestyle=":", color="gray", linewidth=0.5, label=r"PAH $3.3\,{\rm \mu m}$")
-ax1.axvline(4.05, linestyle=":", color="gray", linewidth=0.5, label=r"Br-$\alpha$")
-
-# add plane numbers around PAH 3.3
-sel = np.where( (phot_table["wavelengths"] > 2.6) & (phot_table["wavelengths"] < 3.8) & (~np.isnan(phot_table["aperture_sum_bkgsub"])) )[0]
-[ax1.text(phot_table["wavelengths"][ss] , phot_table["aperture_sum_bkgsub"][ss]*1.1, phot_table["plane"][ss] , fontsize=7, va="bottom", ha="center", rotation=90) for ss in sel]
-
-ax1.legend(fontsize=8)
+# label emission lines directly on the plot
+lines = [
+    (1.28, r"Pa-$\beta$"),
+    (1.60, r"1.6 $\mu$m bump"),
+    (1.87, r"Pa-$\alpha$"),
+    (2.36, r"Br-$\beta$"),
+    (3.30, r"PAH 3.3 $\mu$m"),
+    (4.05, r"Br-$\alpha$"),
+]
+trans = mpl.transforms.blended_transform_factory(ax1.transData, ax1.transAxes)
+for wave, label in lines:
+    ax1.axvline(wave, linestyle=":", color="gray", linewidth=0.5)
+    ax1.text(wave, 0.97, label, transform=trans, fontsize=7,
+             va="top", ha="center", rotation=90, color="gray")
 ax1.tick_params(which="both", axis="both", labelsize=12)
 ax1.set_xlabel(r"Wavelength [$\mu$m]", fontsize=12)
 ax1.set_ylabel(r"Flux [mJy]", fontsize=12)
@@ -354,13 +351,26 @@ plt.show()
 
 ## 6. Create a PAH $3.3\,{\rm \mu m}$ Map
 
-We now create the PAH $3.3\,{\rm \mu m}$ map from the cube. Using the spectrum plot above, we can identify the cube wavelength planes that include the PAH feature (plane 63) and the planes that include the continuum (planes 61/62 on the blue and 64/65 on the red). The emission line map is then simply created by subtracting the continuum planes from the plane containing the emission line.
+We now create the PAH $3.3\,{\rm \mu m}$ map from the cube. The emission line map is created by subtracting the median of nearby continuum planes from the feature plane.
+
+We first identify the relevant planes by matching the observed PAH wavelength to the entries in `wave_tab`. For M101 the redshift is negligible, but the same code works for any target — just set `z` to the known redshift.
+
+```{code-cell} ipython3
+z = 0.0  # redshift of target
+pah_obs = 3.3 * (1 + z)  # observed PAH wavelength in microns
+
+pah_idx = np.argmin(np.abs(wave_tab["wavelengths"] - pah_obs))
+planes_feature = [wave_tab["plane"][pah_idx]]
+planes_continuum = list(wave_tab["plane"][[pah_idx - 2, pah_idx - 1, pah_idx + 1, pah_idx + 2]])
+print(f"Feature plane: {planes_feature}  ({wave_tab['wavelengths'][pah_idx]:.3f} μm)")
+print(f"Continuum planes: {planes_continuum}")
+```
 
-```{warning}
-Note that the plane IDs are 1-indexed (i.e. start from 1). However, the cube is an *ndarray*, which is 0-indexed. If we want to access wavelength plane 63, we will have to chose `cube[63-1,:,:]`.
+```{note}
+For a target at higher redshift, replace `z = 0.0` with the known redshift. The observed PAH wavelength shifts to `3.3 * (1 + z)` microns, and the correct planes are selected automatically from `wave_tab`.
 ```
 
-For the cube creation, we set up a handy function:
+For the map creation, we set up a handy function:
 
 ```{code-cell} ipython3
 def make_map(cube,
@@ -389,10 +399,10 @@ def make_map(cube,
     return(img_map)
 ```
 
-And now we can easily extract the cube with the `planes_feature` and `planes_continuum` defined above.
+With the planes identified, we can now build the map:
 
 ```{code-cell} ipython3
-img_map = make_map(cube_img , planes_feature = [63] , planes_continuum = [61,62,64,65])
+img_map = make_map(cube_img, planes_feature=planes_feature, planes_continuum=planes_continuum)
 ```
 
 Finally, we can plot our PAH $3.3\,{\rm \mu m}$ SPHEREx map!

From 393345250d9f4dfeb08775c3203d7a8422501275 Mon Sep 17 00:00:00 2001
From: Andreas Faisst <afaisst@ipac.caltech.edu>
Date: Wed, 13 May 2026 16:41:11 -0700
Subject: [PATCH 7/7] small changes in spectrum plot

---
 .../spherex/spherex_mosaic/spherex_mosaic.md  | 47 +++++++++++++++++--
 1 file changed, 42 insertions(+), 5 deletions(-)

diff --git a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
index 813b58c3..f1253d4a 100644
--- a/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
+++ b/tutorials/spherex/spherex_mosaic/spherex_mosaic.md
@@ -6,10 +6,10 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.19.1
+    jupytext_version: 1.17.2
 kernelspec:
   name: python3
-  display_name: python3
+  display_name: Python 3 (ipykernel)
   language: python
 ---
 
@@ -337,11 +337,14 @@ lines = [
     (3.30, r"PAH 3.3 $\mu$m"),
     (4.05, r"Br-$\alpha$"),
 ]
+
+# plot spectrum and emission lines
 trans = mpl.transforms.blended_transform_factory(ax1.transData, ax1.transAxes)
 for wave, label in lines:
-    ax1.axvline(wave, linestyle=":", color="gray", linewidth=0.5)
-    ax1.text(wave, 0.97, label, transform=trans, fontsize=7,
-             va="top", ha="center", rotation=90, color="gray")
+    ax1.axvline(wave, ymin=0.3, linestyle=":", color="gray", linewidth=0.5)
+    ax1.text(wave, 0.05, label, transform=trans, fontsize=7,
+             va="bottom", ha="center", rotation=90, color="gray")
+
 ax1.tick_params(which="both", axis="both", labelsize=12)
 ax1.set_xlabel(r"Wavelength [$\mu$m]", fontsize=12)
 ax1.set_ylabel(r"Flux [mJy]", fontsize=12)
@@ -370,6 +373,40 @@ print(f"Continuum planes: {planes_continuum}")
 For a target at higher redshift, replace `z = 0.0` with the known redshift. The observed PAH wavelength shifts to `3.3 * (1 + z)` microns, and the correct planes are selected automatically from `wave_tab`.
 ```
 
+We can also visualize the chosen planes on the spectrum plot to for additional checking.
+
+```{code-cell} ipython3
+## Plot Spectrum (again, for checking planes)
+fig = plt.figure(figsize=(7,4))
+ax1 = fig.add_subplot(1,1,1)
+
+# spectrum
+ax1.plot(phot_table["wavelengths"] , phot_table["aperture_sum_bkgsub"], "o", markersize=3, markerfacecolor="black", markeredgecolor="black")
+
+# plot spectrum and emission lines
+trans = mpl.transforms.blended_transform_factory(ax1.transData, ax1.transAxes)
+for wave, label in lines:
+    ax1.axvline(wave, linestyle=":", color="gray", linewidth=0.5)
+    ax1.text(wave, 0.05, label, transform=trans, fontsize=7,
+             va="bottom", ha="center", rotation=90, color="gray")
+
+# indicated planes
+[ax1.text(phot_table["wavelengths"][ss-1] , phot_table["aperture_sum_bkgsub"][ss-1]*1.1, phot_table["plane"][ss-1] , fontsize=7, va="bottom", ha="center", rotation=90, color="blue") for ss in planes_feature]
+[ax1.text(phot_table["wavelengths"][ss-1] , phot_table["aperture_sum_bkgsub"][ss-1]*1.1, phot_table["plane"][ss-1] , fontsize=7, va="bottom", ha="center", rotation=90, color="red") for ss in planes_continuum]
+ax1.plot(phot_table["wavelengths"][np.asarray(planes_feature)-1] , phot_table["aperture_sum_bkgsub"][np.asarray(planes_feature)-1], "o", markersize=3, markerfacecolor="blue", markeredgecolor="blue")
+ax1.plot(phot_table["wavelengths"][np.asarray(planes_continuum)-1] , phot_table["aperture_sum_bkgsub"][np.asarray(planes_continuum)-1], "o", markersize=3, markerfacecolor="red", markeredgecolor="red")
+
+ax1.tick_params(which="both", axis="both", labelsize=12)
+ax1.set_xlabel(r"Wavelength [$\mu$m]", fontsize=12)
+ax1.set_ylabel(r"Flux [mJy]", fontsize=12)
+
+plt.show()
+```
+
+```{note}
+If you create maps for other lines, it might be worth checking the plane numbers on the spectrum directly to avoid contamination of other line when defining the continuum planes.
+```
+
 For the map creation, we set up a handy function:
 
 ```{code-cell} ipython3