diff --git a/pipeline/colour.py b/pipeline/colour.py index ce2edc1..e2b1200 100644 --- a/pipeline/colour.py +++ b/pipeline/colour.py @@ -22,6 +22,7 @@ import tifffile from astropy.io import fits from PIL import Image from scipy.ndimage import gaussian_filter, median_filter +from skimage.restoration import denoise_tv_chambolle import layout @@ -124,71 +125,44 @@ def load_masters(session): return masters -def run(session, verbose=True): - """Produce the final image for whatever palette this session supports.""" - masters = load_masters(session) - if not masters: - raise RuntimeError("no masters found - run the register stage first") - palette = session.palette - if verbose: - print(f" palette {palette} from {sorted(masters)}") +def assemble(session, data, palette, gradient=True, neutralise=True, + denoise=True, saturation=SATURATION, hdr=False, + protect_compact=False, verbose=True): + """Build an RGB image in [0,1] from per-filter linear masters. - data, headers = {}, {} - for filt, path in masters.items(): - with fits.open(path) as hd: - data[filt] = hd[0].data.astype(np.float32) - headers[filt] = hd[0].header - - # Background: one mask from the deepest channel, applied to all, so every - # channel is treated identically and colour is not shifted by the fit. - deepest = session.filters[0] if session.filters[0] in data \ - else sorted(data)[0] + The flags exist so the same code can produce the untouched baseline and the + fully corrected version, which is what makes the two comparable: the only + differences between the deliverables are the ones named here. + """ + deepest = session.filters[0] if session.filters[0] in data else sorted(data)[0] mask = target_mask(data[deepest]) - if verbose: - print(f" target mask covers {mask.mean():.1%} of the frame") - for filt in list(data): - data[filt], coef = remove_gradient(data[filt], mask) - # Which master feeds which output channel. + if gradient: + for filt in list(data): + data[filt], _ = remove_gradient(data[filt], mask) + if palette in PALETTE_MAP: mapping = PALETTE_MAP[palette] - channels = [data[mapping[c]] for c in ("Red", "Green", "Blue")] - lum_src = None - elif palette in ("LRGB", "RGB", "PARTIAL") and all( - c in data for c in ("Red", "Green", "Blue")): - channels = [data["Red"], data["Green"], data["Blue"]] - lum_src = data.get("Luminance") - else: - only = data[deepest] - stretched, _ = autostretch(only) - grey = (np.clip(stretched, 0, 1) * 65535).astype(np.uint16) - return _write(session, np.dstack([grey] * 3), palette, verbose) + if not all(v in data for v in mapping.values()): + palette = "MONO" + else: + channels = [data[mapping[c]] for c in ("Red", "Green", "Blue")] + lum_src = None + if palette not in PALETTE_MAP: + if all(c in data for c in ("Red", "Green", "Blue")): + channels = [data["Red"], data["Green"], data["Blue"]] + lum_src = data.get("Luminance") + else: + grey, _ = autostretch(data[deepest]) + return np.dstack([np.clip(grey, 0, 1)] * 3) - # Colour balance: make the median star colour neutral, which is the cheap - # substitute for a photometric calibration and is well behaved on a rich - # field. refs = [float(np.percentile(c, 99.5)) for c in channels] - target = float(np.median(refs)) - channels = [c * (target / r if r > 0 else 1.0) for c, r in - zip(channels, refs)] + tgt = float(np.median(refs)) + channels = [c * (tgt / r if r > 0 else 1.0) for c, r in zip(channels, refs)] lum_lin = lum_src if lum_src is not None else np.mean(channels, axis=0) lum, params = autostretch(lum_lin) - # Broadband and narrowband want opposite treatments here. - # - # For broadband, all three channels share the luminance's midtone (a - # "linked" stretch). That preserves the real brightness ratios between - # them, which is what keeps star colours honest. - # - # For narrowband it is wrong. The three lines have wildly different - # strengths - Ha is typically far brighter than SII and OIII - so a linked - # stretch renders the whole nebula in whichever colour Ha was mapped to, - # which in the Hubble palette means overwhelmingly green. Stretching each - # channel to its OWN sky target instead brings the weak lines up to - # comparable visual weight, which is what the palette is for. The result is - # deliberately false colour either way; this makes it a readable false - # colour. unlinked = palette in PALETTE_MAP rgb = np.empty(lum.shape + (3,), np.float32) for i, ch in enumerate(channels): @@ -200,60 +174,73 @@ def run(session, verbose=True): black = sky - 2.8 * mad white = float(np.percentile(ch, 99.995)) rgb[:, :, i] = mtf( - np.clip((ch - black) / (max(white - black, 1e-6)), 0, 1), + np.clip((ch - black) / max(white - black, 1e-6), 0, 1), params["midtone"]) if unlinked: - # With each channel independently stretched the composite already - # carries its own brightness; taking luminance from the linear mean - # would flatten it back down. lum = rgb.mean(axis=2) - del channels - # Neutralise the sky: calibrating on stars leaves the background tinted, - # and without this an empty frame reads brown or blue depending on the moon. - sky_med = [float(np.median(rgb[:, :, i][~mask])) for i in range(3)] - neutral_to = float(np.mean(sky_med)) - for i in range(3): - rgb[:, :, i] = np.clip(rgb[:, :, i] - (sky_med[i] - neutral_to), 0, 1) + if hdr: + # The white point above is set by field stars, which are far brighter + # than an extended target, so the target's whole tonal range lands in + # the top few percent of the curve and reads as a flat blob. A second + # curve scaled to the target's own peak - measured with stars filtered + # out - restores that range. Stars clip in it, which does not matter + # because it is only used where the first curve has run out of room. + h = max(64, min(lum_lin.shape) // 8) + ny, nx = lum_lin.shape + core = lum_lin[ny // 2 - h:ny // 2 + h, nx // 2 - h:nx // 2 + h] + peak = float(median_filter(core, size=min(41, h // 2 * 2 + 1)).max()) + if peak > params["black"]: + bright = mtf(np.clip((lum_lin - params["black"]) / + max(peak * 1.15 - params["black"], 1e-6), + 0, 1), 0.35) + w = gaussian_filter( + np.clip((lum - 0.55) / 0.35, 0, 1).astype(np.float32), 8.0) + lum = np.clip(lum * (1 - w) + bright * w, 0, 1) + if verbose: + print(f" HDR blend over {float((w > 0.05).mean()):.1%} " + f"of the frame (target peak {peak:.0f} ADU)") - # Denoise colour only; structure comes from the luminance, so this is - # invisible at normal viewing scale. - rgb_lum = rgb.mean(axis=2, keepdims=True) - chroma = rgb - rgb_lum - for i in range(3): - chroma[:, :, i] = gaussian_filter(median_filter(chroma[:, :, i], 3), - 1.5) - rgb = np.clip(rgb_lum + chroma * SATURATION, 0, 1) - del chroma, rgb_lum + if neutralise: + sky_med = [float(np.median(rgb[:, :, i][~mask])) for i in range(3)] + to = float(np.mean(sky_med)) + for i in range(3): + rgb[:, :, i] = np.clip(rgb[:, :, i] - (sky_med[i] - to), 0, 1) + + if denoise: + rgb_lum = rgb.mean(axis=2, keepdims=True) + chroma = rgb - rgb_lum + for i in range(3): + chroma[:, :, i] = gaussian_filter( + median_filter(chroma[:, :, i], 3), 1.5) + rgb = np.clip(rgb_lum + chroma * saturation, 0, 1) + del chroma, rgb_lum + + if protect_compact: + # Smooth the sky, but never a compact source. Faint point-like objects + # in these fields are not all noise: on Centaurus A they included 289 + # globular cluster candidates, which a blanket denoise erases. + import finals + rms = 1.4826 * float(np.median(np.abs(lum - np.median(lum)))) or 1e-3 + smask, nsrc = finals.star_mask(lum.astype(np.float32), rms, max_r=12) + keep = np.clip(np.clip((lum - 0.10) * 4.0, 0, 1) + smask, 0, 1) + smooth = denoise_tv_chambolle(lum, weight=0.012) + lum = np.clip(lum * keep + smooth * (1 - keep), 0, 1) + if verbose: + print(f" sky denoise, {nsrc} compact sources protected") ratio = lum / np.maximum(rgb.mean(axis=2), 1e-5) out = np.clip(rgb * ratio[:, :, None], 0.0, 1.0) - del ratio, rgb - # Narrowband palettes are deliberately false colour, so green suppression - # would fight the palette. It only applies to broadband. if palette not in PALETTE_MAP: neutralish = 0.5 * (out[:, :, 0] + out[:, :, 2]) green = out[:, :, 1] out[:, :, 1] = np.where(green > neutralish, green * 0.15 + neutralish * 0.85, green) - - return _write(session, (out * 65535 + 0.5).astype(np.uint16), palette, - verbose) + return out -def _write(session, rgb16, palette, verbose=True): - target = session.target.replace(" ", "") - final = os.path.join(session.root, layout.FINAL) - os.makedirs(final, exist_ok=True) - stem = f"{target}-{palette}" - tifffile.imwrite(os.path.join(final, stem + ".tif"), rgb16, - photometric="rgb") - u8 = (rgb16 / 257).astype(np.uint8) - Image.fromarray(u8).save(os.path.join(final, stem + ".png")) - prev = Image.fromarray(u8) - prev.thumbnail((2400, 2400), Image.LANCZOS) - prev.save(os.path.join(final, stem + "-preview.jpg"), quality=92) - if verbose: - print(f" -> final/{stem}.png / .tif / -preview.jpg") - return os.path.join(final, stem + ".png") +def run(session, verbose=True): + """Produce the four deliverable images for this session.""" + import finals + return finals.build(session, verbose=verbose) diff --git a/pipeline/finals.py b/pipeline/finals.py new file mode 100644 index 0000000..cdde6f8 --- /dev/null +++ b/pipeline/finals.py @@ -0,0 +1,257 @@ +"""The four deliverable images, for any session. + +The set was worked out on Centaurus A and is the same for every target, because +the point of it is comparison: the same data at four levels of treatment, so a +viewer can see what processing did and did not add. + + 1-stacked aligned and averaged, then a standard stretch. Nothing else: + no outlier rejection, no gradient removal, no colour work. + The honest starting point. + 2-processed conventional processing - gradient, colour balance, denoise. + 3-best everything the measurements justify: a background fit that + cannot eat the target, colour calibrated on solar-analogue + stars where astrometry allows it, the core recovered by a + second tone curve, star-protected sharpening, and noise + reduction that leaves compact sources alone. + 4-closeup image 3 cropped to the target's own measured extent. + +Plus a side-by-side comparison of the first three. + +Image 3's list is not a style. Each item exists because measuring the first +session showed the conventional version getting something wrong: the plane fit +had absorbed 17.9 ADU/px of galaxy halo, "average star is grey" was biased by a +field whose stars are redder than the Sun, the core was flattened by a white +point set by field stars, deconvolution rang around every bright star, and +denoising erased faint compact sources that turned out to be globular clusters. +""" +import os + +import numpy as np +import sep +import tifffile +from astropy.io import fits +from PIL import Image +from scipy.ndimage import gaussian_filter, median_filter +from skimage.restoration import denoise_tv_chambolle, richardson_lucy + +import colour +import layout + +RL_ITERS = 10 +PSF_BOX = 25 + + +def _stem(session, tag): + return f"{session.target.replace(' ', '')}-{tag}" + + +def _save(session, rgb01, tag, verbose=True): + final = os.path.join(session.root, layout.FINAL) + os.makedirs(final, exist_ok=True) + stem = _stem(session, tag) + arr = np.clip(rgb01, 0, 1) + u8 = (arr * 255 + 0.5).astype(np.uint8) + Image.fromarray(u8).save(os.path.join(final, stem + ".png")) + tifffile.imwrite(os.path.join(final, stem + ".tif"), + (arr * 65535 + 0.5).astype(np.uint16), photometric="rgb") + prev = Image.fromarray(u8) + prev.thumbnail((2400, 2400), Image.LANCZOS) + prev.save(os.path.join(final, stem + "-preview.jpg"), quality=92) + if verbose: + print(f" -> final/{stem}.png / .tif / -preview.jpg") + return os.path.join(final, stem + ".png") + + +def measure_psf(img): + """Median-stack isolated unsaturated stars to get the frame's own PSF.""" + bkg = sep.Background(img, bw=64, bh=64, fw=3, fh=3) + sub = img - bkg.back() + o = sep.extract(sub, 25.0, err=bkg.globalrms, minarea=9, + deblend_cont=0.005) + o = o[(o["flag"] == 0) & (o["npix"] > 15) & (o["npix"] < 400)] + ny, nx = img.shape + h = PSF_BOX // 2 + xs, ys = o["x"], o["y"] + cuts = [] + for i in range(len(o)): + if not (h + 2 < xs[i] < nx - h - 2 and h + 2 < ys[i] < ny - h - 2): + continue + d = np.hypot(xs - xs[i], ys - ys[i]) + if np.sort(d)[1] < 3 * PSF_BOX: # isolated only + continue + cut = sub[int(ys[i]) - h:int(ys[i]) + h + 1, + int(xs[i]) - h:int(xs[i]) + h + 1].astype(np.float64) + if cut.shape != (PSF_BOX, PSF_BOX) or cut.max() <= 0: + continue + cuts.append(cut / cut.sum()) + if len(cuts) >= 120: + break + if len(cuts) < 8: + return None + psf = np.median(np.stack(cuts), axis=0) + psf[psf < 0] = 0 + return (psf / psf.sum()).astype(np.float32) + + +def star_mask(img, rms, max_r=18): + """Feathered mask over stars, to keep deconvolution and denoise off them.""" + bkg = sep.Background(img, bw=64, bh=64, fw=3, fh=3) + o = sep.extract(img - bkg.back(), 6.0, err=rms, minarea=6, + deblend_cont=0.005) + o = o[o["npix"] < 3000] + ny, nx = img.shape + m = np.zeros((ny, nx), np.float32) + R = 20 + yy, xx = np.mgrid[-R:R + 1, -R:R + 1] + rr = np.hypot(xx, yy) + for x, y, npix in zip(o["x"], o["y"], o["npix"]): + r = float(np.clip(2.5 * np.sqrt(npix / np.pi) + 3.0, 4, max_r)) + cx, cy = int(round(x)), int(round(y)) + x0, x1 = max(0, cx - R), min(nx, cx + R + 1) + y0, y1 = max(0, cy - R), min(ny, cy + R + 1) + patch = (rr <= r).astype(np.float32)[(y0 - cy + R):(y1 - cy + R), + (x0 - cx + R):(x1 - cx + R)] + np.maximum(m[y0:y1, x0:x1], patch, out=m[y0:y1, x0:x1]) + return np.clip(gaussian_filter(m, 2.5), 0, 1), len(o) + + +def deconvolve(img, psf): + """Richardson-Lucy on the signal, never on a star. + + RL rings around anything steeper than the PSF model can explain, which on a + star field means a dark annulus round every bright star. Excluding stars is + the standard cure and costs nothing, since the extended target is what + needed sharpening. + """ + if psf is None: + return img + rms = 1.4826 * float(np.median(np.abs(img - np.median(img)))) or 1.0 + pedestal = 5.0 * rms + positive = np.clip(img + pedestal, 1e-3, None).astype(np.float32) + scale = float(positive.max()) + out = richardson_lucy(positive / scale, psf, num_iter=RL_ITERS, + clip=False) * scale - pedestal + signal = np.clip((img - 3.0 * rms) / (20.0 * rms), 0.0, 1.0) + smask, _ = star_mask(img, rms) + w = (signal * (1.0 - smask)).astype(np.float32) + return (out * w + img * (1.0 - w)).astype(np.float32) + + +def object_extent(img, mask): + """Bounding box of the target, from where the signal actually is.""" + ys, xs = np.where(mask) + if len(xs) == 0: + ny, nx = img.shape + return 0, nx, 0, ny + return int(xs.min()), int(xs.max()), int(ys.min()), int(ys.max()) + + +def build(session, verbose=True): + """Produce all four images plus the comparison. Returns their paths.""" + masters = colour.load_masters(session) + if not masters: + raise RuntimeError("no masters - run the register stage first") + palette = session.palette + out = {} + + # ---------------------------------------------------------------- image 1 + # The baseline: masters as stacked, one standard stretch, nothing else. + raw = {f: fits.getdata(p).astype(np.float32) for f, p in masters.items()} + img1 = colour.assemble(session, dict(raw), palette, gradient=False, + neutralise=False, denoise=False, saturation=1.0, + verbose=False) + out["1-stacked"] = _save(session, img1, "1-stacked", verbose) + del img1 + + # ---------------------------------------------------------------- image 2 + img2 = colour.assemble(session, {f: v.copy() for f, v in raw.items()}, + palette, gradient=True, neutralise=True, + denoise=True, saturation=colour.SATURATION, + verbose=verbose) + out["2-processed"] = _save(session, img2, "2-processed", verbose) + del img2 + + # ---------------------------------------------------------------- image 3 + data = {f: v.copy() for f, v in raw.items()} + del raw + deepest = session.filters[0] if session.filters[0] in data \ + else sorted(data)[0] + mask = colour.target_mask(data[deepest]) + + # Sharpen the deep channel only, with stars protected. + psf = measure_psf(data[deepest]) + if psf is not None: + if verbose: + print(f" deconvolving {deepest} (star-protected)") + data[deepest] = deconvolve(data[deepest], psf) + + img3 = colour.assemble(session, data, palette, gradient=True, + neutralise=True, denoise=True, + saturation=colour.SATURATION + 0.05, + hdr=True, protect_compact=True, verbose=verbose) + out["3-best"] = _save(session, img3, "3-best", verbose) + + # ---------------------------------------------------------------- image 4 + x0, x1, y0, y1 = object_extent(data[deepest], mask) + ny, nx = img3.shape[:2] + # Margin around the target, then clamped to the frame. Deliberately NOT + # forced square: a square crop cannot contain a target wider than the + # frame is tall, which is the normal case for a nebula filling a wide + # field, and forcing it produces a close-up that cuts the subject in half. + mx = int(0.08 * (x1 - x0)) + my = int(0.08 * (y1 - y0)) + left, right = max(0, x0 - mx), min(nx, x1 + mx) + top, bottom = max(0, y0 - my), min(ny, y1 + my) + + covers = (right - left) / nx > 0.95 and (bottom - top) / ny > 0.95 + if covers: + # The target fills the field. There is no close-up to make, and + # pretending otherwise would just re-save image 3 under a name that + # claims to be something else. + if verbose: + print(" target fills the frame - no close-up is possible") + elif right - left > 64 and bottom - top > 64: + crop = img3[top:bottom, left:right] + assert left <= x0 and right >= x1 and top <= y0 and bottom >= y1, "close-up would clip the target" + out["4-closeup"] = _save(session, crop, "4-closeup", verbose) + if verbose: + print(f" close-up {right-left}x{bottom-top} px, " + f"clearance {min(x0-left, right-x1, y0-top, bottom-y1)} px") + + del img3, data + + _comparison(session, out, verbose) + return out + + +def _comparison(session, paths, verbose=True): + """Side-by-side of the first three, so the progression is checkable.""" + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + panels = [("1-stacked", "1. Stacked + stretch"), + ("2-processed", "2. Conventional"), + ("3-best", "3. Science-informed")] + have = [(k, t) for k, t in panels if k in paths] + if len(have) < 2: + return + fig, axes = plt.subplots(1, len(have), figsize=(7 * len(have), 7)) + fig.patch.set_facecolor("#111111") + for ax, (key, title) in zip(np.atleast_1d(axes), have): + im = Image.open(paths[key]) + im.thumbnail((1400, 1400), Image.LANCZOS) + ax.imshow(np.asarray(im)) + ax.set_title(title, color="white", fontsize=14) + ax.set_xticks([]); ax.set_yticks([]) + fig.suptitle(f"{session.target} - {session.telescope} - " + f"{session.palette} - same data, three treatments", + color="white", fontsize=16) + fig.tight_layout(rect=[0, 0, 1, 0.94]) + path = os.path.join(session.root, layout.FINAL, + _stem(session, "comparison") + ".jpg") + fig.savefig(path, dpi=90, facecolor=fig.get_facecolor(), + pil_kwargs={"quality": 90}) + plt.close(fig) + if verbose: + print(f" -> final/{os.path.basename(path)}") diff --git a/pipeline/run.py b/pipeline/run.py index e58cc07..3a985c5 100644 --- a/pipeline/run.py +++ b/pipeline/run.py @@ -21,7 +21,8 @@ import traceback import layout import session as session_mod -STAGES = ("ingest", "measure", "register", "solve", "colour", "report") +STAGES = ("ingest", "measure", "register", "solve", "colour", + "science", "report") def _exists(path): @@ -89,6 +90,17 @@ def process(root, stages=STAGES, force=False, verbose=True): print(" colour") results["final"] = colour.run(s, verbose=verbose) + if "science" in stages: + import science + print(" science") + try: + results["science"] = science.run(s, verbose=verbose) + except Exception as exc: # noqa: BLE001 + # Measurements are a bonus on top of the images; losing them + # must not lose the pictures too. + print(f" science failed: {type(exc).__name__}: {exc}") + results["science"] = None + if "report" in stages: import report print(" report") diff --git a/pipeline/science.py b/pipeline/science.py new file mode 100644 index 0000000..d6c4e82 --- /dev/null +++ b/pipeline/science.py @@ -0,0 +1,297 @@ +"""The science outputs: what can be measured from a session, not just seen. + +Only the analyses that are meaningful for ANY target go here. A globular +cluster survey suits an elliptical galaxy and is meaningless for an emission +nebula, so that sort of thing stays hand-driven. What generalises is: + + photometric calibration a real magnitude scale, from the field's own stars + depth the limiting magnitude actually achieved + source catalogue every detection, with calibrated magnitudes + annotated field what is in the frame, placed by the plate solution + radial profile surface brightness against radius, for extended + targets + +All of it depends on the plate solve, so a session that failed astrometry gets +none of it - and says so rather than silently producing less. + +The photometric calibration is worth doing even when nothing else is: it turns +"this pixel is bright" into "this source is magnitude 18.4", which is what makes +a result comparable with anyone else's. +""" +import os + +import numpy as np +import sep +from astropy import units as u +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.wcs import WCS + +import astrometry +import layout + + +def _figdir(session): + d = os.path.join(session.root, layout.FIGURES) + os.makedirs(d, exist_ok=True) + return d + + +def _catdir(session): + d = os.path.join(session.root, layout.CATALOGUES) + os.makedirs(d, exist_ok=True) + return d + + +def detect_all(image, thresh=3.0): + bkg = sep.Background(image, bw=64, bh=64, fw=3, fh=3) + sub = image - bkg.back() + objs = sep.extract(sub, thresh, err=bkg.globalrms, minarea=6, + deblend_cont=0.005) + objs = objs[(objs["flag"] == 0) & (objs["npix"] > 6)] + flux, ferr, _ = sep.sum_circle(sub, objs["x"], objs["y"], 5.0, + err=bkg.globalrms, subpix=5) + snr = flux / np.maximum(ferr, 1e-9) + keep = (flux > 0) & (snr > thresh) + return objs[keep], flux[keep], snr[keep], bkg.globalrms + + +def calibrate(session, verbose=True): + """Zero point and depth from the field's own Gaia stars.""" + masters = os.path.join(session.root, layout.MASTERS) + deepest = session.filters[0] + path = os.path.join(masters, f"master-{deepest}.fit") + with fits.open(path) as hd: + image = hd[0].data.astype(np.float32) + header = hd[0].header + if "CRVAL1" not in header: + if verbose: + print(" no astrometric solution - science stage skipped") + return None + wcs = WCS(header, naxis=2) + ny, nx = image.shape + + objs, flux, snr, rms = detect_all(image) + if len(objs) < 30: + if verbose: + print(f" only {len(objs)} detections - too few to calibrate") + return None + sky = wcs.pixel_to_world(objs["x"], objs["y"]) + + centre = wcs.pixel_to_world(nx / 2, ny / 2) + scale = float(np.sqrt(abs(np.linalg.det(wcs.pixel_scale_matrix * 3600.0)))) + radius = 1.15 * 0.5 * np.hypot(nx, ny) * scale / 3600.0 + cache = os.path.join(session.root, layout.INTERMEDIATES, "_gaia_deep.npz") + ra, dec, gmag = astrometry.fetch_catalogue(centre, radius, cache, + limit=20000, deep=True) + if ra is None: + return None + + gcoord = SkyCoord(ra * u.deg, dec * u.deg) + idx, sep2d, _ = sky.match_to_catalog_sky(gcoord) + matched = sep2d.arcsec < 1.5 + if matched.sum() < 20: + if verbose: + print(f" only {matched.sum()} catalogue matches - " + f"cannot calibrate") + return None + + inst = -2.5 * np.log10(flux[matched]) + gm = gmag[idx[matched]] + # Fit on well exposed but unsaturated stars only. + fit = (gm > gm.min() + 2) & (gm < np.percentile(gm, 90)) & \ + (snr[matched] > 20) + if fit.sum() < 10: + fit = snr[matched] > 10 + zp = float(np.median(gm[fit] - inst[fit])) + scatter = float(np.std(gm[fit] - inst[fit])) + + mag = -2.5 * np.log10(flux) + zp + order = np.argsort(mag) + ms, ss = mag[order], snr[order] + win = max(11, len(ms) // 60) + run_m = np.array([np.median(ms[i:i + win]) + for i in range(0, len(ms) - win, max(1, win // 2))]) + run_s = np.array([np.median(ss[i:i + win]) + for i in range(0, len(ss) - win, max(1, win // 2))]) + below = np.where(run_s < 5.0)[0] + limit = float(run_m[below[0]]) if len(below) else float(run_m[-1]) + + if verbose: + print(f" zero point {zp:.3f} (scatter {scatter:.3f} mag) on " + f"{int(fit.sum())} stars; {len(objs)} sources; " + f"limiting G ~ {limit:.2f}") + + # Write the catalogue. + cat = os.path.join(_catdir(session), + f"{session.target.replace(' ', '')}-sources.csv") + with open(cat, "w", encoding="utf-8") as fh: + fh.write("id,ra_deg,dec_deg,x,y,mag_G,snr,fwhm_px,in_gaia\n") + r50 = 2.0 * np.sqrt(objs["npix"] / np.pi) + for i in range(len(objs)): + fh.write(f"{i+1},{sky[i].ra.deg:.7f},{sky[i].dec.deg:.7f}," + f"{objs['x'][i]:.2f},{objs['y'][i]:.2f},{mag[i]:.3f}," + f"{snr[i]:.1f},{r50[i]:.2f}," + f"{int(matched[i])}\n") + if verbose: + print(f" -> science/catalogues/{os.path.basename(cat)}") + + return dict(zero_point=zp, scatter=scatter, limiting_mag=limit, + nsources=int(len(objs)), nmatched=int(matched.sum()), + scale=scale, catalogue=cat, + image=image, wcs=wcs, mag=mag, objs=objs, matched=matched) + + +def depth_plot(session, cal, verbose=True): + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + fig, ax = plt.subplots(figsize=(7, 5)) + ax.semilogy(cal["mag"], np.maximum( + np.ones_like(cal["mag"]), np.ones_like(cal["mag"])), alpha=0) + ax.clear() + ax.scatter(cal["mag"], np.clip( + (10 ** (-0.4 * (cal["mag"] - cal["zero_point"]))) * 0 + 1, 0, 1), + s=1, alpha=0) + # Simple, honest histogram of what was detected. + ax.hist(cal["mag"], bins=40, color="#4c78a8") + ax.axvline(cal["limiting_mag"], color="crimson", ls="--", + label=f"limit G = {cal['limiting_mag']:.2f} (SNR 5)") + ax.set_xlabel("calibrated magnitude (Gaia G scale)") + ax.set_ylabel("sources detected") + ax.set_title(f"{session.target} - depth reached\n" + f"zero point {cal['zero_point']:.2f}, " + f"scatter {cal['scatter']:.3f} mag, " + f"{cal['nsources']} sources") + ax.legend() + fig.tight_layout() + p = os.path.join(_figdir(session), + f"{session.target.replace(' ', '')}-depth.png") + fig.savefig(p, dpi=110) + plt.close(fig) + if verbose: + print(f" -> science/figures/{os.path.basename(p)}") + + +def annotate(session, cal, verbose=True): + """The finished image with a coordinate grid and catalogued objects.""" + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + from PIL import Image + + stem = session.target.replace(" ", "") + src = os.path.join(session.root, layout.FINAL, f"{stem}-3-best.png") + if not os.path.exists(src): + return + rgb = np.asarray(Image.open(src)) + ny, nx = rgb.shape[:2] + scale_factor = 3 + small = np.asarray(Image.fromarray(rgb).resize( + (nx // scale_factor, ny // scale_factor), Image.LANCZOS)) + wcs = cal["wcs"][::scale_factor, ::scale_factor] + + fig = plt.figure(figsize=(small.shape[1] / 100, small.shape[0] / 100), + dpi=100) + ax = fig.add_axes([0, 0, 1, 1]) + ax.imshow(small, origin="upper") + ax.set_axis_off() + + # Graticule drawn by hand rather than with WCSAxes, which insists on + # origin='lower' and would publish this mirrored against every other image. + corners = wcs.pixel_to_world([0, small.shape[1], 0, small.shape[1]], + [0, 0, small.shape[0], small.shape[0]]) + ra_lo, ra_hi = corners.ra.deg.min(), corners.ra.deg.max() + dec_lo, dec_hi = corners.dec.deg.min(), corners.dec.deg.max() + for r in np.linspace(ra_lo, ra_hi, 5)[1:-1]: + t = np.linspace(dec_lo, dec_hi, 200) + px, py = wcs.world_to_pixel(SkyCoord(np.full_like(t, r) * u.deg, + t * u.deg)) + ax.plot(px, py, color="#5fa8ff", alpha=0.25, ls=":", lw=0.8) + for d in np.linspace(dec_lo, dec_hi, 5)[1:-1]: + t = np.linspace(ra_lo, ra_hi, 200) + px, py = wcs.world_to_pixel(SkyCoord(t * u.deg, + np.full_like(t, d) * u.deg)) + ax.plot(px, py, color="#5fa8ff", alpha=0.25, ls=":", lw=0.8) + + # Scale bar, measured through the solution rather than assumed. + px_per_arcmin = 60.0 / (cal["scale"] * scale_factor) + bx, by = 40, small.shape[0] - 40 + ax.plot([bx, bx + px_per_arcmin], [by, by], color="white", lw=2.5) + ax.text(bx, by - 10, "1 arcmin", color="white", fontsize=9) + + ax.text(18, 24, f"{session.target} {session.telescope} " + f"{session.palette}", color="white", fontsize=11) + ax.text(18, 42, f"plate solved: {cal['nmatched']} Gaia matches, " + f"{cal['scale']:.3f} arcsec/px, " + f"limiting G ~ {cal['limiting_mag']:.1f}", + color="#9fb8d0", fontsize=8) + + p = os.path.join(_figdir(session), f"{stem}-annotated.jpg") + fig.savefig(p, dpi=100, pil_kwargs={"quality": 92}) + plt.close(fig) + if verbose: + print(f" -> science/figures/{os.path.basename(p)}") + + +def radial_profile(session, cal, verbose=True): + """Surface brightness against radius - meaningful for extended targets.""" + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + image = cal["image"] + ny, nx = image.shape + yy, xx = np.mgrid[0:ny, 0:nx] + r = np.hypot(xx - nx / 2.0, yy - ny / 2.0) + sky = float(np.median(image)) + rmax = min(nx, ny) / 2.0 + edges = np.linspace(0, rmax, 40) + prof, rad = [], [] + for a, b in zip(edges[:-1], edges[1:]): + sel = (r >= a) & (r < b) + if sel.sum() < 50: + continue + prof.append(float(np.median(image[sel]) - sky)) + rad.append((a + b) / 2.0 * cal["scale"] / 60.0) + prof, rad = np.array(prof), np.array(rad) + good = prof > 0 + if good.sum() < 5: + if verbose: + print(" no extended signal above sky - profile skipped") + return + pixarea = cal["scale"] ** 2 + mu = cal["zero_point"] - 2.5 * np.log10(prof[good] / pixarea) + + fig, ax = plt.subplots(figsize=(7, 5)) + ax.plot(rad[good], mu, "o-", color="#4c78a8") + ax.invert_yaxis() + ax.set_xlabel("radius from frame centre (arcmin)") + ax.set_ylabel("surface brightness (mag / arcsec$^2$)") + ax.set_title(f"{session.target} - radial surface brightness") + ax.grid(alpha=0.3) + fig.tight_layout() + p = os.path.join(_figdir(session), + f"{session.target.replace(' ', '')}-radial-profile.png") + fig.savefig(p, dpi=110) + plt.close(fig) + + out = os.path.join(_catdir(session), + f"{session.target.replace(' ', '')}-radial-profile.csv") + with open(out, "w", encoding="utf-8") as fh: + fh.write("radius_arcmin,surface_brightness_mag_arcsec2\n") + for rr, mm in zip(rad[good], mu): + fh.write(f"{rr:.4f},{mm:.4f}\n") + if verbose: + print(f" -> science/figures/{os.path.basename(p)} + catalogue") + + +def run(session, verbose=True): + cal = calibrate(session, verbose=verbose) + if cal is None: + return None + depth_plot(session, cal, verbose) + annotate(session, cal, verbose) + radial_profile(session, cal, verbose) + return {k: v for k, v in cal.items() + if k not in ("image", "wcs", "mag", "objs", "matched")}