This notebook shows how to access the new proper motions available for the SWEEPS field in version 3.1 of the Hubble Source Catalog. Data tables in MAST CasJobs are queried from Python using the mastcasjobs module. Additional information is available on the SWEEPS Proper Motions help page.
This notebook is available for download.
Running the notebook from top to bottom takes about 7 minutes (depending on the speed of your computer).
This notebook requires the use of Python 3.
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:
conda install requests conda install pillow pip install git+git://github.com/dfm/casjobs@master pip install git+git://github.com/rlwastro/mastcasjobs@master pip install fastkde
If you already have an older version of the mastcasjobs
module, you may need to update it:
pip install --upgrade git+git://github.com/rlwastro/mastcasjobs@master
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.
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.
This script prompts for your Casjobs user ID and password during initialization if the environment variables are not defined.
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).
resPath="./" # directory where generated plots are saved
HSCContext= "HSCv3"
%matplotlib inline
import astropy, time, sys, os, requests
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from PIL import Image
from io import BytesIO
# check that version of mastcasjobs is new enough
# we are using some features not in version 0.0.1
from pkg_resources import get_distribution
from distutils.version import StrictVersion as V
assert V(get_distribution("mastcasjobs").version) >= V('0.0.2'), """
A newer version of mastcasjobs is required.
Update mastcasjobs to current version using this command:
pip install --upgrade git+git://github.com/rlwastro/mastcasjobs@master
"""
import mastcasjobs
## For handling ordinary astropy Tables
from astropy.table import Table
from astropy.io import fits, ascii
from fastkde import fastKDE
from scipy.interpolate import RectBivariateSpline
from astropy.modeling import models, fitting
# There are a number of relatively unimportant warnings that
# show up, so for now, suppress them:
import warnings
warnings.filterwarnings("ignore")
# Set page width to fill browser for longer output lines
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
# set width for pprint
astropy.conf.max_width = 150
Set up Casjobs environment.
import getpass
if not os.environ.get('CASJOBS_USERID'):
os.environ['CASJOBS_USERID'] = input('Enter Casjobs UserID:')
if not os.environ.get('CASJOBS_PW'):
os.environ['CASJOBS_PW'] = getpass.getpass('Enter Casjobs password:')
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.
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.
DBtable = "SWEEPS"
jobs = mastcasjobs.MastCasJobs(context="MyDB")
try:
print("Retrieving table MyDB.{} (if it exists)".format(DBtable))
tab = jobs.fast_table(DBtable, verbose=True)
except ValueError:
print("Table MyDB.{} not found, running query to create it".format(DBtable))
# drop table if it already exists
jobs.drop_table_if_exists(DBtable)
#get main information
query = """
select a.ObjID, RA=a.raMean, Dec=a.decMean, RAerr=a.raMeanErr, Decerr=a.decMeanErr,
c.NumFilters, c.NumVisits,
a_f606w=i1.MagMed, a_f606w_n=i1.n, a_f606w_mad=i1.MagMAD,
a_f814w=i2.MagMed, a_f814w_n=i2.n, a_f814w_mad=i2.MagMAD,
bpm=a.pmLat, lpm=a.pmLon, bpmerr=a.pmLatErr, lpmerr=a.pmLonErr,
pmdev=sqrt(pmLonDev*pmLonDev+pmLatDev*pmLatDev),
yr=(a.epochMean - 47892)/365.25+1990,
dT=(a.epochEnd-a.epochStart)/365.25,
yrStart=(a.epochStart - 47892)/365.25+1990,
yrEnd=(a.epochEnd - 47892)/365.25+1990
into mydb.{}
from AstromProperMotions a join AstromSumMagAper2 i1 on
i1.ObjID=a.ObjID and i1.n >=10 and i1.filter ='F606W' and i1.detector='ACS/WFC'
join AstromSumMagAper2 i2 on
i2.ObjID=a.ObjID and i2.n >=10 and i2.filter ='F814W' and i2.detector='ACS/WFC'
join AstromSumPropMagAper2Cat c on a.ObjID=c.ObjID
""".format(DBtable)
t0 = time.time()
jobid = jobs.submit(query, task_name="SWEEPS", context=HSCContext)
print("jobid=",jobid)
results = jobs.monitor(jobid)
print("Completed in {:.1f} sec".format(time.time()-t0))
print(results)
# slower version using CasJobs output queue
# tab = jobs.get_table(DBtable, verbose=True)
# fast version using special MAST Casjobs service
tab = jobs.fast_table(DBtable, verbose=True)
tab
Retrieving table MyDB.SWEEPS (if it exists) 26.2 s: Retrieved 157.86MB table MyDB.SWEEPS 35.6 s: Converted to 443932 row table
ObjID | RA | Dec | RAerr | Decerr | NumFilters | NumVisits | a_f606w | a_f606w_n | a_f606w_mad | a_f814w | a_f814w_n | a_f814w_mad | bpm | lpm | bpmerr | lpmerr | pmdev | yr | dT | yrStart | yrEnd |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
int64 | float64 | float64 | float64 | float64 | int32 | int32 | float64 | int32 | float64 | float64 | int32 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 |
4000709002286 | 269.7911379669984 | -29.206156187411423 | 0.6964818624528099 | 0.2730062330800141 | 2 | 47 | 22.127399444580078 | 47 | 0.021600723266601562 | 21.13010025024414 | 47 | 0.0167999267578125 | 2.087558644949346 | -7.738272329430371 | 0.38854582276386007 | 0.22115673368981237 | 2.887154518133692 | 2013.3007902147058 | 11.371914541371764 | 2003.4361796620085 | 2014.8080942033803 |
4000709002287 | 269.7955922590832 | -29.206151631494986 | 0.24020216760386343 | 0.18524811391217816 | 2 | 47 | 21.508499145507812 | 47 | 0.029998779296875 | 20.69930076599121 | 47 | 0.023900985717773438 | -2.8930568503344967 | -0.7898583846555245 | 0.1316584790053578 | 0.12462185695877996 | 1.474676632663785 | 2013.3007902147058 | 11.371914541371764 | 2003.4361796620085 | 2014.8080942033803 |
4000709002288 | 269.81608933789283 | -29.206155196641195 | 0.3040684131020671 | 0.2850407586200256 | 2 | 47 | 21.654399871826172 | 47 | 0.03650093078613281 | 20.85770034790039 | 47 | 0.017101287841796875 | 4.65866649193795 | -3.2098804580343785 | 0.13931172183651183 | 0.20648097604781626 | 1.9570357322713463 | 2013.3007902147058 | 11.371914541371764 | 2003.4361796620085 | 2014.8080942033803 |
4000709002289 | 269.8259694163096 | -29.20615668840751 | 0.3564325426522067 | 0.39542200297333663 | 2 | 47 | 19.79170036315918 | 47 | 0.028200149536132812 | 19.06909942626953 | 47 | 0.019300460815429688 | -0.45662407290928664 | -2.0909050045433832 | 0.15758175952333653 | 0.2763881282194908 | 2.2415238499377685 | 2013.3007902147058 | 11.371914541371764 | 2003.4361796620085 | 2014.8080942033803 |
4000709002290 | 269.83486415728754 | -29.206155266983643 | 0.16299639839198538 | 0.14062839407811836 | 2 | 46 | 20.566649436950684 | 46 | 0.015949249267578125 | 19.847750663757324 | 46 | 0.014801025390625 | 4.459275526783969 | -2.0433632344343886 | 0.17899727943855331 | 0.18503594468835396 | 1.0091970907052557 | 2013.5152382701992 | 3.0067822490178155 | 2011.8013119543625 | 2014.8080942033803 |
4000709002291 | 269.83512411344606 | -29.2061635244798 | 0.18282583105108072 | 0.2093503650681541 | 2 | 46 | 20.17770004272461 | 46 | 0.028949737548828125 | 19.489749908447266 | 46 | 0.03639984130859375 | 4.090870144734149 | -8.059473158394072 | 0.20446351152189052 | 0.26206212529696193 | 1.293050027027329 | 2013.5152382701992 | 3.0067822490178155 | 2011.8013119543625 | 2014.8080942033803 |
4000709002292 | 269.7964913295107 | -29.20618734483311 | 0.30491102397226527 | 0.26784496777851086 | 2 | 47 | 20.83639907836914 | 47 | 0.02239990234375 | 20.088300704956055 | 47 | 0.02230072021484375 | -1.7001866534338244 | -5.963967462759159 | 0.14814161799844547 | 0.1920517681653374 | 1.9671761489287654 | 2013.3007902147058 | 11.371914541371764 | 2003.4361796620085 | 2014.8080942033803 |
4000709002293 | 269.7872745304419 | -29.206257828852337 | 1.351885504360048 | 1.0460614189471378 | 2 | 46 | 25.99174976348877 | 46 | 0.0852499008178711 | 24.204350471496582 | 46 | 0.05684947967529297 | -2.290436263458843 | -9.596838559623627 | 0.779373480816473 | 0.6485354032580087 | 8.518045574208871 | 2013.3206150152016 | 11.371914541371764 | 2003.4361796620085 | 2014.8080942033803 |
4000709002294 | 269.80888716219647 | -29.206189189626002 | 1.4134596118752103 | 1.2518047802862953 | 2 | 47 | 25.465349197387695 | 46 | 0.08104991912841797 | 23.27630043029785 | 47 | 0.032100677490234375 | 4.381302303073221 | -4.701855876394848 | 0.4784428793001857 | 1.0219363775325352 | 6.469654565147953 | 2013.3007902147058 | 11.371914541371764 | 2003.4361796620085 | 2014.8080942033803 |
4000709002295 | 269.8234425365187 | -29.206188473661626 | 0.35597857562160923 | 1.5877104410672558 | 2 | 36 | 24.13195037841797 | 36 | 0.04940032958984375 | 23.0174503326416 | 36 | 0.038100242614746094 | 4.48015296877619 | -8.000985108228505 | 0.5070876244065738 | 0.7286553380951986 | 4.141157872026995 | 2013.174250263767 | 11.298715372835678 | 2003.4361796620085 | 2014.7348950348442 |
4000709002296 | 269.8288659917366 | -29.206188331342627 | 0.3957945604056192 | 0.5687208923740672 | 2 | 47 | 24.883499145507812 | 47 | 0.05340003967285156 | 23.122600555419922 | 47 | 0.02719879150390625 | -4.005807625507199 | -6.558826318246557 | 0.23540629057177406 | 0.34067420575610496 | 3.4689663362344842 | 2013.3007902147058 | 11.371914541371764 | 2003.4361796620085 | 2014.8080942033803 |
4000709002297 | 269.83003885466655 | -29.206190419891996 | 0.31307078829575036 | 0.25454955949348446 | 2 | 46 | 21.571399688720703 | 46 | 0.02180004119873047 | 20.795000076293945 | 46 | 0.027649879455566406 | -0.8725751200872276 | -7.508546040516253 | 0.43513589539477077 | 0.20855507876805024 | 1.8755318392170417 | 2013.5152382701992 | 3.0067822490178155 | 2011.8013119543625 | 2014.8080942033803 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
4000937456444 | 269.76701171186426 | -29.230781486914456 | 1.580152992881023 | 1.6658550649438768 | 2 | 11 | 24.251300811767578 | 11 | 0.023698806762695312 | 23.03260040283203 | 11 | 0.027099609375 | -3.5849521715767922 | -9.355442553537305 | 2.03043488693805 | 1.7417694248809494 | 5.358128308073762 | 2013.8647109613648 | 2.234081043501644 | 2012.5741947321342 | 2014.8082757756358 |
4000941127104 | 269.7541169187084 | -29.192552431738424 | 15.943917890817005 | 5.288230465807725 | 2 | 13 | 25.18440055847168 | 13 | 0.06609916687011719 | 24.42060089111328 | 13 | 0.07610130310058594 | -1769.4297962839264 | 2184.8990523434904 | 1928.3236885219533 | 1914.9658721714695 | 41.439734192938715 | 2004.1529414597892 | 0.017509189482158703 | 2004.14378291216 | 2004.1612921016422 |
4000941133220 | 269.7226622209822 | -29.189773993145547 | 1.781363174670936 | 1.797325809542357 | 2 | 10 | 23.586899757385254 | 10 | 0.013850212097167969 | 22.529250144958496 | 10 | 0.03670024871826172 | -262.66535678489396 | -109.42569661510447 | 381.1594930877206 | 165.39805763629067 | 5.920443736143895 | 2004.153629846252 | 0.01716819423531972 | 2004.1448605605756 | 2004.1620287548108 |
4000946339500 | 269.6899686412092 | -29.27697845186784 | 1.7581744127758168 | 0.9354091074080774 | 2 | 10 | 23.770350456237793 | 10 | 0.06089973449707031 | 22.6496000289917 | 10 | 0.037400245666503906 | -4.433870501660679 | -1.6701167031929542 | 2.2930714454340286 | 1.3766459346536568 | 4.298430412355294 | 2013.9000436508227 | 2.1913492111623043 | 2012.4990161054518 | 2014.6903653166141 |
4000946380565 | 269.71004704275265 | -29.258318619418723 | 1.4337273570795164 | 1.471655230308599 | 2 | 11 | 23.99720001220703 | 11 | 0.03730010986328125 | 22.51849937438965 | 11 | 0.020000457763671875 | 1.6111158133637464 | -6.745663286831152 | 1.9631691651841359 | 1.720280363564143 | 4.409524084435339 | 2013.9002540765803 | 2.3094412747548714 | 2012.4990161054518 | 2014.8084573802066 |
4000946404892 | 269.67530078168915 | -29.247162734242675 | 1.3763906949106726 | 1.6133894601395353 | 2 | 12 | 24.929399490356445 | 12 | 0.02920055389404297 | 23.21535015106201 | 12 | 0.06190013885498047 | -4.240018447223827 | -5.301395921474695 | 2.3240489286645727 | 2.4646517810785467 | 5.034040969645841 | 2013.3725317320211 | 2.0458630494273624 | 2012.389031376682 | 2014.4348944261094 |
4000946417296 | 269.7016328152704 | -29.246070466593743 | 3.130830982636958 | 3.185144112877607 | 2 | 11 | 25.600500106811523 | 11 | 0.23419952392578125 | 24.028099060058594 | 11 | 0.2503986358642578 | 0.3265864932598421 | -6.491394647992156 | 6.396715459697641 | 2.5782652825923047 | 8.829447098663028 | 2013.7731122440448 | 2.01948620036275 | 2012.6708791162512 | 2014.6903653166141 |
4000949430259 | 269.722986295747 | -29.21197112017265 | 1.361963963773738 | 1.158905612896655 | 3 | 15 | 24.502249717712402 | 14 | 0.06614971160888672 | 23.928850173950195 | 14 | 0.056850433349609375 | 4.745787093281551 | -3.4149224348359475 | 0.9875243981622149 | 0.6006686739150457 | 5.081087460494774 | 2004.567449081762 | 6.212572361916957 | 2004.14378291216 | 2010.356355274077 |
4000949692413 | 269.847657491686 | -29.213464437142868 | 2.3678034671958668 | 2.0224597819540673 | 2 | 10 | 24.422550201416016 | 10 | 0.11170005798339844 | 22.581549644470215 | 10 | 0.0209503173828125 | 1.6834992526734194 | -14.490467506156065 | 3.4740642144508795 | 4.492656398011582 | 5.946113188318635 | 2013.955340339603 | 1.5098539308761485 | 2013.2982402725042 | 2014.8080942033803 |
4000949719295 | 269.8230388680129 | -29.20072120181885 | 6.9256239011552925 | 2.49884610751816 | 2 | 10 | 25.977850914001465 | 10 | 0.10779953002929688 | 24.156999588012695 | 10 | 0.03194999694824219 | -0.009307948981069264 | -9.255580848582802 | 1.8581178329802082 | 1.5241924426465314 | 15.473545552827103 | 2012.4700600740625 | 10.768784408135662 | 2003.4361796620085 | 2014.2049640701443 |
4000979902333 | 269.804403240861 | -29.192813265455396 | 2.3859956776447127 | 1.392848515948042 | 2 | 10 | 25.71024990081787 | 10 | 0.06535053253173828 | 24.112099647521973 | 10 | 0.05470085144042969 | 0.7550276872362829 | -4.7745907025067345 | 3.6257686130953477 | 1.461781567688497 | 5.065034698466592 | 2013.6589915119716 | 2.5308665509739208 | 2012.2040284838704 | 2014.7348950348442 |
4000979908546 | 269.8156665822647 | -29.189854516331526 | 0.9143748593616481 | 2.1238192007085037 | 2 | 12 | 25.4466495513916 | 12 | 0.15910053253173828 | 23.336549758911133 | 12 | 0.06319999694824219 | -7.2055855532918205 | -9.683292833608878 | 2.618264686941967 | 3.0724808097663883 | 5.246267309545055 | 2013.8340336110703 | 1.6017431136587867 | 2013.2063510897215 | 2014.8080942033803 |
4000980227788 | 269.7998513044728 | -29.19707771926519 | 1.8494859814888565 | 1.5991491655370784 | 2 | 12 | 24.400450706481934 | 12 | 0.10389995574951172 | 22.93809986114502 | 12 | 0.02064990997314453 | 0.15278595605010709 | -7.172152518535961 | 0.691934557090257 | 0.49373163797911374 | 5.7619744099400085 | 2012.6568721424214 | 11.371914541371764 | 2003.4361796620085 | 2014.8080942033803 |
x = tab['RA']
y = tab['Dec']
plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,10))
plt.scatter(x, y, s=1)
plt.autoscale(tight=True)
plt.xlabel('RA')
plt.ylabel('Dec')
dc=0.01
plt.xlim(min(x)-dc, max(x)+dc)
plt.ylim(min(y)-dc, max(y)+dc)
plt.gca().invert_xaxis()
plt.text(0.5,0.93,'{:,} stars in SWEEPS'.format(len(x)),
horizontalalignment='left',
transform=plt.gca().transAxes)
Text(0.5, 0.93, '443,932 stars in SWEEPS')
bin = 0.2
hrange = (-20,20)
bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1
plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,10))
plt.hist(tab['lpm'], range=hrange, bins=bincount, label='Longitude',
histtype='step', linewidth=2)
plt.hist(tab['bpm'], range=hrange, bins=bincount, label='Latitude',
histtype='step', linewidth=2)
plt.xlabel('Proper motion [mas/yr]')
plt.ylabel('Number [in {:.2} mas bins]'.format(bin))
plt.legend(loc='upper right')
plt.autoscale(enable=True, axis='x', tight=True)
plt.ylim(0,13500)
plt.title('{:,} stars in SWEEPS'.format(len(tab)))
plt.tight_layout()
plt.savefig('{}sweeps_pmerr_hist.png'.format(resPath))
bin = 0.01
hrange = (0,2)
bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1
plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,10))
plt.hist(tab['lpmerr'], range=hrange, bins=bincount, label='Longitude Error',
histtype='step', cumulative=1, linewidth=2)
plt.hist(tab['bpmerr'], range=hrange, bins=bincount, label='Latitude Error',
histtype='step', cumulative=1, linewidth=2)
plt.xlabel('Proper motion error [mas/yr]')
plt.ylabel('Cumulative number [in {:0.2} mas bins]'.format(bin))
plt.legend(loc='upper right')
plt.autoscale(enable=True, axis='x', tight=True)
plt.ylim(0,500000)
plt.title('{:,} stars in SWEEPS'.format(len(tab)))
plt.tight_layout()
plt.savefig('{}sweeps_pmerr_cumhist.png'.format(resPath))
bin = 0.01
hrange = (0,6)
bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1
plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,10))
plt.hist(tab['lpmerr'], range=hrange, bins=bincount, label='Longitude Error',
histtype='step', linewidth=2)
plt.hist(tab['bpmerr'], range=hrange, bins=bincount, label='Latitude Error',
histtype='step', linewidth=2)
plt.xlabel('Proper motion error [mas/yr]')
plt.ylabel('Number [in {:0.2} mas bins]'.format(bin))
plt.legend(loc='upper right')
plt.yscale('log')
plt.autoscale(enable=True, axis='x', tight=True)
plt.ylim(0,15000)
plt.title('{:,} stars in SWEEPS'.format(len(tab)))
plt.tight_layout()
plt.savefig('{}sweeps_pmerr_loghist.png'.format(resPath))
Exclude objects with dT near zero, and to improve the plotting add a bit of random noise to spread out the quanitized time values.
# restrict to sources with dT > 1 year
dtmin = 1.0
w = np.where(tab['dT']>dtmin)[0]
if ('rw' not in locals()) or len(rw) != len(w):
rw = np.random.random(len(w))
x = np.array(tab['dT'][w]) + 0.5*(rw-0.5)
y = np.log(np.array(tab['lpmerr'][w]))
# Calculate the point density
t0 = time.time()
myPDF,axes = fastKDE.pdf(x,y,numPoints=2**9+1)
print("kde took {:.1f} sec".format(time.time()-t0))
# interpolate to get z values at points
finterp = RectBivariateSpline(axes[1],axes[0],myPDF)
z = finterp(y,x,grid=False)
# Sort the points by density, so that the densest points are plotted last
idx = z.argsort()
xs, ys, zs = x[idx], y[idx], z[idx]
# select a random subset of points in the most crowded regions to speed up plotting
wran = np.where(np.random.random(len(zs))*zs<0.05)[0]
print("Plotting {} of {} points".format(len(wran),len(zs)))
xs = xs[wran]
ys = ys[wran]
zs = zs[wran]
plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,10))
plt.yscale('log')
plt.scatter(xs, np.exp(ys), c=zs, s=2, edgecolor='none', cmap='plasma',
label='Longitude PM error')
plt.autoscale(tight=True, axis='y')
plt.xlim(0.0, max(x)*1.05)
plt.xlabel('Date range [yr]')
plt.ylabel('Proper motion error [mas/yr]')
plt.legend(loc='best')
plt.title('{:,} stars in SWEEPS'.format(len(tab)))
plt.colorbar()
plt.tight_layout()
# plt.savefig('{}sweeps_pmerr_vs_dt.png'.format(resPath))
kde took 3.6 sec Plotting 170401 of 442682 points