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