KDE_mult_overdens.py

import numpy as np
from scipy.ndimage.filters import maximum_filter
from scipy.ndimage.filters import gaussian_filter
from scipy.ndimage.morphology import generate_binary_structure, binary_erosion
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.offsetbox as offsetbox
# from mpl_toolkits.axes_grid1 import make_axes_locatable


def get_coords():
    '''
    Return coordinates from file.
    '''
    # Each sub-list in file_data is a row of the file.
    # file_data = np.loadtxt("TR24_clean.dat")
    file_data = np.loadtxt("CLUSTER.DAT")

    # Extract coordinates and zip them into two lists.
    ra, dec = zip(*file_data)[1], zip(*file_data)[2]

    # ra_rang, dec_rang = max(ra) - min(ra), max(dec) - min(dec)
    # print 'RA range: ', ra_rang
    # print 'DEC range: ', dec_rang

    return ra, dec


def get_2d_histo(x_data, y_data, N_bins, std_dev):
    '''
    Return 2D histogram for the positional data.
    '''
    xmin, xmax = min(x_data), max(x_data)
    ymin, ymax = min(y_data), max(y_data)

    # Calculate the number of bins used.
    x_rang, y_rang = (xmax - xmin), (ymax - ymin)

    # Bin width to create the 2D histogram.
    bin_width = min(x_rang, y_rang) / N_bins

    # Number of bins in x,y given the bin width.
    binsxy = [int(x_rang / bin_width), int(y_rang / bin_width)]

    # hist_2d is the 2D histogram, *edges store the edges of the bins.
    hist_2d, xedges, yedges = np.histogram2d(
        x_data, y_data, range=[[xmin, xmax], [ymin, ymax]], bins=binsxy)

    hist = gaussian_filter(hist_2d, std_dev, mode='nearest')  # 'constant')

    return hist, xedges, yedges


def detect_peaks(image):
    """
    Takes an image and detect the peaks using the local maximum filter.
    Returns a boolean mask of the peaks (i.e. 1 when
    the pixel's value is the neighborhood maximum, 0 otherwise)
    """

    # define an 8-connected neighborhood
    neighborhood = generate_binary_structure(2, 2)

    # apply the local maximum filter; all pixel of maximal value
    # in their neighborhood are set to 1
    local_max = maximum_filter(image, footprint=neighborhood) == image
    # local_max is a mask that contains the peaks we are
    # looking for, but also the background.
    # In order to isolate the peaks we must remove the background from the
    # mask.

    # we create the mask of the background
    background = (image == 0)

    # a little technicality: we must erode the background in order to
    # successfully subtract it form local_max, otherwise a line will
    # appear along the background border (artifact of the local maximum filter)
    eroded_background = binary_erosion(
        background, structure=neighborhood, border_value=1)

    # we obtain the final mask, containing only peaks,
    # by removing the background from the local_max mask
    detected_peaks = local_max - eroded_background

    return detected_peaks


def get_peaks_coords(peaks, ra, dec, histo_edges):
    '''
    Transform peak bin coordinates into actual coordinates.
    '''

    peaks_coords = []
    # for each smoothed chart.
    for i, h_p in enumerate(peaks):
        xedges, yedges = histo_edges[i]
        max_dens_coords = []
        # For each x column in the histogram.
        for x_p_indx, x_col in enumerate(h_p):
            # For each y column in the histogram.
            for y_p_indx, y_p in enumerate(x_col):
                if y_p:
                    x_cent_bin, y_cent_bin = x_p_indx, y_p_indx
                    # x,y coords of the center of the above bin.
                    x_cent_pix = np.average(xedges[x_cent_bin:x_cent_bin + 2])
                    y_cent_pix = np.average(yedges[y_cent_bin:y_cent_bin + 2])
                    max_dens_coords.append([x_cent_pix, y_cent_pix])

        peaks_coords.append(max_dens_coords)

    return peaks_coords


def get_kde_dens(ra, dec, peaks_coords):
    '''
    '''

    # Obtain Gaussian KDE in coordinates.
    values = np.vstack([ra, dec])
    kernel = stats.gaussian_kde(values)

    max_dens = []
    for p in peaks_coords:
        # Evaluate kernel in peak coordinates.
        positions = np.vstack([zip(*p)[0], zip(*p)[1]])
        k_pos = kernel(positions)
        # Store density values in those positions.
        max_dens.append(k_pos)

    return max_dens


def get_filter_dens(peaks_coords, max_dens, thresh=0.75):
    '''
    '''

    max_dens_coords, max_dens_filt = [], []
    for p, md in zip(peaks_coords, max_dens):
        max_d = max(md)
        coords_f, dens_f = [], []
        for i, dens in enumerate(md):

            # Only store those that are larger than a certain threshold.
            if dens > (max_d * thresh):
                dens_f.append(dens)
                coords_f.append(p[i])

        max_dens_coords.append(coords_f)
        max_dens_filt.append(dens_f)

    return max_dens_coords, max_dens_filt


def make_plot(std_dev_lst, clust_histo, ra, dec, peaks_coords, max_dens,
              max_dens_coords, max_dens_filt, thresh):
    '''
    '''

    fig = plt.figure(figsize=(15, 50))
    gs = gridspec.GridSpec(16, 3)

    i, j = 0, 0
    for h, p, m, mc, mf in zip(clust_histo, peaks_coords, max_dens,
                               max_dens_coords, max_dens_filt):
        print 'Plotting: ', i, i + 1, i + 2

        ax1 = plt.subplot(gs[i])
        plt.tick_params(axis='both', which='major', labelsize=8)
        plt.imshow(h.transpose(), origin='lower', interpolation='none')
        # ax.imshow(h.transpose(), origin='lower')
        # Text box.
        ob = offsetbox.AnchoredText(r'$\sigma={}$'.format(std_dev_lst[j]),
                                    loc=2, prop=dict(size=8))
        ob.patch.set(alpha=0.65)
        ax1.add_artist(ob)
        ax1.invert_xaxis()
        j += 1

        # ax2 = plt.subplot(gs[i + 1])
        # cm = plt.cm.gist_earth_r
        # plt.imshow(p.transpose(), origin='lower', cmap=cm)
        # ax2.invert_xaxis()

        siz = 5 + 95 * (np.array(m) / max(m)) ** 2
        ax2 = plt.subplot(gs[i + 1])
        plt.tick_params(axis='both', which='major', labelsize=8)
        plt.xlabel(r'$RA\,(^{\circ})$', fontsize=10)
        plt.ylabel(r'$DEC\,(^{\circ})$', fontsize=10)
        x_c, y_c = zip(*p)[0], zip(*p)[1]
        plt.xlim(min(ra), max(ra))
        plt.ylim(min(dec), max(dec))
        cm = plt.cm.get_cmap('RdYlBu_r')
        plt.scatter(x_c, y_c, c=m, cmap=cm, s=siz, lw=0.25)
        ax2.invert_xaxis()
        ax2.set_aspect('equal')

        # v_min_mp, v_max_mp = min(m), max(m)
        siz = 5 + 95 * (np.array(mf) / max(mf)) ** 6
        ax3 = plt.subplot(gs[i + 2])
        plt.tick_params(axis='both', which='major', labelsize=8)
        plt.xlabel(r'$RA\,(^{\circ})$', fontsize=10)
        plt.ylabel(r'$DEC\,(^{\circ})$', fontsize=10)
        x_c, y_c = zip(*mc)[0], zip(*mc)[1]
        plt.xlim(min(ra), max(ra))
        plt.ylim(min(dec), max(dec))
        cm = plt.cm.get_cmap('RdYlBu_r')
        plt.scatter(x_c, y_c, c=mf, cmap=cm, s=siz, lw=0.25)
        # vmin=v_min_mp, vmax=v_max_mp)
        # Text box.
        ob = offsetbox.AnchoredText('thresh={}'.format(thresh), loc=2,
                                    prop=dict(size=8))
        ob.patch.set(alpha=0.5)
        ax3.add_artist(ob)
        # Invert axis and set aspect.
        ax3.invert_xaxis()
        ax3.set_aspect('equal')

        # Increase counter.
        i += 3

    # Output png file.
    fig.tight_layout()
    plt.savefig('/home/gabriel/Descargas/peak_detect.png', dpi=300)


# Get coordinates.
ra, dec = get_coords()

# Generate 2D histogram of the coordinates, using several different numbers
# of bins and standard deviating values for the Gaussian smoothing.
clust_histo, histo_edges = [], []
std_dev_lst = [2, 2.25, 2.5, 2.75, 3., 3.25, 3.5, 3.75, 4., 4.25, 4.5,
               4.75, 5., 5.25, 5.5, 5.75]
for N_bins in [100]:
    for std_dev in std_dev_lst:
        h, xedges, yedges = get_2d_histo(ra, dec, N_bins, std_dev)
        clust_histo.append(h)
        histo_edges.append([xedges, yedges])

# Detect peaks.
peaks = []
for h in clust_histo:
    peaks.append(detect_peaks(h))

# Transform peaks bin coordinates into actual coordinates.
peaks_coords = get_peaks_coords(peaks, ra, dec, histo_edges)

# Obtain KDE values for each peak in each smoothed cluster.
max_dens = get_kde_dens(ra, dec, peaks_coords)

thresh = 0.9
# Filter found max density coordinates below a certain threshold.
max_dens_coords, max_dens_filt = get_filter_dens(peaks_coords, max_dens,
                                                 thresh)

# Plot results.
make_plot(std_dev_lst, clust_histo, ra, dec, peaks_coords, max_dens,
          max_dens_coords, max_dens_filt, thresh)