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:

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)
SMC 13.18659 -72.8286

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
https://catalogs.mast.stsci.edu/api/v0.1/hsc/v3/summary/magaper2.csv?pagesize=2000000&A_F555W_N.gte=1&A_F814W_N.gte=1&CI.gt=0.9&CI.lt=1.6&MatchDec.gt=-74.3286&MatchDec.lt=-71.3286&MatchRA.gt=8.10582570497584&MatchRA.lt=18.267354295024163&columns=%5BMatchID%2CMatchRA%2CMatchDec%2CCI%2CA_F555W%2CA_F814W%5D
33.5 s: retrieved data and converted to 764094-row astropy table
33.5 s: selected 761835 objects with -1.5 < V-I < 1.5
Out[4]:
Table length=761835
MatchIDMatchRAMatchDecCIA_F555WA_F814WV-I
int64float64float64float64float64float64float64
4266069114.828002-72.1940000.98023.28422.7190.565
4267110412.999241-72.9384521.22920.96921.160-0.191
4277587611.958159-73.3921051.26625.77925.3360.444
4278242711.729423-73.4985761.34224.93124.5470.384
13562113.475491-72.2015341.05722.53522.562-0.027
13563514.703467-72.2016900.99622.66622.5660.100
13642313.194862-73.0384411.19818.99719.569-0.572
14023313.213482-73.0442561.25225.25724.9560.301
15022513.173528-73.0762541.03523.14822.9730.175
15066113.420397-73.3995271.07824.83524.3040.531
.....................
9871277213.903761-72.8642240.99324.12323.8170.306
9879712211.928429-73.3632491.06423.09222.9400.153
9881963411.681767-73.4876541.13625.26025.0630.197
9883448811.927393-73.4931891.04023.47323.2350.238
9886105314.003072-72.8393760.94921.48820.9740.514
9889957012.983860-72.9903621.14025.25225.1380.115
10198442414.613297-72.3051361.02421.55721.578-0.021
10482383812.772684-73.1520881.03523.73523.6110.124
10487932412.393538-73.2787010.95723.63023.3940.236
10804038613.237029-73.0794231.23023.25623.1150.141

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]:
<matplotlib.legend.Legend at 0x1314d2d940>

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")
kde took 7.1 sec
Plotting 211007 of 761835 points
In [ ]: