In [1]:
import os, pdb
import numpy as np
from glob import glob
import fitsio
from astropy.table import vstack, Table, hstack
In [2]:
import seaborn as sns
sns.set(context='talk', style='ticks', font_scale=1.6)
%matplotlib inline
In [4]:
dr8dir = '/global/project/projectdirs/cosmo/data/legacysurvey/dr8'
outdir = '/global/project/projectdirs/desi/users/ioannis/dr8-lslga'
lslgafile = '/global/project/projectdirs/cosmo/staging/largegalaxies/v2.0/LSLGA-v2.0.fits'
In [5]:
lslga = Table(fitsio.read(lslgafile))
#lslga
In [19]:
def build_lslga_dr8(northsouth, overwrite=False):
cols = ['RA', 'DEC', 'TYPE', 'BRICKNAME', 'REF_CAT',
'FRACDEV', 'SHAPEEXP_R', 'SHAPEEXP_E1', 'SHAPEEXP_E2',
'SHAPEDEV_R', 'SHAPEDEV_E1', 'SHAPEDEV_E2', 'REF_ID']
outfile = os.path.join(outdir, 'dr8-lslga-{}.fits'.format(northsouth))
if not os.path.isfile(outfile) or overwrite:
out = []
catfile = glob(os.path.join(dr8dir, northsouth, 'sweep', '8.0', 'sweep-*.fits'))
#catfile = glob(os.path.join(dr8dir, northsouth, 'tractor', '???', 'tractor*.fits'))
for ii, ff in enumerate(catfile):#[:3]):
if ii % 20 == 0:
print('{} / {}'.format(ii, len(catfile)))
cc = fitsio.read(ff, columns=cols, upper=True)
ckeep = np.where(cc['REF_CAT'] == b'L2')[0]
out.append(Table(cc[ckeep]))
out = vstack(out)
# find and remove duplicates
#uval, cnt = np.unique(out['REF_ID'], return_counts=True)
#dup_refid = uval[cnt > 1]
# now match against the LSLGA
match = np.hstack([np.where(refid == lslga['LSLGA_ID'])[0] for refid in out['REF_ID']])
lss = lslga[match]
lss.rename_column('RA', 'RA_LSLGA')
lss.rename_column('DEC', 'DEC_LSLGA')
lss.rename_column('TYPE', 'MORPHTYPE')
out = hstack((lss, out))
print('Writing {}'.format(outfile))
out.write(outfile, overwrite=True)
else:
print('Reading {}'.format(outfile))
out = Table(fitsio.read(outfile))
return out
In [20]:
%time north = build_lslga_dr8('north', overwrite=False)
In [22]:
%time south = build_lslga_dr8('south', overwrite=False)
In [47]:
def merge_north_south(overwrite=False):
mergedfile = os.path.join(outdir, 'dr8-lslga-northsouth.fits')
if not os.path.isfile(mergedfile) or overwrite:
out = []
for northsouth in ('north', 'south'):
catfile = os.path.join(outdir, 'dr8-lslga-{}.fits'.format(northsouth))
cat = Table(fitsio.read(catfile))
cat['REGION'] = northsouth
out.append(cat)
out = vstack(out)
print('Writing {}'.format(mergedfile))
out.write(mergedfile, overwrite=True)
else:
print('Reading {}'.format(mergedfile))
out = Table(fitsio.read(mergedfile))
return out
In [48]:
%time cat = merge_north_south(overwrite=True)
In [49]:
cat
Out[49]:
In [50]:
urefid, uindx = np.unique(cat['REF_ID'], return_index=True)
dup = np.delete(np.arange(len(cat)), uindx)
cat[dup]
cat[cat['REF_ID'] == 462171]
Out[50]:
In [58]:
ww = [tt.strip() == 'REX' for tt in cat['TYPE']]
cat[ww]
Out[58]:
In [73]:
from astrometry.util.fits import fits_table
T = fits_table(os.path.join(outdir, 'dr8-lslga-northsouth.fits'))
In [96]:
def get_ba_pa(e1, e2):
ee = np.hypot(e1, e2)
ba = (1 - ee) / (1 + ee)
#pa = -np.rad2deg(np.arctan2(e2, e1) / 2)
pa = 180 - (-np.rad2deg(np.arctan2(e2, e1) / 2))
return ba, pa
def lslga_model(T):
ba, pa = [], []
for Tone in T:
ttype = Tone.type.strip()
if ttype == 'DEV' or ttype == 'COMP':
e1, e2 = Tone.shapedev_e1, Tone.shapedev_e2
_ba, _pa = get_ba_pa(e1, e2)
elif ttype == 'EXP' or ttype == 'REX':
e1, e2 = Tone.shapeexp_e1, Tone.shapeexp_e2
_ba, _pa = get_ba_pa(e1, e2)
else: # PSF
_ba, _pa = np.nan, np.nan
ba.append(_ba)
pa.append(_pa)
T.ba_model = np.hstack(ba)
T.pa_model = np.hstack(pa)
T.radius_model_arcsec = T.fracdev * T.shapedev_r + (1 - T.fracdev) * T.shapeexp_r
T.radius_arcsec = T.d25 / 2. * 60.
return T
In [97]:
rr = lslga_model(T[:10])
In [98]:
rr.pa, rr.pa_model
rr.ba, rr.ba_model
rr.radius_arcsec, rr.radius_arcsec_model
Out[98]:
In [ ]: