In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.gridspec as gs
import deblend
import astropy.table
import fitsio

Use the skeleton.py code in the willitdeblend repository, (but executed using the code at https://github.com/DarkEnergyScienceCollaboration/WeakLensingDeblending), to create the input image/catalog contained in the following fits files.


In [2]:
fits=fitsio.FITS("pairs.fits")
catalog=astropy.table.Table.read("pair_table.fits")

Some routines to deblend and make plots of blends/deblends.


In [3]:
def img_and_colorbar(fig, grid):
    inner_grid = gs.GridSpecFromSubplotSpec(100, 100, subplot_spec=grid, wspace=0, hspace=0)
    img_ax = fig.add_subplot(inner_grid[:,0:90])
    cbar_ax = fig.add_subplot(inner_grid[:,91:95])
    return img_ax, cbar_ax
    
def plot_blend(fits, catalog, idx):
    hdr = fits[3*idx].read_header()
    blend = fits[3*idx].read()
    sh = blend.shape
    xlim = [0, blend.shape[1]-1]
    ylim = [0, blend.shape[0]-1]
    unblends = [fits[3*idx+1].read(), fits[3*idx+2].read()]
    
    dx = catalog['dx'][2*idx:2*idx+2]*5+2048-hdr['GS_XMIN']-0.5
    dy = catalog['dy'][2*idx:2*idx+2]*5+2048-hdr['GS_YMIN']-0.5
    peaks = [(dx[0]-0.5*sh[1], dy[0]-0.5*sh[0]), (dx[1]-0.5*sh[1], dy[1]-0.5*sh[0])]
    _, _, deblends = deblend.deblend(blend, peaks)    

    center_color='cyan'
    
    #Plotting
    fig = plt.figure(figsize=(12,10))
    grid = gs.GridSpec(3, 17, wspace=1, hspace=0.3)

    #Blend    
    img_ax, cbar_ax = img_and_colorbar(fig, grid[0,0:5])
    img_ax.set_xlim(xlim); img_ax.set_ylim(ylim)
    img_ax.set_title("Blend")
    im = img_ax.imshow(blend)
    vlim = (im.norm.vmin, im.norm.vmax)
    img_ax.plot(dx, dy, marker='+', ls='', color=center_color, ms=10, mew=2)
    plt.colorbar(im, cax=cbar_ax)
    
    #Unblends
    img_ax, cbar_ax = img_and_colorbar(fig, grid[0,6:11])
    img_ax.set_xlim(xlim); img_ax.set_ylim(ylim)
    img_ax.set_title("Unblend[0]")
    im = img_ax.imshow(unblends[0], vmin=vlim[0], vmax=vlim[1])
    img_ax.plot(dx[0], dy[0], marker='+', ls='', color=center_color, ms=10, mew=2)
    plt.colorbar(im, cax=cbar_ax)
    
    img_ax, cbar_ax = img_and_colorbar(fig, grid[0,12:17])
    img_ax.set_xlim(xlim); img_ax.set_ylim(ylim)
    img_ax.set_title("Unblend[1]")
    im = img_ax.imshow(unblends[1], vmin=vlim[0], vmax=vlim[1])
    img_ax.plot(dx[1], dy[1], marker='+', ls='', color=center_color, ms=10, mew=2)
    plt.colorbar(im, cax=cbar_ax)
    
    #Deblends
    img_ax, cbar_ax = img_and_colorbar(fig, grid[1,6:11])
    img_ax.set_xlim(xlim); img_ax.set_ylim(ylim)
    img_ax.set_title("Deblend[0]")
    im = img_ax.imshow(deblends[0], vmin=vlim[0], vmax=vlim[1])
    img_ax.plot(dx[0], dy[0], marker='+', ls='', color=center_color, ms=10, mew=2)
    plt.colorbar(im, cax=cbar_ax)
    
    img_ax, cbar_ax = img_and_colorbar(fig, grid[1,12:17])
    img_ax.set_xlim(xlim); img_ax.set_ylim(ylim)
    img_ax.set_title("Deblend[1]")
    im = img_ax.imshow(deblends[1], vmin=vlim[0], vmax=vlim[1])
    img_ax.plot(dx[1], dy[1], marker='+', ls='', color=center_color, ms=10, mew=2)
    plt.colorbar(im, cax=cbar_ax)    
    
    #Residuals
    img_ax, cbar_ax = img_and_colorbar(fig, grid[2,6:11])
    img_ax.set_xlim(xlim); img_ax.set_ylim(ylim)
    img_ax.set_title("Resid[0]")
    im = img_ax.imshow(unblends[0] - deblends[0])
    img_ax.plot(dx[0], dy[0], marker='+', ls='', color=center_color, ms=10, mew=2)
    plt.colorbar(im, cax=cbar_ax)
    
    img_ax, cbar_ax = img_and_colorbar(fig, grid[2,12:17])
    img_ax.set_xlim(xlim); img_ax.set_ylim(ylim)
    img_ax.set_title("Resid[1]")
    im = img_ax.imshow(unblends[1] - deblends[1])
    img_ax.plot(dx[1], dy[1], marker='+', ls='', color=center_color, ms=10, mew=2)
    plt.colorbar(im, cax=cbar_ax)

And plot!


In [4]:
for i in range(5):
    plot_blend(fits, catalog, i)


deblend.py:20: RuntimeWarning: invalid value encountered in divide
  template_frac = template/template_sum

In [ ]: