Version 0.1
Today our goal is to compute the 1D radial profile of a simulated galaxy. While many software packages exist to do this, we will build this same functionality from scratch. Following this, we will be able to fit a radial profile to a galaxy image and measure its size. In practice, this problem is more complicated, but this example is, nevertheless, instructive.
By M Alpaslan (NYU) & AA Miller (CIERA/Northwestern & Adler)
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
%matplotlib notebook
The simulated galaxy for today's problem comes from the Illustris Simulation, which is a large cosmological simulation of galaxy formation.
Download the fits image for our simulated galaxy: https://northwestern.box.com/s/a8qp2afv6ed9r36v0y0vyc3vf82s1bw0
(Note - The Illustris team has made simulated images of many different galaxies available.)
Problem 1a
Using your favorite tool, display an image of the galaxy using a logarithmic stretch. Do you notice anything unusual about the image?
In [32]:
gal_im = fits.getdata("galaxy.fits")
min_cnts = np.log10(np.percentile(gal_im, 5))
max_cnts = np.log10(np.percentile(gal_im, 95))
plt.imshow(np.log10(gal_im),
vmin = min_cnts, vmax = max_cnts,
cmap='viridis', origin="lower")
plt.colorbar()
plt.tight_layout()
Problem 1b
Plot a histogram of the $\log$ of the pixel values in the simulated image.
Hint - you will likely want more than the matplotlib default number of bins.
In [33]:
plt.hist(np.log10(gal_im.flatten()), bins = 100)
plt.yscale("log")
Problem 1c
Informed by your histogram, replace any highly discrepant pixels within the image with the median pixel value for the image.
Hint - plot the results to confirm this operation occurred as expected.
Note - this procedure is strongly discouraged with real data: if weird things are happening in an image there is a reason, and it is best to understand the origins of the peculiarities. For our current purposes, this step will simplify our analysis.
In [149]:
unusual_pix = np.logical_or(gal_im < 1e-12, gal_im > 10**-2.8)
gal_median_filt = gal_im.copy()
gal_median_filt[unusual_pix] = np.median(gal_im)
plt.imshow(np.log10(gal_median_filt),
cmap='viridis', origin="lower")
plt.colorbar()
plt.tight_layout()
Before we can measure the radial profile of this galaxy, we need to compute the flux intensity at a given radius. By convention, this is done by computing the mean flux along an ellipse centered on the galaxy (where the "radius" in this case is the semi-major axis).
An ellipse with arbitrary centroid and rotation angle can be described as
$$ \frac{((x-h)\cos{\theta} - (y-k)\sin{\theta})^2}{a^2} + \frac{((x-h)\sin{\theta} + (y-k)\cos{\theta})^2}{b^2} = 1,$$where $h, k, a, b, \theta$ are the x centroid, y centroid, semi-major axis, semi-minor axes, and the rotation angle of the elipse, respectively.
This can be re-written as (it's worth checking this on the board or scratch paper):
$$\left(\frac{\sin^2{\theta}}{a^2} + \frac{\cos^2{\theta}}{b^2}\right)y^2 + \left(-2\frac{((x-h)\cos{\theta} + k\sin{\theta})\sin{\theta}}{a^2} + 2\frac{((x-h)\sin{\theta} - k\cos{\theta})\cos{\theta}}{b^2}\right)y + \left(\frac{((x-h)\cos{\theta} + k\sin{\theta})^2}{a^2} + \frac{((x-h)\sin{\theta} - k\cos{\theta})^2}{b^2} - 1\right) = 0$$Problem 2a
Write a function to solve for $y$ given $x$ and the ellipse parameters.
Hint - the function np.roots may be helpful.
In [230]:
def y_ellipse(x, h, k, a, b, theta):
p2 = np.sin(theta)**2/a**2 + np.cos(theta)**2/b**2
dummy1 = ((x-h)*np.cos(theta) + k*np.sin(theta))
dummy2 = ((x-h)*np.sin(theta) - k*np.cos(theta))
p1 = -2*dummy1*np.sin(theta)/a**2 + 2*dummy2*np.cos(theta)/b**2
p0 = dummy1**2/a**2 + dummy2**2/b**2 - 1
y_roots = np.roots([p2, p1, p0])
return y_roots
Problem 2b
Plot an elipse using your equation with $h = 0$, $k = 0$, $a = 4$, $b = 2$, and $\theta = \pi/6$.
Hint - most of this code has been provided for you.
In [289]:
x_grid = np.linspace(-3.605, 3.605, 1000)
y_vals = np.vstack((np.empty_like(x_grid), np.empty_like(x_grid)))
for x_num, x in enumerate(x_grid):
roots = y_ellipse(x, 0, 0, 4, 2, np.pi/6)
y_vals[:,x_num] = np.sort(roots)
plt.plot(np.append(x_grid, x_grid[::-1]), np.append(y_vals[0], y_vals[1][::-1]))
plt.xlim(-4, 4)
plt.ylim(-4, 4)
Out[289]:
Problem 2c
Write a function to measure the average flux of every pixel in the image that is intersected by an ellipse with parameters $k, h, a, b, \theta$.
Hint - np.floor() will be useful for selecting the correct pixel.
Hint 2 - brute force is not elegant, but it is sufficient as long as you throw away imaginary solutions.
In [257]:
def mean_radial_flux(im_data, h, k, a, b, theta):
x_grid = np.arange(im_data.shape[-1])
ellipse_flux = []
for x_num, x in enumerate(x_grid):
roots = y_ellipse(x, h, k, a, b, theta)
if np.isreal(roots).all():
ellipse_flux.append(im_data[int(np.floor(roots[0])), x])
ellipse_flux.append(im_data[int(np.floor(roots[1])), x])
return np.mean(ellipse_flux)
For this and the remaining problems, assume that the galaxy is centered at (x, y) = (128, 138), that the ellipticity (1 - b/a) is fixed at 0.3, and that $\theta = 0.8377$. [In principle we could ask you to determine each of these values, but that is outside the scope of the current problem.]
Problem 2d
Measure the mean "radial" flux of the galaxy at $a = 5$ and $a = 20$. Do these answers make sense?
In [264]:
a5_flux = mean_radial_flux(gal_median_filt, 128, 138, 5, 0.7*5, 0.8377)
a20_flux = mean_radial_flux(gal_median_filt, 128, 138, 20, 0.7*20, 0.8377)
print("The mean flux at 5 pixels is {:.6f}".format(a5_flux))
print("The mean flux at 20 pixels is {:.6f}".format(a20_flux))
Problem 2e
Plot the "mean radial profile" of the galaxy.
Hint - when choosing a grid of radii to measure, recall that the function does not make sub-pixel estimates of the flux.
In [290]:
r_grid = np.arange(2,50, dtype=float)
mean_rad_prof = np.empty_like(r_grid)
for r_num, r in enumerate(r_grid):
mean_rad_prof[r_num] = mean_radial_flux(gal_median_filt, 128, 138, r, 0.7*r, 0.8377)
plt.plot(r_grid, mean_rad_prof)
Out[290]:
The Sersic profile for galaxies can be written as:
$$\ln I(R) = \ln I_0 - kR^{1/n},$$where $I$ is the intensity (or flux), $I_0$ is the intensity at $R = 0$, $k$ is a normalization constant, and $n$ is the "Sersic index", which describes the curvature of the light profile.
We will now use a non-linear optimization routine from scipy.optimize to determine the value of the Sersic index.
Problem 3a
Create a function squared_error() that calculates the sum of the squared difference between $\ln I(R)$ and $\ln I_0 - kR^{1/n}$.
Hint - the 3 unknown parameters ($I_0$, $k$, and $n$) should be passed to the function as a single tuple called params.
In [276]:
def squared_error(params, r, I):
I_0 = params[0]
k = params[1]
n = params[2]
return sum((np.log(I) - (np.log(I_0) - k*r**(1/n)))**2)
Problem 3b
Using the L-BFGS-B method of scipy.optimize.minimize(), determine the value of the Sersic index.
Hint - you will want to provide a reasonable first guess along with bounds for the minimization function.
In [286]:
from scipy.optimize import minimize
res = minimize(squared_error, [4e-4,1,5],
method="L-BFGS-B",
bounds = ((1e-4, 1e-2), (1e-2,1e2), (1,20)),
args=(r_grid, mean_rad_prof))
print("The Sersic index is {:.3f}".format(res.x[2]))
The image that we have been working with reflects the full sensitivity of the Illustris simulation. Ground based images are observed through the atmosphere. Convolve the image with a PSF (assume a circular gaussian with FWHM = 12 pixels) and re-measure the Sersic index. How does your answer compare?
In [ ]:
# complete