Reputation: 1
I am using photutils to perform isophotal fitting on the image of a galaxy.
#First a dummy iteration to get the right centre
geometry = EllipseGeometry(x0=301, y0=307, sma=20, eps=0.8, pa=np.deg2rad(110))
ellipse = Ellipse(masked_img, geometry)
isolist0 = ellipse.fit_image(sma0 = x00, minsma=4, maxsma=210, step=8, maxit=100, minit=10, linear=True, fix_center=False)
#second iteration with center fixed
geometry = EllipseGeometry(x0=np.median(isolist0.x0), y0=np.median(isolist0.y0), sma=20, eps=0.8, pa=np.deg2rad(110))
ellipse = Ellipse(masked_img, geometry)
isolist = ellipse.fit_image(sma0 = x00, minsma=4, maxsma=210, step=8, maxit=200, minit=20, linear=True, nclip=3, integrmode = 'median', fix_center=True)
So far, so good. The result is an isophote list that contains information on the aperture parameters (e.g., radius, center, position angle, ellipticity), as well as information on the flux density profile (specifically the mean intensity and its error, see https://photutils.readthedocs.io/en/stable/api/photutils.isophote.IsophoteList.html#photutils.isophote.IsophoteList).
Now, I have a second image (the same galaxy in a different band), say masked_img_2. I want to obtain the flux profile (isolist.intens and isolist.int_err) for masked_img_2, but without doing the isophotal fitting again; instead, I want to use the same apertures as in my existing isolist object. In other words, I want to do aperture photometry matching the apertures obtained before with the isophotal fitting.
I have tried a few different things but nothing quite works. What came the closest was this:
n_intensities = []
mean_intensities_err = []
# Loop over each isophote in the isolist
for isophote in isolist:
# Extract parameters
sma = isophote.sma # Semi-major axis
eps = isophote.eps # Ellipticity
pa = isophote.pa # Position angle in radians
x0 = isophote.x0 # X center
y0 = isophote.y0 # Y center
# Calculate semi-minor axis
sma_minor = sma * (1 - eps)
# Create the EllipticalAperture for the current isophote
ellipse = EllipticalAperture((x0, y0), sma, sma_minor)
# Calculate points along the ellipse
theta = np.linspace(0, 2 * np.pi, 100) # 100 points around the ellipse
x_coords = x0 + sma * np.cos(theta) * np.cos(pa) - sma_minor * np.sin(theta) * np.sin(pa)
y_coords = y0 + sma * np.cos(theta) * np.sin(pa) + sma_minor * np.sin(theta) * np.cos(pa)
# Clip coordinates to be within image boundaries
x_coords = np.clip(x_coords, 0, masked_img_2.shape[1] - 1).astype(int)
y_coords = np.clip(y_coords, 0, masked_img_2.shape[0] - 1).astype(int)
# Extract pixel values along the elliptical path
ellipse_data = masked_img_2[y_coords, x_coords]
# Calculate mean intensity along the elliptical path
if ellipse_data.size > 0:
mean_intensity = np.nanmean(ellipse_data)
int_err = np.nanstd(ellipse_data) / np.sqrt(ellipse_data.size) # RMS error
else:
mean_intensity = 0 # Or some value to indicate no data
int_err = 0
mean_intensities.append(mean_intensity)
mean_intensities_err.append(int_err)
To test it, I input masked_img instead of masked_img_2, just to see that I can recover the same as when using the isophotal fitting. The results are similar for the intensities but they are not perfect, and can lead to difference of ~0.2 in the surface brightness profiles.
Plot showing that the recovery is not perfect
So, my method does not fully work. Moreover, I would think that there is a much simpler way of using direct photutils to perform matched aperture photometry. Any help is greatly appreciated!
Upvotes: 0
Views: 190