In [1]:
__author__ = 'David Nidever <>, Knut Olsen <>, Robert Nikutta <>, Stephanie Juneau <>' # single string; emails in <>
__version__ = '20180105' # yyyymmdd; version datestamp of this notebook
__datasets__ = ['nsc_dr1']  # enter used datasets by hand

Dwarf galaxies in NSC DR1

In this notebook, we demonstrate the discovery of classical and faint Milky Way dwarf companions in NSC DR1. We query the database around the positions of known dwarfs and apply filtering techniques to reveal the dwarfs as spatial overdensities of filtered sources.

David Nidever, Knut Olsen, Robert Nikutta, Stephanie Juneau & NOAO Data Lab Team

Disclaimer & attribution

If you use this notebook for your published science, please acknowledge the following:

Imports and setup

In [3]:
# Python 2/3 compatibility
from __future__ import print_function # to use print() as a function in Python 2

    input = raw_input # use 'input' function in both Python 2 and 3
except NameError:

# std lib
from getpass import getpass

# 3rd party
import pandas as pd
import numpy as np
import pylab as plt
import matplotlib
from astropy import utils, io, convolution, stats, wcs
from astropy.visualization import make_lupton_rgb
from astropy import units as u
from astropy.stats import median_absolute_deviation as mad
%matplotlib inline

# Data Lab
from dl import authClient as ac, queryClient as qc, storeClient as sc, helpers

#Simple Image Access (SIA) service
from pyvo.dal import sia
svc = sia.SIAService(DEF_ACCESS_URL)

# Quiet the Astropy warnings
import warnings
from astropy.utils.exceptions import AstropyWarning
warnings.simplefilter('ignore', category=AstropyWarning)


In [4]:
# Either get token for anonymous user
token = ac.login('anonymous')

# Authenticated users please uncomment the next line
#token = ac.login(input("Enter user name: "),getpass("Enter password: "))


We're going to select the database and define some functions that we'll use later.

In [5]:
#The NSC DR1 schema
except Exception as e:
Schema: nsc_dr1

      Table Name   Description
      ----------   -----------
            chip   CCD information table
        coverage   Survey coverage table
        exposure   Exposures contributing to the catalog
     nsc_allwise   NSC vs ALLWISE 1-arcsec crossmatch
         nsc_des   NSC vs DES 1-arcsec crossmatch
          object   Primary object table

In [6]:
# A function to retrieve data from a point on the sky
def getData (ra,dec,radius=1.0,columns='*'):

    query_template = \
    """SELECT {0:s} FROM nsc_dr1.object
       WHERE q3c_radial_query(ra,dec,{1:f},{2:f},{3:f})"""

    query = query_template.format(columns,ra,dec,radius)
        result = qc.query(token,sql=query) # by default the result is a CSV formatted string
    except Exception as e:
    df = helpers.convert(result,'pandas')
    return df
In [8]:
# A Mexican-hat convolution filter
def dwarf_filter (ra,dec,fwhm_small=2.0,fwhm_big=20):

    # Based on Koposov et al. (2008).
    # Code by Ken Mighell and Mike Fitzpatrick.
    # Minor edits by Rorbert Nikutta.
    x, y = ra, dec

    # Information about declination (y) [degrees]
    ymean = (y.min() + y.max()) / 2.0
    ydiff_arcmin = (y.max() - y.min()) * 60.0 # convert from degrees to arcmin

    # Information about right ascension (x) [degrees in time]:
    xdiff = x.max() - x.min() # angular separation [degrees (time)] 
    xmean = (x.min() + x.max())/2.0

    # convert from degrees in time to separation in angular degrees:
    xdiff_angular = (x.max() - x.min()) * np.cos(ymean*(np.pi/180.0))

    # convert from degress to arcmin
    xdiff_angular_arcmin = xdiff_angular * 60.0 

    # Get the number of one-arcmin pixels in the X and Y directions:
    nx = np.rint (xdiff_angular_arcmin).astype('int')
    ny = np.rint (ydiff_arcmin).astype('int')
    # Create a two-dimensional histogram of the raw counts:
    Counts, xedges, yedges  = np.histogram2d (x, y, (nx,ny) )
    extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
    raw_hist = np.rot90(Counts).copy() # hack around Pythonic weirdness

    # Make the small and big Gaussian kernels with a standard deviation
    # of the given FWHM in arcmin^2 pixels.
    kernel_small = convolution.Gaussian2DKernel(fwhm_small/2.35,factor=1)
    kernel_big = convolution.Gaussian2DKernel(fwhm_big/2.35,factor=1)

    # Compute the differential convolution kernels.
    conv_big = convolution.convolve (raw_hist, kernel_big)
    conv_small = convolution.convolve (raw_hist, kernel_small)
    conv_delta = conv_small - conv_big
    delta = conv_delta.copy()

    # Compute statistics and the floor
    mean = np.mean (delta, dtype='float64')
    sigma = np.std (delta, dtype='float64')
    sigmaRaw = np.std(raw_hist,dtype='float64')
    median = np.median (delta)                       # not used
    floor = mean 

    # Clip to specified limits.
    clipped = delta.copy()
    clipped[ delta < floor ] = floor

    # Return the computed fields.
    return raw_hist, extent, delta, clipped, sigma
In [9]:
# A little function to download the deepest stacked images
#   adapted from R. Nikutta
def download_deepest_image(ra,dec,fov=0.1,band='g'):
    imgTable =,dec), (fov/np.cos(dec*np.pi/180), fov), verbosity=2).votable.to_table()
    print("The full image list contains", len(imgTable), "entries")
    sel0 = imgTable['obs_bandpass'].astype(str)==band
    sel = sel0 & ((imgTable['proctype'].astype(str)=='Stacked') & (imgTable['prodtype'].astype(str)=='image')) # basic selection
    Table = imgTable[sel] # select
    if (len(Table)>0):
        row = Table[np.argmax(Table['exptime']'float'))] # pick image with longest exposure time
        url = row['access_url'].decode() # get the download URL
        print ('downloading deepest image...')
        image = io.fits.getdata(,cache=True,show_progress=False,timeout=120))

        print ('No image available.')
    return image
In [10]:
# Multi panel image plotter
def plot_images(images,geo=None,panelsize=4,bands=list('gri'),
    n = len(images)
    if geo is None: geo = (n,1)
    fig = plt.figure(figsize=(geo[0]*panelsize,geo[1]*panelsize))
    for j,img in enumerate(images):
        ax = fig.add_subplot(geo[1],geo[0],j+1)
        if img is not None:
            vmin = np.median(img)-2*np.std(img)
            vmax = np.median(img)+2*np.std(img)
            ax.imshow(img,origin='lower',interpolation='none',cmap=cmap,norm=matplotlib.colors.LogNorm(vmin=vmin, vmax=vmax))
            ax.set_title('%s band' % bands[j])

Chapter 1 - Photometry containing dwarfs

We take the positions of the first eight dwarf galaxies and query the database around them.

In [11]:
# Locations of Dwarfs from McConnachie et al. (2012), Bechtol et al. (2015) and others
ra = [86.4,260.1,185.43,56.09]
dec = [34.7,-22.2,-31.99,-43.53]
name = ['Draco','Carina','Hya II','Eri II']
radius = 0.85 # degrees
columns = '''ra,dec,gmag,gerr,imag,rmag,zmag,class_star,fwhm'''
df_dict = {dwarf: pd.DataFrame() for dwarf in name}
In [12]:
from astropy.coordinates import name_resolve
In [13]:
def resolve_coordinates(name):
        coords = name_resolve.get_icrs_coordinates(name)
    except Exception as e:

    ra ='deg').value
    dec ='deg').value      

    return coords, ra, dec        
In [14]:
for j,name_ in enumerate(name):
    coords, ra_, dec_ = resolve_coordinates(name_)
    df = getData(ra_,dec_,radius=radius,columns=columns)
    df_dict[name_] = df
    print(str(len(df))+' objects found')
<SkyCoord (ICRS): (ra, dec) in deg
    ( 260.051625,  57.915361)> 260.051625 57.915361
SELECT ra,dec,gmag,gerr,imag,rmag,zmag,class_star,fwhm FROM nsc_dr1.object
       WHERE q3c_radial_query(ra,dec,260.051625,57.915361,0.850000)
Returning Pandas dataframe
74870 objects found
<SkyCoord (ICRS): (ra, dec) in deg
    ( 100.402888, -50.966196)> 100.402888 -50.966196
SELECT ra,dec,gmag,gerr,imag,rmag,zmag,class_star,fwhm FROM nsc_dr1.object
       WHERE q3c_radial_query(ra,dec,100.402888,-50.966196,0.850000)
Returning Pandas dataframe
429762 objects found
Hya II
<SkyCoord (ICRS): (ra, dec) in deg
    ( 185.42542, -31.98528)> 185.42542 -31.98528
SELECT ra,dec,gmag,gerr,imag,rmag,zmag,class_star,fwhm FROM nsc_dr1.object
       WHERE q3c_radial_query(ra,dec,185.425420,-31.985280,0.850000)
Returning Pandas dataframe
304486 objects found
Eri II
<SkyCoord (ICRS): (ra, dec) in deg
    ( 56.0878, -43.5332)> 56.0878 -43.5332
SELECT ra,dec,gmag,gerr,imag,rmag,zmag,class_star,fwhm FROM nsc_dr1.object
       WHERE q3c_radial_query(ra,dec,56.087800,-43.533200,0.850000)
Returning Pandas dataframe
141685 objects found
CPU times: user 2.4 s, sys: 555 ms, total: 2.96 s
Wall time: 38.7 s

Chapter 2 - Filtering and plotting the dwarfs

We filter the photometry to include objects with S/N>10, that are point-like (CLASS_STAR, and FWHM), and that are relatively blue (g-i < 1.0). We convolve the result with the spatial filter defined above.

In [15]:
# Create a WCS for a tangent plane projection in our region
def get_wcs(ra,dec,image,fov=1.,unit='deg',projection=("RA---TAN","DEC--TAN")):
    npix = image.shape[0]
    crpix = npix/2 + 1
    cdelt = fov/float(npix)
    w = wcs.WCS(naxis=2)
    w.wcs.cunit = (unit,unit)
    w.wcs.crpix = (crpix,crpix)
    w.wcs.cdelt = np.array((-cdelt,cdelt))
    w.wcs.ctype = projection
    w.wcs.crval = (ra,dec),
    return w
In [16]:
# Loop over the dwarfs, filter then and plot a density map
nrow, ncol = 2, 2
fig = plt.figure(figsize=(6*ncol,6*nrow))
for j,name_ in enumerate(df_dict):
    df = df_dict[name_]
    keep = (df['gmag']<90) & (df['rmag']<90) & \
           1.087/df['gerr']>10 & \
           ((df['gmag']-df['rmag'])<1.0) & (df['gmag']>18) & ((df['gmag']-df['rmag'])>(-0.5)) & \
           (np.abs(df['class_star'])>0.8) & (df['fwhm']<2) 
    df2 = df[keep]
    raw_hist, extent, delta, clipped, sigma = dwarf_filter(df2['ra'],df2['dec'])
    w = get_wcs(ra[j],dec[j],clipped,fov=1.)
    ax = fig.add_subplot(nrow,ncol,j+1,projection=w)
    im = plt.imshow(clipped)
Eri II
Hya II

Chapter 3 - Retrieve images

We demonstrate how to retrieve images though the SIA (Simple Image Access) service.

In [17]:
# Get g, r and i stacked image for the third galaxy
rac = ra[2]
decc = dec[2]
band = 'g'
gimage = download_deepest_image(rac, decc, fov=0.1, band=band) # FOV in deg
band = 'r'
rimage = download_deepest_image(rac, decc, fov=0.1, band=band) # FOV in deg
band = 'i'
iimage = download_deepest_image(rac, decc, fov=0.1, band=band) # FOV in deg
The full image list contains 2717 entries
downloading deepest image...
The full image list contains 2717 entries
downloading deepest image...
The full image list contains 2717 entries
downloading deepest image...
In [18]:
# Now plot the image per band and a 3-color image
gimage2 = gimage-np.median(gimage)
rimage2 = rimage-np.median(rimage)
iimage2 = iimage-np.median(iimage)
img = make_lupton_rgb(iimage2, rimage2, gimage2, Q=10, stretch=30)

fig = plt.figure(figsize=[12,12])
ax = fig.add_subplot(1,1,1)

2304.05 39016.5
3063.22 39893.4
2365.8 39021.3