Merge pull request 'solve: prefer the estate's own Astrometry.net over nova for blind solves' (#2) from astrometry-net-fallback into main

The estate now runs Astrometry.net itself (ankh-morpork-infra astrometry/):
solve-field against ~5 GB of local index files, sized from this repo's
TELESCOPES.md so the whole iTelescope fleet's fields of view are covered.

blind.py now tries it before nova.astrometry.net. Speed is the least of the
reasons - a blind solve of a real DSS field returns in 0.7s where nova queues
for minutes. The reasons that matter are that nova requires UPLOADING the
master to a third party and holding an API key, and neither is necessary any
more for the ordinary case. nova remains the fallback, so an estate outage
costs speed and privacy rather than the ability to solve.

pipeline/astrometry_net.py is the client and nothing more: stdlib-only POST,
parses the returned .wcs so the full TAN solution is used rather than a
re-derivation from the summary numbers. It deliberately does NOT decide whether
to trust a solution - blind.py already verifies every blind solve against Gaia
and applies a star-count and residual gate, and two gates that can disagree is
worse than one that is trusted.

The source ('estate' or 'nova') is threaded through the log lines, the returned
method, and the ASTRSOLV card, because a year from now that card is the only
way to tell whether a frame was solved in-house or uploaded.

Also corrected astrometry.py's 'cannot solve without a blind solver' message,
which has been untrue since run.py started falling through to blind.py.

Tested against the live service: blind solve of a DSS2 field with a known
centre returned within ~4 arcsec in 0.7s; the WCS round-trips through astropy;
and an unreachable service falls through to nova instead of raising.

pipeline/astrometry.py

@@ -0,0 +1,424 @@
+"""Stage 4: plate solve the deepest master and copy the WCS to the others.
+
+iTelescope's frames arrive with a PinPoint HISTORY line but no WCS keywords at
+all, so the astrometry has to be redone locally. A blind solve is unnecessary:
+the header gives the pointing to arcminutes and usually the plate scale to four
+figures, so a catalogue can be fetched for the right patch of sky and matched
+directly.
+
+The match is asterism-based, which is invariant to rotation and scale - so the
+camera angle never has to be guessed, and a session at an unknown orientation
+solves as easily as one at a known angle. It is NOT invariant to a mirror flip,
+because a similarity transform cannot include reflection, so both parities are
+tried and the one that matches wins.
+
+The seed match returns only a handful of stars. That is enough to establish an
+approximate solution but too few to trust, so the solution is then used to pair
+every catalogue source with its nearest detection and refit, twice, with a
+shrinking tolerance. The residual that comes out of that is the honest measure
+of whether the solve is real.
+
+Failure here is not fatal to the pipeline: an unsolved session still produces
+images, it just cannot produce positions.
+"""
+import os
+
+import astroalign as aa
+import numpy as np
+from astropy import units as u
+from astropy.coordinates import SkyCoord
+from astropy.io import fits
+from astropy.wcs import WCS
+from astropy.wcs.utils import fit_wcs_from_points
+
+import layout
+import measure as measure_mod
+
+SEED_STARS = 600
+CATALOGUE_FETCH = 6000   # rich fields need far more than the match uses
+
+
+def pointing(header):
+    """Approximate field centre from whatever the header offers."""
+    for ra_key, dec_key in (("OBJCTRA", "OBJCTDEC"), ("RA", "DEC"),
+                            ("CRVAL1", "CRVAL2")):
+        ra, dec = header.get(ra_key), header.get(dec_key)
+        if ra is None or dec is None:
+            continue
+        try:
+            if isinstance(ra, str) and (":" in ra or " " in ra):
+                return SkyCoord(ra, dec, unit=(u.hourangle, u.deg))
+            return SkyCoord(float(ra) * u.deg, float(dec) * u.deg)
+        except Exception:                                     # noqa: BLE001
+            continue
+    return None
+
+
+def plate_scale(header, session=None):
+    """Arcsec per pixel: from the header if present, else from the optics."""
+    val = header.get("HIERARCH iTelescopePlateScaleH")
+    if val:
+        try:
+            return float(val)
+        except (TypeError, ValueError):
+            pass
+    focal = header.get("FOCALLEN")
+    pixel = header.get("XPIXSZ")
+    if focal and pixel:
+        # 206265 arcsec per radian; XPIXSZ is already binned microns.
+        return 206.265 * float(pixel) / float(focal)
+    return session.scale if session else None
+
+
+def fetch_catalogue(centre, radius_deg, cache, limit=SEED_STARS, deep=False):
+    """Gaia sources around the pointing, cached so a re-run costs nothing."""
+    if os.path.exists(cache):
+        z = np.load(cache)
+        return z["ra"], z["dec"], z["g"]
+    ra = dec = g = None
+    mag_cut = 20.5 if deep else 18.0
+
+    # VizieR first, deliberately. It serves the identical DR3 catalogue over a
+    # plain HTTP query and answers in under a second. Gaia's own archive is
+    # excellent but launch_job_async submits to a job QUEUE and polls for the
+    # result, which for a few hundred rows means waiting minutes for work that
+    # takes milliseconds - and it hung for ten minutes on a query of 600 stars
+    # while the network was demonstrably healthy.
+    try:
+        from astroquery.vizier import Vizier
+        # Ask for far more rows than are needed and pick the brightest
+        # locally. VizieR's row_limit truncates the result BEFORE any sort, so
+        # asking for 600 rows in a crowded field returns 600 arbitrary stars
+        # rather than the 600 brightest - and asterism matching only works if
+        # both lists contain the same bright stars. In the LMC field of
+        # NGC 2030 that silently produced no match at all.
+        v = Vizier(columns=["RA_ICRS", "DE_ICRS", "Gmag"],
+                   column_filters={"Gmag": f"<{mag_cut}"}, row_limit=50000)
+        v.TIMEOUT = 60
+        res = v.query_region(centre, radius=radius_deg * u.deg,
+                             catalog="I/355/gaiadr3")
+        if res:
+            t = res[0]
+            ra = np.asarray(t["RA_ICRS"], float)
+            dec = np.asarray(t["DE_ICRS"], float)
+            g = np.asarray(t["Gmag"], float)
+    except Exception as exc:                                  # noqa: BLE001
+        print(f"        VizieR unavailable ({type(exc).__name__})")
+
+    if ra is None or len(ra) == 0:
+        try:
+            from astroquery.gaia import Gaia
+            Gaia.ROW_LIMIT = limit
+            # Synchronous, not async: no job queue for a query this small.
+            job = Gaia.launch_job(f"""
+                SELECT TOP {limit} ra, dec, phot_g_mean_mag
+                FROM gaiadr3.gaia_source
+                WHERE 1 = CONTAINS(POINT('ICRS', ra, dec),
+                         CIRCLE('ICRS', {centre.ra.deg}, {centre.dec.deg},
+                                {radius_deg}))
+                  AND phot_g_mean_mag IS NOT NULL
+                  AND phot_g_mean_mag < {mag_cut}
+                ORDER BY phot_g_mean_mag ASC""")
+            t = job.get_results()
+            ra = np.asarray(t["ra"], float)
+            dec = np.asarray(t["dec"], float)
+            g = np.asarray(t["phot_g_mean_mag"], float)
+        except Exception as exc:                              # noqa: BLE001
+            print(f"        Gaia archive unavailable ({type(exc).__name__})")
+            return None, None, None
+
+    order = np.argsort(g)[:limit]
+    ra, dec, g = ra[order], dec[order], g[order]
+    np.savez_compressed(cache, ra=ra, dec=dec, g=g)
+    return ra, dec, g
+
+
+def project(ra, dec, centre, scale, parity):
+    """Gnomonic projection to pixel-like coordinates for asterism matching."""
+    c = SkyCoord(ra * u.deg, dec * u.deg)
+    dx, dy = centre.spherical_offsets_to(c)
+    return np.column_stack([dx.to_value(u.arcsec) / scale * parity,
+                            dy.to_value(u.arcsec) / scale])
+
+
+
+def _roll_angle(header):
+    """Camera angle from the header, if the telescope recorded it."""
+    for key in ("HIERARCH RollAngle", "RollAngle", "POSANG", "CROTA2"):
+        val = header.get(key)
+        if val is None:
+            continue
+        try:
+            return float(val)
+        except (TypeError, ValueError):
+            continue
+    return None
+
+
+def match_by_pointing(src, ra, dec, centre, scale, roll, shape, verbose=True):
+    """Pair detections with catalogue stars using a known camera angle.
+
+    Asterism matching exists to recover an UNKNOWN orientation. When the header
+    records the roll angle - and iTelescope's does - the whole geometry is
+    already known bar a small pointing error, and a direct search is both more
+    robust and far faster. Rotation, scale and parity are applied, then the
+    residual translation is found by histogramming every detection-to-catalogue
+    offset and taking the peak: a real solution puts thousands of pairs in one
+    bin, and noise spreads flat.
+    """
+    ny, nx = shape
+    # Slide down the catalogue's magnitude ranking as well as trying the four
+    # geometric conventions. In a field as rich as the LMC the brightest few
+    # hundred catalogue stars are far brighter than anything a 60 second
+    # narrowband frame records, so the two "brightest N" lists describe
+    # different populations and overlap barely at all.
+    best = None
+    slices = [(0, 400), (200, 800), (600, 1400), (1200, 2400), (2000, 4000)]
+    for lo, hi in slices:
+        if lo >= len(ra):
+            break
+        ra_s, dec_s = ra[lo:hi], dec[lo:hi]
+        for parity in (1.0, -1.0):
+          for sign in (1.0, -1.0):
+            th = np.radians(sign * roll)
+            cat = project(ra_s, dec_s, centre, scale, parity)
+            rot = np.column_stack([
+                cat[:, 0] * np.cos(th) - cat[:, 1] * np.sin(th),
+                cat[:, 0] * np.sin(th) + cat[:, 1] * np.cos(th)])
+            cat_px = rot + np.array([nx / 2.0, ny / 2.0])
+
+            # Every pairwise offset, histogrammed. Limited to the brightest of
+            # each list to keep this to a few million comparisons.
+            a = src[:400]
+            b = cat_px[:400]
+            dx = (a[:, None, 0] - b[None, :, 0]).ravel()
+            dy = (a[:, None, 1] - b[None, :, 1]).ravel()
+            lim = 0.25 * max(nx, ny)
+            keep = (np.abs(dx) < lim) & (np.abs(dy) < lim)
+            if keep.sum() < 50:
+                continue
+            bins = np.arange(-lim, lim + 12, 12)
+            H, xe, ye = np.histogram2d(dx[keep], dy[keep], bins=[bins, bins])
+            iy, ix = np.unravel_index(np.argmax(H), H.shape)
+            peak = H[iy, ix]
+            if best is None or peak > best[0]:
+                off = (0.5 * (xe[iy] + xe[iy + 1]),
+                       0.5 * (ye[ix] + ye[ix + 1]))
+                best = (peak, parity, sign, off, cat_px, lo, hi)
+    if best is None:
+        return None
+    peak, parity, sign, off, cat_px, lo, hi = best
+    if verbose:
+        print(f"        pointing match: catalogue[{lo}:{hi}], parity "
+              f"{parity:+.0f}, roll {sign:+.0f}, offset "
+              f"({off[0]:+.0f}, {off[1]:+.0f}) px, peak {int(peak)} pairs")
+    if peak < 12:
+        return None
+    shifted = cat_px + np.array(off)
+    idx_base = lo
+    pairs = []
+    for i, (cx, cy) in enumerate(shifted):
+        if not (0 <= cx < nx and 0 <= cy < ny):
+            continue
+        d = np.hypot(src[:, 0] - cx, src[:, 1] - cy)
+        j = int(np.argmin(d))
+        if d[j] < 18:
+            pairs.append((i + idx_base, j))
+    if len(pairs) < 10:
+        return None
+    return np.array([p[0] for p in pairs]), np.array([p[1] for p in pairs])
+
+
+def run(session, verbose=True):
+    """Solve the deepest master; write the WCS into every master."""
+    masters_dir = os.path.join(session.root, layout.MASTERS)
+    deepest = session.filters[0]
+    path = os.path.join(masters_dir, f"master-{deepest}.fit")
+    if not os.path.exists(path):
+        raise RuntimeError("no masters found - run the register stage first")
+
+    with fits.open(path) as hd:
+        image = hd[0].data.astype(np.float32)
+        header = hd[0].header
+    ny, nx = image.shape
+
+    centre = pointing(header)
+    scale = plate_scale(header, session)
+    if centre is None or not scale:
+        # Not fatal, and no longer the end of the road: returning None sends
+        # run.py to blind.py, which needs neither the pointing nor the scale.
+        print("        no pointing or plate scale in the header - "
+              "handing off to the blind solver")
+        return None
+    # Match over the frame's INSCRIBED circle, not its circumscribed one. A
+    # cone big enough to reach the corners also reaches well outside the
+    # frame - on a 4 degree field that is 33 square degrees of catalogue
+    # against 16 of image, so half the catalogue's brightest stars are not in
+    # the picture at all and the two "brightest N" lists barely overlap.
+    # The inscribed circle is inside the frame whatever the camera angle, so
+    # both lists then describe the same patch of sky.
+    radius = 0.95 * 0.5 * min(nx, ny) * scale / 3600.0
+    if verbose:
+        print(f"        pointing {centre.to_string('hmsdms')}, "
+              f"{scale:.4f} arcsec/px, search {radius:.3f} deg")
+
+    xy, _ = measure_mod.load(session)
+    # Detect on the master itself: it is deeper than any single frame.
+    from measure import measure_frame
+    det = measure_frame(path, max_stars=SEED_STARS * 3)
+    src = det["xy"]
+    # Same restriction on the detections, so the two lists cover one region.
+    r_pix = 0.95 * 0.5 * min(nx, ny)
+    inside = np.hypot(src[:, 0] - nx / 2.0, src[:, 1] - ny / 2.0) < r_pix
+    src = src[inside][:SEED_STARS]
+    if len(src) < 12:
+        print(f"        only {len(src)} stars in the master - too few")
+        return None
+
+    cache = os.path.join(session.root, layout.INTERMEDIATES, "_gaia.npz")
+    os.makedirs(os.path.dirname(cache), exist_ok=True)
+    ra, dec, gmag = fetch_catalogue(centre, radius, cache,
+                                    limit=CATALOGUE_FETCH)
+    if ra is None:
+        print("        no catalogue available - skipping the solve")
+        return None
+    if verbose:
+        print(f"        {len(ra)} catalogue stars, {len(src)} detected")
+
+    # If the header recorded the camera angle, use it: a direct geometric
+    # match is far more reliable than recovering the orientation from scratch.
+    roll = _roll_angle(header)
+    seed = None
+    if roll is not None:
+        if verbose:
+            print(f"        header roll angle {roll:.2f} deg")
+        seed = match_by_pointing(src, ra, dec, centre, scale, roll,
+                                 (ny, nx), verbose=verbose)
+    if seed is not None:
+        ci, ii = seed
+        world0 = SkyCoord(ra[ci] * u.deg, dec[ci] * u.deg)
+        wcs = fit_wcs_from_points((src[ii, 0], src[ii, 1]), world0,
+                                  proj_point="center", projection="TAN")
+        best = "pointing"
+    else:
+        best = None
+
+    # Matching needs the two lists to contain the SAME stars, and "brightest"
+    # does not guarantee that. The detector deliberately rejects saturated
+    # cores, so on a well exposed frame the brightest detections are not the
+    # brightest stars in the sky - while the catalogue's brightest are. The
+    # two top-N lists can therefore barely overlap, which is exactly what
+    # happened on the crowded LMC field of NGC 2030.
+    #
+    # So slide a window down the catalogue's magnitude ranking as well as
+    # varying the sample size, and take the best match found.
+    if best == "pointing":
+        pass
+    else:
+      best = None
+      for offset in (0, 15, 40, 80):
+        for n in (120, 200, 60):
+            hi = min(offset + n, len(ra))
+            if hi - offset < 25 or len(src) < 25:
+                continue
+            cat_slice = slice(offset, hi)
+            for parity in (1.0, -1.0):
+                cat_xy = project(ra[cat_slice], dec[cat_slice], centre, scale,
+                                 parity)
+                try:
+                    tform, (cat_m, img_m) = aa.find_transform(
+                        cat_xy, src[:max(n, 120)], max_control_points=60)
+                except Exception:                             # noqa: BLE001
+                    continue
+                if best is None or len(cat_m) > len(best[0]):
+                    # Record which catalogue rows matched, in FULL-list terms.
+                    best = (cat_m, img_m, parity, offset)
+                    if verbose:
+                        print(f"        catalogue[{offset}:{hi}] parity "
+                              f"{parity:+.0f}: matched {len(cat_m)} stars")
+            if best is not None and len(best[0]) >= 8:
+                break
+        if best is not None and len(best[0]) >= 8:
+            break
+
+    if best is None:
+        print("        no match by pointing or asterism - solve failed")
+        return None
+
+    if best != "pointing":
+        cat_m, img_m, parity, _offset = best
+        cat_xy = project(ra, dec, centre, scale, parity)
+        idx = [int(np.argmin(np.hypot(cat_xy[:, 0] - px, cat_xy[:, 1] - py)))
+               for px, py in cat_m]
+        world = SkyCoord(ra[idx] * u.deg, dec[idx] * u.deg)
+        wcs = fit_wcs_from_points((img_m[:, 0], img_m[:, 1]), world,
+                                  proj_point="center", projection="TAN")
+
+    # Refine: the seed match is a handful of stars, so use the approximate
+    # solution to pair every catalogue source with its nearest detection.
+    all_world = SkyCoord(ra * u.deg, dec * u.deg)
+    nmatch, resid = (len(seed[0]) if seed is not None else len(cat_m)), None
+    for tol in (4.0, 2.0):
+        px, py = wcs.world_to_pixel(all_world)
+        pairs = []
+        for i, (cx, cy) in enumerate(zip(px, py)):
+            if not (0 <= cx < nx and 0 <= cy < ny):
+                continue
+            d = np.hypot(src[:, 0] - cx, src[:, 1] - cy)
+            j = int(np.argmin(d))
+            if d[j] <= tol:
+                pairs.append((i, j))
+        if len(pairs) < 20:
+            break
+        ci = np.array([p[0] for p in pairs])
+        ii = np.array([p[1] for p in pairs])
+        wcs = fit_wcs_from_points((src[ii, 0], src[ii, 1]), all_world[ci],
+                                  proj_point="center", projection="TAN")
+        qx, qy = wcs.world_to_pixel(all_world[ci])
+        resid = np.hypot(qx - src[ii, 0], qy - src[ii, 1])
+        nmatch = len(pairs)
+
+    # A quality gate, because a WRONG solution is worse than none: every
+    # position in the science catalogue would inherit the error, silently. A
+    # genuine solve on a field like this matches hundreds of stars and refines
+    # to well under a pixel; a spurious alignment matches a couple of dozen and
+    # stalls at several. Both conditions must hold.
+    MIN_STARS, MAX_RESID_PX = 40, 1.5
+    if resid is None or nmatch < MIN_STARS or np.median(resid) > MAX_RESID_PX:
+        got = "n/a" if resid is None else f"{np.median(resid):.2f} px"
+        print(f"        REJECTED: {nmatch} stars at {got} residual "
+              f"(need >= {MIN_STARS} stars and <= {MAX_RESID_PX} px). "
+              f"Treating as unsolved rather than trusting a doubtful WCS.")
+        return None
+
+    cd = wcs.pixel_scale_matrix * 3600.0
+    solved_scale = float(np.sqrt(abs(np.linalg.det(cd))))
+    rot = float(np.degrees(np.arctan2(cd[0, 1], cd[1, 1])))
+    med_resid = float(np.median(resid)) if resid is not None else float("nan")
+    field = wcs.pixel_to_world(nx / 2.0, ny / 2.0)
+
+    if verbose:
+        print(f"        solved on {nmatch} stars, residual "
+              f"{med_resid:.2f} px ({med_resid * solved_scale:.2f} arcsec)")
+        print(f"        centre {field.to_string('hmsdms')}, "
+              f"{solved_scale:.4f} arcsec/px, PA {rot:.2f} deg, "
+              f"{nx * solved_scale / 60:.1f}' x {ny * solved_scale / 60:.1f}'")
+
+    whdr = wcs.to_header()
+    for filt in session.filters:
+        p = os.path.join(masters_dir, f"master-{filt}.fit")
+        if not os.path.exists(p):
+            continue
+        with fits.open(p, mode="update") as hd:
+            for card in whdr.cards:
+                hd[0].header[card.keyword] = (card.value, card.comment)
+            hd[0].header["ASTRSOLV"] = (
+                f"Gaia DR3 / {nmatch} stars / "
+                f"{med_resid * solved_scale:.2f} arcsec", "local plate solve")
+            hd.flush()
+
+    return dict(nstars=nmatch, residual_px=med_resid,
+                residual_arcsec=med_resid * solved_scale,
+                scale=solved_scale, position_angle=rot,
+                centre=field.to_string("hmsdms"),
+                fov_arcmin=(nx * solved_scale / 60, ny * solved_scale / 60))

pipeline/astrometry_net.py

@@ -0,0 +1,142 @@
+"""Client for the estate's own Astrometry.net service.
+
+    https://astrometry.ankh-morpork.discworld.network
+
+This is solve-field running on docker1/2/3 against ~5 GB of local index files.
+It exists so that blind solving costs seconds instead of minutes and, more to
+the point, so that no image has to leave the house to get solved.
+
+It is used by blind.py, which prefers this service and falls back to the hosted
+nova.astrometry.net when it is unreachable. Everything a solution has to earn -
+verification against Gaia, the star-count and residual gate, writing the WCS
+into the masters - stays in blind.py, shared by both sources. This module only
+gets a WCS out of the estate service; it does not decide whether to believe it.
+
+Only the stdlib is used for the request. A solve is one POST; pulling in
+requests to do it would add a dependency for no gain.
+"""
+import base64
+import io
+import os
+import ssl
+import time
+import urllib.error
+import urllib.request
+
+from astropy.io import fits
+from astropy.wcs import WCS
+
+SERVICE = os.environ.get("ASTROMETRY_SERVICE",
+                         "https://astrometry.ankh-morpork.discworld.network")
+TIMEOUT = int(os.environ.get("ASTROMETRY_TIMEOUT", "180"))
+
+
+def _context():
+    """TLS context for the estate service.
+
+    The service presents a step-ca certificate from the estate's internal CA.
+    On a machine that trusts the Discworld root CA this verifies normally. On
+    one that does not, verification fails - and the right fix is to install the
+    CA, not to turn checking off, so that is what the error message says.
+    ASTROMETRY_CA_BUNDLE points at the root CA if it is not in the system store.
+    ASTROMETRY_INSECURE=1 exists as an explicit, visible escape hatch; it is
+    never the default and it prints when used.
+    """
+    bundle = os.environ.get("ASTROMETRY_CA_BUNDLE")
+    if os.environ.get("ASTROMETRY_INSECURE") == "1":
+        print("        WARNING: TLS verification disabled "
+              "(ASTROMETRY_INSECURE=1)")
+        ctx = ssl.create_default_context()
+        ctx.check_hostname = False
+        ctx.verify_mode = ssl.CERT_NONE
+        return ctx
+    return ssl.create_default_context(cafile=bundle) if bundle else \
+        ssl.create_default_context()
+
+
+def available(timeout=5):
+    """Is the service up and does it have indexes? Cheap enough to ask first."""
+    try:
+        with urllib.request.urlopen(SERVICE + "/healthz", timeout=timeout,
+                                    context=_context()) as resp:
+            return resp.status == 200
+    except Exception:                                         # noqa: BLE001
+        return False
+
+
+def solve(path, centre=None, scale=None, radius_deg=2.0, timeout=TIMEOUT,
+          verbose=True):
+    """Solve one FITS file. Returns (WCS, info dict) or (None, reason).
+
+    `centre` and `scale` are hints, not requirements - pass whatever the header
+    offered. A hinted solve turns a whole-sky search into a cone and is an order
+    of magnitude faster, but the whole point of this fallback is that it still
+    works when there is nothing to hint with.
+    """
+    params = []
+    if centre is not None:
+        params += [f"ra={centre.ra.deg:.6f}", f"dec={centre.dec.deg:.6f}",
+                   f"radius={radius_deg:.3f}"]
+    if scale:
+        # +/-20% around the nominal scale: wide enough to absorb a focal
+        # reducer or a binning the header did not mention, tight enough to
+        # still rule out most index scales.
+        params += [f"scale_low={scale * 0.8:.4f}",
+                   f"scale_high={scale * 1.2:.4f}"]
+    params.append(f"timeout={timeout}")
+    url = f"{SERVICE}/solve?" + "&".join(params)
+
+    with open(path, "rb") as handle:
+        body = handle.read()
+
+    started = time.time()
+    if verbose:
+        hint = "hinted" if (centre is not None or scale) else "blind"
+        print(f"        astrometry.net ({hint}): "
+              f"{len(body) / 1e6:.1f} MB to {SERVICE}")
+    request = urllib.request.Request(
+        url, data=body, method="POST",
+        headers={"Content-Type": "application/octet-stream"})
+    try:
+        with urllib.request.urlopen(request, timeout=timeout + 30,
+                                    context=_context()) as resp:
+            import json
+            payload = json.load(resp)
+    except urllib.error.HTTPError as exc:
+        # 422 is the service saying "no solution", which is an answer, not a
+        # fault. Anything else is the service or the network being broken, and
+        # the two should not read the same in the log.
+        import json
+        try:
+            detail = json.load(exc)
+        except Exception:                                     # noqa: BLE001
+            detail = {}
+        if exc.code == 422:
+            return None, "no solution found"
+        return None, f"service error HTTP {exc.code}: {detail.get('error', '')}"
+    except ssl.SSLCertVerificationError:
+        return None, ("TLS verification failed - install the Discworld root CA "
+                      "or set ASTROMETRY_CA_BUNDLE to it")
+    except Exception as exc:                                  # noqa: BLE001
+        return None, f"{type(exc).__name__}: {exc}"
+
+    if not payload.get("solved"):
+        return None, "no solution found"
+
+    # The service hands back the .wcs file - a bare FITS header - base64
+    # encoded. Parsing it here rather than trusting the summary numbers means
+    # the full distortion-free TAN solution is what gets written to the
+    # masters, not a re-derivation from ra/dec/scale.
+    with fits.open(io.BytesIO(base64.b64decode(payload["wcs_b64"]))) as hdul:
+        wcs = WCS(hdul[0].header)
+
+    info = {k: payload.get(k) for k in
+            ("ra", "dec", "pixscale", "orientation", "parity",
+             "fieldw", "fieldh", "fieldunits")}
+    info["seconds"] = payload.get("seconds", round(time.time() - started, 1))
+    if verbose:
+        print(f"        astrometry.net solved in {info['seconds']}s: "
+              f"centre {info['ra']:.5f} {info['dec']:+.5f}, "
+              f"{info['pixscale']:.3f} arcsec/px, "
+              f"{info['fieldw']:.1f}' x {info['fieldh']:.1f}'")
+    return wcs, info

pipeline/blind.py

@@ -0,0 +1,243 @@
+"""Blind plate solving, for fields the seeded solver cannot handle.
+
+The seeded solver in `astrometry.py` starts from the header's pointing, scale
+and roll angle and matches against Gaia. That works whenever the image and the
+catalogue contain recognisably the same stars. It fails when they do not - a
+single short narrowband frame over a four degree field records a sparse,
+shallow star population that overlaps poorly with any magnitude slice of Gaia,
+and no amount of window-sliding fixes a population mismatch.
+
+A blind solver does not care. It builds geometric hashes from the image's own
+stars and looks them up in pre-built index files, so it needs no pointing, no
+scale and no orientation - only the pixels. astrometry.net is the standard
+implementation, and there are now two of it available:
+
+1. **The estate's own** (`astrometry_net.py`) - solve-field on docker1/2/3
+   against local index files. Tried first: the image never leaves the house, it
+   needs no API key, and it answers in seconds rather than queueing. Its index
+   set covers 5.6 arcmin to 33 degree quads, which is the whole iTelescope
+   fleet's range, but NOT arbitrarily narrow fields.
+2. **The hosted nova.astrometry.net**, through astroquery. The fallback, for
+   when the estate service is unreachable or its index set cannot reach the
+   field. Two things to know about it:
+
+   * **It needs an API key.** Free, from an account at
+     https://nova.astrometry.net (Profile -> API Key). Provide it as the
+     ASTROMETRY_API_KEY environment variable, or in a file named
+     `astrometry.key` beside the session or in c:\\temp\\claudetemp. Without one
+     this module says so and the pipeline carries on unsolved rather than
+     failing.
+   * **It uploads the image.** A stacked master goes to a third-party service.
+     Fine for these targets, and worth knowing before pointing it at anything
+     you would not publish. This is the main reason the estate solver is tried
+     first.
+
+Whichever source answers, the solution is put through the same quality gate as
+the seeded solver - verified against Gaia, with a star count and residual it
+has to meet. A solution has to be good, not merely returned. The source is
+recorded in the ASTRSOLV card and in the returned `method`, so it is always
+possible to tell afterwards whether a frame was solved in-house or uploaded.
+"""
+import os
+
+import numpy as np
+from astropy.io import fits
+from astropy.wcs import WCS
+
+import layout
+
+KEY_FILES = ("astrometry.key",
+             r"c:\temp\claudetemp\astrometry.key")
+MIN_STARS = 40
+MAX_RESID_PX = 1.5
+
+
+def find_key(session=None):
+    """The API key, from the environment or a key file."""
+    key = os.environ.get("ASTROMETRY_API_KEY")
+    if key:
+        return key.strip()
+    candidates = list(KEY_FILES)
+    if session is not None:
+        candidates.insert(0, os.path.join(session.root, "astrometry.key"))
+    for path in candidates:
+        if os.path.exists(path):
+            with open(path, encoding="utf-8") as fh:
+                key = fh.read().strip()
+            if key:
+                return key
+    return None
+
+
+def _solve_via_estate(path, scale_hint=None, verbose=True):
+    """Try the estate's own solver first. Returns a WCS header, or None.
+
+    Preferred over nova for three reasons, in order of how much they matter:
+    the image never leaves the house, it needs no API key, and it answers in
+    seconds rather than sitting in a public queue. If the service is down or
+    unreachable this returns None and the caller goes to nova, so the estate
+    being offline costs speed and privacy, not the ability to solve.
+    """
+    try:
+        import astrometry_net
+    except ImportError:
+        return None
+    if not astrometry_net.available():
+        if verbose:
+            print(f"        estate solver unreachable "
+                  f"({astrometry_net.SERVICE}) - falling back to nova")
+        return None
+
+    wcs, info = astrometry_net.solve(path, scale=scale_hint, verbose=verbose)
+    if wcs is None:
+        # `info` is the reason string here. A genuine "no solution" is worth
+        # reporting but is not a reason to then try nova with the same pixels
+        # and the same index scales - nova's indexes are wider, so it is worth
+        # exactly one more attempt, and the caller makes it.
+        if verbose:
+            print(f"        estate solver: {info}")
+        return None
+    return wcs.to_header()
+
+
+def _solve_via_nova(path, session, scale_hint=None, verbose=True):
+    """The hosted service. Needs an API key and uploads the image."""
+    key = find_key(session)
+    if not key:
+        if verbose:
+            print("        nova blind solve unavailable: no astrometry.net API"
+                  " key.\n          Get one free at https://nova.astrometry.net"
+                  " (Profile -> API Key), then set ASTROMETRY_API_KEY or write"
+                  " it to c:\\temp\\claudetemp\\astrometry.key")
+        return None
+
+    from astroquery.astrometry_net import AstrometryNet
+
+    ast = AstrometryNet()
+    ast.api_key = key
+    ast.TIMEOUT = 600
+
+    # Give it the scale if it is known: a hint turns a search over every
+    # possible scale into a search over one, which is the difference between
+    # a minute and a quarter of an hour.
+    kwargs = dict(solve_timeout=600, publicly_visible="n")
+    if scale_hint:
+        kwargs.update(scale_units="arcsecperpix", scale_type="ev",
+                      scale_est=float(scale_hint), scale_err=20)
+    if verbose:
+        print(f"        blind solving {os.path.basename(path)} via "
+              f"nova.astrometry.net"
+              + (f" (scale hint {scale_hint:.3f} arcsec/px)"
+                 if scale_hint else ""))
+
+    try:
+        wcs_header = ast.solve_from_image(path, force_image_upload=True,
+                                          **kwargs)
+    except Exception as exc:                                  # noqa: BLE001
+        print(f"        nova blind solve failed: {type(exc).__name__}: {exc}")
+        return None
+    if not wcs_header:
+        print("        nova blind solve returned no solution")
+        return None
+    return wcs_header
+
+
+def run(session, scale_hint=None, verbose=True):
+    """Solve the deepest master blind. Returns the same dict as the seeded
+    solver, or None.
+
+    Two possible sources, tried in that order: the estate's own solver, then
+    nova. Whichever answers, the solution then earns its place the same way -
+    verified against Gaia and put through the same gate as the seeded solver.
+    A blind solve that is merely *returned* is not a blind solve that is right.
+    """
+    masters = os.path.join(session.root, layout.MASTERS)
+    deepest = session.filters[0]
+    path = os.path.join(masters, f"master-{deepest}.fit")
+    if not os.path.exists(path):
+        return None
+
+    with fits.open(path) as hd:
+        image = hd[0].data.astype(np.float32)
+    ny, nx = image.shape
+
+    source = "estate"
+    wcs_header = _solve_via_estate(path, scale_hint=scale_hint, verbose=verbose)
+    if wcs_header is None:
+        source = "nova"
+        wcs_header = _solve_via_nova(path, session, scale_hint=scale_hint,
+                                     verbose=verbose)
+    if wcs_header is None:
+        return None
+
+    wcs = WCS(wcs_header)
+    cd = wcs.pixel_scale_matrix * 3600.0
+    solved_scale = float(np.sqrt(abs(np.linalg.det(cd))))
+    rot = float(np.degrees(np.arctan2(cd[0, 1], cd[1, 1])))
+    centre = wcs.pixel_to_world(nx / 2.0, ny / 2.0)
+
+    # Verify against Gaia rather than trusting the service, and put the answer
+    # through the same gate as the seeded solver.
+    import astrometry as seeded
+    from astropy import units as u
+    from astropy.coordinates import SkyCoord
+    from measure import measure_frame
+
+    det = measure_frame(path, max_stars=1500)
+    src = det["xy"]
+    radius = 0.6 * np.hypot(nx, ny) * solved_scale / 3600.0
+    cache = os.path.join(session.root, layout.INTERMEDIATES, "_gaia_blind.npz")
+    ra, dec, _ = seeded.fetch_catalogue(centre, radius, cache, limit=20000,
+                                        deep=True)
+    nmatch, resid = 0, None
+    if ra is not None and len(src) > 10:
+        world = SkyCoord(ra * u.deg, dec * u.deg)
+        px, py = wcs.world_to_pixel(world)
+        pairs = []
+        for i, (cx, cy) in enumerate(zip(px, py)):
+            if not (0 <= cx < nx and 0 <= cy < ny):
+                continue
+            d = np.hypot(src[:, 0] - cx, src[:, 1] - cy)
+            j = int(np.argmin(d))
+            if d[j] < 3.0:
+                pairs.append((i, j))
+        if pairs:
+            ci = np.array([p[0] for p in pairs])
+            ii = np.array([p[1] for p in pairs])
+            qx, qy = wcs.world_to_pixel(world[ci])
+            resid = np.hypot(qx - src[ii, 0], qy - src[ii, 1])
+            nmatch = len(pairs)
+
+    med = float(np.median(resid)) if resid is not None else float("nan")
+    if nmatch < MIN_STARS or not np.isfinite(med) or med > MAX_RESID_PX:
+        print(f"        {source} blind solution REJECTED on verification: "
+              f"{nmatch} Gaia stars at {med:.2f} px "
+              f"(need >= {MIN_STARS} and <= {MAX_RESID_PX})")
+        return None
+
+    if verbose:
+        print(f"        {source} blind solve verified: {nmatch} Gaia stars, "
+              f"{med:.2f} px ({med * solved_scale:.2f} arcsec)")
+        print(f"        centre {centre.to_string('hmsdms')}, "
+              f"{solved_scale:.4f} arcsec/px, PA {rot:.2f} deg")
+
+    for filt in session.filters:
+        p = os.path.join(masters, f"master-{filt}.fit")
+        if not os.path.exists(p):
+            continue
+        with fits.open(p, mode="update") as hd:
+            for card in wcs.to_header().cards:
+                hd[0].header[card.keyword] = (card.value, card.comment)
+            # Name the SOURCE, not just "astrometry.net": a year from now the
+            # only way to tell whether a frame was solved in-house or uploaded
+            # to a public service is this card.
+            hd[0].header["ASTRSOLV"] = (
+                f"astrometry.net blind ({source}) / {nmatch} Gaia stars / "
+                f"{med * solved_scale:.2f} arcsec", "blind plate solve")
+            hd.flush()
+
+    return dict(nstars=nmatch, residual_px=med,
+                residual_arcsec=med * solved_scale, scale=solved_scale,
+                position_angle=rot, centre=centre.to_string("hmsdms"),
+                fov_arcmin=(nx * solved_scale / 60, ny * solved_scale / 60),
+                method=f"astrometry.net blind ({source})")

pipeline/colour.py

@@ -0,0 +1,246 @@
+"""Stage 5: assemble a viewable image from whatever filters the session has.
+
+The palette is chosen from what is present, not assumed:
+
+    LRGB   colour from R/G/B, brightness from the deep luminance
+    RGB    no luminance, so a synthetic one is made from the colour channels
+    SHO    the Hubble palette: SII -> red, Ha -> green, OIII -> blue
+    HOO    Ha -> red, OIII -> green and blue
+    MONO   one filter, rendered as greyscale
+
+Background handling deserves a note, because the obvious approach is wrong in a
+way that is invisible. Fitting a plane to a frame that a large target fills
+makes the plane absorb the target's own outer light - on Centaurus A this was
+measured at -17.9 ADU/px of real halo quietly subtracted away. So the fit
+excludes a central region whose size is derived from where the signal actually
+is, and only a plane is used, never a flexible surface.
+"""
+import os
+
+import numpy as np
+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
+
+BG_TARGET = 0.10          # where the sky sits in the stretched image
+SATURATION = 1.35
+
+PALETTE_MAP = {
+    "SHO": {"Red": "SII", "Green": "Ha", "Blue": "OIII"},
+    "HOO": {"Red": "Ha", "Green": "OIII", "Blue": "OIII"},
+}
+
+
+def mtf(x, midtone):
+    """Midtone transfer function on data already scaled to [0, 1]."""
+    x = np.clip(x, 0.0, 1.0)
+    return ((midtone - 1.0) * x) / ((2.0 * midtone - 1.0) * x - midtone)
+
+
+def autostretch(img, target=BG_TARGET, shadow_sigma=2.8):
+    """Black point just below sky, then an MTF putting sky at `target`."""
+    sky = float(np.median(img))
+    mad = 1.4826 * float(np.median(np.abs(img - sky)))
+    if mad <= 0:
+        mad = max(float(np.std(img)), 1e-6)
+    black = sky - shadow_sigma * mad
+    white = float(np.percentile(img, 99.995))
+    if white <= black:
+        white = black + 1.0
+    norm = np.clip((img - black) / (white - black), 0.0, 1.0)
+    sky_norm = (sky - black) / (white - black)
+    m = ((target - 1.0) * sky_norm) / (2.0 * target * sky_norm - target -
+                                       sky_norm)
+    m = float(np.clip(m, 1e-4, 0.9))
+    return mtf(norm, m), dict(black=black, white=white, midtone=m, sky=sky,
+                              mad=mad)
+
+
+def target_mask(img, frac=0.45):
+    """Ellipse covering the bright central object, sized from the data.
+
+    A fixed radius cannot suit both a galaxy filling the frame and a small
+    nebula. This grows the ellipse until it contains most of the flux above
+    sky, then stops - so the background fit is excluded from wherever the
+    target actually is, whatever its size.
+    """
+    ny, nx = img.shape
+    yy, xx = np.mgrid[0:ny, 0:nx]
+    cy, cx = ny / 2.0, nx / 2.0
+    r = np.hypot((xx - cx) / (nx / 2.0), (yy - cy) / (ny / 2.0))
+    sky = float(np.median(img))
+    signal = np.clip(img - sky, 0, None)
+    total = signal.sum()
+    if total <= 0:
+        return r < 0.5
+    for cut in np.arange(0.25, 0.96, 0.05):
+        if signal[r < cut].sum() >= frac * total:
+            return r < min(cut + 0.15, 0.95)
+    return r < 0.75
+
+
+def remove_gradient(img, mask):
+    """Subtract a plane fitted to tile medians outside `mask`."""
+    ny, nx = img.shape
+    step = max(48, min(ny, nx) // 40)
+    xs, ys, zs = [], [], []
+    for y0 in range(0, ny - step, step):
+        for x0 in range(0, nx - step, step):
+            if mask[y0:y0 + step, x0:x0 + step].any():
+                continue
+            tile = img[y0:y0 + step, x0:x0 + step]
+            zs.append(np.median(tile))
+            xs.append(x0 + step / 2.0)
+            ys.append(y0 + step / 2.0)
+    if len(zs) < 12:
+        return img - float(np.median(img)), None
+    xs, ys, zs = map(np.asarray, (xs, ys, zs))
+    keep = np.ones(len(zs), bool)
+    coef = None
+    for _ in range(3):
+        A = np.column_stack([xs[keep], ys[keep], np.ones(keep.sum())])
+        coef, *_ = np.linalg.lstsq(A, zs[keep], rcond=None)
+        model = coef[0] * xs + coef[1] * ys + coef[2]
+        resid = zs - model
+        s = 1.4826 * np.median(np.abs(resid - np.median(resid)))
+        if s <= 0:
+            break
+        keep = np.abs(resid - np.median(resid)) < 2.5 * s
+    yy, xx = np.mgrid[0:ny, 0:nx]
+    plane = (coef[0] * xx + coef[1] * yy + coef[2]).astype(np.float32)
+    return img - plane, coef
+
+
+def load_masters(session):
+    masters = {}
+    d = os.path.join(session.root, layout.MASTERS)
+    for name in os.listdir(d) if os.path.isdir(d) else []:
+        if name.startswith("master-") and name.endswith(".fit"):
+            filt = name[len("master-"):-len(".fit")]
+            masters[filt] = os.path.join(d, name)
+    return 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.
+
+    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 gradient:
+        for filt in list(data):
+            data[filt], _ = remove_gradient(data[filt], mask)
+
+    if palette in PALETTE_MAP:
+        mapping = PALETTE_MAP[palette]
+        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)
+
+    refs = [float(np.percentile(c, 99.5)) for c in channels]
+    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)
+
+    unlinked = palette in PALETTE_MAP
+    rgb = np.empty(lum.shape + (3,), np.float32)
+    for i, ch in enumerate(channels):
+        if unlinked:
+            rgb[:, :, i], _ = autostretch(ch)
+        else:
+            sky = float(np.median(ch))
+            mad = 1.4826 * float(np.median(np.abs(ch - sky))) or 1.0
+            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),
+                params["midtone"])
+    if unlinked:
+        lum = rgb.mean(axis=2)
+
+    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)")
+
+    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)
+
+    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 out
+
+
+def run(session, verbose=True):
+    """Produce the four deliverable images for this session."""
+    import finals
+    return finals.build(session, verbose=verbose)

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)}")

pipeline/measure.py

@@ -0,0 +1,140 @@
+"""Stage 2: measure every frame and cache its star list.
+
+Two jobs. It produces the quality numbers used to choose a registration
+reference and to weight the stack - sky level, noise, star sharpness - and it
+caches each frame's detected stars, because both registration and the plate
+solve need them and detection costs far more than reading a small array back.
+
+Frames are opened one at a time. A 4096x4096 float32 frame is 67 MB and a
+session can hold ninety of them; there is no reason to have more than one in
+memory at once.
+"""
+import os
+
+import numpy as np
+import sep
+
+from astropy.io import fits
+
+import layout
+
+MAX_STARS = 500          # cached per frame, brightest first
+DETECT_SIGMA = 5.0
+
+
+def measure_frame(path, max_stars=MAX_STARS):
+    """Sky, noise, seeing and a star list for one frame."""
+    with fits.open(path, memmap=False) as hd:
+        data = hd[0].data.astype(np.float32)
+
+    bkg = sep.Background(data, bw=64, bh=64, fw=3, fh=3)
+    sky = float(np.median(bkg.back()))
+    rms = float(bkg.globalrms)
+    sub = data - bkg.back()
+
+    objs = sep.extract(sub, DETECT_SIGMA, err=rms, minarea=5,
+                       deblend_cont=0.005)
+    # Saturated cores and cosmic-ray hits make poor registration anchors and
+    # poor seeing estimates, so they go first.
+    objs = objs[(objs["flag"] == 0) & (objs["npix"] < 3000)]
+
+    # The minimum size has to follow the SAMPLING, not be a fixed number of
+    # pixels. A fixed floor of 12 px assumes a well sampled star: at 0.53
+    # arcsec/px with 5.9 px seeing that is right, but at 3.5 arcsec/px the
+    # stars are undersampled at 1.75 px FWHM and cover only a handful of
+    # pixels each, so the same floor discards nearly every real star and keeps
+    # blends and galaxies instead. That is what made the wide-field sessions
+    # impossible to plate solve: the surviving detections were not the objects
+    # the catalogue lists.
+    if len(objs) > 10:
+        rad0, _ = sep.flux_radius(sub, objs["x"], objs["y"], 6.0 * objs["a"],
+                                  0.5, normflux=objs["flux"])
+        fwhm0 = float(np.median(rad0) * 2.0)
+        min_npix = max(4, int(0.5 * np.pi * (fwhm0 / 2.0) ** 2))
+        objs = objs[objs["npix"] >= min_npix]
+    if len(objs) == 0:
+        del data, sub, bkg
+        return dict(sky=sky, rms=rms, fwhm=np.nan, nstars=0,
+                    xy=np.zeros((0, 2)), flux=np.zeros(0))
+
+    order = np.argsort(objs["flux"])[::-1][:max_stars]
+    objs = objs[order]
+    radius, _ = sep.flux_radius(sub, objs["x"], objs["y"], 6.0 * objs["a"],
+                                0.5, normflux=objs["flux"])
+    fwhm = float(np.median(radius) * 2.0)
+
+    result = dict(
+        sky=sky, rms=rms, fwhm=fwhm, nstars=int(len(objs)),
+        xy=np.column_stack([objs["x"], objs["y"]]).astype(np.float64),
+        flux=objs["flux"].astype(np.float64),
+    )
+    del data, sub, bkg, objs
+    return result
+
+
+def run(session, force=False, verbose=True):
+    """Measure every frame in a session, caching to intermediates/."""
+    cache_path = os.path.join(session.root, layout.INTERMEDIATES,
+                              "_measure.npz")
+    os.makedirs(os.path.dirname(cache_path), exist_ok=True)
+
+    cached = {}
+    if os.path.exists(cache_path) and not force:
+        z = np.load(cache_path, allow_pickle=True)
+        cached = {k: z[k] for k in z.files}
+
+    store, rows = {}, []
+    for frame in session.frames:
+        key = frame.name
+        if f"{key}|xy" in cached and not force:
+            xy = cached[f"{key}|xy"]
+            stats = cached[f"{key}|stats"].item()
+        else:
+            stats = measure_frame(frame.path)
+            xy = stats.pop("xy")
+            flux = stats.pop("flux")
+            store[f"{key}|flux"] = flux
+        store[f"{key}|xy"] = xy
+        store[f"{key}|stats"] = np.array(stats, dtype=object)
+        rows.append((frame, stats))
+        if verbose:
+            print(f"      {frame.filter:10s} {frame.name[-28:]:28s} "
+                  f"sky={stats['sky']:8.1f} rms={stats['rms']:7.2f} "
+                  f"fwhm={stats['fwhm']:5.2f}px stars={stats['nstars']:4d}")
+
+    # Carry forward anything already cached that was not re-measured.
+    for k, v in cached.items():
+        store.setdefault(k, v)
+    np.savez_compressed(cache_path, **store)
+    return rows
+
+
+def load(session):
+    """Star lists and stats for a session, as measured earlier."""
+    path = os.path.join(session.root, layout.INTERMEDIATES, "_measure.npz")
+    if not os.path.exists(path):
+        raise FileNotFoundError(
+            f"{path} not found - run the measure stage first")
+    z = np.load(path, allow_pickle=True)
+    xy = {k.split("|")[0]: z[k] for k in z.files if k.endswith("|xy")}
+    stats = {k.split("|")[0]: z[k].item() for k in z.files
+             if k.endswith("|stats")}
+    return xy, stats
+
+
+def choose_reference(session, stats):
+    """Pick the frame everything else is aligned to.
+
+    The reference sets the output grid, so it should be the sharpest frame of
+    the filter with the most signal - which is also usually the filter with the
+    finest sampling. Choosing a poor reference costs resolution in every other
+    frame, permanently, because registration can only resample onto it.
+    """
+    best_filter = session.filters[0]
+    candidates = [f for f in session.by_filter(best_filter)
+                  if stats.get(f.name, {}).get("nstars", 0) >= 8]
+    if not candidates:
+        candidates = session.by_filter(best_filter)
+    ref = min(candidates,
+              key=lambda f: stats.get(f.name, {}).get("fwhm", np.inf))
+    return ref

pipeline/register.py

@@ -0,0 +1,159 @@
+"""Stage 3: register every frame onto one grid and combine per filter.
+
+Registration matches asterisms between the cached star lists rather than
+cross-correlating pixels: matching a few hundred coordinates is far cheaper
+than comparing 15 megapixel frames, and it copes with the field rotating
+between the east and west sides of the meridian.
+
+**One reference for all filters, not one per filter.** That is what makes the
+masters pixel-aligned, so the colour composite needs no further registration.
+
+The awkward case this has to survive: NGC 6744 was shot with luminance at bin1
+(4096x4096) and colour at bin2 (2048x2048). Frames do not share a shape or a
+pixel scale. Because astroalign returns a similarity transform - rotation,
+translation AND scale - the maths already handles it; what does not handle it
+is assuming the output is the same shape as the input. Every warp is therefore
+given the reference's shape explicitly, and the bin2 frames are resampled up
+onto the bin1 grid.
+"""
+import os
+
+import astroalign as aa
+import numpy as np
+from astropy.io import fits
+from astropy.stats import sigma_clip
+from skimage.transform import warp
+
+import layout
+import measure as measure_mod
+
+aa.MIN_MATCHES_FRACTION = 0.6
+aa.NUM_NEAREST_NEIGHBORS = 8
+
+
+def combine(cube, weights, sigma=3.0):
+    """Weighted sigma-clipped mean, in row blocks to bound peak memory.
+
+    sigma_clip allocates a mask and float64 intermediates; on a full cube that
+    can triple peak usage, which matters on a machine with a few GB free.
+    """
+    n, ny, nx = cube.shape
+    out = np.zeros((ny, nx), np.float32)
+    w = np.asarray(weights, np.float32)[:, None, None]
+    step = 256
+    for y0 in range(0, ny, step):
+        y1 = min(y0 + step, ny)
+        block = cube[:, y0:y1, :]
+        if n >= 3:
+            clipped = sigma_clip(block, sigma=sigma, maxiters=1, axis=0,
+                                 masked=True, copy=True)
+            good = ~clipped.mask
+        else:
+            # With one or two frames there is nothing to reject against, and
+            # clipping would just throw away signal.
+            good = np.ones(block.shape, bool)
+        finite = np.isfinite(block)
+        good &= finite
+        vals = np.where(finite, block, 0.0)
+        wb = np.broadcast_to(w, block.shape) * good
+        denom = wb.sum(axis=0)
+        denom[denom == 0] = np.nan
+        out[y0:y1, :] = np.nansum(vals * wb, axis=0) / denom
+        del block, finite, good, vals, wb, denom
+    return np.nan_to_num(out, nan=0.0)
+
+
+def run(session, verbose=True):
+    """Register and stack every filter; write masters to stacks/masters/."""
+    xy, stats = measure_mod.load(session)
+    ref = measure_mod.choose_reference(session, stats)
+    ref_xy = xy[ref.name]
+    ref_shape = ref.shape
+    if verbose:
+        print(f"      reference {ref.name[-32:]} ({ref.filter}, "
+              f"fwhm {stats[ref.name]['fwhm']:.2f} px, "
+              f"{ref_shape[1]}x{ref_shape[0]})")
+
+    out_dir = os.path.join(session.root, layout.MASTERS)
+    os.makedirs(out_dir, exist_ok=True)
+    written = {}
+
+    for filt in session.filters:
+        frames = session.by_filter(filt)
+        planes, weights, used = [], [], []
+        header0 = None
+        for frame in frames:
+            with fits.open(frame.path, memmap=False) as hd:
+                data = hd[0].data.astype(np.float32)
+                if header0 is None:
+                    header0 = hd[0].header.copy()
+            if frame.name == ref.name:
+                reg = data
+            else:
+                src = xy.get(frame.name)
+                if src is None or len(src) < 8:
+                    print(f"        {frame.name[-28:]}: too few stars - "
+                          f"dropped")
+                    del data
+                    continue
+                try:
+                    tform, _ = aa.find_transform(src, ref_xy)
+                except Exception as exc:                      # noqa: BLE001
+                    print(f"        {frame.name[-28:]}: registration failed "
+                          f"({type(exc).__name__}) - dropped")
+                    del data
+                    continue
+                # output_shape is the whole point: without it a bin2 frame
+                # would be written onto a bin2-sized grid and silently fail to
+                # line up with a bin1 reference.
+                reg = warp(data, inverse_map=tform.inverse,
+                           output_shape=ref_shape, order=3, mode="constant",
+                           cval=np.nan, preserve_range=True).astype(np.float32)
+                del data
+            sky = float(np.nanmedian(reg))
+            reg -= sky
+            planes.append(reg)
+            rms = stats.get(frame.name, {}).get("rms", 1.0) or 1.0
+            weights.append(1.0 / (rms * rms))
+            used.append(frame)
+
+        if not planes:
+            print(f"      {filt}: no frames survived registration - skipped")
+            continue
+
+        cube = np.stack(planes, axis=0)
+        del planes
+        # Normalise for transparency: scale each frame so its bright signal
+        # matches the group, so a frame through thin cloud cannot drag the
+        # mean down.
+        levels = np.array([np.nanpercentile(p, 99.5) for p in cube])
+        ref_level = np.nanmedian(levels)
+        for i, lv in enumerate(levels):
+            if lv > 0:
+                cube[i] *= float(ref_level / lv)
+
+        master = combine(cube, weights)
+        nframes = cube.shape[0]
+        total_exp = sum(f.exptime for f in used)
+        del cube
+
+        hdr = header0
+        for key in ("CBLACK", "CWHITE", "PEDESTAL", "HISTORY"):
+            hdr.remove(key, ignore_missing=True, remove_all=True)
+        hdr["FILTER"] = filt
+        hdr["NCOMBINE"] = (nframes, "frames in this master")
+        hdr["EXPTOTAL"] = (total_exp, "[s] total integration")
+        hdr["STACKREF"] = (ref.name[:60], "registration reference frame")
+        hdr["STACKALG"] = ("sigma-clipped weighted mean" if nframes >= 3
+                           else "weighted mean (too few frames to clip)",
+                           "combine method")
+        hdr["IMAGETYP"] = "Master Light"
+        path = os.path.join(out_dir, f"master-{filt}.fit")
+        fits.PrimaryHDU(master.astype(np.float32), hdr).writeto(
+            path, overwrite=True)
+        written[filt] = path
+        if verbose:
+            print(f"      {filt:10s} {nframes:3d} frames  "
+                  f"{total_exp/60:6.1f} min  -> master-{filt}.fit")
+        del master
+    return written, ref

pipeline/report.py

@@ -0,0 +1,188 @@
+"""Stage 6: write the session's own METHODS.md.
+
+Every session directory gets an account of what was done to it and what came
+out, written for someone who was not there and has no access to the chat that
+produced it. That is the whole point of generating it rather than writing it by
+hand: a hand-written note is written once and then rots, while this is rebuilt
+from the actual results every time the pipeline runs.
+
+It records what was measured, not what was intended - including the things that
+did not work, because a session where the plate solve failed is more useful
+with that stated than with it omitted.
+"""
+import os
+from datetime import datetime, timezone
+
+import numpy as np
+
+import layout
+
+
+def _fmt_exposure(seconds):
+    return f"{seconds / 60:.0f} min" if seconds >= 60 else f"{seconds:.0f} s"
+
+
+def write(session, results, verbose=True):
+    """Compose METHODS.md from whatever the stages actually returned."""
+    stats = results.get("stats") or {}
+    solve = results.get("solve")
+    reference = results.get("reference")
+    palette = session.palette
+    now = datetime.now(timezone.utc).strftime("%Y-%m-%d %H:%M UTC")
+
+    lines = []
+    add = lines.append
+
+    add(f"# {session.target} - {session.telescope}")
+    add("")
+    add(f"Processed {now} by the pipeline in the "
+        "[astrophotography](https://git.discworld.casa/laurence/astrophotography) "
+        "repository (`pipeline/`). This file is generated: re-running the "
+        "pipeline rewrites it from the actual results.")
+    add("")
+
+    # ---------------------------------------------------------- what is here
+    add("## What was captured")
+    add("")
+    add("| Filter | Frames | Exposure | Binning | Total |")
+    add("|---|---|---|---|---|")
+    total = 0.0
+    for filt in session.filters:
+        frames = session.by_filter(filt)
+        exps = sorted({f.exptime for f in frames})
+        bins = sorted({f.binning for f in frames})
+        secs = sum(f.exptime for f in frames)
+        total += secs
+        add(f"| {filt} | {len(frames)} | "
+            f"{'/'.join(f'{e:.0f} s' for e in exps)} | "
+            f"{'/'.join(f'bin{b}' for b in bins)} | {_fmt_exposure(secs)} |")
+    add(f"| **Total** | **{len(session.frames)}** | | | "
+        f"**{_fmt_exposure(total)}** |")
+    add("")
+    add(f"Colour scheme: **{palette}**"
+        + {"SHO": " - the Hubble palette, SII to red, Ha to green, OIII to blue",
+           "HOO": " - Ha to red, OIII to green and blue",
+           "LRGB": " - colour from R/G/B, detail from the deeper luminance",
+           "RGB": " - no luminance channel, so brightness is synthesised from "
+                  "the colour channels",
+           "MONO": " - a single filter, rendered as greyscale"}.get(palette, ""))
+    add("")
+
+    # -------------------------------------------------------- frame quality
+    if stats:
+        add("## Frame quality")
+        add("")
+        add("| Filter | Sky (ADU) | Noise (ADU) | Seeing (px) |")
+        add("|---|---|---|---|")
+        for filt in session.filters:
+            vals = [stats[f.name] for f in session.by_filter(filt)
+                    if f.name in stats]
+            if not vals:
+                continue
+            sky = np.median([v["sky"] for v in vals])
+            rms = np.median([v["rms"] for v in vals])
+            fwhm = np.median([v["fwhm"] for v in vals])
+            add(f"| {filt} | {sky:.1f} | {rms:.2f} | {fwhm:.2f} |")
+        add("")
+
+    # ---------------------------------------------------------- what was done
+    add("## What was done")
+    add("")
+    add("1. **Ingest.** Archives uncompressed into `calibrated/`. Only the "
+        "`calibrated-` set is extracted; the `raw-` archives stay zipped in "
+        "`raw/` as untouched originals, because extracting both would put two "
+        "copies of every frame into the stack.")
+    add("2. **Measure.** Sky, noise and star sharpness per frame, with each "
+        "frame's star list cached for the stages that follow.")
+    if reference:
+        add(f"3. **Register and stack.** Every frame aligned by matching star "
+            f"patterns onto one reference - `{reference.name}` "
+            f"({reference.filter}, the sharpest frame of the deepest filter) - "
+            f"so all masters share a pixel grid and the colour image needs no "
+            f"further alignment. Combined with a "
+            f"{'sigma-clipped ' if len(session.frames) >= 3 else ''}"
+            f"noise-weighted mean.")
+    else:
+        add("3. **Register and stack.**")
+    if solve:
+        add(f"4. **Plate solve.** Matched against Gaia DR3: "
+            f"**{solve['nstars']} stars, {solve['residual_arcsec']:.2f} arcsec "
+            f"residual**. The WCS is written into every master, so any position "
+            f"measured from these images is quotable.")
+    else:
+        add("4. **Plate solve. FAILED for this session** - the star pattern "
+            "match did not converge against Gaia. The images are unaffected, "
+            "but nothing here carries sky coordinates, so positions cannot be "
+            "measured from them.")
+    add(f"5. **Colour.** Assembled as {palette}, with the sky background "
+        "levelled by a plane fitted outside the target and the three channels' "
+        "backgrounds forced neutral.")
+    add("")
+
+    # ------------------------------------------------------------ astrometry
+    if solve:
+        add("## Astrometry")
+        add("")
+        add(f"- Field centre: **{solve['centre']}**")
+        add(f"- Plate scale: **{solve['scale']:.4f} arcsec/pixel**")
+        add(f"- Field of view: {solve['fov_arcmin'][0]:.1f}' x "
+            f"{solve['fov_arcmin'][1]:.1f}'")
+        add(f"- Position angle: {solve['position_angle']:.2f} deg")
+        add(f"- Residual: {solve['residual_px']:.2f} px "
+            f"({solve['residual_arcsec']:.2f} arcsec) on {solve['nstars']} stars")
+        add("")
+
+    # ----------------------------------------------------------- honest limits
+    add("## What limits this data")
+    add("")
+    shortest = min((sum(f.exptime for f in session.by_filter(x))
+                    for x in session.filters), default=0)
+    if len(session.frames) < 6:
+        add(f"- **Very few frames ({len(session.frames)} in total).** With so "
+            "few there is no meaningful outlier rejection, so satellite "
+            "trails and cosmic rays survive into the masters.")
+    if shortest < 600:
+        add(f"- **Short integration** - the thinnest channel has only "
+            f"{_fmt_exposure(shortest)}. No processing recovers signal that "
+            "was never collected.")
+    if palette == "RGB":
+        add("- **No luminance channel**, so brightness is synthesised from the "
+            "colour data and the result is noisier than an LRGB set of the "
+            "same total exposure would be.")
+    if palette == "MONO":
+        add("- **A single filter**, so there is no colour information at all.")
+    if not solve:
+        add("- **No astrometric solution**, so nothing here can be used for "
+            "positions or cross-matched against a catalogue.")
+    fwhms = [v["fwhm"] for v in stats.values() if np.isfinite(v.get("fwhm", np.nan))]
+    if fwhms and solve:
+        seeing = float(np.median(fwhms)) * solve["scale"]
+        add(f"- **Seeing of about {seeing:.1f} arcsec**, which sets the finest "
+            "detail present regardless of how the data is processed.")
+    add("")
+
+    # ------------------------------------------------------------- directories
+    add("## Where things are")
+    add("")
+    add("| Directory | Contents |")
+    add("|---|---|")
+    add("| `raw/` | exactly what the telescope delivered; untouched |")
+    add("| `calibrated/` | the uncompressed calibrated frames |")
+    add("| `stacks/masters/` | one registered, combined master per filter |")
+    add("| `final/` | the finished image |")
+    add("| `intermediates/` | caches; deletable, regenerated by a re-run |")
+    add("")
+    add("Re-run any stage by pointing the pipeline at this directory:")
+    add("")
+    add("```")
+    add(f"set ASTRO_SESSION={session.root}")
+    add("python pipeline/run.py")
+    add("```")
+    add("")
+
+    path = os.path.join(session.root, "METHODS.md")
+    with open(path, "w", encoding="utf-8") as fh:
+        fh.write("\n".join(lines))
+    if verbose:
+        print(f"      -> METHODS.md")
+    return path

pipeline/run.py

@@ -0,0 +1,181 @@
+"""The push-button entry point: process a session, or a whole archive.
+
+    python run.py                          # process $ASTRO_SESSION
+    python run.py PATH                     # process one session
+    python run.py --all PATH               # find and process every session
+    python run.py --stage register PATH    # run one stage
+    python run.py --force PATH             # ignore existing outputs
+
+Stages run in order and each is skipped when its output already exists, so an
+interrupted run continues rather than starting again. A failure in one stage of
+one session does not stop the others: the failure is recorded and the batch
+carries on, because in a batch of twenty the useful outcome is nineteen results
+and one clear error, not nothing.
+"""
+import argparse
+import os
+import sys
+import time
+import traceback
+
+import layout
+import session as session_mod
+
+STAGES = ("ingest", "measure", "register", "solve", "colour",
+          "science", "report")
+
+
+def _exists(path):
+    return os.path.exists(path) and os.path.getsize(path) > 0
+
+
+def process(root, stages=STAGES, force=False, verbose=True):
+    """Run the requested stages for one session. Returns a result summary."""
+    s = session_mod.Session(root)
+    result = {"root": root, "stages": {}, "ok": True}
+    t0 = time.time()
+
+    if "ingest" in stages:
+        s.ingest(verbose=verbose)
+    s.scan(verbose=verbose)
+    if not s.frames:
+        print("    no light frames found - skipping")
+        result["ok"] = False
+        result["error"] = "no frames"
+        return result
+    s.save_manifest()
+    result["target"] = s.target
+    result["telescope"] = s.telescope
+    result["palette"] = s.palette
+    result["frames"] = len(s.frames)
+
+    results = {}
+    try:
+        if "measure" in stages:
+            import measure
+            cache = os.path.join(root, layout.INTERMEDIATES, "_measure.npz")
+            if force or not _exists(cache):
+                print("    measure")
+                measure.run(s, force=force, verbose=False)
+            else:
+                print("    measure (cached)")
+            _, results["stats"] = measure.load(s)
+
+        if "register" in stages:
+            import measure
+            import register
+            masters = os.path.join(root, layout.MASTERS)
+            have = _exists(os.path.join(masters, f"master-{s.filters[0]}.fit"))
+            if force or not have:
+                print("    register")
+                _, ref = register.run(s, verbose=verbose)
+            else:
+                print("    register (cached)")
+                _, st = measure.load(s)
+                ref = measure.choose_reference(s, st)
+            results["reference"] = ref
+
+        if "solve" in stages:
+            import astrometry
+            print("    solve")
+            try:
+                results["solve"] = astrometry.run(s, verbose=verbose)
+            except Exception as exc:                          # noqa: BLE001
+                # An unsolved session still produces images; record and go on.
+                print(f"      solve failed: {type(exc).__name__}: {exc}")
+                results["solve"] = None
+            if not results.get("solve"):
+                # The seeded solver needs the image and the catalogue to hold
+                # recognisably the same stars. When they do not, a blind solve
+                # ignores the pointing entirely and works from the pixels.
+                import blind
+                try:
+                    results["solve"] = blind.run(
+                        s, scale_hint=s.scale, verbose=verbose)
+                except Exception as exc:                      # noqa: BLE001
+                    print(f"      blind solve failed: "
+                          f"{type(exc).__name__}: {exc}")
+
+        if "colour" in stages:
+            import colour
+            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")
+            report.write(s, results, verbose=verbose)
+
+    except Exception as exc:                                  # noqa: BLE001
+        result["ok"] = False
+        result["error"] = f"{type(exc).__name__}: {exc}"
+        print(f"    FAILED: {result['error']}")
+        if verbose:
+            traceback.print_exc(limit=3)
+
+    result["solved"] = bool(results.get("solve"))
+    result["seconds"] = time.time() - t0
+    return result
+
+
+def main(argv=None):
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("path", nargs="?", default=None,
+                    help="session directory, or a tree with --all")
+    ap.add_argument("--all", action="store_true",
+                    help="find and process every session beneath PATH")
+    ap.add_argument("--stage", action="append", choices=STAGES,
+                    help="run only this stage (repeatable)")
+    ap.add_argument("--force", action="store_true",
+                    help="re-run stages even if their output exists")
+    args = ap.parse_args(argv)
+
+    root = args.path or layout.SESSION
+    stages = tuple(args.stage) if args.stage else STAGES
+
+    targets = session_mod.discover(root) if args.all else [root]
+    if not targets:
+        print(f"no sessions found under {root}")
+        return 1
+    print(f"{len(targets)} session(s), stages: {', '.join(stages)}\n")
+
+    results = []
+    for i, t in enumerate(targets, 1):
+        print(f"[{i}/{len(targets)}] {os.path.relpath(t, root)}")
+        results.append(process(t, stages=stages, force=args.force))
+        print()
+
+    print("=" * 72)
+    print(f"{'target':16s} {'scope':16s} {'palette':8s} "
+          f"{'frames':>6s} {'solved':>7s} {'time':>7s}")
+    for r in results:
+        if not r.get("ok") and "target" not in r:
+            print(f"{os.path.basename(r['root']):16s} "
+                  f"FAILED: {r.get('error', 'unknown')}")
+            continue
+        print(f"{r.get('target', '?')[:16]:16s} "
+              f"{r.get('telescope', '?')[:16]:16s} "
+              f"{r.get('palette', '?'):8s} {r.get('frames', 0):6d} "
+              f"{'yes' if r.get('solved') else 'no':>7s} "
+              f"{r.get('seconds', 0):6.0f}s")
+    failed = [r for r in results if not r.get("ok")]
+    unsolved = [r for r in results if r.get("ok") and not r.get("solved")]
+    print(f"\n{len(results) - len(failed)}/{len(results)} processed"
+          + (f", {len(unsolved)} without an astrometric solution" if unsolved
+             else ""))
+    return 1 if failed else 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

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")}

state/DECISIONS.md

@@ -3,6 +3,29 @@
 > A dated, append only log of decisions and their rationale. Newest at the top. Never
 > rewrite past entries; if a decision is reversed, add a new entry that says so.
 
+## 2026-07-21: blind solving goes to the estate's own solver first, nova second
+
+The estate now runs Astrometry.net itself
+(`https://astrometry.ankh-morpork.discworld.network`, config in
+ankh-morpork-infra `astrometry/`): solve-field against ~5 GB of local index
+files, sized from this repo's own TELESCOPES.md so it covers the whole
+iTelescope fleet's fields of view (21.5 arcmin to 10.4 degrees).
+
+`blind.py` now tries it before nova.astrometry.net. The order is not about
+speed, though it is far faster - a blind solve of a real DSS field came back in
+0.7s against nova's minutes-in-a-queue. It is that **nova requires uploading
+the image to a third party and an API key**, and neither is now necessary for
+the ordinary case. nova stays as the fallback for when the estate service is
+unreachable or the field is narrower than the local index set reaches, so
+losing the estate costs speed and privacy, not the ability to solve at all.
+
+What deliberately did NOT change: the seeded solver in `astrometry.py` remains
+the primary path (a good header beats a blind search), and every blind solution
+- from either source - still has to pass the same Gaia verification and
+star-count/residual gate. A wrong WCS is worse than no WCS, and a second solver
+does not change that. The source is recorded in the ASTRSOLV card so it is
+always possible to tell later which one produced a frame's astrometry.
+
 ## 2026-07-21: the pipeline reads headers, not filenames
 
 The first pipeline parsed `calibrated-T32-...-Luminance-BIN2-W-300-001.fit` to

state/TODO.md

@@ -18,8 +18,40 @@ one telescope. Branch `generic-pipeline`.
 filter grouping, palette detection (LRGB / RGB / SHO / HOO / MONO), manifest
 writing. Tested against all five archived sessions.
 
-**Next:** `measure` -> `register/stack` -> `solve` -> `compose` -> `report`,
+**Done since:** `measure.py`, `register.py`, `colour.py`, `astrometry.py`.
-then a `run.py` orchestrator with resumable stages.
+Registration verified on the mixed bin1/bin2 case - upsampled bin2 colour
+aligns to the bin1 luminance master within **0.04 px**. Colour renders LRGB,
+RGB, SHO, HOO and mono.
+
+**Next:** finish `astrometry.py` (see below), then `report.py` (per-session
+METHODS.md) and `run.py` (the one-command orchestrator).
+
+**PLATE SOLVE IS 2 OF 4.** NGC 5128 (by hand) and NGC 2030 solve; NGC 2070 and
+NGC 6744 do not. What is already known, so the next session does not repeat it:
+
+- **Not the plate scale.** Header and FOCALLEN/XPIXSZ agree on both failures,
+  and the session that DOES solve has a header/computed discrepancy
+  (0.527 vs 0.686) yet solves anyway - asterism matching is scale invariant,
+  so scale only sets the catalogue cone radius.
+- **Not the catalogue fetch.** Fixed: VizieR is queried first because Gaia's
+  `launch_job_async` submits to a job queue and hung for ten minutes on 600
+  rows while the network was healthy. Gaia synchronous is the fallback.
+- **Not simple truncation.** Fixed: VizieR's `row_limit` truncates BEFORE
+  sorting, so asking for 600 rows returned 600 arbitrary stars rather than the
+  600 brightest. Now fetches up to 50000 and picks the brightest locally.
+- **Partly saturation.** The detector rejects saturated cores, so the brightest
+  DETECTIONS are not the brightest STARS while the catalogue's are, and the two
+  top-N lists barely overlap. NGC 2030 only solved once a sliding magnitude
+  window was tried, and it matched at `catalogue[15:135]` - skipping the 15
+  brightest. That window search is in place and still is not enough for the
+  other two.
+- **Next things to try:** match on a magnitude-matched subset by estimating the
+  image's own zero point first; loosen `aa.MIN_MATCHES_FRACTION`; raise
+  astroalign's `max_control_points`; or fall back to a genuine blind solver
+  (astrometry.net index files locally, or the nova API) for fields that resist
+  the seeded approach. NGC 2070 is a 4 degree field at 3.5 arcsec/px with only
+  305 detections, and NGC 6744 is a rich field - they may need different
+  handling from each other.
 
 **The blocker before any of that:** the disk is full (see below).