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
# 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
| ra | dec | wavelength | mag | magerr | upper_limit | zeropoint | lam_flam | lam_flam_err |
|---|---|---|---|---|---|---|---|---|
| deg | deg | micron | Jy | erg / (s cm2) | erg / (s cm2) | |||
| float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
| 81.901542 | 34.827953 | 1.235 | 15.454999923706055 | 0.06040000170469284 | 0.0 | 1594.0 | 2.544732245022997e-12 | 1.4156461879395113e-13 |
| 81.901542 | 34.827953 | 1.662 | 14.621999740600586 | 0.05999999865889549 | 0.0 | 1024.0 | 2.6163132694737926e-12 | 1.4458281115775113e-13 |
| 81.901542 | 34.827953 | 2.159 | 14.246000289916992 | 0.06469999998807907 | 0.0 | 666.7 | 1.8539516327037902e-12 | 1.1047865042002163e-13 |
| 81.9014998 | 34.8279561 | 3.4 | 13.37399959564209 | 0.027000000700354576 | 0.0 | 309.54 | 1.2202778434925405e-12 | 3.0345771361898645e-14 |
| 81.9014998 | 34.8279561 | 4.6 | 12.60200023651123 | 0.027000000700354576 | 0.0 | 171.79 | 1.0192029131663901e-12 | 2.5345456151042028e-14 |
| 81.9014998 | 34.8279561 | 12.0 | 9.739999771118164 | 0.04899999871850014 | 0.0 | 31.676 | 1.0054717723991499e-12 | 4.537761137846731e-14 |
| 81.9014998 | 34.8279561 | 22.0 | 7.553999900817871 | 0.1860000044107437 | 0.0 | 8.3635 | 1.0843915457289992e-12 | 1.8576964773270986e-13 |
| 81.90155029515088 | 34.82794419355041 | 0.622 | 18.415817260742188 | 0.005237160580221164 | 0.0 | 3228.75 | 6.694803478683729e-13 | 3.229307516617325e-15 |
| 81.90155029515088 | 34.82794419355041 | 0.511 | 19.706180572509766 | 0.07469302588797017 | 0.0 | 3552.01 | 2.731503451179333e-13 | 1.879132860246366e-14 |
| 81.90155029515088 | 34.82794419355041 | 0.7769 | 17.197399139404297 | 0.012415260738196378 | 0.0 | 2554.95 | 1.3028017437873149e-12 | 1.4897378634343012e-14 |
| 81.9015322 | 34.8279704 | 3.4 | 13.378999710083008 | 0.014999999664723873 | 0.0 | 309.54 | 1.2146710480553427e-12 | 1.6781300313770818e-14 |
| 81.9015322 | 34.8279704 | 4.6 | 12.66100025177002 | 0.010999999940395355 | 0.0 | 171.79 | 9.652963193161815e-13 | 9.779778373796877e-15 |
| 81.90154794582688 | 34.82797707740865 | 3.4 | 4388652343750.0 | 14832420349.121094 | 0.0 | 309.54 | 0.0 | 0.0 |
| 81.90154794582688 | 34.82797707740865 | 4.6 | 8759030273437.5 | 45745952606.20117 | 0.0 | 171.79 | 0.0 | 0.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 welcome message.
# 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")
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")
{'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, there’s no dedicated method in Firefly but you can use this workaround (using low-level API):
fc.dispatch('ImagePlotCntlr.wcsMatch',
payload=dict(matchType='Standard', lockMatch=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 On: 2024-10-17
Contact: irsasupport@ipac.caltech.edu or https://irsa.ipac.caltech.edu/docs/help_desk.html