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

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]:
import hosts
import targeting
import mmthecto

from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table

In [77]:
%matplotlib inline
from matplotlib import pyplot as plt
from matplotlib import rcParams

rcParams['image.interpolation'] = 'none'
rcParams['figure.figsize'] = (16, 10)

This is for AnaK using the Stripe 82 catalog

Ana K


In [17]:
hostlst = hosts.get_saga_hosts_from_google('etollerud@gmail.com', passwd)
anak = [h for h in hostlst if h.name=='AnaK'][0]

In [9]:
# NOTE: this was originally done with a different S82 catalog w/ various issues

In [18]:
s82catfull = fits.getdata('nsa_61945_stripe82coadd.fits')
s82catfull.dtype


Out[18]:
dtype([('ra', '>f8'), ('dec', '>f8'), ('col', '>i2'), ('fieldNum', '>i2'), ('htm15', '>i8'), ('fieldId', '>i4', (5,)), ('E(B-V)', '>f4'), ('aperFlux', '>f4', (5,)), ('aperFluxErr', '>f4', (5,)), ('aperMag', '>f4', (5,)), ('aperMagErr', '>f4', (5,)), ('autoMag', '>f4', (5,)), ('autoMagErr', '>f4', (5,)), ('class_star', '>f4', (5,)), ('flags', '>i2', (5,))])

The original catalog has tons of repeats. Clean it out based on exact matches in ra. New version does'nt have them, though


In [19]:
uras = np.unique(s82catfull['ra'])
msk = np.zeros(len(s82catfull), dtype=bool)
for ura in uras:
    firstura = np.where(ura==s82catfull['ra'])[0][0]
    msk[firstura] = True
s82cat = s82catfull[msk]

In [20]:
len(s82catfull), len(s82cat)


Out[20]:
(9298, 9298)

In [21]:
um,gm,rm,im,zm = s82cat['autoMag'].T
ume,ge,re,ie,ze = s82cat['autoMagErr'].T
s82coo = SkyCoord(s82cat['ra']*u.deg, s82cat['dec']*u.deg)

In [22]:
#now check out the color cuts

colormsk82 = np.ones(len(s82cat), dtype=bool)
for color, (lower, upper) in targeting.bossanova_color_cuts.items():
    c1nm, c2nm = color.split('-')
    c1 = locals()[c1nm+'m']
    c2 = locals()[c2nm+'m']
    c = c1 - c2
    
    colormsk82[c<lower] = False
    colormsk82[c>upper] = False



colormsk82tighter = np.ones(len(s82cat), dtype=bool)
for color, (lower, upper) in targeting.tighter_color_cuts.items():
    c1nm, c2nm = color.split('-')
    c1 = locals()[c1nm+'m']
    c2 = locals()[c2nm+'m']
    c = c1 - c2
    
    colormsk82tighter[c<lower] = False
    colormsk82tighter[c>upper] = False
np.sum(colormsk82), np.sum(colormsk82tighter), len(colormsk82)


Out[22]:
(5634, 3851, 9298)

In [23]:
oldtargs = targeting.select_targets(anak, outercutrad=360*u.kpc, innercutrad=20,
                                 colorcuts=targeting.bossanova_color_cuts,
                                 galvsallcutoff=21, removegalsathighz=True, 
                                 fibermagcut=('r', 22.5))
tcoo = SkyCoord(oldtargs['ra']*u.deg, oldtargs['dec']*u.deg)


Removing 1058 objects at high z w/ good spectra, keeping 17 (total of 85368 objects)

In [24]:
len(s82coo), len(tcoo)


Out[24]:
(9298, 2303)

In [23]:
plt.scatter((s82coo[colormsk82].ra-anak.coords.ra).arcmin, (s82coo[colormsk82].dec-anak.coords.dec).arcmin, s=3,lw=0,alpha=.5, label='S82 good color')
plt.scatter((s82coo[~colormsk82].ra-anak.coords.ra).arcmin, (s82coo[~colormsk82].dec-anak.coords.dec).arcmin, s=3,lw=0,alpha=.5, c='r', label='S82 bad color')
plt.scatter((tcoo.ra-anak.coords.ra).arcmin, (tcoo.dec-anak.coords.dec).arcmin, s=3,lw=0,alpha=1,c='g', label='DR10')
plt.legend(loc=0)


Out[23]:
<matplotlib.legend.Legend at 0x12bf7d9d0>

In [25]:
idx, d2d, d3d = tcoo.match_to_catalog_sky(s82coo[colormsk82])

In [25]:
plt.subplot(1, 2, 1)
plt.hist(d2d.arcmin,bins=100,range=(0,2),histtype='step')
plt.xlabel('arcmin')
plt.subplot(1, 2, 2)
plt.hist(d2d.arcsec,bins=100,range=(0,1),histtype='step')
plt.xlabel('arcsec')
plt.tight_layout()
np.sum(d2d< 0.5*u.arcsec), np.sum(d2d> 0.5*u.arcsec), np.sum(d2d> 1*u.arcmin)


Out[25]:
(1428, 875, 640)

In [26]:
dra = (tcoo.ra - s82coo[colormsk82][idx].ra).arcsec
ddec = (tcoo.dec - s82coo[colormsk82][idx].dec).arcsec
plt.subplot(1,2,1)
plt.scatter(dra, ddec,lw=0,s=5)
plt.xlim(-100, 100)
plt.ylim(-100, 100)
plt.subplot(1,2,2)
plt.scatter(dra, ddec,lw=0,s=5)
plt.scatter([0], [0], c='r')
plt.xlim(-0.5, 0.5)
plt.ylim(-0.5, 0.5)
plt.tight_layout()



In [27]:
#compare mags for those that do match
matchmsk = d2d < 0.5*u.arcsec
dr = oldtargs[matchmsk]['r'] - s82cat['autoMag'][colormsk82][idx[matchmsk], 2]
plt.scatter(oldtargs[matchmsk]['r'], dr, lw=0, s=8, alpha=.8, c='r')
plt.axhline(0, color='k',ls=':')
plt.xlabel('$r$',fontsize=48)
plt.ylabel('$r_{\\rm DR10} - r_{\\rm S82}$',fontsize=48)
plt.ylim(-.5,0.5)
plt.xlim(16,21.1)
plt.tight_layout()



In [28]:
idx2, d2d2, d3d2 = s82coo[colormsk82].match_to_catalog_sky(tcoo)

plt.subplot(1, 2, 1)
plt.hist(d2d2.arcmin,bins=100,range=(0,2),histtype='step')
plt.xlabel('arcmin')
plt.subplot(1, 2, 2)
plt.hist(d2d2.arcsec,bins=100,range=(0,1),histtype='step')
plt.xlabel('arcsec')
plt.tight_layout()
np.sum(d2d2< 0.5*u.arcsec), np.sum(d2d2> 0.5*u.arcsec), np.sum(d2d2> 1*u.arcmin)


Out[28]:
(1572, 4062, 794)

In [29]:
idxself, d2dself, d3dself = s82coo.match_to_catalog_sky(s82coo, 2)

plt.subplot(1, 2, 1)
plt.hist(d2dself.arcmin,bins=100,range=(0,2),histtype='step')
plt.xlabel('arcmin')
plt.subplot(1, 2, 2)
plt.hist(d2dself.arcsec,bins=100,range=(0,1),histtype='step')
plt.xlabel('arcsec')
plt.tight_layout()
np.sum(d2dself< 0.5*u.arcsec), np.sum(d2dself> 0.5*u.arcsec), np.sum(d2dself> 1*u.arcmin)


Out[29]:
(1766, 7532, 2)

In [30]:
DEFAULT_URL='http://skyservice.pha.jhu.edu/DR10/ImgCutout/getjpeg.aspx?ra={ra}&dec={dec}&width={w}&height={h}&opt={flags}&scale={scale}'
def view_coo(coo, flags='', baseurl=DEFAULT_URL, size=(512, 512), scale=1):
    from IPython import display
    
    if not coo.isscalar:
        raise TypeError('Can only show scalar coordinates')
    pxscale = 0.396127*scale
    url = baseurl.format(ra=coo.ra.deg, dec=coo.dec.deg, flags=flags, w=size[0], h=size[1], scale=pxscale)
    return display.Image(url=url, format='jpeg')

There are some "quasi-duplicates" that are quite close

Check non-match properties


In [105]:
c = (gm - im)[colormsk82]
mag = rm[colormsk82]

matches = d2d< 0.5*u.arcsec

plt.scatter(c[matches], mag[matches], c='b', lw=0, s=2)
plt.scatter(c[~matches], mag[~matches], c='r', lw=0, s=4)
plt.xlim(-1,2)
plt.ylim(23,15)


Out[105]:
(23, 15)

In [107]:
view_coo(s82coo[d2d2> 0.5*u.arcsec][100], 'G',scale=.2)


Out[107]:

Just try making a catalog anyway


In [27]:
anadist = s82coo.separation(anak.coords)
distmsk = (anadist>anak.physical_to_projected(20*u.kpc))&(anadist < anak.environsarcmin*u.arcmin)
starmsk = s82cat['class_star'][:,2]<0.5
magmsk = s82cat['autoMag'][:,2]<21.5

In [28]:
targs = Table(data=[(np.arange(len(s82coo))+1), s82coo.ra.deg, s82coo.dec.deg, 
                    s82cat['autoMag'][:,1], s82cat['autoMag'][:,2], s82cat['autoMag'][:,3],
                    s82cat['autoMagErr'][:,1], s82cat['autoMagErr'][:,2], s82cat['autoMagErr'][:,3],
                    ],
              names=['objID', 'ra', 'dec', 'g', 'r', 'i', 'g_err', 'r_err', 'i_err'])

targranks = np.zeros(len(targs), dtype=int)
targranks[distmsk&starmsk&colormsk82] = 2
targranks[distmsk&starmsk&colormsk82tighter] = 4

msk = targranks > 0
targs = targs[msk]
targranks = targranks[msk]

In [29]:
#now add in all the stuff from dr10 that does *not* match s82
idx, d2d, d3d = tcoo.match_to_catalog_sky(s82coo)
nomatch = d2d > 0.5*u.arcsec

targstoadd = oldtargs[nomatch]
targrankstoadd = np.ones(len(targstoadd), dtype=int) * 3
gmr = oldtargs[nomatch]['g']-oldtargs[nomatch]['r']
rmi = oldtargs[nomatch]['r']-oldtargs[nomatch]['i']
targrankstoadd[(gmr<targeting.tighter_color_cuts['g-r'][1])&
               (rmi<targeting.tighter_color_cuts['r-i'][1])] +=2

ids = np.concatenate((targs['objID'], targstoadd['objID']))
ras = np.concatenate((targs['ra'], targstoadd['ra']))
decs = np.concatenate((targs['dec'], targstoadd['dec']))
gs = np.concatenate((targs['g'], targstoadd['g']))
rs = np.concatenate((targs['r'], targstoadd['r']))
is_ = np.concatenate((targs['i'], targstoadd['i']))
g_errs = np.concatenate((targs['g_err'], targstoadd['g_err']))
r_errs = np.concatenate((targs['r_err'], targstoadd['r_err']))
i_errs = np.concatenate((targs['i_err'], targstoadd['i_err']))
targs = Table(data=[ids, ras, decs, gs, rs, is_, g_errs, r_errs, i_errs], 
              names=['objID', 'ra', 'dec', 'g', 'r', 'i', 'g_err', 'r_err', 'i_err'])
targranks = np.concatenate((targranks, targrankstoadd))

In [112]:
host = anak
tab = mmthecto.generate_catalog(host, targs, targranks, fluxrank=1, repeatflux=10, 
                                fnout='mmthecto/{0}_aug2014.cat'.format(host.name),
                                fluxfnout='mmthecto/{0}_aug2014.fluxstars'.format(host.name))


Including 2519 targets
Found 175 Flux stars
Found 704 guide stars

In [113]:
for r in np.unique(tab['rank']):
    print r, np.sum(tab['rank']==r)


 704
1 1750
2 914
3 45
4 1486
5 74

In [117]:
#check out the calib stars (opens a browser window)
targeting.sampled_imagelist(tab[tab['rank']=='1'], None)
None

In [118]:
#check out the primary targets from s82 (opens a browser window)
targeting.sampled_imagelist(tab[tab['rank']=='2'], None)
None

In [328]:
#check out the secondary targets from s82 (opens a browser window)
targeting.sampled_imagelist(tab[tab['rank']=='4'], None)
None

In [ ]:
#check out the secondary targets from dr10 (opens a browser window)
targeting.sampled_imagelist(tab[tab['rank']=='3'], None)
None

In [ ]:
#check out the secondary targets from dr10 (opens a browser window)
targeting.sampled_imagelist(tab[tab['rank']=='5'], None)
None

In [49]:
targeting.bossanova_color_cuts, targeting.tighter_color_cuts


Out[49]:
({'g-r': (None, 1.3), 'r-i': (None, 0.7)},
 {'g-r': (None, 1.0), 'r-i': (None, 0.5)})

A catalog that removes already-observed


In [4]:
already_observed = fits.getdata('mmthecto/AnaK_obs_asofAug14.fits.gz')
successful_obs = already_observed[already_observed['ZQUALITY']>2]

In [30]:
obssc = SkyCoord(successful_obs['RA']*u.deg, successful_obs['DEC']*u.deg)
targssc = SkyCoord(targs['ra']*u.deg, targs['dec']*u.deg)

In [31]:
idx, d2d, d3d = targssc.match_to_catalog_sky(obssc)

In [32]:
notprevobsmsk = d2d > 0.5*u.arcsec
targsnotprevobs = targs[notprevobsmsk]
targsnotprevobssc = targssc[notprevobsmsk]
np.sum(notprevobsmsk), np.sum(~notprevobsmsk), len(d2d)


Out[32]:
(2447, 777, 3224)

In [80]:
plt.scatter(successful_obs['RA'], successful_obs['DEC'],lw=0,c='r',s=5)
plt.scatter(targssc.ra.deg[notprevobsmsk], targssc.dec.deg[notprevobsmsk],lw=0,c='g',s=5)
plt.scatter(targssc.ra.deg[~notprevobsmsk], targssc.dec.deg[~notprevobsmsk],lw=0,c='k',s=5)
plt.scatter(s82cat['ra'], s82cat['dec'],lw=0,c='c',s=5)
plt.tight_layout()