Trying out methods to create profiles, i.e. assigning each bin in 2D images to a bin in a 1D profile and then computing aggregate values for that 1D bin.
An example would be:
distance
; Call np.digitize
with a d
binning larger than the pixel size. This creates a label image label
assigning each image pixel to a profile bin.counts
, background
, exposure
as well as exess excess = counts - background
and flux = excess / exposure
.import numpy as np
import pandas as pd
from astropy.io import fits
import scipy.ndimage
import tevpy.image
%pylab inline
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 5))
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
<matplotlib.figure.Figure at 0x106a69f90>
<matplotlib.figure.Figure at 0x106a69f90>
counts = fits.getdata('counts.fits')
background = fits.getdata('background.fits')
exposure = fits.getdata('exposure.fits')
excess = counts - background
flux = excess / exposure
# This is just an example of some exclusion and distance image.
# Typically this will be a pre-computed input.
significance = fits.getdata('significance_0.1.fits')
significance = np.nan_to_num(significance)
significance = scipy.ndimage.gaussian_filter(significance, 3)
exclusion = significance < 5
bin_size = 0.02
distance = bin_size * tevpy.utils.image.exclusion_distance(exclusion)
y, x = np.indices(significance.shape, dtype='int32')
plt.imshow(significance, origin='lower', vmax=20)
plt.colorbar();
plt.imshow(exclusion, origin='lower', cmap='gray')
plt.colorbar();
plt.imshow(distance, origin='lower')
plt.colorbar();
# Compute labels (in this case in distance)
# TODO: implement alternative methods to bin (e.g. equal numbers of pixels per bin)
n_bins = 5
bins = np.arange(n_bins)
distance_bins = np.linspace(distance.min(), distance.max(), n_bins)
bin_size = 0.02 # deg
distance_bins_deg = bin_size * distance_bins[:-1]
labels = np.digitize(distance.flat, distance_bins, right=True).reshape(distance.shape)
plt.imshow(labels, origin='lower')
plt.colorbar();
bin_entries = np.histogram(labels.flat, bins=bins)[0]
bin_sum = np.histogram(labels.flat, weights=significance.flat, bins=bins)[0]
#plt.hist(labels.flat, weights=significance.flat)
bin_average = bin_sum / bin_entries
# TODO: check binning! Only two pixels in the first bin???
#print distance_bins
#print bin_entries
#print distance.min()
print labels.min(), labels.max()
print bin_entries
0 19 [ 2 4240 10122 13893 17589 21212 23927 22644 22466 22178 22044 20726 18719 16433 11638 6838 3236 1604 1120]
# Check plot
plt.errorbar(distance_bins_deg, bin_entries, fmt='o');
plt.xlabel('Distance (deg)')
plt.ylabel('Number of pixels in bin');
plt.errorbar(distance_bins_deg, bin_average, xerr=0.5 * bin_size * np.diff(distance_bins), fmt='o');
plt.xlabel('Distance (deg)')
plt.ylabel('Average value')
plt.grid()
plt.ylim(-10, 20)
(-10, 20)
# We split the binning from the flux point computation because it's
# more flexible like this.
distance_bin_edges = compute_binning(distance, n_bins=5, method='same width')
flux_profile = FluxProfile(x_edges=distance_bin_edges, x=distance,
counts=counts, background=background, exposure=exposure)
flux_profile.compute()
figure(); flux_profile.plot('n_entries')
figure(); flux_profile.plot('counts_mean')
#flux_profile.save('profile.fits')
print flux_profile.data['label'].data
# It's easy to get back images from the 1D flux profile DataFrames
labels = flux_profile.data['label'].reshape(flux_profile.shape)
plt.imshow(labels, origin='lower')
plt.colorbar();