
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


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:')


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


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))


# 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


bin = 0.01
hrange = (0,6)
bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1

tsplit = 6
dmaglim = 0.05

plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,12))
plt.subplot(211)
wmag = np.where((tab['a_f606w_mad']<dmaglim) & (tab['a_f814w_mad']<dmaglim))[0]
w = wmag[tab['dT'][wmag]<=tsplit]
plt.hist(tab['lpmerr'][w], range=hrange, bins=bincount, label='Longitude Error', 
           histtype='step', linewidth=2)
plt.hist(tab['bpmerr'][w], 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 with dT < {} yrs, dmag < {}'.format(len(w),tsplit,dmaglim))
plt.tight_layout()

plt.subplot(212)
w = wmag[tab['dT'][wmag]>tsplit]
plt.hist(tab['lpmerr'][w], range=hrange, bins=bincount, label='Longitude Error', 
           histtype='step', linewidth=2)
plt.hist(tab['bpmerr'][w], 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 with dT > {} yrs, dmag < {}'.format(len(w),tsplit,dmaglim))
plt.tight_layout()

plt.savefig('{}sweeps_pmerr_loghist2.png'.format(resPath))


bin = 1
hrange = (0,130)
bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1
plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,10))
plt.hist(tab['NumVisits'], range=hrange, bins=bincount, label='Number of visits ', 
           histtype='step', linewidth=2)
plt.xlabel('Number of visits')
plt.ylabel('Number of objects')
plt.autoscale(enable=True, axis='x', tight=True)
plt.ylim(0,200000)
plt.title('{:,} stars in SWEEPS'.format(len(tab)))
plt.tight_layout()
plt.savefig('{}sweeps_numvisits_hist.png'.format(resPath))


bin = 1
hrange = (2000, 2020)
bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1
plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,10))
plt.hist(tab['yr'], range=hrange, bins=bincount, label='year ', histtype='step', linewidth=2)
plt.xlabel('mean detection epoch (year)')
plt.ylabel('Number of objects')
plt.autoscale(enable=True, axis='x', tight=True)
plt.ylim(0,300000)
plt.title('{:,} stars in SWEEPS'.format(len(tab)))
plt.tight_layout()
plt.savefig('{}sweeps_year_hist.png'.format(resPath))


bin = 0.25
hrange = (0, 15)
bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1
plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,10))
plt.hist(tab['dT'], range=hrange, bins=bincount, label='year ', histtype='step', linewidth=2)
plt.xlabel('time span (years)')
plt.ylabel('Number of objects')
plt.autoscale(enable=True, axis='x', tight=True)
plt.yscale('log')
plt.title('{:,} stars in SWEEPS'.format(len(tab)))
plt.tight_layout()
plt.savefig('{}sweeps_year_hist.png'.format(resPath))


bin = 0.025
hrange = (16,28)
bincount = int((hrange[1]-hrange[0])/bin + 0.5) + 1
plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,10))
plt.hist(tab['a_f606w'], range=hrange, bins=bincount, label='F606W', 
           histtype='step',  linewidth=2)
plt.hist(tab['a_f814w'], range=hrange, bins=bincount, label='F814W', 
           histtype='step', linewidth=2)
plt.xlabel('magnitude')
plt.ylabel('Number of stars[in {:0.2} magnitude bins]'.format(bin))
plt.legend(loc='upper right')
plt.autoscale(enable=True, axis='x', tight=True)
plt.title('{:,} stars in SWEEPS'.format(len(tab)))
plt.tight_layout()
plt.savefig('{}sweeps_mag_hist.png'.format(resPath))


bin = 0.001
hrange = (0,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['a_f606w_mad'], range=hrange, bins=bincount, label='F606W', 
           histtype='step',  linewidth=2)
plt.hist(tab['a_f814w_mad'], range=hrange, bins=bincount, label='F814W', 
           histtype='step', linewidth=2)
plt.xlabel('magnitude error (median absolute deviation)')
plt.ylabel('Number of stars[in {:0.2} magnitude bins]'.format(bin))
plt.legend(loc='upper right')
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_magerr_hist.png'.format(resPath))


f606w = tab['a_f606w']
f814w = tab['a_f814w']
RminusI = f606w-f814w

# Calculate the point density
w = np.where((RminusI > -1) & (RminusI < 4))[0]
x = np.array(RminusI[w])
y = np.array(f606w[w])
t0 = time.time()
myPDF,axes = fastKDE.pdf(x,y,numPoints=2**10+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.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')
plt.autoscale(tight=True)
plt.xlabel('A_F606W - A_F814W')
plt.ylabel('A_F606W')
plt.gca().invert_yaxis()
plt.xlim(-1,4)
plt.ylim(27.5,17.5)
plt.colorbar()
plt.text(.93,.93,'{:,} stars in SWEEPS'.format(len(x)),
       horizontalalignment='right',
       transform=plt.gca().transAxes)
plt.savefig("{}sweeps_colormag1.png".format(resPath))

kde took 4.1 sec
Plotting 147230 of 443502 points


# define function
def positions(Obj, jobs=None):
    """
    input parameter Obj is the value of the ObjID 
    optional jobs parameter re-uses casjobs jobs variable
    output plots change in (lon, lat) as a function of time
    overplots proper motion fit
    provides number of objects and magnitude/color information
    """
    if not jobs:
        jobs = mastcasjobs.MastCasJobs(context=HSCContext)

    # execute these as "system" queries so they don't fill up your Casjobs history

    # get the measured positions as a function of time
    query = """SELECT dT, dLon, dLat 
        from AstromSourcePositions where ObjID={}
        order by dT
        """.format(Obj)
    pos = jobs.quick(query, context=HSCContext, task_name="SWEEPS/Microlensing",
                     astropy=True, system=True)
    
    # get the PM fit parameters
    query = """SELECT pmlon, pmlonerr, pmlat, pmlaterr
        from AstromProperMotions where ObjID={}
        """.format(Obj)
    pm = jobs.quick(query, context=HSCContext, task_name="SWEEPS/Microlensing",
                    astropy=True, system=True)
    
    lpm = pm['pmlon'][0]
    bpm = pm['pmlat'][0]
    
    # get the intercept for the proper motion fit referenced to the start time
    # time between mean epoch and zero (ref) epoch (years)

    # get median magnitudes and colors for labeling
    query = """SELECT a_f606w=i1.MagMed, a_f606_m_f814w=i1.MagMed-i2.MagMed
        from AstromSumMagAper2 i1 
        join AstromSumMagAper2 i2 on i1.ObjID=i2.ObjID 
        where i1.ObjID={} and i1.filter='f606w' and i2.filter='f814w' 
        """.format(Obj)
    phot = jobs.quick(query, context=HSCContext, task_name="SWEEPS/Microlensing",
                      astropy=True, system=True)
    f606w = phot['a_f606w'][0]
    f606wmf814w = phot['a_f606_m_f814w'][0]

    x = pos['dT']
    y = pos['dLon']
    plt.rcParams.update({'font.size':10})
    plt.figure(1,(6,3))
    plt.subplot(121)
    plt.scatter(x, y, s=10)
    # xpm = np.linspace(0, max(x), 10)
    xpm = np.array([x.min(),x.max()])
    ypm = lpm*xpm
    plt.plot(xpm, ypm, '-r')
    plt.xlabel('dT (yrs)')
    plt.ylabel('dLon (mas)')
    y = pos['dLat']
    plt.subplot(122)
    plt.scatter(x, y, s=10)
    ypm = bpm*xpm
    plt.plot(xpm, ypm, '-r')
    plt.xlabel('dT (yrs)')
    plt.ylabel('dLat (mas)')
    plt.suptitle("""ObjID {0}
{1} detections,  (lpm, bpm) = ({2:.1f}, {3:.1f}) mas/yr
(f606w, f606w-f814w) = ({4:.1f}, {5:.1f})""".format(Obj, len(x), lpm, bpm, f606w, f606wmf814w),
                  size=10)
    plt.tight_layout(rect=[0, 0.0, 1, 0.88])
    plt.show()
    plt.close()


n = tab['NumVisits']
dev = tab['pmdev']
objid = tab['ObjID']
lpmerr0 = np.array(tab['lpmerr'])
bpmerr0 = np.array(tab['bpmerr'])
wi = np.where( (dev < 1.5) & (n > 90) & (np.sqrt(bpmerr0**2+lpmerr0**2) < 1.0))[0]
print("Plotting {} objects".format(len(wi)))
for o in objid[wi]:
    positions(o, jobs=jobs)

Plotting 21 objects


f606w = tab['a_f606w']
f814w = tab['a_f814w']
RminusI = f606w-f814w

w = np.where((RminusI > -1) & (RminusI < 4))[0]
f606w_mad = tab['a_f606w_mad']
f814w_mad = tab['a_f814w_mad']

t0=time.time()
# Calculate the point density
x1 = np.array(f606w[w])
y1 = np.array(f606w_mad[w])
y1log = np.log(y1)
myPDF1,axes1 = fastKDE.pdf(x1,y1log,numPoints=2**10+1)
finterp = RectBivariateSpline(axes1[1],axes1[0],myPDF1)
z1 = finterp(y1log,x1,grid=False)
# Sort the points by density, so that the densest points are plotted last
idx = z1.argsort()
xs1, ys1, zs1 = x1[idx], y1[idx], z1[idx]

# select a random subset of points in the most crowded regions to speed up plotting
wran = np.where(np.random.random(len(zs1))*zs1<0.05)[0]
print("Plotting {} of {} points".format(len(wran),len(zs1)))
xs1 = xs1[wran]
ys1 = ys1[wran]
zs1 = zs1[wran]

x2 = np.array(f814w[w])
y2 = np.array(f814w_mad[w])
y2log = np.log(y2)
myPDF2,axes2 = fastKDE.pdf(x2,y2log,numPoints=2**10+1)
finterp = RectBivariateSpline(axes2[1],axes2[0],myPDF2)
z2 = finterp(y2log,x2,grid=False)
idx = z2.argsort()
xs2, ys2, zs2 = x2[idx], y2[idx], z2[idx]
print("{:.1f} s: completed kde".format(time.time()-t0))

# select a random subset of points in the most crowded regions to speed up plotting
wran = np.where(np.random.random(len(zs2))*zs2<0.05)[0]
print("Plotting {} of {} points".format(len(wran),len(zs2)))
xs2 = xs2[wran]
ys2 = ys2[wran]
zs2 = zs2[wran]

xr = (18,27)
xx = np.arange(501)*(xr[1]-xr[0])/500.0 + xr[0]
xcut1 = 24.2
xnorm1 = 0.03
xcut2 = 23.0
xnorm2 = 0.03

# only plot a subset of the points to speed things up
qsel = 3
xs1 = xs1[::qsel]
ys1 = ys1[::qsel]
zs1 = zs1[::qsel]
xs2 = xs2[::qsel]
ys2 = ys2[::qsel]
zs2 = zs2[::qsel]

plt.figure(1,(15,8))
plt.subplot(121)
plt.yscale('log')
plt.scatter(xs1,ys1,c=zs1,s=2,edgecolor='none',cmap='plasma')
plt.autoscale(tight=True)
# overplot an error limit that varies with magnitude of the form listed below
plt.plot(xx,xnorm1 * (1. + 10.**(0.4*(xx-xcut1))),linewidth=2.0,
           label='$%.2f (1+10^{0.4(M-%.1f)})$' % (xnorm1,xcut1))
plt.legend(loc='upper left')
plt.xlabel('F606W')
plt.ylabel('F606W_MAD')

plt.subplot(122)
plt.yscale('log')
plt.scatter(xs2,ys2,c=zs2,s=2,edgecolor='none',cmap='plasma')
plt.autoscale(tight=True)
# overplot an error limit that varies with magnitude of the form listed below
plt.plot(xx,xnorm2 * (1. + 10.**(0.4*(xx-xcut2))),linewidth=2.0,
           label='$%.2f (1+10^{0.4(M-%.1f)})$' % (xnorm2,xcut2))
plt.legend(loc='upper left')
plt.xlabel('F814W')
plt.ylabel('F814W_MAD')
plt.tight_layout()

Plotting 274457 of 443502 points
8.5 s: completed kde
Plotting 250807 of 443502 points


def noisecut(tab, factor=1.0):
    """Return boolean array with noise cut in f606w and f814w using model
    factor is normalization factor to use (>1 means allow more noise)
    """
    f606w = tab['a_f606w']
    f814w = tab['a_f814w']
    f606w_mad = tab['a_f606w_mad']
    f814w_mad = tab['a_f814w_mad']
    
    # noise model computed above
    xcut_f606w = 24.2
    xnorm_f606w = 0.03 * factor
    xcut_f814w = 23.0
    xnorm_f814w = 0.03 * factor
    return ((f606w_mad < xnorm_f606w*(1+10.0**(0.4*(f606w-xcut_f606w))))
          & (f814w_mad < xnorm_f814w*(1+10.0**(0.4*(f814w-xcut_f814w))))
           )

# low-noise objects
good = noisecut(tab,factor=1.0)

# Calculate the point density
w = np.where((RminusI > 0) & (RminusI < 3) & good)[0]
x = np.array(RminusI[w])
y = np.array(f606w[w])
t0 = time.time()
myPDF,axes = fastKDE.pdf(x,y,numPoints=2**10+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.075)[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.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')
plt.xlabel('A_F606W - A_F814W')
plt.ylabel('A_F606W')
plt.gca().invert_yaxis()
plt.xlim(-1,4)
plt.ylim(27.5,17.5)
plt.colorbar()
plt.text(.93,.93,'{:,} stars in SWEEPS'.format(len(x)),
       horizontalalignment='right',
       transform=plt.gca().transAxes)
plt.savefig("{}sweeps_colormag2.png".format(resPath))

kde took 3.1 sec
Plotting 129892 of 333525 points


# good defined above
f606w = tab['a_f606w']
f814w = tab['a_f814w']
RminusI = f606w-f814w
w = np.where((RminusI > 0) & (RminusI < 3) & good)[0]
lpm = np.array(tab['lpm'][w])
bpm = np.array(tab['bpm'][w])
x = np.array(RminusI[w])
y = np.array(f606w[w])

nbins = 50
count2d, yedge, xedge = np.histogram2d(y, x, bins=nbins)
lpm_sum = np.histogram2d(y, x, bins=nbins, weights=lpm-lpm.mean())[0]
bpm_sum = np.histogram2d(y, x, bins=nbins, weights=bpm-bpm.mean())[0]
lpm_sumsq = np.histogram2d(y, x, bins=nbins, weights=(lpm-lpm.mean())**2)[0]
bpm_sumsq = np.histogram2d(y, x, bins=nbins, weights=(bpm-bpm.mean())**2)[0]

ccount = count2d.clip(1) 
lpm_mean = lpm_sum/ccount
bpm_mean = bpm_sum/ccount
lpm_rms = np.sqrt(lpm_sumsq/ccount-lpm_mean**2)
bpm_rms = np.sqrt(bpm_sumsq/ccount-bpm_mean**2)
lpm_msigma = lpm_rms/np.sqrt(ccount)
bpm_msigma = bpm_rms/np.sqrt(ccount)

plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,10))
ww = np.where(count2d > 100)
yy, xx = np.mgrid[:nbins,:nbins]
xx = (0.5*(xedge[1:]+xedge[:-1]))[xx]
yy = (0.5*(yedge[1:]+yedge[:-1]))[yy]
# plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')
Q = plt.quiver(xx[ww],yy[ww],lpm_mean[ww],bpm_mean[ww],color='red',width=0.0015)
qlength = 5
plt.quiverkey(Q,0.8,0.97,qlength,'{} mas/yr'.format(qlength),coordinates='axes',labelpos='W')
plt.gca().invert_yaxis()
plt.autoscale(tight=True)
plt.xlim((xedge[0],xedge[-1]))
plt.ylim((yedge[-1],yedge[0]))
plt.xlabel('A_F606W - A_F814W')
plt.ylabel('A_F606W')
plt.savefig('{}sweeps_vecmean.png'.format(resPath))


plt.rcParams.update({'font.size': 16})
plt.figure(1,(20,6))

sub1 = plt.subplot(131)
plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')
plt.xlim((xedge[0],xedge[-1]))
plt.ylim((yedge[0],yedge[-1]))
plt.gca().invert_yaxis()
plt.xlabel('A_F606W - A_F814W')
plt.ylabel('A_F606W')
plt.text(0.97, 0.97,
       'Color-magnitude\n{:,} stars\nin SWEEPS'.format(len(x)),
       horizontalalignment='right', verticalalignment='top',
       transform=plt.gca().transAxes)
plt.colorbar()

sub2 = plt.subplot(132)
plt.xlim((xedge[0],xedge[-1]))
plt.ylim((yedge[0],yedge[-1]))
mask = (lpm_msigma <= 1.0) & (count2d > 10)
im = (lpm_mean+lpm.mean())*mask
im[~mask] = np.nan
vmax = np.nanmax(np.abs(im))
plt.imshow(im,cmap='RdYlGn',aspect="auto",origin="lower",
            extent=(xedge[0],xedge[-1],yedge[0],yedge[-1]))
plt.gca().invert_yaxis()
plt.xlabel('A_F606W - A_F814W')
plt.text(0.97, 0.97,'Mean Longitude PM\n$\sigma(\mathrm{PM}) \leq 1$ and $N > 10$',
       horizontalalignment='right', verticalalignment='top',
       transform=plt.gca().transAxes)
cbar = plt.colorbar()
cbar.ax.set_ylabel('mas/yr', rotation=270)

sub2 = plt.subplot(133)
plt.xlim((xedge[0],xedge[-1]))
plt.ylim((yedge[0],yedge[-1]))
im = lpm_rms*(count2d>10)
plt.imshow(im,cmap='magma',aspect="auto",origin="lower",
            extent=(xedge[0],xedge[-1],yedge[0],yedge[-1]),
            norm=LogNorm(vmin=im[im>0].min(), vmax=im.max()))
plt.gca().invert_yaxis()
plt.xlabel('A_F606W - A_F814W')
plt.text(0.97, 0.97,'RMS in Longitude PM\nN > 10',
       horizontalalignment='right', verticalalignment='top',
       transform=plt.gca().transAxes)
cbar = plt.colorbar()
cbar.ax.set_ylabel('mas/yr', rotation=270)

plt.tight_layout()
plt.savefig("{}sweeps_PMlon.png".format(resPath))


plt.rcParams.update({'font.size': 16})
plt.figure(1,(20,6))

sub1 = plt.subplot(131)
plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')
plt.xlim((xedge[0],xedge[-1]))
plt.ylim((yedge[0],yedge[-1]))
plt.gca().invert_yaxis()
plt.xlabel('A_F606W - A_F814W')
plt.ylabel('A_F606W')
plt.text(0.97,0.97,
       'Color-magnitude\n{:,} stars\nin SWEEPS'.format(len(x)),
       horizontalalignment='right', verticalalignment='top',
       transform=plt.gca().transAxes)
plt.colorbar()

sub2 = plt.subplot(132)
plt.xlim((xedge[0],xedge[-1]))
plt.ylim((yedge[0],yedge[-1]))
mask = (bpm_msigma <= 1.0) & (count2d > 10)
im = (bpm_mean+bpm.mean())*mask
im[~mask] = np.nan
vmax = np.nanmax(np.abs(im))
plt.imshow(im,cmap='RdYlGn',aspect="auto",origin="lower",
            extent=(xedge[0],xedge[-1],yedge[0],yedge[-1]))
plt.gca().invert_yaxis()
plt.xlabel('A_F606W - A_F814W')
plt.text(.97,.97,'Mean Latitude PM\n$\sigma(\mathrm{PM}) \leq 1$ and $N > 10$',
       horizontalalignment='right', verticalalignment='top',
       transform=plt.gca().transAxes)
cbar = plt.colorbar()
cbar.ax.set_ylabel('mas/yr', rotation=270)

sub2 = plt.subplot(133)
plt.xlim((xedge[0],xedge[-1]))
plt.ylim((yedge[0],yedge[-1]))
im = bpm_rms*(count2d>10)
plt.imshow(im,cmap='magma',aspect="auto",origin="lower",
            extent=(xedge[0],xedge[-1],yedge[0],yedge[-1]),
            norm=LogNorm(vmin=im[im>0].min(), vmax=im.max()))
plt.gca().invert_yaxis()
plt.xlabel('A_F606W - A_F814W')
plt.text(.97,.97,'RMS in Latitude PM\nN > 10',
       horizontalalignment='right', verticalalignment='top',
       transform=plt.gca().transAxes)
cbar = plt.colorbar()
cbar.ax.set_ylabel('mas/yr', rotation=270)

plt.tight_layout()
plt.savefig("{}sweeps_PMlat.png".format(resPath))


# locate ridge line
iridgex = np.argmax(myPDF,axis=1)
pdfx = myPDF[np.arange(len(iridgex),dtype=int),iridgex]
# pdfx = myPDF.max(axis=1)
wx = np.where(pdfx > 0.1)[0]
iridgex = iridgex[wx]
# use weighted sum of 2*hw+1 points around peak
hw = 10
pridgex = 0.0
pdenom = 0.0
for k in range(-hw,hw+1):
    wt = myPDF[wx,iridgex+k]
    pridgex = pridgex + k*wt
    pdenom = pdenom + wt
pridgex = iridgex + pridgex/pdenom
ridgex = np.interp(pridgex, np.arange(len(axes[0])), axes[0])

# Fit the data using a polynomial model
x0 = axes[1][wx].min()
x1 = axes[1][wx].max()
p_init = models.Polynomial1D(9)
fit_p = fitting.LinearLSQFitter()
xx = (axes[1][wx]-x0)/(x1-x0)
yoff = 0.65
yy = np.log(ridgex - yoff)
p = fit_p(p_init, xx, yy)

# define useful functions for the ridge line

def ridge_color(f606w, function=p, yoff=yoff, x0=x0, x1=x1):
    """Return R-I position of ridge line as a function of f606w magnitude
    
    function, yoff, x0, x1 are from polynomial fit above
    """
    return yoff + np.exp(p((f606w-x0)/(x1-x0)))

# calculate grid of function values for approximate inversion
rxgrid = axes[1][wx]
rygrid = ridge_color(rxgrid)
color_domain = [rygrid[0],rygrid[-1]]
mag_domain = [axes[1][wx[0]], axes[1][wx[-1]]]
print("color_domain {}".format(color_domain))
print("mag_domain   {}".format(mag_domain))

def ridge_mag(RminusI, xgrid=rxgrid, ygrid=rygrid):
    """Return f606w position of ridge line as a function of R-I color
    
    Uses simple linear interpolation to get approximate value
    """
    f606w = np.interp(RminusI,ygrid,xgrid)
    f606w[(RminusI < ygrid[0]) | (RminusI > ygrid[-1])] = np.nan
    return f606w

color_domain [0.6883374859498668, 2.156051675823446]
mag_domain   [19.356585323810577, 26.480549454689026]


ridgexf = yoff + np.exp(p(xx))

plt.figure(1,(12,6))
plt.subplot(121)
plt.plot(axes[1][wx], ridgex, 'bo')
plt.plot(axes[1][wx], ridgexf, color='red')
plt.ylabel('R-I')
plt.xlabel('R')

# check the derivative plot to see if it stays positive
deriv = np.exp(p(xx))*models.Polynomial1D.horner(xx, 
                            (p.parameters * np.arange(len(p.parameters)))[1:])
plt.subplot(122)
plt.semilogy(axes[1][wx], np.exp(p(xx))*models.Polynomial1D.horner(xx, 
                            (p.parameters * np.arange(len(p.parameters)))[1:]),color='red')
plt.xlabel('R')
plt.ylabel('Fit derivative')
plt.title('Min deriv {:.6f}'.format(deriv.min()))
plt.tight_layout()


plt.rcParams.update({'font.size': 16})
plt.figure(1,(12,10))
plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')
plt.xlim((xedge[0],xedge[-1]))
plt.ylim((yedge[0],yedge[-1]))
plt.gca().invert_yaxis()

# overplot ridge line
plt.plot(ridge_color(axes[1][wx]), axes[1][wx], color='red')
plt.plot(axes[0], ridge_mag(axes[0]), color='green')

plt.xlabel('A_F606W - A_F814W')
plt.ylabel('A_F606W')
plt.text(.97,.97,
       'Color-magnitude\n{:,} stars\nin SWEEPS'.format(len(x)),
       horizontalalignment='right', verticalalignment='top',
       transform=plt.gca().transAxes)
plt.colorbar()
# plt.savefig('{}sweeps_ridgeline.png'.format(resPath))

<matplotlib.colorbar.Colorbar at 0x7f8ef8541e20>


yloc = ridge_mag(x)
x1 = x[np.isfinite(yloc)]
wy = np.where((axes[0] >= x1.min()) & (axes[0] <= x1.max()))[0]
ridgey = ridge_mag(axes[0][wy])

# Weighted histogram
dmagmin = -2.0
dmagmax = 1.0
xmax = axes[0][wy[-1]]
# xmin = axes[0][wy[0]]
xmin = 1.0
wsel = np.where((y-yloc >= dmagmin) & (y-yloc <= dmagmax) & (x >= xmin) & (x <= xmax))[0]

x2 = y[wsel]-yloc[wsel]
y2 = lpm[wsel]-lpm.mean()
hrange = (dmagmin, dmagmax)
hbins = 50
count1d, xedge1d = np.histogram(x2,range=hrange,bins=hbins)
lpm_sum1d = np.histogram(x2,range=hrange,bins=hbins,weights=y2)[0]
lpm_sumsq1d = np.histogram(x2,range=hrange,bins=hbins,weights=y2**2)[0]

ccount1d = count1d.clip(1)
lpm_mean1d = lpm_sum1d/ccount1d
lpm_rms1d = np.sqrt(lpm_sumsq1d/ccount1d-lpm_mean1d**2)
lpm_msigma1d = lpm_rms1d/np.sqrt(ccount1d)
lpm_mean1d += lpm.mean()

x1d = 0.5*(xedge1d[1:]+xedge1d[:-1])

plt.rcParams.update({'font.size': 16})
plt.figure(1,(14,6))

sub1 = plt.subplot(121)
plt.scatter(xs, ys, c=zs, s=2, edgecolor='none', cmap='plasma')
plt.xlim((xedge[0],xedge[-1]))
plt.ylim((yedge[0],yedge[-1]))
plt.gca().invert_yaxis()
plt.xlabel('R - I (A_F606W - A_F814W)')
plt.ylabel('R (A_F606W)')
plt.text(.97,.97,
       'Color-magnitude\n{:,} stars\nin SWEEPS'.format(len(x)),
       horizontalalignment='right', verticalalignment='top',
       transform=plt.gca().transAxes)
xboundary = np.hstack((axes[0][wy],axes[0][wy[::-1]]))
yboundary = np.hstack((ridgey+dmagmax,ridgey[::-1]+dmagmin))
wb = np.where((xboundary >= xmin) & (xboundary <= xmax))
xboundary = xboundary[wb]
yboundary = yboundary[wb]
print(xboundary[0],yboundary[0],xboundary[-1],yboundary[-1])
print(xboundary.shape,yboundary.shape)
xboundary = np.append(xboundary,xboundary[0])
yboundary = np.append(yboundary,yboundary[0])
plt.plot(xboundary, yboundary, color='red')

sub1 = plt.subplot(122)
# don't plot huge error points
wp = np.where(lpm_msigma1d < 1)
plt.errorbar(x1d[wp], lpm_mean1d[wp], xerr=(xedge1d[1]-xedge1d[0])/2.0, 
               yerr=lpm_msigma1d[wp], linestyle='')
plt.autoscale(tight=True)
plt.xlabel('Distance from ridge line [R mag]')
plt.ylabel('Mean Longitude PM [mas/yr]')
plt.gca().invert_xaxis()
plt.text(.03,.97,'{:,} stars in SWEEPS\n1-sigma error bars on mean PM'.format(len(x2)),
       horizontalalignment='left', verticalalignment='top',
       transform=plt.gca().transAxes)
plt.tight_layout()

1.0020693112164736 23.784969024110666 1.0020693112164736 20.784969024110666
(394,) (394,)


w = np.where((RminusI > 0) & (RminusI < 3) & good)[0]

# Calculate the point density
x = np.array(RminusI[w])
y = np.array(f606w[w])
myPDF,axes = fastKDE.pdf(x,y,numPoints=2**10+1)
finterp = RectBivariateSpline(axes[1],axes[0],myPDF)
z = finterp(y,x,grid=False)
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.1)[0]
print("Plotting {} of {} points".format(len(wran),len(zs)))
xs = xs[wran]
ys = ys[wran]
zs = zs[wran]

# locate ridge line in magnitude as a function of color
xloc = ridge_color(y)
ridgex = ridge_color(axes[1][wx])

# locate ridge line in color as function of magnitude
yloc = ridge_mag(x)
x1 = x[np.isfinite(yloc)]
wy = np.where((axes[0] >= x1.min()) & (axes[0] <= x1.max()))[0]
ridgey = ridge_mag(axes[0][wy])

# low-noise objects
print("Selected {:,} low-noise objects".format(len(w)))

# red objects
ylim = yloc - 1.5 - (yloc - 25.0).clip(0)/(1+10.0**(-0.4*(yloc-26.0)))
wred = np.where((y<25) & (y > 19.5) & (y < ylim) & (x-xloc > 0.35))[0]
#wred = np.where((y<25) & (y > 19.5) & ((y-yloc) < -1.5)
#                   & (x-xloc > 0.3))[0]
#                   & (x > 1.1) & (x < 2.5) & (x-xloc > 0.2))[0]
print("Selected {:,} red objects".format(len(wred)))
# main sequence objects
wmain = np.where((y>21) & (y<22.4) & (np.abs(x-xloc) < 0.1))[0]
print("Initially selected {:,} MS objects".format(len(wmain)))
# sort by distance from ridge and select the closest
wmain = wmain[np.argsort(np.abs(x[wmain]-xloc[wmain]))]
wmain = wmain[:len(wred)]
print("Selected {:,} MS objects closest to ridge".format(len(wmain)))

plt.rcParams.update({'font.size': 14})

plt.figure(1,(12,8))
plt.subplot2grid((2,3), (0,0), rowspan=2, colspan=2)
plt.scatter(xs,ys,c=zs,s=2,cmap='plasma',edgecolor='none')
plt.scatter(x[wred],y[wred],c='red',s=2,edgecolor='none')
plt.scatter(x[wmain],y[wmain],c='blue',s=2,edgecolor='none')
plt.xlim(0,3)
plt.ylim(18,27.5)
plt.xlabel('F606W-F814W [mag]')
plt.ylabel('F606W [mag]')
plt.gca().invert_yaxis()
plt.title('{:,} stars in SWEEPS'.format(len(x)))

lrange = (-20,5)
brange = (-12.5,12.5)

# plot MS and red points in random order
wsel = w[np.append(wmain,wred)]
colors = np.array(['blue']*len(wsel))
colors[len(wmain):] = 'red'
irs = np.argsort(np.random.random(len(wsel)))
wsel = wsel[irs]
colors = colors[irs]

plt.subplot2grid((2,3), (0,2), rowspan=1, colspan=1)
plt.scatter(tab['lpm'][w],tab['bpm'][w],c='darkgray',s=0.1)
plt.scatter(tab['lpm'][wsel],tab['bpm'][wsel],c=colors,s=2)
plt.xlim(*lrange)
plt.ylim(*brange)
plt.ylabel('Latitude PM [mas/yr]')

plt.subplot2grid((2,3), (1,2), rowspan=1, colspan=1)
bin = 0.5
bincount = int((lrange[1]-lrange[0])/bin + 0.5) + 1
plt.hist(tab['lpm'][w[wmain]], range=lrange, bins=bincount, label='MS', color='blue', 
           histtype='step')
plt.hist(tab['lpm'][w[wred]], range=lrange, bins=bincount, label='Red', color='red', 
           histtype='step')
plt.xlim(*lrange)
plt.xlabel('Longitude PM [mas/yr]')
plt.ylabel('Number [in {:.2} mas bins]'.format(bin))
plt.legend(loc='upper left')
plt.tight_layout()
plt.savefig('{}sweeps_calamida.png'.format(resPath))

Plotting 156351 of 333525 points
Selected 333,525 low-noise objects
Selected 2,797 red objects
Initially selected 59,605 MS objects
Selected 2,797 MS objects closest to ridge

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

Hubble Source Catalog SWEEPS Proper Motion Notebook¶

June 2019, Steve Lubow and Rick White¶

Instructions:¶

Table of Contents¶

Initialization ¶

Install Python modules¶

Set up your CasJobs account information¶

Other optional configuration¶

Create table in MyDB with selected SWEEPS objects¶

Properties of Full Catalog ¶

Sky Coverage ¶

Proper Motion Histograms ¶

Proper motion histograms for lon and lat¶

Proper motion error cumulative histogram for lon and lat¶

Proper motion error log histogram for lon and lat¶

Proper motion error as a function of dT¶

Proper motion error log histogram for lon and lat¶

Number of Visits Histogram ¶

Time Histograms ¶

Magnitude Histograms ¶

Aper2 magnitude histograms for F606W and F814W¶

Aper2 magnitude error histograms for F606W and F814W¶

Color-Magnitude Diagram ¶

Color-magnitude diagram¶

Detection Positions ¶

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.¶

Good Photometric Objects ¶

Look at photometric error distribution to pick out good photometry objects as a function of magnitude ¶

Define function to apply noise cut and plot color-magnitude diagram with cut¶

Science Applications ¶

Proper Motions of Good Objects ¶

Average proper motion in color-magnitude bins¶

RMS in longitude PM as a function of color/magnitude¶

Mean longitude PM as image¶

Mean latitude PM as image¶

Proper Motions in Bulge and Disk ¶

Fit a smooth function to the main ridgeline of color-magnitude diagram¶

Plot the ridgeline on the CMD¶

Binned distribution of PM(Long) vs magnitude offset from ridge line¶

Reproduce Figure 1 from Calamida et al. 2014 ¶

Hubble Source Catalog SWEEPS Proper Motion Notebook¶

June 2019, Steve Lubow and Rick White¶

Instructions:¶

Table of Contents¶

Initialization ¶

Install Python modules¶

Set up your CasJobs account information¶

Other optional configuration¶

Create table in MyDB with selected SWEEPS objects¶

Properties of Full Catalog ¶

Sky Coverage ¶

Proper Motion Histograms ¶

Proper motion histograms for lon and lat¶

Proper motion error cumulative histogram for lon and lat¶

Proper motion error log histogram for lon and lat¶

Proper motion error as a function of dT¶

Proper motion error log histogram for lon and lat¶

Number of Visits Histogram ¶

Time Histograms ¶

Magnitude Histograms ¶

Aper2 magnitude histograms for F606W and F814W¶

Aper2 magnitude error histograms for F606W and F814W¶

Color-Magnitude Diagram ¶

Color-magnitude diagram¶

Detection Positions ¶

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.¶

Good Photometric Objects ¶

Look at photometric error distribution to pick out good photometry objects as a function of magnitude ¶

Define function to apply noise cut and plot color-magnitude diagram with cut¶

Science Applications ¶

Proper Motions of Good Objects ¶

Average proper motion in color-magnitude bins¶

RMS in longitude PM as a function of color/magnitude¶

Mean longitude PM as image¶

Mean latitude PM as image¶

Proper Motions in Bulge and Disk ¶

Fit a smooth function to the main ridgeline of color-magnitude diagram¶

Plot the ridgeline on the CMD¶

Binned distribution of PM(Long) vs magnitude offset from ridge line¶

Reproduce Figure 1 from Calamida et al. 2014¶

Reproduce Figure 1 from Calamida et al. 2014 ¶