The purpose of this notebook is to build a group catalog (using a simple friends-of-friends algorithm) from a diameter-limited (D25>5 arcsec) parent sample of galaxies defined and documented as part of the Legacy Survey Large Galaxy Atlas.
In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt
In [2]:
import astropy.units as u
from astropy.table import Table, Column
from astropy.coordinates import SkyCoord
In [3]:
import fitsio
from pydl.pydlutils.spheregroup import spheregroup
In [4]:
%matplotlib inline
In [5]:
LSLGAdir = os.getenv('LSLGA_DIR')
In [6]:
mindiameter = 0.25 # [arcmin]
linking_length = 2.5 # [arcmin]
In [7]:
suffix = '0.05'
In [8]:
ledafile = os.path.join(LSLGAdir, 'sample', 'leda-logd25-{}.fits'.format(suffix))
leda = Table.read(ledafile)
keep = (np.char.strip(leda['OBJTYPE']) != 'g') * (leda['D25'] / 60 > mindiameter)
leda = leda[keep]
keep = ['SDSS' not in gg and '2MAS' not in gg for gg in leda['GALAXY']]
#keep = np.logical_and( (np.char.strip(leda['OBJTYPE']) != 'g'), ~((leda['D25'] / 60 > 2.5) * (leda['BMAG'] > 16)) )
leda = leda[keep]
leda
Out[8]:
In [9]:
fig, ax = plt.subplots()
ax.scatter(leda['RA'], leda['DEC'], s=1, alpha=0.5)
Out[9]:
In [10]:
fig, ax = plt.subplots()
ax.hexbin(leda['BMAG'], leda['D25'] / 60, extent=(5, 20, 0, 20),
mincnt=1)
ax.set_xlabel('B mag')
ax.set_ylabel('D(25) (arcmin)')
Out[10]:
In [11]:
if False:
these = (leda['RA'] > 200) * (leda['RA'] < 210) * (leda['DEC'] > 5) * (leda['DEC'] < 10.0)
leda = leda[these]
print(np.sum(these))
In [12]:
%time grp, mult, frst, nxt = spheregroup(leda['RA'], leda['DEC'], linking_length / 60.0)
In [13]:
npergrp, _ = np.histogram(grp, bins=len(grp), range=(0, len(grp)))
nbiggrp = np.sum(npergrp > 1).astype('int')
nsmallgrp = np.sum(npergrp == 1).astype('int')
ngrp = nbiggrp + nsmallgrp
In [14]:
print('Found {} total groups, including:'.format(ngrp))
print(' {} groups with 1 member'.format(nsmallgrp))
print(' {} groups with 2-5 members'.format(np.sum( (npergrp > 1)*(npergrp <= 5) ).astype('int')))
print(' {} groups with 5-10 members'.format(np.sum( (npergrp > 5)*(npergrp <= 10) ).astype('int')))
print(' {} groups with >10 members'.format(np.sum( (npergrp > 10) ).astype('int')))
In [15]:
groupcat = Table()
groupcat.add_column(Column(name='GROUPID', dtype='i4', length=ngrp, data=np.arange(ngrp))) # unique ID number
groupcat.add_column(Column(name='GALAXY', dtype='S1000', length=ngrp))
groupcat.add_column(Column(name='NMEMBERS', dtype='i4', length=ngrp))
groupcat.add_column(Column(name='RA', dtype='f8', length=ngrp)) # average RA
groupcat.add_column(Column(name='DEC', dtype='f8', length=ngrp)) # average Dec
groupcat.add_column(Column(name='DIAMETER', dtype='f4', length=ngrp))
groupcat.add_column(Column(name='D25MAX', dtype='f4', length=ngrp))
groupcat.add_column(Column(name='D25MIN', dtype='f4', length=ngrp))
In [16]:
leda_groupid = leda.copy()
leda_groupid.add_column(Column(name='GROUPID', dtype='i4', length=len(leda)))
leda_groupid
Out[16]:
In [17]:
smallindx = np.arange(nsmallgrp)
In [18]:
ledaindx = np.where(npergrp == 1)[0]
groupcat['RA'][smallindx] = leda['RA'][ledaindx]
groupcat['DEC'][smallindx] = leda['DEC'][ledaindx]
groupcat['NMEMBERS'][smallindx] = 1
groupcat['GALAXY'][smallindx] = np.char.strip(leda['GALAXY'][ledaindx])
groupcat['DIAMETER'][smallindx] = leda['D25'][ledaindx] # [arcsec]
groupcat['D25MAX'][smallindx] = leda['D25'][ledaindx] # [arcsec]
groupcat['D25MIN'][smallindx] = leda['D25'][ledaindx] # [arcsec]
leda_groupid['GROUPID'][ledaindx] = groupcat['GROUPID'][smallindx]
In [19]:
bigindx = np.arange(nbiggrp) + nsmallgrp
In [20]:
coord = SkyCoord(ra=leda['RA']*u.degree, dec=leda['DEC']*u.degree)
In [21]:
def biggroups():
for grpindx, indx in zip(bigindx, np.where(npergrp > 1)[0]):
ledaindx = np.where(grp == indx)[0]
_ra, _dec = np.mean(leda['RA'][ledaindx]), np.mean(leda['DEC'][ledaindx])
d25min, d25max = np.min(leda['D25'][ledaindx]), np.max(leda['D25'][ledaindx])
groupcat['RA'][grpindx] = _ra
groupcat['DEC'][grpindx] = _dec
groupcat['D25MAX'][grpindx] = d25max
groupcat['D25MIN'][grpindx] = d25min
groupcat['NMEMBERS'][grpindx] = len(ledaindx)
groupcat['GALAXY'][grpindx] = ','.join(np.char.strip(leda['GALAXY'][ledaindx]))
leda_groupid['GROUPID'][ledaindx] = groupcat['GROUPID'][grpindx]
# Get the distance of each object from the group center.
cc = SkyCoord(ra=_ra*u.degree, dec=_dec*u.degree)
diameter = 2 * coord[ledaindx].separation(cc).arcsec.max()
groupcat['DIAMETER'][grpindx] = np.max( (diameter*1.02, d25max) )
In [22]:
%time biggroups()
In [23]:
leda_groupid
Out[23]:
In [24]:
groupcat
Out[24]:
In [25]:
ww = np.where(groupcat['NMEMBERS'] >= 2)[0]
fig, ax = plt.subplots()
ax.scatter(groupcat['RA'][ww], groupcat['DEC'][ww], s=1, alpha=0.5)
Out[25]:
In [26]:
groupfile = os.path.join(LSLGAdir, 'sample', 'leda-logd25-{}-groupcat.fits'.format(suffix))
print('Writing {}'.format(groupfile))
groupcat.write(groupfile, overwrite=True)
In [27]:
ledafile_groupid = os.path.join(LSLGAdir, 'sample', 'leda-logd25-{}-groupid.fits'.format(suffix))
print('Writing {}'.format(ledafile_groupid))
leda_groupid.write(ledafile_groupid, overwrite=True)