{ "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']tsplit]\n", "plt.hist(tab['lpmerr'][w], range=hrange, bins=bincount, label='Longitude Error', \n", " histtype='step', linewidth=2)\n", "plt.hist(tab['bpmerr'][w], 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 with dT > {} yrs, dmag < {}'.format(len(w),tsplit,dmaglim))\n", "plt.tight_layout()\n", "\n", "plt.savefig('{}sweeps_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", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(12,10))\n", "plt.hist(tab['NumVisits'], range=hrange, bins=bincount, label='Number of visits ', \n", " histtype='step', linewidth=2)\n", "plt.xlabel('Number of visits')\n", "plt.ylabel('Number of objects')\n", "plt.autoscale(enable=True, axis='x', tight=True)\n", "plt.ylim(0,200000)\n", "plt.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "plt.tight_layout()\n", "plt.savefig('{}sweeps_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", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(12,10))\n", "plt.hist(tab['yr'], range=hrange, bins=bincount, label='year ', histtype='step', linewidth=2)\n", "plt.xlabel('mean detection epoch (year)')\n", "plt.ylabel('Number of objects')\n", "plt.autoscale(enable=True, axis='x', tight=True)\n", "plt.ylim(0,300000)\n", "plt.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "plt.tight_layout()\n", "plt.savefig('{}sweeps_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": {}, "outputs": [], "source": [ "bin = 0.25\n", "hrange = (0, 15)\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['dT'], range=hrange, bins=bincount, label='year ', histtype='step', linewidth=2)\n", "plt.xlabel('time span (years)')\n", "plt.ylabel('Number of objects')\n", "plt.autoscale(enable=True, axis='x', tight=True)\n", "plt.yscale('log')\n", "plt.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "plt.tight_layout()\n", "plt.savefig('{}sweeps_year_hist.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Magnitude Histograms " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aper2 magnitude histograms for F606W and F814W" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bin = 0.025\n", "hrange = (16,28)\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['a_f606w'], range=hrange, bins=bincount, label='F606W', \n", " histtype='step', linewidth=2)\n", "plt.hist(tab['a_f814w'], range=hrange, bins=bincount, label='F814W', \n", " histtype='step', linewidth=2)\n", "plt.xlabel('magnitude')\n", "plt.ylabel('Number of stars[in {:0.2} magnitude bins]'.format(bin))\n", "plt.legend(loc='upper right')\n", "plt.autoscale(enable=True, axis='x', tight=True)\n", "plt.title('{:,} stars in SWEEPS'.format(len(tab)))\n", "plt.tight_layout()\n", "plt.savefig('{}sweeps_mag_hist.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aper2 magnitude error histograms for F606W and F814W" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bin = 0.001\n", "hrange = (0,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['a_f606w_mad'], range=hrange, bins=bincount, label='F606W', \n", " histtype='step', linewidth=2)\n", "plt.hist(tab['a_f814w_mad'], range=hrange, bins=bincount, label='F814W', \n", " histtype='step', linewidth=2)\n", "plt.xlabel('magnitude error (median absolute deviation)')\n", "plt.ylabel('Number of stars[in {:0.2} magnitude bins]'.format(bin))\n", "plt.legend(loc='upper right')\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_magerr_hist.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Color-Magnitude Diagram " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Color-magnitude diagram\n", "\n", "Plot the color-magnitude diagram for the ~440k points retrieved from the database. This uses fastkde to compute the kernel density estimate for the crowded plot, which is very fast. See https://pypi.org/project/fastkde/ for instructions -- or just do\n", "
\n", "pip install fastkde" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f606w = tab['a_f606w']\n", "f814w = tab['a_f814w']\n", "RminusI = f606w-f814w\n", "\n", "# Calculate the point density\n", "w = np.where((RminusI > -1) & (RminusI < 4))[0]\n", "x = np.array(RminusI[w])\n", "y = np.array(f606w[w])\n", "t0 = time.time()\n", "myPDF,axes = fastKDE.pdf(x,y,numPoints=2**10+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.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')\n", "plt.autoscale(tight=True)\n", "plt.xlabel('A_F606W - A_F814W')\n", "plt.ylabel('A_F606W')\n", "plt.gca().invert_yaxis()\n", "plt.xlim(-1,4)\n", "plt.ylim(27.5,17.5)\n", "plt.colorbar()\n", "plt.text(.93,.93,'{:,} stars in SWEEPS'.format(len(x)),\n", " horizontalalignment='right',\n", " transform=plt.gca().transAxes)\n", "plt.savefig(\"{}sweeps_colormag1.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, jobs=None):\n", " \"\"\"\n", " input parameter Obj is the value of the ObjID \n", " optional jobs parameter re-uses casjobs jobs variable\n", " output plots change in (lon, lat) as a function of time\n", " overplots proper motion fit\n", " provides number of objects and magnitude/color information\n", " \"\"\"\n", " if not jobs:\n", " jobs = mastcasjobs.MastCasJobs(context=HSCContext)\n", "\n", " # execute these as \"system\" queries so they don't fill up your Casjobs history\n", "\n", " # get the measured positions as a function of time\n", " query = \"\"\"SELECT dT, dLon, dLat \n", " from AstromSourcePositions where ObjID={}\n", " order by dT\n", " \"\"\".format(Obj)\n", " pos = jobs.quick(query, context=HSCContext, task_name=\"SWEEPS/Microlensing\",\n", " astropy=True, system=True)\n", " \n", " # get the PM fit parameters\n", " query = \"\"\"SELECT pmlon, pmlonerr, pmlat, pmlaterr\n", " from AstromProperMotions where ObjID={}\n", " \"\"\".format(Obj)\n", " pm = jobs.quick(query, context=HSCContext, task_name=\"SWEEPS/Microlensing\",\n", " astropy=True, system=True)\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", " # get median magnitudes and colors for labeling\n", " query = \"\"\"SELECT a_f606w=i1.MagMed, a_f606_m_f814w=i1.MagMed-i2.MagMed\n", " from AstromSumMagAper2 i1 \n", " join AstromSumMagAper2 i2 on i1.ObjID=i2.ObjID \n", " where i1.ObjID={} and i1.filter='f606w' and i2.filter='f814w' \n", " \"\"\".format(Obj)\n", " phot = jobs.quick(query, context=HSCContext, task_name=\"SWEEPS/Microlensing\",\n", " astropy=True, system=True)\n", " f606w = phot['a_f606w'][0]\n", " f606wmf814w = phot['a_f606_m_f814w'][0]\n", "\n", " x = pos['dT']\n", " y = pos['dLon']\n", " plt.rcParams.update({'font.size':10})\n", " plt.figure(1,(6,3))\n", " plt.subplot(121)\n", " plt.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", " plt.plot(xpm, ypm, '-r')\n", " plt.xlabel('dT (yrs)')\n", " plt.ylabel('dLon (mas)')\n", " y = pos['dLat']\n", " plt.subplot(122)\n", " plt.scatter(x, y, s=10)\n", " ypm = bpm*xpm\n", " plt.plot(xpm, ypm, '-r')\n", " plt.xlabel('dT (yrs)')\n", " plt.ylabel('dLat (mas)')\n", " plt.suptitle(\"\"\"ObjID {0}\n", "{1} detections, (lpm, bpm) = ({2:.1f}, {3:.1f}) mas/yr\n", "(f606w, f606w-f814w) = ({4:.1f}, {5:.1f})\"\"\".format(Obj, len(x), lpm, bpm, f606w, f606wmf814w),\n", " size=10)\n", " plt.tight_layout(rect=[0, 0.0, 1, 0.88])\n", " plt.show()\n", " plt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot positions of 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, jobs=jobs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Good Photometric Objects " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Look at photometric error distribution to pick out good photometry objects as a function of magnitude \n", "\n", "The photometric error is mainly a function of magnitude. We make a cut slightly above the typical error to exclude objects that have poor photometry. (In the SWEEPS field, that most often is the result of blending and crowding.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f606w = tab['a_f606w']\n", "f814w = tab['a_f814w']\n", "RminusI = f606w-f814w\n", "\n", "w = np.where((RminusI > -1) & (RminusI < 4))[0]\n", "f606w_mad = tab['a_f606w_mad']\n", "f814w_mad = tab['a_f814w_mad']\n", "\n", "t0=time.time()\n", "# Calculate the point density\n", "x1 = np.array(f606w[w])\n", "y1 = np.array(f606w_mad[w])\n", "y1log = np.log(y1)\n", "myPDF1,axes1 = fastKDE.pdf(x1,y1log,numPoints=2**10+1)\n", "finterp = RectBivariateSpline(axes1[1],axes1[0],myPDF1)\n", "z1 = finterp(y1log,x1,grid=False)\n", "# Sort the points by density, so that the densest points are plotted last\n", "idx = z1.argsort()\n", "xs1, ys1, zs1 = x1[idx], y1[idx], z1[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(zs1))*zs1<0.05)[0]\n", "print(\"Plotting {} of {} points\".format(len(wran),len(zs1)))\n", "xs1 = xs1[wran]\n", "ys1 = ys1[wran]\n", "zs1 = zs1[wran]\n", "\n", "x2 = np.array(f814w[w])\n", "y2 = np.array(f814w_mad[w])\n", "y2log = np.log(y2)\n", "myPDF2,axes2 = fastKDE.pdf(x2,y2log,numPoints=2**10+1)\n", "finterp = RectBivariateSpline(axes2[1],axes2[0],myPDF2)\n", "z2 = finterp(y2log,x2,grid=False)\n", "idx = z2.argsort()\n", "xs2, ys2, zs2 = x2[idx], y2[idx], z2[idx]\n", "print(\"{:.1f} s: completed kde\".format(time.time()-t0))\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(zs2))*zs2<0.05)[0]\n", "print(\"Plotting {} of {} points\".format(len(wran),len(zs2)))\n", "xs2 = xs2[wran]\n", "ys2 = ys2[wran]\n", "zs2 = zs2[wran]\n", "\n", "xr = (18,27)\n", "xx = np.arange(501)*(xr[1]-xr[0])/500.0 + xr[0]\n", "xcut1 = 24.2\n", "xnorm1 = 0.03\n", "xcut2 = 23.0\n", "xnorm2 = 0.03\n", "\n", "# only plot a subset of the points to speed things up\n", "qsel = 3\n", "xs1 = xs1[::qsel]\n", "ys1 = ys1[::qsel]\n", "zs1 = zs1[::qsel]\n", "xs2 = xs2[::qsel]\n", "ys2 = ys2[::qsel]\n", "zs2 = zs2[::qsel]\n", "\n", "plt.figure(1,(15,8))\n", "plt.subplot(121)\n", "plt.yscale('log')\n", "plt.scatter(xs1,ys1,c=zs1,s=2,edgecolor='none',cmap='plasma')\n", "plt.autoscale(tight=True)\n", "# overplot an error limit that varies with magnitude of the form listed below\n", "plt.plot(xx,xnorm1 * (1. + 10.**(0.4*(xx-xcut1))),linewidth=2.0,\n", " label='$%.2f (1+10^{0.4(M-%.1f)})$' % (xnorm1,xcut1))\n", "plt.legend(loc='upper left')\n", "plt.xlabel('F606W')\n", "plt.ylabel('F606W_MAD')\n", "\n", "plt.subplot(122)\n", "plt.yscale('log')\n", "plt.scatter(xs2,ys2,c=zs2,s=2,edgecolor='none',cmap='plasma')\n", "plt.autoscale(tight=True)\n", "# overplot an error limit that varies with magnitude of the form listed below\n", "plt.plot(xx,xnorm2 * (1. + 10.**(0.4*(xx-xcut2))),linewidth=2.0,\n", " label='$%.2f (1+10^{0.4(M-%.1f)})$' % (xnorm2,xcut2))\n", "plt.legend(loc='upper left')\n", "plt.xlabel('F814W')\n", "plt.ylabel('F814W_MAD')\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define function to apply noise cut and plot color-magnitude diagram with cut\n", "Note that we reduce the R-I range to 0-3 here because there are very few objects left bluer than R-I = 0 or redder than R-I = 3." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def noisecut(tab, factor=1.0):\n", " \"\"\"Return boolean array with noise cut in f606w and f814w using model\n", " factor is normalization factor to use (>1 means allow more noise)\n", " \"\"\"\n", " f606w = tab['a_f606w']\n", " f814w = tab['a_f814w']\n", " f606w_mad = tab['a_f606w_mad']\n", " f814w_mad = tab['a_f814w_mad']\n", " \n", " # noise model computed above\n", " xcut_f606w = 24.2\n", " xnorm_f606w = 0.03 * factor\n", " xcut_f814w = 23.0\n", " xnorm_f814w = 0.03 * factor\n", " return ((f606w_mad < xnorm_f606w*(1+10.0**(0.4*(f606w-xcut_f606w))))\n", " & (f814w_mad < xnorm_f814w*(1+10.0**(0.4*(f814w-xcut_f814w))))\n", " )\n", "\n", "# low-noise objects\n", "good = noisecut(tab,factor=1.0)\n", "\n", "# Calculate the point density\n", "w = np.where((RminusI > 0) & (RminusI < 3) & good)[0]\n", "x = np.array(RminusI[w])\n", "y = np.array(f606w[w])\n", "t0 = time.time()\n", "myPDF,axes = fastKDE.pdf(x,y,numPoints=2**10+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.075)[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.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')\n", "plt.xlabel('A_F606W - A_F814W')\n", "plt.ylabel('A_F606W')\n", "plt.gca().invert_yaxis()\n", "plt.xlim(-1,4)\n", "plt.ylim(27.5,17.5)\n", "plt.colorbar()\n", "plt.text(.93,.93,'{:,} stars in SWEEPS'.format(len(x)),\n", " horizontalalignment='right',\n", " transform=plt.gca().transAxes)\n", "plt.savefig(\"{}sweeps_colormag2.png\".format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Science Applications " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Proper Motions of Good Objects " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Average proper motion in color-magnitude bins " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# good defined above\n", "f606w = tab['a_f606w']\n", "f814w = tab['a_f814w']\n", "RminusI = f606w-f814w\n", "w = np.where((RminusI > 0) & (RminusI < 3) & good)[0]\n", "lpm = np.array(tab['lpm'][w])\n", "bpm = np.array(tab['bpm'][w])\n", "x = np.array(RminusI[w])\n", "y = np.array(f606w[w])\n", "\n", "nbins = 50\n", "count2d, yedge, xedge = np.histogram2d(y, x, bins=nbins)\n", "lpm_sum = np.histogram2d(y, x, bins=nbins, weights=lpm-lpm.mean())[0]\n", "bpm_sum = np.histogram2d(y, x, bins=nbins, weights=bpm-bpm.mean())[0]\n", "lpm_sumsq = np.histogram2d(y, x, bins=nbins, weights=(lpm-lpm.mean())**2)[0]\n", "bpm_sumsq = np.histogram2d(y, x, bins=nbins, weights=(bpm-bpm.mean())**2)[0]\n", "\n", "ccount = count2d.clip(1) \n", "lpm_mean = lpm_sum/ccount\n", "bpm_mean = bpm_sum/ccount\n", "lpm_rms = np.sqrt(lpm_sumsq/ccount-lpm_mean**2)\n", "bpm_rms = np.sqrt(bpm_sumsq/ccount-bpm_mean**2)\n", "lpm_msigma = lpm_rms/np.sqrt(ccount)\n", "bpm_msigma = bpm_rms/np.sqrt(ccount)\n", "\n", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(12,10))\n", "ww = np.where(count2d > 100)\n", "yy, xx = np.mgrid[:nbins,:nbins]\n", "xx = (0.5*(xedge[1:]+xedge[:-1]))[xx]\n", "yy = (0.5*(yedge[1:]+yedge[:-1]))[yy]\n", "# plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')\n", "Q = plt.quiver(xx[ww],yy[ww],lpm_mean[ww],bpm_mean[ww],color='red',width=0.0015)\n", "qlength = 5\n", "plt.quiverkey(Q,0.8,0.97,qlength,'{} mas/yr'.format(qlength),coordinates='axes',labelpos='W')\n", "plt.gca().invert_yaxis()\n", "plt.autoscale(tight=True)\n", "plt.xlim((xedge[0],xedge[-1]))\n", "plt.ylim((yedge[-1],yedge[0]))\n", "plt.xlabel('A_F606W - A_F814W')\n", "plt.ylabel('A_F606W')\n", "plt.savefig('{}sweeps_vecmean.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## RMS in longitude PM as a function of color/magnitude" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Mean longitude PM as image" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(20,6))\n", "\n", "sub1 = plt.subplot(131)\n", "plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')\n", "plt.xlim((xedge[0],xedge[-1]))\n", "plt.ylim((yedge[0],yedge[-1]))\n", "plt.gca().invert_yaxis()\n", "plt.xlabel('A_F606W - A_F814W')\n", "plt.ylabel('A_F606W')\n", "plt.text(0.97, 0.97,\n", " 'Color-magnitude\\n{:,} stars\\nin SWEEPS'.format(len(x)),\n", " horizontalalignment='right', verticalalignment='top',\n", " transform=plt.gca().transAxes)\n", "plt.colorbar()\n", "\n", "sub2 = plt.subplot(132)\n", "plt.xlim((xedge[0],xedge[-1]))\n", "plt.ylim((yedge[0],yedge[-1]))\n", "mask = (lpm_msigma <= 1.0) & (count2d > 10)\n", "im = (lpm_mean+lpm.mean())*mask\n", "im[~mask] = np.nan\n", "vmax = np.nanmax(np.abs(im))\n", "plt.imshow(im,cmap='RdYlGn',aspect=\"auto\",origin=\"lower\",\n", " extent=(xedge[0],xedge[-1],yedge[0],yedge[-1]))\n", "plt.gca().invert_yaxis()\n", "plt.xlabel('A_F606W - A_F814W')\n", "plt.text(0.97, 0.97,'Mean Longitude PM\\n$\\sigma(\\mathrm{PM}) \\leq 1$ and $N > 10$',\n", " horizontalalignment='right', verticalalignment='top',\n", " transform=plt.gca().transAxes)\n", "cbar = plt.colorbar()\n", "cbar.ax.set_ylabel('mas/yr', rotation=270)\n", "\n", "sub2 = plt.subplot(133)\n", "plt.xlim((xedge[0],xedge[-1]))\n", "plt.ylim((yedge[0],yedge[-1]))\n", "im = lpm_rms*(count2d>10)\n", "plt.imshow(im,cmap='magma',aspect=\"auto\",origin=\"lower\",\n", " extent=(xedge[0],xedge[-1],yedge[0],yedge[-1]),\n", " norm=LogNorm(vmin=im[im>0].min(), vmax=im.max()))\n", "plt.gca().invert_yaxis()\n", "plt.xlabel('A_F606W - A_F814W')\n", "plt.text(0.97, 0.97,'RMS in Longitude PM\\nN > 10',\n", " horizontalalignment='right', verticalalignment='top',\n", " transform=plt.gca().transAxes)\n", "cbar = plt.colorbar()\n", "cbar.ax.set_ylabel('mas/yr', rotation=270)\n", "\n", "plt.tight_layout()\n", "plt.savefig(\"{}sweeps_PMlon.png\".format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mean latitude PM as image" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(20,6))\n", "\n", "sub1 = plt.subplot(131)\n", "plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')\n", "plt.xlim((xedge[0],xedge[-1]))\n", "plt.ylim((yedge[0],yedge[-1]))\n", "plt.gca().invert_yaxis()\n", "plt.xlabel('A_F606W - A_F814W')\n", "plt.ylabel('A_F606W')\n", "plt.text(0.97,0.97,\n", " 'Color-magnitude\\n{:,} stars\\nin SWEEPS'.format(len(x)),\n", " horizontalalignment='right', verticalalignment='top',\n", " transform=plt.gca().transAxes)\n", "plt.colorbar()\n", "\n", "sub2 = plt.subplot(132)\n", "plt.xlim((xedge[0],xedge[-1]))\n", "plt.ylim((yedge[0],yedge[-1]))\n", "mask = (bpm_msigma <= 1.0) & (count2d > 10)\n", "im = (bpm_mean+bpm.mean())*mask\n", "im[~mask] = np.nan\n", "vmax = np.nanmax(np.abs(im))\n", "plt.imshow(im,cmap='RdYlGn',aspect=\"auto\",origin=\"lower\",\n", " extent=(xedge[0],xedge[-1],yedge[0],yedge[-1]))\n", "plt.gca().invert_yaxis()\n", "plt.xlabel('A_F606W - A_F814W')\n", "plt.text(.97,.97,'Mean Latitude PM\\n$\\sigma(\\mathrm{PM}) \\leq 1$ and $N > 10$',\n", " horizontalalignment='right', verticalalignment='top',\n", " transform=plt.gca().transAxes)\n", "cbar = plt.colorbar()\n", "cbar.ax.set_ylabel('mas/yr', rotation=270)\n", "\n", "sub2 = plt.subplot(133)\n", "plt.xlim((xedge[0],xedge[-1]))\n", "plt.ylim((yedge[0],yedge[-1]))\n", "im = bpm_rms*(count2d>10)\n", "plt.imshow(im,cmap='magma',aspect=\"auto\",origin=\"lower\",\n", " extent=(xedge[0],xedge[-1],yedge[0],yedge[-1]),\n", " norm=LogNorm(vmin=im[im>0].min(), vmax=im.max()))\n", "plt.gca().invert_yaxis()\n", "plt.xlabel('A_F606W - A_F814W')\n", "plt.text(.97,.97,'RMS in Latitude PM\\nN > 10',\n", " horizontalalignment='right', verticalalignment='top',\n", " transform=plt.gca().transAxes)\n", "cbar = plt.colorbar()\n", "cbar.ax.set_ylabel('mas/yr', rotation=270)\n", "\n", "plt.tight_layout()\n", "plt.savefig(\"{}sweeps_PMlat.png\".format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Proper Motions in Bulge and Disk " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fit a smooth function to the main ridgeline of color-magnitude diagram\n", "Fit the R-I vs R values, but force the function to increase montonically with R. We use a log transform of the y coordinate to help." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# locate ridge line\n", "iridgex = np.argmax(myPDF,axis=1)\n", "pdfx = myPDF[np.arange(len(iridgex),dtype=int),iridgex]\n", "# pdfx = myPDF.max(axis=1)\n", "wx = np.where(pdfx > 0.1)[0]\n", "iridgex = iridgex[wx]\n", "# use weighted sum of 2*hw+1 points around peak\n", "hw = 10\n", "pridgex = 0.0\n", "pdenom = 0.0\n", "for k in range(-hw,hw+1):\n", " wt = myPDF[wx,iridgex+k]\n", " pridgex = pridgex + k*wt\n", " pdenom = pdenom + wt\n", "pridgex = iridgex + pridgex/pdenom\n", "ridgex = np.interp(pridgex, np.arange(len(axes[0])), axes[0])\n", "\n", "# Fit the data using a polynomial model\n", "x0 = axes[1][wx].min()\n", "x1 = axes[1][wx].max()\n", "p_init = models.Polynomial1D(9)\n", "fit_p = fitting.LinearLSQFitter()\n", "xx = (axes[1][wx]-x0)/(x1-x0)\n", "yoff = 0.65\n", "yy = np.log(ridgex - yoff)\n", "p = fit_p(p_init, xx, yy)\n", "\n", "# define useful functions for the ridge line\n", "\n", "def ridge_color(f606w, function=p, yoff=yoff, x0=x0, x1=x1):\n", " \"\"\"Return R-I position of ridge line as a function of f606w magnitude\n", " \n", " function, yoff, x0, x1 are from polynomial fit above\n", " \"\"\"\n", " return yoff + np.exp(p((f606w-x0)/(x1-x0)))\n", "\n", "# calculate grid of function values for approximate inversion\n", "rxgrid = axes[1][wx]\n", "rygrid = ridge_color(rxgrid)\n", "color_domain = [rygrid[0],rygrid[-1]]\n", "mag_domain = [axes[1][wx[0]], axes[1][wx[-1]]]\n", "print(\"color_domain {}\".format(color_domain))\n", "print(\"mag_domain {}\".format(mag_domain))\n", "\n", "def ridge_mag(RminusI, xgrid=rxgrid, ygrid=rygrid):\n", " \"\"\"Return f606w position of ridge line as a function of R-I color\n", " \n", " Uses simple linear interpolation to get approximate value\n", " \"\"\"\n", " f606w = np.interp(RminusI,ygrid,xgrid)\n", " f606w[(RminusI < ygrid[0]) | (RminusI > ygrid[-1])] = np.nan\n", " return f606w" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the results to check that they look reasonable." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ridgexf = yoff + np.exp(p(xx))\n", "\n", "plt.figure(1,(12,6))\n", "plt.subplot(121)\n", "plt.plot(axes[1][wx], ridgex, 'bo')\n", "plt.plot(axes[1][wx], ridgexf, color='red')\n", "plt.ylabel('R-I')\n", "plt.xlabel('R')\n", "\n", "# check the derivative plot to see if it stays positive\n", "deriv = np.exp(p(xx))*models.Polynomial1D.horner(xx, \n", " (p.parameters * np.arange(len(p.parameters)))[1:])\n", "plt.subplot(122)\n", "plt.semilogy(axes[1][wx], np.exp(p(xx))*models.Polynomial1D.horner(xx, \n", " (p.parameters * np.arange(len(p.parameters)))[1:]),color='red')\n", "plt.xlabel('R')\n", "plt.ylabel('Fit derivative')\n", "plt.title('Min deriv {:.6f}'.format(deriv.min()))\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot the ridgeline on the CMD" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(12,10))\n", "plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')\n", "plt.xlim((xedge[0],xedge[-1]))\n", "plt.ylim((yedge[0],yedge[-1]))\n", "plt.gca().invert_yaxis()\n", "\n", "# overplot ridge line\n", "plt.plot(ridge_color(axes[1][wx]), axes[1][wx], color='red')\n", "plt.plot(axes[0], ridge_mag(axes[0]), color='green')\n", "\n", "plt.xlabel('A_F606W - A_F814W')\n", "plt.ylabel('A_F606W')\n", "plt.text(.97,.97,\n", " 'Color-magnitude\\n{:,} stars\\nin SWEEPS'.format(len(x)),\n", " horizontalalignment='right', verticalalignment='top',\n", " transform=plt.gca().transAxes)\n", "plt.colorbar()\n", "# plt.savefig('{}sweeps_ridgeline.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Binned distribution of PM(Long) vs magnitude offset from ridge line" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "yloc = ridge_mag(x)\n", "x1 = x[np.isfinite(yloc)]\n", "wy = np.where((axes[0] >= x1.min()) & (axes[0] <= x1.max()))[0]\n", "ridgey = ridge_mag(axes[0][wy])\n", "\n", "# Weighted histogram\n", "dmagmin = -2.0\n", "dmagmax = 1.0\n", "xmax = axes[0][wy[-1]]\n", "# xmin = axes[0][wy[0]]\n", "xmin = 1.0\n", "wsel = np.where((y-yloc >= dmagmin) & (y-yloc <= dmagmax) & (x >= xmin) & (x <= xmax))[0]\n", "\n", "x2 = y[wsel]-yloc[wsel]\n", "y2 = lpm[wsel]-lpm.mean()\n", "hrange = (dmagmin, dmagmax)\n", "hbins = 50\n", "count1d, xedge1d = np.histogram(x2,range=hrange,bins=hbins)\n", "lpm_sum1d = np.histogram(x2,range=hrange,bins=hbins,weights=y2)[0]\n", "lpm_sumsq1d = np.histogram(x2,range=hrange,bins=hbins,weights=y2**2)[0]\n", "\n", "ccount1d = count1d.clip(1)\n", "lpm_mean1d = lpm_sum1d/ccount1d\n", "lpm_rms1d = np.sqrt(lpm_sumsq1d/ccount1d-lpm_mean1d**2)\n", "lpm_msigma1d = lpm_rms1d/np.sqrt(ccount1d)\n", "lpm_mean1d += lpm.mean()\n", "\n", "x1d = 0.5*(xedge1d[1:]+xedge1d[:-1])\n", "\n", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(14,6))\n", "\n", "sub1 = plt.subplot(121)\n", "plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')\n", "plt.xlim((xedge[0],xedge[-1]))\n", "plt.ylim((yedge[0],yedge[-1]))\n", "plt.gca().invert_yaxis()\n", "plt.xlabel('R - I (A_F606W - A_F814W)')\n", "plt.ylabel('R (A_F606W)')\n", "plt.text(.97,.97,\n", " 'Color-magnitude\\n{:,} stars\\nin SWEEPS'.format(len(x)),\n", " horizontalalignment='right', verticalalignment='top',\n", " transform=plt.gca().transAxes)\n", "xboundary = np.hstack((axes[0][wy],axes[0][wy[::-1]]))\n", "yboundary = np.hstack((ridgey+dmagmax,ridgey[::-1]+dmagmin))\n", "wb = np.where((xboundary >= xmin) & (xboundary <= xmax))\n", "xboundary = xboundary[wb]\n", "yboundary = yboundary[wb]\n", "print(xboundary[0],yboundary[0],xboundary[-1],yboundary[-1])\n", "print(xboundary.shape,yboundary.shape)\n", "xboundary = np.append(xboundary,xboundary[0])\n", "yboundary = np.append(yboundary,yboundary[0])\n", "plt.plot(xboundary, yboundary, color='red')\n", "\n", "sub1 = plt.subplot(122)\n", "# don't plot huge error points\n", "wp = np.where(lpm_msigma1d < 1)\n", "plt.errorbar(x1d[wp], lpm_mean1d[wp], xerr=(xedge1d[1]-xedge1d[0])/2.0, \n", " yerr=lpm_msigma1d[wp], linestyle='')\n", "plt.autoscale(tight=True)\n", "plt.xlabel('Distance from ridge line [R mag]')\n", "plt.ylabel('Mean Longitude PM [mas/yr]')\n", "plt.gca().invert_xaxis()\n", "plt.text(.03,.97,'{:,} stars in SWEEPS\\n1-sigma error bars on mean PM'.format(len(x2)),\n", " horizontalalignment='left', verticalalignment='top',\n", " transform=plt.gca().transAxes)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reproduce Figure 1 from Calamida et al. 2014\n", "\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "w = np.where((RminusI > 0) & (RminusI < 3) & good)[0]\n", "\n", "# Calculate the point density\n", "x = np.array(RminusI[w])\n", "y = np.array(f606w[w])\n", "myPDF,axes = fastKDE.pdf(x,y,numPoints=2**10+1)\n", "finterp = RectBivariateSpline(axes[1],axes[0],myPDF)\n", "z = finterp(y,x,grid=False)\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.1)[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", "# locate ridge line in magnitude as a function of color\n", "xloc = ridge_color(y)\n", "ridgex = ridge_color(axes[1][wx])\n", "\n", "# locate ridge line in color as function of magnitude\n", "yloc = ridge_mag(x)\n", "x1 = x[np.isfinite(yloc)]\n", "wy = np.where((axes[0] >= x1.min()) & (axes[0] <= x1.max()))[0]\n", "ridgey = ridge_mag(axes[0][wy])\n", "\n", "# low-noise objects\n", "print(\"Selected {:,} low-noise objects\".format(len(w)))\n", "\n", "# red objects\n", "ylim = yloc - 1.5 - (yloc - 25.0).clip(0)/(1+10.0**(-0.4*(yloc-26.0)))\n", "wred = np.where((y<25) & (y > 19.5) & (y < ylim) & (x-xloc > 0.35))[0]\n", "#wred = np.where((y<25) & (y > 19.5) & ((y-yloc) < -1.5)\n", "# & (x-xloc > 0.3))[0]\n", "# & (x > 1.1) & (x < 2.5) & (x-xloc > 0.2))[0]\n", "print(\"Selected {:,} red objects\".format(len(wred)))\n", "# main sequence objects\n", "wmain = np.where((y>21) & (y<22.4) & (np.abs(x-xloc) < 0.1))[0]\n", "print(\"Initially selected {:,} MS objects\".format(len(wmain)))\n", "# sort by distance from ridge and select the closest\n", "wmain = wmain[np.argsort(np.abs(x[wmain]-xloc[wmain]))]\n", "wmain = wmain[:len(wred)]\n", "print(\"Selected {:,} MS objects closest to ridge\".format(len(wmain)))\n", "\n", "plt.rcParams.update({'font.size': 14})\n", "\n", "plt.figure(1,(12,8))\n", "plt.subplot2grid((2,3), (0,0), rowspan=2, colspan=2)\n", "plt.scatter(xs,ys,c=zs,s=2,cmap='plasma',edgecolor='none')\n", "plt.scatter(x[wred],y[wred],c='red',s=2,edgecolor='none')\n", "plt.scatter(x[wmain],y[wmain],c='blue',s=2,edgecolor='none')\n", "plt.xlim(0,3)\n", "plt.ylim(18,27.5)\n", "plt.xlabel('F606W-F814W [mag]')\n", "plt.ylabel('F606W [mag]')\n", "plt.gca().invert_yaxis()\n", "plt.title('{:,} stars in SWEEPS'.format(len(x)))\n", "\n", "lrange = (-20,5)\n", "brange = (-12.5,12.5)\n", "\n", "# plot MS and red points in random order\n", "wsel = w[np.append(wmain,wred)]\n", "colors = np.array(['blue']*len(wsel))\n", "colors[len(wmain):] = 'red'\n", "irs = np.argsort(np.random.random(len(wsel)))\n", "wsel = wsel[irs]\n", "colors = colors[irs]\n", "\n", "plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1)\n", "plt.scatter(tab['lpm'][w],tab['bpm'][w],c='darkgray',s=0.1)\n", "plt.scatter(tab['lpm'][wsel],tab['bpm'][wsel],c=colors,s=2)\n", "plt.xlim(*lrange)\n", "plt.ylim(*brange)\n", "plt.ylabel('Latitude PM [mas/yr]')\n", "\n", "plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1)\n", "bin = 0.5\n", "bincount = int((lrange[1]-lrange[0])/bin + 0.5) + 1\n", "plt.hist(tab['lpm'][w[wmain]], range=lrange, bins=bincount, label='MS', color='blue', \n", " histtype='step')\n", "plt.hist(tab['lpm'][w[wred]], range=lrange, bins=bincount, label='Red', color='red', \n", " histtype='step')\n", "plt.xlim(*lrange)\n", "plt.xlabel('Longitude PM [mas/yr]')\n", "plt.ylabel('Number [in {:.2} mas bins]'.format(bin))\n", "plt.legend(loc='upper left')\n", "plt.tight_layout()\n", "plt.savefig('{}sweeps_calamida.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Mean and median proper motions of bulge stars compared with SgrA* (-6.379 $\\pm$ 0.026, -0.202 $\\pm$ 0.019) mas/yr Reid and Brunthaler (2004)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lpmmain=np.mean(tab['lpm'][w[wmain]])\n", "bpmmain=np.mean(tab['bpm'][w[wmain]])\n", "\n", "norm=1.0/np.sqrt(len(tab['lpm'][w[wmain]]))\n", "print(\"Bulge stars mean PM longitude {:.2f} +- {:.3f} latitude {:.2f} +- {:.3f}\".format(\n", " np.mean(tab['lpm'][w[wmain]]), norm*np.std(tab['lpm'][w[wmain]]),\n", " np.mean(tab['bpm'][w[wmain]]), norm*np.std(tab['bpm'][w[wmain]])\n", " ))\n", "\n", "norm=1.2533/np.sqrt(len(tab['lpm'][w[wmain]]))\n", "print(\"Bulge stars median PM longitude {:.2f} +- {:.3f} latitude {:.2f} +- {:.3f}\".format(\n", " np.median(tab['lpm'][w[wmain]]), norm*np.std(tab['lpm'][w[wmain]]),\n", " np.median(tab['bpm'][w[wmain]]), norm*np.std(tab['bpm'][w[wmain]])\n", " ))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# White Dwarfs " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "w = np.where((RminusI > 0) & (RminusI < 3) & good)[0]\n", "\n", "wwd = np.where((RminusI < 0.7) \n", " & (RminusI > 0) \n", " & (f606w > 22.5) \n", " & good \n", " & (tab['lpm'] < 5) \n", " & (tab['lpm'] > -20))[0]\n", "xwd = np.array(RminusI[wwd])\n", "ywd = np.array(f606w[wwd])\n", "\n", "# Calculate the point density\n", "x = np.array(RminusI[w])\n", "y = np.array(f606w[w])\n", "\n", "myPDF,axes = fastKDE.pdf(x,y,numPoints=2**10+1)\n", "finterp = RectBivariateSpline(axes[1],axes[0],myPDF)\n", "z = finterp(y,x,grid=False)\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.1)[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", "# locate ridge line in magnitude as a function of color\n", "xloc = ridge_color(y)\n", "ridgex = ridge_color(axes[1][wx])\n", "\n", "# locate ridge line in color as function of magnitude\n", "yloc = ridge_mag(x)\n", "x1 = x[np.isfinite(yloc)]\n", "wy = np.where((axes[0] >= x1.min()) & (axes[0] <= x1.max()))[0]\n", "ridgey = ridge_mag(axes[0][wy])\n", "\n", "# low-noise objects\n", "print(\"Selected {:,} low-noise objects\".format(len(w)))\n", "\n", "# red objects\n", "ylim = yloc - 1.5 - (yloc - 25.0).clip(0)/(1+10.0**(-0.4*(yloc-26.0)))\n", "wred = np.where((y<25) & (y > 19.5) & (y < ylim) & (x-xloc > 0.35))[0]\n", "#wred = np.where((y<25) & (y > 19.5) & ((y-yloc) < -1.5)\n", "# & (x-xloc > 0.3))[0]\n", "# & (x > 1.1) & (x < 2.5) & (x-xloc > 0.2))[0]\n", "print(\"Selected {:,} red objects\".format(len(wred)))\n", "# main sequence objects\n", "wmain = np.where((y>21) & (y<22.4) & (np.abs(x-xloc) < 0.1))[0]\n", "print(\"Initially selected {:,} MS objects\".format(len(wmain)))\n", "# sort by distance from ridge and select the closest\n", "wmain = wmain[np.argsort(np.abs(x[wmain]-xloc[wmain]))]\n", "wmain = wmain[:len(wred)]\n", "print(\"Selected {:,} MS objects closest to ridge\".format(len(wmain)))\n", "\n", "plt.rcParams.update({'font.size': 14})\n", "\n", "plt.figure(1,(12,8))\n", "plt.subplot2grid((2,3), (0,0), rowspan=2, colspan=2)\n", "plt.scatter(xs,ys,c=zs,s=2,cmap='plasma',edgecolor='none')\n", "plt.scatter(x[wred],y[wred],c='red',s=2,edgecolor='none')\n", "plt.scatter(x[wmain],y[wmain],c='blue',s=2,edgecolor='none')\n", "plt.scatter(xwd,ywd,c='green',s=10,edgecolor='none')\n", "plt.xlim(0,3)\n", "plt.ylim(18,27.5)\n", "plt.xlabel('F606W-F814W [mag]')\n", "plt.ylabel('F606W [mag]')\n", "plt.gca().invert_yaxis()\n", "plt.title('{:,} stars in SWEEPS'.format(len(x)))\n", "\n", "lrange = (-20,5)\n", "brange = (-12.5,12.5)\n", "\n", "# plot MS and red points in random order\n", "wsel = w[np.append(wmain,wred)]\n", "colors = np.array(['blue']*len(wsel))\n", "colors[len(wmain):] = 'red'\n", "irs = np.argsort(np.random.random(len(wsel)))\n", "wsel = wsel[irs]\n", "colors = colors[irs]\n", "\n", "plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1)\n", "plt.scatter(tab['lpm'][w],tab['bpm'][w],c='darkgray',s=0.1)\n", "plt.scatter(tab['lpm'][wsel],tab['bpm'][wsel],c=colors,s=2)\n", "plt.scatter(tab['lpm'][wwd],tab['bpm'][wwd],c='green',s=2)\n", "plt.xlim(*lrange)\n", "plt.ylim(*brange)\n", "plt.ylabel('Latitude PM [mas/yr]')\n", "\n", "plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1)\n", "bin = 0.5\n", "bincount = int((lrange[1]-lrange[0])/bin + 0.5) + 1\n", "plt.hist(tab['lpm'][wwd], range=lrange, bins=bincount, label='WD', color='green', \n", " histtype='step')\n", "plt.hist(tab['lpm'][w[wmain]], range=lrange, bins=bincount, label='MS', color='blue', \n", " histtype='step')\n", "plt.hist(tab['lpm'][w[wred]], range=lrange, bins=bincount, label='Red', color='red', \n", " histtype='step')\n", "plt.xlim(*lrange)\n", "plt.xlabel('Longitude PM [mas/yr]')\n", "plt.ylabel('Number [in {:.2} mas bins]'.format(bin))\n", "plt.legend(loc='upper left')\n", "plt.tight_layout()\n", "# plt.savefig('{}sweeps_calamida_wd.png'.format(resPath))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## White dwarf mean PM and uncertainity" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "norm=1.0/np.sqrt(len(tab['lpm'][wwd]))\n", "print(\"WD PM mean +- stdev longitude {:.2f} +- {:.3f} latitude {:.2f} +- {:.3f}\".format(\n", " np.mean(tab['lpm'][wwd]), norm*np.std(tab['lpm'][wwd]),\n", " np.mean(tab['bpm'][wwd]), norm*np.std(tab['bpm'][wwd])\n", " ))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## WDs generally closer to bulge (MS) PM distribution" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bin = 0.5\n", "bincount = int((lrange[1]-lrange[0])/bin + 0.5) + 1\n", "plt.hist(tab['lpm'][wwd], range=lrange, bins=bincount, label='WD', color='green', \n", " histtype='step', density=True, cumulative=1)\n", "plt.hist(tab['lpm'][w[wmain]], range=lrange, bins=bincount, label='MS', color='blue', \n", " histtype='step', density=True, cumulative=1)\n", "plt.hist(tab['lpm'][w[wred]], range=lrange, bins=bincount, label='Red', color='red', \n", " histtype='step', density=True, cumulative=1)\n", "plt.xlabel('Longitude PM [mas/yr]')\n", "plt.ylabel('Cumulative fraction'.format(bin))\n", "plt.legend(loc='upper left')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Look for quasar candidates (low PM blue stars) \n", "Note this includes all objects, not just \"good\" objects with low mag noise, because quasars might be variable too." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "wqso1 = np.where((RminusI < 0.5) \n", " & (np.sqrt(tab['bpm']**2+tab['lpm']**2) < 1.0) \n", " & (tab['NumFilters'] > 2))[0]\n", "wqso1 = wqso1[np.argsort(f606w[wqso1])]\n", "\n", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(10,10))\n", "plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')\n", "plt.plot(RminusI[wqso1],f606w[wqso1],'rs',markersize=10,fillstyle='none')\n", "plt.xlim((xedge[0],xedge[-1]))\n", "plt.ylim((yedge[0],yedge[-1]))\n", "plt.gca().invert_yaxis()\n", "plt.xlabel('R - I (A_F606W - A_F814W)')\n", "plt.ylabel('R (A_F606W)')\n", "plt.text(0.97,0.97, \"\"\"{:,} stars in SWEEPS\n", "{} low PM blue stars\n", "with NumFilters$\\\\geq3$\"\"\".format(len(x),len(wqso1)),\n", " horizontalalignment='right', verticalalignment='top',\n", " transform=plt.gca().transAxes)\n", "tab[wqso1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "objid = tab['ObjID']\n", "print(\"Plotting {} objects\".format(len(wqso1)))\n", "for o in objid[wqso1]:\n", " positions(o,jobs=jobs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Try again with different criteria on proper motion fit quality" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wqso2 = np.where((RminusI < 0.5) \n", " & (np.sqrt(tab['bpm']**2+tab['lpm']**2) < 2.0)\n", " & (np.sqrt(tab['bpmerr']**2+tab['lpmerr']**2) < 2.0)\n", " & (tab['pmdev'] < 4.0)\n", " )[0]\n", "wqso2 = wqso2[np.argsort(f606w[wqso2])]\n", "\n", "plt.rcParams.update({'font.size': 16})\n", "plt.figure(1,(10,10))\n", "plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')\n", "plt.plot(RminusI[wqso2],f606w[wqso2],'rs',markersize=10,fillstyle='none')\n", "plt.xlim((xedge[0],xedge[-1]))\n", "plt.ylim((yedge[0],yedge[-1]))\n", "plt.gca().invert_yaxis()\n", "plt.xlabel('R - I (A_F606W - A_F814W)')\n", "plt.ylabel('R (A_F606W)')\n", "plt.text(0.97,0.97, \"\"\"{:,} stars in SWEEPS\n", "{} low PM blue stars\n", "with NumFilters$\\\\geq3$\"\"\".format(len(x),len(wqso2)),\n", " horizontalalignment='right', verticalalignment='top',\n", " transform=plt.gca().transAxes)\n", "tab[wqso2]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "objid = tab['ObjID']\n", "print(\"Plotting {} objects\".format(len(wqso2)))\n", "for o in objid[wqso2]:\n", " positions(o,jobs=jobs)" ] }, { "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 main sequence sample (Galactic center approximately)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f606w = tab['a_f606w']\n", "f814w = tab['a_f814w']\n", "RminusI = f606w-f814w\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-bpmmain)**2+(lpm0-lpmmain)**2)\n", "pmerr0 = np.sqrt(bpmerr0**2+lpmerr0**2)\n", "\n", "# sort samples by decreasing PM\n", "wpml = np.where((pmtot0 > 12) & (pmtot0 < 15) & (pmerr0 < 1.0) & good)[0]\n", "wpml = wpml[np.argsort(-pmtot0[wpml])]\n", "xpml = np.array(RminusI[wpml])\n", "ypml = np.array(f606w[wpml])\n", "\n", "wpmh = np.where((pmtot0 > 15) & (pmerr0 < 1.0) & good)[0]\n", "wpmh = wpmh[np.argsort(-pmtot0[wpmh])]\n", "xpmh = np.array(RminusI[wpmh])\n", "ypmh = np.array(f606w[wpmh])\n", "\n", "# Calculate the point density\n", "w = np.where((RminusI > -1) & (RminusI < 4) & good)[0]\n", "x = np.array(RminusI[w])\n", "y = np.array(f606w[w])\n", "t0 = time.time()\n", "myPDF,axes = fastKDE.pdf(x,y,numPoints=2**10+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.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')\n", "plt.autoscale(tight=True)\n", "plt.xlabel('A_F606W - A_F814W')\n", "plt.ylabel('A_F606W')\n", "plt.gca().invert_yaxis()\n", "plt.xlim(-1,4)\n", "plt.ylim(27.5,17.0)\n", "plt.colorbar()\n", "plt.scatter(xpml, ypml, s=40, c=\"green\", marker =\"^\", label='10 mas/yr < PM < 15 mas/yr')\n", "plt.scatter(xpmh, ypmh, s=80, c=\"red\", marker =\"^\", label='15 mas/yr < PM')\n", "plt.text(0.,0.95,'High proper motions relative to mean of main sequence sample (bulge) \\n \\n {:,} stars in SWEEPS'.format(len(x)),\n", " horizontalalignment='left',\n", " transform=plt.gca().transAxes)\n", "plt.legend(loc='upper right')\n", "# plt.savefig(\"{}sweeps_colormag_highpm.png\".format(resPath))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "print(\"Plotting {} objects\".format(len(wpmh)))\n", "for o in tab[\"ObjID\"][wpmh]:\n", " positions(o,jobs=jobs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Very red high proper motion objects" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "wpmred = np.where((pmtot0 > 12) & (pmerr0 < 1.0) & (RminusI > 2.6) & good)[0]\n", "print(\"Plotting {} objects\".format(len(wpmred)))\n", "for o in tab[\"ObjID\"][wpmred]:\n", " positions(o,jobs=jobs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Very blue high proper motion objects" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "wpmblue = np.where((pmtot0 > 8) & (pmerr0 < 1.0) & (RminusI < 0.5) & good)[0]\n", "print(\"Plotting {} objects\".format(len(wpmblue)))\n", "for o in tab[\"ObjID\"][wpmblue]:\n", " positions(o,jobs=jobs)" ] }, { "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", "# wsel = wpmred\n", "# wsel = wpmblue\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", "plt.rcParams.update({\"font.size\":14})\n", "plt.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", " plt.subplot(nrows,ncols,2*jim+1)\n", " plt.imshow(im1,origin=\"upper\")\n", " plt.plot(xcross,ycross,'g')\n", " plt.title(hlatab['StartTime'][k])\n", " plt.ylabel(\"ObjID {}\".format(tab[\"ObjID\"][i]))\n", "\n", " k = -1\n", " im2 = get_image(hlatab['URL'][k])\n", " plt.subplot(nrows,ncols,2*jim+2)\n", " plt.imshow(im2,origin=\"upper\")\n", " plt.plot(xcross,ycross,'g')\n", " plt.title(hlatab['StartTime'][k])" ] }, { "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','a_f606w','a_f814w','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", "plt.rcParams.update({\"font.size\":11})\n", "plt.figure(1,(20, (20/ncols)*nrows))\n", "t0 = time.time()\n", "for k in range(nim):\n", " im1 = get_image(hlatab['URL'][k])\n", " plt.subplot(nrows,ncols,k+1)\n", " plt.imshow(im1,origin=\"upper\")\n", " plt.plot(xcross,ycross,'g')\n", " plt.title(hlatab['StartTime'][k])\n", " if ((k+1) % 10)==0:\n", " print(\"{:.1f} s: finished {} of {}\".format(time.time()-t0,k+1,nim))\n", "plt.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 (ipykernel)", "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.9.7" } }, "nbformat": 4, "nbformat_minor": 2 }