diff --git a/deploy_to_fornax_manifest.in b/deploy_to_fornax_manifest.in index 57db50d2..3fbcbc4f 100644 --- 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 index ce50448b..44c37db9 100644 --- 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,7 +171,7 @@ 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") @@ -179,6 +179,11 @@ def process_cutout(row, ra, dec, cache): 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 index 48f6e4e8..1bad9159 100644 --- 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 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 index c96188f1..f629d475 100644 --- 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,11 +73,12 @@ 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 @@ -85,6 +86,9 @@ from photutils.psf import ImagePSF # 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() diff --git a/tutorials/spherex/spherex_source_discovery/README.md b/tutorials/spherex/spherex_source_discovery/README.md deleted file mode 100644 index e8ec7a09..00000000 --- a/tutorials/spherex/spherex_source_discovery/README.md +++ /dev/null @@ -1,62 +0,0 @@ -# SPHEREx Source Discovery Tool - -## Overview - -The SPHEREx Source Discovery Tool is the Python package `spherex_source_discovery_tool` (included in this directory), which is used to discover and extract sources from SPHEREx Spectral Images and visualize their spectra. -The notebook [](spherex_source_discovery_tool_demo.md) demonstrates how to use it. - -### Directory Contents - -- `spherex_source_discovery_tool_demo.md`: Demo notebook. -- `conda-spherex_sdt.yml`: Conda environment definition file. -- `spherex_source_discovery_tool/`: Python package, including functions and configuration files to extract and visualize sources. - - Python modules - - `aperture_photometry.py`: Functions adapted from Zafar Rustamkulov to perform aperture photometry. - - `bokeh_viz.py`: Functions for dynamic `bokeh` visualizations. - - `firefly_viz.py`: Functions for `Firefly` visualizations. - - `source_extraction.py`: Functions to work with SExtractor. - - `sdt_utils.py`: Miscellaneous utility functions for working with IRSA tables and SPHEREx image headers. - - SExtractor files - - `default.conv`: Convolution mask used as detection filter. - - `default.nnw`: Table of neural-network weights for star/galaxy separation. - - `default_sdt.sex`: SExtractor default config file. Specifies catalog, extraction, photometry, etc. options. - - `default_sdt.param`: SExtractor default parameters file. Specifies columns to include in output SExtractor catalog. - -## Setup - -### On local machine - -All necessary packages and tools (SExtractor) are listed in the `conda-spherex_sdt.yml` file. - -To create a conda environment with the dependencies on your local machine use: - ``` - conda env create --file conda-spherex_sdt.yml - ``` - -Then, make the environment available in the list of kernels for your JupyterLab: - ``` - conda activate spherex_sdt - python -m ipykernel install --user --name=spherex_sdt - ``` - -To use the environment in your Jupyter notebooks, either start JupyterLab in that environment by typing - -``` -conda activate spherex_sdt -jupyter-lab -``` - -or select the environment `spherex_sdt` in your Jupyter Notebook using the dropdown on the upper left. - -### On Fornax - -Installing the conda environment on the [NASA Fornax Science Console](https://science.nasa.gov/astrophysics/programs/physics-of-the-cosmos/community/the-fornax-initiative/) needs slightly different steps. These can be reviewed in the documentation [create a new environment](https://docs.fornax.sciencecloud.nasa.gov/compute-environments/#create-new-env). - -In order to install this specific conda environment on Fornax, the file name of the `yml` file specifically needs to be in the following format `conda-*.yml`. The `yml` file distributed here is already in that format (`conda-spherex_sdt.yml`). Once this is set, open a new terminal (click on the large "+" button right under the "File" menu tab) and type to following command _inside_ the same directory where the `yml` file is located: - ``` - setup-conda-env --user - ``` - -Note that we use the `--user` option here, which will keep the environment available for subsequent Fornax sessions. - -To use the environment in your Jupyter notebooks on Fornax, directly select the environment `spherex_sdt` in your Jupyter Notebook using the dropdown on the upper left. diff --git a/tutorials/spherex/spherex_source_discovery/spherex_source_discovery_tool_demo.md b/tutorials/spherex/spherex_source_discovery/spherex_source_discovery_tool_demo.md index dd7df3f5..442b36c3 100644 --- a/tutorials/spherex/spherex_source_discovery/spherex_source_discovery_tool_demo.md +++ b/tutorials/spherex/spherex_source_discovery/spherex_source_discovery_tool_demo.md @@ -53,6 +53,23 @@ More information is available in the [SPHEREx Explanatory Supplement](https://ir ## 3. Requirements +This directory contains the following files: + +- `spherex_source_discovery_tool_demo.md`: This notebook. +- `conda-spherex_sdt.yml`: Conda environment definition file. +- `spherex_source_discovery_tool/`: Python package, including functions and configuration files to extract and visualize sources. + - Python modules + - `aperture_photometry.py`: Functions adapted from Zafar Rustamkulov to perform aperture photometry. + - `bokeh_viz.py`: Functions for dynamic `bokeh` visualizations. + - `firefly_viz.py`: Functions for `Firefly` visualizations. + - `source_extraction.py`: Functions to work with SExtractor. + - `sdt_utils.py`: Miscellaneous utility functions for working with IRSA tables and SPHEREx image headers. + - SExtractor files + - `default.conv`: Convolution mask used as detection filter. + - `default.nnw`: Table of neural-network weights for star/galaxy separation. + - `default_sdt.sex`: SExtractor default config file. Specifies catalog, extraction, photometry, etc. options. + - `default_sdt.param`: SExtractor default parameters file. Specifies columns to include in output SExtractor catalog. + +++ ### 3.1 On local machines