In [1]:
# This changes the current directory to the base saga directory - make sure to run this first!
# This is necessary to be able to import the py files and use the right directories,
# while keeping all the notebooks in their own directory.
import os
import sys
import numpy as np
if 'saga_base_dir' not in locals():
saga_base_dir = os.path.abspath('..')
if saga_base_dir not in sys.path:
os.chdir(saga_base_dir)
In [2]:
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table
from astropy.io import fits
In [3]:
%matplotlib inline
from matplotlib import rcParams, style
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
rcParams['image.cmap'] = 'viridis'
rcParams['image.origin'] = 'lower'
rcParams['axes.prop_cycle'] = style.library['seaborn-deep']['axes.prop_cycle']
rcParams['figure.figsize'] = (14, 8)
rcParams['axes.titlesize'] = rcParams['axes.labelsize'] = 16
rcParams['xtick.labelsize'] = rcParams['ytick.labelsize'] = 14
# this final line just shows the first color that plt.plot uses.
# plt.bar and plt.scatter don't use this so we'll use it explicitly below.
next(iter(rcParams['axes.prop_cycle']))['color']
Out[3]:
In [4]:
dirtytab = Table.read('SAGADropbox/data/saga_spectra_dirty.fits.gz')
In [5]:
print len(dirtytab)
dirtytab[:5].show_in_notebook()
Out[5]:
Lets check what ZQUALITY values are present
In [6]:
np.unique(dirtytab['ZQUALITY'])
Out[6]:
In [7]:
np.sum(dirtytab['ZQUALITY'].mask)
Out[7]:
~1000 of the of the ZQUALITY entries are masked (that's what the "--" means in the np.unique call). Lets check their value...
In [8]:
zqmsked = dirtytab['ZQUALITY'][dirtytab['ZQUALITY'].mask]
zqmsked.mask = False
np.unique(zqmsked)
Out[8]:
Hmm, all of the ZQUALITY values that are 1 are masked. That's probably not intended, so we'll un-mask them
In [9]:
dirtytab['ZQUALITY'].mask = False
np.unique(dirtytab['ZQUALITY'])
Out[9]:
In [10]:
for zq in np.unique(dirtytab['ZQUALITY']):
print('ZQUALITY={} is {:%}'.format(zq, np.sum(zq==dirtytab['ZQUALITY'])/len(dirtytab)))
Lets just have a quick look at the CMD
In [11]:
plt.scatter(dirtytab['g'] - dirtytab['r'], dirtytab['r'], c='k', s=1, lw=0, alpha=.25)
plt.xlim(-1,4)
plt.ylim(24,15)
Out[11]:
In [12]:
okspec_msk = dirtytab['ZQUALITY']>2
gmr = dirtytab['g'] - dirtytab['r']
gmr_okspec = gmr[okspec_msk]
gmr_badspec = gmr[~okspec_msk]
In [13]:
# not that these are the bin *edges*, so even though there are 11 elements, we end up with 10 bins
completeness_bins = np.linspace(0,2.5, 11)
completeness_bin_centers = (completeness_bins[1:] + completeness_bins[:-1])/2
n_okspec_inbin = np.histogram(gmr_okspec, completeness_bins)[0]
n_badspec_inbin = np.histogram(gmr_badspec, completeness_bins)[0]
n_total_inbin = n_okspec_inbin + n_badspec_inbin
In [14]:
# note that bar takes the *left* edge of a bin, which is why we use completeness_bins[:-1]
plt.bar(completeness_bins[:-1], n_okspec_inbin/n_total_inbin, width=np.diff(completeness_bins), color='#4C72B0')
plt.xlabel('$g-r$')
plt.ylabel('completeness fraction')
plt.tight_layout()
In [15]:
plt.scatter(completeness_bin_centers, n_okspec_inbin/n_total_inbin, color='#4C72B0')
plt.ylim(0,1)
plt.xlabel('$g-r$')
plt.ylabel('completeness fraction')
plt.tight_layout()