Using Firefly visualization tools in Python to vet SEDs

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 (https://doi.org/10.3847/1538-3881/ace32f).

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

  • jupytext - for working with Markdown notebooks. It must be installed before your Jupyterlab session is started. See these instructions for helpful information.

To run the notebook by itself, run the cell below (skip if you already have the requirements installed):

# Install libraries if necessary by uncommenting the following
# !pip install -r requirements.txt

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"]
\[[2.5447322 \times 10^{-12},~2.6163133 \times 10^{-12},~1.8539516 \times 10^{-12},~1.2202778 \times 10^{-12},~1.0192029 \times 10^{-12},~1.0054718 \times 10^{-12},~1.0843915 \times 10^{-12},~6.6948035 \times 10^{-13},~2.7315035 \times 10^{-13},~1.3028017 \times 10^{-12},~1.214671 \times 10^{-12},~9.6529632 \times 10^{-13},~0,~0] \; \mathrm{\frac{erg}{s\,cm^{2}}}\]
phot_tbl["lam_flam_err"]
\[[1.4156462 \times 10^{-13},~1.4458281 \times 10^{-13},~1.1047865 \times 10^{-13},~3.0345771 \times 10^{-14},~2.5345456 \times 10^{-14},~4.5377611 \times 10^{-14},~1.8576965 \times 10^{-13},~3.2293075 \times 10^{-15},~1.8791329 \times 10^{-14},~1.4897379 \times 10^{-14},~1.67813 \times 10^{-14},~9.7797784 \times 10^{-15},~0,~0] \; \mathrm{\frac{erg}{s\,cm^{2}}}\]

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

phot_tbl
QTable length=14
radecwavelengthmagmagerrupper_limitzeropointlam_flamlam_flam_err
degdegmicronJyerg / (s cm2)erg / (s cm2)
float64float64float64float64float64float64float64float64float64
81.90154234.8279531.23515.4549999237060550.060400001704692840.01594.02.544732245022997e-121.4156461879395113e-13
81.90154234.8279531.66214.6219997406005860.059999998658895490.01024.02.6163132694737926e-121.4458281115775113e-13
81.90154234.8279532.15914.2460002899169920.064699999988079070.0666.71.8539516327037902e-121.1047865042002163e-13
81.901499834.82795613.413.373999595642090.0270000007003545760.0309.541.2202778434925405e-123.0345771361898645e-14
81.901499834.82795614.612.602000236511230.0270000007003545760.0171.791.0192029131663901e-122.5345456151042028e-14
81.901499834.827956112.09.7399997711181640.048999998718500140.031.6761.0054717723991499e-124.537761137846731e-14
81.901499834.827956122.07.5539999008178710.18600000441074370.08.36351.0843915457289992e-121.8576964773270986e-13
81.9015502951508834.827944193550410.62218.4158172607421880.0052371605802211640.03228.756.694803478683729e-133.229307516617325e-15
81.9015502951508834.827944193550410.51119.7061805725097660.074693025887970170.03552.012.731503451179333e-131.879132860246366e-14
81.9015502951508834.827944193550410.776917.1973991394042970.0124152607381963780.02554.951.3028017437873149e-121.4897378634343012e-14
81.901532234.82797043.413.3789997100830080.0149999996647238730.0309.541.2146710480553427e-121.6781300313770818e-14
81.901532234.82797044.612.661000251770020.0109999999403953550.0171.799.652963193161815e-139.779778373796877e-15
81.9015479458268834.827977077408653.44388652343750.014832420349.1210940.0309.540.00.0
81.9015479458268834.827977077408654.68759030273437.545745952606.201170.0171.790.00.0

Upload to Firefly#

This notebook assumes that jupyter_firefly_extensions has been installed, and that the Firefly server to use has been specified before Jupyterlab was started.

If everything has been properly configured, executing the next cell will bring up a Firefly tab in Jupyterlab with the message “Firefly Ready”:

# Uncomment to use the jupyter_firefly_extensions
#fc = FireflyClient.make_lab_client()

# Uncomment to use a separate web browser tab
fc = FireflyClient.make_client(url="https://irsa.ipac.caltech.edu/irsaviewer")
Open your web browser to this link

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:

# Temporary work-around to Firefly error for Slate and/or Internet connection problems
fc.display_url()
Open your web browser to this link

Write the photometry table to a FITS file:

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

Upload the table 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")
fc.show_fits(fc.upload_file(h_fname), plot_id="2MASS H")
fc.show_fits(fc.upload_file(k_fname), plot_id="2MASS K")

{'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}

In the future, the user should be able to implement turning on “Align and Lock by WCS” using the Python API as part of this notebook.

The current alternative is to use the following link to launch IRSA Finder Chart in a separate browser and turn on that function there: https://irsa.ipac.caltech.edu/applications/finderchart/servlet/api?mode=getResult&locstr=270.599026+-24.147018+Equ+J2000&subsize=0.0167&sources=2MASS&DoSearch=true

In that Firefly display, it should be possible to turn on “Align and Lock by WCS” in the image toolbar.

About This Notebook#

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

Updated On: 2024-07-31

Contact: irsasupport@ipac.caltech.edu or https://irsa.ipac.caltech.edu/docs/help_desk.html

Contents