import pyfits
from scipy import ndimage
from scipy import signal
import numpy as np
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
Read in a few images:
j = pyfits.getdata('/astro/ferguson1/ferguson/python_course/hlsp_hudf09_hst_wfc3ir_hudf09_F125W_v1_sci.fits')
h = pyfits.getdata('/astro/ferguson1/ferguson/python_course/hlsp_hudf09_hst_wfc3ir_hudf09_F160W_v1_sci.fits')
weight = pyfits.getdata('/astro/ferguson1/ferguson/python_course/hlsp_hudf09_hst_wfc3ir_hudf09_F160W_v1_wht.fits')
fig, ax = subplots(figsize=(8, 8))
ax.imshow(h, vmin=-0.05, vmax=0.1, cmap=cm.gray) # Display the H-band image
<matplotlib.image.AxesImage at 0x107a18fd0>
Determine a 3-sigma-clipped mean H-band background and standard deviation in pixels that have data (weight > 0):
filt = weight > 0
for n in range(10): # 10 clipping iterations
bkgd = h[filt].mean()
rms = h[filt].std()
filt = filt & (h > bkgd - 3. * rms) & (h < bkgd + 3. * rms)
print bkgd, rms
0.00207353405252 0.0443872555883 0.00123125074417 0.00675409262619 0.000567707646256 0.00226487980137 0.000309228836302 0.00129066125838 0.0002034082969 0.00100020279111 0.000158265381271 0.000900917552984 0.000139326139565 0.000863149701725 0.000131785046243 0.000848352850905 0.000128567800048 0.000842378567371 0.000127299835581 0.000839969679175
Smooth the image, and find connected regions above some threshold in the smoothed image:
Unsharp mask the image and keep positive peaks
hsmooth = ndimage.gaussian_filter(h, 2.5) # Smooth with a 2.5 pixel gaussian
hsmoother = ndimage.gaussian_filter(h, 9.5) # Smooth with a bigger gaussian
snr = (hsmooth - bkgd) / rms
diff = hsmooth - hsmoother # Unsharp mask the image
peaks = np.choose((snr > 5) & (diff > 0), (0., diff)) # Keep just the positive peaks
fig, (ax1, ax2) = subplots(1, 2, sharex=True, sharey=True, figsize=(16, 8))
ax1.imshow(h, vmin=-0.005, vmax=0.01, cmap=cm.gray)
ax2.imshow(peaks, vmin=-0.01, vmax=0.01, cmap=cm.gray)
<matplotlib.image.AxesImage at 0x14f3e7950>
Segment this image, and display this segmentation map along with the original:
labels, n = ndimage.measurements.label(peaks) # Segment the image
print n
fig, ax = subplots(figsize=(8, 8))
ax.imshow(labels)
1082
<matplotlib.image.AxesImage at 0x107697150>
Find the centroids of all the sources in J and H.
positions_h = ndimage.measurements.center_of_mass(h, labels, range(n))
positions_j = ndimage.measurements.center_of_mass(j, labels, range(n))
n = len(positions_j)
print n
1082
Calculate the position offsets between the two bands:
pj, ph = np.array(positions_j), np.array(positions_h)
xj, xh = pj[:,0], ph[:,0]
yj, yh = pj[:,1], ph[:,1]
dx = xh - xj
dy = yh - yj
dr = np.sqrt(dx ** 2 + dy ** 2)
Plot offsets as a vector field.
fig, (ax1, ax2) = subplots(1, 2, sharex=True, sharey=True, figsize=(16, 8))
ax1.imshow(h, vmin=-0.005, vmax=0.01, cmap=cm.gray)
ax2.imshow(h, vmin=-0.005, vmax=0.01, cmap=cm.gray)
ax2.quiver(yh, xh, dx, dy, color='y', units='xy', scale=0.01) # Note x,y are reversed on the image display
<matplotlib.quiver.Quiver at 0x14f475890>