Finding Galaxies within 100 kpc of Q1532+0613

  |   Source

This will demonstrate how scikit's image processing library scikit-image is useful to my data for Q1532+0613.

The data was taken by the Gemini North during the 2015A observing run (GN-2015A-Q77). The data can be accessed here

In [76]:
# Import matplotlib (plotting), skimage (image processing) and interact (user interfaces)
# This enables their use in the Notebook.
%matplotlib inline
import astropy.io.fits as fits

#from skimage import data
from skimage.feature import blob_doh
from skimage.color import rgb2gray

from ipywidgets import interact, fixed

# Extract the first 500px square of the Hubble Deep Field.
#path = 'blog/posts/Gemini_Data/Imaging'
#image = path+
#image = data.hubble_deep_field()[0:500, 0:500]
image, header = fits.getdata('/home/ec2-user/blog/posts/Gemini_Data/Imaging/obj_stacked_r.fits', 0, header=True)
image = image/image.max()
image_gray = rgb2gray(image)
def plot_blobs(max_sigma=30, threshold=0.1, gray=False):
    """
    Plot the image and the blobs that have been found.
    """
    blobs = blob_doh(image_gray, max_sigma=max_sigma, threshold=threshold)
    
    fig, ax = plt.subplots(figsize=(8,8))
    ax.set_title('Galaxies in the Hubble Deep Field')
    
    if gray:
        ax.imshow(image_gray, interpolation='nearest', cmap='gray_r')
        circle_color = 'red'
    else:
        ax.imshow(image, interpolation='nearest')
        circle_color = 'yellow'
    for blob in blobs:
        y, x, r = blob
        c = plt.Circle((x, y), r, color=circle_color, linewidth=2, fill=False)
        ax.add_patch(c)
In [74]:
# Use interact to explore the galaxy detection algorithm.
interact(plot_blobs, max_sigma=(10, 40, 2), threshold=(0.005, 0.02, 0.001));
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
<ipython-input-73-29b83ec9d469> in plot_blobs(max_sigma, threshold, gray)
     21     Plot the image and the blobs that have been found.
     22     """
---> 23     blobs = blob_doh(image_gray, max_sigma=max_sigma, threshold=threshold)
     24 
     25     fig, ax = plt.subplots(figsize=(8,8))

/home/ec2-user/.pyenv/versions/3.5.1/lib/python3.5/site-packages/skimage/feature/blob.py in blob_doh(image, min_sigma, max_sigma, num_sigma, threshold, overlap, log_scale)
    406 
    407     hessian_images = [_hessian_matrix_det(image, s) for s in sigma_list]
--> 408     image_cube = np.dstack(hessian_images)
    409 
    410     local_maxima = peak_local_max(image_cube, threshold_abs=threshold,

/home/ec2-user/.pyenv/versions/3.5.1/lib/python3.5/site-packages/numpy/lib/shape_base.py in dstack(tup)
    366 
    367     """
--> 368     return _nx.concatenate([atleast_3d(_m) for _m in tup], 2)
    369 
    370 def _replace_zero_by_x_arrays(sub_arys):

MemoryError: 
In [71]:
print(image_gem.dtype, image.size)
>f4 7160832
In [79]:
t=image_gem/image_gem.max()
In [81]:
t.ndim
Out[81]:
2
In [ ]:
 
In [ ]: