Hubble Source Catalog API Notebook: SMC Color-Magnitude Diagram¶
August 2019, Rick White¶
A new MAST interface supports queries to the current and previous versions of the Hubble Source Catalog. It allows searches of the summary table (with multi-filter mean photometry) and the detailed table (with all the multi-epoch measurements). It also has an associated API, which is used in this notebook.
This is based on part of HSC Use Case #2.
- It searches the HSC for point-like objects in the Small Magellanic Cloud (SMC) with ACS/WFC V and I band measurements,
- selects a subset of those objects in a V-I color range,
- plots the positions of the objects on the sky, and
- plots the color-magnitude diagram for the selected objects.
The whole process takes only about 2 minutes to complete.
This notebook is available for download. Another simple notebook demonstrates other search capabilities of the API to find variable objects and plot their light curves. A more complex notebook that shows how to access the proper motion tables using the HSC API is also available.
Instructions:¶
- Complete the initialization steps described below.
- Run the notebook.
Running the notebook from top to bottom takes about 2 minutes.
Table of Contents¶
Initialization ¶
Install Python modules¶
This notebook requires the use of Python 3.
This needs the requests and fastkde modules in addition to the common requirements of astropy, numpy and scipy. For anaconda versions of Python the installation commands are:
conda install requests pip install fastkde
In [1]:
%matplotlib inline
import astropy, pylab, time, sys, os, requests, json
import numpy as np
from astropy.table import Table
from astropy.io import ascii
from fastkde import fastKDE
from scipy.interpolate import RectBivariateSpline
from astropy.modeling import models, fitting
# Set page width to fill browser for longer output lines
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
# set width for pprint
astropy.conf.max_width = 150
Useful functions¶
Execute HSC searches and resolve names using MAST query.
In [2]:
hscapiurl = "https://catalogs.mast.stsci.edu/api/v0.1/hsc"
def hsccone(ra,dec,radius,table="summary",release="v3",format="csv",magtype="magaper2",
columns=None, baseurl=hscapiurl, verbose=False,
**kw):
"""Do a cone search of the HSC catalog
Parameters
----------
ra (float): (degrees) J2000 Right Ascension
dec (float): (degrees) J2000 Declination
radius (float): (degrees) Search radius (<= 0.5 degrees)
table (string): summary, detailed, propermotions, or sourcepositions
release (string): v3 or v2
magtype (string): magaper2 or magauto (only applies to summary table)
format: csv, votable, json
columns: list of column names to include (None means use defaults)
baseurl: base URL for the request
verbose: print info about request
**kw: other parameters (e.g., 'numimages.gte':2)
"""
data = kw.copy()
data['ra'] = ra
data['dec'] = dec
data['radius'] = radius
return hscsearch(table=table,release=release,format=format,magtype=magtype,
columns=columns,baseurl=baseurl,verbose=verbose,**data)
def hscsearch(table="summary",release="v3",magtype="magaper2",format="csv",
columns=None, baseurl=hscapiurl, verbose=False,
**kw):
"""Do a general search of the HSC catalog (possibly without ra/dec/radius)
Parameters
----------
table (string): summary, detailed, propermotions, or sourcepositions
release (string): v3 or v2
magtype (string): magaper2 or magauto (only applies to summary table)
format: csv, votable, json
columns: list of column names to include (None means use defaults)
baseurl: base URL for the request
verbose: print info about request
**kw: other parameters (e.g., 'numimages.gte':2). Note this is required!
"""
data = kw.copy()
if not data:
raise ValueError("You must specify some parameters for search")
if format not in ("csv","votable","json"):
raise ValueError("Bad value for format")
url = "{}.{}".format(cat2url(table,release,magtype,baseurl=baseurl),format)
if columns:
# check that column values are legal
# create a dictionary to speed this up
dcols = {}
for col in hscmetadata(table,release,magtype)['name']:
dcols[col.lower()] = 1
badcols = []
for col in columns:
if col.lower().strip() not in dcols:
badcols.append(col)
if badcols:
raise ValueError('Some columns not found in table: {}'.format(', '.join(badcols)))
# two different ways to specify a list of column values in the API
# data['columns'] = columns
data['columns'] = '[{}]'.format(','.join(columns))
# either get or post works
# r = requests.post(url, data=data)
r = requests.get(url, params=data)
if verbose:
print(r.url)
r.raise_for_status()
if format == "json":
return r.json()
else:
return r.text
def hscmetadata(table="summary",release="v3",magtype="magaper2",baseurl=hscapiurl):
"""Return metadata for the specified catalog and table
Parameters
----------
table (string): summary, detailed, propermotions, or sourcepositions
release (string): v3 or v2
magtype (string): magaper2 or magauto (only applies to summary table)
baseurl: base URL for the request
Returns an astropy table with columns name, type, description
"""
url = "{}/metadata".format(cat2url(table,release,magtype,baseurl=baseurl))
r = requests.get(url)
r.raise_for_status()
v = r.json()
# convert to astropy table
tab = Table(rows=[(x['name'],x['type'],x['description']) for x in v],
names=('name','type','description'))
return tab
def cat2url(table="summary",release="v3",magtype="magaper2",baseurl=hscapiurl):
"""Return URL for the specified catalog and table
Parameters
----------
table (string): summary, detailed, propermotions, or sourcepositions
release (string): v3 or v2
magtype (string): magaper2 or magauto (only applies to summary table)
baseurl: base URL for the request
Returns a string with the base URL for this request
"""
checklegal(table,release,magtype)
if table == "summary":
url = "{baseurl}/{release}/{table}/{magtype}".format(**locals())
else:
url = "{baseurl}/{release}/{table}".format(**locals())
return url
def checklegal(table,release,magtype):
"""Checks if this combination of table, release and magtype is acceptable
Raises a ValueError exception if there is problem
"""
releaselist = ("v2", "v3")
if release not in releaselist:
raise ValueError("Bad value for release (must be one of {})".format(
', '.join(releaselist)))
if release=="v2":
tablelist = ("summary", "detailed")
else:
tablelist = ("summary", "detailed", "propermotions", "sourcepositions")
if table not in tablelist:
raise ValueError("Bad value for table (for {} must be one of {})".format(
release, ", ".join(tablelist)))
if table == "summary":
magtypelist = ("magaper2", "magauto")
if magtype not in magtypelist:
raise ValueError("Bad value for magtype (must be one of {})".format(
", ".join(magtypelist)))
def mastQuery(request, url='https://mast.stsci.edu/api/v0/invoke'):
"""Perform a MAST query.
Parameters
----------
request (dictionary): The MAST request json object
url (string): The service URL
Returns the returned data content
"""
# Encoding the request as a json string
requestString = json.dumps(request)
r = requests.post(url, data={'request': requestString})
r.raise_for_status()
return r.text
def resolve(name):
"""Get the RA and Dec for an object using the MAST name resolver
Parameters
----------
name (str): Name of object
Returns RA, Dec tuple with position
"""
resolverRequest = {'service':'Mast.Name.Lookup',
'params':{'input':name,
'format':'json'
},
}
resolvedObjectString = mastQuery(resolverRequest)
resolvedObject = json.loads(resolvedObjectString)
# The resolver returns a variety of information about the resolved object,
# however for our purposes all we need are the RA and Dec
try:
objRa = resolvedObject['resolvedCoordinate'][0]['ra']
objDec = resolvedObject['resolvedCoordinate'][0]['decl']
except IndexError as e:
raise ValueError("Unknown object '{}'".format(name))
return (objRa, objDec)
Find objects in the SMC ¶
This is based on HSC Use Case #2, which includes an example of creating a color-magnitude diagram for the SMC using MAST CasJobs. This is simple to do using the HSC API.
Use MAST name resolver to get position of the SMC ¶
In [3]:
target = 'SMC'
ra, dec = resolve(target)
print(target,ra,dec)
Select objects with the desired magnitudes and colors near the SMC ¶
This searches the summary table for objects in a 3x3 degree box centered on the galaxy that have measurements in both ACS F555W and F814W. It computes the V-I color and selects only objects in the range -1.5 < V-I < 1.5. This large query ultimately returns more than 700,000 objects and takes about a minute to complete.
In [4]:
# save typing a quoted list of columns
columns = """MatchID,MatchRA,MatchDec,CI,A_F555W,A_F814W""".split(",")
columns = [x.strip() for x in columns]
columns = [x for x in columns if x and not x.startswith('#')]
# select objects with at least one ACS F555W and ACS F814W measurement
# and with concentration index 0.9 < CI < 1.6, consistent with point sources
# search a large 3x3 degree box in RA and Dec centered on the SMC
ddec = 1.5
dra = ddec/np.cos(np.radians(dec))
constraints = {'A_F555W_N.gte': 1, 'A_F814W_N.gte': 1, 'CI.gt':0.9, 'CI.lt':1.6,
'MatchDec.gt': dec-ddec, 'MatchDec.lt': dec+ddec,
'MatchRA.gt': ra-dra, 'MatchRA.lt': ra+dra}
# do a search with a large number of rows allowed
t0 = time.time()
tab = ascii.read(hscsearch(table="summary",release='v3',
columns=columns,verbose=True,pagesize=2000000,**constraints))
print("{:.1f} s: retrieved data and converted to {}-row astropy table".format(time.time()-t0, len(tab)))
# compute color column and select for objects in more limited color range
tab['V-I'] = tab['A_F555W'] - tab['A_F814W']
tab = tab[(tab['V-I'] < 1.5) & (tab['V-I'] > -1.5)]
print("{:.1f} s: selected {} objects with -1.5 < V-I < 1.5".format(time.time()-t0, len(tab)))
# clean up the output format
tab['A_F555W'].format = "{:.3f}"
tab['A_F814W'].format = "{:.3f}"
tab['V-I'].format = "{:.3f}"
tab['CI'].format = "{:.3f}"
tab['MatchRA'].format = "{:.6f}"
tab['MatchDec'].format = "{:.6f}"
tab
Out[4]:
Plot object positions on the sky ¶
We mark the galaxy center as well. These fields are sprinkled all over the galaxy (as determined by the HST proposals).
In [5]:
pylab.rcParams.update({'font.size': 16})
pylab.figure(1,(10,10))
pylab.plot(tab['MatchRA'], tab['MatchDec'], 'bo', markersize=1,
label='{} HSC measurements'.format(len(tab)))
pylab.plot(ra,dec,'rx',label=target,markersize=10)
pylab.gca().invert_xaxis()
pylab.gca().set_aspect(1.0/np.cos(np.radians(dec)))
pylab.xlabel('RA [deg]')
pylab.ylabel('Dec [deg]')
pylab.legend(loc='best')
Out[5]:
Plot the color-magnitude diagram ¶
This uses the fastkde module to get a kernel density estimate in order to plot a dense scatterplot.
In [7]:
# Calculate the point density
t0 = time.time()
x = tab['V-I']
y = tab['A_F555W']
myPDF,axes = fastKDE.pdf(x,y,numPoints=2**9+1)
print("kde took {:.1f} sec".format(time.time()-t0))
# interpolate to get z values at points
finterp = RectBivariateSpline(axes[1],axes[0],myPDF)
z = finterp(y,x,grid=False)
# Sort the points by density, so that the densest points are plotted last
idx = z.argsort()
xs, ys, zs = x[idx], y[idx], z[idx]
# select a random subset of points in the most crowded regions to speed up plotting
wran = np.where(np.random.random(len(zs))*zs<0.05)[0]
print("Plotting {} of {} points".format(len(wran),len(zs)))
xs = xs[wran]
ys = ys[wran]
zs = zs[wran]
pylab.rcParams.update({'font.size': 16})
pylab.figure(1,(12,10))
pylab.scatter(xs, ys, c=zs, s=2, edgecolor='', cmap='plasma')
pylab.ylabel('V [mag]')
pylab.xlabel('V - I [mag]')
pylab.xlim(-1.5,1.5)
pylab.ylim(14,27)
pylab.gca().invert_yaxis()
pylab.title('{:,} stars in the Small Magellanic Cloud'.format(len(tab)))
pylab.colorbar()
pylab.tight_layout()
pylab.savefig("smc_cmd.png")
In [ ]: