From 8bbc5ce3d1c6fce7dc8b8639cd295b6c98b0ae83 Mon Sep 17 00:00:00 2001 From: Caden Gobat <36030084+cgobat@users.noreply.github.com> Date: Tue, 17 Mar 2026 18:24:54 -0700 Subject: [PATCH 1/9] Use astropy modules for FITS I/O and WCS math --- core/postflight/bayer_photometry.py | 80 ++++++----------------------- 1 file changed, 16 insertions(+), 64 deletions(-) diff --git a/core/postflight/bayer_photometry.py b/core/postflight/bayer_photometry.py index d9822ae..d4b1d8a 100644 --- a/core/postflight/bayer_photometry.py +++ b/core/postflight/bayer_photometry.py @@ -22,6 +22,10 @@ from pathlib import Path from typing import Optional, Tuple +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.wcs import WCS + logger = logging.getLogger("seevar.bayer_photometry") # IMX585 Bayer pattern — GRBG @@ -97,86 +101,34 @@ def aperture_flux( return net_flux, sky_median, snr -# --------------------------------------------------------------------------- -# FITS loader — pure numpy, no astropy dependency -# --------------------------------------------------------------------------- - class BayerFITS: """ - Lightweight raw FITS reader. Parses header and pixel array. + Lightweight FITS reader. Parses header and pixel array. Validates Bayer pattern and saturation before science extraction. """ def __init__(self, fits_path: Path): self.path = Path(fits_path) - self.header: dict = {} + self.header: fits.Header = None self.array: Optional[np.ndarray] = None - self._wcs: dict = {} + self._wcs: WCS = None def load(self) -> bool: try: - with open(self.path, "rb") as f: - raw = f.read() - except OSError as e: - logger.error("Cannot open %s: %s", self.path, e) + with fits.open(self.path) as f: + hdu = f[0].copy() + except Exception as e: + logger.error("Failed to read FITS file %s: %s", self.path, e) return False - header = {} - header_blocks = 0 - found_end = False - - for block_start in range(0, len(raw), 2880): - block = raw[block_start: block_start + 2880] - header_blocks += 1 - for i in range(0, 2880, 80): - rec = block[i:i+80].decode("ascii", errors="replace") - key = rec[:8].strip() - if key == "END": - found_end = True - break - if "=" in rec[:30]: - k, _, rest = rec.partition("=") - val_str = rest.split("/")[0].strip().strip("'").strip() - try: - if "." in val_str: - header[k.strip()] = float(val_str) - elif val_str in ("T", "F"): - header[k.strip()] = (val_str == "T") - else: - header[k.strip()] = int(val_str) - except ValueError: - header[k.strip()] = val_str - if found_end: - break - - self.header = header - bitpix = int(header.get("BITPIX", -32)) - naxis1 = int(header.get("NAXIS1", 0)) - naxis2 = int(header.get("NAXIS2", 0)) - - data_start = header_blocks * 2880 - n_bytes = abs(bitpix) // 8 * naxis1 * naxis2 - dt = {8: ">u1", 16: ">u2", -16: ">u2", 32: ">i4", - -32: ">f4", -64: ">f8"}.get(bitpix, ">f4") - - self.array = np.frombuffer( - raw[data_start: data_start + n_bytes], dtype=dt - ).reshape(naxis2, naxis1) - - self._wcs = { - "crval1": float(header.get("CRVAL1", 0)), - "crval2": float(header.get("CRVAL2", 0)), - "crpix1": float(header.get("CRPIX1", naxis1 / 2)), - "crpix2": float(header.get("CRPIX2", naxis2 / 2)), - "cdelt1": float(header.get("CDELT1", -0.001042)), - "cdelt2": float(header.get("CDELT2", 0.001042)), - } + self.header = hdu.header + self.array = hdu.data + self._wcs = WCS(hdu.header) return True def world_to_pixel(self, ra: float, dec: float) -> Tuple[int, int]: - w = self._wcs - px = w["crpix1"] + (w["crval1"] - ra) / abs(w["cdelt1"]) - py = w["crpix2"] + (dec - w["crval2"]) / abs(w["cdelt2"]) + coord = SkyCoord(ra=ra, dec=dec, unit="deg") + px, py = coord.to_pixel(self._wcs) return int(round(px)), int(round(py)) def is_saturated(self, cx: int, cy: int, box: int = 5) -> Tuple[bool, float]: From 29be96ee9f9291d8a4796a55c1f5ec7474f47112 Mon Sep 17 00:00:00 2001 From: Caden Gobat <36030084+cgobat@users.noreply.github.com> Date: Tue, 17 Mar 2026 18:35:22 -0700 Subject: [PATCH 2/9] Read bayer pattern from FITS header rather than hardcoding it Doesn't seem like this is even used anywhere currently, but this is more robust regardless --- core/postflight/bayer_photometry.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/core/postflight/bayer_photometry.py b/core/postflight/bayer_photometry.py index d4b1d8a..9d53a8d 100644 --- a/core/postflight/bayer_photometry.py +++ b/core/postflight/bayer_photometry.py @@ -28,11 +28,6 @@ logger = logging.getLogger("seevar.bayer_photometry") -# IMX585 Bayer pattern — GRBG -# Row Even: G, R, G, R (col 0,1,2,3...) -# Row Odd: B, G, B, G -BAYER_PATTERN = "GRBG" - # Default aperture geometry (pixels) — used as fallback if PSF fit fails R_AP_DEFAULT = 8 R_SKY_IN_DEFAULT = 12 @@ -112,6 +107,7 @@ def __init__(self, fits_path: Path): self.header: fits.Header = None self.array: Optional[np.ndarray] = None self._wcs: WCS = None + self.bayer_pattern: str = None def load(self) -> bool: try: @@ -124,6 +120,12 @@ def load(self) -> bool: self.header = hdu.header self.array = hdu.data self._wcs = WCS(hdu.header) + + self.bayer_pattern = self.header.get("BAYERPAT", "GRBG") + # For IMX585 (S30 Pro), we expect GRBG, which looks like + # Even rows: G, R, G, R (col 0,1,2,3...) + # Odd rows: B, G, B, G + return True def world_to_pixel(self, ra: float, dec: float) -> Tuple[int, int]: From 55c15d3f32b912f2175e17f5cee7c261c692623a Mon Sep 17 00:00:00 2001 From: Caden Gobat <36030084+cgobat@users.noreply.github.com> Date: Tue, 17 Mar 2026 19:01:51 -0700 Subject: [PATCH 3/9] Use astropy to construct and write FITS files --- core/flight/pilot.py | 54 +++++++------------------------------------- 1 file changed, 8 insertions(+), 46 deletions(-) diff --git a/core/flight/pilot.py b/core/flight/pilot.py index f8876bb..69adb4a 100644 --- a/core/flight/pilot.py +++ b/core/flight/pilot.py @@ -18,6 +18,7 @@ from typing import Optional, Tuple import numpy as np +from astropy.io import fits from astropy.coordinates import EarthLocation, AltAz, SkyCoord, get_body from astropy.time import Time import astropy.units as u @@ -405,14 +406,6 @@ def sovereign_stamp( # --- Assemble header ----------------------------------------------------- h = { - # --- Mandatory FITS structural keys (priority order preserved by write_fits) --- - "SIMPLE": True, - "BITPIX": 16, - "NAXIS": 2, - "NAXIS1": width, - "NAXIS2": height, - "BZERO": 32768.0, # CRITICAL: unsigned uint16 → signed int16 FITS convention - "BSCALE": 1.0, # --- Target identity --- "OBJECT": target.name, "OBJCTRA": _hours_to_hms(target.ra_hours), @@ -474,48 +467,17 @@ def sovereign_stamp( def write_fits(array: np.ndarray, header_dict: dict, output_path: Path) -> bool: output_path.parent.mkdir(parents=True, exist_ok=True) - # 1. Offset the unsigned 16-bit IMX585 data into signed 16-bit space for FITS standard - array_offset = array.astype(np.int32) - 32768 - array_signed = array_offset.astype(np.int16) - - if array_signed.dtype.byteorder not in (">",): - array_signed = array_signed.byteswap().view(array_signed.dtype.newbyteorder(">")) - - def card(key: str, value, comment: str = "") -> str: - key = key.upper()[:8].ljust(8) - if isinstance(value, bool): val_str = f"{'T' if value else 'F':>20}" - elif isinstance(value, int): val_str = f"{value:>20}" - elif isinstance(value, float): val_str = f"{value:>20.10G}" - elif isinstance(value, str): - val_str = f"'{value.replace('\'', '\'\''):<8}'" - val_str = f"{val_str:<20}" - else: - val_str = f"'{str(value):<8}'" - val_str = f"{val_str:<20}" - c = f" / {comment}" if comment else "" - return f"{key}= {val_str}{c}"[:80].ljust(80) - - # BZERO and BSCALE must appear immediately after NAXIS properties - priority_keys = ["SIMPLE", "BITPIX", "NAXIS", "NAXIS1", "NAXIS2", "BZERO", "BSCALE"] - records = [card(k, header_dict[k]) for k in priority_keys if k in header_dict] - records.extend([card(k, v) for k, v in header_dict.items() if k not in priority_keys]) - - records.append("COMMENT SeeVar v5.1.0 M2 -- BZERO Signed-Integer Protected".ljust(80)) - records.append("END".ljust(80)) - - while (len(records) * 80) % 2880 != 0: records.append(" " * 80) - header_bytes = "".join(records).encode("ascii") + header = fits.Header() + header.update(header_dict) + header["COMMENT"] = f"{SWCREATE} M2" - data_bytes = array_signed.tobytes() - remainder = len(data_bytes) % 2880 - if remainder: data_bytes += b"\x00" * (2880 - remainder) + hdu = fits.PrimaryHDU(data=array.astype(np.uint16), header=header) try: - with open(output_path, "wb") as f: - f.write(header_bytes) - f.write(data_bytes) + hdu.writeto(output_path, overwrite=True) return True - except OSError: return False + except OSError: + return False class DiamondSequence: """Sovereign acquisition loop — STATE_MACHINE.md DiamondSequence. From 6d5b4e76d93763f2391e52c39ad5dc9d1f94ca2e Mon Sep 17 00:00:00 2001 From: Caden Gobat <36030084+cgobat@users.noreply.github.com> Date: Tue, 17 Mar 2026 19:02:47 -0700 Subject: [PATCH 4/9] Consistently handle time keywords using astropy `Time` object --- core/flight/pilot.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/core/flight/pilot.py b/core/flight/pilot.py index 69adb4a..6aefbd0 100644 --- a/core/flight/pilot.py +++ b/core/flight/pilot.py @@ -420,7 +420,9 @@ def sovereign_stamp( "CTYPE1": "RA---TAN", "CTYPE2": "DEC--TAN", # --- Timing --- - "DATE-OBS": utc_obs.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-3], + "DATE-OBS": t_astropy.isot, # ISO 8601 format + "MJD-OBS": t_astropy.mjd, # Modified Julian Date + "JD": round(t_astropy.jd, 6), # Julian Date "EXPTIME": target.exp_ms / 1000.0, # --- Instrument --- "INSTRUME": INSTRUMENT, @@ -458,10 +460,6 @@ def sovereign_stamp( if target.auid: h["AUID"] = target.auid - # JD — Julian Date of exposure start. - # Explicit key for AAVSO photometry tools and reporter cross-reference. - h["JD"] = round(t_astropy.jd, 6) - return h def write_fits(array: np.ndarray, header_dict: dict, output_path: Path) -> bool: From c8ec7749bd4a100c108c2067b7376eb7d8a44c6a Mon Sep 17 00:00:00 2001 From: Caden Gobat <36030084+cgobat@users.noreply.github.com> Date: Tue, 17 Mar 2026 19:03:38 -0700 Subject: [PATCH 5/9] Remove fake WCS --- core/flight/pilot.py | 9 --------- 1 file changed, 9 deletions(-) diff --git a/core/flight/pilot.py b/core/flight/pilot.py index 6aefbd0..5095d08 100644 --- a/core/flight/pilot.py +++ b/core/flight/pilot.py @@ -410,15 +410,6 @@ def sovereign_stamp( "OBJECT": target.name, "OBJCTRA": _hours_to_hms(target.ra_hours), "OBJCTDEC": _deg_to_dms(target.dec_deg), - # --- WCS --- - "CRVAL1": ra_deg, - "CRVAL2": target.dec_deg, - "CRPIX1": width / 2.0, - "CRPIX2": height / 2.0, - "CDELT1": -0.001042, - "CDELT2": 0.001042, - "CTYPE1": "RA---TAN", - "CTYPE2": "DEC--TAN", # --- Timing --- "DATE-OBS": t_astropy.isot, # ISO 8601 format "MJD-OBS": t_astropy.mjd, # Modified Julian Date From a98c0282c2deed165109c0da6caea22bfd68fe02 Mon Sep 17 00:00:00 2001 From: Caden Gobat <36030084+cgobat@users.noreply.github.com> Date: Tue, 17 Mar 2026 20:08:37 -0700 Subject: [PATCH 6/9] Clarify that RPC JSON should be terminated by CR+LF, not just newline --- dev/logic/API_PROTOCOL.MD | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dev/logic/API_PROTOCOL.MD b/dev/logic/API_PROTOCOL.MD index e88635a..1ab7f1d 100644 --- a/dev/logic/API_PROTOCOL.MD +++ b/dev/logic/API_PROTOCOL.MD @@ -25,7 +25,7 @@ S30-Pro's network address. Never hardcoded. ## JSON-RPC Wire Format (port 4700) -All messages are newline-terminated JSON, sent as UTF-8: +All messages are CRLF-terminated JSON, sent as UTF-8: ```python msg = {"id": , "method": "", "params": } From f3508ca00f80f2ddb51a318e92ec6a0ca371f0cf Mon Sep 17 00:00:00 2001 From: Caden Gobat <36030084+cgobat@users.noreply.github.com> Date: Tue, 17 Mar 2026 21:58:16 -0700 Subject: [PATCH 7/9] Use photutils for aperture photometry --- core/postflight/bayer_photometry.py | 37 ++++++++++++++++------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/core/postflight/bayer_photometry.py b/core/postflight/bayer_photometry.py index 9d53a8d..413d96b 100644 --- a/core/postflight/bayer_photometry.py +++ b/core/postflight/bayer_photometry.py @@ -25,6 +25,7 @@ from astropy.coordinates import SkyCoord from astropy.io import fits from astropy.wcs import WCS +from photutils.aperture import ApertureStats, CircularAnnulus, CircularAperture logger = logging.getLogger("seevar.bayer_photometry") @@ -67,29 +68,31 @@ def aperture_flux( Returns (net_flux, sky_median, snr). """ - y_idx, x_idx = np.ogrid[-cy:image.shape[0]-cy, -cx:image.shape[1]-cx] - abs_y = y_idx + cy - abs_x = x_idx + cx + data = np.asanyarray(image, dtype=np.float64) + yy, xx = np.indices(data.shape) - if bayer_channel == "G": - b_mask = (abs_y % 2) == (abs_x % 2) + if bayer_channel == "ALL": + b_mask = np.ones_like(yy, dtype=bool) + elif bayer_channel == "G": + b_mask = (yy % 2) == (xx % 2) elif bayer_channel == "R": - b_mask = (abs_y % 2 == 0) & (abs_x % 2 == 1) + b_mask = (yy % 2 == 0) & (xx % 2 == 1) elif bayer_channel == "B": - b_mask = (abs_y % 2 == 1) & (abs_x % 2 == 0) + b_mask = (yy % 2 == 1) & (xx % 2 == 0) else: - b_mask = np.ones_like(abs_y, dtype=bool) + raise ValueError(f"Invalid bayer channel: {bayer_channel!r}") - r2 = (x_idx**2 + y_idx**2).astype(np.float64) - ap_mask = (r2 <= r_ap ** 2) & b_mask - sky_mask = (r2 >= r_sky_in ** 2) & (r2 <= r_sky_out ** 2) & b_mask + exclusion_mask = ~b_mask + aperture = CircularAperture((cx, cy), r=r_ap) + sky_annulus = CircularAnnulus((cx, cy), r_in=r_sky_in, r_out=r_sky_out) - sky_vals = image[sky_mask].astype(np.float64) - sky_median = float(np.median(sky_vals)) if len(sky_vals) > 0 else 0.0 - sky_std = float(sky_vals.std()) if len(sky_vals) > 0 else 1.0 + sky_stats = ApertureStats(data, sky_annulus, mask=exclusion_mask) + sky_median = sky_stats.median + sky_std = sky_stats.std - ap_sum = float(image[ap_mask].astype(np.float64).sum()) - n_ap = int(ap_mask.sum()) + ap_stats = ApertureStats(data, aperture, mask=exclusion_mask) + ap_sum = ap_stats.sum + n_ap = ap_stats.sum_aper_area net_flux = ap_sum - sky_median * n_ap snr = net_flux / (sky_std * math.sqrt(n_ap)) if sky_std > 0 and n_ap > 0 else 0.0 @@ -224,6 +227,8 @@ def measure_star( "cx": cx, "cy": cy, "peak": peak, "saturated": saturated, "r_ap": r_ap, + "r_sky_in": r_sky_in_used, + "r_sky_out": r_sky_out_used, } for ch in ("G", "R", "B", "ALL"): From 1304768dd3f5eb473be0fc46b10a8c782cb62ec1 Mon Sep 17 00:00:00 2001 From: Caden Gobat <36030084+cgobat@users.noreply.github.com> Date: Tue, 17 Mar 2026 22:24:03 -0700 Subject: [PATCH 8/9] Use `np.*` instead of `math.*` functions for one less import --- core/postflight/bayer_photometry.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/core/postflight/bayer_photometry.py b/core/postflight/bayer_photometry.py index 413d96b..a2c7b94 100644 --- a/core/postflight/bayer_photometry.py +++ b/core/postflight/bayer_photometry.py @@ -17,7 +17,6 @@ """ import logging -import math import numpy as np from pathlib import Path from typing import Optional, Tuple @@ -94,7 +93,7 @@ def aperture_flux( ap_sum = ap_stats.sum n_ap = ap_stats.sum_aper_area net_flux = ap_sum - sky_median * n_ap - snr = net_flux / (sky_std * math.sqrt(n_ap)) if sky_std > 0 and n_ap > 0 else 0.0 + snr = net_flux / (sky_std * np.sqrt(n_ap)) if sky_std > 0 and n_ap > 0 else 0.0 return net_flux, sky_median, snr @@ -328,7 +327,7 @@ def differential_magnitude( if comp_snr < MIN_COMP_SNR: continue - zp = v_mag + 2.5 * math.log10(m[flux_key]) + zp = v_mag + 2.5 * np.log10(m[flux_key]) zero_points.append(zp) weights.append(comp_snr ** 2) # SNR² weighting @@ -347,11 +346,11 @@ def differential_magnitude( np.sqrt(np.sum(w_arr * (zp_arr - avg_zp) ** 2) / w_sum) ) - magnitude = avg_zp - 2.5 * math.log10(target_flux) + magnitude = avg_zp - 2.5 * np.log10(target_flux) # Photometric error: quadrature sum of weighted ZP scatter and SNR noise snr_err = 1.0857 / target_snr if target_snr > 0 else 9.99 - total_err = round(math.sqrt(zp_std ** 2 + snr_err ** 2), 3) + total_err = round(np.sqrt(zp_std ** 2 + snr_err ** 2), 3) return { "status": "ok", From b25eb0ccda0663e50b7350cc5c58870d18a44f14 Mon Sep 17 00:00:00 2001 From: Caden Gobat <36030084+cgobat@users.noreply.github.com> Date: Fri, 20 Mar 2026 09:04:28 -0700 Subject: [PATCH 9/9] Remove unused width & height args to `sovereign_stamp()` function --- core/flight/orchestrator.py | 2 +- core/flight/pilot.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/core/flight/orchestrator.py b/core/flight/orchestrator.py index 649a341..bc3ad65 100644 --- a/core/flight/orchestrator.py +++ b/core/flight/orchestrator.py @@ -132,7 +132,7 @@ def step(tag, msg): timestamp = utc_obs.strftime("%Y%m%dT%H%M%S") out_path = LOCAL_BUFFER / f"SIM_{safe_name}_{timestamp}_Raw.fits" - header = sovereign_stamp(target, utc_obs, width, height) + header = sovereign_stamp(target, utc_obs) write_fits(array, header, out_path) return FrameResult(success=True, path=out_path, width=width, height=height, elapsed_s=4.0) diff --git a/core/flight/pilot.py b/core/flight/pilot.py index 8ee23b7..76e7480 100644 --- a/core/flight/pilot.py +++ b/core/flight/pilot.py @@ -257,7 +257,7 @@ def _read_gps_ram() -> dict: return {"lat": lat, "lon": lon, "elevation": elev} except Exception: return {"lat": 0.0, "lon": 0.0, "elevation": 0.0} -def sovereign_stamp(target: AcquisitionTarget, utc_obs: datetime, width: int, height: int, ccd_temp: Optional[float] = None) -> dict: +def sovereign_stamp(target: AcquisitionTarget, utc_obs: datetime, ccd_temp: Optional[float] = None) -> dict: ra_deg = target.ra_hours * 15.0 t_astropy = Time(utc_obs) @@ -407,7 +407,7 @@ def notify(step, msg): safe_name = target.name.replace(" ", "_").replace("/", "-") out_path = LOCAL_BUFFER / f"{safe_name}_{utc_obs.strftime('%Y%m%dT%H%M%S')}_Raw.fits" - header = sovereign_stamp(target, utc_obs, width, height, ccd_temp=ccd_temp) + header = sovereign_stamp(target, utc_obs, ccd_temp=ccd_temp) ok = write_fits(array, header, out_path) elapsed = time.monotonic() - t_start