302.2. Stellar multi-color photometry

302.2. Stellar multi-color photometry¶

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: 2026-03-19
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20

Learning objective: Explore stellar color-magnitude and color-color diagrams in the Fornax dSph field.

LSST data products: Object

Packages: lsst.rsp.get_tap_service

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¶

This tutorial demonstrates the analysis of Rubin DP1 stellar photometry within the Fornax dSph field, utilizing color-magnitude diagrams (CMDs) and color-color Diagrams (CCDs) to characterize resolved stellar populations and their spatial distributions.

Constructing a Hertzsprung–Russell (H–R) diagram, by plotting stellar luminosity or absolute magnitude against surface temperature or spectral type, reveals distinct evolutionary phases, ranging from the hydrogen-burning main sequence to the evolved giant branches. While these fundamental parameters are traditionally derived via spectroscopy, astronomers leverage multi-band photometry to approximate these physical quantities through observed magnitudes and colors. Analyzing CMDs and CCDs enables the efficient characterization of large, resolved stellar populations and provides a robust framework for age, metallicity, and line-of-sight dust determination.

The Fornax dSph is a known satellite of the Milky Way. The LSSTComCam field targeting this galaxy is centered at (RA, Dec) = (40.080, -34.450) degrees and is representative of relatively crowded regions that LSST will encounter in the Milky Way.

Related tutorials: See the 300-level DP1 tutorials on stellar PSF photometry, the Fornax dSph field exploration, and the Milky Way foreground dust correction.

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 pacakge (pipelines.lsst.io), import the RSP Table Access Protocol (TAP) service and the DustValues class from the rubin_sim package (rubin_sim).

In [1]:

import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from matplotlib.patches import Polygon, Rectangle, Ellipse, Patch
from matplotlib.path import Path

from lsst.rsp import get_tap_service
from rubin_sim.phot_utils import DustValues

1.2. Define parameters and functions¶

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

In [2]:

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

Define a circular region with a 1-degree radius centered on the Fornax dSph field. Coordinates are in degrees.

In [3]:

ra_cen = 40.080
dec_cen = -34.450
radius = 1.0

Set the environment variable RUBIN_SIM_DATA_DIR to /rubin/rubin_sim_data to make the current rubin_sim throughput data available.

In [4]:

os.environ['RUBIN_SIM_DATA_DIR'] = '/rubin/rubin_sim_data'

List bands that are used in this notebook.

In [5]:

bands = "gri"

Set font sizes globally for all subsequent matplotlib figures in this notebook.

In [6]:

plt.rcParams["font.size"] = 15
plt.rcParams["figure.titlesize"] = 15
plt.rcParams["legend.fontsize"] = 12
plt.rcParams["axes.titlesize"] = 12

2. Construct CMD¶

Retrieve PSF photometry for stars in the field, apply Galactic dust extinction corrections, and evaluate the impact on ($g-i$, $i$) CMDs.

2.1. Query for point-like objects¶

Query the Object table in the Fornax dSph field for coordinates, PSF photometry in $gri$, and $E(B-V)$ values. Restrict to point sources (refExtendedness = 0) to only retrieve stars.

In [7]:

query = """
        SELECT coord_ra, coord_dec, g_psfMag, r_psfMag, i_psfMag, ebv
        FROM dp1.Object
        WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec),
              CIRCLE('ICRS', {}, {}, {})) = 1
        AND refExtendedness = 0
        """.format(ra_cen, dec_cen, radius)

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

Job phase is COMPLETED

Fetch the results, store them as a table.

In [9]:

assert job.phase == 'COMPLETED'
df = job.fetch_result().to_table()

2.2. Apply dust correction¶

The DustValues().r_x attribute provides the extinction coefficients for the Rubin filters. These coefficients are derived from the CCM89 extinction law (Cardelli, Clayton, & Mathis 1989), assuming a flat spectral energy distribution and a default total-to-selective extinction ratio of 3.1, representative of the diffuse Milky Way interstellar medium. Apply dust correction in each filter and save the dust-corrected magnitudes as new columns in the df table.

In [10]:

R_band = DustValues().r_x
A_band = {band: R_band[band] * df['ebv'] for band in bands}
for band in bands:
    df[f"{band}_psfMag0"] = df[f"{band}_psfMag"] - A_band[band]

Compute the average dust correction in $i$ and $g-i$.

In [11]:

shift_in_i = np.ma.median(A_band['i'])
shift_in_gi = np.ma.median(A_band['g'] - A_band['i'])

print(f"Average extinction in i: {shift_in_i: .3f} mag")
print(f"Average reddening in g-i: {shift_in_gi: .3f} mag")

Average extinction in i:  0.050 mag
Average reddening in g-i:  0.039 mag

2.3. Dust correction and its impact on the CMD¶

Construct observed and de-reddened ($g−i$, $i$) CMDs, then subtract the former from the latter to visualize the impact of dust correction on the stellar distribution.

In [12]:

grid_size = 200
extent = [-1, 2, 15, 27]

counts_obs, xedges, yedges = np.histogram2d(
    df['g_psfMag'] - df['i_psfMag'], df['i_psfMag'],
    bins=grid_size, range=[[extent[0], extent[1]], [extent[2], extent[3]]]
)

counts_corr, _, _ = np.histogram2d(
    df['g_psfMag0'] - df['i_psfMag0'], df['i_psfMag0'],
    bins=grid_size, range=[[extent[0], extent[1]], [extent[2], extent[3]]]
)

diff = counts_corr.T - counts_obs.T

fig, ax = plt.subplots(figsize=(6, 5))

v_limit = np.percentile(np.abs(diff), 99.5)
im = ax.imshow(diff, extent=extent, aspect='auto', origin='lower',
               cmap='RdBu_r', vmin=-v_limit, vmax=v_limit, interpolation='nearest')

ax.set_xlim(-1, 2)
ax.set_ylim(27, 15)
ax.set_xlabel(r'$g-i$')
ax.set_ylabel(r'$i$')
ax.arrow(1.5, 17, -shift_in_gi*5, -shift_in_i*5, color='black',
         width=0.01, head_width=0.05, length_includes_head=True)
ax.text(1.5, 16.5, 'Correction\nDirection', color='black', ha='center')

cbar = fig.colorbar(im, ax=ax, label='De-reddened CMD - Observed CMD')
plt.show()

Figure 1: Residual 2-dimensional histogram illustrating the shift in the ($g-i$, $i$) CMD of the Fornax dSph after applying Galactic dust corrections. The color scale represents the change in stellar counts per bin ($\Delta$N = N$_{corrected}$ - N$_{observed}$), where blue (negative) regions indicate the original observed positions of stars and red (positive) regions represent their intrinsic, de-reddened positions. The black vector denotes the direction and magnitude of the correction (magnified by x5 for clarity), demonstrating the systematic recovery of the stellar locus toward higher effective temperatures and luminosities.

3. Classification of stellar evolutionary phases¶

CMDs provide a powerful diagnostic for disentangling the complex star formation history of a stellar system. By defining geometric selection criteria in an CMD, it is possible to isolate distinct stellar populations in the stellar system. In this section, we apply these selections to the de-reddened ($g−i$, $i$) CMD of the Fornax dSph to classify individual stars by evolutionary phase, facilitating a comparative study of their spatial and photometric distributions.

3.1. Identify key evolutionary phases¶

Classify key evolutionary phases, including the Main-Sequence Turn-Off (MSTO; low-mass hydrogen-burning stars), Red Giant Branch (RGB; evolved stars with helium cores and hydrogen-burning shells), Horizontal Branch (HB; core helium-burning stars), and Red Clump (RC; metal-rich counterparts to HB stars), within the CMD. The presence of a Blue Plume (BP) suggests recent star formation within the galaxy, although these stars may also be interpreted as a population of blue stragglers.

In [13]:

phase_colors = {
    'MSTO': '#0072B2',
    'BP': '#CC79A7',
    'RGB': '#D55E00',
    'HB': '#009E73',
    'RC': '#E69F00'
}

regions = [
    {'type': 'poly', 'label': 'MSTO',
     'xy': [(-0.1, 24.0), (-0.1, 24.5), (0.5, 24.5), (0.65, 23.5), (0.4, 23.5), (0.25, 24.0)],
     'edgecolor': phase_colors['MSTO']},
    {'type': 'rect', 'label': 'BP',
     'xy': (-0.6, 22.0), 'width': 0.5, 'height': 2.0, 'edgecolor': phase_colors['BP']},
    {'type': 'poly', 'label': 'RGB',
     'xy': [(1.7, 17.5), (1.45, 17.5), (0.8, 20.0), (0.45, 23.0), (0.75, 23.0), (1.0, 20.0)],
     'edgecolor': phase_colors['RGB']},
    {'type': 'rect', 'label': 'RC',
     'xy': (0.65, 20.3), 'width': 0.3, 'height': 0.8, 'edgecolor': phase_colors['RC']},
    {'type': 'rect', 'label': 'HB',
     'xy': (-0.4, 20.5), 'width': 1.0, 'height': 1.3, 'edgecolor': phase_colors['HB']}
]

In [14]:

fig, ax = plt.subplots(figsize=(6, 5))
hb = ax.hexbin(df['g_psfMag0'] - df['i_psfMag0'], df['i_psfMag0'], gridsize=150,
               extent=[-1, 2, 15, 27], cmap='gray_r', norm=LogNorm(vmax=70))
ax.set_xlim(-1, 2)
ax.set_ylim(27, 15)
ax.set_xlabel(r'$g_0 - i_0$')
ax.set_ylabel(r'$i_0$')

legend_elements = []
for reg in regions:
    if reg['type'] == 'rect':
        patch = Rectangle(reg['xy'], reg['width'], reg['height'],
                          fill=False, lw=2, edgecolor=reg['edgecolor'], label=reg['label'])
    elif reg['type'] == 'poly':
        patch = Polygon(reg['xy'], fill=False, lw=2, edgecolor=reg['edgecolor'], label=reg['label'])

    ax.add_patch(patch)
    legend_elements.append(patch)

ax.legend(handles=legend_elements, loc='upper left')

plt.show()

Figure 2: ($g−i$, $i$) CMD of the Fornax dSph. Overlaid geometric regions (polygons) define the selection criteria used to isolate distinct evolutionary phases.

3.2. Tag stars with evolutionary phase¶

Define functions to filter stars using the selection masks established in the previous section.

In [15]:

def get_poly_mask(color, mag, vertices):
    """
    Generate a boolean mask for points falling within a specified 2D polygon.

    Parameters
    ----------
    color : array-like
        1D array of stellar color indices (e.g., g - i).
    mag : array-like
        1D array of stellar magnitudes (e.g., i).
    vertices : list of tuples
        A list of (x, y) coordinates defining the polygon vertices in the 
        (color, mag) plane.

    Returns
    -------
    mask : numpy.ndarray
        A boolean array where True indicates the point is inside the polygon.
    """
    path = Path(vertices)
    points = np.column_stack((color, mag))

    return path.contains_points(points)

def get_rect_mask(color, mag, x_min, y_min, width, height):
    """
    Generate a boolean mask for points within a rectangular region.

    Parameters
    ----------
    color : array-like
        1D array of stellar color indices.
    mag : array-like
        1D array of stellar magnitudes.
    x_min : float
        The minimum color value.
    y_min : float
        The minimum magnitude value.
    width : float
        The extent of the box along the color axis.
    height : float
        The extent of the box along the magnitude axis.

    Returns
    -------
    mask : numpy.ndarray
        A boolean array where True indicates the point is inside the rectangle.
    """
    return (color >= x_min) & (color <= x_min + width) & \
           (mag >= y_min) & (mag <= y_min + height)

Apply the selection masks, assign labels to each population, and add a categorical column (phase_label) to the df table.

In [16]:

df['phase_label'] = 'unclassified'

for reg in regions:
    if reg['type'] == 'rect':
        mask = get_rect_mask(df['g_psfMag0'] - df['i_psfMag0'], df['i_psfMag0'],
                             reg['xy'][0], reg['xy'][1], reg['width'], reg['height'])
    elif reg['type'] == 'poly':
        mask = get_poly_mask(df['g_psfMag0'] - df['i_psfMag0'], df['i_psfMag0'], reg['xy'])

    df['phase_label'][mask] = reg['label']

4. Multi-Diagnostic Visualization¶

Examine the distribution of each labeled population in both color-color and spatial planes. Using these combined diagnostics has great potential to identify population-specific trends, such as spatial gradients, and to ensure the purity of the isolated evolutionary phases.

4.1. CCD diagnostics¶

Visualize how these CMD-selected populations behave in the color-color space ($g-i$ vs. $r-i$).

In [17]:

fig, ax = plt.subplots(figsize=(6, 5))

gi = df['g_psfMag0'] - df['i_psfMag0']
ri = df['r_psfMag0'] - df['i_psfMag0']

msto_mask = df['phase_label'] == 'MSTO'
h = ax.hist2d(gi[msto_mask], ri[msto_mask], bins=200,
              range=((-1, 2), (-1, 2)), norm=LogNorm(), cmap='Blues')
msto_patch = Patch(color=phase_colors['MSTO'], label=f"MSTO (N={msto_mask.sum()})")

for phase, color in phase_colors.items():
    subset = df['phase_label'] == phase
    if phase == 'MSTO':
        continue
    ax.scatter(gi[subset], ri[subset], s=1, color=color,
               label=f"{phase} (N={len(df[subset])})")
ax.set_xlim(-1, 2)
ax.set_ylim(-1, 2)
ax.set_xlabel('$g_0 - i_0$')
ax.set_ylabel('$r_0 - i_0$')

handles, labels = ax.get_legend_handles_labels()
handles.append(msto_patch)
labels.append(msto_patch.get_label())
ax.legend(handles=handles, markerscale=3)

plt.tight_layout()
plt.show()

Figure 3: Distribution of each labeled population in the ($g-i$ vs. $r-i$) CCD. Plot the MSTO as a 2D histogram due to its high stellar density, while representing the remaining populations with scatter plots. Note that while each population occupies a distinct region of the CCD with partial overlap, they collectively define the stellar locus.

4.2. Spatial distribution by populations¶

Explore spatial distribution of distinct populations within the field footprint. Note that the lack of stars in the inner region is an artifact of the LSST Science Pipelines skipping the deblending process in crowded areas.

Define the structural parameters of the Fornax dSph, including the center, position angle, and ellipticity, to construct its elliptical isophotes at the elliptical radii of 0.4, 0.8, and 1.0 degrees from the center.

In [18]:

ra0, dec0 = 39.997, -34.449
pa = 42.0
ellipticity = 0.3
e_radii = [0.4, 0.8, 1.0]

Examine the spatial distribution of each population across the field.

In [19]:

n_phases = len(phase_colors)
n_cols = 3
n_rows = int(np.ceil(n_phases / n_cols))

fig, axes = plt.subplots(n_rows, n_cols, figsize=(10, 3 * n_rows), sharex=True, sharey=True)
axes = axes.flatten()

bg_mask = df['phase_label'] == 'unclassified'
ra_bg, dec_bg = df['coord_ra'][bg_mask], df['coord_dec'][bg_mask]

for i, (phase, color) in enumerate(phase_colors.items()):
    ax = axes[i]
    ax.scatter(ra_bg, dec_bg, s=1, color='lightgrey', alpha=0.3, rasterized=True)
    subset = df[df['phase_label'] == phase]
    ax.scatter(subset['coord_ra'], subset['coord_dec'], s=3, color=color, alpha=0.3)
    ax.set_title(f"Phase: {phase}", fontsize=12, fontweight='bold')
    ax.invert_xaxis()

    for er in e_radii:
        width = 2 * er
        height = 2 * er * (1 - ellipticity)

        ellipse_patch = Ellipse((ra0, dec0), width, height, angle=pa,
                                edgecolor='k', facecolor='none', lw=2, ls='--')
        ax.add_patch(ellipse_patch)

    if i >= n_phases - n_cols:
        ax.set_xlabel('RA [deg]')
    if i % n_cols == 0:
        ax.set_ylabel('Dec [deg]')

for j in range(i + 1, len(axes)):
    axes[j].axis('off')

plt.tight_layout()
plt.show()

Figure 4: Spatial distribution of five distinct stellar populations. The underlying gray density map represents the remaining stars in the Fornax dSph field that fall outside the defined selection masks. Black dashed lines denote elliptical isophotes with semi-major axes of 0.4, 0.8, and 1.0 degrees from the center. It is difficult to identify clear spatial trends as a function of stellar population, primarily due to the limited radial coverage and incomplete PSF photometry in the DP1 dataset for this field.

4.2.1. Known stellar population variation with radius¶

The stellar populations of the Fornax dSph are known to vary with elliptical radius. While the galaxy is dominated by intermediate-age (1-10 Gyr) populations, it also contains ancient (10-14 Gyr) and young ($<$1 Gyr) components. The galaxy displays a prominent radial age gradient, where younger, more metal-rich populations are more centrally concentrated.

The figure below displays the ($B-V$, $V$) CMDs for each elliptical radial bin (Fig. 5 in de Boer et al. (2012)). The DP1-based CMD appears qualitatively similar to those representing the two outermost radii.

In [ ]: