303.2. Galaxy Shapes#

303_2_Galaxy_Shapes.r29_2_0

303.2. Galaxy Shapes in DP1¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 1
Container Size: Large
LSST Science Pipelines version: r29.2.0
Last verified to run: 2025-09-02
Repository: github.com/lsst/tutorial-notebooks

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

LSST data products: objectTable, deepCoadd

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

Credit: Originally developed by Christina Williams and the Rubin Community Science team. This notebook is partly based on DP0.2 notebook 18 on galaxy photometry, and benefitted from earlier exploration of simulated data and notebook development by Melissa Graham and Dan Taranu, and helpful discussions with Jim Bosch.

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 deepCoadd images and appear in the objectTable as part of the LSST pipelines data products. The focus will be on galaxies. Data products related to shapes for the purpose of cosmological will be demonstrated elsewhere.

This notebook partly is based on a DP0.2 notebook (18) about galaxy photometry. Please note that the objectTable columns and contents may have changed between DP0.2 and DP1.

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).

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.

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.

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 lsst.geom as geom

from lsst.rsp import get_tap_service
from lsst.rsp.utils import get_pyvo_auth
from lsst.rsp.service import get_siav2_service

import lsst.afw.display as afwDisplay
from lsst.afw.image import ExposureF
from lsst.afw.fits import MemFileManager
import lsst.afw.geom.ellipses as ellipses

from pyvo.dal.adhoc import DatalinkResults, SodaQuery

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

from photutils.aperture import SkyEllipticalAperture

import galsim as gs
from lsst.gauss2d import Ellipse, EllipseMajor, Covariance

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 DP1 tutorial notebook 103.4 that demonstrates the Rubin image cutout service.

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
    -------
    mem : 'object'
        An lsst.afw.fits._fits.MemFileManager object containing the
        cutout image bytes returned by the cutout tool. The
        cutout image bytes can be accessed as an LSST awf image
        exposure class format such as: cutout = ExposureF(mem)
    """
    service = get_siav2_service("dp1")
    spherePoint = geom.SpherePoint(ra*geom.degrees, dec*geom.degrees)

    circle = (ra, dec, cutout_size)
    results = service.search(pos=circle, calib_level=3)
    table = results.to_table()
    tx = np.where((table['dataproduct_subtype'] == 'lsst.deep_coadd')
                  & (table['lsst_band'] == band))[0]

    datalink_url = results[tx[0].item()].access_url
    dl = DatalinkResults.from_result_url(datalink_url, session=get_pyvo_auth())

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

    sq.circle = (spherePoint.getRa().asDegrees() * u.deg,
                 spherePoint.getDec().asDegrees() * u.deg,
                 cutout_size * u.deg)
    cutout_bytes = sq.execute_stream().read()
    mem = MemFileManager(len(cutout_bytes))
    mem.setData(cutout_bytes, len(cutout_bytes))

    return mem

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')

Get an instance of the TAP service, and assert that it exists.

In [5]:

service = get_tap_service("tap")
assert service is not None

2. Find galaxies in ECDFS¶

The deepest DP1 dataset is in the Extended Chandra Deep Field South (ECDFS). This section will identify a galaxy sample in ECDFS to inspect morphologies.

In [6]:

target_ra = 53.195
target_dec = -27.703

Here, query the DP1 objectTable for galaxies in that field, identified as being extended (i_extendedness = 1). Additionally, check that their i_kronFlux_flag_*, sersic_no_data_flag, and shape_flag are all = 0, indicating there is not a problem with modeling their shapes and light profiles. Finally, make sure to identify galaxies that are well detected, with signal to noise (i_sersicFlux / i_sersicFluxErr > 20).

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_sersicFlux, obj.i_bdFluxB, obj.i_bdFluxD, " + \
        "obj.i_bdE1, obj.i_bdE2,obj.i_bdReB, obj.i_bdReD, " + \
        "obj.i_kronRad, obj.i_kronFlux_flag, obj.sersic_no_data_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.sersic_index, obj.sersic_reff_x, " + \
        "obj.sersic_reff_y, obj.sersic_rho " + \
        "FROM dp1.Object AS obj " + \
        "WHERE (obj.i_sersicFlux/obj.i_sersicFluxErr > 20) AND " + \
        "(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.01)) = 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=45

objectId	coord_ra	coord_dec	i_extendedness	refBand	i_kronFlux	i_sersicFlux	i_bdFluxB	i_bdFluxD	i_bdE1	i_bdE2	i_bdReB	i_bdReD	i_kronRad	i_kronFlux_flag	sersic_no_data_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	sersic_index	sersic_reff_x	sersic_reff_y	sersic_rho
	deg	deg			nJy	nJy	nJy	nJy	pix2	pix2	pix2	pix2						pix2	pix2	pix2		pix2	pix2	pix2		pix	pix
int64	float64	float64	float32	str1	float32	float32	float32	float32	float32	float32	float32	float32	float32	bool	bool	bool	bool	float32	float32	float32	bool	float32	float32	float32	float32	float32	float32	float32
611255759837093254	53.199663576429884	-27.69672467773119	1.0	i	9167.93	8797.69	10215.9	8741.72	3.22591	0.133742	2.23074	1.9165	4.30207	False	False	False	False	7.88048	0.418324	7.83588	False	7.88048	0.418324	7.83588	1.53463	1.26448	1.26892	0.121253
611255759837093257	53.2013027127527	-27.69803467521063	1.0	g	1011.95	847.968	821.272	735.639	5.28265	0.516451	1.79221	1.70645	5.69154	False	False	False	False	10.2205	1.00729	6.91665	False	10.5533	2.55631	6.03135	5.23214	1.77874	0.891193	0.224274
611255759837093256	53.199323537028725	-27.697616273641266	1.0	i	1433.41	1478.56	1254.47	1172.93	1.65818	0.506255	1.19915	1.32481	4.78868	False	False	False	False	7.2517	1.31589	5.66316	False	7.2517	1.31589	5.66316	5.43043	1.56729	1.11669	0.266155
611255759837093255	53.20079988330049	-27.69771002778915	1.0	i	2263.07	2051.1	2711.6	2145.78	4.65769	-0.270504	3.47921	2.45467	4.90205	False	False	False	False	9.15616	-1.22335	10.538	False	9.15616	-1.22335	10.538	0.54759	1.85628	1.88304	-0.532391
611255759837091299	53.189932625363795	-27.70316352393926	1.0	i	1235.88	1262.88	1127.9	933.024	9.50846	-0.0530412	3.0772	2.45304	6.63315	False	False	False	False	10.871	-0.691796	8.93709	False	10.871	-0.691796	8.93709	5.12773	2.81826	2.33318	-0.0995091
611255759837091240	53.19576102156065	-27.712042926357224	1.0	i	831.393	806.686	998.433	838.846	7.41561	0.708229	2.76345	2.1667	5.54082	False	False	False	False	11.1099	4.43951	8.95225	False	11.1099	4.43951	8.95225	0.588713	2.36954	1.40759	0.807632
611255759837091245	53.19825585131647	-27.708810591163374	1.0	i	1701.99	1567.08	1240.8	1013.14	16.469	-0.0973306	4.00825	3.04401	9.02746	False	False	False	False	14.6072	-1.90237	12.1894	False	14.6072	-1.90237	12.1894	4.39	4.23944	4.02601	-0.100785
611255759837091244	53.1939968968357	-27.705295910713442	1.0	i	2584.96	2761.34	2573.75	2044.85	11.7688	0.0275334	3.42156	2.47198	6.38867	False	False	False	False	10.5125	0.272488	9.15349	False	10.5125	0.272488	9.15349	4.06828	2.8193	2.57095	-0.199753
611255759837091270	53.19321526934131	-27.70681303285501	1.0	i	1802.15	1504.65	1209.3	1001.39	9.96979	-0.0680188	3.01336	2.39373	7.93978	False	False	False	False	11.0714	-0.808657	8.42779	False	11.0714	-0.808657	8.42779	4.51188	2.70914	3.67353	-0.0963647
611255759837091205	53.20166007949987	-27.701457020887304	1.0	i	2228.36	2155.86	2103.92	2033.06	0.433463	-0.105842	0.932089	1.18611	4.24703	False	False	False	False	6.1122	-0.182979	5.48326	False	6.1122	-0.182979	5.48326	5.73555	0.678497	0.574767	-0.0677738
611255759837091207	53.20403115758294	-27.700924965271117	1.0	i	2812.03	2877.05	3404.89	2624.16	13.7445	-0.274521	3.7557	2.54334	4.98968	False	False	False	False	11.7781	-1.19316	8.38531	False	11.7781	-1.19316	8.38531	1.62787	2.52214	1.69002	-0.379808
611255759837091208	53.20216324196284	-27.70787640593279	1.0	i	2750.79	2653.0	2942.86	2389.01	1.59789	-0.0413736	2.88718	2.22512	5.77308	False	False	False	False	6.57994	-0.145237	14.7968	False	6.57994	-0.145237	14.7968	2.61329	0.888863	3.05821	-0.358154
611255759837091211	53.191206335750934	-27.70212017723286	1.0	i	2635.29	2546.95	2835.09	2473.43	6.92824	0.128161	1.93105	1.74826	4.38752	False	False	False	False	10.3695	0.403928	5.58871	False	10.3695	0.403928	5.58871	1.32442	1.90727	0.807838	0.19795
611255759837091202	53.18817754735505	-27.696451896379912	1.0	i	2644.55	2699.1	2384.24	2190.11	0.962495	-0.100553	1.35535	1.42809	4.95363	False	False	False	False	6.39635	-0.25593	6.741	False	6.39635	-0.25593	6.741	5.5419	1.24572	1.25692	-0.220287
611255759837091210	53.20077589572504	-27.708360254856945	1.0	i	4023.21	5116.39	3932.03	2815.29	25.5751	-0.729051	5.24711	3.09745	8.27158	False	False	False	False	12.7095	-6.20462	16.4689	False	12.7095	-6.20462	16.4689	5.17289	5.64379	7.37058	-0.613854
611255759837091209	53.193020348877866	-27.711318143993562	1.0	i	2101.38	2264.62	1660.17	1600.93	1.64673	-0.385424	0.885352	1.1249	5.30923	False	False	False	False	6.9008	-0.602611	4.8183	False	6.9008	-0.602611	4.8183	5.9926	1.85742	1.34364	-0.180917
611255759837091194	53.19799954320962	-27.704597266052158	1.0	i	3737.29	3825.56	3891.94	3142.44	20.4077	-0.380813	3.03848	2.29698	5.6742	False	False	False	False	12.9795	-1.69005	7.06923	False	12.9795	-1.69005	7.06923	3.63462	2.98284	1.51729	-0.475803
611255759837091193	53.19976362399971	-27.700561117315264	1.0	i	3147.25	2916.28	3409.5	2944.75	5.19667	-0.33668	2.12727	1.87133	4.79487	False	False	False	False	8.78622	-0.904973	6.91885	False	8.78622	-0.904973	6.91885	1.65124	1.46056	0.962814	-0.520299
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
611255759837091216	53.19507671290735	-27.711217496615912	1.0	i	1305.49	1435.28	1172.62	1128.07	0.396632	0.00323104	0.895857	1.1162	4.62004	False	False	False	False	5.91342	0.350328	5.20046	False	5.91342	0.350328	5.20046	5.91035	1.18428	1.12088	0.00308204
611255759837091215	53.19634762555406	-27.7107203625018	1.0	i	1566.39	1769.24	1537.82	1371.31	3.66132	-0.132667	1.82193	1.7389	5.2346	False	False	False	False	8.80965	-0.442735	6.27896	False	8.80965	-0.442735	6.27896	4.93977	2.00956	1.50471	0.202759
611255759837091225	53.203033312936526	-27.702650825367492	1.0	i	1206.18	1124.58	1166.45	1068.11	1.33577	0.136562	1.54009	1.60282	5.21439	False	False	False	False	6.38554	0.543231	7.69054	False	6.38554	0.543231	7.69054	4.40214	0.780963	1.10066	0.0919333
611255759837091223	53.19133905504191	-27.69734067771905	1.0	i	1753.01	1820.87	1719.5	1426.47	8.02068	-0.411969	2.85118	2.27556	6.1105	False	False	False	False	9.25644	-1.70977	9.39643	False	9.25644	-1.70977	9.39643	3.32987	2.33693	2.31436	-0.350048
611255759837091222	53.19158483124395	-27.697974457019285	1.0	i	1763.13	1727.69	1967.83	1595.72	8.65264	-0.253494	3.1391	2.37408	4.97156	False	False	False	False	9.49625	-1.07268	9.41265	False	9.49625	-1.07268	9.41265	2.9407	1.6089	1.74039	-0.179497
611255759837091221	53.19088191344716	-27.70520162742041	1.0	i	1666.35	1564.06	1395.24	1195.52	4.24345	-0.00805663	2.18096	1.88244	6.34685	False	False	False	False	8.01557	-0.401021	7.63052	False	8.01557	-0.401021	7.63052	5.61203	1.97106	1.67376	0.0392122
611255759837091220	53.196577692737144	-27.700127927591033	1.0	i	5374.03	5975.73	4726.01	3480.51	31.3068	0.278838	4.66122	2.92068	7.81654	False	False	False	False	14.0152	1.46378	9.8381	False	14.0152	1.46378	9.8381	5.67441	5.22565	4.30789	0.204635
611255759837091183	53.19052691345921	-27.70369636345784	1.0	i	9130.47	9303.03	8096.21	6980.75	2.87209	-0.362452	2.10526	1.87771	5.73355	False	False	False	False	7.63359	-1.01733	8.84075	False	7.63359	-1.01733	8.84075	5.32051	1.67143	2.38697	-0.494772
611255759837091182	53.19607892095075	-27.697162722948054	1.0	i	8269.5	8083.6	8967.91	7491.05	2.23681	0.614162	2.25172	1.83268	5.27515	False	False	False	False	6.92381	3.49071	13.9592	False	6.92381	3.49071	13.9592	1.70049	1.29916	2.70297	0.772744
611255759837091181	53.186052360179254	-27.70077226911401	1.0	i	9066.08	8476.79	9428.47	8317.47	4.65767	0.393145	1.82218	1.69977	4.61619	False	False	False	False	9.0416	1.098	6.29213	False	9.0416	1.098	6.29213	1.13808	1.6107	1.01877	0.381051
611255759837091179	53.19957726312809	-27.709189233462414	1.0	i	16309.7	16529.3	20388.8	14616.0	28.0596	-0.296493	5.21768	3.07209	6.01311	False	False	False	False	14.0238	-1.68532	10.7514	False	14.0238	-1.68532	10.7514	2.09223	2.79632	2.27321	-0.284873
611255759837091190	53.18594169632371	-27.70744166959267	1.0	i	2984.23	2919.6	2953.71	2804.97	0.862144	-0.53797	1.05166	1.25984	3.70287	False	False	False	False	6.16282	-0.980656	6.35813	False	6.16282	-0.980656	6.35813	4.33012	0.659557	0.877332	-0.7304
611255759837091189	53.185571541256586	-27.705933603375957	1.0	i	2980.03	2885.69	2916.61	2896.18	-0.00286265	0.0678388	0.548499	0.812523	3.72316	False	False	False	False	5.07739	0.0436112	4.83607	False	5.07739	0.0436112	4.83607	5.4931	0.321514	0.297506	-0.138061
611255759837091188	53.19394763466494	-27.70471031560458	1.0	i	6716.69	6972.46	5779.77	5243.5	1.54367	-0.0102744	1.47309	1.50218	5.23287	False	False	False	False	7.04307	-0.00160825	6.65124	False	7.04307	-0.00160825	6.65124	5.95321	1.42767	1.66664	-0.0900605
611255759837091185	53.20523369884081	-27.704637091875746	1.0	i	5709.96	5497.5	5892.75	5569.11	0.643645	0.626311	1.01278	1.1938	3.8653	False	False	False	False	5.94876	1.15018	6.73427	False	5.94876	1.15018	6.73427	0.962335	0.861126	1.06969	0.617976
611255759837091178	53.190008149385	-27.69768546880768	1.0	i	15401.8	16105.1	15066.3	12206.6	6.95458	-0.268255	2.93516	2.2201	5.89136	False	False	False	False	8.75571	-0.88721	9.326	False	8.75571	-0.88721	9.326	4.48714	2.23955	2.34583	-0.125619
611255759837091177	53.197279286215014	-27.69603331458278	1.0	i	15319.0	14810.5	14363.7	14265.7	0.123316	0.0583554	0.300789	0.436867	3.51895	False	False	False	False	5.10946	0.052575	3.96201	False	5.10946	0.052575	3.96201	5.99511	0.462065	0.0479186	0.18464
611255759837091175	53.194881936221606	-27.703425462481054	1.0	i	1545340.0	1707590.0	2289690.0	1268050.0	559.808	0.766072	19.1187	7.37653	12.6523	False	False	False	False	47.8407	29.8852	48.1299	False	47.8407	29.8852	48.1299	2.31913	10.006	10.1266	0.791414
611255759837073869	53.19386252386268	-27.695489347891403	1.0	i	1073.07	1250.06	948.937	862.971	2.09955	0.804063	1.38	1.40022	5.18558	False	False	False	False	7.58669	4.22837	10.0594	False	7.58669	4.22837	10.0594	5.79553	1.8468	3.37274	0.929968

3. Visualize galaxy sizes/shapes¶

This section will make an image cutout of a 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 objectTable.

Below, use a galaxy that has been pre-selected in ECDFS to have large size to visualize, and define the SkyCoord for its location.

In [10]:

wh = np.where((tab['objectId'] == 611255759837091175))[0]
indx = 0

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

The sizes in the LSST objectTable 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 objectTable 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 deepCoadd 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 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'][wh][indx],
                                         tab['shape_yy'][wh][indx],
                                         tab['shape_xy'][wh][indx]))
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"][wh][indx],
                                  sigma_y_sq=tab["shape_yy"][wh][indx],
                                  cov_xy=tab["shape_xy"][wh][indx]))

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

print(ellipse_Rmaj, Rmaj)
print(ellipse_Rmin, Rmin)
print(ellipse_theta, theta)

8.824446180695961 8.824446180695961
4.254380064377689 4.254380064377689
0.7878174326516222 0.7878174326516222

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 objectTable 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 objectTable 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'][wh][indx]

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 shape¶

Convert the Sersic morphological parameters stored in the objectTable (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'][wh][indx],
                         sigma_y=tab['sersic_reff_y'][wh][indx],
                         rho=tab['sersic_rho'][wh][indx]), 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

mem = make_image_cutout(tab['coord_ra'][wh][indx],
                        tab['coord_dec'][wh][indx], cutout_size=cutout_size)

cutout = ExposureF(mem)

Plot the cutout and overplot the various shapes and apertures.

In [20]:

plt.subplot(projection=WCS(cutout.getWcs().getFitsMetadata()))
extent = (cutout.getBBox().beginX, cutout.getBBox().endX,
          cutout.getBBox().beginY, cutout.getBBox().endY)

plt.imshow(cutout.image.array, extent=extent, origin='lower', cmap='gray',
           vmin=0.1, vmax=500, norm='asinh')

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

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

sersic_pix_ellipse = sersic_ellipse.to_pixel(WCS(cutout.getWcs().getFitsMetadata()))
sersic_pix_ellipse.plot(color=colors[2], lw=2, label='Sersic Half-light Profile')
plt.legend()

Out[20]:

<matplotlib.legend.Legend at 0x7b289992b5f0>

No description has been provided for this image

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 objectTable.

A straitforward 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 objectTable, indicating (in this case) a good fit.

In [21]:

sersic = gs.Sersic(n=tab['sersic_index'][wh][indx], half_light_radius=r_major * np.sqrt(ba),
                   flux=tab['i_sersicFlux'][wh][indx]).shear(q=ba, beta=pa * gs.degrees)

img = sersic.drawImage(nx=500, ny=500, method="real_space", scale=arcsec_per_pix).array

print('Total flux of model = ', np.sum(img), 'nJy, very close to flux from objectTable = ',
      tab['i_sersicFlux'][wh][indx], ' nJy')

fig, ax = plt.subplots()
im = ax.imshow(img, 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 =  1707512.4 nJy, very close to flux from objectTable =  1707590.0  nJy

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

5. Known Issues¶

5.1 BD Models¶

BD models (where B and D refers to bulge and disk) produce measurements similar to those generated from modeling a single Sersic profile. The models are then integrated to estimate the total flux <f>_bdFluxB and <f>_bdFluxD. But, rather than using the best fit sersic index obtained from modeling the 6 images in each filter, the BD fluxes are measured from a model with the sersic index n fixed to either 1 (exponential disk in the case of D) or 4 (de Vaucouleurs bulge in the case of B). Note that this is not an equal decomposition of the two component sersic modeling as part of cModel.

There are known issues related to these parameters:

The units in the DP1 schema are wrong for <f>_bdReD and <f>_bdReB (where <f>_ indicates a measurement exists for each of the 6 filters). These are the effective radii measured on each filter image and the units are pixel.
The units in the DP1 schema are wrong for ellipticity parameters <f>_bdE1 and <f>_bdE2.
A bug currently affects the objectTable values for <f>_bdE1. In the meantime this parameter should not be used.
A note here for clarity: BD model shape parameters are not cataloged in the same way as the metrics for Sersic models. In the future, galsim.Shear can be used to convert the ellipicity parameters<f>_bdE1 and <f>_bdE2 to extract the shape parameters, but awaits the bug fix to <f>_bdE1 before it can be used.

It is currently recommended to not use any BD parameters for science with DP1.

5.2 Kron¶

A known issue for Kron shapes exists.

The units in the DP1 schema is wrong for Kron radius. The unit is pixel.