Quick Start

Loading Example Data

The following data is a cutout of a group of bright galaxies in Abell 2744. The original data was acquired by the Hubble Frontier Fields team via the WFC3 instrument in the F105W filter and can be directly downloaded from the Mikulski Archive for Space Telescopes.

Loading Image

You can use astropy’s CCDData to load the example data and visualize it through matplotlib.

from astropy.nddata import CCDData

image = CCDData.read('data/abell_2744_dwarf_galaxy_f105w.fits.gz')

import numpy as np
from matplotlib import pyplot as plt

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['image.origin'] = 'lower'

vmax = 0.05 # Use the image std as max and min of all plots
vmin = - vmax

plt.imshow(image.data, vmin=vmin, vmax=vmax)
plt.title("Galaxy in Abell 2744")
plt.xlabel("Pixels")
plt.ylabel("Pixels")
plt.show()

Making Cutouts

Use astropy’s Cutout2D function to make cutouts of sources:

from astropy.nddata import Cutout2D

cutout_image = Cutout2D(image, position=(100,100), size=40)


plt.imshow(cutout_image.data, vmin=vmin, vmax=vmax)
plt.title("Cutout Galaxy")
plt.xlabel("Pixels")
plt.ylabel("Pixels")
plt.show()

Defining AstroPy Sérsic Models

You can define a 2D Sérsic model using astropy as follows. PetroFit provides a helper function get_default_sersic_bounds provides a Python dictionary with default parameter bounds, which are useful when fitting (used to constrain the parameter space).

bounds = {
    'amplitude': (0., None),
    'r_eff': (0, None),
    'n': (0, 10),
    'ellip': (0, 1),
    'theta': (-2 * np.pi, 2 * np.pi),
}

from astropy.modeling import models
from petrofit.modeling import get_default_sersic_bounds

sersic_model = models.Sersic2D(

        amplitude=0.1, # Intensity at r_eff
        r_eff=10, # Effective or half-lilght radius
        n=1, # Sersic index
        x_0=20, # center of model in the x direction
        y_0=20, # center of model in the y direction
        ellip=0.1, # Ellipticity
        theta=0.0, # Rotation angle in radians, counterclockwise from the positive x-axis.

        bounds=get_default_sersic_bounds(), # Parameter bounds
)

To add x_0 and y_0 bounds to the default bounds, you can update the dictionary as you would a regular Python dictionary:

bound_dict = get_default_sersic_bounds()

bound_dict.update( {'x_0': (10, 30),  'y_0': (10, 30)} )

bound_dict

{'amplitude': (0.0, None),
 'r_eff': (0, None),
 'n': (0, 10),
 'ellip': (0, 1),
 'theta': (-6.283185307179586, 6.283185307179586),
 'x_0': (10, 30),
 'y_0': (10, 30)}

You can update the model bounds as follows:

sersic_model.bounds.update(bound_dict)

If you need to do this once, you can pass your updates to the get_default_sersic_bounds function as follows:

bound_dict = get_default_sersic_bounds( {'x_0': (10, 30),  'y_0': (10, 30)} )
bound_dict

{'amplitude': (0.0, None),
 'r_eff': (0, None),
 'n': (0, 10),
 'ellip': (0, 1),
 'theta': (-6.283185307179586, 6.283185307179586),
 'x_0': (10, 30),
 'y_0': (10, 30)}

Making Compound Models (Combining Models)

You can combine multiple models to form a compound model by adding, subtracting, multiplying, and dividing individual models. For example we add the Sersic model from the last section to itself to form a two component sersic model (notice that the number of parameters double):

compound_sersic_model = sersic_model + sersic_model

from petrofit.modeling import print_model_params

print_model_params(compound_sersic_model)

0.1000  amplitude_0
10.0000 r_eff_0
1.0000  n_0
20.0000 x_0_0
20.0000 y_0_0
0.1000  ellip_0
0.0000  theta_0
0.1000  amplitude_1
10.0000 r_eff_1
1.0000  n_1
20.0000 x_0_1
20.0000 y_0_1
0.1000  ellip_1
0.0000  theta_1

Making a PSF Convolved Model

The petrofit PSFConvolvedModel2D is a Fittable2DModel that adds PSF convolution and model to image sampling to astropy core models. PSFConvolvedModel2D makes an image of the underlying model and samples it onto a grid. The model image is then convolved with a PSF if one is provided. Since PSFConvolvedModel2D is a Fittable2DModel, it can be used to fit model images to data images. For example, we wrap an astropy Sersic2D model in this doc with PSFConvolvedModel2D, which produces an oversampled and PSF convolved version of the Sersic profile at each iteration of the Levenberg-Marquardt fitting algorithm. Note that ``PSFModel`` is deprecated and replaced by ``PSFConvolvedModel2D``.

############
# Load PSF #
############

from astropy.io import fits

# Load PSF image (2D array)
PSF = fits.getdata('data/f105w_psf.fits.gz')

# Normalize PSF
PSF = PSF / PSF.sum()

# Note that the PSF shape is odd on all sides
print("PSF Shape = {}".format(PSF.shape))

# Plot PSF and use vmax and vmin to show difraction spikes
plt.imshow(PSF, vmin=0, vmax=PSF.std()/10)
plt.show()

PSF Shape = (51, 51)

from petrofit.modeling import PSFConvolvedModel2D

psf_sersic_model = PSFConvolvedModel2D(sersic_model, psf=PSF, oversample=4)

Converting Models to Images

To convert any model to an image use the model_to_image function.

from petrofit.modeling import model_to_image

# Size of model image
size = 40

# sersic model image
model_image = model_to_image(model=sersic_model, size=size)

# PSF convolved model image
psf_model_image = model_to_image(model=psf_sersic_model, size=size)

Plot model image

fig, axs = plt.subplots(1,2, figsize=(15,7.5))

plt.sca(axs[0])
plt.imshow(model_image, vmin=vmin, vmax=vmax)
plt.title('Sersic Model Image')

plt.sca(axs[1])
plt.imshow(psf_model_image, vmin=vmin, vmax=vmax)
plt.title('PSF Convolved Model Image')

plt.show()

Fitting Model to Image

Use the fit_model function to fit 2D models to images as follows:

%%time

from petrofit.modeling import fit_model

fitted_model, fit_info = fit_model(
    cutout_image.data, psf_sersic_model,
    maxiter=10000,
    epsilon=1.4901161193847656e-08,
    acc=1e-09,
)

CPU times: user 401 ms, sys: 88 µs, total: 401 ms
Wall time: 401 ms

Convert the fitted model into an image

from petrofit.modeling import plot_fit

axs, model_image, residual_image = plot_fit(fitted_model, cutout_image, return_images=True,
                                            vmax=vmax, vmin=vmin, figsize=[24, 12])

for ax in axs:
    ax.set_xlabel('Pixels')
    ax.set_ylabel('Pixels')

plt.show()

Fitting Multiple Sources

If the locations of the sources are known, one can fit all sources at the same time by creating a compound model. Note that x_0 and y_0 are known beforehand using photometric centroids. Below, a compound model of 3 Sérsic components is defined and the original image is fit (i.e not the cutout we have been working with).

# Center elliptical galaxy we have been fitting:
galaxy_model_1 = models.Sersic2D(

        amplitude=0.1, # Intensity at r_eff
        r_eff=10, # Effective or half-lilght radius
        n=1.7384901, # Sersic index
        x_0=99.97722657736085, # center of model in the x direction
        y_0=99.12324178530918, # center of model in the y direction
        ellip=0.1, # Ellipticity
        theta=0.0, # Rotation angle in radians, counterclockwise from the positive x-axis.

        bounds=get_default_sersic_bounds(), # Parameter bounds
)

# Football shaped galaxy
galaxy_model_2 = models.Sersic2D(

        amplitude=0.1, # Intensity at r_eff
        r_eff=10, # Effective or half-lilght radius
        n=1, # Sersic index
        x_0=138.56315299695075, # center of model in the x direction
        y_0=89.27757468116197, # center of model in the y direction
        ellip=0.7, # Ellipticity
        theta=0.7, # Rotation angle in radians, counterclockwise from the positive x-axis.

        bounds=get_default_sersic_bounds(), # Parameter bounds
)

# Large galaxy near the bottom corner
galaxy_model_3 = models.Sersic2D(

        amplitude=0.1, # Intensity at r_eff
        r_eff=10, # Effective or half-lilght radius
        n=1, # Sersic index
        x_0=178.72302596615611, # center of model in the x direction
        y_0=63.506754312433046      , # center of model in the y direction
        ellip=0.2, # Ellipticity
        theta=0.0, # Rotation angle in radians, counterclockwise from the positive x-axis.

        bounds=get_default_sersic_bounds(), # Parameter bounds
)

Make compound PSF model as follows:

all_galaxies_model = galaxy_model_1 + galaxy_model_2 + galaxy_model_3

all_galaxies_psf_model = PSFConvolvedModel2D(all_galaxies_model, psf=PSF)

Fit the model

%%time

from petrofit.modeling import fit_model

all_galaxies_fitted_model, fit_info = fit_model(
    image.data, all_galaxies_psf_model,
    maxiter=10000,
    epsilon=1.4901161193847656e-08,
    acc=1e-09,
)

/home/docs/checkouts/readthedocs.org/user_builds/petrofit/conda/v0.4.0/lib/python3.10/site-packages/astropy/modeling/functional_models.py:3060: RuntimeWarning: divide by zero encountered in true_divide
  z = np.sqrt((x_maj / a) ** 2 + (x_min / b) ** 2)

CPU times: user 7.37 s, sys: 7.89 ms, total: 7.38 s
Wall time: 7.38 s

plot_fit(all_galaxies_fitted_model, image, return_images=False,
         vmax=vmax, vmin=vmin, figsize=[24, 12])
plt.show()

Looks like the bottom corner galaxy is a spiral, let us add another component for the spiral and fit again:

%%time

# Redefine model with an extra component for galaxy 3
all_galaxies_model = galaxy_model_1 + galaxy_model_2 + galaxy_model_3 + galaxy_model_3

# PSF model
all_galaxies_psf_model = PSFConvolvedModel2D(all_galaxies_model, psf=PSF)

# Fit the model
all_galaxies_fitted_model, fit_info = fit_model(
    image.data, all_galaxies_psf_model,
    maxiter=10000,
    epsilon=1.4901161193847656e-08,
    acc=1e-09,
)

# Plot the fit
plot_fit(all_galaxies_fitted_model, image, return_images=False,
         vmax=vmax, vmin=vmin, figsize=[24, 12])
plt.show()

/home/docs/checkouts/readthedocs.org/user_builds/petrofit/conda/v0.4.0/lib/python3.10/site-packages/astropy/modeling/functional_models.py:3062: RuntimeWarning: divide by zero encountered in true_divide
  return amplitude * np.exp(-bn * (z ** (1 / n) - 1))

CPU times: user 34.8 s, sys: 131 ms, total: 34.9 s
Wall time: 34.8 s

Fitting Image Backgrounds

The fit_background function can be used to fit the background pixels using a 2D plane. It will sigma clip the pixels (sigma value provided by the user) and fit a 2D plane to the clipped image. Users can also provide their own 2D models.

from petrofit.modeling import fit_background

bg_model, fit_info = fit_background(image, sigma=3.0)
bg_image = model_to_image(bg_model, size=(image.shape[1], image.shape[0]))

WARNING: Model is linear in parameters; consider using linear fitting methods. [astropy.modeling.fitting]

Plot backround

fig, axs = plt.subplots(1,2)

plt.sca(axs[0])
plt.imshow(bg_image)
plt.title("Background Image")
plt.xlabel("Pixels")
plt.ylabel("Pixels")

plt.sca(axs[1])
plt.imshow(image.data - bg_image, vmin=vmin, vmax=vmax)
plt.title("Background subtracted image")
plt.xlabel("Pixels")
plt.ylabel("Pixels")
plt.show()