{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Hubble Source Catalog SWEEPS Proper Motion Notebook\n", "### August 2019, Steve Lubow and 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", "The [web-based user interface](https://catalogs.mast.stsci.edu/hsc) to the HSC does not currently include access to the new proper motions available for the [SWEEPS](https://media.stsci.edu/news_release/news/2011-16) field in version 3.1 of the [Hubble Source Catalog](https://archive.stsci.edu/hst/hsc). However those tables are accessible via the API. This notebook shows how to use them.\n", "\n", "This notebook is similar to the [previously released version](../sweeps/sweeps_hscv3p1.html) that uses CasJobs rather than the new API. The Casjobs interface is definitely more powerful and flexible, but the API is simpler to use for users who are not already experienced Casjobs users. Currently the API does not include easy access to the colors and magnitudes of the SWEEPS objects, but they will be added soon.\n", "\n", "Additional information is available on the \n", "[SWEEPS Proper Motions help page](https://archive.stsci.edu/hst/hsc/help/sweeps/hsc_sweeps_pm.html).\n", "\n", "This notebook is available for [download](sweeps_hscv3p1_api.ipynb). A simpler notebook that provides a quick introduction to the HSC API is also [available](hscv3_api.html). Another [simple notebook](hscv3_smc_api.html) generates a color-magnitude diagram for the Small Magellanic Cloud in only a couple of minutes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Instructions: \n", "* Complete the initialization steps [described below](#Initialization).\n", "* Run the notebook to completion. \n", "* Modify and rerun any sections of the Table of Contents below.\n", "\n", "Running the notebook from top to bottom takes about \n", "4 minutes (depending on the speed of your computer).\n", "\n", "\n", "# Table of Contents\n", "* [Initialization](#Initialization)\n", "* [Properties of Full Catalog](#fullcat)\n", " * [Sky Coverage](#SkyCoverage)\n", " * [Proper Motion Distributions](#pmhist)\n", " * [Visit Distribution](#visitshist)\n", " * [Time Coverage Distributions](#timehist)\n", " * [Detection Positions](#detpos)\n", " * [Positions for a Sample With Good PMs](#good_sample)\n", "* [Science Applications](#sciap)\n", " * [High Proper Motion Objects](#hpm)\n", " * [HLA Cutout Images for Selected Objects](#cutouts)" ] }, { "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 some special 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",
    "conda install pillow\n",
    "pip install fastkde\n",
    "
\n", "\n", "#### Other optional configuration\n", "\n", "If desired, you can set `resPath`, the output directory, in the next code block (the default location is the current working directory, which is probably the same directory as this script)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "resPath=\"./\" # directory where generated plots are saved\n", "\n", "%matplotlib inline\n", "import astropy, pylab, time, sys, os, requests, json\n", "import numpy as np\n", "from matplotlib.colors import LogNorm\n", "\n", "## For handling ordinary astropy Tables\n", "from astropy.table import Table\n", "from astropy.io import fits, ascii\n", "\n", "from PIL import Image\n", "from io import BytesIO\n", "\n", "from fastkde import fastKDE\n", "from scipy.interpolate import RectBivariateSpline\n", "from astropy.modeling import models, fitting\n", "\n", "# There are a number of relatively unimportant warnings that \n", "# show up, so for now, suppress them:\n", "import warnings\n", "warnings.filterwarnings(\"ignore\")\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": [ "### Use metadata query to get information on available columns\n", "\n", "This query works for any of the tables in the API (summary, detailed, propermotions, sourcepositions)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "meta = hscmetadata(\"propermotions\")\n", "print(' '.join(meta['name']))\n", "meta[:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Retrieve data on selected SWEEPS objects\n", "\n", "This makes a single large request to the HSC search interface to the get the contents of the proper motions table. Despite its large size (460K rows), the query is relatively efficient: it takes about 25 seconds to retrieve the results from the server, and then another 20 seconds to convert it to an astropy table. The speed of the table conversion will depend on your computer.\n", "\n", "Note that the query restricts the sample to objects with at least 20 images total spread over at least 10 different visits.\n", "These constraints can be modified depending on your science goals." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "columns = \"\"\"ObjID,raMean,decMean,raMeanErr,decMeanErr,NumFilters,NumVisits,\n", " pmLat,pmLon,pmLatErr,pmLonErr,pmLatDev,pmLonDev,epochMean,epochStart,epochEnd\"\"\".split(\",\")\n", "columns = [x.strip() for x in columns]\n", "columns = [x for x in columns if x and not x.startswith('#')]\n", "\n", "# missing -- impossible with current data I think\n", "# MagMed, n, MagMAD\n", "\n", "constraints = {'NumFilters.gt':1, 'NumVisits.gt':9, 'NumImages.gt':19}\n", "\n", "# note the pagesize parameter, which allows retrieving very large results\n", "# the default pagesize is 50000 rows\n", "\n", "t0 = time.time()\n", "results = hscsearch(table=\"propermotions\",release='v3',columns=columns,verbose=True,\n", " pagesize=500000, **constraints)\n", "print(\"{:.1f} s: retrieved data\".format(time.time()-t0))\n", "tab = ascii.read(results)\n", "print(\"{:.1f} s: converted to astropy table\".format(time.time()-t0))\n", "\n", "# change some column names for consistency with the Casjobs version of this notebook\n", "\n", "tab.rename_column(\"raMean\",\"RA\")\n", "tab.rename_column(\"decMean\",\"Dec\")\n", "tab.rename_column(\"raMeanErr\",\"RAerr\")\n", "tab.rename_column(\"decMeanErr\",\"Decerr\")\n", "tab.rename_column(\"pmLat\",\"bpm\")\n", "tab.rename_column(\"pmLon\",\"lpm\")\n", "tab.rename_column(\"pmLatErr\",\"bpmerr\")\n", "tab.rename_column(\"pmLonErr\",\"lpmerr\")\n", "\n", "# compute some additional columns\n", "\n", "tab['pmdev'] = np.sqrt(tab['pmLonDev']**2+tab['pmLatDev']**2)\n", "tab['yr'] = (tab['epochMean'] - 47892)/365.25+1990\n", "tab['dT'] = (tab['epochEnd']-tab['epochStart'])/365.25\n", "tab['yrStart'] = (tab['epochStart'] - 47892)/365.25+1990\n", "tab['yrEnd'] = (tab['epochEnd'] - 47892)/365.25+1990\n", "\n", "# delete some columns that are not needed after the computations\n", "del tab['pmLonDev'], tab['pmLatDev'], tab['epochEnd'], tab['epochStart'], tab['epochMean']\n", "\n", "tab" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Properties of Full Catalog " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sky Coverage " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = tab['RA']\n", "y = tab['Dec']\n", "\n", "pylab.rcParams.update({'font.size': 16})\n", "pylab.figure(1,(12,10))\n", "pylab.scatter(x, y, s=1)\n", "pylab.autoscale(tight=True)\n", "pylab.xlabel('RA')\n", "pylab.ylabel('Dec')\n", "dc=0.01\n", "pylab.xlim(min(x)-dc, max(x)+dc)\n", "pylab.ylim(min(y)-dc, max(y)+dc)\n", "pylab.gca().invert_xaxis()\n", "pylab.text(0.5,0.93,'{:,} stars in SWEEPS'.format(len(x)),\n", " horizontalalignment='left',\n", " transform=pylab.gca().transAxes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Proper Motion Histograms " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Proper motion histograms for lon and lat " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bin = 0.2\n", "hrange = (-20,20)\n", "bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1\n", "pylab.rcParams.update({'font.size': 16})\n", "pylab.figure(1,(12,10))\n", "pylab.hist(tab['lpm'], range=hrange, bins=bincount, label='Longitude', \n", " histtype='step', linewidth=2)\n", "pylab.hist(tab['bpm'], range=hrange, bins=bincount, label='Latitude', \n", " histtype='step', linewidth=2)\n", "pylab.xlabel('Proper motion [mas/yr]')\n", "pylab.ylabel('Number [in {:.2} mas bins]'.format(bin))\n", "pylab.legend(loc='upper right')\n", "pylab.autoscale(enable=True, axis='x', tight=True)\n", "pylab.ylim(0,13500)\n", "pylab.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "pylab.tight_layout()\n", "pylab.savefig('{}sweeps_api_pmerr_hist.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Proper motion error cumulative histogram for lon and lat" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bin = 0.01\n", "hrange = (0,2)\n", "bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1\n", "pylab.rcParams.update({'font.size': 16})\n", "pylab.figure(1,(12,10))\n", "pylab.hist(tab['lpmerr'], range=hrange, bins=bincount, label='Longitude Error', \n", " histtype='step', cumulative=1, linewidth=2)\n", "pylab.hist(tab['bpmerr'], range=hrange, bins=bincount, label='Latitude Error', \n", " histtype='step', cumulative=1, linewidth=2)\n", "pylab.xlabel('Proper motion error [mas/yr]')\n", "pylab.ylabel('Cumulative number [in {:0.2} mas bins]'.format(bin))\n", "pylab.legend(loc='upper right')\n", "pylab.autoscale(enable=True, axis='x', tight=True)\n", "pylab.ylim(0,500000)\n", "pylab.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "pylab.tight_layout()\n", "pylab.savefig('{}sweeps_api_pmerr_cumhist.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Proper motion error log histogram for lon and lat " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bin = 0.01\n", "hrange = (0,6)\n", "bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1\n", "pylab.rcParams.update({'font.size': 16})\n", "pylab.figure(1,(12,10))\n", "pylab.hist(tab['lpmerr'], range=hrange, bins=bincount, label='Longitude Error', \n", " histtype='step', linewidth=2)\n", "pylab.hist(tab['bpmerr'], range=hrange, bins=bincount, label='Latitude Error', \n", " histtype='step', linewidth=2)\n", "pylab.xlabel('Proper motion error [mas/yr]')\n", "pylab.ylabel('Number [in {:0.2} mas bins]'.format(bin))\n", "pylab.legend(loc='upper right')\n", "pylab.yscale('log')\n", "pylab.autoscale(enable=True, axis='x', tight=True)\n", "pylab.ylim(0,15000)\n", "pylab.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "pylab.tight_layout()\n", "pylab.savefig('{}sweeps_api_pmerr_loghist.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Proper motion error as a function of dT\n", "\n", "Exclude objects with dT near zero, and to improve the plotting add a bit of random noise to spread out the quanitized time values." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# restrict to sources with dT > 1 year\n", "dtmin = 1.0\n", "w = np.where(tab['dT']>dtmin)[0]\n", "if ('rw' not in locals()) or len(rw) != len(w):\n", " rw = np.random.random(len(w))\n", "x = np.array(tab['dT'][w]) + 0.5*(rw-0.5)\n", "y = np.log(np.array(tab['lpmerr'][w]))\n", "\n", "# Calculate the point density\n", "t0 = time.time()\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.yscale('log')\n", "pylab.scatter(xs, np.exp(ys), c=zs, s=2, edgecolor='', cmap='plasma', \n", " label='Longitude PM error')\n", "pylab.autoscale(tight=True, axis='y')\n", "pylab.xlim(0.0, max(x)*1.05)\n", "pylab.xlabel('Date range [yr]')\n", "pylab.ylabel('Proper motion error [mas/yr]')\n", "pylab.legend(loc='best')\n", "pylab.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "pylab.colorbar()\n", "pylab.tight_layout()\n", "pylab.savefig('{}sweeps_api_pmerr_vs_dt.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Proper motion error log histogram for lon and lat \n", "\n", "Divide sample into points with $<6$ years of data and points with more than 6 years of data." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "bin = 0.01\n", "hrange = (0,6)\n", "bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1\n", "\n", "tsplit = 6\n", "\n", "pylab.rcParams.update({'font.size': 16})\n", "pylab.figure(1,(12,12))\n", "pylab.subplot(211)\n", "w = np.where(tab['dT']<=tsplit)[0]\n", "pylab.hist(tab['lpmerr'][w], range=hrange, bins=bincount, label='Longitude Error', \n", " histtype='step', linewidth=2)\n", "pylab.hist(tab['bpmerr'][w], range=hrange, bins=bincount, label='Latitude Error', \n", " histtype='step', linewidth=2)\n", "pylab.xlabel('Proper motion error [mas/yr]')\n", "pylab.ylabel('Number [in {:0.2} mas bins]'.format(bin))\n", "pylab.legend(loc='upper right')\n", "pylab.yscale('log')\n", "pylab.autoscale(enable=True, axis='x', tight=True)\n", "pylab.ylim(0,15000)\n", "pylab.title('{:,} stars in SWEEPS with dT < {} yrs'.format(len(w),tsplit))\n", "pylab.tight_layout()\n", "\n", "pylab.subplot(212)\n", "w = np.where(tab['dT']>tsplit)\n", "pylab.hist(tab['lpmerr'][w], range=hrange, bins=bincount, label='Longitude Error', \n", " histtype='step', linewidth=2)\n", "pylab.hist(tab['bpmerr'][w], range=hrange, bins=bincount, label='Latitude Error', \n", " histtype='step', linewidth=2)\n", "pylab.xlabel('Proper motion error [mas/yr]')\n", "pylab.ylabel('Number [in {:0.2} mas bins]'.format(bin))\n", "pylab.legend(loc='upper right')\n", "pylab.yscale('log')\n", "pylab.autoscale(enable=True, axis='x', tight=True)\n", "pylab.ylim(0,15000)\n", "pylab.title('{:,} stars in SWEEPS with dT > {} yrs'.format(len(w),tsplit))\n", "pylab.tight_layout()\n", "\n", "pylab.savefig('{}sweeps_api_pmerr_loghist2.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Number of Visits Histogram " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bin = 1\n", "hrange = (0,130)\n", "bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1\n", "pylab.rcParams.update({'font.size': 16})\n", "pylab.figure(1,(12,10))\n", "pylab.hist(tab['NumVisits'], range=hrange, bins=bincount, label='Number of visits ', \n", " histtype='step', linewidth=2)\n", "pylab.xlabel('Number of visits')\n", "pylab.ylabel('Number of objects')\n", "pylab.autoscale(enable=True, axis='x', tight=True)\n", "pylab.ylim(0,200000)\n", "pylab.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "pylab.tight_layout()\n", "pylab.savefig('{}sweeps_api_numvisits_hist.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Time Histograms \n", "\n", "First plot histogram of observation dates." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bin = 1\n", "hrange = (2000, 2020)\n", "bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1\n", "pylab.rcParams.update({'font.size': 16})\n", "pylab.figure(1,(12,10))\n", "pylab.hist(tab['yr'], range=hrange, bins=bincount, label='year ', histtype='step', linewidth=2)\n", "pylab.xlabel('mean detection epoch (year)')\n", "pylab.ylabel('Number of objects')\n", "pylab.autoscale(enable=True, axis='x', tight=True)\n", "pylab.ylim(0,300000)\n", "pylab.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "pylab.tight_layout()\n", "pylab.savefig('{}sweeps_api_year_hist.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then plot histogram of observation duration for the objects." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "bin = 0.25\n", "hrange = (0, 15)\n", "bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1\n", "pylab.rcParams.update({'font.size': 16})\n", "pylab.figure(1,(12,10))\n", "pylab.hist(tab['dT'], range=hrange, bins=bincount, label='year ', histtype='step', linewidth=2)\n", "pylab.xlabel('time span (years)')\n", "pylab.ylabel('Number of objects')\n", "pylab.autoscale(enable=True, axis='x', tight=True)\n", "pylab.yscale('log')\n", "pylab.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "pylab.tight_layout()\n", "pylab.savefig('{}sweeps_api_year_hist.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Detection Positions \n", "\n", "Define a function to plot the PM fit for an object." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# define function\n", "def positions(Obj):\n", " \"\"\"\n", " input parameter Obj is the value of the ObjID \n", " output plots change in (lon, lat) as a function of time\n", " overplots proper motion fit\n", " provides ObjID and proper motion information in labels\n", " \"\"\"\n", "\n", " # get the measured positions as a function of time\n", " pos = ascii.read(hscsearch(\"sourcepositions\", columns=\"dT,dLon,dLat\".split(','), objid=Obj))\n", " pos.sort('dT')\n", " \n", " # get the PM fit parameters\n", " pm = ascii.read(hscsearch(\"propermotions\", columns=\"pmlon,pmlonerr,pmlat,pmlaterr\".split(','), objid=Obj))\n", " \n", " lpm = pm['pmlon'][0]\n", " bpm = pm['pmlat'][0]\n", " \n", " # get the intercept for the proper motion fit referenced to the start time\n", " # time between mean epoch and zero (ref) epoch (years)\n", "\n", " x = pos['dT']\n", " y = pos['dLon']\n", " pylab.rcParams.update({'font.size':10})\n", " pylab.figure(1,(6,3))\n", " pylab.subplot(121)\n", " pylab.scatter(x, y, s=10)\n", " # xpm = np.linspace(0, max(x), 10)\n", " xpm = np.array([x.min(),x.max()])\n", " ypm = lpm*xpm\n", " pylab.plot(xpm, ypm, '-r')\n", " pylab.xlabel('dT (yrs)')\n", " pylab.ylabel('dLon (mas)')\n", " y = pos['dLat']\n", " pylab.subplot(122)\n", " pylab.scatter(x, y, s=10)\n", " ypm = bpm*xpm\n", " pylab.plot(xpm, ypm, '-r')\n", " pylab.xlabel('dT (yrs)')\n", " pylab.ylabel('dLat (mas)')\n", " pylab.suptitle(\"\"\"ObjID {}\n", "{} detections, (lpm, bpm) = ({:.1f}, {:.1f}) mas/yr\"\"\".format(Obj, len(x), lpm, bpm),\n", " size=10)\n", " pylab.tight_layout(rect=[0, 0.0, 1, 0.88])\n", " pylab.show()\n", " pylab.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot positions for a sample of objects with good proper motions \n", "\n", "This selects objects that are detected in more than 90 visits with a median absolute deviation from the fit of less than 1.5 mas and proper motion error less than 1.0 mas/yr." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "n = tab['NumVisits']\n", "dev = tab['pmdev']\n", "objid = tab['ObjID']\n", "lpmerr0 = np.array(tab['lpmerr'])\n", "bpmerr0 = np.array(tab['bpmerr'])\n", "wi = np.where( (dev < 1.5) & (n > 90) & (np.sqrt(bpmerr0**2+lpmerr0**2) < 1.0))[0]\n", "print(\"Plotting {} objects\".format(len(wi)))\n", "for o in objid[wi]:\n", " positions(o)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Science Applications " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## High Proper Motion Objects \n", "\n", "Get a list of objects with high, accurately measured proper motions.\n", "Proper motions are measured relative to the Galactic center." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "lpm_sgra = -6.379 # +- 0.026\n", "bpm_sgra = -0.202 # +- 0.019\n", "\n", "lpm0 = np.array(tab['lpm'])\n", "bpm0 = np.array(tab['bpm'])\n", "lpmerr0 = np.array(tab['lpmerr'])\n", "bpmerr0 = np.array(tab['bpmerr'])\n", "pmtot0 = np.sqrt((bpm0-bpm_sgra)**2+(lpm0-lpm_sgra)**2)\n", "pmerr0 = np.sqrt(bpmerr0**2+lpmerr0**2)\n", "dev = tab['pmdev']\n", "\n", "# sort samples by decreasing PM\n", "wpmh = np.where((pmtot0 > 15) & (pmerr0 < 1.0) & (dev < 5))[0]\n", "wpmh = wpmh[np.argsort(-pmtot0[wpmh])]\n", "\n", "print(\"Plotting {} objects\".format(len(wpmh)))\n", "for o in tab[\"ObjID\"][wpmh]:\n", " positions(o)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Get HLA cutout images for selected objects \n", "\n", "Get HLA color cutout images for the high-PM objects. The `query_hla` function gets a table of all the color images that are available at a given position using the f814w+f606w filters. The `get_image` function reads a single cutout image (as a JPEG color image) and returns a PIL image object.\n", "\n", "See the documentation on [HLA VO services](http://hla.stsci.edu/hla_help.html#services) and the [fitscut image cutout service](http://hla.stsci.edu/fitscutcgi_interface.html) for more information on the web services being used." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "def query_hla(ra,dec,size=0.0,imagetype=\"color\",inst=\"ACS\",format=\"image/jpeg\",\n", " spectral_elt=(\"f814w\",\"f606w\"),autoscale=95.0,asinh=1,\n", " naxis=33):\n", " # convert a list of filters to a comma-separated string\n", " if not isinstance(spectral_elt,str):\n", " spectral_elt = \",\".join(spectral_elt)\n", " siapurl = (\"https://hla.stsci.edu/cgi-bin/hlaSIAP.cgi?\"\n", " \"pos={ra},{dec}&size={size}&imagetype={imagetype}&inst={inst}\"\n", " \"&format={format}&spectral_elt={spectral_elt}\"\n", " \"&autoscale={autoscale}&asinh={asinh}\"\n", " \"&naxis={naxis}\").format(**locals())\n", " votable = Table.read(siapurl,format=\"votable\")\n", " return votable\n", "\n", "def get_image(url):\n", " \n", " \"\"\"Get image from a URL\"\"\"\n", " \n", " r = requests.get(url)\n", " im = Image.open(BytesIO(r.content))\n", " return im\n", "\n", "# display earliest and latest images side-by-side\n", "# top 10 highest PM objects\n", "wsel = wpmh[:10]\n", "nim = len(wsel)\n", "icols = 1 # objects per row\n", "ncols = 2*icols # two images for each object\n", "nrows = (nim+icols-1)//icols\n", "\n", "imsize = 33\n", "xcross = np.array([-1,1,0,0,0])*2 + imsize/2\n", "ycross = np.array([0,0,0,-1,1])*2 + imsize/2\n", "\n", "pylab.rcParams.update({\"font.size\":14})\n", "pylab.figure(1,(12, (12/ncols)*nrows))\n", "for jim, i in enumerate(wsel):\n", " hlatab = query_hla(tab[\"RA\"][i],tab[\"Dec\"][i],naxis=imsize)\n", " # sort by observation date\n", " hlatab = hlatab[np.argsort(hlatab['StartTime'])]\n", " k = 0\n", " im1 = get_image(hlatab['URL'][k])\n", " pylab.subplot(nrows,ncols,2*jim+1)\n", " pylab.imshow(im1,origin=\"upper\")\n", " pylab.plot(xcross,ycross,'g')\n", " pylab.title(hlatab['StartTime'][k].decode('utf-8'))\n", " pylab.ylabel(\"ObjID {}\".format(tab[\"ObjID\"][i]))\n", "\n", " k = -1\n", " im2 = get_image(hlatab['URL'][k])\n", " pylab.subplot(nrows,ncols,2*jim+2)\n", " pylab.imshow(im2,origin=\"upper\")\n", " pylab.plot(xcross,ycross,'g')\n", " pylab.title(hlatab['StartTime'][k].decode('utf-8'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Look at the entire collection of images for the highest PM object" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "i = wpmh[0]\n", "print(tab['ObjID','RA','Dec','bpm','lpm','yr','dT'][i])\n", "imsize = 33\n", "hlatab = query_hla(tab[\"RA\"][i],tab[\"Dec\"][i],naxis=imsize)\n", "# sort by observation date\n", "hlatab = hlatab[np.argsort(hlatab['StartTime'])]\n", "\n", "nim = len(hlatab)\n", "ncols = 8\n", "nrows = (nim+ncols-1)//ncols\n", "\n", "xcross = np.array([-1,1,0,0,0])*2 + imsize/2\n", "ycross = np.array([0,0,0,-1,1])*2 + imsize/2\n", "\n", "pylab.rcParams.update({\"font.size\":11})\n", "pylab.figure(1,(20, (20/ncols)*nrows))\n", "t0 = time.time()\n", "for k in range(nim):\n", " im1 = get_image(hlatab['URL'][k])\n", " pylab.subplot(nrows,ncols,k+1)\n", " pylab.imshow(im1,origin=\"upper\")\n", " pylab.plot(xcross,ycross,'g')\n", " pylab.title(hlatab['StartTime'][k].decode('utf-8'))\n", " if ((k+1) % 10)==0:\n", " print(\"{:.1f} s: finished {} of {}\".format(time.time()-t0,k+1,nim))\n", "pylab.tight_layout()\n", "print(\"{:.1f} s: finished {}\".format(time.time()-t0,nim))" ] }, { "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 }