{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Hubble Source Catalog SWEEPS Proper Motion Notebook\n", "### June 2019, Steve Lubow and Rick White\n", "\n", "This notebook shows how to access 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). Data tables in [MAST CasJobs](https://mastweb.stsci.edu/hcasjobs) are queried from Python using the [mastcasjobs](https://github.com/rlwastro/mastcasjobs) module. 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.ipynb)." ] }, { "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 7 minutes (depending on the speed of your computer).\n", "\n", "\n", "# Table of Contents\n", "* [Intialization](#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", " * [Magnitude Distributions](#maghist) \n", " * [Color Magnitude Diagram](#cmdall)\n", " * [Detection Positions](#detpos)\n", "* [Good Photometric Objects](#goodphot)\n", "* [Science Applications](#sciap)\n", " * [Proper Motions on the CMD](#goodpm)\n", " * [Proper Motions in Bulge Versus Disk](#bulgedisk)\n", " * [White Dwarfs](#wd)\n", " * [QSO Candidates](#qsocan)\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 git+git://github.com/dfm/casjobs@master\n", "pip install git+git://github.com/rlwastro/mastcasjobs@master\n", "pip install fastkde\n", "\n", "\n", "If you already have an older version of the `mastcasjobs` module, you may need to update it:\n", "\n", "
\n", "pip install --upgrade git+git://github.com/rlwastro/mastcasjobs@master\n", "\n", "\n", "### Set up your CasJobs account information\n", "\n", "You must have a MAST Casjobs account (see https://mastweb.stsci.edu/hcasjobs to create one). Note that MAST Casjobs accounts are independent of SDSS Casjobs accounts.\n", "\n", "For easy startup, you can optionally set the environment variables `CASJOBS_USERID` and/or `CASJOBS_PW` with your Casjobs account information. The Casjobs user ID and password are what you enter when logging into Casjobs.\n", "\n", "This script prompts for your Casjobs user ID and password during initialization if the environment variables are not defined.\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", "HSCContext= \"HSCv3\"\n", "\n", "%matplotlib inline\n", "import astropy, time, sys, os, requests\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from matplotlib.colors import LogNorm\n", "\n", "from PIL import Image\n", "from io import BytesIO\n", "\n", "# check that version of mastcasjobs is new enough\n", "# we are using some features not in version 0.0.1\n", "from pkg_resources import get_distribution\n", "from distutils.version import StrictVersion as V\n", "assert V(get_distribution(\"mastcasjobs\").version) >= V('0.0.2'), \"\"\"\n", "A newer version of mastcasjobs is required.\n", "Update mastcasjobs to current version using this command:\n", "pip install --upgrade git+git://github.com/rlwastro/mastcasjobs@master\n", "\"\"\"\n", "\n", "import mastcasjobs\n", "\n", "## For handling ordinary astropy Tables\n", "from astropy.table import Table\n", "from astropy.io import fits, ascii\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": [ "Set up Casjobs environment." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import getpass\n", "if not os.environ.get('CASJOBS_USERID'):\n", " os.environ['CASJOBS_USERID'] = input('Enter Casjobs UserID:')\n", "if not os.environ.get('CASJOBS_PW'):\n", " os.environ['CASJOBS_PW'] = getpass.getpass('Enter Casjobs password:')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create table in MyDB with selected SWEEPS objects\n", "\n", "Note that the query restricts the sample to matches with at least 10 detections in each of F606W and F814W. This can be modified depending on the science goals.\n", "\n", "This uses an existing MyDB.SWEEPS table if it already exists in your CasJobs account. If you want to change the query, either change the name of the output table or drop the table to force it to be recreated. Usually the query completes in about a minute, but if the database server is heavily loaded then it can take much longer." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "DBtable = \"SWEEPS\"\n", "jobs = mastcasjobs.MastCasJobs(context=\"MyDB\")\n", "\n", "try:\n", " print(\"Retrieving table MyDB.{} (if it exists)\".format(DBtable))\n", " tab = jobs.fast_table(DBtable, verbose=True)\n", "except ValueError:\n", " print(\"Table MyDB.{} not found, running query to create it\".format(DBtable))\n", "\n", " # drop table if it already exists\n", " jobs.drop_table_if_exists(DBtable)\n", "\n", " #get main information\n", " query = \"\"\"\n", " select a.ObjID, RA=a.raMean, Dec=a.decMean, RAerr=a.raMeanErr, Decerr=a.decMeanErr,\n", " c.NumFilters, c.NumVisits,\n", " a_f606w=i1.MagMed, a_f606w_n=i1.n, a_f606w_mad=i1.MagMAD,\n", " a_f814w=i2.MagMed, a_f814w_n=i2.n, a_f814w_mad=i2.MagMAD,\n", " bpm=a.pmLat, lpm=a.pmLon, bpmerr=a.pmLatErr, lpmerr=a.pmLonErr,\n", " pmdev=sqrt(pmLonDev*pmLonDev+pmLatDev*pmLatDev),\n", " yr=(a.epochMean - 47892)/365.25+1990, \n", " dT=(a.epochEnd-a.epochStart)/365.25,\n", " yrStart=(a.epochStart - 47892)/365.25+1990,\n", " yrEnd=(a.epochEnd - 47892)/365.25+1990\n", " into mydb.{}\n", " from AstromProperMotions a join AstromSumMagAper2 i1 on \n", " i1.ObjID=a.ObjID and i1.n >=10 and i1.filter ='F606W' and i1.detector='ACS/WFC'\n", " join AstromSumMagAper2 i2 on \n", " i2.ObjID=a.ObjID and i2.n >=10 and i2.filter ='F814W' and i2.detector='ACS/WFC'\n", " join AstromSumPropMagAper2Cat c on a.ObjID=c.ObjID\n", " \"\"\".format(DBtable)\n", "\n", " t0 = time.time()\n", " jobid = jobs.submit(query, task_name=\"SWEEPS\", context=HSCContext)\n", " print(\"jobid=\",jobid)\n", " results = jobs.monitor(jobid)\n", " print(\"Completed in {:.1f} sec\".format(time.time()-t0))\n", " print(results)\n", "\n", " # slower version using CasJobs output queue\n", " # tab = jobs.get_table(DBtable, verbose=True)\n", " \n", " # fast version using special MAST Casjobs service\n", " tab = jobs.fast_table(DBtable, verbose=True)\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", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(12,10))\n", "plt.scatter(x, y, s=1)\n", "plt.autoscale(tight=True)\n", "plt.xlabel('RA')\n", "plt.ylabel('Dec')\n", "dc=0.01\n", "plt.xlim(min(x)-dc, max(x)+dc)\n", "plt.ylim(min(y)-dc, max(y)+dc)\n", "plt.gca().invert_xaxis()\n", "plt.text(0.5,0.93,'{:,} stars in SWEEPS'.format(len(x)),\n", " horizontalalignment='left',\n", " transform=plt.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", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(12,10))\n", "plt.hist(tab['lpm'], range=hrange, bins=bincount, label='Longitude', \n", " histtype='step', linewidth=2)\n", "plt.hist(tab['bpm'], range=hrange, bins=bincount, label='Latitude', \n", " histtype='step', linewidth=2)\n", "plt.xlabel('Proper motion [mas/yr]')\n", "plt.ylabel('Number [in {:.2} mas bins]'.format(bin))\n", "plt.legend(loc='upper right')\n", "plt.autoscale(enable=True, axis='x', tight=True)\n", "plt.ylim(0,13500)\n", "plt.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "plt.tight_layout()\n", "plt.savefig('{}sweeps_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", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(12,10))\n", "plt.hist(tab['lpmerr'], range=hrange, bins=bincount, label='Longitude Error', \n", " histtype='step', cumulative=1, linewidth=2)\n", "plt.hist(tab['bpmerr'], range=hrange, bins=bincount, label='Latitude Error', \n", " histtype='step', cumulative=1, linewidth=2)\n", "plt.xlabel('Proper motion error [mas/yr]')\n", "plt.ylabel('Cumulative number [in {:0.2} mas bins]'.format(bin))\n", "plt.legend(loc='upper right')\n", "plt.autoscale(enable=True, axis='x', tight=True)\n", "plt.ylim(0,500000)\n", "plt.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "plt.tight_layout()\n", "plt.savefig('{}sweeps_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", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(12,10))\n", "plt.hist(tab['lpmerr'], range=hrange, bins=bincount, label='Longitude Error', \n", " histtype='step', linewidth=2)\n", "plt.hist(tab['bpmerr'], range=hrange, bins=bincount, label='Latitude Error', \n", " histtype='step', linewidth=2)\n", "plt.xlabel('Proper motion error [mas/yr]')\n", "plt.ylabel('Number [in {:0.2} mas bins]'.format(bin))\n", "plt.legend(loc='upper right')\n", "plt.yscale('log')\n", "plt.autoscale(enable=True, axis='x', tight=True)\n", "plt.ylim(0,15000)\n", "plt.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "plt.tight_layout()\n", "plt.savefig('{}sweeps_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": { "scrolled": false }, "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", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(12,10))\n", "plt.yscale('log')\n", "plt.scatter(xs, np.exp(ys), c=zs, s=2, edgecolor='none', cmap='plasma', \n", " label='Longitude PM error')\n", "plt.autoscale(tight=True, axis='y')\n", "plt.xlim(0.0, max(x)*1.05)\n", "plt.xlabel('Date range [yr]')\n", "plt.ylabel('Proper motion error [mas/yr]')\n", "plt.legend(loc='best')\n", "plt.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "plt.colorbar()\n", "plt.tight_layout()\n", "# plt.savefig('{}sweeps_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": {}, "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", "dmaglim = 0.05\n", "\n", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(12,12))\n", "plt.subplot(211)\n", "wmag = np.where((tab['a_f606w_mad']