In [1]:
__author__ = 'Stephanie Juneau, NOAO Data Lab Team'
__version__ = '20170602' # yyyymmdd

Star/Galaxy/QSO Classification in DECaLS

by Stéphanie Juneau, Robert Nikutta, Knut Olsen and the NOAO Data Lab Team

In this notebook, we investigate the optical and infrared colors of astronomical sources detected in the DECam Legacy Survey (DECaLS), which comprises stars, galaxies and quasars (or QSOs, which stands for Quasi-Stellar Objects). We will plot color combinations, and take into account the source "type" as defined from the light profile shape.

DECaLS will cover ~9500 deg² in the g, r, z bands to depths of g=24.7, r=23.9, z=23.0. Infrared WISE data are also extracted at the location of the DECaLS sources.

The third data release (DR3) includes a portion of DECaLS, covering a disjoint footprint with 4300 deg² in g-band, 4600 deg² in r-band and 8100 deg² in z-band, of which 4200 deg² has been observed in all three optical filters.

You can read more about DECaLS on the Data Lab survey page (here) and on the Legacy Survey team website (here).

Import

In [2]:
# Import packages including some Data Lab (dl) specific ones.
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.ticker import NullFormatter
from astropy.utils.data import download_file  #import file from URL
from scipy.stats import binned_statistic_2d
%matplotlib inline

from cStringIO import StringIO  #C script to handle string format
from astropy.table import Table

from dl import authClient as ac, queryClient as qc
print 'Done importing'
Done importing

Authenticate

In [3]:
# To save to virtual space, you need to log in your account (not anonymous)
token = ac.login('anonymous')
print 'Done getting token'
Done getting token

Introduction

Ellipticals are red; Spirals are blue.

An imaging survey includes a zoo of different astronomical objects. There are foreground stars from our own Milky Way galaxies, and background galaxies at various distances, including QSOs (quasi-stellar objects) with actively accreting supermassive black holes.

**Figure 1.** Small section of the DECaLS DR3 image from the Legacy Survey [viewer](http://legacysurvey.org/viewer) around RA=253.3, Dec=29.2 degrees.

In preparation for DESI (Dark Energy Spectroscopic Experiment), the imaging Legacy Survey (including DECaLS) will be used to select over 30 million targets for the spectroscopy campain. The main DESI sample will comprise emission-line galaxies (ELGs), luminous red galaxies (LRGs), and QSOs. There will also be a Bright Galaxy Survey (BGS) targetting a magnitude-limited sample (r<19.5), including stars within the Milky Way.

**Figure 2.** Same section of the DECaLS DR3 image as in Figure 1 but with labels from a preliminary target selection algorithm applied to DR2 observations, and selecting for ELGs (emission-line galaxies), LRGs (luminous red galaxies), and QSOs.

Object Shapes/Types

The object shape (2D light profile) is modeled by the Tractor (Lang, Hogg & Mykytyn) as part of the procedure to compute model photometry.

Possible shapes for DECaLS DR3:

  • PSF (point spread function: size will vary with the seeing of the observations)
  • SIMP (“simple” galaxies: round, exponential profile with 0.45″ effective radius)
  • EXP (exponential profile; spiral galaxies)
  • DEV (deVaucouleurs profile; elliptical galaxies)
  • COMP (composite deVaucouleurs+exponential at same centroid)
**Figure 3.** Images of galaxies including a nearby elliptical galaxy, a nearby spiral galaxy, and a QSO.

Summary of shape fitting and assignment:

  1. Fit all sources with PSF and SIMP; keep if significant (5 sigma); otherwise discard.
  2. If SIMP is better model, fit again with EXP and DEV; keep if 3 sigma improvement
  3. Fit again with COMP (composite of EXP and DEV); keep if further 3 sigma improvement

Magnitudes and Colors

Magnitudes are obtained through a set of filters similar to the u,g,r,i,z set used for SDSS. In this work, we use:

  • *g,r,z* from the Dark Energy Camera (DECam)
  • *W1, W2* forced photometry in [WISE](http://www.nasa.gov/mission_pages/WISE/main/) channels 1 & 2 (3.4 and 4.6 microns)

Colors are defined as a difference between magnitudes in two bands. A "redder" color means that the object is comparatively brighter in the redder (i.e., longer wavelength) band. Conversely, a "bluer" color means that the object is comparatively brighther in the bluer (i.e., shorter wavelenght) band.

We will use the following colors:

  • *g-r*
  • *r-z*
  • *z-W1*
  • *W1-W2*

Query DECaLS Tractor Photometry Catalog

The photometry is derived from Tractor modeling of sources, and the database includes model photometry, type (shape), as well as aperture photometry in various aperture sizes. In this work, we will use Tractor model magnitudes.

In [25]:
# Write query statement (adql)
query = """
        SELECT dered_mag_g as gmag, dered_mag_r as rmag, dered_mag_z as zmag, 
               dered_mag_w1 as w1mag, dered_mag_w2 as w2mag, type,
               snr_g, snr_r, snr_z, ra, dec 
        FROM ls_dr3.tractor_primary
        WHERE (snr_g>3 and snr_r>3 and snr_z>3)
        LIMIT 100000"""

# type             = object type (PSF, SIMP, EXP, DEV, COMP)
# dered_mag_g,r,z  = AB magnitudes in DECam g,r,z bands corrected for Galactic reddening
# dered_mag_w1,w2  = AB magnitudes in WISE bands W1 & W2 corrected for Galactic reddening
#
# WHERE: requirement that S/N>3 in each DECaLS band

print query
        SELECT dered_mag_g as gmag, dered_mag_r as rmag, dered_mag_z as zmag, 
               dered_mag_w1 as w1mag, dered_mag_w2 as w2mag, type,
               snr_g, snr_r, snr_z, ra, dec 
        FROM ls_dr3.tractor_primary
        WHERE (snr_g>3 and snr_r>3 and snr_z>3)
        LIMIT 100000
In [35]:
%%time
# If the query is LARGE (e.g., LIMIT > 100,000 or NO LIMIT),
# run with Async=True

# Execute the Query synchronously if short
response = qc.query(token, adql=query, fmt='csv')

print 'Time'
Time
CPU times: user 94 ms, sys: 14 ms, total: 108 ms
Wall time: 2.59 s
In [28]:
# Reformat output into a table
result = Table.read(StringIO(response), format='csv')  #dictionary

# Print a few rows from the result table
print result[:5]
len(result)
   gmag      rmag      zmag   ...   snr_z         ra           dec      
--------- --------- --------- ... --------- ------------- --------------
21.861643 21.573879 21.473604 ... 25.698452 316.858107514 -1.10833708696
21.805428 20.300669  18.60038 ... 270.78387 316.855405345 -1.10466866907
21.766214 20.849237 20.280699 ... 42.201897 316.856820721 -1.10317664752
24.301088  23.40452 23.168093 ...    3.7672 316.860981255 -1.10567611692
22.818901 21.761492  21.12874 ... 34.124718 316.864666304 -1.12261113057
Out[28]:
100000

Optical Color-Color Diagram

In [29]:
# Select range of interest
thres = 5.   #threshold value for S/N (here, making it more stringent than query)
keep = (result['snr_g']>thres)&(result['snr_r']>thres)&(result['snr_z']>thres)

# Colors
g_r   = result['gmag'][keep] - result['rmag'][keep]
r_z   = result['rmag'][keep] - result['zmag'][keep]
z_w1  = result['zmag'][keep] - result['w1mag'][keep]
w1_w2 = result['w1mag'][keep] - result['w2mag'][keep]

print len(g_r)


# Classification per object type
objtype = result['type'][keep]
75186
In [30]:
col0 = r_z   #r-z color
col1 = g_r   #g-r color

# 2D-histogram of objects
fig, ax1 = plt.subplots(1, 1, figsize=(9, 8))
im1 = ax1.hexbin(col0, col1, bins='log', cmap=plt.cm.viridis,
               mincnt=1, extent=(-1., 3, -1., 3))
ax1.set_ylabel('g-r',fontsize=20)
ax1.set_xlabel('r-z',fontsize=20)

#color bar
cb = plt.colorbar(im1,label='log(N)')
In [31]:
col0 = r_z   #r-z color
col1 = g_r   #g-r color

# List of the types
typeList = ['All','PSF','SIMP','EXP','DEV','COMP']

fig, axes = plt.subplots(2, 3, figsize=(15, 10), sharex='all', sharey='all')
axes = axes.flatten()

for i,typ in enumerate(typeList):
    ax = axes[i]
    if typ=='All':
        selec = np.ones(len(col0),dtype='bool')   #for All objects: array of 1's
    else:
        selec = (objtype==typ)    #for each type, select on the type parameter
    im = ax.hexbin(col0[selec], col1[selec], bins='log', cmap=plt.cm.viridis,
               mincnt=1, extent=(-1., 3, -1., 3))
    if i>=3: ax.set_xlabel('r-z',fontsize=20)
    if i%3==0: ax.set_ylabel('g-r',fontsize=20)
    ax.text(0.1,0.9,typ,transform=ax.transAxes,fontsize=14,color='red',backgroundcolor='white')

plt.subplots_adjust(wspace=0., hspace=0.)

plt.show()

Optical-Infrared Colors

Combine information from optical/near-infrared DECam observations and from WISE mid-infrared observations. The latter were extracted using "forced photometry" at the position of DECaLS objects.

In [32]:
# First, try z-W1 versus infrared color W1-W2
col0 = z_w1
col1 = w1_w2

# List
typeList = ['All','PSF','SIMP','EXP','DEV','COMP']

fig, axes = plt.subplots(2, 3, figsize=(15, 10), sharex='all', sharey='all')
axes = axes.flatten()

for i,typ in enumerate(typeList):
    ax = axes[i]
    if typ=='All':
        selec = np.ones(len(objtype),dtype='bool')
    else:
        selec = (objtype==typ)
    im = ax.hexbin(col0[selec], col1[selec], bins='log', cmap=plt.cm.viridis,
               mincnt=1, extent=(-3, 3, -2., 1.5))
    if i>=3: ax.set_xlabel('z-w1',fontsize=20)
    if i%3==0: ax.set_ylabel('w1-w2',fontsize=20)
    ax.text(0.8,0.9,typ,transform=ax.transAxes,fontsize=14,color='red',backgroundcolor='white')

plt.subplots_adjust(wspace=0.02, hspace=0.02)

plt.show()