312.3. Painting and measuring a tidal feature#
312.3. Painting and measuring a tidal feature¶
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-07-06
Repository: github.com/lsst/tutorial-notebooks
Learning objective: Load a DP1 coadd stamp, interactively paint a mask over a visible tidal feature using the %matplotlib widget backend, and measure the mean and median surface brightness of the painted region in $\mathrm{mag\ arcsec}^{-2}$.
LSST data products: deep_coadd, skyMap
Packages: lsst.daf.butler, lsst.geom, numpy, matplotlib, scipy, astropy
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¶
Tidal features, which include streams, tidal tails, and shells made up from stars stripped from interacting galaxies, are typically too irregular in shape for automated aperture photometry. Visual identification and manual masking is a common first step, and is the basis of dedicated inspection tools such as JAFAR, which uses expert and citizen-science classifications to build catalogs of low-surface-brightness features in wide-field survey images.
The target galaxy used throughout this notebook is taken from the Extended Chandra Deep Field South (ECDFS) imaged by DP1; the same target appears in tutorials 312.1 and 312.2.
The TidalMaskPainter class defined in Section 1.2 provides a simple interactive brush using the %matplotlib widget backend. Clicking and dragging marks a region; the resulting boolean mask is used in Section 5 to measure the mean and median surface brightness in AB magnitudes per square arcsecond ($\mathrm{mag\ arcsec}^{-2}$).
Related tutorials: Tutorial 312.1 introduces image display stretches; 312.2 covers noise-suppression techniques including pixel binning and adaptive smoothing.
1.1. Import packages¶
lsst.daf.butler and lsst.geom access the data and handle sky-coordinate and pixel-geometry operations. astropy.visualization provides ImageNormalize and AsinhStretch, which apply a non-linear stretch that compresses the bright galaxy nucleus while preserving faint diffuse structure. Standard scientific Python packages (numpy, matplotlib) handle array operations and figures.
import numpy as np
import matplotlib.pyplot as plt
from astropy.visualization import AsinhStretch, ImageNormalize
import lsst.geom as geom
from lsst.daf.butler import Butler
1.2. Define parameters and functions¶
Set plotting defaults, then define TARGET, BAND, and STAMP_HALF: the (RA, Dec) sky coordinate of the target galaxy in degrees, the imaging band, and the postage-stamp half-width in pixels (450 pixels corresponds to $\pm 1.5$ arcmin at the LSST pixel scale). PIXEL_SCALE and PIXEL_AREA are derived from the coadd WCS in Section 2.
plt.rcParams.update(
{"figure.dpi": 150, "image.origin": "lower", "axes.titlesize": 10})
TARGET = (53.010, -28.350)
BAND = "i"
STAMP_HALF = 450
build_sky_arr masks detected sources and detector artifacts by setting flagged pixels to NaN, preventing bright compact objects from biasing sky statistics. sky_stats returns display limits and the per-pixel sky noise $\hat{\sigma}$ (the same quantity used throughout Tutorial 312.1).
Sky noise is the pixel-to-pixel flux variation on blank patches of sky, measured after masking sources and artifacts. It sets the fundamental detection floor: any structure fainter than a few $\hat{\sigma}$ per pixel is invisible without additional noise suppression. The noise is estimated robustly as $\sigma_\mathrm{sky} = 1.4826 \times \mathrm{MAD}$, where the factor 1.4826 makes the MAD-based estimate consistent with the standard deviation of a Gaussian distribution.
def build_sky_arr(
masked_image,
arr,
bad_planes=("SAT", "NO_DATA", "BAD", "CR", "BRIGHT_OBJECT"),
source_planes=("DETECTED",),
):
"""Return arr with artifact and source pixels set to NaN.
Parameters
----------
masked_image : lsst.afw.image.MaskedImage
Masked image whose mask plane is used to identify bad pixels.
arr : ndarray, shape (ny, nx)
Pixel array to mask, typically the image or variance plane.
bad_planes : tuple of str, optional
Mask plane names treated as artifacts.
source_planes : tuple of str, optional
Mask plane names treated as detected sources.
Returns
-------
masked_arr : ndarray, shape (ny, nx)
Copy of arr with artifact, source, and non-finite pixels set to NaN.
"""
mask_arr = masked_image.getMask().getArray()
mask_pd = masked_image.getMask().getMaskPlaneDict()
exclude = np.zeros(arr.shape, dtype=bool)
for plane in list(bad_planes) + list(source_planes):
if plane in mask_pd:
exclude |= (mask_arr & (1 << mask_pd[plane])) != 0
return np.where(exclude | ~np.isfinite(arr), np.nan, arr)
def sky_stats(arr, lo_sigma=2.0, hi_sigma=8.0):
"""Return (vmin, vmax, sigma) anchored to the sky background.
Uses the robust median as sky level and 1.4826 x MAD as pixel noise.
Pass a source-masked array (e.g. from build_sky_arr) so that detected
sources do not bias the estimates.
Parameters
----------
arr : ndarray
Pixel array, typically source-masked (NaN at bad/source pixels).
lo_sigma : float, optional
Lower stretch limit in units of sky noise below the sky median.
hi_sigma : float, optional
Upper stretch limit in units of sky noise above the sky median.
Returns
-------
vmin : float
Lower display limit (sky - lo_sigma * noise).
vmax : float
Upper display limit (sky + hi_sigma * noise).
sigma : float
Per-pixel noise estimate (1.4826 * MAD).
"""
flat = arr[np.isfinite(arr)]
sky = np.median(flat)
sigma = 1.4826 * np.median(np.abs(flat - sky))
return sky - lo_sigma * sigma, sky + hi_sigma * sigma, sigma
sky_norm wraps sky_stats into an ImageNormalize object for direct use with imshow.
def sky_norm(arr, lo_sigma=2.0, hi_sigma=8.0):
"""Return an ImageNormalize anchored to the sky background.
Wraps sky_stats. Pass a source-masked array so that detected objects do
not bias the sky median or noise estimate.
Parameters
----------
arr : ndarray
Pixel array, typically source-masked (NaN at bad/source pixels).
lo_sigma : float, optional
Lower stretch limit in units of sky noise below the sky median.
hi_sigma : float, optional
Upper stretch limit in units of sky noise above the sky median.
Returns
-------
norm : astropy.visualization.ImageNormalize
Normalization with asinh stretch and sky-anchored limits.
"""
vmin, vmax, _ = sky_stats(arr, lo_sigma=lo_sigma, hi_sigma=hi_sigma)
return ImageNormalize(vmin=vmin, vmax=vmax, stretch=AsinhStretch())
nJy_per_pix_to_mu converts pixel flux in nJy to AB surface brightness in mag arcsec$^{-2}$, using PIXEL_AREA.
PIXEL_SCALE and PIXEL_AREA are not passed to the function as arguments: they are assigned as global variables in Section 2, immediately after the coadd is loaded. Python looks up a name used inside a function body in the global namespace at the time the function is called, not when it is defined, so nJy_per_pix_to_mu can be defined here, before PIXEL_AREA exists, and still work correctly as long as the notebook is run top to bottom before the function is used.
def nJy_per_pix_to_mu(flux_nJy):
"""Convert nJy pixel**-1 to AB mag arcsec**-2 for the LSST pixel scale.
Parameters
----------
flux_nJy : float or array_like
Pixel flux in nJy pixel**-1.
Returns
-------
mu : float or array_like
Surface brightness in AB mag arcsec**-2.
"""
return -2.5 * np.log10(np.abs(flux_nJy) * 1e-9 / PIXEL_AREA / 3631.0)
TidalMaskPainter creates an interactive painting canvas using the %matplotlib widget backend. Clicking and dragging on the displayed image applies a circular brush at the cursor position, blending an orange overlay onto the canvas in real time. The boolean mask recording which pixels have been painted is available via the .mask attribute.
The class connects three matplotlib canvas events: button_press_event, motion_notify_event, and button_release_event. The brush footprint is a filled disc of radius brush pixels; the _paint method updates both the boolean mask and the displayed canvas.
_ORANGE_RGB = np.array([254.0, 97.0, 0.0])
class TidalMaskPainter:
"""Interactive brush painter for marking tidal features.
Click and drag on the displayed image to paint; retrieve the boolean
mask via the `.mask` attribute.
Parameters
----------
arr : ndarray, shape (ny, nx)
Pixel array to display as the painting canvas.
norm : astropy.visualization.ImageNormalize
Normalization applied to arr for display.
brush : int, optional
Brush radius in pixels.
Attributes
----------
mask : ndarray of bool, shape (ny, nx)
Boolean mask; True where the user has painted.
"""
def __init__(self, arr, norm, brush=8):
ny, nx = arr.shape
self.mask = np.zeros((ny, nx), dtype=bool)
self._brush = brush
self._active = False
self._last_drawn = (-1, -1)
arr_disp = np.where(np.isfinite(arr), arr, 0.0)
normed = np.clip(np.asarray(norm(arr_disp)), 0.0, 1.0)
self._base = (plt.cm.Greys_r(normed) * 255).astype(np.uint8)
self._canvas = self._base.copy()
self.fig, self.ax = plt.subplots(figsize=(5, 5))
self._im = self.ax.imshow(
self._canvas,
origin="lower",
extent=[-0.5, nx - 0.5, -0.5, ny - 0.5],
)
self.ax.set_title(f"Click and drag to paint -- brush = {brush} px")
self.fig.canvas.mpl_connect("button_press_event", self._press)
self.fig.canvas.mpl_connect("motion_notify_event", self._motion)
self.fig.canvas.mpl_connect("button_release_event", self._release)
def _paint(self, bx, by, force=False):
"""Paint a circular brush footprint at pixel (bx, by)."""
r = self._brush
ny, nx = self.mask.shape
y0, y1 = max(0, by - r), min(ny, by + r + 1)
x0, x1 = max(0, bx - r), min(nx, bx + r + 1)
yy, xx = np.ogrid[y0:y1, x0:x1]
circle = (yy - by) ** 2 + (xx - bx) ** 2 <= r**2
self.mask[y0:y1, x0:x1] |= circle
base_s = self._base[y0:y1, x0:x1, : 3].astype(np.float32)
canvas_s = self._canvas[y0:y1, x0:x1]
blended = np.clip(0.55 * base_s + 0.45 * _ORANGE_RGB,
0, 255).astype(np.uint8)
canvas_s[circle, : 3] = blended[circle]
ldx = bx - self._last_drawn[0]
ldy = by - self._last_drawn[1]
if force or ldx * ldx + ldy * ldy >= r * r:
self._im.set_data(self._canvas)
self.fig.canvas.draw_idle()
self._last_drawn = (bx, by)
def _press(self, ev):
"""Activate painting and draw at the click position."""
if ev.inaxes == self.ax and ev.button == 1:
self._active = True
self._paint(int(round(ev.xdata)), int(round(ev.ydata)), force=True)
def _motion(self, ev):
"""Continue painting as the cursor moves."""
if self._active and ev.inaxes == self.ax:
self._paint(int(round(ev.xdata)), int(round(ev.ydata)))
def _release(self, ev):
"""Deactivate painting and finalise the canvas."""
self._active = False
self._im.set_data(self._canvas)
self.fig.canvas.draw_idle()
2. Load a coadd postage stamp¶
Instantiate the Butler with the DP1 repository and load a $\pm 1.5$ arcmin stamp centered on the target at RA = 53.010, Dec = -28.350.
Resolve the target sky coordinate to its tract and patch in the lsst_cells_v1 tessellation. The tract WCS and outer bounding box are retrieved here for use in the next step.
TARGET_RA, TARGET_DEC = TARGET
butler = Butler("dp1", collections="LSSTComCam/DP1")
skymap_name = "lsst_cells_v1"
skymap = butler.get("skyMap", skymap=skymap_name)
sky_point = geom.SpherePoint(TARGET_RA * geom.degrees, TARGET_DEC * geom.degrees)
tract_info = skymap.findTract(sky_point)
patch_info = tract_info.findPatch(sky_point)
tract_id = tract_info.getId()
patch_id = patch_info.getSequentialIndex()
patch_wcs = tract_info.getWcs()
patch_bbox = patch_info.getOuterBBox()
data_id = {"band": BAND, "tract": tract_id, "patch": patch_id}
Convert the sky coordinate to a pixel position within the patch WCS, then construct a square bounding box of STAMP_HALF pixels on each side centred on the host galaxy. The box is clipped to the patch boundary to avoid requesting pixels outside the coadd.
host_pix = patch_wcs.skyToPixel(sky_point)
stamp_bbox = geom.Box2I(
geom.Point2I(int(host_pix.getX()) - STAMP_HALF,
int(host_pix.getY()) - STAMP_HALF),
geom.Extent2I(2 * STAMP_HALF, 2 * STAMP_HALF),
)
stamp_bbox.clip(patch_bbox)
Pass the bounding box to the Butler as a parameter to retrieve only the stamp cutout from the deep_coadd dataset. The pixel scale and pixel area are derived from the stamp WCS.
stamp = butler.get("deep_coadd", dataId=data_id,
parameters={"bbox": stamp_bbox})
stamp_arr = stamp.image.array.copy()
host_x = int(host_pix.getX()) - stamp_bbox.getMinX()
host_y = int(host_pix.getY()) - stamp_bbox.getMinY()
PIXEL_SCALE = stamp.getWcs().getPixelScale().asArcseconds()
PIXEL_AREA = PIXEL_SCALE**2
arcmin = 2 * STAMP_HALF * PIXEL_SCALE / 60
print(
f"Loaded {BAND}-band stamp: {stamp_arr.shape} "
f"({arcmin: .1f}' x {arcmin: .1f}')"
)
print(f"Host galaxy at stamp pixel: ({host_x}, {host_y})")
print(f"Pixel scale: {PIXEL_SCALE: .4f} arcsec/pixel")
Loaded i-band stamp: (847, 900) ( 3.0' x 3.0') Host galaxy at stamp pixel: (450, 450) Pixel scale: 0.2000 arcsec/pixel
3. Inspect the image and build the sky mask¶
Mask detected sources and detector artifacts using build_sky_arr, leaving only empty-sky pixels. The sky noise $\hat{\sigma} = 1.4826 \times \mathrm{MAD}$ of these pixels sets the detection floor for the surface brightness measurement in Section 5. The stamp is displayed at native resolution with an arcsinh stretch from $-2$ to $+40\,\hat{\sigma}$; the orange cross marks the target position.
stamp_sky = build_sky_arr(stamp.getMaskedImage(), stamp_arr)
sky_val = np.nanmedian(stamp_sky[np.isfinite(stamp_sky)])
_, _, sky_sigma = sky_stats(stamp_sky, lo_sigma=2, hi_sigma=40)
norm = sky_norm(stamp_sky, lo_sigma=2, hi_sigma=40)
mu_sky = nJy_per_pix_to_mu(sky_val)
mu_sigma = nJy_per_pix_to_mu(sky_sigma)
print(f"Sky level: {sky_val: .1f} nJy/pix = {mu_sky: .1f} mag/arcsec^2")
print(f"Per-pixel 1-sigma noise: {sky_sigma: .1f} nJy/pix = {mu_sigma: .1f} mag/arcsec^2")
fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(stamp_arr, norm=norm, cmap="Greys_r", origin="lower")
ax.plot(host_x, host_y, "+", color="#FE6100",
ms=14, mew=2.5, label="target galaxy")
ax.legend(fontsize=9)
Sky level: -0.1 nJy/pix = 30.0 mag/arcsec^2 Per-pixel 1-sigma noise: 2.8 nJy/pix = 26.8 mag/arcsec^2
<matplotlib.legend.Legend at 0x7f774feccd70>
Figure 1: i-band stamp centered on the target galaxy (orange cross), displayed with an arcsinh stretch from $-2$ to $+40\,\sigma$. The stretch is tuned to reveal faint diffuse emission around the central source.
4. Paint the tidal feature¶
Run the cell below to open the interactive painting figure. Click and drag to paint over a candidate tidal feature; the painted region is shown in orange. The brush radius is set to 15 pixels ($3''$). When finished painting, run Section 5 to compute the photometry.
For the best result, avoid the bright host galaxy core and any obvious foreground stars or background galaxies — only detector artifacts are masked, so all astronomical sources contribute to the measurement. The median surface brightness is more robust than the mean if the painted region contains residual compact sources.
To start over, re-run this cell to reset the painter, then run Section 5 again.
%matplotlib widget
painter = TidalMaskPainter(stamp_arr, norm, brush=15)
5. Measure surface brightness¶
Run the cell below after painting. stamp_artifacts retains detected sources but excludes detector artifacts (saturated pixels, bad columns, cosmic rays), so the painted tidal feature remains measurable. The painted mask is intersected with this array to exclude bad pixels, and the mean and median surface brightness are computed in $\mathrm{mag\ arcsec}^{-2}$. The quoted uncertainty on the mean is $(2.5/\ln(10)) \times \hat{\sigma}/\sqrt{N}$, a statistical lower bound that does not account for pixel-to-pixel correlations in the coadd.
Note: A more accurate uncertainty can be obtained by placing many randomly drawn apertures of the same pixel count on blank-sky pixels (stamp_sky) and measuring the standard deviation of their mean fluxes. This approach accounts for pixel-to-pixel correlations automatically.
stamp_artifacts = build_sky_arr(
stamp.getMaskedImage(), stamp_arr, source_planes=())
plt.close("all")
native_mask = painter.mask
pixels = mean_flux = med_flux = None
mu_mean = mu_med = depth_1s = depth_3s = None
if native_mask.sum() == 0:
print("No pixels painted. Run the painting cell and draw over a feature first.")
else:
valid = native_mask & np.isfinite(stamp_artifacts)
n_total = int(native_mask.sum())
n_valid = int(valid.sum())
print(
f"Painted region: {n_total} pixels ({n_total * PIXEL_AREA: .0f} arcsec2)")
print(
f"After masking: {n_valid} usable pixels "
f"({100 * n_valid / max(n_total, 1): .0f}%)"
)
print()
if n_valid < 10:
print("Too few usable pixels -- try a larger aperture.")
else:
pixels = stamp_arr[valid]
mean_flux = float(np.nanmean(pixels))
med_flux = float(np.nanmedian(pixels))
depth_1s = nJy_per_pix_to_mu(sky_sigma)
depth_3s = nJy_per_pix_to_mu(3 * sky_sigma)
if mean_flux > 0:
mu_mean = nJy_per_pix_to_mu(mean_flux)
sigma_mu = (2.5 / np.log(10)) * sky_sigma / (np.sqrt(n_valid) * mean_flux)
print(
f"Mean SB: {mu_mean: .2f} +/- {sigma_mu: .2f} mag/arcsec2"
f" (stat. lower bound)"
)
if med_flux > 0:
mu_med = nJy_per_pix_to_mu(med_flux)
print(f"Median SB: {mu_med: .2f} mag/arcsec2")
print()
print(f"1-sigma depth: {depth_1s: .2f} mag/arcsec2")
print(f"3-sigma depth: {depth_3s: .2f} mag/arcsec2")
print(f"SNR (median): {med_flux / sky_sigma: .1f}")
Painted region: 7485 pixels ( 299 arcsec2) After masking: 7413 usable pixels ( 99%) Mean SB: 24.77 +/- 0.00 mag/arcsec2 (stat. lower bound) Median SB: 24.95 mag/arcsec2 1-sigma depth: 26.77 mag/arcsec2 3-sigma depth: 25.58 mag/arcsec2 SNR (median): 5.3
Plot the painted region overlaid on the stamp and the surface brightness distribution of above-sky pixels in the painted region.
if pixels is None:
print("Run Section 5 first with a painted region of at least 10 usable pixels.")
else:
ny_m, nx_m = native_mask.shape
fig, axes = plt.subplots(1, 2, figsize=(8, 4))
axes[0].imshow(stamp_arr, norm=norm, cmap="Greys_r", origin="lower")
rgba = np.zeros((ny_m, nx_m, 4), dtype=float)
rgba[native_mask] = [0.996, 0.380, 0.0, 0.5]
axes[0].imshow(
rgba, origin="lower", extent=[-0.5, nx_m - 0.5, -0.5, ny_m - 0.5]
)
axes[0].set_title("Painted region")
pos_pix = pixels[pixels > 0]
mu_pos = nJy_per_pix_to_mu(pos_pix)
axes[1].hist(
mu_pos, bins=60, color="#FE6100", edgecolor="k", linewidth=0.4, alpha=0.85
)
if mean_flux > 0:
axes[1].axvline(
mu_mean, color="k", ls="--", lw=2, label=f"mean = {mu_mean: .1f}"
)
if med_flux > 0:
axes[1].axvline(
mu_med, color="k", ls="-", lw=2, label=f"median = {mu_med: .1f}"
)
axes[1].axvline(
depth_1s, color="0.5", ls=":", lw=1.5, label=f"1-sigma = {depth_1s: .1f}"
)
axes[1].set_xlim(right=depth_1s + 2)
axes[1].set_xlabel("Surface brightness ($\\mathrm{mag\\ arcsec}^{-2}$)")
axes[1].set_ylabel("Count")
axes[1].set_title("Surface brightness distribution")
axes[1].legend(fontsize=8)
Figure 2: Left -- painted region (orange overlay) on the native i-band image. Right -- surface brightness distribution ($\mathrm{mag\ arcsec}^{-2}$) of above-sky pixels in the painted region; vertical lines mark the mean (black dashed), median (black solid), and 1-sigma depth (grey dotted).
6. Summary¶
- Tidal features and other irregular structures may be best measured by defining a freehand region directly on the image; the resulting pixel sample gives direct access to the flux distribution of the selected area without assumptions about shape or symmetry.
- The median is the preferred surface brightness estimator for diffuse emission because it is robust to unmasked compact sources within the aperture, unlike the mean which is pulled upward by any residual point sources.
- Comparing the measured surface brightness against the per-pixel sky noise floor shows directly whether a candidate feature is a genuine detection or consistent with sky noise..