Compare photometry in the new and old Stripe82 catalogs

to Gaia DR2 photometry and derive corrections for systematics

Link to helper tools

Link to data reading

Link to data analysis


In [24]:
%matplotlib inline
from astropy.table import Table
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.table import hstack
import matplotlib.pyplot as plt 
import numpy as np
from astroML.plotting import hist
# for astroML installation see https://www.astroml.org/user_guide/installation.html

Helper Tools


In [230]:
### selection tools and numerical analysis 

# robust standard deviation
def sigG(arr):
    return 0.741*(np.quantile(arr, 0.75)-np.quantile(arr, 0.25))

def checkNobs(d, band):
    str1 = band + '_Nobs_old'
    arr1 = d[str1]
    print(band, 'band, OLD:')
    printStats(arr1)
    str1 = band + '_Nobs_new'
    arr1 = d[str1]
    print('        NEW:')
    printStats(arr1)
    print('      DIFF:')
    str2 = 'd' + band
    arr2 = d[str2]
    printStats(arr2)
    return
  
def printStats(arr):
    print('           ', np.min(arr), np.mean(arr), np.median(arr), np.max(arr), np.size(arr)) 
    return 

# given vectors x and y, fit medians in bins from xMin to xMax, with Nbin steps,
# and return xBin, medianBin, medianErrBin 
def fitMedians(x, y, xMin, xMax, Nbin, verbose=1): 

    # first generate bins
    xEdge = np.linspace(xMin, xMax, (Nbin+1)) 
    xBin = np.linspace(0, 1, Nbin)
    nPts = 0*np.linspace(0, 1, Nbin)
    medianBin = 0*np.linspace(0, 1, Nbin)
    sigGbin = -1+0*np.linspace(0, 1, Nbin) 
    for i in range(0, Nbin): 
        xBin[i] = 0.5*(xEdge[i]+xEdge[i+1]) 
        yAux = y[(x>xEdge[i])&(x<=xEdge[i+1])]
        if (yAux.size > 0):
            nPts[i] = yAux.size
            medianBin[i] = np.median(yAux)
            # robust estimate of standard deviation: 0.741*(q75-q25)
            sigmaG = 0.741*(np.percentile(yAux,75)-np.percentile(yAux,25))
            # uncertainty of the median: sqrt(pi/2)*st.dev/sqrt(N)
            sigGbin[i] = np.sqrt(np.pi/2)*sigmaG/np.sqrt(nPts[i])
        else:
            nPts[i] = 0
            medianBin[i] = 0
            sigGbin[i] = 0
        
    if (verbose):
        print('median:', np.median(medianBin[nPts>0]), 'std.dev:', np.std(medianBin[nPts>0]))

    return xBin, nPts, medianBin, sigGbin

In [23]:
## for Gaia-related analysis

from scipy import stats
from scipy import optimize

# this function computes polynomial models given some data x
# and parameters theta
def polynomial_fit(theta, x):
    """Polynomial model of degree (len(theta) - 1)"""
    return sum(t * x ** n for (n, t) in enumerate(theta))

# compute the data log-likelihood given a model
def logL(data, theta, model=polynomial_fit):
    """Gaussian log-likelihood of the model at theta"""
    x, y, sigma_y = data
    y_fit = model(theta, x)
    return sum(stats.norm.logpdf(*args)
               for args in zip(y, y_fit, sigma_y))

# a direct optimization approach is used to get best model 
# parameters (which minimize -logL)
def best_theta(data, degree, model=polynomial_fit):
    theta_0 = (degree + 1) * [0]
    neg_logL = lambda theta: -logL(data, theta, model)
    return optimize.fmin_bfgs(neg_logL, theta_0, disp=False)

 
# this function computes a linear combination of 4 functions
# given parameters theta
def linear_fit(coeffs, x, w, y, z):
    ffit = coeffs[0]*x + coeffs[1]*w + coeffs[2]*y + coeffs[3]*z 
    return ffit

# compute the data log-likelihood given a model
def logL2(dataL, coeffs, model=linear_fit):
    """Gaussian log-likelihood of the model at theta"""
    x, w, y, z, f, sigma_f = dataL
    f_fit = model(coeffs, x, w, y, z)
    return sum(stats.norm.logpdf(*args)
               for args in zip(f, f_fit, sigma_f))

# a direct optimization approach is used to get best model 
# parameters (which minimize -logL)
def best_lintheta(dataL, degree=4, model=linear_fit):
    coeffs_0 = degree * [0]
    neg_logL = lambda coeffs: -logL2(dataL, coeffs, model)
    return optimize.fmin_bfgs(neg_logL, coeffs_0, disp=False)

In [297]:
### plots 

# quick plot 
def qp(d, Xstr, Xmin, Xmax, Ystr, Ymin, Ymax):
    ax = plt.axes()
    ax.set_xlabel(Xstr)
    ax.set_ylabel(Ystr)
    ax.scatter(d[Xstr], d[Ystr], s=0.01, c='blue')  
    ax.set_xlim(Xmin, Xmax)
    ax.set_ylim(Ymin, Ymax)
    plt.show()
    return

# quick plot - compare three subsamples
def qp3(d1, d2, d3, Xstr, Xmin, Xmax, Ystr, Ymin, Ymax):
    ax = plt.axes()
    ax.set_xlabel(Xstr)
    ax.set_ylabel(Ystr)
    ax.scatter(d1[Xstr], d1[Ystr], s=0.01, c='green') 
    ax.scatter(d2[Xstr], d2[Ystr], s=0.01, c='red') 
    ax.scatter(d3[Xstr], d3[Ystr], s=0.01, c='blue') 
    ax.set_xlim(Xmin, Xmax)
    ax.set_ylim(Ymin, Ymax)
    plt.show()
    return

# quick plot - binned median
def qpBM(d, Xstr, Xmin, Xmax, Ystr, Ymin, Ymax, nBin, Nsigma=3, offset=0.01):
         
    print('medianAll:', np.median(d[Ystr]), 'std.dev.All:', sigG(d[Ystr]))
    print('N=', np.size(d[Ystr]), 'min=', np.min(d[Ystr]), 'max=', np.max(d[Ystr]))

    ax = plt.axes()
    ax.scatter(d[Xstr], d[Ystr], s=0.01, c='black') 
    # binning
    xBinM, nPtsM, medianBinM, sigGbinM = fitMedians(d[Xstr], d[Ystr], Xmin, Xmax, nBin, 1)
    # plotting
    ax.scatter(xBinM, medianBinM, s=30.0, c='black', alpha=0.8)
    ax.scatter(xBinM, medianBinM, s=15.0, c='yellow', alpha=0.3)
    #
    TwoSigP = medianBinM + Nsigma*sigGbinM
    TwoSigM = medianBinM - Nsigma*sigGbinM 
    ax.plot(xBinM, TwoSigP, c='yellow')
    ax.plot(xBinM, TwoSigM, c='yellow')
    #
    rmsBin = np.sqrt(nPtsM) / np.sqrt(np.pi/2) * sigGbinM
    rmsP = medianBinM + rmsBin
    rmsM = medianBinM - rmsBin
    ax.plot(xBinM, rmsP, c='cyan')
    ax.plot(xBinM, rmsM, c='cyan')
    # 
    xL = np.linspace(-100,100)
    ax.plot(xL, 0*xL, c='red')
    ax.plot(xL, 0*xL+offset, '--', c='red')
    ax.plot(xL, 0*xL-offset, '--', c='red')
    # 
    ax.set_xlabel(Xstr)
    ax.set_ylabel(Ystr)
    ax.set_xlim(Xmin, Xmax)
    ax.set_ylim(Ymin, Ymax)
    plt.show()
    return
 
def qphist(arr, xMin, xMax, xLabel, verbose = False):
    ax = plt.axes()
    hist(arr, bins='knuth', ax=ax, histtype='stepfilled', ec='k', fc='#AAAAAA')
    ax.set_xlabel(xLabel)
    ax.set_ylabel('n')
    ax.set_xlim(xMin, xMax)
    plt.show()
    if (verbose):
        print('Min, max: ', np.min(arr),np.max(arr)) 
        print('Mean, median: ', np.mean(arr),np.median(arr)) 
        print('sigG, st.dev.: ', sigG(arr),np.std(arr)) 
    return 

def qpH0(arr, xMin, xMax, xLabel, nBins=0, verbose = False):
    ax = plt.axes()
    if (nBins>0):
        hist, bins = np.histogram(arr, bins=nBins)
        center = (bins[:-1]+bins[1:])/2
        ax.plot(center, hist, drawstyle='steps', c='blue')   
    else:
        plt.hist(arr, bins='auto', histtype='stepfilled', ec='k', fc='red') 
 
    ax.set_xlabel(xLabel)
    ax.set_ylabel('n')
    ax.set_xlim(xMin, xMax)
    ax.plot([-1000, 1000], [0, 0], '--k')
    plt.show()
    if (verbose):
        print('Min, max: ', np.min(arr),np.max(arr)) 
        print('Mean, median: ', np.mean(arr),np.median(arr)) 
        print('sigG, st.dev.: ', sigG(arr),np.std(arr))  
    return 

def qpHdm(d, band, dmMax, xMin, xMax, nBins=50, verbose=False):
    str = 'd' + band
    dm = 1000*d[str]
    dmOK = dm[np.abs(dm)<(1000*dmMax)]
    xLabel = str + ' (milimag)'
    qpH0(dmOK, 1000*xMin, 1000*xMax, xLabel, nBins, verbose)
    return np.mean(dmOK), np.median(dmOK), sigG(dmOK)


def qp2hist(arr1, arr2, xMin, xMax, xLabel, nBins=0, verbose = False):
    ax = plt.axes()
    if (nBins>0):
        hist, bins = np.histogram(arr1, bins=nBins)
        center = (bins[:-1]+bins[1:])/2
        ax.plot(center, hist, drawstyle='steps', c='red')   
        hist2, bins2 = np.histogram(arr2, bins=nBins)
        center2 = (bins2[:-1]+bins2[1:])/2
        ax.plot(center2, hist2, drawstyle='steps', c='blue')   
    else:
        plt.hist(arr1, bins='auto', histtype='stepfilled', ec='k', fc='yellow')
        plt.hist(arr2, bins='auto', histtype='stepfilled', ec='red', fc='blue')
 
    ax.set_xlabel(xLabel)
    ax.set_ylabel('n')
    ax.set_xlim(xMin, xMax)
    plt.show()
    if (verbose):
        print('Min: ', np.min(arr1),np.min(arr2)) 
        print('Median: ', np.median(arr1),np.median(arr2)) 
        print('sigG: ', sigG(arr1),sigG(arr2)) 
        print('Max: ', np.max(arr1),np.max(arr2)) 
    return

In [ ]:
## plots involving Gaia data

Define paths and catalogs


In [4]:
ZIdataDir = "/Users/ivezic/Work/Science/CalibrationV2/Data"
# the original SDSS catalog from 2007
sdssOldCat = ZIdataDir + "/" + "stripe82calibStars_v2.6.dat"
# Karun's new catalog from 2020
sdssNewCat = ZIdataDir + "/" + "NEW_stripe82calibStars_v0.dat" 
# Gaia DR2 
GaiaDR2Cat = ZIdataDir + "/" + "Stripe82_GaiaDR2.dat"
# for testing, 1% of the sample 
GaiaDR2Cat1perc = ZIdataDir + "/" + "Stripe82_GaiaDR2_1percent.dat"

In [5]:
# both new and old files use identical data structure
colnamesSDSS = ['calib_fla', 'ra', 'dec', 'raRMS', 'decRMS', 'nEpochs', 'AR_val', 
                'u_Nobs', 'u_mMed', 'u_mMean', 'u_mErr', 'u_rms_scatt', 'u_chi2',
                'g_Nobs', 'g_mMed', 'g_mMean', 'g_mErr', 'g_rms_scatt', 'g_chi2',
                'r_Nobs', 'r_mMed', 'r_mMean', 'r_mErr', 'r_rms_scatt', 'r_chi2',
                'i_Nobs', 'i_mMed', 'i_mMean', 'i_mErr', 'i_rms_scatt', 'i_chi2',
                'z_Nobs', 'z_mMed', 'z_mMean', 'z_mErr', 'z_rms_scatt', 'z_chi2']

In [6]:
%%time
# old
sdssOld = Table.read(sdssOldCat, format='ascii', names=colnamesSDSS) 
np.size(sdssOld)


CPU times: user 23.5 s, sys: 5.33 s, total: 28.8 s
Wall time: 28.9 s
Out[6]:
1006849

In [7]:
%%time
# new 
sdssNewRaw = Table.read(sdssNewCat, format='csv', names=colnamesSDSS)
# rejects objects with bad Declination, sample size from 1007136 to 1005470 
sdssNew = sdssNewRaw[sdssNewRaw['dec']<90]
np.size(sdssNew)


CPU times: user 8.6 s, sys: 2.73 s, total: 11.3 s
Wall time: 9.78 s
Out[7]:
1005470

Simple positional match using ra/dec


In [8]:
sdssOld_coords = SkyCoord(ra = sdssOld['ra']*u.degree, dec= sdssOld['dec']*u.degree) 
sdssNew_coords = SkyCoord(ra = sdssNew['ra']*u.degree, dec= sdssNew['dec']*u.degree) 
# this is matching sdssNew to sdssOld, so that indices are into sdssNew catalog
# makes sense in this case since the sdssOld catalog is (a little bit) bigger 
# than sdssNew (1006849 vs 1005470)
idx, d2d, d3d = sdssNew_coords.match_to_catalog_sky(sdssOld_coords)

In [11]:
# object separation is an object with units, 
# I add that as a column so that one can 
# select based on separation to the nearest matching object
new_old = hstack([sdssNew, sdssOld[idx]], table_names = ['new', 'old'])
new_old['sep_2d_arcsec'] = d2d.arcsec
# good matches between the old and new catalogs
MAX_DISTANCE_ARCSEC = 0.5
sdss = new_old[(new_old['sep_2d_arcsec'] < MAX_DISTANCE_ARCSEC)]
print(np.size(sdss))


1003204

now match to Gaia DR2...


In [16]:
colnamesGaia = ['ra', 'dec', 'nObs', 'Gmag', 'flux', 'fluxErr', 'pmra', 'pmdec']
if (1):
    gaia = Table.read(GaiaDR2Cat, format='ascii', names=colnamesGaia)
else: 
    gaia = Table.read(GaiaDR2Cat1perc, format='ascii', names=colnamesGaia)
gaia['raG'] = gaia['ra']
gaia['decG'] = gaia['dec']

In [67]:
sdss_coords = SkyCoord(ra = sdss['ra_old']*u.degree, dec= sdss['dec_old']*u.degree) 
gaia_coords = SkyCoord(ra = gaia['raG']*u.degree, dec= gaia['decG']*u.degree) 

# this is matching gaia to sdss, so that indices are into sdss catalog
# makes sense in this case since the sdss catalog is bigger than gaia
idxG, d2dG, d3dG = gaia_coords.match_to_catalog_sky(sdss_coords)  

# object separation is an object with units, 
# I add that as a column so that one can 
# select based on separation to the nearest matching object
gaia_sdss = hstack([gaia, sdss[idxG]], table_names = ['gaia', 'sdss'])
gaia_sdss['sepSG_2d_arcsec'] = d2dG.arcsec

In [78]:
### code for generating new quantities, such as dra, ddec, colors, differences in mags, etc
def derivedColumns(matches):
    matches['dra'] = (matches['ra_new']-matches['ra_old'])*3600
    matches['ddec'] = (matches['dec_new']-matches['dec_old'])*3600
    matches['ra'] = matches['ra_old']
    ra = matches['ra'] 
    matches['raW'] = np.where(ra > 180, ra-360, ra) 
    matches['dec'] = matches['dec_old']
    matches['u'] = matches['u_mMed_old']
    matches['g'] = matches['g_mMed_old']
    matches['r'] = matches['r_mMed_old']
    matches['i'] = matches['i_mMed_old']
    matches['z'] = matches['z_mMed_old']
    matches['ug'] = matches['u_mMed_old'] - matches['g_mMed_old']
    matches['gr'] = matches['g_mMed_old'] - matches['r_mMed_old']
    matches['ri'] = matches['r_mMed_old'] - matches['i_mMed_old']
    matches['gi'] = matches['g_mMed_old'] - matches['i_mMed_old']
    matches['du'] = matches['u_mMed_old'] - matches['u_mMed_new']
    matches['dg'] = matches['g_mMed_old'] - matches['g_mMed_new']
    matches['dr'] = matches['r_mMed_old'] - matches['r_mMed_new']
    matches['di'] = matches['i_mMed_old'] - matches['i_mMed_new']
    matches['dz'] = matches['z_mMed_old'] - matches['z_mMed_new']
    # Gaia 
    matches['draGold'] = -3600*(matches['ra_old'] - matches['raG']) 
    matches['draGnew'] = -3600*(matches['ra_new'] - matches['raG']) 
    matches['ddecGold'] = -3600*(matches['dec_old'] - matches['decG']) 
    matches['ddecGnew'] = -3600*(matches['dec_new'] - matches['decG']) 
    # photometric
    matches['gGr_old'] = matches['Gmag'] - matches['r_mMed_old']
    matches['gGr_new'] = matches['Gmag'] - matches['r_mMed_new']
    return

In [79]:
derivedColumns(gaia_sdss)

Select good matches and compare both catalogs to Gaia DR2


In [295]:
## first limit astrometric distance and 
## require at least 4 epochs as in the old catalog
MAX_DISTANCE_ARCSEC = 0.5
m1 = gaia_sdss[(gaia_sdss['sepSG_2d_arcsec'] < MAX_DISTANCE_ARCSEC)]
a1 = m1['g_Nobs_new']
a2 = m1['r_Nobs_new']
a3 = m1['i_Nobs_new']
mOK = m1[(a1>3)&(a2>3)&(a3>3)]
print(len(new_old))
print(len(m1))
print(len(mOK))


1005470
892987
892196

In [81]:
## additional constraints for the u and z band:
## require at least 4 epochs both in the old and new catalogs 
mOKu = mOK[(mOK['u_Nobs_old']>3)&(mOK['u_Nobs_new']>3)]
mOKz = mOK[(mOK['z_Nobs_old']>3)&(mOK['u_Nobs_new']>3)]
print(len(mOK))
print(len(mOKu))
print(len(mOKz))


892196
499584
891601

In [83]:
qp(mOK, 'draGold', -0.7, 0.7, 'pmra', -50, 50)



In [84]:
qp(mOK, 'ddecGold', -0.7, 0.7, 'pmdec', -50, 50)



In [85]:
qp(mOK, 'draGnew', -0.7, 0.7, 'pmra', -50, 50)



In [86]:
qp(mOK, 'ddecGnew', -0.7, 0.7, 'pmdec', -50, 50)



In [87]:
qpBM(mOK, 'raW', -51.5, 60, 'ddecGnew', -0.6, 0.6, 120)


median: -0.06282995037513173 std.dev: 0.0148788419801274

In [88]:
qpBM(mOK, 'decG', -1.3, 1.3, 'ddecGnew', -0.6, 0.6, 120)


median: -0.06621597938480939 std.dev: 0.005257749440549402

In [89]:
qpBM(mOK, 'raW', -51.5, 60, 'ddecGold', -0.6, 0.6, 120)


median: -0.07346659294063806 std.dev: 0.015496055371269446

In [94]:
qpBM(mOK, 'raW', -51.5, 60, 'draGnew', -0.6, 0.6, 120)


medianAll: 0.011811212948487082 std.dev.All: 0.08298736520212664
N= 892196 min= -0.6685039625608624 max= -0.6685039625608624
median: 0.034776713221162936 std.dev: 0.028315364402572804

In [95]:
qpBM(mOK, 'raW', -51.5, 60, 'draGold', -0.6, 0.6, 120)


medianAll: 0.04640065991594611 std.dev.All: 0.14529865150637988
N= 892196 min= -0.4992510178453813 max= -0.4992510178453813
median: 0.05769705885541043 std.dev: 0.02619740837997715

In [294]:
# doGaiaAll(mOK)
def doGaiaAll(d, Cstr):
    # Cstr = 'gGr_old' or 'gGr_new'  
    gi = d['gi']
    Gr = d[Cstr]
    Gmag = d['Gmag']
    xBin, nPts, medianBin, sigGbin = fitMedians(gi, Gr, -0.7, 4.0, 47, 0)
    data = np.array([xBin, medianBin, sigGbin])
    Ndata = xBin.size
    # get best-fit parameters  
    thetaCloc = best_theta(data,7)
    # generate best fit lines on a fine grid 
    xfit = np.linspace(-1.1, 4.3, 1000)
    yfit = polynomial_fit(thetaCloc, xfit) 

    d['gGrFit'] = polynomial_fit(thetaCloc, gi)
    d['dgGr'] = d[Cstr] - d['gGrFit']  
    Dc = d[(d['gi']>0.4)&(d['gi']<3.0)]
    
    DcB = Dc[(Dc['Gmag']>14.5)&(Dc['Gmag']<20.0)]
    DcB['GrResid'] = DcB['dgGr'] - np.median(DcB['dgGr'])
    DcBok = DcB[np.abs(DcB['dgGr'])<0.1]
    
    qpBM(DcBok, 'dec', -1.3, 1.3, 'GrResid', -0.03, 0.03, 260) 
    qpBM(DcBok, 'raW', -51.5, 60, 'GrResid', -0.03, 0.03, 112) 
    
    return thetaC, DcBok

In [298]:
thetaC, DcBok_old = doGaiaAll(mOK, 'gGr_old')


medianAll: 6.945092618722629e-05 std.dev.All: 0.02004531973622611
N= 522222 min= -0.09426535645927858 max= 0.10572221954619065
median: 0.0008714020577174544 std.dev: 0.007166684347508925
medianAll: 6.945092618722629e-05 std.dev.All: 0.02004531973622611
N= 522222 min= -0.09426535645927858 max= 0.10572221954619065
median: -0.00035407516604597095 std.dev: 0.003458964256715476

In [299]:
thetaC, DcBok_new = doGaiaAll(mOK, 'gGr_new')


medianAll: -2.4249865089840533e-05 std.dev.All: 0.019276353963555082
N= 521594 min= -0.09528747200368677 max= 0.10470693719856021
median: 0.0003502347122703881 std.dev: 0.00632767845062075
medianAll: -2.4249865089840533e-05 std.dev.All: 0.019276353963555082
N= 521594 min= -0.09528747200368677 max= 0.10470693719856021
median: -0.0005412109918604639 std.dev: 0.003418639846584066

In [103]:
xBin, nPts, medianBin, sigGbin = fitMedians(mOK['gi'], mOK['gGr_new'], -0.7, 4.0, 47, 0)
xBin2, nPts2, medianBin2, sigGbin2 = fitMedians(mOK['gi'], mOK['gGr_old'], -0.7, 4.0, 47, 0)

In [127]:
data = np.array([xBin, medianBin, sigGbin])
Ndata = xBin.size
# get best-fit parameters for linear, quadratic and cubic models
theta1 = best_theta(data,3)
theta2 = best_theta(data,5)
theta3 = best_theta(data,7)
# generate best fit lines on a fine grid 
xfit = np.linspace(-1.1, 4.3, 1000)
yfit1 = polynomial_fit(theta1, xfit)
yfit2 = polynomial_fit(theta2, xfit)
yfit3 = polynomial_fit(theta3, xfit)

In [128]:
# Plot a (gaiaG - r)  vs (g-i)  for photometric transformation
# plot
fig,ax = plt.subplots(1,1,figsize=(12,8))
ax.scatter(mOK['gi'], mOK['gGr_old'], s=0.01, c='red') 
ax.scatter(mOK['gi'], mOK['gGr_new'], s=0.01, c='blue') 
# medians
ax.scatter(xBin, medianBin, s=30.0, c='yellow', alpha=0.5)
ax.scatter(xBin2, medianBin2, s=30.0, c='cyan', alpha=0.5)
ax.set_xlim(-1,4)
ax.set_ylim(-1.5,0.5)
ax.set_xlabel('SDSS(g-i)')
ax.set_ylabel('Gaia G - SDSS r')
# ax.errorbar(x, y, sigma_y, fmt='ok', ecolor='gray')
ax.plot(xfit, yfit1, label='best P3 model')
ax.plot(xfit, yfit2, label='best P5 model')
ax.plot(xfit, yfit3, label='best P7 model')
ax.legend(loc='best', fontsize=14)
plt.show()



In [136]:
theta3


Out[136]:
array([-0.01295894,  0.1512374 , -0.12672826, -0.06557411,  0.23939858,
       -0.17380899,  0.04648134, -0.00423948])

In [133]:
## use 7-th degree polynomial and 0.4 < g-i < 3.0 for recalibration
thetaC = theta3
mOK['gGrFit'] = polynomial_fit(thetaC, mOK['gi'])
mOK['dgGr'] = mOK['gGr_new'] - mOK['gGrFit']
mOKc = mOK[(mOK['gi']>0.4)&(mOK['gi']<3.0)]

In [134]:
xBinM, nPtsM, medianBinM, sigGbinM = fitMedians(mOK['gi'], mOK['dgGr'], -0.7, 4.0, 47, 0)

In [137]:
fig,ax = plt.subplots(1,1,figsize=(12,8))
ax.scatter(mOK['gi'], mOK['dgGr'], s=0.01, c='red')
ax.scatter(mOKc['gi'], mOKc['dgGr'], s=0.01, c='blue')
# medians
ax.scatter(xBinM, medianBinM, s=30.0, c='black', alpha=0.8)
ax.scatter(xBinM, medianBinM, s=20.0, c='yellow', alpha=0.3)
TwoSigP = medianBinM + 2*sigGbinM
TwoSigM = medianBinM - 2*sigGbinM 
ax.plot(xBinM, TwoSigP, c='yellow')
ax.plot(xBinM, TwoSigM, c='yellow')
rmsBin = np.sqrt(nPtsM) / np.sqrt(np.pi/2) * sigGbinM
rmsP = medianBinM + rmsBin
rmsM = medianBinM - rmsBin
ax.plot(xBinM, rmsP, c='cyan')
ax.plot(xBinM, rmsM, c='cyan')
ax.set_xlim(-1,4.2)
ax.set_ylim(-0.14,0.14)
ax.set_xlabel('g-i')
ax.set_ylabel('residuals for (Gaia G - SDSS r)')
xL = np.linspace(-10,10)
ax.plot(xL, 0*xL+0.00, c='red')
ax.plot(xL, 0*xL+0.01, c='red')
ax.plot(xL, 0*xL-0.01, c='red')
ax.plot(0*xL+0.4, xL, c='yellow')
ax.plot(0*xL+3.0, xL, c='yellow')


Out[137]:
[<matplotlib.lines.Line2D at 0x6420de7d0>]

In [138]:
xBin3, nPts3, medianBin3, sigGbin3 = fitMedians(mOKc['Gmag'], mOKc['dgGr'], 14.0, 21.0, 70, 1)


median: -0.006547697202531096 std.dev: 0.010185923312915577

In [140]:
# Plot a (gaiaG - r)  vs (g-i)  for photometric transformation
# plot
fig,ax = plt.subplots(1,1,figsize=(12,8))
ax.scatter(mOK['Gmag'], mOK['dgGr'] , s=0.01, c='red') 
ax.scatter(mOKc['Gmag'], mOKc['dgGr'] , s=0.01, c='blue')  
# medians
ax.scatter(xBin3, medianBin3, s=30.0, c='yellow', alpha=0.5)
ax.set_xlim(13.5,21.5)
ax.set_ylim(-0.2,0.2)
ax.set_xlabel('Gaia G')
ax.set_ylabel('residuals for (Gaia G - SDSS r)')
# 
offset = 0.01
xL = np.linspace(-100,100)
ax.plot(xL, 0*xL, c='red')
ax.plot(xL, 0*xL+offset, c='red')
ax.plot(xL, 0*xL-offset, c='red')
plt.show()



In [169]:
xBin3x, nPts3x, medianBin3x, sigGbin3x = fitMedians(mOKc['Gmag'], mOKc['dgGr'], 14.0, 21.0, 140, 1)


median: -0.006669137553808449 std.dev: 0.0102064345051115

In [172]:
# Plot a (gaiaG - r)  vs (g-i)  for photometric transformation
# plot
fig,ax = plt.subplots(1,1,figsize=(12,8))
ax.scatter(mOKc['Gmag'], mOKc['dgGr']+0.0065 , s=0.01, c='blue')  
# medians
ax.scatter(xBin3x, medianBin3x+0.0065, s=30.0, c='black', alpha=0.95)
ax.set_xlim(15.0,16.5)
ax.set_ylim(-0.05,0.05)
ax.set_xlabel('Gaia G')
ax.set_ylabel('residuals for (Gaia G - SDSS r)')
# 
offset = 0.01
xL = np.linspace(-100,100)
ax.plot(xL, 0*xL, c='red')
ax.plot(xL, 0*xL+offset, c='red')
ax.plot(xL, 0*xL-offset, c='red')
#ax.plot(xL, 0*xL-0.0065, c='blue')
#ax.plot(xL, 0*xL-0.009, c='blue')
ax.plot(xL, 0*xL-0.0025, '--', c='red')


plt.show()



In [ ]:


In [146]:
## for zero point calibration, take 14.5 < G < 20.0 
mOKcB = mOKc[(mOKc['Gmag']>14.5)&(mOKc['Gmag']<20.0)]
mOKcB['GrResid'] = mOKcB['dgGr'] - np.median(mOKcB['dgGr'])
mOKcBok = mOKcB[np.abs(mOKcB['dgGr'])<0.1]
print(np.size(mOKcB), np.size(mOKcBok))


539281 521594

In [150]:
sigG(mOKcBok['GrResid'])


Out[150]:
0.019276353963555082

In [156]:
qpBM(mOKcBok, 'dec', -1.3, 1.3, 'GrResid', -0.03, 0.03, 260)


medianAll: -2.4249865089840533e-05 std.dev.All: 0.019276353963555082
N= 521594 min= -0.09528747200368677 max= -0.09528747200368677
median: 0.0003502347122703881 std.dev: 0.00632767845062075

In [174]:
qpBM(mOKcBok, 'raW', -51.5, 60, 'GrResid', -0.03, 0.03, 112)


medianAll: -2.4249865089840533e-05 std.dev.All: 0.019276353963555082
N= 521594 min= -0.09528747200368677 max= -0.09528747200368677
median: -0.0005412109918604639 std.dev: 0.003418639846584066

Recalibrate R.A. residuals


In [312]:
RAbin, RAnPts, RAmedianBin, RAsigGbin = fitMedians(mOKcBok['raW'], mOKcBok['GrResid'], -51.5, 60.0, 112, 1)


median: -0.0005412109918604639 std.dev: 0.003418639846584066

Recalibrate Dec residuals


In [216]:
def correctZeropoint(X, Z, thetaAll):
    xLeft = -1.3
    for i in range(0,12):
        xRight = -1.2655 + (i+1)*0.2109
        theta = thetaAll[i]
        xB = X[(X>xLeft)&(X<xRight)]
        zB = Z[(X>xLeft)&(X<xRight)]
        Zfit = polynomial_fit(theta, xB)
        xLeft = xRight
        if (i==0):
            Xcorr = xB
            Zcorr = zB - Zfit
        else: 
            Xcorr = np.concatenate((Xcorr, xB), axis=0)
            Zcorr = np.concatenate((Zcorr, zB-Zfit), axis=0)

    return Xcorr, Zcorr

In [180]:
aprint(np.size(mOKcBok), np.size(mOKcBok2))


521594 505506

In [256]:
decOK = mOKcBok2['dec_new']
GrResid =  mOKcBok2['GrResid']
fig,ax = plt.subplots(1,1,figsize=(8,6))
ax.scatter(decOK, GrResid, s=0.01, c='blue')
ax.set_xlim(-1.3,1.3)
ax.set_ylim(-0.06,0.06)
ax.set_xlabel('Declination (deg)')
ax.set_ylabel('Gaia G - SDSS G')
xBin, nPts, medianBin, sigGbin = fitMedians(decOK, GrResid, -1.266, 1.264, 252, 0)
ax.scatter(xBin, medianBin, s=30.0, c='black', alpha=0.9)
ax.scatter(xBin, medianBin, s=15.0, c='yellow', alpha=0.5)
TwoSigP = medianBin + 2*sigGbin
TwoSigM = medianBin - 2*sigGbin 
ax.plot(xBin, TwoSigP, c='yellow')
ax.plot(xBin, TwoSigM, c='yellow')
xL = np.linspace(-100,100)
ax.plot(xL, 0*xL+0.00, c='yellow')
ax.plot(xL, 0*xL+0.01, c='red')
ax.plot(xL, 0*xL-0.01, c='red')
dCleft = -1.3
ax.plot(0*xL+dCleft, xL, c='red')
alltheta = []
for i in range(0,12):
    decCol = -1.2655 + (i+1)*0.2109
    ax.plot(0*xL+decCol, xL, c='red')
    xR = xBin[(xBin>dCleft)&(xBin<decCol)]
    yR = medianBin[(xBin>dCleft)&(xBin<decCol)]
    dyR = sigGbin[(xBin>dCleft)&(xBin<decCol)]
    data = np.array([xR, yR, dyR])
    theta2 = best_theta(data,5)
    alltheta.append(theta2)
    yfit = polynomial_fit(theta2, xR)
    ax.plot(xR, yfit, c='cyan', lw=2)
    dCleft = decCol
    rrr = yR - yfit
    # print(i, np.median(rrr), np.std(rrr))  # 2 milimag scatter 
    # print(i, theta2)



In [217]:
# derive corrections as a function of Dec 
decArr = np.array(mOKcBok2['dec_new'])
GrResidArr =  np.array(mOKcBok2['GrResid'])
decCorr, GrResidCorr = correctZeropoint(decArr, GrResidArr, alltheta)


Out[217]:
array([-0.00513725,  0.00434107, -0.01067382, ..., -0.08163946,
        0.02329312,  0.00564562])

In [218]:
decOK = decCorr
GrResid = GrResidCorr
fig,ax = plt.subplots(1,1,figsize=(8,6))
ax.scatter(decOK, GrResid, s=0.01, c='blue')
ax.set_xlim(-1.3,1.3)
ax.set_ylim(-0.06,0.06)
ax.set_xlabel('Declination (deg)')
ax.set_ylabel('Gaia G - SDSS G')
xBin, nPts, medianBin, sigGbin = fitMedians(decOK, GrResid, -1.266, 1.264, 252, 0)
ax.scatter(xBin, medianBin, s=30.0, c='black', alpha=0.9)
ax.scatter(xBin, medianBin, s=15.0, c='yellow', alpha=0.5)
TwoSigP = medianBin + 2*sigGbin
TwoSigM = medianBin - 2*sigGbin 
ax.plot(xBin, TwoSigP, c='yellow')
ax.plot(xBin, TwoSigM, c='yellow')
xL = np.linspace(-100,100)
ax.plot(xL, 0*xL+0.00, c='yellow')
ax.plot(xL, 0*xL+0.01, c='red')
ax.plot(xL, 0*xL-0.01, c='red')
dCleft = -1.3
ax.plot(0*xL+dCleft, xL, c='red')
for i in range(0,12):
    decCol = -1.2655 + (i+1)*0.2109
    ax.plot(0*xL+decCol, xL, c='red')
    xR = xBin[(xBin>dCleft)&(xBin<decCol)]
    yR = medianBin[(xBin>dCleft)&(xBin<decCol)]
    dyR = sigGbin[(xBin>dCleft)&(xBin<decCol)]
    data = np.array([xR, yR, dyR])
    theta2 = best_theta(data,3)
    yfit = polynomial_fit(theta2, xR)
    ax.plot(xR, yfit, c='cyan', lw=2)
    dCleft = decCol
    rrr = yR - yfit
    print(i, np.median(rrr), np.std(rrr))  # 2 milimag scatter 
    # print(i, theta2)


0 -6.328293341071765e-05 0.0008897334930961753
1 3.868832921149882e-05 0.0007555771868818445
2 -2.3085644507424372e-05 0.001326947358511655
3 0.00011695199956686733 0.0006424957583309187
4 -0.00010382885183050049 0.000989211765621928
5 -7.529043206195634e-05 0.00043862247407898397
6 0.00013627982081246303 0.0010755662228955484
7 -2.0934097579194955e-05 0.00037227220392653624
8 0.00010158709132892033 0.0006089783558551822
9 2.0550015710149583e-05 0.00045569048111683883
10 5.1899683392996054e-05 0.0007653208176340183
11 -0.0001359306520898329 0.0005977744384741594

In [250]:
# let's now correct all mags with this correction
# derive corrections as a function of Dec 
thetaRecalib = alltheta
mOKtest = mOK[mOK['r_Nobs_new']>3]
decArray = np.array(mOKtest['dec_new'])
mOKtest['GrResid'] = mOKtest['dgGr'] - np.median(mOKtest['dgGr'])
correctionArr =  np.array(mOKtest['GrResid'])
decCorr, corrections = correctZeropoint(decArray, correctionArr, thetaRecalib)
printStats(corrections)
# matches['gGr_new'] = matches['Gmag'] - matches['r_mMed_new']


            -6.449147447645764 13.46681232616975 -0.0002827079476654219 10021.08671872744 892196

In [272]:
decLeft = -1.3
for i in range(0,12):
    decRight = -1.2655 + (i+1)*0.2109
    decArr = np.linspace(decLeft, decRight, 100)
    thetaBin = thetaRecalib[i] 
    ZPfit = polynomial_fit(thetaBin, decArr)
    if (i==0):
        decCorrGrid = decArr
        ZPcorr = ZPfit
    else: 
        decCorrGrid = np.concatenate([decCorrGrid, decArr]) 
        ZPcorr = np.concatenate([ZPcorr, ZPfit])
    decLeft = decRight

In [313]:
mOKtest = mOK[mOK['r_Nobs_new']>3]
decGrid2correct = mOKtest['dec']
ZPcorrectionsDec = np.interp(decGrid2correct, decCorrGrid, ZPcorr)
raWGrid2correct = mOKtest['raW'] 
ZPcorrectionsDec = np.interp(raWGrid2correct, RAbin, RAmedianBin)

In [314]:
fig,ax = plt.subplots(1,1,figsize=(8,6))
ax.scatter(decGrid2correct, ZPcorrectionsDec, s=0.01, c='blue')
ax.plot(decCorrGrid, ZPcorr, c='red')
ax.set_xlim(-1.3,1.3)
ax.set_ylim(-0.02,0.02)
ax.set_xlabel('Declination (deg)')
ax.set_ylabel('Correction')


Out[314]:
Text(0, 0.5, 'Correction')

In [304]:
mOKtest = mOK[mOK['r_Nobs_new']>3]

In [305]:
thetaCtest, DcBokTest_new = doGaiaAll(mOKtest, 'gGr_new')


medianAll: -2.4249865089840533e-05 std.dev.All: 0.019276353963555082
N= 521594 min= -0.09528747200368677 max= 0.10470693719856021
median: 0.0003502347122703881 std.dev: 0.00632767845062075
medianAll: -2.4249865089840533e-05 std.dev.All: 0.019276353963555082
N= 521594 min= -0.09528747200368677 max= 0.10470693719856021
median: -0.0005412109918604639 std.dev: 0.003418639846584066

In [306]:
mOKtest['ZPcorrections'] = ZPcorrections
mOKtest['r_mMed_new'] = mOKtest['r_mMed_new'] + mOKtest['ZPcorrections']
mOKtest['gGr_new'] = mOKtest['Gmag'] - mOKtest['r_mMed_new']
mOKtest['gGrFit'] = polynomial_fit(thetaC, mOKtest['gi'])
mOKtest['dgGr'] = mOKtest['gGr_new'] - mOKtest['gGrFit']

In [307]:
thetaCtest, DcBokTest_new = doGaiaAll(mOKtest, 'gGr_new')


medianAll: -1.740482951925401e-05 std.dev.All: 0.01726075645822507
N= 521724 min= -0.09563441220592167 max= 0.10436005687996558
median: -4.2614853139463435e-05 std.dev: 0.0007056170351647916
medianAll: -1.740482951925401e-05 std.dev.All: 0.01726075645822507
N= 521724 min= -0.09563441220592167 max= 0.10436005687996558
median: -0.0002705818761659425 std.dev: 0.0034037088581967954

In [308]:
thetaCtest, DcBokTest_new = doGaiaAll(mOK, 'gGr_new')


medianAll: -2.4249865089840533e-05 std.dev.All: 0.019276353963555082
N= 521594 min= -0.09528747200368677 max= 0.10470693719856021
median: 0.0003502347122703881 std.dev: 0.00632767845062075
medianAll: -2.4249865089840533e-05 std.dev.All: 0.019276353963555082
N= 521594 min= -0.09528747200368677 max= 0.10470693719856021
median: -0.0005412109918604639 std.dev: 0.003418639846584066

In [ ]: