303.2. Galaxy shapes

303.2. Galaxy shapes in DP2¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
Container Size: Large
LSST Science Pipelines version: r30.0.11
Last verified to run: 2026-09-04
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20

Learning objective: Explore the available measurements of galaxy shapes produced by the LSST pipelines and their applications.

LSST data products: Object, deep_coadd

Packages: lsst.afw, lsst.rsp, lsst.geom, lsst.gauss2d, astropy, photutils, galsim

Credit: Originally developed by the Rubin Community Science team. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.

Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

The LSST Science Pipelines make a variety of automated shape and morphology measurements for extended sources that are useful for galaxy evolution science. This notebook will teach the user about these measurements. They are performed on the deep_coadd images and appear in the Object table as part of the LSST pipelines data products. The focus will be on galaxies. Data products related to shapes for the purpose of cosmological analyses will be demonstrated elsewhere.

1.1. Import packages¶

Import numpy, a fundamental package for scientific computing with arrays in Python (numpy.org), and matplotlib, a comprehensive library for data visualization (matplotlib.org; matplotlib gallery). The io package provides input and output functions for handling data transfer.

From the pyvo package, import some functions that will enable using the image cutout tool. From astropy and photutils import packages to enable plotting images and drawing shapes on images with WCS information.

From the lsst package, import modules for accessing the Table Access Protocol (TAP) service, and image display functions from the LSST Science Pipelines (pipelines.lsst.io). Also import some geometric functions to help plot photometric apertures.

Finally, import galsim which is the galaxy simulation package, which was used to build the DP0.2 simulated images. This package and the lsst package gauss2d are useful for reconstructing the sersic profiles that are used to model the galaxies in Rubin images.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import io

from pyvo.dal.adhoc import SodaQuery

from astropy.coordinates import SkyCoord
import astropy.units as u

from photutils.aperture import SkyEllipticalAperture

from lsst.rsp import RSPDiscovery
from lsst.rsp.utils import get_pyvo_auth
import lsst.afw.display as afwDisplay
import lsst.afw.geom.ellipses as ellipses
from lsst.gauss2d import Ellipse, EllipseMajor, Covariance
from lsst.images.serialization import read_archive

import galsim as gs

1.2. Define parameters and functions¶

Define a function to generate an image cutout, using the Rubin image cutout service. Further information about the cutout tool can be found in DP2 tutorial notebook 103.6 that demonstrates the Rubin image cutout service.

Since more than one deep_coadd can overlap any galaxy, and sometimes objects can be too close to the deep_coadd edge to generate a usable cutout, this function also performs a check for which overlapping deep_coadd is the best (the object is closest to its center; indicated by best_index).

In [2]:

def make_image_cutout(ra, dec, cutout_size=0.01, band='i'):
    """
    Wrapper function to generate a cutout using the cutout tool.
    Default is to show cutout in i-band.

    Parameters
    ----------
    ra, dec : 'float'
        the ra and dec of the cutout center
    cutout_size : 'float', optional
        radial edge length in degrees of the cutout

    Returns
    -------
    cutout : 'lsst.images object'
    """

    sia_client = discovery.get_sia_client()

    eff_wl = 622.1e-09
    circle = (ra, dec, cutout_size)

    results = sia_client.search(pos=circle, calib_level=3,
                                dpsubtype='lsst.deep_coadd',
                                band=eff_wl)

    if len(results) == 0:
        raise ValueError(f"No images found for RA={ra}, Dec={dec}")

    table = results.to_table()
    distances = ((table['s_ra'] - ra)**2 + (table['s_dec'] - dec)**2)**0.5
    best_index = int(distances.argmin())
    dl_result = discovery.get_datalink_results(results[best_index])

    f"Datalink status: {dl_result.status}."

    sq = SodaQuery.from_resource(dl_result,
                                 dl_result.get_adhocservice_by_id("cutout-sync-exposure"),
                                 session=get_pyvo_auth())

    cutout_ra = ra * u.deg
    cutout_dec = dec * u.deg
    Radius = cutout_size * u.deg

    sq.circle = (cutout_ra, cutout_dec, Radius)
    cutout_bytes = sq.execute_stream().read()
    sq.raise_if_error()

    cutout = read_archive(io.BytesIO(cutout_bytes))
    return cutout

Define parameters to use colorblind-friendly colors with matplotlib.

In [3]:

plt.style.use('seaborn-v0_8-colorblind')
prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']

Set the afwDisplay backend to matplotlib to enable use of various packages to overplot on images.

In [4]:

afwDisplay.setDefaultBackend('matplotlib')

Instantiate RSPDiscovery with the DP2 release, create an instance of the TAP service, and assert that it exists.

In [5]:

discovery = RSPDiscovery("dp2")
service = discovery.get_tap_client()

2. Find galaxies in ECDFS¶

One of the deep drilling fields in the DP2 is the Extended Chandra Deep Field South (ECDFS). This section uses a galaxy that was previously identified in ECDFS that is known to be large and edge-on.

In [6]:

target_ra = 53.19488194
target_dec = -27.70342546

Here, query the DP2 Object table for the morphological parameters for this galaxy. Check that the identified galaxy is extended (i_extendedness = 1), that its i_kronFlux_flag_*, sersic_no_data_flag, and shape_flag are all = 0, indicating there is not a problem with modeling the shape and light profile.

Pick a very small search radius of 0.72 arcseconds to ensure the return is the previously specified large, edge-on galaxy.

The query will retrieve a number of shape measurements and flags that will be explained in Section 3.

In [7]:

query = "SELECT obj.objectId, obj.coord_ra, obj.coord_dec, " + \
        "obj.i_extendedness, obj.refBand, " + \
        "obj.i_kronFlux, obj.i_kronRad, obj.i_kronFlux_flag, " + \
        "obj.i_kronFlux_flag_small_radius, obj.i_kronFlux_flag_bad_radius, " + \
        "obj.shape_xx, obj.shape_xy, obj.shape_yy, obj.shape_flag, " + \
        "obj.i_ixx, obj.i_ixy, obj.i_iyy, " + \
        "obj.i_sersicFlux, obj.i_exponentialFlux,  obj.sersic_no_data_flag, " + \
        "obj.sersic_index, obj.sersic_reff_x, " + \
        "obj.sersic_reff_y, obj.sersic_rho, " + \
        "obj.sersic_reff_major, obj.sersic_reff_minor, obj.sersic_theta, " + \
        "obj.i_cModel_dev_reff_major, obj.i_cModel_dev_reff_minor, obj.i_cModel_dev_theta, " + \
        "obj.i_cModel_exp_reff_major, obj.i_cModel_exp_reff_minor, obj.i_cModel_exp_theta, " + \
        "obj.i_cModelFlux, obj.i_cModel_fracDev, " + \
        "obj.exponential_reff_major, obj.exponential_reff_minor, obj.exponential_theta " + \
        "FROM dp2.Object AS obj " + \
        "WHERE (obj.i_extendedness = 1) AND (obj.shape_flag = 0) AND " + \
        "(obj.i_kronFlux_flag_small_radius = 0) AND " + \
        "(obj.i_kronFlux_flag_bad_radius = 0) AND " + \
        "(obj.i_kronFlux_flag = 0) AND (obj.sersic_no_data_flag = 0) AND " + \
        "CONTAINS(POINT('ICRS', obj.coord_ra, obj.coord_dec), " + \
        "CIRCLE('ICRS',"+str(target_ra)+","+str(target_dec)+", 0.0002)) = 1 "

In [8]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)

Job phase is COMPLETED

In [9]:

results = job.fetch_result()
tab = results.to_table()
tab

Out[9]:

Table length=1

objectId	coord_ra	coord_dec	i_extendedness	refBand	i_kronFlux	i_kronRad	i_kronFlux_flag	i_kronFlux_flag_small_radius	i_kronFlux_flag_bad_radius	shape_xx	shape_xy	shape_yy	shape_flag	i_ixx	i_ixy	i_iyy	i_sersicFlux	i_exponentialFlux	sersic_no_data_flag	sersic_index	sersic_reff_x	sersic_reff_y	sersic_rho	sersic_reff_major	sersic_reff_minor	sersic_theta	i_cModel_dev_reff_major	i_cModel_dev_reff_minor	i_cModel_dev_theta	i_cModel_exp_reff_major	i_cModel_exp_reff_minor	i_cModel_exp_theta	i_cModelFlux	i_cModel_fracDev	exponential_reff_major	exponential_reff_minor	exponential_theta
	deg	deg			nJy	pix				pix2	pix2	pix2		pix2	pix2	pix2	nJy	nJy			pix	pix		arcsec	arcsec		arcsec	arcsec		arcsec	arcsec		nJy		arcsec	arcsec
int64	float64	float64	float32	str1	float32	float32	bool	bool	bool	float32	float32	float32	bool	float32	float32	float32	float32	float32	bool	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32
755370947912955099	53.194873939824234	-27.703425524958117	1.0	i	1.52285e+06	12.5642	False	False	False	49.6591	30.993	50.6145	False	49.6591	30.993	50.6145	1.66207e+06	1.39646e+06	False	2.20258	9.38052	9.53709	0.785649	2.52803	0.875735	-45.5537	6.28992	2.26001	-45.6413	2.37127	0.931405	-45.4914	2.29825e+06	1.0	2.07659	0.75108	-45.5496

3. Visualize galaxy morphology¶

This section will make an image cutout of the large galaxy, and compare the shapes of the different photometric and morphological measurements. This section demonstrates how to reconstruct the Gaussian Ellipse, Kron ellipse, and Sersic half-light shape, using the corresponding shape parameters in the Object table.

Below, define the SkyCoord for the galaxy's location.

In [10]:

coord = SkyCoord(ra=tab['coord_ra'][0]*u.degree,
                 dec=tab['coord_dec'][0]*u.degree, frame='icrs')

The sizes in the LSST Object table are in units of pixel. Note that the pixel scale of the LSST data is 0.2 arcseconds per pixel. Define the conversion in the cell below.

In [11]:

arcsec_per_pix = 0.2

3.1. Gaussian ellipse¶

In the next cell, use Object table shape parameters to reconstruct the galaxy shape, approximated as a 2D Gaussian. The LSST pipelines measures this shape, parameterized by three parameters or "moments" measured on the deep_coadd in each band individually: <band>_ixx, <band>_iyy, <band>_ixy (weak lensing experts may recognize these come from the measured re-Gaussianization method of Hirata & Seljak 2003, implemented by Mandelbaum et al. 2005, and called HSM moments. The corresponding moments measured on the reference band or refBand are also stored as shape_xx, shape_yy, and shape_xy. These moments have not been corrected for the point spread function (PSF; thus not applicable for weak lensing) but PSF effects are small for galaxy sizes much larger than the PSF. These shape parameters can be converted to more commonly used set of morphological parameters using the LSST package ellipses.

The cell below will demonstrate how to extract these more commonly used parameters (e.g. semi-major radius A or Rmaj, semi-minor radius B or Rmin, which can be used to obtain the axis ratio ba defined as Rmin/Rmaj, and the position angle theta. theta is defined in radians counterclockwise from the x-axis.

The semi-major and minor diameters are the "sigmas" of the 2D gaussian model approximation to the light profile. These sizes are unrelated to the Kron photometry or shape, but note that the Kron measurements use the axis ratio and position angle from this measurement.

The 2D Gaussian radius (or equally the Sersic radius) are good metrics of size, and better than kronRad, which can be come uncertain or fail for very small galaxies. Size metrics including Sersic size or gaussian sizes should be used for small galaxies. The kronRad and Kron aperture perform better for large galaxies (i.e. much larger than the PSF).

The example below converts these parameters using the ellipses package in lsst.afw.geom.

In [12]:

axes = ellipses.Axes(ellipses.Quadrupole(tab['shape_xx'][0],
                                         tab['shape_yy'][0],
                                         tab['shape_xy'][0]))
Rmaj = axes.getA()
Rmin = axes.getB()
theta = axes.getTheta()

Equivalently, this could be done using a different set of packages (the cell below demonstrates instead how to use the ellipse package in lsst.gauss2d). The cell below also validates that the same answer is retrieved with both methods.

In [13]:

ellipse = EllipseMajor(Covariance(sigma_x_sq=tab["shape_xx"][0],
                                  sigma_y_sq=tab["shape_yy"][0],
                                  cov_xy=tab["shape_xy"][0]))

ellipse_Rmaj = ellipse.r_major
ellipse_Rmin = ellipse.axrat * ellipse.r_major
ellipse_theta = ellipse.angle

print(ellipse_Rmaj, Rmaj, ' arcsec')
print(ellipse_Rmin, Rmin, ' arcsec')
print(ellipse_theta, theta, ' radians')

9.0074125081613 9.0074125081613  arcsec
4.374942013913388 4.374942013913388  arcsec
0.7931041204932631 0.7931041204932631  radians

Use photutils to define an aperture which enables visualizing the shape parameters on an image cutout. Note that photutils convention for pa or theta is 90 degrees or pi/2 offset from that of the LSST pipeline functions.

In [14]:

gaussell_ellipse = SkyEllipticalAperture(coord,
                                         Rmaj * arcsec_per_pix * u.arcsec,
                                         Rmin * arcsec_per_pix * u.arcsec,
                                         theta=(np.pi/2 + theta) * u.rad)

3.2. Kron ellipse¶

In the next cell, reconstruct the Kron aperture. The method of the Kron implementation in the LSST pipelines is to use the Object table parameter <band>_kronRad, in combination with the same axis-ratio and rotation angle that come from the Gaussian shape parameters from section 3.1.

First, define R_mom, which is the (circularized) moment radius defined as sqrt(Rmaj * Rmin) where Rmaj and Rmin came from the Gaussian ellipse measurement. Together these define the axis ratio. This will be combined with kronRad to obtain the semin-major and minor radii used for the kron aperture. kronRad in the `Object' table has units of pixel. A known issue is that the schema implies kronRad is unitless, but this is incorrect and the units listed in the schema will be fixed in DP2.

In [15]:

kronRad = tab['i_kronRad'][0]

R_mom = np.sqrt(Rmaj * Rmin)

kron_aperture_Rmaj = 2.5 * Rmaj * kronRad * arcsec_per_pix / R_mom
kron_aperture_Rmin = 2.5 * Rmin * kronRad * arcsec_per_pix / R_mom

Use photutils to define the Kron aperture.

In [16]:

kron_aperture = SkyEllipticalAperture(coord,
                                      kron_aperture_Rmaj * u.arcsec,
                                      kron_aperture_Rmin * u.arcsec,
                                      theta=(np.pi/2 + theta) * u.rad)

3.3. Sersic profiles¶

3.3.1. Sersic parameters in x-y¶

Convert the Sersic morphological parameters stored in the Object table (the effective Sersic radii in x and y directions, sersic_reff_x, sersic_reff_y, and the correlation coefficient from the multiband Sersic model fit sersic_rho that is related to orientation angle) into the more conventional sersic parameters (axis ratio ba defined as the ratio of the semi-minor half-light radius divided to semi-major half-light radius, the position angle pa, and the semi-major half-light radius r_major). The position angle (pa) convention is counter-clockwise relative to the x-axis. In this case, return the pa in units of degrees instead of the previous example in radians, to demonstrate its use with photutils.

In [17]:

x = EllipseMajor(Ellipse(sigma_x=tab['sersic_reff_x'][0],
                         sigma_y=tab['sersic_reff_y'][0],
                         rho=tab['sersic_rho'][0]), degrees=True)
ba = x.axrat
pa = x.angle
r_major = x.r_major * arcsec_per_pix

Use photutils to define the Sersic half-light ellipse.

In [18]:

sersic_ellipse = SkyEllipticalAperture(coord,
                                       r_major * u.arcsec,
                                       r_major * ba * u.arcsec,
                                       theta=(pa+90) * u.deg)

Generate an image cutout with edge size 0.008 degrees (~29 arcseconds).

In [19]:

cutout_size = 0.008

cutout = make_image_cutout(tab['coord_ra'][0], tab['coord_dec'][0],
                           cutout_size=cutout_size)

Plot the cutout and overplot the various shapes and apertures. This requires some manipulation of the WCS in order for lsst.images and astropy to plot consistently.

First, extract the WCS once using the lsst.images syntax. This replaces the old WCS(cutout.getWcs().getFitsMetadata()) from DP1. Then create a subplot with the extracted WCS.

In [20]:

astropy_wcs = cutout.sky_projection.as_fits_wcs(bbox=cutout.bbox)

plt.subplot(projection=astropy_wcs)

plt.imshow(cutout.image.array, origin='lower', cmap='gray', vmin=1,
           vmax=1000, norm='asinh', aspect='equal')

gaussell_pix_ellipse = gaussell_ellipse.to_pixel(astropy_wcs)
gaussell_pix_ellipse.plot(color=colors[0], lw=3, label='Gaussian (shape) ellipse')

kron_pix_aperture = kron_aperture.to_pixel(astropy_wcs)
kron_pix_aperture.plot(color=colors[1], lw=3, label='Kron Aperture')

sersic_pix_ellipse = sersic_ellipse.to_pixel(astropy_wcs)
sersic_pix_ellipse.plot(color=colors[2], lw=2, label='Sersic Half-light Profile')

plt.legend()

Out[20]:

<matplotlib.legend.Legend at 0x7be399602ba0>

Figure 1: An i-band image cutout of an example galaxy (grayscale). The Kron aperture used for Kron photometry is overplotted (green), as well as the half-light profile from the 2D Gaussian model (blue) and Sersic model (orange).

4. Sersic reconstruction¶

For some science applications it is useful to compare the on-sky image to the Sersic model, for example to assess goodness of fit. The cells below demonstrate how to reconstruct the best-fitting Sersic model that produced the <band>_sersicFlux, and represents the morphological parameters stored in the Object table.

A straightforward way to do this is using the galsim, which is the galaxy simulation package used to build the DP0.2 simulated images. Note that the algorithm used to measure the Sersic fluxes, multiprofit, is a Gaussian mixture model that approximates a Sersic model. Renderings of the galaxy using as a pure Sersic model using galsim is simpler to plot, but may not be perfectly identical to that built with Gaussian mixture.

It is also good to be aware of objects with Sersic index n that are measured to be suspiciously close to 1 (e.g. 1+/- 1e-6). This may happen for poorly fit objects because n=1 is the starting guess, and the optimizer tends to give up if it can't improve within 5-10 iterations.

Below, build the Sersic model using the galsim function gs.Sersic using the Sersic parameters from the cell above. Note that the total flux of the model very closely approximates the sersicFlux from the `Object' table, indicating (in this case) a good fit.

Note that models that exhibit broad wings in the light profile may require a larger cutout to reproduce the measured flux when summing model pixels (e.g. for cModel in Section 5.2).

4.1. Standard Sersic parameters¶

New in DP2, the Object table now contains standard Sersic parameters (semi-major and minor half-light radii, position angle) which can be used in place of the x-y parameters that were needed in DP1. For this section, as a demonstration, replace those calculated in Section 3.3.1 using the x-y Sersic parameters with the ones stored in the Object table.

In [21]:

r_major = tab['sersic_reff_major']
r_minor = tab['sersic_reff_minor']
ba = tab['sersic_reff_minor'] / tab['sersic_reff_major']
pa = tab['sersic_theta'][0] + 90

As above, generate the best fit 2D Sersic model of the galaxy. Note that galsim takes as input the circularized effective radius, which is defined as Reff = r_major * sqrt(ba).

In [22]:

sersic = gs.Sersic(n=tab['sersic_index'][0], half_light_radius=r_major * np.sqrt(ba),
                   flux=tab['i_sersicFlux'][0]).shear(q=ba, beta=pa * gs.degrees)
img_h, img_w = cutout.image.array.shape

sersic_model = sersic.drawImage(nx=img_w, ny=img_h, scale=arcsec_per_pix).array

print('Total flux of model = ', np.sum(sersic_model),
      'nJy, very close to flux from Object table = ',
      tab['i_sersicFlux'][0], ' nJy')

fig, ax = plt.subplots()
im = ax.imshow(sersic_model, origin='lower', interpolation='nearest',
               vmin=0, vmax=1000, norm='asinh')
cbar = fig.colorbar(im, ax=ax)
cbar.set_label('Flux [nJy/pixel]', rotation=270, labelpad=25)
cbar.set_ticks([-0.1, 10, 30, 60, 100, 1000])
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Total flux of model =  1.6613828e+06 nJy, very close to flux from Object table =  1.66207e+06  nJy

Figure 2: The best fit Sersic model of the galaxy in Figure 1 that was used to measure sersicFlux of the galaxy.

5. Fixed bulge and disk shapes¶

DP2 contains new shape measurements that can be used to reconstruct the light profiles that model galaxies with fixed Sersic index to n=1 (<f>__exponential_*) and the cModel (composite, or two component model). This is the sum of a bulge: deVaucouleurs profile with n=4 and a disk: exponential profile n=1.

5.1. Exponential profile¶

Below, generate the model for the best fitting exponential model (n=1). Note that since the best fit single Sersic index sersic_index n=2.2 in the Object table, fitting with a different Sersic index n=1 results in a more extended light distribution (compared to n=2.2, which is more concentrated), and as a result, a different measurement of total flux from that when Sersic index is left as a free parameter.

In [23]:

ba = tab['exponential_reff_minor'][0]/tab['exponential_reff_major'][0]
use_beta = (tab['exponential_theta'][0] + 90) * gs.degrees
sm = gs.Sersic(n=1, half_light_radius=tab['exponential_reff_major'][0] * np.sqrt(ba),
               flux=tab['i_exponentialFlux'][0]).shear(q=ba, beta=use_beta)

exp_model = sm.drawImage(nx=img_w, ny=img_h, scale=arcsec_per_pix).array

print('Total flux of model = ', np.sum(exp_model), 'nJy, very close to flux from Object table = ',
      tab['i_exponentialFlux'][0], ' nJy')

fig, ax = plt.subplots()
im = ax.imshow(exp_model, origin='lower', interpolation='nearest', vmin=0, vmax=1000, norm='asinh')
cbar = fig.colorbar(im, ax=ax)
cbar.set_label('Flux [nJy/pixel]', rotation=270, labelpad=25)
cbar.set_ticks([-0.1, 10, 30, 60, 100, 1000])
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Total flux of model =  1.3964601e+06 nJy, very close to flux from Object table =  1.39646e+06  nJy

Figure 3: The best fit exponential Sersic model of the galaxy (i.e. Sersic index n is fixed to one) for the galaxy in Figure 1. This model was used to measure exponentialFlux of the galaxy.

5.2. cModel Reconstruction¶

The cModel fluxes are measured using the meas_modelfit algorithm. It does an initial guess to identify the centroid and preliminary shape moments, and then models the galaxy as an exponential disk with n=1 (to generate exp_reff and expFlux), then models the galaxy as a bulge with a deVaucouleurs profile using n=4 (to generate dev_r, devFlux, then fits a model as a linear combination of those, where the amplitudes (or fraction of flux, fracDev) are left free, with the resulting output stored as cModel fluxes.

Thus, <f>_cModel_expFlux  and <f>_cModel_devFlux are total fluxes and do not simply add up to <f>_cModelFlux. One can use the <f>_cModel_fracDev (standing for fraction which is deVaucouleurs, of <f>_cModel_devFlux) to reconstruct the cModel profile. The fraction of flux in the exponential disk component is thus <f>_cModel_fracDev.

First, store the parameters of the two Sersic components for cModel.

In [24]:

frac_dev = tab['i_cModel_fracDev'][0]
total_flux = tab['i_cModelFlux'][0]

bulge_flux = total_flux * frac_dev
bulge_ba = tab['i_cModel_dev_reff_minor']/tab['i_cModel_dev_reff_major']
bulge_reff = tab['i_cModel_dev_reff_major'][0] * np.sqrt(bulge_ba)
bulge_pa = tab['i_cModel_dev_theta'][0] + 90

disk_flux = total_flux * (1.0 - frac_dev)
disk_ba = tab['i_cModel_exp_reff_minor']/tab['i_cModel_exp_reff_major']
disk_reff = tab['i_cModel_exp_reff_major'][0] * np.sqrt(disk_ba)
disk_pa = tab['i_cModel_exp_theta'][0] + 90

The cell below demonstrates how to reconstruct the intrinsic cModel light profile of the galaxy using these parameters.

In [25]:

disk = gs.Exponential(flux=disk_flux, half_light_radius=disk_reff)
disk = disk.shear(q=disk_ba, beta=disk_pa * gs.degrees)

bulge = gs.DeVaucouleurs(flux=bulge_flux, half_light_radius=bulge_reff)
bulge = bulge.shear(q=bulge_ba, beta=bulge_pa * gs.degrees)

cmodel_galaxy = disk + bulge

Next, have galsim generate an image of the best cModel shape.

Warning: The following cell produces a pink warning message that is safe to ignore. It appears because the de Vaucouleurs profiles have sharp central cusps and also large wings which ideally requires a large grid of frequencies for the Fast Fourier Transform (FFT). Because this demonstration is only visualizing the light profile, and is not concerned with re-calculating photometry with high accuracy, this warning can be safely ignored.

In [26]:

cmodel = cmodel_galaxy.drawImage(nx=img_w, ny=img_h, scale=arcsec_per_pix).array

/opt/lsst/software/stack/conda/envs/lsst-scipipe-12.3.0-exact/lib/python3.13/site-packages/galsim/errors.py:441: GalSimFFTSizeWarning: drawFFT requires a very large FFT.
The required FFT size would be 11390 x 11390, which requires 2.90 GB of memory.
If you can handle the large FFT and want to suppress this warning,
you may update gsparams.maximum_fft_size.
  warnings.warn(GalSimFFTSizeWarning(message, size))

Finally, visualize the model of the observations.

In [27]:

plt.figure(figsize=(6, 6))
plt.imshow(cmodel, origin='lower', interpolation='nearest', vmin=0,
           vmax=1000, norm='asinh')

print('Total flux of model = ', np.sum(cmodel),
      'nJy, very close to flux from Object table = ',
      tab['i_cModelFlux'][0], ' nJy')

plt.title(f"Reconstructed cModel (fracDev = {frac_dev})")
plt.colorbar(label='Pixel Flux')
plt.show()

Total flux of model =  2.1509728e+06 nJy, very close to flux from Object table =  2.29825e+06  nJy

Figure 4: Reconstructed light profile using the cModel bulge plus disk 2 component model.

6. Apply PSF convolution¶

Sections 4-5 visualized the intrinsic light profile (without the instrumental effect due to the PSF). Below, reconstruct the PSF at the location of the galaxy using the metadata stored in the deep_coadd image cutout, and convolve the intrinsic light profile with the PSF.

In [28]:

target_x = cutout.bbox.x.min + (cutout.bbox.x.max - cutout.bbox.x.min) / 2
target_y = cutout.bbox.y.min + (cutout.bbox.y.max - cutout.bbox.y.min) / 2

target_point = (target_x, target_y)

local_psf_image = cutout.psf.compute_kernel_image(x=target_x, y=target_y)
psf_array = local_psf_image.array

galsim_psf = gs.InterpolatedImage(gs.Image(psf_array), scale=arcsec_per_pix)

final_profile = gs.Convolve([cmodel_galaxy, galsim_psf])

Finally, visualize the intrinsic light profile of the galaxy after PSF convolution. Note that due to the larger wings of the disk component compared to other models, the relatively small box chosen for good visualization means that the model sum misses about ~6% of the flux stored in the Object table.

In [29]:

img_w_psf = final_profile.drawImage(nx=img_w, ny=img_h, scale=arcsec_per_pix).array

plt.figure(figsize=(6, 6))
plt.imshow(img_w_psf, origin='lower', interpolation='nearest', vmin=0, vmax=1000, norm='asinh')

print('Total flux of model = ', np.sum(img_w_psf), 'nJy, very close to flux from Object table = ',
      tab['i_cModelFlux'][0], ' nJy')

plt.title("Reconstructed cModel (convolved with PSF)")
plt.colorbar(label='Pixel Flux')
plt.show()

Total flux of model =  2.150725e+06 nJy, very close to flux from Object table =  2.29825e+06  nJy

Figure 5: cModel light profile convolved with the PSF evaluated at the location of galaxy in the deep_coadd.