303.1. Galaxy photometry

303.1. Galaxy photometry 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 photometry produced by the LSST pipelines and their applications.

LSST data products: Object table

Packages: lsst.rsp

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¶

Photometry is the measurement of how much light is apparent from astronomical sources. The amount of light arriving on the telescope from the object is typically referred to as the flux density, or apparent magnitude (depending on units). Flux density is defined as the amount of energy arriving on the telescope per unit area, per unit time, per unit frequency (or wavelength) of the light.

The LSST Science Pipelines makes a variety of photometric measurements for point-like and extended sources. This notebook will teach the user about the automated photometry measurements for extended sources that are measured on the deep_coadd images and appear in the object table as part of the LSST pipelines data products.

The photometry measurements in the catalogs are flux densities in units of nano-Jansky [nJy]. 1 Jy = 10^{-23} ergs/s/cm^2/Hz.

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 the RSPDiscovery module for accessing the Table Access Protocol (TAP) service.

Finally, import scipy.stats package.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from lsst.rsp import RSPDiscovery
from scipy.stats import binned_statistic

1.2. Define parameters and functions¶

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

In [2]:

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

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

2. Types of photometry¶

This section will explore photometry measurements produced by the LSST pipelines, and provide some guidance for which are optimal for various applications for science with galaxies.

2.1. Explore the schema¶

Numerous photometry measurements are produced by the LSST Pipelines. Two types of photometry are in the object table. The first are total fluxes (see Section 3), which aim to approximate (or model) all of the light coming from object. The second class of fluxes are measured inside an on-sky aperture but not corrected for flux that may fall outside: thus they are apparent fluxes but do not recover the intrisic (total) flux (see Section 4 and 5). The apparent fluxes are optimized for other purposes, such as for measuring accurate light profiles or accurate colors.

Descriptions for all the measurements in the object table can be found in the DP2 schema; the schema can also be obtained programmatically as demonstrated in the following code cells.

First, see what is available in the object table by querying the tap_schema columns, and printing all the parameters available related to "Flux" measured in the i-band (as an example). For clarity, the return also omits errors and flags associated with the photometric measurements outlined in section 1.1.

In [4]:

query = "SELECT column_name, datatype, description, unit " \
        "FROM tap_schema.columns " \
        "WHERE table_name = 'dp2.Object'"
results = service.search(query).to_table()

In [5]:

search_string = 'Flux'
band = 'i_'
exclude1 = 'Err'
exclude2 = 'flag'
for cname in results['column_name']:
    if cname.find(search_string) > -1 and cname.find(band) > -1 and \
       cname.find(exclude1) == -1 and cname.find(exclude2) == -1:
        print(cname)

i_psfFlux
i_free_psfFlux
i_cModel_devFlux
i_cModel_expFlux
i_gaapPsfFlux
i_gaap0p7Flux
i_gaap1p0Flux
i_kronFlux
i_ap03Flux
i_ap06Flux
i_ap09Flux
i_ap12Flux
i_ap17Flux
i_ap25Flux
i_ap35Flux
i_ap50Flux
i_ap70Flux
i_cModelFlux
i_cModelFlux_inner
i_free_cModelFlux
i_free_cModelFlux_inner
i_free_cModel_devFlux
i_free_cModel_expFlux
i_deblend_zeroFlux
i_psfFlux_area
i_exponentialFlux
i_sersicFlux

The object catalog also has pre-computed AB magnitudes (Mag columns) for cModel and psf fluxes. Query the tap_schema columns, and print all the parameters available related to Mag measured in the i-band.

In [6]:

search_string = 'Mag'
band = 'i_'
exclude1 = 'Err'
for cname in results['column_name']:
    if cname.find(search_string) > -1 and cname.find(band) > -1 and \
       cname.find(exclude1) == -1:

        print(cname)

i_psfMag
i_cModelMag

2.2. Select a galaxy sample¶

Query the DP2 object table for a selection of photometric measurements.

Limit the search to contain galaxies using the i_extendedness flag. This will exclude many point-like objects, but not all. While some contamination by point-like objects is tolerable for the purposes of this demonstration, users should consider the variety of extendedness measurements when deciding how to constrain samples for scientific analysis (see, e.g., the flag usage guidance).

Further, identify objects which have been detected at high signal to noise > 20 and whose photometric measurements have not been flagged as having an issue (i_kronFlux_flag or i_cModel_flag or sersic_no_data_flag = 0 means the photometry is ok). Also exclude very bright galaxies (i-band magnitude < 20).

Search for the sample using the DP2 imaging obtained in the Extended Chandra Deep Field South (ECDFS; center ra, dec = 53.2, -28.1 in degrees).

In [7]:

target_ra = 53.2
target_dec = -28.1
query = "SELECT obj.objectId, obj.coord_ra, obj.coord_dec, " + \
        "obj.detect_fromBlend, obj.detect_isIsolated, " + \
        "obj.i_blendedness, obj.i_extendedness, " + \
        "obj.i_kronFlux, obj.i_kronFluxErr, obj.i_kronRad, " + \
        "obj.i_cModelFlux, obj.i_cModelFluxErr, " + \
        "obj.i_cModel_devFlux, obj.i_cModel_expFlux, obj.i_cModel_fracDev, " + \
        "obj.i_gaap1p0Flux, obj.r_gaap1p0Flux, " + \
        "obj.sersic_index, " + \
        "obj.i_sersicFlux, obj.i_exponentialFlux, " + \
        "obj.i_kronFlux_flag, obj.i_cModel_flag, obj.sersic_no_data_flag, " + \
        "obj.i_ap03Flux, obj.i_ap06Flux, obj.i_ap09Flux, obj.i_ap12Flux, " + \
        "obj.i_ap17Flux, obj.i_ap25Flux, obj.i_ap35Flux, obj.i_ap50Flux " + \
        "FROM dp2.Object AS obj " + \
        "WHERE (obj.i_cModelFlux/obj.i_cModelFluxErr > 20) AND " + \
        "(obj.i_extendedness = 1) AND (obj.sersic_no_data_flag = 0) AND " + \
        "(obj.i_kronFlux_flag = 0) AND (obj.i_cModel_flag = 0) AND " + \
        "(scisql_nanojanskyToAbMag(obj.i_cModelFlux) > 20) AND " + \
        "CONTAINS(POINT('ICRS', obj.coord_ra, obj.coord_dec), " + \
        "CIRCLE('ICRS',"+str(target_ra)+","+str(target_dec)+", 0.1)) = 1 "

In [8]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
    job.raise_if_error()
assert job.phase == 'COMPLETED'
tab = job.fetch_result().to_table()

Job phase is COMPLETED

Print the results of the search query.

In [9]:

tab

Out[9]:

Table length=2448

objectId	coord_ra	coord_dec	detect_fromBlend	detect_isIsolated	i_blendedness	i_extendedness	i_kronFlux	i_kronFluxErr	i_kronRad	i_cModelFlux	i_cModelFluxErr	i_cModel_devFlux	i_cModel_expFlux	i_cModel_fracDev	i_gaap1p0Flux	r_gaap1p0Flux	sersic_index	i_sersicFlux	i_exponentialFlux	i_kronFlux_flag	i_cModel_flag	sersic_no_data_flag	i_ap03Flux	i_ap06Flux	i_ap09Flux	i_ap12Flux	i_ap17Flux	i_ap25Flux	i_ap35Flux	i_ap50Flux
	deg	deg					nJy	nJy	pix	nJy	nJy	nJy	nJy		nJy	nJy		nJy	nJy				nJy	nJy	nJy	nJy	nJy	nJy	nJy	nJy
int64	float64	float64	bool	bool	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	bool	bool	bool	float32	float32	float32	float32	float32	float32	float32	float32
755369504803917291	53.29995352431074	-28.0547425924061	False	True	0.0	1.0	1414.72	72.0451	4.94271	1403.65	38.7482	1647.6	1409.58	0.0	1150.62	803.32	0.998802	1405.98	1405.63	False	False	False	486.864	1016.58	1320.04	1336.43	1356.0	1537.43	2091.29	2839.97
755369504803917615	53.27852121568102	-28.041728904182452	False	True	0.0	1.0	3233.27	57.1106	3.88724	3356.16	37.6656	3742.26	3369.76	0.0	2883.62	2008.29	0.560724	3277.18	3333.87	False	False	False	1234.44	2631.58	3077.74	3254.31	3200.47	3216.93	3563.03	4032.17
755369504803917532	53.293283291734404	-28.04488854033139	False	True	0.0	1.0	2862.44	59.0203	4.03617	2737.52	32.9685	2807.54	2715.79	0.318602	2570.97	1629.27	1.37526	2725.83	2709.26	False	False	False	1244.24	2295.36	2630.83	2842.01	2848.95	2853.79	3625.39	4788.06
755369504803917500	53.28234741164852	-28.046407171944857	False	True	0.0	1.0	1322.73	67.146	4.42771	1232.63	32.7592	1264.39	1238.41	0.0	1182.68	559.705	1.0	1246.6	1246.6	False	False	False	589.831	1028.35	1232.82	1313.65	1401.64	1460.33	1725.96	2961.38
755369504803917506	53.28618008377495	-28.046274702255996	False	True	0.0	1.0	620.809	67.0139	4.51427	632.483	30.597	631.973	626.249	1.0	593.758	497.226	1.0	604.945	604.945	False	False	False	303.533	545.726	552.365	598.198	718.719	1037.21	1287.13	1584.13
755369573523419082	53.11010031308377	-28.056637623159684	True	False	-2.78063	1.0	563.165	36.124	2.7085	624.258	29.8603	643.001	626.239	0.0	562.47	415.416	1.0	581.709	581.709	False	False	False	276.948	504.797	552.537	607.75	541.965	673.904	392.952	780.458
755369573523419084	53.27540362944091	-28.056449925729464	True	False	-4.2554	1.0	1029.83	52.9737	3.42355	1019.68	33.4041	1037.58	1022.55	0.0	981.126	727.301	1.0	1052.26	1049.51	False	False	False	486.334	878.368	999.324	997.932	1069.73	1203.73	1532.87	2718.87
755369573523419086	53.264253287024864	-28.05643489446994	True	False	-7.21747	1.0	1302.45	227.223	14.5546	1166.24	33.0592	1176.58	1169.59	0.0	1093.95	658.694	1.0	1137.71	1137.71	False	False	False	562.651	1001.03	1065.43	986.574	992.808	1111.53	1216.97	1824.89
755369573523419072	53.26825865443503	-28.057639907755988	True	False	-5.03204	1.0	1057.88	52.191	3.36517	1063.47	33.556	1062.44	1040.46	1.0	1007.23	806.643	1.0	1044.11	1044.11	False	False	False	511.261	897.725	1024.83	1054.55	1000.83	1126.17	1602.54	2289.83
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
755368886328627624	53.210566565711	-28.188256343563943	False	True	0.0	1.0	816.844	69.6818	4.47262	814.582	36.89	866.124	781.559	0.429168	723.926	641.206	1.0	754.747	754.747	False	False	False	344.887	639.556	777.827	780.634	811.567	1026.0	1290.44	2384.96
755368886328627633	53.2346745806509	-28.18745262077458	False	True	0.0	1.0	8953.93	73.7362	4.68576	8822.35	42.4398	11111.3	8857.92	0.0	7263.55	3894.11	0.858384	8806.48	8903.12	False	False	False	3123.0	6690.56	8248.61	8671.76	8923.72	9156.67	9125.59	9784.35
755368886328627643	53.17988648404575	-28.187288276715655	False	True	0.0	1.0	4999.43	95.9461	6.08847	5586.85	65.5502	7906.7	5007.96	0.285344	3079.07	2655.66	1.0356	5040.16	5013.77	False	False	False	1070.05	2692.59	3902.09	4609.6	4800.91	4786.93	5343.2	5701.76
755368886328627476	53.20910825829933	-28.19454588576937	False	True	0.0	1.0	864.557	57.0972	3.69539	797.156	33.3584	814.509	793.016	0.307192	768.142	416.252	1.0	752.568	752.568	False	False	False	386.161	671.055	834.975	854.457	941.633	1005.54	1425.33	2756.37
755368886328627524	53.20838526719028	-28.192398129060187	False	True	0.0	1.0	2095.75	69.2284	4.45958	2127.89	37.5518	2350.58	2137.45	0.0	1845.79	1228.81	0.574232	2079.57	2106.74	False	False	False	842.635	1709.91	1933.09	2024.0	2097.3	2204.54	2447.99	3175.77
755368886328627515	53.211117287822674	-28.192754826530244	False	True	0.0	1.0	4278.46	80.579	5.15763	4399.4	53.8613	6290.01	4415.84	0.0	2780.08	973.613	0.505871	4165.57	4405.97	False	False	False	952.226	2484.3	3501.72	3958.39	4215.32	4356.83	4528.19	5175.47
755368886328627522	53.21756223930313	-28.192313405372573	False	True	0.0	1.0	2620.04	79.9185	5.08378	2674.47	51.7611	3803.18	2684.46	0.0	1824.85	808.034	0.592978	2613.27	2754.99	False	False	False	663.561	1650.63	2208.42	2482.83	2528.47	2781.84	3291.19	4010.59
755369642242894084	53.086708739954794	-28.10107353300232	True	False	-4.20846	1.0	1536.13	82.0766	5.477	1379.96	39.023	1452.1	1290.84	0.58103	1160.57	902.409	1.63947	1294.6	1268.82	False	False	False	508.727	1012.82	1253.1	1435.79	1450.24	1669.58	2068.82	2898.41
755369642242894228	53.086865385546545	-28.098331173104896	True	False	-28.6942	1.0	7911.47	69.861	4.63945	7865.07	41.0404	9612.65	7877.11	0.0	6601.18	3083.45	0.788995	7831.06	7941.01	False	False	False	2716.71	5977.5	7327.75	7666.57	7839.84	7984.3	8270.61	8535.5

Below, convert the fluxes (units $nJy$) extracted from the object table into AB magnitudes using: $m_{AB} = -2.5log( f_{nJy}) + 31.4$ (see e.g. AB Magnitudes Wikipedia page for the conversion between flux in Jy and AB magnitudes).

Warning: The following cell will produce warnings for invalid value encountered in log10, which happens if the source flux is negative. This occasionally happens from aperture photometry if the included pixels inside the aperture have negative values and can be safely ignored for this example. The log10 will return a NaN which can be filtered out later.

In [10]:

cmodel_mag = -2.50 * np.log10(tab['i_cModelFlux']) + 31.4
sersic_mag = -2.50 * np.log10(tab['i_sersicFlux']) + 31.4
cmodel_exp_mag = -2.50 * np.log10(tab['i_cModel_expFlux']) + 31.4
cmodel_dev_mag = -2.50 * np.log10(tab['i_cModel_devFlux']) + 31.4
exponential_mag = -2.50 * np.log10(tab['i_exponentialFlux']) + 31.4
ap06_mag = -2.50 * np.log10(tab['i_ap06Flux']) + 31.4
ap09_mag = -2.50 * np.log10(tab['i_ap09Flux']) + 31.4
ap17_mag = -2.50 * np.log10(tab['i_ap17Flux']) + 31.4
ap12_mag = -2.50 * np.log10(tab['i_ap12Flux']) + 31.4
ap35_mag = -2.50 * np.log10(tab['i_ap35Flux']) + 31.4
kron_mag = -2.50 * np.log10(tab['i_kronFlux']) + 31.4
gaap_mag = -2.50 * np.log10(tab['i_gaap1p0Flux']) + 31.4

/tmp/ipykernel_1047/3036620627.py:10: RuntimeWarning: invalid value encountered in log10
  ap35_mag = -2.50 * np.log10(tab['i_ap35Flux']) + 31.4
/tmp/ipykernel_1047/3036620627.py:11: RuntimeWarning: invalid value encountered in log10
  kron_mag = -2.50 * np.log10(tab['i_kronFlux']) + 31.4

3. Total fluxes¶

The following total flux measurements are available for galaxies from entries in the object table: total fluxes from a Sersic model (with shape parameters left free), from cModel (sum of bulge and disk Sersic components fitted to the galaxy); and bulge (de Vaucouleurs) and/or exponential disk fluxes (with Sersic index fixed to either n=4 or 1).

Sersic fluxes¶

This photometric measurement models all galaxies as a single Sersic profile and calculates its total flux according to the best fitting model. See the Sersic profile Wikipedia page for more information.

The best fit Sersic model is evaluated based on a model to all filters, and then the flux in each filter is calculated by integrating the galaxy light from the best fit model, assuming the best-fit Sersic shape parameters.

<f>_sersicFlux    : Flux from the final sersic fit. Forced on <f>-band.
<f>_sersicFluxErr : Uncertainty of <f>_sersicFlux
sersic_no_data_flag   : Failure flag for <f>_sersicFlux

The LSST pipeline package responsible for Sersic fluxes is called multiprofit. Visit the multiprofit documentation.

Exponential fluxes¶

This photometric measurement models all galaxies as a single Sersic profile with Sersic index n fixed to 1 (exponential disk) and calculates its total flux according to the best fitting model. See the Sersic profile Wikipedia page for more information.

The best fit exponential model is evaluated based on a model to all filters, and then the flux in each filter is calculated by integrating the galaxy light from the best fit model, assuming the best-fit Sersic shape parameters when n is fixed to 1.

<f>_exponentialFlux    : Flux from the final Sersic fit assuming n=1. Forced on <f>-band.
<f>_exponentialFluxErr : Uncertainty of <f>_exponentialFlux

Composite Model (CModel) fluxes¶

Similar in nature to the SDSS cModel photometry, these will be familiar to SDSS users. In short, it is the linear combination of the best fit exponential (disk or D; Sersic index n = 1) and de Vaucouleurs (bulge or B; Sersic index n = 4) profiles. Thus, cModel is a good compromise between the exponential fluxes that assume either a disk (n=1) or a bulge (n=4).

<f>_cModelFlux    : Flux from the final cmodel fit. Forced on <f>-band.
<f>_cModelFluxErr : Uncertainty of <f>_cModelFlux
<f>_cModel_flag   : Failure flag for <f>_cModelFlux

For most cases the "fixed" cModel photometry (i.e. the catalog entries listed above) are preferred to that measured with more degrees of freedom labeled <f>_free_cModelFlux. The difference is that the fixed ones above uses a reference band (recorded as refBand in the schema) where the galaxy is well detected to determine the other parameters, which are then fixed when fitting for the flux in the other bands. The _free_cModelFlux measurement allows all parameters to be free and independent in each filter. The fixed _cModelFlux measurements are generally recommended for galaxy science applications where total flux measurements are needed (e.g. for intrinsic luminosity or mass).

In the object table, a pre-computed AB magnitude (<f>_cModelMag) also exists for this measurement.

Exponential disk fluxes (assuming n=1)¶

These fluxes assume an exponential profile, and is measured on the reference band with a fixed n=1 Sersic index, using the meas_modelfit package. It is run as an initialization for the cModel output.

<f>_cModel_expFlux    : Flux from the cmodel fit with n=1. Forced on <f>-band.
<f>_cModel_expFluxErr : Uncertainty of <f>_cModel_expFlux

De Vaucouleurs bulge fluxes (assuming n=4)¶

These fluxes assume an de Vaucouleurs profile, and is measured on the reference band with a fixed n=4 Sersic index, using the meas_modelfit package. It is run as an initialization for the cModel output.

<f>_cModel_devFlux    : Flux from the cmodel fit with n=4. Forced on <f>-band.
<f>_cModel_devFluxErr : Uncertainty of <f>_cModel_devFlux

Below, store the bulge and disc components of the cModel 2-Sersic-component flux. These were previously stored in the DP1 schema as bdFluxB (bulge, or n=4) and bdFluxD (disk, or n=1), but are deprecated in DP2. They can be reconstructed using the i_cModel_fracDev value, which holds the fraction of the cModel flux from the de Vaucouleurs (n=4) component, to decompose.

In [11]:

tab['i_bdFluxB'] = tab['i_cModelFlux'] * tab['i_cModel_fracDev']
tab['i_bdFluxD'] = tab['i_cModelFlux'] * (1 - tab['i_cModel_fracDev'])
bdFluxD_mag = -2.50 * np.log10(tab['i_bdFluxD']) + 31.4
bdFluxB_mag = -2.50 * np.log10(tab['i_bdFluxB']) + 31.4

/tmp/ipykernel_1047/1438055475.py:3: RuntimeWarning: divide by zero encountered in log10
  bdFluxD_mag = -2.50 * np.log10(tab['i_bdFluxD']) + 31.4
/tmp/ipykernel_1047/1438055475.py:4: RuntimeWarning: divide by zero encountered in log10
  bdFluxB_mag = -2.50 * np.log10(tab['i_bdFluxB']) + 31.4

3.1. Comparing total fluxes¶

This section will make several plots that show how the different photometric measurements compare.

First, compare cModelFlux (two component bulge+disk Sersic model flux) to sersicFlux (single Sersic model flux with shape parameters free), cModel_devFlux (Sersic with n=4) and cModel_expFlux (Sersic with n=1).

In [12]:

fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1,
                              width_ratios=[0.8, 0.2], figsize=(10, 6))

bins = np.arange(16, 27, 1)

ylims = [-0.5, 0.5]

ax.plot(cmodel_mag, (cmodel_mag-cmodel_exp_mag), 's', alpha=.3,
        label='cModel-cModel_exp', color=colors[0])
x = cmodel_mag
y = (cmodel_mag-exponential_mag)
bin_mean, bin_edge, binnum = binned_statistic(x, y,
                                              statistic='median', bins=bins)
binctr = bin_edge[:-1] + (bin_edge[1:]-bin_edge[:-1])/2
ax.plot(binctr, bin_mean, color='k', lw=3, label=None, zorder=10)
ax.plot(binctr, bin_mean, color=colors[0],
        lw=2, label='bin median', zorder=11)

ax.plot(cmodel_mag, (cmodel_mag-cmodel_dev_mag), '^', alpha=.3,
        label='cModel-cModel_dev', color=colors[1])
x = cmodel_mag
y = (cmodel_mag-cmodel_dev_mag)
bin_mean, bin_edge, binnum = binned_statistic(x, y,
                                              statistic='median', bins=bins)
binctr = bin_edge[:-1] + (bin_edge[1:]-bin_edge[:-1])/2
ax.plot(binctr, bin_mean, color='k', lw=3, label=None, zorder=10)
ax.plot(binctr, bin_mean, color=colors[1], lw=2, label='bin median', zorder=11)


ax.plot(cmodel_mag, (cmodel_mag - sersic_mag), 'o', alpha=.3,
        label='cModel - sersic', color=colors[2])

x = cmodel_mag
y = (cmodel_mag-sersic_mag)
bin_mean, bin_edge, binnum = binned_statistic(x, y,
                                              statistic='median', bins=bins)
binctr = bin_edge[:-1] + (bin_edge[1:]-bin_edge[:-1])/2
ax.plot(binctr, bin_mean, color='k', lw=3, label=None, zorder=10)
ax.plot(binctr, bin_mean, color=colors[2], lw=2, label='bin median', zorder=11)

ax.axhline(0, linestyle='--')
ax.set_xlabel('cModel Magnitude')
ax.set_ylabel('cModel mag - other mag')
ax.set_ylim([-1, 1])
ax.legend()

ax2.hist((cmodel_mag-exponential_mag), orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 40), align='mid',
         histtype="step", color=colors[0], stacked=True, fill=False, label='exp')
ax2.hist((cmodel_mag-cmodel_dev_mag), orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 40), align='mid',
         histtype="step", color=colors[1], stacked=True, fill=False, label='dev')
ax2.hist((cmodel_mag-sersic_mag), orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 40), align='mid',
         histtype="step", color=colors[2], stacked=True, fill=False, label='Sersic')

ax2.set_ylim(ylims)
ax2.axhline(0, linestyle='--', color='k')
ax.axhline(0, linestyle='--', color='k')
ax.set_xlabel('i-band cModel magnitude [AB magnitude]')
ax.set_ylabel('Difference btwn cModel mag - Other mag [AB magnitude]')
ax.set_ylim(ylims)

ax.legend()
ax2.legend()

Out[12]:

<matplotlib.legend.Legend at 0x7e24e64c9590>

Figure 1: Left panel: comparison between cModel flux with three other total flux measurements: Sersic flux with Sersic index n left free (orange), cModel_exp flux with n fixed to 1 (exponential disk; blue) and cModel_dev with n fixed to 4 (de Vaucouleurs or bulge profile; green). Right panel: histogram of values where values near 0 indicate the measured values are comparable. Sersic and cModel fluxes are more comparable than with n=1 or 4 fixed which have more scatter and disagree more for brighter galaxies (which also tend to be larger on-sky).

3.2. Exponential disk fluxes¶

Two estimates of total flux based on modeling objects as exponential disks (Sersic model with Sersic index n = 1) are in the object table. They are calculated using different algorithms, but the plot below demonstrates that the fluxes are very similar.

In [13]:

fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(8, 6))

one2one = np.arange(0, 35000, 1)

ax.plot(tab['i_exponentialFlux'], tab['i_cModel_expFlux'], 'o', alpha=.3)

ax.plot(one2one, one2one, linestyle='--', color='k', label='1:1')
ax.set_xlabel('i_exponentialFlux')
ax.set_ylabel('i_cModel_expFlux')
plt.legend()

Out[13]:

<matplotlib.legend.Legend at 0x7e24e63116d0>

Figure 2: A comparison of the galaxy flux (modeled as an exponential disk with n=1) using the cModel algorithm (cModel_expFlux; y-axis) vs the multiprofit (exponentialFlux; x-axis). The results are very similar with little scatter indicating that the two methods provide similar fluxes. Note that Figure 1 shows that fixing the Sersic index to n=1 can underestimate the flux relative to the two component cModelFlux and the sersicFlux when the index is a free parameter, and may better agree with the actual light profile.

Below, compare the cModelFlux with the sum of its decomposed components (bulge and disk) that were calculated using the cModel_devFrac value stored in the object table.

In [14]:

plt.plot(tab['i_cModelFlux'], tab['i_bdFluxB'] + tab['i_bdFluxD'], '.', alpha=.3,
         color='blue')

plt.xlabel('i_cModelFlux')
plt.ylabel('cModel Bulge + Disk components')

Out[14]:

Text(0, 0.5, 'cModel Bulge + Disk components')

Figure 3: A comparison of the sum of the decomposed cModel bulge and disk components calculated using the fracDev parameter (y-axis) vs the composite cModel flux (x-axis) stored in the object table. The values are identical.

Below, store two shape parameters that will be useful for interpreting photometric apertures. First, the Kron Radius, which is a good proxy for the size of the galaxy light profile. Then, the Sersic index, which describes the shape of the light profile of galaxies.

In [15]:

i_kronRad = tab['i_kronRad']
sersic_index = tab['sersic_index']

Plot a comparison of the flux measurements when sersic index is fixed, as a function of measured sersic index when sersic index is left free.

In [16]:

fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1,
                              width_ratios=[0.8, 0.2], figsize=(10, 6))

bins = np.arange(0, 7, 0.25)

ylims = [-1.2, 1.2]

ax.plot(sersic_index, (sersic_mag-cmodel_exp_mag), 's', alpha=.3,
        label='sersic-cModel_exp', color=colors[0])
x = sersic_index
y = (sersic_mag-cmodel_exp_mag)
bin_mean, bin_edge, binnum = binned_statistic(x, y,
                                              statistic='median', bins=bins)
binctr = bin_edge[:-1] + (bin_edge[1:]-bin_edge[:-1])/2
ax.plot(binctr, bin_mean, color='k', lw=3, label=None, zorder=10)
ax.plot(binctr, bin_mean, color=colors[0], lw=2, label='bin median', zorder=11)

ax.plot(sersic_index, (sersic_mag-cmodel_dev_mag), '^', alpha=.3,
        label='sersic-cModel_dev', color=colors[2])
x = sersic_index
y = (sersic_mag-cmodel_dev_mag)
bin_mean, bin_edge, binnum = binned_statistic(x, y,
                                              statistic='median', bins=bins)
binctr = bin_edge[:-1] + (bin_edge[1:]-bin_edge[:-1])/2
ax.plot(binctr, bin_mean, color='k', lw=3, label=None, zorder=10)
ax.plot(binctr, bin_mean, color=colors[2], lw=2, label='bin median', zorder=11)

ax.axhline(0, linestyle='--')

ax.set_ylabel('Sersic mag - B or D mag')
ax.set_ylim([-1, 1])
ax.legend()


ax2.hist((sersic_mag-cmodel_exp_mag), orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 40), align='mid',
         histtype="step", color=colors[0], stacked=True, fill=False, label='exp')
ax2.hist((sersic_mag-cmodel_dev_mag), orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 40), align='mid',
         histtype="step", color=colors[2], stacked=True, fill=False, label='dev')

ax2.set_ylim(ylims)
ax2.axhline(0, linestyle='--', color='k')
ax.axhline(0, linestyle='--', color='k')
ax.set_xlabel('Sersic Index n')
ax.set_ylabel('Sersic mag - Other mag')
ax.set_ylim(ylims)
ax.set_xlim([0, 8])
ax.legend()
ax2.legend()

Out[16]:

<matplotlib.legend.Legend at 0x7e24e60820d0>

Figure 4: Comparison of the Sersic mag (flux measured with Sersic parameters left free) with the flux that is measured with Sersic index fixed to either n=1 (exp) or n=4 (dev) as a function of Sersic index n. Fixing the Sersic index can inject scatter in the flux relative to leaving it free, but the Sersic flux converges to that measured when fixed to a disk(bulge) when the measured sersic index is n=1(n=4).

4. Apparent fluxes¶

This section explores three types of apparent fluxes: Kron (elliptical apertures that typically includes more than 90% of intrinsic light), and aperture photometry (flux measured inside circles of varying size).

Here, apparent fluxes refers to photometry that are not corrected for flux outside of the measurement aperture (i.e. not corrected to be a total flux). These measurements have their own applications but should not be used to measure mass or luminosity.

Kron fluxes¶

A decent summary of Kron fluxes in the NED documentation. The aperture used for the fluxes is 2.5 x R1 where R1 is the luminosity weighted radius (also called "first moment"; Kron et al. 1980).

<f>_kronFlux      : Flux from Kron Flux algorithm. Measured on <f>-band.
<f>_kronFluxErr   : Uncertainty of <f>_kronFlux.
<f>_kronFlux_flag : Failure flag for <f>_kronFlux.

The Kron radius, <f>_kronRad, is also available. In this case of LSST pipeline output, the Kron flux is not corrected for light that is emitted outside of the Kron aperture. While in many cases it will collect the majority of light, it will not be as accurate as the cModel for science cases requiring total flux.

Aperture fluxes¶

This contains the enclosed flux inside a given aperture (and not corrected to total fluxes using an aperture correction that accounts for the flux falling outside the aperture). Fixed aperture size refers to the aperture radius in pixels.

<f>_ap<pix>Flux     : Flux within <pix>-pixel aperture. Forced on <f>-band.
<f>_ap<pix>FluxErr  : Uncertainty of <f>_ap<pix>Flux.
<f>_ap<pix>FluxFlag : Failure flag for <f>_ap<pix>Flux.

The apertures are 3, 6, 9, 12, 17, 25, 35, 50, and 70 pixels. In the column name, apertures are 03, 06, 09, 12, and so on. While aperture fluxes are not corrected for the loss outside the aperture, if the aperture size is much larger than the galaxy size then it will approximate the total flux of the galaxy. The general application of these measurements are for measuring radial profiles (see Section 4.2 below).

4.1. Comparing total to apparent fluxes¶

This section will make several plots that compare the cModel flux (which we take as the fiducial total flux, as commonly used for SDSS) to some of the apparent flux measurements. Kron, and aperture photometry are all measures of light within a fixed aperture.

Generally, magnitudes measured using aperture photometry in the LSST pipeline are fainter than those measured from cModel, because the fixed circular aperture systematically underestimates the flux in the galaxy wings (and the lost flux increases as the intrinsic size of the galaxy increases, e.g. as traced by the Kron radius).

First, compare cModel to Kron which should typically enclose 90% of the light.

In [17]:

fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1,
                              width_ratios=[0.8, 0.2], figsize=(10, 6))

bins = np.arange(16, 27, 1)

ylims = [-1.2, 1.2]

ax.plot(cmodel_mag, (cmodel_mag - kron_mag), 's', alpha=.3,
        label='cModel-Kron', color=colors[0])
x = cmodel_mag
y = (cmodel_mag-kron_mag)
bin_mean, bin_edge, binnum = binned_statistic(x, y,
                                              statistic='median', bins=bins)
binctr = bin_edge[:-1] + (bin_edge[1:]-bin_edge[:-1]) / 2
ax.plot(binctr, bin_mean, color='k', lw=3, label=None, zorder=10)
ax.plot(binctr, bin_mean, color=colors[0], lw=2, label='bin median', zorder=11)

ax.plot(cmodel_mag, (cmodel_mag-sersic_mag), 'o', alpha=.3,
        label='cModel - Sersic', color=colors[2])

x = cmodel_mag
y = (cmodel_mag-sersic_mag)
bin_mean, bin_edge, binnum = binned_statistic(x, y,
                                              statistic='median', bins=bins)
binctr = bin_edge[:-1] + (bin_edge[1:]-bin_edge[:-1])/2
ax.plot(binctr, bin_mean, color='k', lw=3, label=None, zorder=10)
ax.plot(binctr, bin_mean, color=colors[2], lw=2, label='bin median', zorder=11)

ax.axhline(0, linestyle='--')
ax.set_xlabel('cModel Magnitude')
ax.set_ylabel('cModel mag - Sersic mag')
ax.set_ylim([-1, 1])
ax.legend()


ax2.hist((cmodel_mag-sersic_mag), orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 50), align='mid',
         histtype="step", color=colors[2], stacked=True, fill=False, label='Sersic')
ax2.hist((cmodel_mag-kron_mag), orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 50), align='mid',
         histtype="step", color=colors[0], stacked=True, fill=False, label='Kron')

ax2.set_ylim(ylims)
ax2.axhline(0, linestyle='--', color='k')
ax.axhline(0, linestyle='--', color='k')
ax.set_xlabel('Kron Radius [i-band; pixels]')
ax.set_ylabel('cModel mag - Other mag')
ax.set_ylim(ylims)

ax.legend()
ax2.legend()

Out[17]:

<matplotlib.legend.Legend at 0x7e24e6167750>

Figure 5: A figure comparing the difference between cModel and Kron magnitudes, compared to the difference between cModel and Sersic magnitudes (left panel). Generally, both Kron and Sersic the measurements are in comparable agreement with cModel. The right panel shows the histogram of magnitude differences, demonstrating that there are not systematic offsets but slightly higher scatter from Kron with respect to cModel.

Below, explore how circular aperture photometry compares to cModel. Generally, magnitudes measured using aperture photometry in the LSST pipeline are fainter than those measured from cModel, because the fixed circular aperture systematically underestimates the flux in the galaxy wings (and the lost flux increases as the intrinsic size of the galaxy increases, e.g. as traced by the Kron radius).

In [18]:

fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1,
                              width_ratios=[0.8, 0.2], figsize=(10, 6))

bins = np.arange(2, 15, 1)
ylims = [-1.5, 1.5]

ax.plot(i_kronRad, (cmodel_mag-ap06_mag), '^', alpha=.3,
        label='6-pixel aperture', color=colors[0])
ax.plot(i_kronRad, (cmodel_mag-ap09_mag), 's', alpha=.3,
        label='9-pixel aperture', color=colors[1])
ax.plot(i_kronRad, (cmodel_mag-ap12_mag), 'o', alpha=.3,
        label='12-pixel aperture', color=colors[2])
ax.plot(i_kronRad, (cmodel_mag-ap17_mag), '.', alpha=.3,
        label='17-pixel aperture', color=colors[3])

ax2.hist((cmodel_mag-ap17_mag), edgecolor=colors[3], orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 40), align='mid',
         histtype="step", stacked=True, fill=False)
ax2.hist((cmodel_mag-ap12_mag), edgecolor=colors[2], orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 40), align='mid',
         histtype="step", stacked=True, fill=False)
ax2.hist((cmodel_mag-ap09_mag), edgecolor=colors[1], orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 40), align='mid',
         histtype="step", stacked=True, fill=False)
ax2.hist((cmodel_mag-ap06_mag), edgecolor=colors[0], orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 40), align='mid',
         histtype="step", stacked=True, fill=False)

ax2.set_ylim(ylims)
ax.axhline(0, linestyle='--', color='k')
ax2.axhline(0, linestyle='--', color='k')
ax.set_xlabel('Kron Radius [i-band; pixels]')
ax.set_ylabel('cModel mag - Aperture mag')
ax.set_ylim(ylims)
ax.set_xlim([2, 12])
ax.legend()

Out[18]:

<matplotlib.legend.Legend at 0x7e24e6017b10>

Figure 6: A comparison of the difference between cModel photometry and aperture photometry measured by the LSST pipelines for four different aperture sizes as a function of galaxy size (as measured using the Kron radius). The left panel shows the scatter plot of difference in photometry vs Kron radius, and the right panel shows a histogram of these values that demonstrate that larger aperture sizes have photometry that is closer to the cModel. The right panel shows histograms of the data in the left panel, where the colors indicate for the same data in each panel.

In [19]:

fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(7, 6))
ylims = [-1.2, 1.2]
bins = np.arange(16, 27, 1)

ax.plot(cmodel_mag, (cmodel_mag-ap06_mag), '^', alpha=.1,
        label='cModel 6-pix aperture', color=colors[0])
x = cmodel_mag
y = (cmodel_mag-ap06_mag)
bin_mean, bin_edge, binnum = binned_statistic(x, y,
                                              statistic='median', bins=bins)
binctr = bin_edge[:-1] + (bin_edge[1:]-bin_edge[:-1])/2
ax.plot(binctr, bin_mean, color=colors[0], lw=2, label='bin median', zorder=10)

ax.plot(cmodel_mag, (cmodel_mag-ap09_mag), 's', alpha=.1,
        label='cModel 9-pix aperture', color=colors[1])
x = cmodel_mag
y = (cmodel_mag-ap09_mag)
bin_mean, bin_edge, binnum = binned_statistic(x, y,
                                              statistic='median', bins=bins)
binctr = bin_edge[:-1] + (bin_edge[1:]-bin_edge[:-1])/2
ax.plot(binctr, bin_mean, color=colors[1], lw=2,
        label='bin median', zorder=10)

ax.plot(cmodel_mag, (cmodel_mag-ap12_mag), 'o', alpha=.1,
        label='cModel 12-pix aperture', color=colors[2])
x = cmodel_mag
y = (cmodel_mag-ap12_mag)
bin_mean, bin_edge, binnum = binned_statistic(x, y,
                                              statistic='median', bins=bins)
binctr = bin_edge[:-1] + (bin_edge[1:]-bin_edge[:-1])/2
ax.plot(binctr, bin_mean, color=colors[2], lw=2,
        label='bin median', zorder=10)

ax.plot(cmodel_mag, (cmodel_mag-ap17_mag), '.', alpha=.1,
        label='cModel 17-pix aperture', color=colors[3])
x = cmodel_mag
y = (cmodel_mag-ap17_mag)
bin_mean, bin_edge, binnum = binned_statistic(x, y,
                                              statistic='median', bins=bins)
binctr = bin_edge[:-1] + (bin_edge[1:]-bin_edge[:-1])/2
ax.plot(binctr, bin_mean, color=colors[3], lw=2,
        label='bin median', zorder=10)

ax.axhline(0, linestyle='--')
ax.set_xlabel('cModel Magnitude')
ax.set_ylabel('cModel mag - Aperture mag')
ax.set_ylim([-1, 1])
ax.set_xlim([20, 25])
ax.legend()

Out[19]:

<matplotlib.legend.Legend at 0x7e24e5edb390>

Figure 7: A similar comparison of the difference between cModel photometry and aperture photometry measured by the LSST pipelines for four different aperture sizes, this time as a function of cModel magnitude. Running median is included.

These two figures show that the aperture photometry typically under-estimates the flux relative to the total flux estimated using cModel. As expected, there is a general trend for larger apertures to get closer to the total flux from cModel for large galaxies (i.e. whose Kron Radius is larger). There is also a general trend for the aperture photometry to be less discrepant at fainter magnitudes, since faint galaxies tend to be small.

4.2. Application of aperture photometry: radial profile¶

A science application for the aperture photometry is easy visualization of the radial profile of galaxies. In the cell below, make this plot for both a large galaxy (first) and a smaller galaxy of similar brightness (second). The query looks for bright galaxies whose cModel magnitude ~ 20 ABmag. Dividing the aperture flux by the surface area of the aperture yields the surface brightness, which can be plotted as a function of radius from the center of the galaxy to compare radial light profiles.

In [20]:

wh = np.where((tab['i_kronRad'] > 20) & (cmodel_mag > 20)
              & (cmodel_mag < 21.5))[0]

indx = 0
arcsec_per_pix = 0.2

rad = np.array([3, 6, 9, 12, 17, 25, 35, 50]) * arcsec_per_pix
area = np.pi * rad**2
profile = np.array([tab['i_ap03Flux'][wh][indx], tab['i_ap06Flux'][wh][indx],
                    tab['i_ap09Flux'][wh][indx], tab['i_ap12Flux'][wh][indx],
                    tab['i_ap17Flux'][wh][indx], tab['i_ap25Flux'][wh][indx],
                    tab['i_ap35Flux'][wh][indx], tab['i_ap50Flux'][wh][indx]]) / area

plt.plot(rad, profile, linestyle=':',
         label='Large Radius R='
         + str(np.round(i_kronRad[wh][indx]*arcsec_per_pix, 2)))
plt.xlabel('Aperture Radius [arcsec]')
plt.ylabel(r'Surface Brightness [nJy arcsec$^{-2}$]')

wh2 = np.where((tab['i_kronRad'] < 8) & (tab['i_kronRad'] > 5)
               & (cmodel_mag > 20) & (cmodel_mag < 21.5))[0]

indx = 0
print("large galaxy mag = ", cmodel_mag[wh][indx], " small galaxy mag = ", cmodel_mag[wh2][indx])

profile = np.array([tab['i_ap03Flux'][wh2][indx], tab['i_ap06Flux'][wh2][indx],
                    tab['i_ap09Flux'][wh2][indx], tab['i_ap12Flux'][wh2][indx],
                    tab['i_ap17Flux'][wh2][indx], tab['i_ap25Flux'][wh2][indx],
                    tab['i_ap35Flux'][wh2][indx], tab['i_ap50Flux'][wh2][indx]])/area

plt.plot(rad, profile,
         label='Small Radius R='
         + str(round(i_kronRad[wh2][indx]*arcsec_per_pix, 2)))
plt.legend()
plt.yscale('log')

large galaxy mag =  20.943514251708983  small galaxy mag =  20.850939893722533

Figure 8: Plot demonstrating the use of aperture photometry to plot the surface brightness profile (as a function of aperture radius) for a galaxy with large Kron radius (green solid) and small Kron radius (blue dotted).

5. Photometry for color¶

This section will explore GaaP fluxes (which are optimized for measuring accurate colors between bands).

GaaP fluxes¶

These are the Gaussian-aperture-and-PSF flux that is defined in Kuijken et al. 2008. The main goal of this method is to measure accurate colors while accounting for the different spatial resolution between filters. This is sometimes achieved in other datasets by convolving all images to the largest PSF, but this process of PSF-matching is computationally very time consuming for large images, thus motivating GaaP as a faster alternative. It is not a measure of total flux in a filter. Several measurement apertures are available.

Aperture

<f>_gaap<ap>Flux    : GaaP flux with <ap> aperture after multiplying the seeing aperture. Forced on <f>-band.
<f>_gaap<ap>FluxErr : Uncertainty of <f>_gaap<ap>Flux.

Where the measurement apertures are 0.7 and 1.0 arcseconds. In the column name <ap> appears as 0p7 and 1p0. Multiplying by the "seeing aperture" refers to convolving the PSF with a kernel so that the PSF is as if the seeing were 1.15 arcseconds. This has the effect of smearing the images of all filters consistently so that the colors are accurate.

For photometric redshifts, and other analysis where accurate colors are important, it is recommended to start with the GaaP fluxes with 1.0 aperture, which was found to have better overall performance compared to other aperture sizes. Experiment yourself to see how it works for your science case.

5.1. Kron and GaaP comparison¶

In the next cell, compare the cModel instead to the Kron and GaaP measures.

In [21]:

fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1,
                              width_ratios=[0.8, 0.2], figsize=(10, 6))
ylims = [-2, 2]

ax.plot(i_kronRad, (cmodel_mag-gaap_mag), 's', alpha=.3,
        label='gaap1p0', color=colors[0])

ax.plot(i_kronRad, (cmodel_mag-kron_mag), 'o', alpha=.3,
        label='Kron', color=colors[2])

ax2.hist((cmodel_mag-gaap_mag), orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 40), align='mid',
         histtype="step", color=colors[0], stacked=True, fill=False, label='gaap1p0')
ax2.hist((cmodel_mag-kron_mag), orientation="horizontal",
         bins=np.linspace(ylims[0], ylims[1], 40), align='mid',
         histtype="step", color=colors[2], stacked=True, fill=False, label='Kron')

ax2.set_ylim(ylims)
ax.axhline(0, linestyle='--', color='k')
ax2.axhline(0, linestyle='--', color='k')
ax.set_xlabel('Kron Radius [i-band; pixels]')
ax.set_ylabel('cModel mag - Other mag')
ax.set_ylim(ylims)
ax.set_xlim([2, 15])
ax.legend()
ax2.legend()

Out[21]:

<matplotlib.legend.Legend at 0x7e24e5bd4410>

Figure 9: The left panel figure shows the i-band magnitude difference between cModel and Kron (orange circles) and between cModel and gaap1p0 (blue squares) vs the Kron radius (a proxy for galaxy size) for the galaxies in the query. The dashed line indicates where the two magnitudes would have the same value. The gaap1p0 magnitude always underestimates the flux, but the offset becomes worse for larger galaxies (relative to the fixed aperture). The right panel shows the histogram of the magnitude differences in the left panel, illustrating that while cModel - Kron magnitudes are similar on average (blue histogram) the gaap1p0 systematically underestimates the flux relative to cModel (orange histogram).

5.2. CMD with GaaP¶

This section demonstrates using GaaP photometry to calculate accurate galaxy colors to identify different types of galaxies. First, define magnitudes from g, r, and i band photometry, then compare the colors of galaxies that overlap the galaxy cluster with that in the field. In clusters, galaxies tend to be old, red elliptical galaxies and thus exhibit a well defined red sequence in color space.

The earlier query in Section 2 returned signal-to-noise ratio $>20$ galaxies from a blank, "field" location (the ECDFS) and stored them in the tab table. These will be dominated by bluer star forming galaxies which are most common in field environments.

Add a new query near a known galaxy cluster PSZ2 G309.43-72.86 at redshift z=0.35, from the ELAIS-S1 field, and use a signal-to-noise ratio $>50$ instead of $20$.

In [22]:

cluster_ra = 10.2082
cluster_dec = -44.1307
query = "SELECT obj.objectId, obj.coord_ra, obj.coord_dec, " + \
        "obj.detect_fromBlend, obj.detect_isIsolated, " + \
        "obj.i_blendedness, obj.i_extendedness, " + \
        "obj.i_kronFlux, obj.i_kronFluxErr, obj.i_kronRad, " + \
        "obj.i_cModelFlux, obj.i_cModelFluxErr, obj.i_gaap1p0Flux, " + \
        "obj.r_gaap1p0Flux, " + \
        "obj.i_kronFlux_flag, obj.i_cModel_flag " + \
        "FROM dp2.Object AS obj " + \
        "WHERE (obj.i_cModelFlux/obj.i_cModelFluxErr > 50) AND " + \
        "(obj.i_extendedness = 1) AND " + \
        "(obj.i_kronFlux_flag = 0) AND (obj.i_cModel_flag = 0) AND " + \
        "CONTAINS(POINT('ICRS', obj.coord_ra, obj.coord_dec), " + \
        "CIRCLE('ICRS',"+str(cluster_ra)+","+str(cluster_dec)+", 0.2)) = 1 "

In [23]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
    job.raise_if_error()
assert job.phase == 'COMPLETED'
tab2 = job.fetch_result().to_table()

Job phase is COMPLETED

First, calculate the magnitudes of galaxies in the 'field' location. Then, calculate the magnitudes for the other filters near the galaxy cluster from the query performed in the cell above. This will enable plotting their colors.

Warning: Like in Section 2, the following cell will produce warnings for invalid value encountered in log10, which happens if the source flux is negative. This happens for a small number of objects and since the goal of the plot is to see the distribution of the majority of sources, the warning can be safely ignored.

In [24]:

r_field_gaap_mag = -2.50 * np.log10(tab['r_gaap1p0Flux']) + 31.4
i_field_gaap_mag = -2.50 * np.log10(tab['i_gaap1p0Flux']) + 31.4
i_field_cmodel_mag = -2.50 * np.log10(tab['i_cModelFlux']) + 31.4
i_cluster_gaap_mag = -2.50 * np.log10(tab2['i_gaap1p0Flux']) + 31.4
r_cluster_gaap_mag = -2.50 * np.log10(tab2['r_gaap1p0Flux']) + 31.4
i_cluster_cmodel_mag = -2.50 * np.log10(tab2['i_cModelFlux']) + 31.4

In [25]:

fig, (ax, ax1) = plt.subplots(ncols=1, nrows=2, figsize=(10, 6), sharex=True)

ax.plot(i_field_cmodel_mag, (r_field_gaap_mag-i_field_gaap_mag),
        '.', alpha=.1, color='blue', label='Field Galaxies (ECDFS)')
ax.set_xlabel('i-band Magnitude [cModel]')
ax.set_ylabel('r-i color')
ax.set_ylim([-1, 2])
ax.legend()

ax1.plot(i_cluster_cmodel_mag, (r_cluster_gaap_mag-i_cluster_gaap_mag),
         '.', alpha=.1, color='r', label='Cluster Galaxies (PSZ2 G309.43-72.86; ELAIS-S1)')
ax1.set_xlabel('i-band Magnitude [cModel]')
ax1.set_ylabel('r-i color')
ax1.set_ylim([-1, 2])
ax1.legend()

Out[25]:

<matplotlib.legend.Legend at 0x7e24e5f57ed0>

Figure 10: The r − i vs. i color-magnitude diagram for galaxies selected in the queries. Top panel shows the SNR$>20$ galaxies selected from in a random field that does not contain a galaxy cluster (ECDFS). The bottom panel shows the SNR$>50$ galaxies from a field with a galaxy cluster, PSZ2 G309.43-72.86. The cluster galaxies appear as a "red sequence" with red r-i colors, because the Balmer / 4000 Angstrom break spectral feature that traces older stars sits between the bands. The faint-end cutoff is different for the two subsets due to the different SNR used in the query selection.

A very nice red sequence appears from the red, old galaxies in the cluster!

6. Exercise for the learner¶

Compare the <f>_free_cModelFlux measurements to <f>_cModelFlux in the filters that are not the reference band where the <f>_cModelFlux was measured (i.e. refBand). Investigate how leaving the cModel measurements free differs from the one measured with parameters fixed to the refBand, as a function of decreasing signal to noise. As an additional exercise, check how the signal to noise in colors measured using <f>_gaap1p0Flux values compare to those measured with the larger 3.0" aperture, <f>_gaap3p0Flux, where the larger aperture may increase the noise.