This notebook examines the photometric reliability of the Bok data using repeat measurements of SDSS standard stars. For each reference star the mean and rms magnitudes are calculated from all individual measurements. The reference stars are then grouped into bins of magnitude. The distribution of the rms values (the photometric scatter) within the magnitude bins gives an indication of the photometric reliability.
First, some preliminaries:
In [2]:
%pylab inline
from matplotlib import ticker
import os,sys
sys.path.append('../bok')
from collections import OrderedDict
from astropy.table import Table,join
from astropy.stats import sigma_clip
import bokrmphot
_ = os.environ.setdefault('BOKRMDIR','../bok')
Load a photometric statistics table. This table contains the statistics of the photometry within magnitude bins. The "season" keyword restricts the data to a single observing season, so that the statistics can be compared between seasons. The resulting table is saved to a FITS file so that it can be reloaded later.
In [3]:
def load_pstable(psTabFn,photRun=None,season='2014',**kwargs):
if not os.path.exists(psTabFn):
psTab = bokrmphot.binned_phot_stats(photRun,season=season,**kwargs)
psTab.write(psTabFn)
psTab = Table.read(psTabFn)
for c in psTab.colnames:
if c.startswith('sig'):
psTab[c].format = "{:.4f}"
elif c=='outlierFrac':
psTab[c].format = "{:.5f}"
return psTab
Plot the results from the photometric statistics table. A central line shows the median rms within magnitude bins, and a shaded region shows the inter-quartile range of the rms. I.e., half of all reference stars at a given magnitude had an rms value less than the median line.
In [40]:
def plot_compare_scatters(psTabs,labels,colors='rgbcymk',filt='g',ax=None,units='mag',sfx=None):
if sfx is None:
sfx = ['','']
def _append_arr(arr):
return arr
# used this for drawstyle=steps-post, but no equiv. for fill_between
#return np.concatenate([arr,[arr[-1]]])
if units=='mag':
mscl = 1.0
elif units=='mmag':
mscl = 1e3
if ax is None:
figure()
ax = subplot(111)
for tab,l,c,_sfx in zip(psTabs,labels,colors,sfx):
ax.fill_between(tab['mbins'],_append_arr(mscl*tab['sig%s25'%_sfx]),
_append_arr(mscl*tab['sig%s75'%_sfx]),
edgecolor='none',color=c,alpha=0.5)
ax.plot(tab['mbins'],_append_arr(mscl*tab['sig%s50'%_sfx]),color=c,lw=1.5,label=l)
legend(loc='upper left')
xlabel('%s magnitude' % filt)
ylabel('per-object std [%s]'%units)
xlim(17,19)
The next set of photometric statistics come from some early reductions of the first-year Bok data. Nov2015g was an early run of bokpipe. Note that the first round of difference imaging was based on an IDL pipeline; the associated catalogs are versioned 20140905 and have an asymptotic scattter of ~2.5% in g band, much worse than the bokpipe runs. Jan2017 and Feb2017 were subsequent runs with incremental improvements; the changes can be recovered from the git repository.
In [3]:
psNov2015g = load_pstable('../phot_stats_Nov2015g.fits')
print psNov2015g['mbins','sig50','sig50Clip','outlierFrac']
In [4]:
psJan2017g = load_pstable('../phot_stats_Jan2017g.fits')
print psJan2017g['mbins','sig50','sig50Clip','outlierFrac']
In [5]:
psFeb2017g = load_pstable('../phot_stats_Feb2017g.fits')
print psFeb2017g['mbins','sig50','sig50Clip','outlierFrac']
In [6]:
psJan2017i = load_pstable('../phot_stats_Jan2017i.fits')
print psJan2017i['mbins','sig50','sig50Clip','outlierFrac']
In [7]:
psFeb2017i = load_pstable('../phot_stats_Feb2017i.fits')
print psFeb2017i['mbins','sig50','sig50Clip','outlierFrac']
The improvement from the Nov2015g to Jan2017g versions is significant.
In [9]:
plot_compare_scatters([psNov2015g,psJan2017g],['Nov15g','Jan17g'],'gb')
Not much improvement in the Feb2017g version.
In [10]:
plot_compare_scatters([psJan2017g,psFeb2017g],['Jan17g','Feb17g'],'gb')
The scatter increases significantly if fainter stars are included in the zeropoint. Going from mag$<19.5$ to mag$<20.5$ more than doubles the scatter. Do not understand why this is...
In [11]:
plot_compare_scatters([psJan2017i,psFeb2017i],['Jan17i','Feb17i'],'gb')
In [14]:
for psTab,l in zip([psNov2015g,psJan2017g,psFeb2017g],['Nov15g','Jan17g','Feb17g']):
ii = where(psTab['outlierFrac']>0)[0]
plot(psTab['mbins'][ii],log10(psTab['outlierFrac'][ii]),label=l)
legend()
Out[14]:
In [15]:
for psTab,l in zip([psJan2017i,psFeb2017i],['Jan17i','Feb17i']):
ii = where(psTab['outlierFrac']>0)[0]
plot(psTab['mbins'][ii],log10(psTab['outlierFrac'][ii]),label=l)
legend()
Out[15]:
In [5]:
procJul2017 = OrderedDict()
for season in ['2014','2015','2016']:
procJul2017[season] = load_pstable('../archive/summer_2017/phot_stats_%sg_v2.fits'%season,
'cleanstars',season,catdir='../archive/summer_2017/')
print season
print procJul2017[season]['mbins','sig50','outlierFrac']
Comparing the results for the three observing seasons. The 2014 data has less scatter, possibly because there are more data and the sky flats are more accurate.
In [7]:
plot_compare_scatters(procJul2017.values(),procJul2017.keys(),'gbr')
In [20]:
procOct2017 = OrderedDict()
for season in ['2014','2015','2016','2017']:
procOct2017[season] = load_pstable('../archive/october_2017/phot_stats_%sg_v2.fits'%season,
'cleanstars',season,catdir='../archive/october_2017/')
print season
print procOct2017[season]['mbins','sig50','outlierFrac']
This compares the Oct2017 processing to the Jul2017 processing. The main improvement was in cleaning up some gain correction errors, and fixing a long-standing bug where the zeropoints were based on an aperture (15 pixel) that subsequently had an aperture correction applied. A change was made to apply the correction to the zeropoint. Now the 2014 data reach sub-1% repeatability at $g<17.7$!
In [21]:
figure(figsize=(14,4))
subplots_adjust(0.05,0.1,0.99,0.99,0.22)
for pnum,season in enumerate(['2014','2015','2016'],start=1):
ax = subplot(1,3,pnum)
plot_compare_scatters([procJul2017[season],procOct2017[season]],
[season+'-Jul2017',season+'-Oct2017'],'gb',ax=ax)
ax.set_ylim(0,0.04)
ax.grid()
Again comparing the four seasons. The 2017 data mainly have 300s exposure times, whereas previous years mainly had 150s exposures.
In [15]:
plot_compare_scatters(procOct2017.values(),procOct2017.keys(),'gbrm')
And now the i-band:
In [22]:
procOct2017i = OrderedDict()
for season in ['2014','2015','2016','2017']:
procOct2017i[season] = load_pstable('../archive/october_2017/phot_stats_%si_v2.fits'%season,
'cleanstars',season,band='i',
catdir='../archive/october_2017/')
print season
print procOct2017i[season]['mbins','sig50','outlierFrac']
In [23]:
plot_compare_scatters(procOct2017i.values(),procOct2017.keys(),'gbrm')
So the 2016 i-band data has some issues...
As a shortcut to solving the full ubercalibration problem (even with a large number of repeats the Bok data have been found to be somewhat noisy for recovering ubercal parameters), a selfcal() routine was implemented that assumes the extinction parameters from SDSS and uses them to derive per-image zeropoints from the mean magnitude offsets. The above results were all based on zeropoints obtained from external calibration using SDSS; the following results come from the selfcal procedure which provides zeropoints on an internal Bok system.
In [16]:
class NewBinnedStats(object):
def __init__(self,tabf):
self.tab = Table.read(tabf)
def get(self,season,filt):
i = where((self.tab['season']==season)&(self.tab['filter']==filt))[0][0]
rv = Table()
rv['mbins'] = array(self.tab.meta['MBINS'].split(',')).astype(float)
for k in self.tab.colnames:
if k.startswith('sig'):
rv[k] = self.tab[k][i]
return rv
In [28]:
oct17noselfcal = NewBinnedStats('../archive/october_2017/phot_stats_oct17noselfcal.fits')
oct17selfcal = NewBinnedStats('../archive/october_2017/phot_stats_oct17.fits')
In [35]:
figure(figsize=(14,7))
subplots_adjust(0.05,0.1,0.99,0.99,0.22)
pnum = 1
for b in 'gi':
for season in ['2014','2015','2016','2017']:
ax = subplot(2,4,pnum)
plot_compare_scatters([oct17noselfcal.get(season,b),oct17selfcal.get(season,b)],
[season+'-Oct2017orig',season+'-Oct2017selfcal'],
'gb',filt=b,ax=ax,units='mmag')
if b=='g':
ax.set_ylim(0,40)
else:
ax.set_ylim(0,80)
ax.grid()
pnum += 1
selfcal has a small effect on g-band photometry (maybe slight improvement in 2015). But huge effect on i-band, and in particular addresses the weird scatter in the 2016 data. Now asymptotically reaching 10 mmag scatter at g,i=17.
Check the improvement by restricting to only photometric frames:
In [36]:
oct17selfcalphoto = NewBinnedStats('../archive/october_2017/phot_stats_oct17_photo.fits')
In [38]:
figure(figsize=(14,7))
subplots_adjust(0.05,0.1,0.99,0.99,0.22)
pnum = 1
for b in 'gi':
for season in ['2014','2015','2016','2017']:
ax = subplot(2,4,pnum)
plot_compare_scatters([oct17selfcal.get(season,b),oct17selfcalphoto.get(season,b)],
[season+'-Oct2017selfcal',season+'-Oct2017selfcalphoto'],
'br',filt=b,ax=ax,units='mmag')
if b=='g':
ax.set_ylim(0,40)
else:
ax.set_ylim(0,80)
ax.grid()
pnum += 1
Finally, compare the internal calib with the external rms (Bok-SDSS):
In [41]:
figure(figsize=(14,7))
subplots_adjust(0.05,0.1,0.99,0.99,0.22)
pnum = 1
for b in 'gi':
for season in ['2014','2015','2016','2017']:
ax = subplot(2,4,pnum)
plot_compare_scatters([oct17selfcal.get(season,b),oct17selfcal.get(season,b)],
[season+'-Oct2017selfcal-int',season+'-Oct2017selfcal-ext'],
['C0','C1'],filt=b,ax=ax,units='mmag',sfx=['','Ext'])
if b=='g':
ax.set_ylim(0,40)
else:
ax.set_ylim(0,80)
ax.grid()
pnum += 1
In [ ]: