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