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)
In [ ]: