<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Parts adapted from DECALS low-SB_completeness AnaK overlap.ipynb
Get AnaK Info and download
Step1:  basic photometric addtions to the catalogs
Step2:  Basic residual comparisons
Step3:  Inspect objects with failed residuals
Step4:  ???? Why are so many of the residuals NaN/infs ???
Step5:  Inspect objects with high residuals vs. flux
Note that this includes only those with r<22 to ensure there's not a flux effect
Step6:  X-matching to SDSS
Step8:  now inspect those that are in the upper-left of that plot
<ASSISTANT_TASK:>
Python Code:
hsts = hosts.get_saga_hosts_from_google(clientsecretjsonorfn='client_secrets.json', useobservingsummary=False)
anak  = [h for h in hsts if h.name=='AnaK']
assert len(anak)==1
anak = anak[0]
bricknames = []
with open('decals_dr3/anakbricks') as f:
    for l in f:
        l = l.strip()
        if l != '':
            bricknames.append(l)
print(bricknames)
base_url = 'http://portal.nersc.gov/project/cosmo/data/legacysurvey/dr3/tractor/{first3}/tractor-{brickname}.fits'
for brickname in ProgressBar(bricknames, ipython_widget=True):
    url = base_url.format(brickname=brickname, first3=brickname[:3])
    target = os.path.join('decals_dr3/catalogs/', url.split('/')[-1])
    if not os.path.isfile(target):
        !wget $url -O $target
    else:
        print(target, 'already exists, not downloading')
bricks = Table.read('decals_dr3/survey-bricks.fits.gz')
bricksdr3 = Table.read('decals_dr3/survey-bricks-dr3.fits.gz')
catalog_fns = ['decals_dr3/catalogs/tractor-{}.fits'.format(bnm) for bnm in bricknames]
decals_catalogs = [Table.read(fn) for fn in catalog_fns]
dcatall = table.vstack(decals_catalogs, metadata_conflicts='silent')
sdss_catalog = Table.read('catalogs/base_sql_nsa{}.fits.gz'.format(anak.nsaid))
Explanation: Parts adapted from DECALS low-SB_completeness AnaK overlap.ipynb
Get AnaK Info and download
End of explanation
for dcat in [dcatall]:
    for magnm, idx in zip('grz', [1, 2, 4]):
        mag, mag_err = fluxivar_to_mag_magerr(dcat['decam_flux'][:, idx], dcat['decam_flux_ivar'][:, idx])
        dcat[magnm] = mag
        dcat[magnm + '_err'] = mag_err
    
    dcat['sb_r_0.5'] = compute_sb(0.5*u.arcsec, dcat['decam_apflux'][:, 2, :])
    dcat['sb_r_0.75'] = compute_sb(0.75*u.arcsec, dcat['decam_apflux'][:, 2, :])
    dcat['sb_r_1'] = compute_sb(1.0*u.arcsec, dcat['decam_apflux'][:, 2, :])
    dcat['sb_r_2'] = compute_sb(2.0*u.arcsec, dcat['decam_apflux'][:, 2, :])
Explanation: basic photometric addtions to the catalogs
End of explanation
DECALS_AP_SIZES
apmag, apmagerr = fluxivar_to_mag_magerr(dcatall['decam_apflux'], dcatall['decam_apflux_ivar'])
apmagres, _ = fluxivar_to_mag_magerr(dcatall['decam_apflux_resid'], dcatall['decam_apflux_ivar'])
apdiff = subselect_aperture(apmagres - apmag, 'r')
apcolor = subselect_aperture(apmag, 'g') - subselect(apmag, 'r')
apmagx = subselect_aperture(apmag, 'r')
good = ~np.isnan(apdiff)&~np.isnan(apcolor)&~np.isnan(apmagx)
plt.scatter(apcolor[good], apmagx[good], c=apdiff[good], cmap='viridis', alpha=.1, lw=0, s=1)
plt.colorbar()
plt.xlim(-2,5)
plt.ylim(26,15)
Explanation: Basic residual comparisons
End of explanation
for band in 'ugrizy':
    reses = subselect_aperture(apmagres, band, None)
    print('Band', band)
    for ap, res in zip(DECALS_AP_SIZES, reses.T):
        print('Aperture', ap,'has', 100*np.sum(np.isfinite(res))/len(res),'% good')
Explanation: Inspect objects with failed residuals
End of explanation
rs = subselect_aperture(apmagres, 'r')
catnotfin = dcatall[~np.isfinite(rs)]
catnotfin['apmagres_allaps'] = subselect_aperture(apmagres, 'r', None)[~np.isfinite(rs)]
make_cutout_comparison_table(catnotfin[np.random.permutation(len(catnotfin))[:10]], 
                             inclres=True, inclmod=True, inclsdss=False, doprint=False,
                             add_annotation=['apmagres_allaps'])
Explanation: ???? Why are so many of the residuals NaN/infs ???
End of explanation
dmag_of_ap_distr = {}
for ap in DECALS_AP_SIZES:
    rs = subselect_aperture(apmag, 'r', ap)
    rres = subselect_aperture(apmagres, 'r', ap)
    dmag_of_ap_distr[ap] = dmag = rs - rres
    plt.hist(dmag[np.isfinite(dmag)&(rs<22*u.mag)], bins=100, histtype='step', label=str(ap), normed=True)
    
plt.legend(loc=0)
plt.xlabel('r_flux - r_resid')
ap = 1.0*u.arcsec
rs = subselect_aperture(apmag, 'r', ap)
rres = subselect_aperture(apmagres, 'r', ap)
dmag = rs - rres
perc = 95
p = np.percentile(dmag[np.isfinite(dmag)&(rs<22*u.mag)], perc)
print('nobjs in', perc,'percentile:', np.sum(dmag.value>p), 'cutoff is', p)
msk = np.isfinite(dmag)&(dmag.value>p)&(rs<22*u.mag)
dcatbadres = dcatall[msk]
dcatbadres['dmag'] = dmag[msk]
dcatbadres['r'] = rs[msk]
make_cutout_comparison_table(dcatbadres[:10], 
                             inclres=True, inclmod=True, inclsdss=False, doprint=False,
                             add_annotation=['dmag', 'r'])
Explanation: Inspect objects with high residuals vs. flux
Note that this includes only those with r<22 to ensure there's not a flux effect
End of explanation
#cut out the non-overlap region
dsc = SkyCoord(dcatall['ra'], dcatall['dec'], unit=u.deg)
dcutall = dcatall[dsc.separation(anak.coords) < 1*u.deg]
dsc = SkyCoord(dcutall['ra'], dcutall['dec'], unit=u.deg)
ssc = SkyCoord(sdss_catalog['ra'], sdss_catalog['dec'], unit=u.deg)
threshold = 1*u.arcsec
idx, d2d, _ = ssc.match_to_catalog_sky(dsc)
plt.hist(d2d.arcsec, bins=100, range=(0, 3),histtype='step', log=True)
plt.axvline(threshold.to(u.arcsec).value, c='k')
None
dmatchmsk = idx[d2d<threshold]
dmatch = dcutall[dmatchmsk]
smatch = sdss_catalog[d2d<threshold]
idx, d2d, _ = dsc.match_to_catalog_sky(ssc)
dnomatchmsk = d2d>threshold
dnomatch = dcutall[dnomatchmsk]
plt.figure(figsize=(12, 10))
xnm = 'r'
ynm = 'sb_r_0.5'
ap = 1*u.arcsec
apmag, apmagerr = fluxivar_to_mag_magerr(dnomatch['decam_apflux'], dnomatch['decam_apflux_ivar'])
apmagres, _ = fluxivar_to_mag_magerr(dnomatch['decam_apflux_resid'], dnomatch['decam_apflux_ivar'])
rs = subselect_aperture(apmag, xnm, ap)
rres = subselect_aperture(apmagres, xnm, ap)
dmag = rs - rres
dnstar = dnomatch['type']=='PSF '
dnoext = -2.5*np.log10(dnomatch['decam_mw_transmission'][:, 2])
r0 = (dnomatch[xnm] - dnoext)
sb = dnomatch[ynm] - dnoext
plt.scatter(r0[~dnstar], sb[~dnstar], 
            c=dmag[~dnstar], lw=0, alpha=1, s=3, label='Glx in DECALS, not in SDSS', vmax=0, vmin=-5,
            cmap='viridis_r')
plt.colorbar().set_label('r_ap - r_res [{}]'.format(ap))
plt.axvline(20.75, color='k', ls=':')
plt.xlim(17, 23)
plt.ylim(18, 28)
plt.xlabel(r'$r_{0, {\rm DECaLS}}$', fontsize=28)
plt.ylabel(r'$SB_{0.5^{\prime \prime}, {\rm DECaLS}}$', fontsize=28)
plt.xticks(fontsize=24)
plt.yticks(fontsize=24)
plt.legend(loc='lower right', fontsize=20)
Explanation: X-matching to SDSS
End of explanation
msk = (r0[~dnstar]<20.75)&(sb[~dnstar]>24)
cat = dnomatch[~dnstar][msk]
cat['dmag'] = dmag[~dnstar][msk]
p = np.percentile(cat['dmag'][np.isfinite(cat['dmag'])], 10)
catlower = cat[cat['dmag']<p]
print(len(catlower))
make_cutout_comparison_table(catlower[np.random.permutation(len(catupper))[:10]], 
                             inclres=True, inclmod=True, inclsdss=False, doprint=False,
                             add_annotation=['dmag', 'r'])
p = np.percentile(cat['dmag'][np.isfinite(cat['dmag'])], 90)
catupper = cat[cat['dmag']>p]
print(len(catupper))
make_cutout_comparison_table(catupper[np.random.permutation(len(catupper))[:10]], 
                             inclres=True, inclmod=True, inclsdss=False, doprint=False,
                             add_annotation=['dmag', 'r'])
# things we want to identify automatically
disky_things_from_marla = 
mag ra dec unk
19.4785    354.157853072    0.102747198705    false
19.3311    354.255069696    0.619242808225    false
19.1127    354.284237813    0.160859730123    false
18.6914    354.069400831    -0.115904590235    false
19.3392    354.415038652    -0.0645766794736    false
19.4624    354.114525096    0.00532483801292    false
19.1087    354.534841322    0.436958919955    false
19.0242    354.447705125    0.266811924681    false
19.0354    354.136706143    0.691858529943    false
19.0534    354.53888635    0.236989192453    false
19.2452    354.568916976    0.837946572935    false
19.3481    353.473924073    -0.0912764847886    false
19.268    354.011615043    -0.445983276629    false
19.3429    354.408722681    0.904017267882    true
19.1914    354.851596507    0.401264434888    true
19.2529    353.412613626    0.640755638979    false
19.3848    354.878706758    0.364835196656    false
disky_things_from_marla = Table.read([disky_things_from_marla], format='ascii')
sc_marla = SkyCoord(disky_things_from_marla['ra'], disky_things_from_marla['dec'], unit=u.deg)
idx, d2d, _ = sc_marla.match_to_catalog_sky(SkyCoord(dcatall['ra'], dcatall['dec'], unit=u.deg))
np.sum(d2d < 1*u.arcsec)/len(d2d)
matchcat = dcatall[idx]
apmag, apmagerr = fluxivar_to_mag_magerr(matchcat['decam_apflux'], matchcat['decam_apflux_ivar'])
apmagres, _ = fluxivar_to_mag_magerr(matchcat['decam_apflux_resid'], matchcat['decam_apflux_ivar'])
rs = subselect_aperture(apmag, xnm, None)
rres = subselect_aperture(apmagres, xnm, None)
matchcat['dmag'] = rs - rres
dmagsigs = []
for dmag, dmagdistr in zip(matchcat['dmag'].T, dmag_of_ap_distr.values()):
    msk = np.isfinite(dmagdistr)
    dmagsigs.append(np.mean(dmagdistr[msk].value)-dmag/np.std(dmagdistr[msk].value))
    print(np.mean(dmagdistr[msk].value).shape, dmag.shape, np.std(dmagdistr[msk].value).shape, dmagsigs[-1].shape)
matchcat['dmag_sig'] = np.array(dmagsigs).T
make_cutout_comparison_table(matchcat[np.random.permutation(len(matchcat))[:10]], 
                             inclres=True, inclmod=True, inclsdss=False, doprint=False,
                             add_annotation=['dmag', 'r', 'dmag_sig'])
Explanation: now inspect those that are in the upper-left of that plot
End of explanation
<END_TASK>