Off-centered weighted numpy histogram2d?

Question

I'm attempting to generate a model PSF from a set of observed stars. I'm following the great example provided by ali_m in this answer (MCVE below)

The 5 stars I'm using look like this:

where the center (peak intensity) is at bins [9, 9]. The results of their combination via numpy's hitsogram2d is this:

showing a peak density at bins [8, 8]. To center it at [9, 9], I have to obtain the centroids (see below) as:

cx, cy = np.array([1.] * len(stars)), np.array([1.] * len(stars))

instead. Why is this?

---

import numpy as np
from matplotlib import pyplot as plt

stars = # Uploaded here: http://pastebin.com/tjLqM9gQ

fig, ax = plt.subplots(2, 3, figsize=(5, 5))
for i in range(5):
    ax.flat[i].imshow(
        stars[i], cmap=plt.cm.viridis, interpolation='nearest',
        origin='lower', vmin=0.)
    ax.flat[i].axhline(9., ls='--', lw=2, c='w')
    ax.flat[i].axvline(9., ls='--', lw=2, c='w')
fig.tight_layout()

# (nstars, ny, nx) pixel coordinates relative to each centroid
# pixel coordinates (integer)
x, y = np.mgrid[:20, :20]
# centroids (float)
cx, cy = np.array([0.] * len(stars)), np.array([0.] * len(stars))
dx = cx[:, None, None] + x[None, ...]
dy = cy[:, None, None] + y[None, ...]

# 2D weighted histogram
bins = np.linspace(0., 20., 20)
h, xe, ye = np.histogram2d(dx.ravel(), dy.ravel(), bins=bins,
                           weights=stars.ravel())
fig, ax = plt.subplots(1, 1, subplot_kw={'aspect': 'equal'})
ax.hold(True)
ax.imshow(h, cmap=plt.cm.viridis, interpolation='nearest',
          origin='lower', vmin=0.)
ax.axhline(8., ls='--', lw=2, c='w')
ax.axvline(8., ls='--', lw=2, c='w')

plt.show()

Answer 1

The reason, the histogram is not centered at the point (9,9) where the single star intensity distribution is centered, is that the code to generate it shifts around the bins of the histogram.

As I already suggested in the comments, keep things simple. E.g. we do not need plots to see the problem. Also, I do not understand what those dx dy are, so let's avoid them.

We can then calculate the histogram by

import numpy as np

stars = # Uploaded here: http://pastebin.com/tjLqM9gQ

# The argmax of a single star results in (9,9)
single_star_argmax  = np.unravel_index(np.argmax(stars[0]), stars[0].shape)

# Create a meshgrid of coordinates (0,1,...,19) times (0,1,...,19)
y,x = np.mgrid[:len(stars[0,:,0]), :len(stars[0,0,:])]
# duplicating the grids
xcoord, ycoord = np.array([x]*len(stars)),  np.array([y]*len(stars)) 
# compute histogram with coordinates as x,y
# and [20,20] bins    
h, xe, ye = np.histogram2d(xcoord.ravel(), ycoord.ravel(), 
                           bins=[len(stars[0,0,:]), len(stars[0,:,0])],
                           weights=stars.ravel())

# The argmax of the combined stars results in (9,9)
combined_star_argmax =  np.unravel_index(np.argmax(h), h.shape)

print single_star_argmax
print combined_star_argmax
print single_star_argmax == combined_star_argmax
# prints:
# (9, 9)
# (9, 9)
# True

The only problem in the original code really was the line bins = np.linspace(0., 20., 20) which creates 20 points between 0 and 20,
0. 1.05263158 2.10526316 ... 18.94736842 20.
This scales the bin size to ~1.05 and lets your argmax occur already "earlier" then expected.
What you really want are 20 points between 0 and 19, np.linspace(0,19,20) or np.arange(0,20)
To avoid such mistakes, one can simply give the length of the original array as argument, bins=20.