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
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
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)
Out[6]:
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)
Out[7]:
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))
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)
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))
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))
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)
In [88]:
qpBM(mOK, 'decG', -1.3, 1.3, 'ddecGnew', -0.6, 0.6, 120)
In [89]:
qpBM(mOK, 'raW', -51.5, 60, 'ddecGold', -0.6, 0.6, 120)
In [94]:
qpBM(mOK, 'raW', -51.5, 60, 'draGnew', -0.6, 0.6, 120)
In [95]:
qpBM(mOK, 'raW', -51.5, 60, 'draGold', -0.6, 0.6, 120)
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')
In [299]:
thetaC, DcBok_new = doGaiaAll(mOK, 'gGr_new')
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]:
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]:
In [138]:
xBin3, nPts3, medianBin3, sigGbin3 = fitMedians(mOKc['Gmag'], mOKc['dgGr'], 14.0, 21.0, 70, 1)
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)
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))
In [150]:
sigG(mOKcBok['GrResid'])
Out[150]:
In [156]:
qpBM(mOKcBok, 'dec', -1.3, 1.3, 'GrResid', -0.03, 0.03, 260)
In [174]:
qpBM(mOKcBok, 'raW', -51.5, 60, 'GrResid', -0.03, 0.03, 112)
In [312]:
RAbin, RAnPts, RAmedianBin, RAsigGbin = fitMedians(mOKcBok['raW'], mOKcBok['GrResid'], -51.5, 60.0, 112, 1)
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))
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]:
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)
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']
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]:
In [304]:
mOKtest = mOK[mOK['r_Nobs_new']>3]
In [305]:
thetaCtest, DcBokTest_new = doGaiaAll(mOKtest, 'gGr_new')
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')
In [308]:
thetaCtest, DcBokTest_new = doGaiaAll(mOK, 'gGr_new')
In [ ]: