In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage as ndi
from astropy.visualization import simple_norm
from astropy.modeling import models
import photutils
import time
import statmorph
%matplotlib inline
In [2]:
ny, nx = 240, 240
y, x = np.mgrid[0:ny, 0:nx]
sersic_model = models.Sersic2D(amplitude=1, r_eff=20, n=1.5, x_0=0.5*nx, y_0=0.4*ny,
ellip=0.5, theta=0.5)
image = sersic_model(x, y)
plt.imshow(image, cmap='gray', origin='lower',
norm=simple_norm(image, stretch='log', log_a=10000))
Out[2]:
In [3]:
size = 20 # on each side from the center
sigma_psf = 2.0
y, x = np.mgrid[-size:size+1, -size:size+1]
psf = np.exp(-(x**2 + y**2)/(2.0*sigma_psf**2))
psf /= np.sum(psf)
plt.imshow(psf, origin='lower', cmap='gray')
Out[3]:
Now we convolve the image with the PSF.
In [4]:
image = ndi.convolve(image, psf)
plt.imshow(image, cmap='gray', origin='lower',
norm=simple_norm(image, stretch='log', log_a=10000))
Out[4]:
In [5]:
np.random.seed(1)
snp = 100.0
image += (1.0 / snp) * np.random.standard_normal(size=(ny, nx))
plt.imshow(image, cmap='gray', origin='lower',
norm=simple_norm(image, stretch='log', log_a=10000))
Out[5]:
The code will ask for one of two input arguments: (1) a weight map, which is a 2D array (of the same size as the input image) representing one standard deviation at each pixel value, or (2) the gain, a scalar that can be multiplied by the science image to obtain the number of electron counts per pixel. This is used internally by statmorph to calculate the weight map.
Here we assume, rather arbitrarily, that there is an average of 1000 electron counts/pixel at the effective radius (where we defined the amplitude as 1.0), so that the gain is 1000.0.
In [6]:
gain = 1000.0
In [7]:
threshold = photutils.detect_threshold(image, 1.5)
npixels = 5 # minimum number of connected pixels
segm = photutils.detect_sources(image, threshold, npixels)
Although statmorph is designed to process all the sources labeled by the segmentation map, in this example we only focus on the main (largest) source found in the image.
In [8]:
# Keep only the largest segment
label = np.argmax(segm.areas) + 1
segmap = segm.data == label
plt.imshow(segmap, origin='lower', cmap='gray')
Out[8]:
We regularize a bit the shape of the segmentation map:
In [9]:
segmap_float = ndi.uniform_filter(np.float64(segmap), size=10)
segmap = segmap_float > 0.5
plt.imshow(segmap, origin='lower', cmap='gray')
Out[9]:
Now that we have all the required data, we are ready to measure the morphology of the source just created. Note that we include the PSF as a keyword argument. In principle, this results in more correct Sersic profile fits, although it can also make the code run slower, depending on the size of the PSF.
In [10]:
start = time.time()
source_morphs = statmorph.source_morphology(image, segmap, gain=gain, psf=psf)
print('Time: %g s.' % (time.time() - start))
In general, source_morphs is a list of objects, each corresponding to a labeled source in the image. Here we only focus on the first labeled source.
In [11]:
morph = source_morphs[0]
Now we print some of the morphological properties just calculated:
In [12]:
print('xc_centroid =', morph.xc_centroid)
print('yc_centroid =', morph.yc_centroid)
print('ellipticity_centroid =', morph.ellipticity_centroid)
print('elongation_centroid =', morph.elongation_centroid)
print('orientation_centroid =', morph.orientation_centroid)
print('xc_asymmetry =', morph.xc_asymmetry)
print('yc_asymmetry =', morph.yc_asymmetry)
print('ellipticity_asymmetry =', morph.ellipticity_asymmetry)
print('elongation_asymmetry =', morph.elongation_asymmetry)
print('orientation_asymmetry =', morph.orientation_asymmetry)
print('rpetro_circ =', morph.rpetro_circ)
print('rpetro_ellip =', morph.rpetro_ellip)
print('rhalf_circ =', morph.rhalf_circ)
print('rhalf_ellip =', morph.rhalf_ellip)
print('r20 =', morph.r20)
print('r80 =', morph.r80)
print('Gini =', morph.gini)
print('M20 =', morph.m20)
print('F(G, M20) =', morph.gini_m20_bulge)
print('S(G, M20) =', morph.gini_m20_merger)
print('sn_per_pixel =', morph.sn_per_pixel)
print('C =', morph.concentration)
print('A =', morph.asymmetry)
print('S =', morph.smoothness)
print('sersic_amplitude =', morph.sersic_amplitude)
print('sersic_rhalf =', morph.sersic_rhalf)
print('sersic_n =', morph.sersic_n)
print('sersic_xc =', morph.sersic_xc)
print('sersic_yc =', morph.sersic_yc)
print('sersic_ellip =', morph.sersic_ellip)
print('sersic_theta =', morph.sersic_theta)
print('sky_mean =', morph.sky_mean)
print('sky_median =', morph.sky_median)
print('sky_sigma =', morph.sky_sigma)
print('flag =', morph.flag)
print('flag_sersic =', morph.flag_sersic)
Note that the fitted Sersic profile is in pretty good agreement with the "true" Sersic profile that we originally defined (n=1.5, r_eff=20, etc.). However, such agreement tends to deteriorate somewhat at higher noise levels and larger Sersic indices (not to mention that real galaxies are not always well described by Sersic profiles). Other morphological measurements that are more robust to noise, which are also calculated by statmorph, include the Gini-M20 (Lotz et al. 2004), CAS (Conselice 2003) and MID (Freeman et al. 2013) statistics, as well as the outer asymmetry (Wen et al. 2014) and shape asymmetry (Pawlik et al. 2016).
Also note that statmorph calculates two different "bad measurement" flags (where 0 means good measurement and 1 means bad):
flag : indicates a problem with the basic morphological measurements.
flag_sersic : indicates if there was a problem/warning during the Sersic profile fitting.
In general, flag==0 should always be enforced, while flag_sersic==0 should only be used when interested in Sersic fits (which might fail for merging galaxies and other "irregular" objects).
In [13]:
ny, nx = image.shape
y, x = np.mgrid[0:ny, 0:nx]
fitted_model = statmorph.ConvolvedSersic2D(
amplitude=morph.sersic_amplitude,
r_eff=morph.sersic_rhalf,
n=morph.sersic_n,
x_0=morph.sersic_xc,
y_0=morph.sersic_yc,
ellip=morph.sersic_ellip,
theta=morph.sersic_theta)
fitted_model.set_psf(psf) # always required when using ConvolvedSersic2D
image_model = fitted_model(x, y)
bg_noise = (1.0 / snp) * np.random.standard_normal(size=(ny, nx))
fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(131)
ax.imshow(image, cmap='gray', origin='lower',
norm=simple_norm(image, stretch='log', log_a=10000))
ax.set_title('Original image')
ax = fig.add_subplot(132)
ax.imshow(image_model + bg_noise, cmap='gray', origin='lower',
norm=simple_norm(image, stretch='log', log_a=10000))
ax.set_title('Fitted model')
ax = fig.add_subplot(133)
residual = image - image_model
ax.imshow(residual, cmap='gray', origin='lower',
norm=simple_norm(residual, stretch='linear'))
ax.set_title('Residual')
Out[13]:
In [14]:
from statmorph.utils.image_diagnostics import make_figure
fig = make_figure(morph)
In [15]:
fig.savefig('tutorial.png', dpi=150)