Cutout stacking

Stack continuum map stamps centred on a prior source catalogue.

This example shows how to use stack() to produce a weighted-average stacked map from a set of source positions. No real data is needed — we build a synthetic ContMap with known injected sources and then stack on their positions.

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.coordinates import SkyCoord
from astropy.modeling import models
from astropy.nddata import StdDevUncertainty
from astropy.stats.funcs import gaussian_fwhm_to_sigma
from astropy.table import Table
from astropy.wcs import WCS
from photutils.datasets import make_model_image

from nikamap import ContBeam, ContMap

rng = np.random.default_rng(42)

Build a synthetic map

We create a 128 × 128 pixel map at 2″/pixel with a 12.5″ FWHM beam and inject 30 point sources at uniform random positions.

shape = (256, 256)
pixscale = 2 * u.arcsec
fwhm = 12.5 * u.arcsec
noise_level = 1e-3  # Jy/beam

wcs = WCS(naxis=2)
wcs.wcs.crpix = [shape[1] / 2, shape[0] / 2]
wcs.wcs.cdelt = [-pixscale.to("deg").value, pixscale.to("deg").value]
wcs.wcs.crval = [0.0, 0.0]
wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]

# Random source positions (pixel, integer for exact centering)
nsources = 30
margin = 15  # pixels — keep sources away from edges
x_pix = rng.integers(margin, shape[1] - margin, nsources)
y_pix = rng.integers(margin, shape[0] - margin, nsources)
peak_flux = 5e-3  # Jy/beam

beam_std_pix = (fwhm / pixscale).decompose().value * gaussian_fwhm_to_sigma
source_table = Table(
    {
        "x_mean": x_pix.astype(float),
        "y_mean": y_pix.astype(float),
        "amplitude": np.full(nsources, peak_flux),
        "x_stddev": np.full(nsources, beam_std_pix),
        "y_stddev": np.full(nsources, beam_std_pix),
        "theta": np.zeros(nsources),
    }
)

source_map = (
    make_model_image(shape, models.Gaussian2D(), source_table,
                     model_shape=shape, x_name="x_mean", y_name="y_mean")
)
noise_map = rng.normal(0, noise_level, size=shape)

cm = ContMap(
    source_map + noise_map,
    uncertainty=StdDevUncertainty(np.full(shape, noise_level)),
    wcs=wcs,
    unit=u.Jy / u.beam,
    beam=ContBeam(major=fwhm, pixscale=pixscale),
)

Recover the sky coordinates of the injected sources

world = wcs.pixel_to_world(x_pix, y_pix)
coords = SkyCoord(world)

print(f"Number of stacked positions: {len(coords)}")

Number of stacked positions: 30

Stack the cutouts

stack() extracts a stamp of angular size size around each coordinate and returns an inverse-variance weighted average. The uncertainty of the output map is derived from the combined weights.

stacked = cm.stack(coords, size=40 * u.arcsec)

print(stacked)

INFO: overwriting masked ndarray's current mask with specified mask. [astropy.nddata.nddata]
[[-1.51963863e-04  3.01123592e-04  2.58660627e-04  2.55794530e-04
   4.66767218e-05 -1.98243604e-05 -1.56301076e-04  2.93453373e-04
   1.65482361e-04  3.28288485e-04  7.38877421e-05  2.47519912e-04
   1.10860890e-04 -2.36705098e-05 -2.29370940e-04  5.99733890e-05
   2.69996980e-04 -1.19418136e-05 -2.20768494e-04  2.74542051e-04
   5.19814098e-04]
 [ 8.38169836e-05  3.29675334e-04  3.48521201e-04  1.93159030e-04
   1.12831741e-04  1.04671105e-05  1.53086555e-04  2.44054614e-04
   1.47548428e-04  3.42936615e-04 -3.66731003e-05  4.73143857e-04
  -3.72565956e-05  2.36255645e-04  3.60704828e-05 -2.18046026e-04
   2.40959596e-04  8.79053606e-05  1.51360775e-05  1.37600085e-05
   3.75653268e-04]
 [ 9.33126650e-05  2.71851347e-04  9.40303434e-05  4.69523479e-05
   1.30065418e-04  1.39675308e-04 -1.13069401e-04 -4.62287028e-05
   2.53615043e-05  1.26991505e-05  1.39267920e-04 -2.23712008e-05
   1.00324137e-04 -6.69839509e-05 -9.24852891e-05 -1.09992220e-04
   9.13496323e-06  2.84245149e-04  5.05086494e-06  1.42877385e-04
  -1.60283223e-04]
 [ 3.59744135e-04  1.34196277e-04  5.24724362e-05  2.25753622e-04
   3.36398619e-04  6.19259587e-05 -2.09692417e-04  1.10253312e-04
   3.24891449e-04  3.79612771e-04  1.03582831e-04  6.24726205e-06
  -1.52462402e-04  2.16645317e-04  3.30293301e-04  7.25813818e-06
  -1.67490366e-04  1.88568415e-04  4.68457858e-04  2.84564432e-04
   3.22481499e-04]
 [-2.51862612e-04  2.26748516e-05  6.92990005e-05  1.48595592e-04
   1.80008992e-04 -7.41595408e-05  3.68935682e-05  2.42851298e-04
   3.83714600e-04  2.88782556e-04  5.14342707e-04  3.65165521e-04
   3.74796891e-04  3.53739164e-04 -7.02023592e-05 -1.17866225e-05
   1.11522224e-04  2.65324088e-04  1.53967077e-04 -5.19864922e-05
   2.72951563e-04]
 [ 1.40403503e-04  1.51832166e-04  6.38644499e-05 -7.44453587e-05
   2.58459161e-04  4.82889814e-04  3.43361149e-04  4.20688359e-04
   5.72545288e-04  1.01000305e-03  5.63206523e-04  8.99641816e-04
   4.66935346e-04  3.34058076e-04 -7.49298506e-05 -9.95622993e-05
   6.73094875e-05  4.43632768e-04  2.04159291e-04  3.61180119e-04
   1.94120118e-04]
 [ 4.44374051e-05  4.72072318e-04  5.66234467e-04  2.00066995e-04
   3.41454807e-05  4.89524679e-04  6.63170294e-04  6.31857745e-04
   1.25568455e-03  1.51553978e-03  1.47693668e-03  1.43171060e-03
   1.13705795e-03  6.00117782e-04  2.86827848e-04  1.26966952e-04
   1.06881033e-04  2.20362094e-04  8.47860344e-05  4.19054420e-04
   4.33551934e-04]
 [ 5.01079172e-05  1.76323774e-04 -1.06049265e-04  1.45238216e-04
   3.85048742e-04  5.25728017e-04  8.90178066e-04  1.42342469e-03
   1.81144251e-03  2.46753713e-03  2.74438439e-03  2.58284798e-03
   1.79852716e-03  1.30734511e-03  6.78531001e-04  5.70537039e-04
   2.56964152e-04  3.08693560e-05  3.07132592e-04 -6.99000713e-05
   6.03764463e-04]
 [ 3.64119808e-04  4.33611418e-04  2.74146135e-05  8.10516473e-05
   3.76101962e-04  5.78651942e-04  1.28816210e-03  2.13031224e-03
   2.84369862e-03  3.64944907e-03  4.05363471e-03  3.50935781e-03
   2.31993739e-03  2.22726784e-03  1.51242952e-03  6.22849330e-04
   3.85820434e-04  2.70319723e-04 -1.88178602e-04  2.54457107e-04
   6.78448553e-05]
 [ 3.27365549e-04  3.21709691e-04  2.48144738e-04  1.29144373e-04
   6.40889033e-04  9.60050091e-04  1.50364515e-03  2.33062836e-03
   3.45120158e-03  4.45507636e-03  4.44824069e-03  4.13782168e-03
   3.23233592e-03  2.55885253e-03  1.39212396e-03  8.36765168e-04
   4.20279050e-04  1.60394176e-04  2.12160353e-04  6.41986869e-04
   6.96211306e-04]
 [ 2.18412473e-04  2.55346490e-04  3.67064262e-04  2.62819648e-04
   4.04251242e-04  7.58170791e-04  1.68425305e-03  2.68589251e-03
   3.78695639e-03  4.67772403e-03  4.97888639e-03  4.67516736e-03
   3.85678471e-03  2.53566800e-03  1.69951962e-03  9.93360139e-04
   6.06114375e-04  5.80343445e-04  3.32269760e-04  5.98476522e-04
   3.81179685e-04]
 [ 2.55350359e-04  3.56189360e-04  2.02788314e-04  6.35174044e-04
   2.25394455e-04  8.33188914e-04  1.49099957e-03  2.57133719e-03
   3.44758891e-03  4.19926408e-03  4.56727641e-03  3.95999106e-03
   3.68008045e-03  2.36510748e-03  1.40332040e-03  9.31793604e-04
   5.14954807e-04  5.82508888e-04  5.11497231e-04  2.42150831e-04
   5.90880585e-04]
 [ 6.96992713e-04  6.44022199e-06  1.18339862e-04  5.81225032e-04
   2.39936325e-04  7.20875127e-04  1.06606039e-03  2.08108680e-03
   2.74726781e-03  3.23323305e-03  3.82577416e-03  3.24831654e-03
   2.73487883e-03  2.12804381e-03  1.41960431e-03  7.97435343e-04
   1.81887644e-04  3.18115394e-04  4.01439009e-04  9.09823433e-05
   4.09553231e-04]
 [ 4.01798915e-04  9.62353002e-05  3.68858578e-04  4.45521816e-04
   4.61929973e-04  4.87499524e-04  1.19892711e-03  1.48373686e-03
   2.05008628e-03  2.40502028e-03  2.58260244e-03  2.59894275e-03
   1.95818904e-03  1.44035814e-03  9.83426482e-04  2.65169984e-04
   1.53688102e-04  4.82152727e-05  4.12207221e-04  2.47690379e-04
  -5.41158845e-05]
 [ 4.02214353e-04  3.63384351e-04  2.29313210e-04  2.92319518e-04
  -1.70268482e-04  4.28838558e-04  5.47130442e-04  8.63383184e-04
   1.09195608e-03  1.29507140e-03  1.29391523e-03  1.26252206e-03
   1.00523540e-03  4.91159922e-04  5.02322246e-04  4.07899781e-04
   1.73470354e-04  1.61111502e-04  1.75678250e-04  4.92751461e-05
   1.93837101e-04]
 [-6.52626334e-05  1.18720947e-05  7.93358976e-05  1.38950417e-04
   2.45970186e-04  2.16020993e-04  2.33320978e-04  4.07028570e-04
   7.58975633e-04  1.13875068e-03  6.32750697e-04  6.86959974e-04
   9.76057105e-04  4.20921340e-04  2.74488425e-04  2.50272882e-04
   1.53752238e-04  2.43104561e-04  4.24720921e-04  2.17713594e-04
   3.24289315e-04]
 [-1.01094530e-04  1.21455491e-04 -3.22835708e-06  3.40056538e-05
   7.32050583e-05 -9.25091721e-06  2.92980884e-04  1.13681335e-04
   4.20659613e-04  6.93907793e-04  2.88786184e-04  3.29071811e-04
   3.50449749e-04  2.00010659e-04  3.81698081e-04  3.27964437e-04
  -3.42853815e-04  3.58096900e-04  9.53595090e-05  2.72044082e-04
  -1.43899556e-04]
 [ 6.83455218e-04 -9.80025276e-05 -9.36983886e-05  7.13324067e-05
   8.25356297e-05  1.90871904e-04  3.99881473e-05  5.60284286e-04
   2.30053635e-04  1.53531812e-04  4.60848184e-04  5.61485100e-04
   2.09795074e-04  3.27912493e-04  1.68297357e-06 -1.46183708e-05
  -1.12421669e-04  1.64822971e-04  3.19767498e-04  6.86736752e-04
   5.01636362e-04]
 [ 4.25385374e-04  7.48597015e-05 -2.83781900e-04  3.36268715e-07
  -3.09983752e-04 -2.50055780e-05  7.41819690e-05  1.84400466e-04
   6.45288430e-05 -3.19193761e-05  2.47172616e-04  1.47371606e-04
  -1.11010485e-04 -8.73709324e-05  4.94847491e-04  3.72435521e-04
  -2.74663240e-05 -5.73908101e-05  2.01993179e-04  3.11780568e-04
  -1.14275184e-04]
 [ 1.76061970e-04  1.36737795e-04  2.84853138e-05 -2.33908795e-04
   7.04721484e-05  1.90351714e-04 -1.08376715e-04  5.96965584e-05
   2.80357163e-07  1.94586564e-04  1.50760261e-04  2.33206618e-05
   9.05643251e-05  1.38773937e-04  2.17906198e-04 -2.62699823e-05
  -1.46218736e-04 -1.74526958e-04 -1.40594497e-05  2.60125738e-04
  -1.27919780e-04]
 [-2.47651253e-04  4.58509752e-04  3.31634948e-04  1.66019071e-04
  -1.48785062e-04 -1.99780134e-04 -7.74666736e-05  2.37745788e-04
   1.12161282e-04 -7.56067594e-06  5.09613341e-04  2.20074696e-04
   4.64650048e-04  1.78954775e-05 -2.38509069e-04  1.65637732e-04
   1.49193734e-04 -1.72373637e-04  1.78704015e-04  1.42814218e-04
   3.52591582e-04]] Jy / beam

The peak of the stacked map should recover the injected peak flux

peak = float(np.nanmax(stacked.data))
print(f"Injected peak flux : {peak_flux * 1e3:.1f} mJy/beam")
print(f"Stacked peak flux  : {peak * 1e3:.1f} mJy/beam")

Injected peak flux : 5.0 mJy/beam
Stacked peak flux  : 5.0 mJy/beam

Plot

fig, axes = plt.subplots(1, 2, figsize=(10, 4),
                         subplot_kw={"projection": cm.wcs})
cm.plot(ax=axes[0], cbar=True)
axes[0].set_title("Synthetic map")
axes[0].scatter([c.ra.deg for c in coords],
                [c.dec.deg for c in coords],
                transform=axes[0].get_transform("world"),
                s=10, c="red", label="sources")
axes[0].legend(loc="upper right")

stacked.plot(ax=axes[1], cbar=True)
axes[1].set_title("Stacked stamp")
plt.tight_layout()