Using Firefly visualization tools in Python to vet SEDs

Learning Goals¶

By the end of this tutorial, the following goals can be accomplished:

Select a position and categories of data for search of possible nearby sources
Pull photometry for different bands via TAP queries on chosen IRSA catalogs
Add data columns from result to table, including energy density values (converted from magnitudes)
Create a simple SED plot using compiled photometric data
Display matching image cutouts from IRSA data sets for selected position

Introduction¶

The purpose of this notebook is to show how to make and vet Spectral Energy Distributions (SEDs) using Firefly. All examples use the python package pyvo to search for catalogs and images from IRSA.

Detailed specifications for this notebook come from a paper led by Luisa Rebull on the selection of possible Young Stellar Objects (YSOs) in IC 417 (Rebull et al. (2023)).

In this case, the SEDs generated by photometry data (from several different survey catalogs) near a given position can be compared with corresponding images to verify the location of a point source that could also be considered as a YSO if the graph shape is sound.

Imports¶

Prior to starting Jupyterlab, you should specify the Firefly server by setting the FIREFLY_URL environment variable (use the value https://irsa.ipac.caltech.edu/irsaviewer if you’re not sure), or by one of the other ways listed on the jupyter_firefly_extension page.

The packages needed for this notebook are in the requirements.txt file. They can be installed with pip install -r requirements.txt before you start your Jupyterlab session.

astropy.units - for specifying and manipulating units
astropy.constants - for constants
astropy.coordinates.Skycoord - for representing coordinates
astropy.table.QTable - for tables with units
firefly_client.FireflyClient - Python API to Firefly for displaying tables, images and charts
numpy - for working with arrays
pyvo - for queries to Virtual Observatory services at the archives

# Uncomment the next line to install dependencies if needed.
# !pip install numpy astropy pyvo firefly_client

Call the imports below:

import astropy.units as u
from astropy import constants as const
from astropy.coordinates import SkyCoord
from astropy.table import QTable
from astropy.utils.data import download_file
from firefly_client import FireflyClient
import numpy as np
import pyvo

Choose Target as Example¶

From Figure 10 of the referenced paper, pick the source in the upper right corner (J052736.37+344940.6):

target = SkyCoord(ra="05h27m36.37s", dec="+34d49m40.6s")

Get Photometry from IRSA¶

The available data at IRSA include the following catalogs (with matching radius at 1 arcsec except where noted by reference paper):

Gaia DR3
2MASS
Spitzer-IRAC (GLIMPSE? SEIP Source List?)
Spitzer-MIPS-24 (GLIMPSE? SEIP Source List?)
AllWISE
CatWISE
unWISE
MSX
Akari
Herschel-PACS (70um and 160um)

The NITARP Wiki page on central wavelengths may be helpful for finding zero-point information for specific catalogs, as well as the more general astroquery service on SVO filter data.

Whlie this notebook will focus on a subset of the above list (2MASS, Gaia DR3, and WISE data), note that there is a section where additional data can be included as well.

The variable phot_tbl is the table containing the photometry for the target and all the derived quantities needed to make SED plots:

Right ascension (in degrees)
Declination (in degrees)
Wavelength for the bandpass (in microns)
Photometry measurement
Photometry uncertainty
Is the measurement an upper limit? (Note: 1 means yes, 0 means no)
Zeropoint for the bandpass in units of Janskys (needed for converting photometric measurements to energy densities)

Set the variables for the photmetry table and each quantity in that data structure:

phot_tbl = QTable(names=['ra', 'dec', 'wavelength', 'mag', 'magerr', 'upper_limit', 'zeropoint'],
                  units=[u.deg, u.deg, u.micron, None, None, None, u.Jy])

We query catalog data using IRSA’s TAP (VO Table Access Protocol) service. It’s sufficient to set up the service once, at the beginning of the session:

irsa_tap = pyvo.dal.TAPService("https://irsa.ipac.caltech.edu/TAP")

Retrieve 2MASS photometry¶

To start, define a function to retrieve 2MASS photometry via a TAP query to IRSA that includes upper limits (note that only the ph_qual column on photometric quality of source is used and not the rd_flg column on “read flag”):

def get_2mass_phot(target, radius=1*u.arcsec):
    """Get photometry from 2MASS

    Parameters:
    -----------
    target: astropy.coordinates.SkyCoord
          coordinates of the target(s)

    radius: astropy.units.Quantity
          matching radius, default 1*u.arcsec

    Returns:
    --------
    rows: list
         lists of photometry points to add to table
    """

    query = f"""
    SELECT ra,dec,j_m,j_msigcom,h_m,h_msigcom,k_m,k_msigcom,ph_qual,rd_flg
    FROM fp_psc
    WHERE CONTAINS(POINT('ICRS', ra, dec),
    CIRCLE('ICRS', {target.ra.deg}, {target.dec.deg},
    {radius.to('deg').value}))=1
    """

    results = irsa_tap.search(query)

    rowj = [results["ra"]*u.deg,
            results["dec"]*u.deg,
            1.235*u.micron,
            results["j_m"][0],
            results["j_msigcom"][0],
            1 if results["ph_qual"] == "U" else 0,
            1594*u.Jy]

    rowh = [results["ra"]*u.deg,
            results["dec"]*u.deg,
            1.662*u.micron,
            results["h_m"][0],
            results["h_msigcom"][0],
            1 if results["ph_qual"] == "U" else 0,
            1024*u.Jy]

    rowks = [results["ra"]*u.deg,
            results["dec"]*u.deg,
            2.159*u.micron,
            results["k_m"][0],
            results["k_msigcom"][0],
            1 if results["ph_qual"] == "U" else 0,
            666.7*u.Jy]

    return([rowj, rowh, rowks])

Consider that while the search radius is currently set at 1 degree from the position, that value can be adjusted by changing what the number on the radius line is.

Now gather the rows of 2MASS data and add each of them to the photometry table:

rows = get_2mass_phot(target)

for row in rows:
    phot_tbl.add_row(row)

Retrieve AllWISE Photometry¶

According to the paper referenced above, the detections from the AllWISE catalog were retained if the data quality flags were A, B, or C; if the data quality flag was Z, then the data were provisionally retained with a very large error bar, 30% larger than what appears in the catalog.

Based on that information, define a function to release AllWISE photometry via a TAP query to IRSA that includes the needed data quality flags and retains upper limits:

def get_allwise_phot(target, radius=1*u.arcsec):
    """Get photometry from AllWISE

    Parameters:
    -----------
    target: astropy.coordinates.SkyCoord
          coordinates of the target(s)

    radius: astropy.units.Quantity
          matching radius, defaults to 1*u.arcsec

    Returns:
    --------
    rows: list
         lists of photometry points to add to table
    """

    query = f"""
    SELECT ra,dec,w1mpro,w1sigmpro,w2mpro,w2sigmpro,w3mpro,
          w3sigmpro,w4mpro,w4sigmpro,ph_qual
    FROM allwise_p3as_psd
    WHERE CONTAINS(POINT('ICRS', ra, dec),
    CIRCLE('ICRS', {target.ra.deg}, {target.dec.deg},
    {radius.to('deg').value}))=1
    """

    results = irsa_tap.search(query)

    roww1 = [results["ra"]*u.deg,
             results["dec"]*u.deg,
             3.4*u.micron,
             results["w1mpro"][0],
             results["w1sigmpro"][0],
             1 if results["ph_qual"] == "U" else 0,
             309.54*u.Jy]

    roww2 = [results["ra"]*u.deg,
             results["dec"]*u.deg,
             4.6*u.micron,
             results["w2mpro"][0],
             results["w2sigmpro"][0],
             1 if results["ph_qual"] == "U" else 0,
             171.79*u.Jy]

    roww3 = [results["ra"]*u.deg,
             results["dec"]*u.deg,
             12*u.micron,
             results["w3mpro"][0],
             results["w3sigmpro"][0],
             1 if results["ph_qual"] == "U" else 0,
             31.676*u.Jy]

    roww4 = [results["ra"]*u.deg,
             results["dec"]*u.deg,
             22*u.micron,
             results["w4mpro"][0],
             results["w4sigmpro"][0],
             1 if results["ph_qual"] == "U" else 0,
             8.3635*u.Jy]

    return([roww1, roww2, roww3, roww4])

rows = get_allwise_phot(target)

for row in rows:
    phot_tbl.add_row(row)

Retrieve Gaia DR3 Photometry¶

Define a function to retrieve Gaia DR3 photometry via a TAP query to IRSA (not considering upper limits as Gaia DR3 ony provides “detections or not” info so actual limits are not given):

def get_gaia_phot(target, radius=1*u.arcsec):
    """Get photometry from Gaia DR3

    Parameters:
    -----------
    target: astropy.coordinates.SkyCoord
          coordinates of the target(s)

    radius: astropy.units.Quantity
          matching radius, default is 1*u.arcsec

    Returns:
    --------
    rows: list
         lists of photometry points to add to table
    """

    query = f"""
    SELECT designation,ra,ra_error,dec,dec_error,ra_dec_corr,phot_g_mean_flux_over_error,phot_g_mean_mag,
       phot_bp_mean_flux_over_error,phot_bp_mean_mag,phot_rp_mean_flux_over_error,phot_rp_mean_mag
    FROM gaia_dr3_source
    WHERE CONTAINS(POINT('ICRS', ra, dec),
    CIRCLE('ICRS', {target.ra.deg}, {target.dec.deg},
    {radius.to('deg').value}))=1
    """

    results = irsa_tap.search(query)

    rowg = [results["ra"]*u.deg,
            results["dec"]*u.deg,
            0.622*u.micron,
            results["phot_g_mean_mag"][0],
            2.5/(np.log(10)*results["phot_g_mean_flux_over_error"][0]),
            0,
            3228.75*u.Jy]

    rowbp = [results["ra"]*u.deg,
            results["dec"]*u.deg,
            0.511*u.micron,
            results["phot_bp_mean_mag"][0],
            2.5/(np.log(10)*results["phot_bp_mean_flux_over_error"][0]),
            0,
            3552.01*u.Jy]

    rowrp = [results["ra"]*u.deg,
            results["dec"]*u.deg,
            0.7769*u.micron,
            results["phot_rp_mean_mag"][0],
            2.5/(np.log(10)*results["phot_rp_mean_flux_over_error"][0]),
            0,
            2554.95*u.Jy]

    return([rowg, rowbp, rowrp])

rows = get_gaia_phot(target)

for row in rows:
    phot_tbl.add_row(row)

Retrieve CatWISE Photometry¶

Define a function to retrieve CatWISE photometry via a TAP query to IRSA that includes upper limits (note that CatWISE2020 Catalog is being used rather than CatWISE Preliminary Catalog):

def get_catwise_phot(target, radius=1*u.arcsec):
    """Get photometry from CatWISE

    Parameters:
    -----------
    target: astropy.coordinates.SkyCoord
          coordinates of the target(s)

    radius: astropy.units.Quantity
          matching radius, defaults to 1*u.arcsec

    Returns:
    --------
    rows: list
         lists of photometry points to add to table
    """

    query = f"""
    SELECT ra,dec,w1mpro,w1sigmpro,w2mpro,w2sigmpro
    FROM catwise_2020
    WHERE CONTAINS(POINT('ICRS', ra, dec),
    CIRCLE('ICRS', {target.ra.deg}, {target.dec.deg},
    {radius.to('deg').value}))=1
    """

    results = irsa_tap.search(query)

    roww1 = [results["ra"]*u.deg,
             results["dec"]*u.deg,
             3.4*u.micron,
             results["w1mpro"][0],
             results["w1sigmpro"][0],
             1 if results["w1sigmpro"] is None else 0,
             309.54*u.Jy]

    roww2 = [results["ra"]*u.deg,
             results["dec"]*u.deg,
             4.6*u.micron,
             results["w2mpro"][0],
             results["w2sigmpro"][0],
             1 if results["w2sigmpro"] is None else 0,
             171.79*u.Jy]

    return([roww1, roww2])

rows = get_catwise_phot(target)

for row in rows:
    phot_tbl.add_row(row)

Retrieve unWISE Photometry¶

Define a function to retrieve unWISE photometry via a TAP query to IRSA(note that flux units are given in nano-mag, zero-point is carried over from AllWISE, and no flag for upper limits is given):

def get_unwise_phot(target, radius=1*u.arcsec):
    """Get photometry from CatWISE

    Parameters:
    -----------
    target: astropy.coordinates.SkyCoord
          coordinates of the target(s)

    radius: astropy.units.Quantity
          matching radius, defaults to 1*u.arcsec

    Returns:
    --------
    rows: list
         lists of photometry points to add to table
    """

    query = f"""
    SELECT ra,dec,flux_1,dflux_1,flux_2,dflux_2,fluxlbs_1
    FROM unwise_2019
    WHERE CONTAINS(POINT('ICRS', ra, dec),
    CIRCLE('ICRS', {target.ra.deg}, {target.dec.deg},
    {radius.to('deg').value}))=1
    """

    results = irsa_tap.search(query)

    roww1 = [results["ra"]*u.deg,
             results["dec"]*u.deg,
             3.4*u.micron,
             pow(10,9)*results["flux_1"][0],
             pow(10,9)*results["dflux_1"][0],
             0,
             309.54*u.Jy]

    roww2 = [results["ra"]*u.deg,
             results["dec"]*u.deg,
             4.6*u.micron,
             pow(10,9)*results["flux_2"][0],
             pow(10,9)*results["dflux_2"][0],
             0,
             171.79*u.Jy]

    return([roww1, roww2])

rows = get_unwise_phot(target)

for row in rows:
    phot_tbl.add_row(row)

Retrieve Optional Photometry¶

Insert more data retrievals here, following the above formats.

Note: After more columns are added in the conversion to spectral energy density, the phot_tbl.add_row functions will thrown an exception that number of values doesn’t match the number of columns.

# SEIP
# GLIMPSE

# MSX, AKARI, Herschel PACS

# Pan-STARRS, UKIDSS?

Convert Magnitudes to Energy Densities¶

Use the zeropoint information embedded by the catalog retrieval functions to compute energy densities:

factor = (const.c*phot_tbl["wavelength"]**-1).to(u.Hz)

phot_tbl["lam_flam"] = ((factor*phot_tbl["zeropoint"] *
                            10**(-0.4*phot_tbl["mag"]))
                       ).to(u.erg*u.s**-1*u.cm**-2)

As the above formula was for the photometry measurements given, compute the photometry uncertainty using a similar formula:

phot_tbl["lam_flam_err"] = phot_tbl["lam_flam"]*phot_tbl["magerr"]*np.log(10)/2.5

Check the energy density values for each converted value below:

phot_tbl["lam_flam"]

phot_tbl["lam_flam_err"]

Run the following cell to view the full photometry table and its values for each column (ex. zeropoint, upper limit):

phot_tbl

Upload to Firefly¶

There are two ways to initialize a Firefly client from Python, depending on whether you’re running the notebook in JupyterLab or not. Assuming you have jupyter-firefly-extensions set up in your environment as explained here, you can use make_lab_client() in JupyterLab, which will open the Firefly viewer in a new tab within the Lab. Otherwise, you can use make_client() in a Jupyter Notebook (or even a Python shell), which will open the Firefly viewer in a new web browser tab.

You also need a Firefly server to communicate with your Firefly Python client. In this notebook, we use a public Firefly server: the IRSA Viewer (https://irsa.ipac.caltech.edu/irsaviewer). However, you can also run a local Firefly server via a Firefly Docker image and access it at http://localhost:8080/firefly. The URL of the Firefly server is read by both make_client() and make_lab_client() through the environment variable FIREFLY_URL. However, make_client() also allows you to pass the URL directly as the url parameter.

# Uncomment when using within Jupyter Lab with jupyter_firefly_extensions installed
# fc = FireflyClient.make_lab_client()

# Uncomment for contexts other than above 
fc = FireflyClient.make_client(url="https://irsa.ipac.caltech.edu/irsaviewer")

fc.reinit_viewer() # to clean the state, if this cell ran earlier

In the event that there are problems with the tab opened above, run the below command to obtain a web link that can be opened in a browser directly:

fc.display_url()

Write the photometry table to a FITS file:

tblpath = "./phot_tbl.fits"
phot_tbl.write(tblpath, format="fits", overwrite=True)

Upload the table (FITS file) to Firefly:

tval = fc.upload_file(tblpath)

Show the table in Firefly:

fc.show_table(tval, tbl_id="photometry", title="SED Photometry")

{'success': True}

Specify the columns and markers for the SED plot:

trace1 = dict(tbl_id="photometry", x="tables::wavelength", y="tables::lam_flam", mode="markers",
            type="scatter", error_y=dict(array="tables::lam_flam_err"),
            marker=dict(size=4))
trace_data = [trace1]

Specify the title and axes for the SED plot:

layout_s = dict(title="SED", xaxis=dict(title="Wavelength (microns)", type="log"),
                yaxis=dict(title="Energy density (erg/s/cm2)", type="log"))

Display the SED plot in Firefly:

fc.show_chart(layout=layout_s, data=trace_data)

{'success': True,
 'warning': 'group_id unnecessary and ignored in triview mode'}

Retrieve and Display Images¶

Use queries directed at Simple Image Access services to retrieve images:

irsa_sia = pyvo.dal.SIAService("https://irsa.ipac.caltech.edu/cgi-bin/2MASS/IM/nph-im_sia?type=at&ds=asky&")

Get 2MASS Images¶

def get_2mass_images(target, search_size=1*u.arcsec):
    """Retrieve 2MASS images

    Parameters:
    -----------
    target: astropy.coordinates.SkyCoord
        coordinates of the target(s)

    search_size: astropy.units.Quantity
        matching radius, defaults to 1*u.arcsec

    Returns:
    --------
    j_filename: `str`
         path to J-band image in FITS format

    h_filename: `str`
         path to H-band image in FITS format

    k_filename: `str`
         path to K-band image in FITS format
    """
    im_table = irsa_sia.search(pos=target, size=search_size)
    url_j = url_h = url_k = None
    # Get URL for the first instance of each band
    for i in range(len(im_table)):
        if url_j is None and im_table[i]["band"] == "J":
            url_j = im_table[i].getdataurl()
        if url_h is None and im_table[i]["band"] == "H":
            url_h = im_table[i].getdataurl()
        if url_k is None and im_table[i]["band"] == "K":
            url_k = im_table[i].getdataurl()
    j_filename = download_file(url_j, cache=True)
    h_filename = download_file(url_h, cache=True)
    k_filename = download_file(url_k, cache=True)
    return (j_filename, h_filename, k_filename)

j_fname, h_fname, k_fname = get_2mass_images(target)

fc.show_fits(fc.upload_file(j_fname), plot_id="2MASS J", title="2MASS J")
fc.show_fits(fc.upload_file(h_fname), plot_id="2MASS H", title="2MASS H")
fc.show_fits(fc.upload_file(k_fname), plot_id="2MASS K", title="2MASS K")

WARNING: show_fits() is deprecated. Use show_fits_image() instead.

WARNING: show_fits() is deprecated. Use show_fits_image() instead.

WARNING: show_fits() is deprecated. Use show_fits_image() instead.

{'success': True}

fc.set_zoom(plot_id="2MASS J", factor=4.0)
fc.set_pan(plot_id="2MASS J", x=target.ra.deg,
           y=target.dec.deg, coord="j2000")
fc.set_zoom(plot_id="2MASS H", factor=4.0)
fc.set_pan(plot_id="2MASS H", x=target.ra.deg,
           y=target.dec.deg, coord="j2000")
fc.set_zoom(plot_id="2MASS K", factor=4.0)
fc.set_pan(plot_id="2MASS K", x=target.ra.deg,
           y=target.dec.deg, coord="j2000")

{'success': True}

For turning on “Align and Lock by WCS” in these images, use:

fc.align_images(lock_match=True)

{'success': True}

Note that zooming or panning one image will do the same in the other images as well.

About This Notebook¶

Author: David Shupe (IRSA Scientist) and Joyce Kim (IRSA Scientist) in conjunction with the IRSA Science Team

Updated: 2024-12-19

Contact: the IRSA Helpdesk with questions or reporting problems.

References¶

Rebull, L. M., Anderson, R. L., Hall, G., Kirkpatrick, J. D., Koenig, X., Odden, C. E., Rodriguez, B., Sanchez, R., Senson, B., Urbanowski, V., Austin, M., Blood, K., Kerman, E., Long, J., & Roosa, N. (2023). Young Stellar Object Candidates in IC 417. The Astronomical Journal, 166(3), 87. 10.3847/1538-3881/ace32f

Special Topics

Access HATS Collections Using LSDB: Euclid Q1 and ZTF DR23

Special Topics

Parallelizing image convolution