Note

This page is generated from the Jupyter notebook notebooks/05_noise_and_saturation.ipynb in the repository. You can download it and run it interactively.

Noise and saturation

wcc_sim applies the wcc_etc noise model exactly: sky + dark, per-frame full-well / ADC saturation, Poisson noise, read noise, and digitization to ADU. This notebook walks through it with a controlled star field: a brightness ladder from G = 6 to 16 at +2 waves defocus.

Requirements: a kernel with wcc_etc and wcc_sim installed. Fully offline (synthetic catalog).

import numpy as np
import matplotlib.pyplot as plt

import wcc_etc
import wcc_sim

from wcc_sim.plotting import use_style, TEAL, RED, AMBER
use_style()

print(f'wcc_sim {wcc_sim.__version__} | wcc_etc {wcc_etc.__version__}')

wcc_sim 0.1.0 | wcc_etc 0.1.0

1. A brightness ladder

Six stars, G = 6 to 16, placed along a column at x = 400. (The left half of the frame is not star-free: the G = 6 star’s scattered-light halo covers the whole array, see the scattered-light page of the user guide and field.plot_psf_profile() in notebook 10, so the background statistics below come from a source-free run.) Coordinates are generated through the same TAN WCS the simulator will build.

from astropy.table import Table
from wcc_sim import simulate_field
from wcc_sim.detectors import get_geometry
from wcc_sim.wcsutil import build_wcs

RA0, DEC0 = 150.1, 2.2
SHAPE = (512, 512)
EXPTIME = 90.0
GMAGS = np.array([6.0, 8.0, 10.0, 12.0, 14.0, 16.0])

g = get_geometry('zwo:r')
w = build_wcs(RA0, DEC0, g.plate_scale_mas, 0.0, SHAPE)
ys = np.linspace(60, 460, len(GMAGS))
ras, decs = w.pixel_to_world_values(np.full_like(ys, 400.0), ys)
cat = Table({
    'source_id': np.arange(len(GMAGS)), 'ra': ras, 'dec': decs,
    'phot_g_mean_mag': GMAGS,
    'phot_bp_mean_mag': np.full(len(GMAGS), np.nan),  # -> G2V
    'phot_rp_mean_mag': np.full(len(GMAGS), np.nan),
})

field = simulate_field(RA0, DEC0, sensorfilter='zwo:r', focus=2, exptime=EXPTIME,
                       seed=42, catalog=cat, shape=SHAPE)
p = field.params
print(f"sky {p['sky_e_s']:.3f} + dark {p['dark_e_s']:.4f} e-/s/px | "
      f"read noise {p['read_noise']:.2f} e- | gain {p['gain']:.3f} e-/ADU | "
      f"well depth {p['well_depth']} e-")
print(f'saturated pixels: {field.saturation_mask.sum()}')

/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: divide by zero encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: overflow encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: invalid value encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: divide by zero encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: overflow encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: invalid value encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total

sky 0.014 + dark 0.0031 e-/s/px | read noise 3.04 e- | gain 0.261 e-/ADU | well depth 16274.929434794496 e-
saturated pixels: 20524

ax = field.plot(satmask=True, percent=99.5)
for gmag, y in zip(GMAGS, ys):
    ax.annotate(f'G = {gmag:.0f}', (400, y), xytext=(310, y - 4), color=RED)
ax.set_title('+2 waves defocus, 90 s, saturated pixels outlined in red')
plt.show()

Note the shape of the saturated regions: the defocused PSF saturates in its bright ring first. For stars near the threshold the central pixel stays below full well. This is why the per-star saturated flag checks a 32 px window (wcc_sim.render.star_saturated), not just the central pixel.

sat_flags = np.asarray(field.catalog['saturated'])
for gmag, flag in zip(GMAGS, sat_flags):
    print(f'G = {gmag:4.1f}  saturated = {bool(flag)}')
assert sat_flags[GMAGS <= 8].all() and not sat_flags[GMAGS >= 14].any()

G =  6.0  saturated = True
G =  8.0  saturated = True
G = 10.0  saturated = True
G = 12.0  saturated = True
G = 14.0  saturated = False
G = 16.0  saturated = False

2. The background noise budget

Without sources the pixel variance should be Poisson on (sky + dark) × exptime plus the read-noise variance. Compare the measured standard deviation of image_e (electrons, before digitization) against that expectation.

sky = simulate_field(RA0, DEC0, sensorfilter='zwo:r', focus=2, exptime=EXPTIME,
                     seed=43, catalog=cat[:0], shape=SHAPE)   # no sources
bkg = sky.image_e
level_e = (p['sky_e_s'] + p['dark_e_s']) * EXPTIME
expected = np.sqrt(level_e + p['read_noise'] ** 2)
measured = bkg.std()
print(f'background level: {bkg.mean():8.2f} e-  (expected {level_e:8.2f})')
print(f'background std:   {measured:8.2f} e-  (expected {expected:8.2f})')
assert abs(measured / expected - 1) < 0.05

fig, ax = plt.subplots(figsize=(7.5, 4.2))
ax.hist(bkg.ravel(), bins=80, density=True, color=TEAL, alpha=0.6,
        label='pixels (source-free run)')
xx = np.linspace(bkg.min(), bkg.max(), 300)
ax.plot(xx, np.exp(-0.5 * ((xx - level_e) / expected) ** 2)
        / (expected * np.sqrt(2 * np.pi)), color=RED,
        label='Gaussian, ETC budget')
ax.set_xlabel('pixel value [e$^-$]'); ax.set_ylabel('density')
ax.legend()
plt.show()

/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: divide by zero encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: overflow encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: invalid value encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: divide by zero encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: overflow encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: invalid value encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total

background level:     1.51 e-  (expected     1.50)
background std:       3.27 e-  (expected     3.28)

3. n_reads: trading saturation for read noise

n_reads follows the ETC convention of coadded frames: the exposure is split into equal frames that are summed. Saturation is evaluated per frame, so more reads keep bright stars unsaturated, at the cost of read-noise variance growing as n_reads.

field9 = simulate_field(RA0, DEC0, sensorfilter='zwo:r', focus=2, exptime=EXPTIME,
                        n_reads=9, seed=42, catalog=cat, shape=SHAPE)
for f, label in [(field, 'n_reads = 1'), (field9, 'n_reads = 9')]:
    n_sat = int(f.saturation_mask.sum())
    stars = ', '.join(f'G={m:.0f}' for m, s in
                      zip(GMAGS, np.asarray(f.catalog['saturated'])) if s)
    bstd = (f.image_e - f.image_clean)[:, :200].std()   # noise realization only
    print(f'{label}: {n_sat:6d} saturated px | saturated stars: {stars or "none"} | '
          f'noise std {bstd:6.2f} e-')
assert field9.saturation_mask.sum() < field.saturation_mask.sum()

/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: divide by zero encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: overflow encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: invalid value encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: divide by zero encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: overflow encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: invalid value encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total

n_reads = 1:  20524 saturated px | saturated stars: G=6, G=8, G=10, G=12 | noise std   4.85 e-
n_reads = 9:   8764 saturated px | saturated stars: G=6, G=8, G=10 | noise std   9.87 e-

4. Noiseless closure

With add_noise=False the realisation is the expectation image itself, so image_e == image_clean wherever nothing saturates. The one exception is deliberate: the delivered frame respects the per-frame well ceiling with or without noise, so the saturated pixels of the G = 6 and G = 8 stars are capped in image_e and left at their true expectation in image_clean. This is the basis for photometric closure tests.

clean = simulate_field(RA0, DEC0, sensorfilter='zwo:r', focus=2, exptime=EXPTIME,
                       catalog=cat, shape=SHAPE, add_noise=False)
unsat = ~clean.saturation_mask
diff = np.abs(clean.image_e - clean.image_clean)
print(f'max |image_e - image_clean| off the saturated pixels: {diff[unsat].max()}')
print(f'saturated pixels: {int((~unsat).sum())}, capped at {clean.params["saturation_e"]:.0f} e- '
      f'(true expectation up to {clean.image_clean[~unsat].max():.3g} e-)')
assert diff[unsat].max() == 0.0
assert np.all(clean.image_e[~unsat] == clean.params['saturation_e'])

/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: divide by zero encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: overflow encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:35: RuntimeWarning: invalid value encountered in matmul
  cx = (psf.sum(axis=0) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: divide by zero encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: overflow encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total
/Users/gudmundurstefansson/Dropbox/mypylib/notebooks/GIT/wcc-sim/src/wcc_sim/psf.py:36: RuntimeWarning: invalid value encountered in matmul
  cy = (psf.sum(axis=1) @ i) / total

max |image_e - image_clean| off the saturated pixels: 0.0
saturated pixels: 20524, capped at 16275 e- (true expectation up to 6.07e+06 e-)

Summary

• Noise model = sky + dark, per-frame full-well/ADC cap, Poisson, read noise (variance × n_reads), gain + bias + clip to ADU, with the same semantics as the ETC.

• Defocused stars saturate in the ring first; the per-star flag uses a 32 px window.

• Background statistics match the ETC noise budget to a few percent.

• More reads → less per-frame saturation, more read noise; add_noise=False gives the exact expectation image.