311.1. Milky Way Dust Extinction Correction

311.1. Milky Way Dust Extinction Correction¶

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

Learning objective: To learn how to correct for Milky Way foreground dust extinction.

LSST data products: object

Packages: lsst.rsp.RSPDiscovery

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¶

Interstellar dust removes flux by absorbing photons and scattering them out of the observer's line of sight, a process known as "extinction". The amount of extinction depends on wavelength: dust affects shorter (bluer) wavelengths more strongly than longer (redder) ones, leading to a "reddening" of the light. The wavelength dependence is described by an extinction curve, from which the size distribution and composition of interstellar dust grains can be inferred. Such curves are used to correct observations for dust effects. Without these corrections, the derived properties of celestial objects, such as intrinsic luminosity and color, will be inaccurate.

This notebook demonstrates how to correct for the effects of Galactic foreground dust along the line of sight using the ebv column in the Object table, which provides $E(B-V)$ values at given RA/Dec coordinates per Schlegel, Finkbeiner & Davis (1998) (SFD98 hereafter). Interstellar reddening is quantified by the color excess, defined as the difference between an object’s observed and intrinsic (i.e., dust-free) color indices: $E(B-V) = (B-V)_{observed} - (B-V)_{intrinsic}$. E(B-V) measures how much interstellar dust reddens incoming flux and serves as a direct proxy for the amount of the line-of-sight extinction. The total-to-selective extinction ratio at wavelength $\lambda$ is defined as $R_\lambda = \frac{A_\lambda}{E(B-V)}$, which characterizes the degree to which dust dims light at that wavelength.

Related tutorials: See the 200-level and 300-level DP2 tutorials for guidance on the Object table.

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 module for accessing the Table Access Protocol (TAP) service, and the DustValues class from the rubin_sim package (rubin_sim).

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from astropy import units as u
from astropy.coordinates import SkyCoord
import os

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

1.2. Define parameters and functions¶

Instantiate RSPDiscovery with the DP2 release, create an instance of the TAP service.

In [2]:

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

Define a 1.0-degree radius around the coordinates (RA, Dec) = (225 deg, -40 deg).

In [3]:

ra_cen = 225
dec_cen = -40
radius = 1.0

List bands available in the target region.

In [4]:

bands = "ugrizy"

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

In [5]:

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

2. Explore reddening property¶

The Object table contains forced photometric measurements on the deep coadded images at the locations with signal-to-noise ratio > 5 in a deep_coadd of any filter, as well as corresponding $E(B-V)$ values from the SFD98 dust map (ebv column).

2.1. Query Object table¶

Query the Object table for coordinates, photometric measurements, extendedness parameter, and $E(B-V)$ values.

In [6]:

query = """
        SELECT coord_ra, coord_dec,
            u_psfMag, g_psfMag, r_psfMag, i_psfMag, z_psfMag, y_psfMag,
            u_cModelMag, g_cModelMag, r_cModelMag, i_cModelMag, z_cModelMag, y_cModelMag,
            refExtendedness, ebv
        FROM dp2.Object
        WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec),
            CIRCLE('ICRS', {}, {}, {})) = 1
        """.format(ra_cen, dec_cen, radius)

In [7]:

job = rsp_tap.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 and store them as a table.

In [8]:

table = job.fetch_result().to_table()
print(f"The query returned {len(table)} objects.")

The query returned 1101019 objects.

Option to display the table of results.

In [9]:

# table

2.2. Display maps¶

Display the object count map, median $E(B-V)$ map, and 1-D $E(B-V)$ histogram for the target region.

In [10]:

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 4))

hb = ax1.hexbin(table['coord_ra'], table['coord_dec'],
                gridsize=100, bins='log')
ax1.set_xlabel('RA [deg]')
ax1.set_ylabel('Dec [deg]')
ax1.invert_xaxis()
cb = fig.colorbar(hb, ax=ax1, label='Counts of all objects per bin')
ax1.annotate('', xy=(0.15, 0.95), xytext=(0.15, 0.85),
             arrowprops=dict(facecolor='m', width=2, headwidth=5),
             xycoords='axes fraction')
ax1.text(0.15, 0.95, 'N', transform=ax1.transAxes,
         ha='center', va='bottom', fontsize=10)
ax1.annotate('', xy=(0.05, 0.85), xytext=(0.15, 0.85),
             arrowprops=dict(facecolor='m', width=2, headwidth=5),
             xycoords='axes fraction')
ax1.text(0.04, 0.85, 'E', transform=ax1.transAxes,
         ha='right', va='center', fontsize=10)

hb2 = ax2.hexbin(table['coord_ra'], table['coord_dec'], table['ebv'],
                 reduce_C_function=np.ma.median, gridsize=100)
ax2.set_xlabel('RA [deg]')
ax2.set_ylabel('Dec [deg]')
ax2.invert_xaxis()
cb = fig.colorbar(hb2, label="Median E(B-V)")

ax3.hist(table['ebv'], bins=100)
ax3.set_xlabel('E(B-V)')
ax3.set_ylabel('Number of objects')

plt.tight_layout()
plt.show()

Figure 1: The left panel shows the spatial distribution of all objects in the field, revealing clumps and filaments. The middle panel presents the median reddening per spatial bin, which correlates spatially with the stellar density map. The right panel shows the distribution of E(B-V) values in the field. This region exhibits significantly high dust reddening, with a median of $E(B-V)$ = 0.09.

3. Correct the effects of dust¶

In rubin_sim.phot_utils, DustValues().r_x provides the $R_{\lambda}$ values for the Rubin filters. These values are calculated using the CCM89 extinction curve (Cardelli, Clayton, and Mathis 1989), assuming a flat spectral energy distribution. The default $R_{V}$ is 3.1 (for the Milky Way diffuse ISM), but it is possible to specify a different value with the keyword argument r_v.

Print out $R_{\lambda}$ values for $ugrizy$.

In [11]:

R_band = DustValues().r_x
print(R_band)

{'u': np.float64(4.757217815396922), 'g': np.float64(3.6605664439892616), 'r': np.float64(2.70136780871597), 'i': np.float64(2.0536599130965882), 'z': np.float64(1.5900964472616756), 'y': np.float64(1.3077049588254708)}

3.1. Compute extinction¶

$E(B−V)$ sets the scale of dust along the line of sight, while $R_\lambda$, combined with an extinction curve, specifies the wavelength dependance. Together, extinction in a given band can be computed as: $A_\lambda = R_\lambda\times\,E(B-V)$. Apply this relation to compute extinction in each band.

In [12]:

A_band = {band: R_band[band] * table['ebv'] for band in bands}

3.2. Apply extinction corrections¶

Apply extinction corrections to the observed magnitudes in each band ($m_{\lambda,0} = m_\lambda - A_\lambda$), and add the corrected magnitudes back into the table.

In [13]:

for band in bands:
    table[f"{band}_psfMag0"] = table[f"{band}_psfMag"] - A_band[band]
    table[f"{band}_cModelMag0"] = table[f"{band}_cModelMag"] - A_band[band]

3.3. Compare magnitudes¶

Compare magnitudes before and after dust correction both for stars (i.e., refExtendedness = 0) and galaxies (i.e., refExtendedness = 1).

In [14]:

is_star = table['refExtendedness'] == 0
is_galaxy = table['refExtendedness'] == 1

bins = np.linspace(0, 1.0, 100)

fig, axes = plt.subplots(6, 2, figsize=(10, 18), sharex=True)
plt.subplots_adjust(wspace=0.25, hspace=0.35)

axes[0, 0].set_title("PSF (stars)")
axes[0, 1].set_title("cModel (galaxies)")

for i, b in enumerate(bands):
    m = table[f"{b}_psfMag"]
    m0 = table[f"{b}_psfMag0"]
    ok = is_star & np.isfinite(m) & np.isfinite(m0)
    dmag = np.asarray(m[ok]) - np.asarray(m0[ok])
    ax = axes[i, 0]
    ax.hist(dmag, bins=bins)
    med = np.median(dmag) if dmag.size else np.nan
    if np.isfinite(med):
        ax.axvline(med, c="r", linestyle="--")
        ax.text(0.98, 0.95, f"median = {med: .3f} mag",
                transform=ax.transAxes, ha="right", va="top")
    ymin, ymax = ax.get_ylim()
    ax.set_ylim(ymin, ymax * 1.2)
    ax.set_ylabel(f"{b}-band count")

    m = table[f"{b}_cModelMag"]
    m0 = table[f"{b}_cModelMag0"]
    ok = is_galaxy & np.isfinite(m) & np.isfinite(m0)
    dmag = np.asarray(m[ok]) - np.asarray(m0[ok])
    ax = axes[i, 1]
    ax.hist(dmag, bins=bins)
    med = np.median(dmag) if dmag.size else np.nan
    if np.isfinite(med):
        ax.axvline(med, c="r", linestyle="--")
        ax.text(0.98, 0.95, f"median = {med: .3f} mag",
                transform=ax.transAxes, ha="right", va="top")
    ymin, ymax = ax.get_ylim()
    ax.set_ylim(ymin, ymax * 1.2)

axes[-1, 0].set_xlabel(r"m - m$_{0}$ [mag]")
axes[-1, 1].set_xlabel(r"m - m$_{0}$ [mag]")

plt.tight_layout()
plt.show()

Figure 2: Histograms of magnitude differences between observed and dust-corrected values (i.e., extinction) for stars (left panels, PSF magnitudes) and galaxies (right panels, cModel magnitudes) in the $ugrizy$ bands. The median of each distribution is marked with a red dashed line, and its value is printed inside each panel. The extinction values are largest in the bluest $u$ band and decrease toward longer wavelengths, with the reddest $y$ band showing the smallest dust effect, independent of the photometric measurement (PSF vs. cModel).

4. Caveat for using the SFD98 map¶

The SFD98 map is known to be uncertain and unreliable along sightlines with Galactic latitude |b| < 5$^{\circ}$ due to contamination in the 100 $\mu$m dust emission map.

The Trifid Nebula (M20) lies within this zone. It consists of open star cluster, emission nebular, and reflection nebular with hot, structured dust, conditions under which the SFD98 map is known to be biased. In addition, the map provides the full dust column to infinity, whereas much of the dust toward the Trifid Nebula is at only ~1.5 kpc. In short, applying the SFD98 map to foreground stars or nearby H II regions leads to over-correction.

4.1. Retrieve Nebula data¶

Define the coordinates of the Trifid Nebula and confirm that its Galactic latitude lies within 5$^{\circ}$ of the Galactic planewithin (|b| < 5$^{\circ}$).

In [15]:

ra_cen_neb = 270.675451
dec_cen_neb = -22.971902

c = SkyCoord(ra=ra_cen_neb*u.degree, dec=dec_cen_neb*u.degree,
             frame='icrs')
print(c.galactic)

<SkyCoord (Galactic): (l, b) in deg
    (7.08589368, -0.28780377)>

Retrieve only the columns required to construct a ($g-r,g$) color–magnitude diagram for stars in the Trifid Nebula field.

In [16]:

query = """
        SELECT objectId, coord_ra, coord_dec, g_psfMag, r_psfMag, ebv
        FROM dp2.Object
        WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec),
            CIRCLE('ICRS', {}, {}, {})) = 1
        """.format(ra_cen_neb, dec_cen_neb, radius)

In [17]:

job = rsp_tap.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 and save it as an astropy table.

In [18]:

table_neb = job.fetch_result().to_table()

4.2. Apply extinction correction¶

Repeat the procedures from Sections 3.1 and 3.2 to apply dust corrections to the Trifid Nebula field photometry, limiting the correction to the retrieved $g$ and $r$ bands.

In [19]:

bands_neb = "gr"
A_band = {band: R_band[band] * table_neb['ebv'] for band in bands_neb}

for band in bands_neb:
    table_neb[f"{band}_psfMag0"] = (
        table_neb[f"{band}_psfMag"] - A_band[band]
    )

Compare the ($g-r, g$) color-magnitude diagram constructed from the observed photometry with that built from the dust-corrected photometry.

In [20]:

g_ok = np.isfinite(table_neb["g_psfMag"])
r_ok = np.isfinite(table_neb["r_psfMag"])
ok = g_ok & r_ok

g = table_neb["g_psfMag"][ok]
r = table_neb["r_psfMag"][ok]
g0 = table_neb["g_psfMag0"][ok]
r0 = table_neb["r_psfMag0"][ok]

fig, axes = plt.subplots(1, 2, figsize=(8, 4))

h1 = axes[0].hist2d(
    g-r, g,
    range=[(-2, 6), (15, 35)],
    bins=200,
    norm=LogNorm()
)
axes[0].set_xlabel("g − r")
axes[0].set_ylabel("g")
axes[0].set_title("Observed")
axes[0].invert_yaxis()

h2 = axes[1].hist2d(
    g0-r0, g0,
    range=[(-40, 6), (-100, 15)],
    bins=200,
    norm=LogNorm(vmax=h1[0].max())
)
axes[1].set_xlabel("g − r")
axes[1].set_title("Dust-corrected")
axes[1].invert_yaxis()

plt.tight_layout()
plt.show()

Figure 3: $(g-r,\, g)$ color–magnitude diagrams of the Trifid Nebula region using observed magnitudes (left) and dust-corrected magnitudes (right). The unrealistic color and magnitude ranges after dust correction demonstrate that the SFD98 map does not provide reliable E(B-V) values for this low Galactic latitude field.

5. Exercise for the learner¶

Apply dust correction with a different $R_{V}$ than the default 3.1 using the keyword argument r_v in DustValues(). Compare these photometry with those corrected with $R_{V}$ = 3.1 in Section 3.

In [ ]: