diff --git a/pipeline/blind.py b/pipeline/blind.py new file mode 100644 index 0000000..26812a2 --- /dev/null +++ b/pipeline/blind.py @@ -0,0 +1,180 @@ +"""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 this uses the hosted nova.astrometry.net service through +astroquery. + +Two things to know: + +* **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 reports that clearly and the + pipeline carries on unsolved rather than failing. +* **It uploads the image.** A stacked master goes to a third-party service. + That is fine for these targets and worth knowing before pointing it at + anything you would not publish. + +Results are still put through the same quality gate as the seeded solver: a +solution has to be good, not merely returned. +""" +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 run(session, scale_hint=None, verbose=True): + """Solve the deepest master blind. Returns the same dict as the seeded + solver, or None.""" + key = find_key(session) + if not key: + if verbose: + print(" 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 + + 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) + header = hd[0].header + ny, nx = image.shape + + 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" blind solve failed: {type(exc).__name__}: {exc}") + return None + if not wcs_header: + print(" blind solve returned no solution") + 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" blind solution REJECTED on verification: {nmatch} " + f"Gaia stars at {med:.2f} px " + f"(need >= {MIN_STARS} and <= {MAX_RESID_PX})") + return None + + if verbose: + print(f" 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) + hd[0].header["ASTRSOLV"] = ( + f"astrometry.net blind / {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="astrometry.net blind") diff --git a/pipeline/run.py b/pipeline/run.py index 3a985c5..ea9610b 100644 --- a/pipeline/run.py +++ b/pipeline/run.py @@ -84,6 +84,17 @@ def process(root, stages=STAGES, force=False, verbose=True): # 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