In [3]:
import os, time
import fitsio
import numpy as np
import matplotlib.pyplot as plt
from glob import glob
from astropy.table import vstack, Table, hstack
In [4]:
import seaborn as sns
sns.set(context='talk', style='ticks', font_scale=1.6)
%matplotlib inline
In [5]:
outdir = '/Users/ioannis/work/legacysurvey/dr8'
#outdir = '/global/project/projectdirs/desi/users/ioannis/dr8-lslga'
In [6]:
#lslgadir = '/global/project/projectdirs/cosmo/staging/largegalaxies/v2.0'
lslgadir = '/Users/ioannis/research/projects/LSLGA/sample/v2.0'
lslgafile = os.path.join(lslgadir, 'LSLGA-v2.0.fits')
lslga = Table.read(lslgafile)
print('Read {} galaxies from {}'.format(len(lslga), lslgafile))
In [7]:
cut = lslga['IN_DESI'] * (lslga['D25'] > 1.5)
len(lslga[cut]) * 7 / 14000
Out[7]:
In [123]:
MASKBITS = dict(
NPRIMARY = 0x1, # not PRIMARY
BRIGHT = 0x2,
SATUR_G = 0x4,
SATUR_R = 0x8,
SATUR_Z = 0x10,
ALLMASK_G = 0x20,
ALLMASK_R = 0x40,
ALLMASK_Z = 0x80,
WISEM1 = 0x100, # WISE masked
WISEM2 = 0x200,
BAILOUT = 0x400, # bailed out of processing
MEDIUM = 0x800, # medium-bright star
GALAXY = 0x1000, # LSLGA large galaxy
CLUSTER = 0x2000, # Cluster catalog source
)
In [148]:
randomfile = os.path.join(outdir, 'randoms-inside-dr8-0.31.0-1.fits')
ff = fitsio.FITS(randomfile)
nrows = ff[1].get_nrows()
nrand = 10000000
factor = nrows / nrand / 5000
%time rand = ff[1].read(rows=np.arange(nrand))
print('Read {}/{} random positions from {}'.format(len(rand), nrows, randomfile))
In [149]:
print(rand.dtype)
In [150]:
infoot = (rand['NOBS_G'] > 0) * (rand['NOBS_R'] > 0) * (rand['NOBS_Z'] > 0) * ((rand['MASKBITS'] & MASKBITS['NPRIMARY']) == 0)
In [151]:
totarea = np.sum(infoot) * factor
print('Total area = {:.3f}'.format(totarea))
for bit in MASKBITS.keys():
area = np.sum(infoot * ((rand['MASKBITS'] & MASKBITS[bit]) != 0) * factor)
frac = 100 * area / totarea
print(' {}: {:.3f} deg2, {:.3f}%'.format(bit, area, frac))
In [6]:
def read_dr8_lslga(region='north', clobber=False):
sweepdir = '/global/project/projectdirs/cosmo/work/legacysurvey/dr8/{}/sweep/8.0/'.format(region)
outfile = os.path.join(outdir, 'dr8-lslga-{}.fits'.format(region))
if os.path.isfile(outfile) and not clobber:
print('Reading {}'.format(outfile))
out = Table.read(outfile)
else:
out = []
catfile = glob(os.path.join(sweepdir, 'sweep*.fits'))
for ii, ff in enumerate(catfile):
if ii % 50 == 0:
print('{} / {}'.format(ii, len(catfile)))
cc = Table(fitsio.read(ff, columns=['REF_CAT', 'REF_ID', 'RA', 'DEC', 'TYPE', 'FRACDEV',
'SHAPEEXP_R', 'SHAPEEXP_E1', 'SHAPEEXP_E2',
'SHAPEDEV_R', 'SHAPEDEV_E1', 'SHAPEDEV_E2'], upper=True))
ckeep = np.where(cc['REF_CAT'] == 'L2')[0]
#keep = np.where([rcat.decode('utf-8').strip() == 'L2' for rcat in cc['REF_CAT']])[0]
if len(ckeep) > 0:
cc = cc[ckeep]
cc = cc[np.argsort(cc['REF_ID'])] # sort!
_, uu = np.unique(cc['REF_ID'], return_index=True)
if len(uu) != len(cc):
print('Duplicate large galaxy in {}'.format(os.path.basename(ff)))
cc = cc[uu]
lkeep = np.where(np.isin(lslga['LSLGA_ID'], cc['REF_ID']))[0]
if len(lkeep) != len(cc):
print('Still a problem with duplicates!')
lss = lslga[lkeep]
assert(np.all(lss['LSLGA_ID'] == cc['REF_ID']))
lss.rename_column('RA', 'RA_LSLGA')
lss.rename_column('DEC', 'DEC_LSLGA')
lss.rename_column('TYPE', 'MORPHTYPE')
out.append(hstack((lss, cc)))
#import pdb ; pdb.set_trace()
out = vstack(out)
out.write(outfile, overwrite=True)
return out
In [12]:
%time north = read_dr8_lslga(region='north', clobber=False)
%time south = read_dr8_lslga(region='south', clobber=False)
In [23]:
cat = vstack((north, south))
print('Read {} large galaxies.'.format(len(cat)))
cat[:5]
Out[23]:
In [15]:
def qa_radec(png=None):
dra, ddec = (cat['RA'] - cat['RA_LSLGA']).data * 3600, (cat['DEC'] - cat['DEC_LSLGA']).data * 3600
print(np.median(dra), np.std(dra), np.median(ddec), np.std(ddec))
fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(dra, ddec, s=50)
ax.set_xlim(-40, 40)
ax.set_ylim(-40, 40)
ax.axhline(y=0, ls='-', alpha=0.8, color='k')
ax.axvline(x=0, ls='-', alpha=0.8, color='k')
ax.set_xlabel(r'$\Delta\,$(RA) (arcsec)')
ax.set_ylabel(r'$\Delta\,$(Dec) (arcsec)')
plt.subplots_adjust(bottom=0.15, left=0.15)
if png:
pngfile = os.path.join(outdir, png)
print('Writing {}'.format(pngfile))
fig.savefig(pngfile)
In [16]:
qa_radec(png='qa-radec.png')
In [20]:
def get_e1e2(ba, phi):
ab = 1. / ba
e = (ab - 1) / (ab + 1)
ee = -np.log(1 - e)
angle = np.deg2rad(2. * (-phi))
ee1 = ee * np.cos(angle)
ee2 = ee * np.sin(angle)
return ee1, ee2
def type2properties(cat, objtype):
this = np.where([objtype == tt.strip() for tt in cat['TYPE']])[0]
if objtype == 'EXP' or objtype == 'REX':
reff, e1, e2 = cat['SHAPEEXP_R'][this], cat['SHAPEEXP_E1'][this], cat['SHAPEEXP_E2'][this]
elif objtype == 'DEV':
reff, e1, e2 = cat['SHAPEDEV_R'][this], cat['SHAPEDEV_E1'][this], cat['SHAPEDEV_E2'][this]
lslga_rad = cat['D25'][this] / 2 * 60 / 2 # [arcmin]
lslga_e1, lslga_e2 = get_e1e2(cat['BA'][this], 180-cat['PA'][this])
return reff, e1, e2, lslga_rad, lslga_e1, lslga_e2, cat[this]
In [21]:
def qa_morph(png=None):
fig, ax = plt.subplots(1, 3, figsize=(18, 5))
#for objtype in ('COMP', 'DEV', 'EXP', 'PSF', 'REX'):
for objtype, marker in zip(('DEV', 'EXP', 'REX'), ('s', 'o', '*')):
reff, e1, e2, lslga_rad, lslga_e1, lslga_e2, _ = type2properties(cat, objtype)
ax[0].scatter(reff, lslga_rad, label=objtype, s=20, alpha=0.7, marker=marker)
if objtype != 'REX':
ax[1].scatter(e1, lslga_e1, s=20, alpha=0.7, marker=marker)
ax[2].scatter(e2, lslga_e2, s=20, alpha=0.7, marker=marker)
ax[0].set_xlim(0, 40)
ax[0].set_ylim(0, 40)
ax[0].axhline(y=5, color='k')
ax[0].set_xlabel(r'r$_{eff}$ (Tractor, arcsec)')
ax[0].set_ylabel(r'0.5 * R(25) (LSLGA, arcsec)')
ax[1].set_xlabel(r'e$_1$ (Tractor)')
ax[1].set_ylabel(r'e$_1$ (LSLGA)')
ax[2].set_xlabel(r'e$_2$ (Tractor)')
ax[2].set_ylabel(r'e$_2$ (LSLGA)')
for xx in ax[1:]:
xx.set_xlim(-1, 1)
xx.set_ylim(-1, 1)
leg = ax[0].legend(loc='upper left', frameon=True, fontsize=16, markerscale=2)
for l in leg.get_lines():
l.set_alpha(1)
for xx in ax:
xx.plot(xx.get_xlim(), xx.get_ylim(), color='k', ls='--')
plt.subplots_adjust(wspace=0.35, bottom=0.2)
if png:
pngfile = os.path.join(outdir, png)
print('Writing {}'.format(pngfile))
fig.savefig(pngfile)
In [22]:
qa_morph(png='qa-lslga-morph.png')
In [95]:
def qa_maskbits(png=None):
from astropy.io import fits
from astropy.visualization import AsinhStretch as Stretch
from astropy.visualization import ImageNormalize
from astropy.visualization import ZScaleInterval as Interval
from astrometry.util.util import Tan
brick = '1940p412'
imgfile = os.path.join(outdir, 'legacysurvey-{}-image-r.fits.fz'.format(brick))
img = fitsio.read(imgfile)
msk = fitsio.read(os.path.join(outdir, 'legacysurvey-{}-maskbits.fits.fz'.format(brick)))
hdr = fits.getheader(imgfile, ext=1)
wcs = Tan(hdr['CRVAL1'], hdr['CRVAL2'], hdr['CRPIX1'], hdr['CRPIX2'], hdr['CD1_1'],
hdr['CD1_2'], hdr['CD2_1'], hdr['CD2_2'], hdr['NAXIS1'], hdr['NAXIS2'])
cmap = plt.cm.viridis
stretch = Stretch(a=0.9)
interval = Interval(contrast=0.5, nsamples=10000)
norm = ImageNormalize(img, interval=interval, stretch=stretch)
inchperax = 4
fig, ax = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(inchperax*3, 3))
ax[0].imshow(img, origin='lower', norm=norm, cmap=cmap, interpolation='none')
ax[1].imshow((msk & 2**12) != 0, origin='lower', interpolation='none')
ax[2].imshow((msk & 2**1) != 0, origin='lower', interpolation='none')
for xx in ax:
xx.get_xaxis().set_visible(False)
xx.get_yaxis().set_visible(False)
xx.set_aspect('equal')
xx.axis('off')
xx.autoscale(False)
fig.subplots_adjust(wspace=0)
plt.tight_layout(w_pad=0)
if png:
pngfile = os.path.join(outdir, png)
print('Writing {}'.format(pngfile))
fig.savefig(pngfile, bbox_inches='tight')#, pad_inches=0)
In [20]:
def qa_PGC3087924(png=None):
from astropy.io import fits
from astropy.visualization import AsinhStretch as Stretch
from astropy.visualization import ImageNormalize
from astropy.visualization import ZScaleInterval as Interval
from astrometry.util.util import Tan
brick = '1940p412'
imgfile = os.path.join(outdir, 'legacysurvey-{}-image-r.fits.fz'.format(brick))
img = fitsio.read(imgfile)
msk = fitsio.read(os.path.join(outdir, 'legacysurvey-{}-maskbits.fits.fz'.format(brick)))
hdr = fits.getheader(imgfile, ext=1)
wcs = Tan(hdr['CRVAL1'], hdr['CRVAL2'], hdr['CRPIX1'], hdr['CRPIX2'], hdr['CD1_1'],
hdr['CD1_2'], hdr['CD2_1'], hdr['CD2_2'], hdr['NAXIS1'], hdr['NAXIS2'])
cmap = plt.cm.viridis
stretch = Stretch(a=0.9)
interval = Interval(contrast=0.5, nsamples=10000)
norm = ImageNormalize(img, interval=interval, stretch=stretch)
rad = 250
_, x0, y0 = wcs.radec2pixelxy(194.0090, 41.2926)
x0, y0 = np.int(x0), np.int(y0)
imgcut = img[y0-rad:y0+rad, x0-rad:x0+rad]
mskcut = msk[y0-rad:y0+rad, x0-rad:x0+rad]
fig, ax = plt.subplots()
#ax.imshow(imgcut, origin='lower', norm=norm, cmap=cmap, interpolation='none')
ax.imshow(mskcut, origin='lower')
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
ax.set_aspect('equal')
ax.axis('off')
ax.autoscale(False)
fig.subplots_adjust(wspace=0)
plt.tight_layout(w_pad=0)
if png:
pngfile = os.path.join(outdir, png)
print('Writing {}'.format(pngfile))
fig.savefig(pngfile, bbox_inches='tight', pad_inches=0)
In [21]:
qa_PGC3087924(png='PGC3087924-maskbits.jpg')
In [96]:
qa_maskbits(png='qa-maskbits.png')
In [ ]:
stop
In [ ]:
reff, e1, e2, lslga_rad, lslga_e1, lslga_e2, devcat = type2properties(cat, 'DEV')
big = np.where(reff > 25)[0]
for cc in devcat[big][:20]:
size = np.round(cc['D25'] * 1.5 * 60 / 0.262).astype(int)
montagefile = 'montage/{}.jpg'.format(cc['GALAXY'].lower())
jpgfile = []
#print('http://legacysurvey.org/viewer-dev?ra={}&dec={}&zoom=14&layer=dr8b-decam&lslga'.format(cc['RA'], cc['DEC']))
for ii, imtype in enumerate(('', '-model', '-resid')):
jpgfile.append('jpg/{}{}.jpg'.format(cc['GALAXY'].lower(), imtype))
url = '"http://legacysurvey.org/viewer-dev/jpeg-cutout?ra={}&dec={}&size={}&layer=dr8b-decam{}"'.format(
cc['RA'], cc['DEC'], size, imtype)
if not os.path.exists(jpgfile[ii]):
cmd = 'wget --continue -O {} {}'.format(jpgfile[ii], url)
print(cmd)
os.system(cmd)
time.sleep(1)
cmd = 'montage -bordercolor white -borderwidth 1 -tile 3x1 -geometry +0+0 '
cmd = cmd+' '.join(ff for ff in jpgfile)
cmd = cmd+' {}'.format(montagefile)
print(cmd)
os.system(cmd)
In [ ]: