laurence/astrophotography
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Compare commits

base: laurence:main
Branches Tags
laurence:main
laurence:astrometry-net-fallback
laurence:generic-pipeline
..
compare: laurence:generic-pipeline
Branches Tags
laurence:main
laurence:astrometry-net-fallback
laurence:generic-pipeline

No commits in common. "main" and "generic-pipeline" have entirely different histories.
main ... generic-pi

12 changed files with 2 additions and 2334 deletions
Download patch file Download diff file Expand all files Collapse all files

424
pipeline/astrometry.py
View file

@ -1,424 +0,0 @@
"""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))

142
pipeline/astrometry_net.py
View file

@ -1,142 +0,0 @@
"""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

243
pipeline/blind.py
View file

@ -1,243 +0,0 @@
"""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})")

246
pipeline/colour.py
View file

@ -1,246 +0,0 @@
"""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)

257
pipeline/finals.py
View file

@ -1,257 +0,0 @@
"""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)}")

140
pipeline/measure.py
View file

@ -1,140 +0,0 @@
"""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

159
pipeline/register.py
View file

@ -1,159 +0,0 @@
"""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

188
pipeline/report.py
View file

@ -1,188 +0,0 @@
"""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

181
pipeline/run.py
View file

@ -1,181 +0,0 @@
"""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())

297
pipeline/science.py
View file

@ -1,297 +0,0 @@
"""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")}

23
state/DECISIONS.md
View file

@ -3,29 +3,6 @@
> A dated, append only log of decisions and their rationale. Newest at the top. Never > 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. > 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 ## 2026-07-21: the pipeline reads headers, not filenames
The first pipeline parsed `calibrated-T32-...-Luminance-BIN2-W-300-001.fit` to The first pipeline parsed `calibrated-T32-...-Luminance-BIN2-W-300-001.fit` to

36
state/TODO.md
View file

@ -18,40 +18,8 @@ one telescope. Branch `generic-pipeline`.
filter grouping, palette detection (LRGB / RGB / SHO / HOO / MONO), manifest filter grouping, palette detection (LRGB / RGB / SHO / HOO / MONO), manifest
writing. Tested against all five archived sessions. writing. Tested against all five archived sessions.
**Done since:** `measure.py`, `register.py`, `colour.py`, `astrometry.py`. **Next:** `measure` -> `register/stack` -> `solve` -> `compose` -> `report`,
Registration verified on the mixed bin1/bin2 case - upsampled bin2 colour then a `run.py` orchestrator with resumable stages.
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). **The blocker before any of that:** the disk is full (see below).