In [1]:
from __future__ import print_function, division
In [2]:
# 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
from time import time
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 [3]:
import targeting
import casjobs
import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table
from astropy.io import fits
from astropy.utils.console import ProgressBar
from astropy.utils import data
In [4]:
%matplotlib inline
from matplotlib import style, pyplot as plt
plt.style.use('seaborn-deep')
plt.rcParams['image.cmap'] = 'viridis'
plt.rcParams['image.origin'] = 'lower'
plt.rcParams['figure.figsize'] = (14, 8)
plt.rcParams['axes.titlesize'] = plt.rcParams['axes.labelsize'] = 16
plt.rcParams['xtick.labelsize'] = plt.rcParams['ytick.labelsize'] = 14
In [5]:
from IPython import display
In [6]:
# from the DECALS low-SB_brick selection and data download notebook
bricknames = ['1181m012', '2208m005']
In [7]:
catalog_fns = ['decals_dr3/catalogs/tractor-{}.fits'.format(bnm) for bnm in bricknames]
decals_catalogs = [Table.read(fn) for fn in catalog_fns]
decals_catalogs[0]
Out[7]:
In [8]:
sdss_fns = ['decals_dr3/catalogs/sdss_comparison_{}.csv'.format(bnm) for bnm in bricknames]
sdss_catalogs = [Table.read(fn) for fn in sdss_fns]
sdss_catalogs[0]
Out[8]:
In [9]:
bricks = Table.read('decals_dr3/survey-bricks.fits.gz')
bricksdr3 = Table.read('decals_dr3/survey-bricks-dr3.fits.gz')
In [10]:
A0p5 = 2.5*np.log10(np.pi*(0.5)**2)
A0p75 = 2.5*np.log10(np.pi*(0.75)**2)
A1 = 2.5*np.log10(np.pi*(1.)**2)
for dcat in decals_catalogs:
dcat['g'] = np.array(22.5 - 2.5*np.log10(dcat['decam_flux'][:, 1]))*u.mag
dcat['r'] = np.array(22.5 - 2.5*np.log10(dcat['decam_flux'][:, 2]))*u.mag
dcat['z'] = np.array(22.5 - 2.5*np.log10(dcat['decam_flux'][:, 4]))*u.mag
dcat['sb_r_0.5'] = np.array(22.5 - 2.5*np.log10(dcat['decam_apflux'][:, 2, 0]) + A0p5)*u.mag * u.arcsec**-2
dcat['sb_r_0.75'] = np.array(22.5 - 2.5*np.log10(dcat['decam_apflux'][:, 2, 1]) + A0p75)*u.mag * u.arcsec**-2
dcat['sb_r_1'] = np.array(22.5 - 2.5*np.log10(dcat['decam_apflux'][:, 2, 2]) + A1)*u.mag * u.arcsec**-2
In [11]:
dcat = decals_catalogs[0]
scat = sdss_catalogs[0]
bricknm = bricknames[0]
bricksdr3[bricknm==bricksdr3['brickname']]
Out[11]:
In [12]:
plt.scatter(dcat['ra'], dcat['dec'], lw=0, c='g')
plt.scatter(scat['RA'], scat['DEC'], lw=0, c='r')
Out[12]:
For unclear reasons, all of the SDSS fields near the deep-r area seem to be problematic somehow. Forging ahead, but probably not a good idea to dig into SDSS areal completeness too much here
In [13]:
plt.hist(scat['r'], bins=50, histtype='step', range=(15, 25), label='SDSS (Nx3)', weights=[3]*len(scat))
plt.hist(dcat['r'], bins=50, histtype='step', range=(15, 25), label='DECALS')
plt.legend(loc=0)
plt.xlabel('mag')
Out[13]:
Note that the relative numbers are boosted here for SDSS because of the areal coverage effect visible above. But it's clear that the magnitude limit is at least 3 deeper for DECALS.
In [14]:
points_depth_mag = -2.5*(np.log10(5*dcat['decam_depth']**-0.5)-9)
plt.hist(points_depth_mag[:, 2], bins=100, range=(22, 25.5), histtype='step', label='r')
plt.hist(points_depth_mag[:, 1], bins=100, range=(22, 25.5), histtype='step', label='g')
plt.xlabel(r'DECAM $5\sigma$ depth [mag]')
plt.legend()
Out[14]:
And this confirms that the typical depth for DECAM is ~2-3 mags below the SDSS
In [15]:
ssc = SkyCoord(scat['RA'], scat['DEC'], unit=u.deg)
dsc = SkyCoord(dcat['ra'], dcat['dec'], unit=u.deg)
idx, d2d, _ = dsc.match_to_catalog_sky(ssc)
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.hist(d2d.arcsec, bins=100, histtype='step', range=(0, 3), log=True)
ax2.hist(d2d.arcsec, bins=100, histtype='step', range=(0, 30), log=True)
idx, d2d, _ = ssc.match_to_catalog_sky(dsc)
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.hist(d2d.arcsec, bins=100, histtype='step', range=(0, 3), log=True)
ax2.hist(d2d.arcsec, bins=100, histtype='step', range=(0, 30), log=True)
None
OK, looks like all the w/i ~1" matches are real, and essentially everything has a match from the SDSS
In [16]:
dmatch = dcat[idx[d2d<1*u.arcsec]]
smatch = scat[d2d<1*u.arcsec]
plt.axhline(0, c='k')
plt.scatter(dmatch['r'], smatch['r']- dmatch['r'], lw=0)
plt.ylim(-1, 1)
plt.xlabel(r'$r_{\rm DECALS}$', fontsize=18)
plt.ylabel(r'$r_{\rm SDSS} - r_{\rm DECALS}$', fontsize=18)
plt.xlim(10, 24)
Out[16]:
OK, so there's a small flux offset, probably due to different flux definitions, but probably not a big systematic effect. Lets press on.
In [17]:
starmsk = smatch['PHOT_SG']=='STAR'
plt.scatter(smatch['r'][starmsk], smatch['SB_PETRO_R'][starmsk], lw=0, c='r', alpha=.5, label='star')
plt.scatter(smatch['r'][~starmsk], smatch['SB_PETRO_R'][~starmsk], lw=0, c='b', alpha=.5, label='gal')
plt.legend(loc='upper left')
plt.xlim(10, 24)
plt.ylim(15, 28)
plt.axvline(21, c='k')
plt.xlabel(r'$r_{\rm SDSS}$', fontsize=18)
plt.ylabel(r'$SB_{\rm SDSS}$', fontsize=18)
Out[17]:
The vertical line is the typical SAGA cutoff of r~21
In [18]:
fig, (ax1,ax2) = plt.subplots(2, 1)
starmsk_match = dmatch['type']=='PSF '
ax1.scatter(dmatch['r'][starmsk_match], dmatch['sb_r_0.5'][starmsk_match], lw=0, c='r', alpha=.5, label='star')
ax1.scatter(dmatch['r'][~starmsk_match], dmatch['sb_r_0.5'][~starmsk_match], lw=0, c='b', alpha=.5, label='gal')
ax1.legend(loc='lower right')
starmsk = dcat['type']=='PSF '
plt.scatter(dcat['r'][starmsk], dcat['sb_r_0.5'][starmsk], lw=0, c='r', alpha=.5)
plt.scatter(dcat['r'][~starmsk], dcat['sb_r_0.5'][~starmsk], lw=0, c='b', alpha=.5)
for ax in (ax1, ax2):
ax.set_xlim(13, 25)
ax.set_ylim(15, 28)
ax.axvline(21, c='k')
ax.axhline(24.5, c='k', ls=':')
ax.set_xlabel(r'$r_{\rm DECALS}$', fontsize=18)
ax1.set_ylabel(r'$SB_{\rm DECALS}$, in SDSS', fontsize=18)
ax2.set_ylabel(r'$SB_{\rm DECALS}$, all', fontsize=18)
Out[18]:
The horizontal line here is to guide the eye: it's a limit above which there's nothing in the SDSS that DECALS says is that low SB (at least to our r<21 cutoff).
But there are some things at very low DECAM SB that aren't showing up in the SDSS X-match... lets see what they are.
In [19]:
to_check = (dcat['r']<21)&(dcat['sb_r_0.5']>24.5)
de_cutout_url = 'http://legacysurvey.org/viewer/jpeg-cutout/?ra={0.ra.deg}&dec={0.dec.deg}&layer=decals-dr3&pixscale=0.1&bands=grz'
sd_cutout_url = 'http://legacysurvey.org/viewer/jpeg-cutout/?ra={0.ra.deg}&dec={0.dec.deg}&layer=sdssco&pixscale=0.1&bands=gri'
tabrows = []
for row in dcat[to_check]:
sc = SkyCoord(row['ra'], row['dec'], unit=u.deg)
objstr = '{}_{}<br>RA={:.4f}<br>Dec={:.4f}<br>r={:.2f}<br>sb={:.2f}'.format(row['brickid'], row['objid'], row['ra'], row['dec'], row['r'], row['sb_r_0.5'])
deimg = '<img src="{}">'.format(de_cutout_url.format(sc))
sdimg = '<img src="{}">'.format(sd_cutout_url.format(sc))
tabrows.append('<tr><td>{}</td><td>{}</td><td>{}</td></tr>'.format(objstr, deimg, sdimg))
display.HTML("""
<table>
<tr>
<th>obj</th>
<th>DECALS</th>
<th>SDSS</th>
</tr>
{}
</table>
""".format('\n'.join(tabrows)))
Out[19]:
Aha - looks like they are all tractor/catalog failures rather than real low-SB objects.
The lone suspicious ones are:
In [20]:
dcat = decals_catalogs[1]
scat = sdss_catalogs[1]
bricknm = bricknames[1]
bricksdr3[bricknm==bricksdr3['brickname']]
Out[20]:
In [21]:
plt.scatter(dcat['ra'], dcat['dec'], lw=0, c='g')
plt.scatter(scat['RA'], scat['DEC'], lw=0, c='r')
Out[21]:
For unclear reasons, all of the SDSS fields near the deep-r area seem to be problematic somehow. Forging ahead, but probably not a good idea to dig into SDSS areal completeness too much here
In [22]:
plt.hist(scat['r'], bins=50, histtype='step', range=(15, 25), label='SDSS', weights=[1]*len(scat))
plt.hist(dcat['r'], bins=50, histtype='step', range=(15, 25), label='DECALS')
plt.legend(loc=0)
plt.xlabel('mag')
Out[22]:
Unlike deep_r, here there is no areal twiddle-factor because the coverage is almost the same (there's a tiny bit more in DECALS for unclear reasons, but that's mostly in the noise). That confirms the result before that the magnitude limit is quite a bit deeper for DECALS.
In [23]:
points_depth_mag = -2.5*(np.log10(5*dcat['decam_depth']**-0.5)-9)
plt.hist(points_depth_mag[:, 2], bins=100, range=(22, 25.5), histtype='step', label='r')
plt.hist(points_depth_mag[:, 1], bins=100, range=(22, 25.5), histtype='step', label='g')
plt.xlabel(r'DECAM $5\sigma$ depth [mag]')
plt.legend()
Out[23]:
And this confirms that the typical depth for DECAM is ~1-2 mags below the SDSS
In [24]:
ssc = SkyCoord(scat['RA'], scat['DEC'], unit=u.deg)
dsc = SkyCoord(dcat['ra'], dcat['dec'], unit=u.deg)
idx, d2d, _ = dsc.match_to_catalog_sky(ssc)
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.hist(d2d.arcsec, bins=100, histtype='step', range=(0, 3), log=True)
ax2.hist(d2d.arcsec, bins=100, histtype='step', range=(0, 30), log=True)
idx, d2d, _ = ssc.match_to_catalog_sky(dsc)
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.hist(d2d.arcsec, bins=100, histtype='step', range=(0, 3), log=True)
ax2.hist(d2d.arcsec, bins=100, histtype='step', range=(0, 30), log=True)
None
OK, looks like all the w/i ~1" matches are real, and essentially everything has a match from the SDSS
In [25]:
dmatch = dcat[idx[d2d<1*u.arcsec]]
smatch = scat[d2d<1*u.arcsec]
plt.axhline(0, c='k')
plt.scatter(dmatch['r'], smatch['r']- dmatch['r'], lw=0)
plt.ylim(-1, 1)
plt.xlabel(r'$r_{\rm DECALS}$', fontsize=18)
plt.ylabel(r'$r_{\rm SDSS} - r_{\rm DECALS}$', fontsize=18)
plt.xlim(10, 24)
Out[25]:
OK, so there's a small flux offset, probably due to different flux definitions, but probably not a big systematic effect. Lets press on.
In [26]:
starmsk = smatch['PHOT_SG']=='STAR'
plt.scatter(smatch['r'][starmsk], smatch['SB_PETRO_R'][starmsk], lw=0, c='r', alpha=.5, label='star')
plt.scatter(smatch['r'][~starmsk], smatch['SB_PETRO_R'][~starmsk], lw=0, c='b', alpha=.5, label='gal')
plt.legend(loc='upper left')
plt.xlim(10, 24)
plt.ylim(15, 28)
plt.axvline(21, c='k')
plt.xlabel(r'$r_{\rm SDSS}$', fontsize=18)
plt.ylabel(r'$SB_{\rm SDSS}$', fontsize=18)
Out[26]:
The vertical line is the typical SAGA cutoff of r~21
In [27]:
fig, (ax1,ax2) = plt.subplots(2, 1)
starmsk_match = dmatch['type']=='PSF '
ax1.scatter(dmatch['r'][starmsk_match], dmatch['sb_r_0.5'][starmsk_match], lw=0, c='r', alpha=.5, label='star')
ax1.scatter(dmatch['r'][~starmsk_match], dmatch['sb_r_0.5'][~starmsk_match], lw=0, c='b', alpha=.5, label='gal')
ax1.legend(loc='lower right')
starmsk = dcat['type']=='PSF '
plt.scatter(dcat['r'][starmsk], dcat['sb_r_0.5'][starmsk], lw=0, c='r', alpha=.5)
plt.scatter(dcat['r'][~starmsk], dcat['sb_r_0.5'][~starmsk], lw=0, c='b', alpha=.5)
for ax in (ax1, ax2):
ax.set_xlim(13, 25)
ax.set_ylim(15, 28)
ax.axvline(21, c='k')
ax.axhline(24.5, c='k', ls=':')
ax.set_xlabel(r'$r_{\rm DECALS}$', fontsize=18)
ax1.set_ylabel(r'$SB_{\rm DECALS}$, in SDSS', fontsize=18)
ax2.set_ylabel(r'$SB_{\rm DECALS}$, all', fontsize=18)
Out[27]:
The horizontal line here is to guide the eye: it's a limit above which there's nothing in the SDSS that DECALS says is that low SB (at least to our r<21 cutoff).
But there are some things at very low DECAM SB that aren't showing up in the SDSS X-match... lets see what they are.
In [28]:
to_check = (dcat['r']<21)&(dcat['sb_r_0.5']>24.5)
de_cutout_url = 'http://legacysurvey.org/viewer/jpeg-cutout/?ra={0.ra.deg}&dec={0.dec.deg}&layer=decals-dr3&pixscale=0.1&bands=grz'
sd_cutout_url = 'http://legacysurvey.org/viewer/jpeg-cutout/?ra={0.ra.deg}&dec={0.dec.deg}&layer=sdssco&pixscale=0.1&bands=gri'
tabrows = []
for row in dcat[to_check]:
sc = SkyCoord(row['ra'], row['dec'], unit=u.deg)
objstr = '{}_{}<br>RA={:.4f}<br>Dec={:.4f}<br>r={:.2f}<br>sb={:.2f}'.format(row['brickid'], row['objid'], row['ra'], row['dec'], row['r'], row['sb_r_0.5'])
deimg = '<img src="{}">'.format(de_cutout_url.format(sc))
sdimg = '<img src="{}">'.format(sd_cutout_url.format(sc))
tabrows.append('<tr><td>{}</td><td>{}</td><td>{}</td></tr>'.format(objstr, deimg, sdimg))
display.HTML("""
<table>
<tr>
<th>obj</th>
<th>DECALS</th>
<th>SDSS</th>
</tr>
{}
</table>
""".format('\n'.join(tabrows)))
Out[28]:
One possibly real object: 328371_785. But it's suspiciously near a bright star...
In [29]:
# this yields [0.5,0.75,1.0,1.5,2.0,3.5,5.0,7.0] arcsec aperture mags
-2.5*np.log10(dcat[dcat['objid']==785]['decam_apflux'][0][2])+22.5
Out[29]:
Even being really conservative, it's r>21
In [30]:
for scat, bnm in zip(sdss_catalogs, bricknames):
plt.hist(scat['EXTINCTION_R'], histtype='step', label=bnm)
plt.legend(loc=0)
None
It looks like there are not low surface brightness objects that DECALS sees but SDSS does not.
While the DECALS catalog does show a few such objects, they're almost entirely (maybe entirely entirely?) artifacts from nearby bright stars.