Voronoi tesselation¶

This notebook shows how to perform voronoi tesselation using the vorbinpackage. The following cell is the mandatory data loading.

In [1]:

import astropy.units as u
from xsb_fluc.data.cluster import Cluster

cluster = Cluster(
    imglink='data/A2142/mosaic_a2142.fits.gz',
    explink='data/A2142/mosaic_a2142_expo.fits.gz',
    bkglink='data/A2142/mosaic_a2142_bkg.fits.gz',
    reglink='data/A2142/src_ps.reg',
    nhlink='data/A2142/A2142_nh.fits',
    ra=239.58615,
    dec=27.229433,
    r_500=1.403*u.Mpc, 
    redshift=0.09,
)

cluster = cluster.reduce_to_r500(1).rebin(5)
cluster.plot_cluster()

import astropy.units as u from xsb_fluc.data.cluster import Cluster cluster = Cluster( imglink='data/A2142/mosaic_a2142.fits.gz', explink='data/A2142/mosaic_a2142_expo.fits.gz', bkglink='data/A2142/mosaic_a2142_bkg.fits.gz', reglink='data/A2142/src_ps.reg', nhlink='data/A2142/A2142_nh.fits', ra=239.58615, dec=27.229433, r_500=1.403*u.Mpc, redshift=0.09, ) cluster = cluster.reduce_to_r500(1).rebin(5) cluster.plot_cluster()

WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Stellar reference frame 
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: EQUINOX = '2000.0 ' / Coordinate system equinox 
a floating-point value was expected. [astropy.wcs.wcs]

Using this cluster, we can perform a Voronoi tesselation using the implementation from the vorbinpackage, which requires a signal-to-noise function. Here, I do not target a specific SNR but rather a given number of counts (700, which is far greater then what you usually do). The vorbin package works with 1D inputs, so we flatten the various data from our cluster. We also filter for pixels that are too far away from the center and that are not exposed.

In [2]:

import numpy as np
from vorbin.voronoi_2d_binning import voronoi_2d_binning

def sn_func(index, signal, noise):

    return np.sum(signal[index])

X,Y = np.indices(cluster.shape)
valid = (cluster.exp>0)&(cluster.coords.separation(cluster.center) < cluster.t_500)

x = X[valid].flatten()
y = Y[valid].flatten()
s = cluster.img[valid].flatten()
n = np.ones_like(s)

bin_number, xBin, yBin, xBar, yBar, sn, nPixels, scale = voronoi_2d_binning(x, y, s, n, 700, pixelsize=1, quiet = True, sn_func=sn_func, plot=False)

np.savetxt('data/A2142/voronoi.txt', bin_number)

import numpy as np from vorbin.voronoi_2d_binning import voronoi_2d_binning def sn_func(index, signal, noise): return np.sum(signal[index]) X,Y = np.indices(cluster.shape) valid = (cluster.exp>0)&(cluster.coords.separation(cluster.center) < cluster.t_500) x = X[valid].flatten() y = Y[valid].flatten() s = cluster.img[valid].flatten() n = np.ones_like(s) bin_number, xBin, yBin, xBar, yBar, sn, nPixels, scale = voronoi_2d_binning(x, y, s, n, 700, pixelsize=1, quiet = True, sn_func=sn_func, plot=False) np.savetxt('data/A2142/voronoi.txt', bin_number)

Plotting the results¶

Here I just plot the result. In practice, the image is flattened and reduced to a 1D vector, but this script is made to unflatten it and plot the corresponding Voronoi regions.

In [3]:

import matplotlib.pyplot as plt
import cmasher as cmr
from matplotlib.colors import LogNorm

indexes = (cluster.exp>0)&(cluster.coords.separation(cluster.center) < cluster.t_500)
y_ref, x_ref = cluster.y_ref[indexes], cluster.x_ref[indexes]
img = np.array(cluster.img[indexes].astype(int))
exp = np.array(cluster.exp[indexes].astype(np.float32))

voronoi_img = np.zeros_like(cluster.img)*np.nan
bin_number_img = np.zeros_like(cluster.img)*np.nan
voronoi_img_flat = voronoi_img[indexes]
bin_number_img_flat = np.zeros_like(voronoi_img_flat)

unique_bin = np.unique(bin_number)

for i, number in enumerate(unique_bin):
    bin_index = (bin_number == number)
    voronoi_img_flat[bin_index] = np.sum(img[bin_index])
    bin_number_img_flat[bin_index] = number
    
voronoi_img[indexes] = voronoi_img_flat
bin_number_img[indexes] = bin_number_img_flat

fig, axs = plt.subplots(
    figsize=(12, 5),
    nrows=1,
    ncols=2,
    subplot_kw={'projection': cluster.wcs}
)

img_plot = axs[0].imshow(voronoi_img, norm=LogNorm(vmin=500), cmap=cmr.cosmic)

plt.colorbar(
    mappable=img_plot,
    ax=axs[0],
    location='bottom',
    label="Image (Photons)"
)

number_plot = axs[1].imshow(bin_number_img, cmap='tab20')

plt.colorbar(
    mappable=number_plot,
    ax=axs[1],
    location='bottom',
    label="Bin number"
)

plt.show();

import matplotlib.pyplot as plt import cmasher as cmr from matplotlib.colors import LogNorm indexes = (cluster.exp>0)&(cluster.coords.separation(cluster.center) < cluster.t_500) y_ref, x_ref = cluster.y_ref[indexes], cluster.x_ref[indexes] img = np.array(cluster.img[indexes].astype(int)) exp = np.array(cluster.exp[indexes].astype(np.float32)) voronoi_img = np.zeros_like(cluster.img)*np.nan bin_number_img = np.zeros_like(cluster.img)*np.nan voronoi_img_flat = voronoi_img[indexes] bin_number_img_flat = np.zeros_like(voronoi_img_flat) unique_bin = np.unique(bin_number) for i, number in enumerate(unique_bin): bin_index = (bin_number == number) voronoi_img_flat[bin_index] = np.sum(img[bin_index]) bin_number_img_flat[bin_index] = number voronoi_img[indexes] = voronoi_img_flat bin_number_img[indexes] = bin_number_img_flat fig, axs = plt.subplots( figsize=(12, 5), nrows=1, ncols=2, subplot_kw={'projection': cluster.wcs} ) img_plot = axs[0].imshow(voronoi_img, norm=LogNorm(vmin=500), cmap=cmr.cosmic) plt.colorbar( mappable=img_plot, ax=axs[0], location='bottom', label="Image (Photons)" ) number_plot = axs[1].imshow(bin_number_img, cmap='tab20') plt.colorbar( mappable=number_plot, ax=axs[1], location='bottom', label="Bin number" ) plt.show();