<jupyter_start><jupyter_text>Identifying junk in the large-galaxy atlasThe goal of this notebook is to investigate the junk in the parent (Hyperleda) large-galaxy sample, so it can be culled (in `LSLGA-build-parent`).<jupyter_code>import os
import numpy as np
import seaborn as sns
sns.set(style='ticks', font_scale=1.8, palette='Set2')
import matplotlib.pyplot as plt
from astrometry.libkd.spherematch import tree_build_radec, tree_search_radec
import LSLGA.io
import LSLGA.misc
%matplotlib inline
pp = LSLGA.io.read_hyperleda(verbose=True, version='v2.0')
pp = pp[np.isfinite(pp['MAG'])]
pp<jupyter_output>Read 1436176 objects from /Users/ioannis/research/projects/LSLGA/sample/v2.0/hyperleda-d25min10-18nov14.fits
Read 1436176 objects from /Users/ioannis/research/projects/LSLGA/sample/v2.0/hyperleda-d25min10-18nov14-allwise.fits
  Identified 1379531/1436176 (96.06%) objects with AllWISE photometry.<jupyter_text>Compare optical and AllWISE diametersTo do this comparison we require W1RSEMI to be measured, which is only done if the WISE source is within 2 arcsec of a 2MASS/XSC galaxy.  See:* http://wise2.ipac.caltech.edu/docs/release/allwise/expsup/sec5_3bvii.htmlxsc* http://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec4_4c.htmlxsc<jupyter_code>w1 = np.where( (pp['CNTR'] != -1) )[0]
w2 = np.where( np.isfinite(pp['W1RSEMI'][w1]) )[0]
print('{:.2f}% galaxies in AllWISE, {:.2f}% with measured WISE/2MASS diameters.'.format(
    100*len(w1)/len(pp), 100*len(w2)/len(w1)))<jupyter_output>96.58% galaxies in AllWISE, 24.16% with measured WISE/2MASS diameters.<jupyter_text>This agreement isn't crazy. Basically the D(25) diameters are much larger than the isophotal IR diameters.  However, we should probably look into the sources that are small in the optical and big in the IR.<jupyter_code>fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(pp['D25'][w1][w2], pp['W1RSEMI'][w1][w2]/60, s=1)
ax.plot([-1, 5], [-1, 5], lw=2, ls='-', color='k', alpha=0.7)
ax.set_xlim(0, 3)
ax.set_ylim(0, 3)<jupyter_output><empty_output><jupyter_text>Find galaxies that are large in the optical but unresolved in AllWISE/2MASSWe know the catalog contains spuriously large "galaxies", e.g., due to filter reflections in the SDSS imaging.  The strategy here is to find "large" objects that WISE says are unresolved.<jupyter_code>j1 = np.where( (pp['CNTR'] != -1) )[0]
#j2 = np.where( np.isnan(pp['W1RSEMI'][j1]) )[0] # this depends on being close to a 2MASS source
j2 = np.where( pp['EXT_FLG'][j1] == 0 )[0] # this depends on the WISE PSF only
print('{:.2f}% galaxies in AllWISE, {:.2f}% are unresolved.'.format(100*len(j1)/len(pp), 100*len(j2)/len(j1)))
fig, ax = plt.subplots(figsize=(8, 6))
_ = ax.hist(pp['D25'][j1][j2], bins=100, range=(0, 2), log=True)
            #cumulative=True, histtype='step', lw=2, density=True)
ax.axhline(y=20, lw=2, color='k')
ax.set_xlabel('D(25) (arcmin)')
_ = ax.set_title('Galaxies unresolved by AllWISE')
j3 = np.where( (pp['D25'][j1][j2] > 1.0) )[0]
#pp['GALAXY', 'PGC', 'RA', 'DEC', 'TYPE', 'D25', 'DIAM_REF', 'MAG', 'MAG_REF', 'SB_D25'][j1][j2][j3]
fig, ax = plt.subplots(figsize=(8, 6))
_ = ax.hist(pp['SB_D25'][j1][j2][j3], bins=100, log=True)#, range=(0, 3), log=True)
ax.axvline(x=27, lw=2, color='k')
ax.set_xlabel('Mean Surface Brightness (mag/arcsec2)')
_ = ax.set_title('Galaxies unresolved by AllWISE and D(25)>1 arcmin')<jupyter_output><empty_output><jupyter_text>The figure below is interesting:1. Galaxies with D(25)<10 arcsec are almost exclusively from the SDSS: they have half-light radii but not isophotal radii. We will likely reject these from the large-galaxy sample because Tractor does reasonably well on them.2. There's a spur of galaxies with D(25)~0.98 arcmin. Based on previous visual inspection, these are almost all from the SDSS and are cases where the Petrosian fitting failed (e.g., because of a filter reflection).  Reject these.3. I think all the galaxies with D(25)>0.95 arcmin and SB>27 mag/arcsec2 are also spurious or at least suspect.  This cut will likely reject bona fide low surface brightness dwarf galaxies, but I don't think we want them in the sample anyway (e.g., as large-galaxy priors).<jupyter_code>def plotbox(ax):
    ax.plot([0.95, 1000], [27, 27], color='k', lw=1)
    ax.plot([0.95, 0.95], [27, 40], color='k', lw=1)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))#, sharey=True)
if False:
    ax1.scatter(pp['D25'][j1][j2], pp['SB_D25'][j1][j2], s=2)
    ax2.scatter(pp['D25'][j1][j2], pp['SB_D25'][j1][j2], s=2)
else:
    ax1.scatter(pp['D25'], pp['SB_D25'], s=2)
    ax2.scatter(pp['D25'], pp['SB_D25'], s=2)
ax1.set_xlim(0.05, 10)
ax1.set_ylim(18, 36)
ax1.set_xscale('log')
ax1.set_ylabel('Mean Surface Brightness (mag/arcsec2)')
for ax in (ax1, ax2):
    ax.axvline(x=10/60, lw=2, ls='--', color='k')
    ax.set_xlabel('D(25) (arcmin)')
    plotbox(ax)

ax2.set_xlim(0.5, 5)
ax2.set_ylim(25, 35)
if False:
    fig.suptitle('Galaxies unresolved by AllWISE')
else:
    fig.suptitle('All galaxies')

fig.subplots_adjust(wspace=0.1)<jupyter_output><empty_output><jupyter_text>Here are the galaxies culled by our $D(25)>0.95$ arcmin & $SB_{25}>27$ mag/arcsec2 cut:Notice the sample is only about 500 objects.<jupyter_code>cut = np.where( (pp['D25'] > 0.95) * (pp['SB_D25'] > 27) )[0]
#print(pp[cut]['GALAXY'].data)
pp['GALAXY', 'PGC', 'RA', 'DEC', 'TYPE', 'D25', 'DIAM_REF', 'MAG', 'MAG_REF', 'SB_D25', 'EXT_FLG', 'IN_ALLWISE'][cut]<jupyter_output><empty_output><jupyter_text>Plot the cumulative fraction of galaxies above various angular diameter cuts.<jupyter_code>keep = np.delete( np.arange(len(pp)), cut )
#srt = np.argsort(pp['D25'][keep])
for d25cut in (0.5, 0.75, 1.0, 1.5):
    nn = np.sum( (pp['D25'][keep] > d25cut) )
    print('Number (fraction) of galaxies with D(25)>{:.2f} arcmin: {}/{} ({:.3f}%)'.format(
        d25cut, nn, len(keep), 100 * nn / len(keep) ))

fig, ax = plt.subplots(figsize=(8, 6))
_ = ax.hist(pp['D25'][keep], bins=100, range=(0, 2), #log=True,
            cumulative=True, histtype='step', lw=2, density=True)
ax.set_xlabel('D(25) (arcmin)')
_ = ax.set_title('All Galaxies')<jupyter_output>Number (fraction) of galaxies with D(25)>0.50 arcmin: 174901/1422394 (12.296%)
Number (fraction) of galaxies with D(25)>0.75 arcmin: 52007/1422394 (3.656%)
Number (fraction) of galaxies with D(25)>1.00 arcmin: 21182/1422394 (1.489%)
Number (fraction) of galaxies with D(25)>1.50 arcmin: 6156/1422394 (0.433%)<jupyter_text>Inspect all the "SDSS" and "2MASS" galaxies, especially the largest ones.<jupyter_code>tiles = LSLGA.io.read_desi_tiles()
%time indesi = LSLGA.misc.is_point_in_desi(tiles, pp['RA'], pp['DEC']).astype(bool)
lsbcut = (pp['SB_D25'] > 27) * (pp['D25'] > 0.95)
diamcut = pp['D25'] > 20/60
sample = pp[~lsbcut * diamcut * indesi]
print(len(sample))
%time kdsample = tree_build_radec(sample['RA'], sample['DEC'])
issdss = np.array(['SDSS' in gg.decode('utf-8') for gg in sample['GALAXY'].data])
is2mass = np.array(['MASS' in gg.decode('utf-8') for gg in sample['GALAXY'].data])
print(np.sum(issdss), np.sum(is2mass))<jupyter_output>30221 4<jupyter_text>Here it's likely that many (all?) of the largest "SDSS" galaxies are spurious; check them.<jupyter_code>rest = np.logical_or( ~issdss, ~is2mass )
fig, ax = plt.subplots(figsize=(8, 6))
_ = ax.hist(sample['D25'][rest], bins=100, range=(0, 5), log=True, label='All')
_ = ax.hist(sample['D25'][issdss], bins=100, range=(0, 5), log=True, label='SDSS')
_ = ax.hist(sample['D25'][is2mass], bins=100, range=(0, 5), log=True, label='2MASS')
ax.legend(loc='upper right')
#ax.axhline(y=20, lw=2, color='k')
_ = ax.set_xlabel('D(25) (arcmin)')
_ = ax.set_title('SDSS and 2MASS-named Galaxies')
ww = issdss * (sample['D25'] > 1)
print(np.sum(ww))
sdssgal = [gg.decode('utf-8') for gg in sample['GALAXY'][ww].data]
#print(sdssgal)
out = sample['RA', 'DEC', 'GALAXY'][ww]
out.rename_column('GALAXY', 'NAME')
#out.write('/Users/ioannis/Downloads/inspect.fits', overwrite=True)<jupyter_output><empty_output><jupyter_text>Flag all the SDSS galaxies inside other large galaxies.<jupyter_code>verb = False
galingal = np.zeros( len(sample), dtype=bool)
for ss in sample:
    I = tree_search_radec(kdsample, ss['RA'], ss['DEC'], ss['D25'] / 2 / 60)
    if len(I) > 1:
        isdss = np.array(['SDSS' in gg for gg in sample['GALAXY'][I]])
        if np.sum(isdss) > 0:
            galingal[I[isdss]] = True
            if verb:
                print(sample[I]['GALAXY', 'RA', 'DEC', 'D25'])
                print('http://legacysurvey.org/viewer-dev?ra={:.5f}&dec={:.4f}&zoom=15&layer=ls-dr67&lslga'.format(
                    sample[I[0]]['RA'], sample[I[0]]['DEC']))
                import pdb ; pdb.set_trace()
print('  Found {}/{} ({:.2f}%) SDSS galaxies inside another large galaxy.'.format(
    np.sum(galingal), len(sample), 100*np.sum(galingal)/len(sample)))
sample[galingal]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1814702240234155, 'fraction_numerical': 0.061127997298210064, 'mean_word_length': 5.287331917905166, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '422140', 'license': 'MIT', 'n_tokens_mistral': 3883, 'n_tokens_neox': 3546, 'n_words': 832, 'path': 'py/SGA/sandbox/sample-find-junk.ipynb', 'repo_name': 'moustakas/SGA', 'size': '8883'}