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()