jkrick · Oct 28, 2026 · Oct 22, 2026 · Oct 22, 2026 · Oct 22, 2026 · Oct 22, 2026
diff --git a/deploy_to_fornax_manifest.in b/deploy_to_fornax_manifest.in
@@ -20,7 +20,6 @@ tutorials/spherex/spherex_cutouts.md
 tutorials/spherex/spherex_intro.md
 tutorials/spherex/spherex_psf.md
 tutorials/spherex/spherex_source_discovery/conda-spherex_sdt.yml
-tutorials/spherex/spherex_source_discovery/README.md
 tutorials/spherex/spherex_source_discovery/spherex_source_discovery_tool_demo.md
 tutorials/spherex/spherex_source_discovery/spherex_source_discovery_tool/__init__.py
 tutorials/spherex/spherex_source_discovery/spherex_source_discovery_tool/aperture_photometry.py

diff --git a/tutorials/spherex/spherex_cutouts.md b/tutorials/spherex/spherex_cutouts.md
@@ -46,7 +46,7 @@ The following packages must be installed to run this notebook.
 
 ```{code-cell} ipython3
 # Uncomment the next line to install dependencies if needed.
-# !pip install numpy astropy pyvo matplotlib
+# !pip install numpy astropy astroquery matplotlib
 ```
 
 ```{code-cell} ipython3
@@ -58,11 +58,12 @@ import urllib.error
 import astropy.units as u
 import matplotlib.pyplot as plt
 import numpy as np
-import pyvo
 from astropy.coordinates import SkyCoord
 from astropy.io import fits
 from astropy.table import Table
+from astropy.visualization import quantity_support
 from astropy.wcs import WCS
+from astroquery.ipac.irsa import Irsa
 
 # Suppress logging temporarily to prevent astropy
 # from repeatedly printing out warning notices related to alternate WCSs
@@ -73,6 +74,9 @@ logging.getLogger('astropy').setLevel(logging.ERROR)
 # astropy's default timeout limit. Increase it.
 from astropy.utils.data import conf
 conf.remote_timeout = 120
+
+# Let matplotlib recognize astropy Quantities and auto-label axes with their units.
+quantity_support()
 ```
 
 ## 4. Specify inputs and outputs
@@ -83,8 +87,7 @@ In this example, we are creating cutouts of the Pinwheel galaxy (M101) for the S
 
 ```{code-cell} ipython3
 # Choose a position.
-ra = 210.80227 * u.degree
-dec = 54.34895 * u.degree
+coord = SkyCoord(ra=210.80227 * u.degree, dec=54.34895 * u.degree, frame="icrs")
 
 # Choose a cutout size.
 size = 0.1 * u.degree
@@ -98,34 +101,30 @@ output_filename = 'spherex_cutouts_mef.fits'
 
 ## 5. Query IRSA for a list of cutouts that satisfy the criteria specified above.
 
-Here we show how to use the `pyvo` TAP SQL query to retrieve all images that overlap with the position defined above.
+Here we show how to use `astroquery`'s IRSA TAP interface to retrieve all images that overlap with the position defined above.
 This query will retrieve a table of URLs that link to the MEF cutouts.
 Each row in the table corresponds to a single cutout and includes the data access URL and an observation timestamp.
 The results are sorted from oldest to newest.
 
 ```{code-cell} ipython3
-# Define the service endpoint for IRSA's Table Access Protocol (TAP)
-# so that we can query SPHEREx metadata tables.
-service = pyvo.dal.TAPService("https://irsa.ipac.caltech.edu/TAP")
-
 # Define a query that will search the appropriate SPHEREx metadata tables
 # for spectral images that cover the chosen coordinates and match the
 # specified bandpass. Return the cutout data access URL and the time of observation.
 # Sort by observation time.
 query = f"""
 SELECT
-    'https://irsa.ipac.caltech.edu/' || a.uri || '?center={ra.value},{dec.value}d&size={size.value}' AS uri,
+    'https://irsa.ipac.caltech.edu/' || a.uri || '?center={coord.ra.deg},{coord.dec.deg}d&size={size.value}' AS uri,
     p.time_bounds_lower
 FROM spherex.artifact a
 JOIN spherex.plane p ON a.planeid = p.planeid
-WHERE 1 = CONTAINS(POINT('ICRS', {ra.value}, {dec.value}), p.poly)
+WHERE 1 = CONTAINS(POINT('ICRS', {coord.ra.deg}, {coord.dec.deg}), p.poly)
         AND p.energy_bandpassname = '{bandpass}'
 ORDER BY p.time_bounds_lower
 """
 
 # Execute the query and return as an astropy Table.
 t1 = time.time()
-results = service.search(query)
+results = Irsa.query_tap(query)
 print("Time to do TAP query: {:2.2f} seconds.".format(time.time() - t1))
 print("Number of images found: {}".format(len(results)))
 ```
@@ -141,12 +140,13 @@ This function takes in a row of the catalog that we created above and does the f
 
 - It downloads the cutout
 - It computes the wavelength of the center pixel of the cutout (in micro-meters)
+- It records the cutout's processing version (the `VERSION` keyword in its primary header), so that it is not lost when many cutouts are later combined into one MEF
 - It combines the image HDUs into a new HDU and adds it to the table row.
 
 Note that the values of the rows are being added in place.
 
 ```{code-cell} ipython3
-def process_cutout(row, ra, dec, cache):
+def process_cutout(row, coord, cache):
     '''
     Downloads the cutouts given in a row of the table including all SPHEREx images overlapping with a position.
 
@@ -156,8 +156,8 @@ def process_cutout(row, ra, dec, cache):
     row : astropy.table row
         Row of a table that will be changed in place by this function. The table
         is created by the SQL TAP query.
-    ra,dec : coordinates (astropy units)
-        Ra and Dec coordinates (same as used for the TAP query) with attached astropy units
+    coord : astropy.coordinates.SkyCoord
+        Coordinates (same as used for the TAP query).
     cache : bool
         If set to `True`, the output of cached and the cutout processing will run faster next time.
         Turn this feature off by setting ``cache=False``.
@@ -171,14 +171,19 @@ def process_cutout(row, ra, dec, cache):
 
         # Compute pixel coordinates corresponding to cutout position.
         spatial_wcs = WCS(header)
-        x, y = spatial_wcs.world_to_pixel(SkyCoord(ra=ra, dec=dec, unit="deg", frame="icrs"))
+        x, y = spatial_wcs.world_to_pixel(coord)
 
         # Compute wavelength at cutout position.
         spectral_wcs = WCS(header, fobj=hdulist, key="W")
         spectral_wcs.sip = None
         wavelength, bandwidth = spectral_wcs.pixel_to_world(x, y)
         row["central_wavelength"] = wavelength.to(u.micrometer).value
 
+        # The FITS primary header carries the processing version (VERSION), which isn't repeated
+        # in the extension headers. Keep it so it isn't silently dropped when we build the
+        # combined MEF below.
+        row["version"] = hdulist[0].header.get("VERSION", "")
+
         # Collect the HDUs for this cutout and append the row's cutout_index to the EXTNAME.
         hdus = []
         for hdu in hdulist[1:]:  # skip the primary header
@@ -192,7 +197,7 @@ It catches transient read errors and simply skips those cutouts, which is suffic
 For science use cases, users may instead want to implement their own retry logic or error handling strategy.
 
 ```{code-cell} ipython3
-def process_cutout_with_error_handling(row, ra, dec, cache):
+def process_cutout_with_error_handling(row, coord, cache):
     '''
     Call `process_cutout` while catching transient read errors.
 
@@ -202,15 +207,16 @@ def process_cutout_with_error_handling(row, ra, dec, cache):
     row : astropy.table row
         Row of a table that will be changed in place by this function. The table
         is created by the SQL TAP query.
-    ra,dec : coordinates (astropy units)
-        Ra and Dec coordinates (same as used for the TAP query) with attached astropy units
+    coord : astropy.coordinates.SkyCoord
+        Coordinates (same as used for the TAP query).
     cache : bool
         If set to `True`, the output of cached and the cutout processing will run faster next time.
         Turn this feature off by setting ``cache=False``.
     '''
     try:
-        process_cutout(row, ra, dec, cache=cache)
-    except (TimeoutError, urllib.error.HTTPError, http.client.IncompleteRead):
+        process_cutout(row, coord, cache=cache)
+    except (TimeoutError, urllib.error.HTTPError, http.client.IncompleteRead,
+            http.client.RemoteDisconnected):
         # Transient read errors. Skip this cutout.
         row["hdus"] = None
 ```
@@ -247,11 +253,12 @@ To learn more about the cache please read the [astropy cache management document
 results_table_serial = results.to_table()
 results_table_serial["cutout_index"] = range(1, len(results_table_serial) + 1)
 results_table_serial["central_wavelength"] = np.full(len(results_table_serial), np.nan)
+results_table_serial["version"] = np.full(len(results_table_serial), "", dtype="U32")
 results_table_serial["hdus"] = np.full(len(results_table_serial), None)
 
 t1 = time.time()
 for row in results_table_serial:
-    process_cutout_with_error_handling(row, ra, dec, cache=False)
+    process_cutout_with_error_handling(row, coord, cache=False)
 print("Time to create cutouts in serial mode: {:2.2f} minutes.".format((time.time() - t1) / 60))
 
 # Drop rows that failed to download.
@@ -268,16 +275,7 @@ The choice of this value depends on the number of cores but also on the number o
 
 ```{tip}
 A good value for the maximum number of workers is between 7 and 12 for a machine with 8 cores.
-```
-
-```{tip}
-The astropy `fits.open()` supports a caching argument.
-This can be passed through in the `process_cutout()` function.
-If `cache=True` is set, the images are cached and the cutout creation is sped up next time the code is run (even if the Jupyter kernel is restarted!).
-The downside is that the images are saved on the machine where this notebook is run (usually in `~/.astropy/cache/`).
-If many cutouts are created, this can sum up to a large cached data volume, in which case `cache=False` is preferred.
-
-To learn more about the cache please read the [astropy cache management documentation](https://docs.astropy.org/en/stable/utils/data.html#cache-management).
+See the caching tip in section 7.1 above -- the same `cache` argument applies here.
 ```
 
 Again, before running the cutout processing we define some place holders.
@@ -286,11 +284,12 @@ Again, before running the cutout processing we define some place holders.
 results_table_parallel = results.to_table()
 results_table_parallel["cutout_index"] = range(1, len(results_table_parallel) + 1)
 results_table_parallel["central_wavelength"] = np.full(len(results_table_parallel), np.nan)
+results_table_parallel["version"] = np.full(len(results_table_parallel), "", dtype="U32")
 results_table_parallel["hdus"] = np.full(len(results_table_parallel), None)
 
 t1 = time.time()
 with concurrent.futures.ThreadPoolExecutor(max_workers=10) as executor:
-    futures = [executor.submit(process_cutout_with_error_handling, row, ra, dec, False)
+    futures = [executor.submit(process_cutout_with_error_handling, row, coord, False)
                for row in results_table_parallel]
     concurrent.futures.wait(futures)
 print("Time to create cutouts in parallel mode: {:2.2f} minutes.".format((time.time() - t1) / 60))
@@ -315,6 +314,7 @@ cols = fits.ColDefs([
     fits.Column(name="cutout_index", format="J", array=results_table_parallel["cutout_index"], unit=""),
     fits.Column(name="observation_date", format="D", array=results_table_parallel["time_bounds_lower"], unit="d"),
     fits.Column(name="central_wavelength", format="D", array=results_table_parallel["central_wavelength"], unit="um"),
+    fits.Column(name="version", format="A32", array=results_table_parallel["version"], unit=""),
     fits.Column(name="access_url", format="A200", array=results_table_parallel["uri"], unit=""),
 ])
 table_hdu = fits.BinTableHDU.from_columns(cols)
@@ -324,6 +324,7 @@ table_hdu.header["EXTNAME"] = "CUTOUT_INFO"
 ## 9. Create the final MEF
 
 Now, we create a primary HDU and combine it with the summary table HDU and the cutout image and wavelength HDUs to a final HDU list.
+The primary header is left blank because the cutouts generally come from different processing versions, so no single `VERSION` would be correct for all of them; the `version` column in `CUTOUT_INFO` above records each cutout's version.
 
 ```{code-cell} ipython3
 primary_hdu = fits.PrimaryHDU()
@@ -335,10 +336,11 @@ combined_hdulist = fits.HDUList(hdulist_list)
 ## 10. Write the final MEF
 
 Finally, we save the full HDU list to a multi-extension FITS file to disk.
+We use `overwrite=False` so that re-running this cell does not silently clobber an existing output file; delete the file first (or change `output_filename` above) if you want to regenerate it.
 
 ```{code-cell} ipython3
 # Write the final MEF
-combined_hdulist.writeto(output_filename, overwrite=True)
+combined_hdulist.writeto(output_filename, overwrite=False)
 ```
 
 ## 11. Test and visualize the final result
@@ -363,14 +365,18 @@ with fits.open(output_filename) as hdul:
 ```
 
 Plot the images with the wavelength corresponding to their central pixel.
+We read `central_wavelength` as an astropy `Quantity` (via the column's `.quantity` property, which picks up the `um` unit written into the FITS header) rather than a bare float, so the unit travels with the value instead of being assumed.
 
 ```{code-cell} ipython3
+central_wavelength = summary_table["central_wavelength"].quantity
+
 fig = plt.figure(figsize=(15, 5))
 axs = [fig.add_subplot(2, 5, ii + 1) for ii in range(10)]
 
 for ii, img in enumerate(imgs):
     axs[ii].imshow(imgs[ii], norm="log", origin="lower")
-    axs[ii].text(0.05, 0.05, r"$\lambda_{\rm center} = %2.4g \,{\rm \mu m}$" % summary_table["central_wavelength"][ii],
+    wl = central_wavelength[ii].to(u.micrometer)
+    axs[ii].text(0.05, 0.05, rf"$\lambda_{{\rm center}} = {wl.value:2.4g} \,{{\rm \mu m}}$",
                  va="bottom", ha="left", color="white", transform=axs[ii].transAxes)
 
 plt.show()

diff --git a/tutorials/spherex/spherex_intro.md b/tutorials/spherex/spherex_intro.md
@@ -164,8 +164,11 @@ Transient read errors occur sometimes, so we'll catch those and retry a few time
 max_retries = 3
 for attempt in range(max_retries):
     try:
-        # Load the data.
-        hdulist = fits.open(spectral_image_url)
+        # Load the data, reading every extension into memory so that
+        # hdulist stays usable after the file is closed.
+        with fits.open(spectral_image_url, memmap=False) as hdulist:
+            for hdu in hdulist:
+                hdu.data
-        # Load the data, reading every extension into memory so that
-        # hdulist stays usable after the file is closed.
-        with fits.open(spectral_image_url, memmap=False) as hdulist:
-            for hdu in hdulist:
-                hdu.data
+        # Load the data, reading every extension into memory so that
+        # hdulist stays usable after the file is closed.
+        with fits.open(spectral_image_url, memmap=False) as hdulist:
+            for hdu in hdulist:
+                hdu_data = hdu.data.copy()
-        # Load the data, reading every extension into memory so that
-        # hdulist stays usable after the file is closed.
-        with fits.open(spectral_image_url, memmap=False) as hdulist:
-            for hdu in hdulist:
-                hdu.data
+        # Load the data, reading every extension into memory so that
+        # hdulist stays usable after the file is closed.
+        with fits.open(spectral_image_url, memmap=False) as hdulist:
+            for hdu in hdulist:
+                hdu_data = hdu.data.copy()
         break
     except (TimeoutError, urllib.error.HTTPError, http.client.IncompleteRead):
         if attempt == max_retries - 1:

diff --git a/tutorials/spherex/spherex_psf.md b/tutorials/spherex/spherex_psf.md
@@ -58,7 +58,7 @@ The following packages must be installed to run this notebook.
 
 ```{code-cell} ipython3
 # Uncomment the next line to install dependencies if needed.
-# !pip install astropy matplotlib numpy pyvo photutils
+# !pip install astropy matplotlib numpy astroquery photutils
 ```
 
 ```{code-cell} ipython3
@@ -73,18 +73,22 @@ from packaging.version import Version
 import astropy.units as u
 import matplotlib.pyplot as plt
 import numpy as np
-import pyvo
 from astropy.coordinates import SkyCoord
 from astropy.io import fits
 from astropy.table import Table
+from astropy.visualization import quantity_support
 from astropy.wcs import WCS
+from astroquery.ipac.irsa import Irsa
 
 from photutils.psf import ImagePSF
 
 # The time it takes to read SPHEREx files can exceed
 # astropy's default timeout limit. Increase it.
 from astropy.utils.data import conf
 conf.remote_timeout = 120
+
+# Let matplotlib recognize astropy Quantities and auto-label axes with their units.
+quantity_support()
 ```
 
 ## 4. Get SPHEREx Cutout and Obtain Position on Parent Frame
@@ -104,35 +108,30 @@ To learn more about how to access SPHEREx spectral images and how to download cu
 ```
 
 ```{code-cell} ipython3
-ra = 305.59875000000005 * u.degree
-dec = 41.14888888888889 * u.degree
+coord = SkyCoord(ra=305.59875000000005 * u.degree, dec=41.14888888888889 * u.degree, frame="icrs")
 size = 0.01 * u.degree
 ```
 
 Once we defined the coordinates of interest and the size of the cutout, we run a TAP query to gather all SPHEREx spectral images that cover the coordinates.
 
 ```{code-cell} ipython3
-# Define the service endpoint for IRSA's Table Access Protocol (TAP)
-# so that we can query SPHEREx metadata tables.
-service = pyvo.dal.TAPService("https://irsa.ipac.caltech.edu/TAP")
-
 # Define a query that will search the appropriate SPHEREx metadata tables
 # for spectral images that cover the chosen coordinates of interest.
 # Return the cutout data access URL and the time of observation.
 # Sort by observation time.
 query = f"""
 SELECT
-    'https://irsa.ipac.caltech.edu/' || a.uri || '?center={ra.value},{dec.value}d&size={size.value}' AS uri,
+    'https://irsa.ipac.caltech.edu/' || a.uri || '?center={coord.ra.deg},{coord.dec.deg}d&size={size.value}' AS uri,
     p.time_bounds_lower, p.energy_bandpassname, p.time_bounds_lower, p.time_bounds_upper, p.provenance_version
 FROM spherex.artifact a
 JOIN spherex.plane p ON a.planeid = p.planeid
-WHERE 1 = CONTAINS(POINT('ICRS', {ra.value}, {dec.value}), p.poly)
+WHERE 1 = CONTAINS(POINT('ICRS', {coord.ra.deg}, {coord.dec.deg}), p.poly)
 ORDER BY p.time_bounds_lower
 """
 
 # Execute the query and return as an astropy Table.
 t1 = time.time()
-results = service.search(query)
+results = Irsa.query_tap(query)
 results = results.to_table()
 print("Time to do TAP query: {:2.2f} seconds.".format(time.time() - t1))
 print("Number of images found: {}".format(len(results)))
@@ -170,7 +169,8 @@ for attempt in range(max_retries):
             cutout_header = hdul['IMAGE'].header
             cutout = hdul['IMAGE'].data
         break
-    except (TimeoutError, urllib.error.HTTPError, http.client.IncompleteRead):
+    except (TimeoutError, urllib.error.HTTPError, http.client.IncompleteRead,
+            http.client.RemoteDisconnected):
         if attempt == max_retries - 1:
             raise
         time.sleep(10 * (attempt + 1))
@@ -220,7 +220,7 @@ We do this by first determining the pixel $(x,y)$ coordinates of our coordinates
 
 ```{code-cell} ipython3
 wcs = WCS(cutout_header)
-xpix_cutout, ypix_cutout = wcs.world_to_pixel(SkyCoord(ra=ra, dec=dec))
+xpix_cutout, ypix_cutout = wcs.world_to_pixel(coord)
 
 print(f"Pixel values of coordinates of interest on cutout image: x = {xpix_cutout}, y = {ypix_cutout}")
 ```
@@ -476,11 +476,13 @@ ax2.set_ylabel("$y$ PSF zones")
 plt.colorbar(im2 , ax=ax2, label=r"$N_{\rm eff}$ (PSF effective area)")
 
 ## As a function of wavelength
-x = epsf_bintable["CWAVE"] # wavelength at the center of the PSF zone
+# CWAVE has no unit metadata in the FITS table, so attach the known unit (microns)
+# explicitly rather than relying on a hardcoded axis label string.
+wave = epsf_bintable["CWAVE"] * u.micrometer  # wavelength at the center of the PSF zone
 y = epsf_bintable["NEFF_MEAN"]
-ax3.scatter(x , y, c=z , marker="o", s=5)
+ax3.scatter(wave, y, c=z , marker="o", s=5)
 
-ax3.set_xlabel(r"Wavelength [$\mu m$]")
+ax3.set_xlabel(f"Wavelength [{wave.unit.to_string('latex_inline')}]")
 ax3.set_ylabel(r"$N_{\rm eff}$ (PSF effective area)")
 
 plt.show()