completeness_2d_hist.py

import numpy as np
import astropy.wcs as wcs
import astropy.io.fits as fits
import matplotlib.pyplot as plt
from matplotlib import cm
import scipy.stats as sps
from astropy.table import Table, Column, vstack
from astropy.coordinates import SkyCoord, Distance, Angle
from astropy import units as u
from matplotlib.ticker import NullFormatter


def compute_completeness(data_file, orig_cat, realobj_sep=2.5, mag_tol=1, make_plots=True):
    '''
    input: data_file: results of fake object insertion: RA, Dec, X(arcmin), Y(arxmin), x(pix),y(pix),mag_in,mag_out
           orig_cat: RA/Dec of objects in real catalog
           realobj_sep: minimum distance from a real object for a fake obj to be considered detected (arcsec)
           mag_tol: absolute value of magnitude difference to be considered matched.

    output: 
    '''
    art_results = Table.read(data_file, format='ascii.commented_header')
    orig_cat = Table.read(orig_cat, format='ascii.no_header')

    rgc = np.sqrt(art_results['Major']**2+art_results['Minor']**2)
    art_results.add_column(Column(rgc, name = 'Rgc'))

    # find objects in artificial results list that are close to real objects in image
    matched_art_orig = filter_cat(art_results['RA'],art_results['Dec'], orig_cat['col1'],orig_cat['col2'], sep=realobj_sep)
    n_pos_filt = len(matched_art_orig)
    print '{} objects filtered with separation {:.2f}'.format(n_pos_filt, realobj_sep)

    # find objects in artificial results list whose magnitudes are too different from input
    nmag_filt = np.logical_and(art_results['mag_out']>0, np.abs(art_results['mag_in'] - art_results['mag_out'])>=mag_tol).sum()
    print '{} objects filtered with delta mag {:.2f}'.format(nmag_filt, mag_tol)

    # now make a detection mask
    art_results['mag_out'][matched_art_orig] = 0 # mark the too-close ones as undetected
    detected = np.logical_and(art_results['mag_out']>0, np.abs(art_results['mag_in'] - art_results['mag_out'])<mag_tol)
    art_results.add_column(Column(detected, name = 'detected'))

    if make_plots:
        xedge,yedge,top, bot = compute_histos(art_results, 'mag_in','Rgc', 'detected')
        fancy_2d_plot(xedge,yedge,top,bot)

    return(art_results)


def fancy_2d_plot(xedge,yedge, top_hist, bot_hist):
    '''plots 2-D and associated 1-D completeness histograms
       (ratios of top_hist/bot_hist per bin)
    xedge: bin edges in x-direction
    yedge: bin edges in y-direction
    top_hist: numerator histogram
    bot_hist: denominator histogram
    '''

    nullfmt   = NullFormatter()         # no labels
    
    # definitions for the axes
    left, width = 0.1, 0.65
    bottom, height = 0.1, 0.65
    bottom_h = left_h = left+width+0.02
    
    rect_scatter = [left, bottom, width, height]
    rect_histx = [left, bottom_h, width, 0.2]
    rect_histy = [left_h, bottom, 0.15, height]
    
    # start with a rectangular Figure
    f = plt.figure(figsize=(8,8))
    
    ax2d = plt.axes(rect_scatter)
    axHistx = plt.axes(rect_histx)
    axHisty = plt.axes(rect_histy)
    
    # no labels
    axHistx.xaxis.set_major_formatter(nullfmt)
    axHisty.yaxis.set_major_formatter(nullfmt)
    
    # the 2D plot:
    df = top_hist/bot_hist
    # note the all-important transpose method!
    plot2d = ax2d.imshow(df.T, interpolation='none', origin='low',\
      extent=[xedge[0], xedge[-1], yedge[0], yedge[-1]],aspect='auto',cmap='gray')
    ax2d.set_xlabel('Input magnitude')
    ax2d.set_ylabel('Galactocentric distance [arcmin]')
    plt.colorbar(plot2d)
    
    # the histograms
    xhist = top_hist.sum(axis=1)/bot_hist.sum(axis=1)
    axHistx.bar(left=xedge[:-1],height=xhist, width = xedge[1:]-xedge[:-1])
    axHistx.set_xlim( ax2d.get_xlim() )
    axHistx.set_ylabel('Detection fraction')
    axHistx.set_yticks([0.1, 0.4, 0.7, 1.0])

    yhist = top_hist.sum(axis=0)/bot_hist.sum(axis=0)
    axHisty.barh(bottom=yedge[:-1],width=yhist, height = yedge[1:]-yedge[:-1])
    axHisty.set_ylim( ax2d.get_ylim() )
    axHisty.set_xlim(0,0.8)
    axHisty.set_xlabel('Detection fraction')
    axHisty.set_xticks([0.2, 0.5, 0.8])
    plt.show()
    return

def compute_histos(art_results, col1, col2, mask):
    ''' computes two 2-D histograms, using the bins from the first one for the second
        art_results: astropy table with relevant data
        col1: name of column for 'X' quantity in the 2-D histogram
        col2: name of column for 'Y' quantity in the 2-D histogram
        mask: name of mask column for making the second histogram
        '''
    hist_2d_all, xedge, yedge = np.histogram2d(art_results[col1], art_results[col2], bins=8)
    hist_2d_det, xj, yj = np.histogram2d(art_results[col1][art_results[mask]], art_results[col2][art_results[mask]],\
       bins=[xedge, yedge])

    return xedge, yedge, hist_2d_det, hist_2d_all

def filter_cat(ra1, dec1, ra2, dec2, sep=2.5):
    '''find objects in cat1 that match positions in cat2
       input: RA/Dec pairs for catalog1 and catalog2

       output: indices into cat1 for objects that are near objects in cat2
    '''

    cat1_coords = SkyCoord(ra=ra1*u.deg, dec=dec1*u.deg)
    cat2_coords = SkyCoord(ra=ra2*u.deg, dec=dec2*u.deg)

    # find the objects in cat1 which are near objects in cat2
    # match: see http://astropy.readthedocs.org/en/latest/coordinates/matchsep.html#searching-around-coordinates
    idx1,idx2,sep,dist = cat2_coords.search_around_sky(cat1_coords, sep*u.arcsec)
    
    return(idx1)