Use astropy for FITS I/O and WCS transformations + photutils for aperture photometry

core/flight/orchestrator.py

@@ -96,7 +96,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)
 
         step("T7", f"write_fits — SIM payload saved to {out_path}")

core/flight/pilot.py

@@ -30,6 +30,7 @@
 from urllib.parse import urlencode
 
 import numpy as np
+from astropy.io import fits
 from astropy.coordinates import AltAz, EarthLocation, SkyCoord, get_body
 from astropy.time import Time
 import astropy.units as u
@@ -552,10 +553,8 @@ def _read_gps_ram() -> dict:
     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:
-    ra_deg    = target.ra_hours * 15.0
+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)
     gps       = _read_gps_ram()
     site_lat, site_lon, site_elev = gps["lat"], gps["lon"], gps["elevation"]
@@ -585,18 +584,13 @@ def sovereign_stamp(target: AcquisitionTarget, utc_obs: datetime,
             pass
 
     h = {
-        "SIMPLE":   True,   "BITPIX": 16, "NAXIS": 2,
-        "NAXIS1":   width,  "NAXIS2": height,
-        "BZERO":    32768.0, "BSCALE": 1.0,
         "OBJECT":   target.name,
         "OBJCTRA":  _hours_to_hms(target.ra_hours),
         "OBJCTDEC": _deg_to_dms(target.dec_deg),
-        "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",
-        "DATE-OBS": utc_obs.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-3],
-        "EXPTIME":  target.exp_ms / 1000.0,
+        "DATE-OBS": t_astropy.isot,          # ISO 8601 format
+        "MJD-OBS":  round(t_astropy.mjd, 6), # Modified Julian Date
+        "JD":       round(t_astropy.jd, 6),  # Julian Date
+        "EXPTIME":  target.exp_ms / 1000.0,  # Exposure time (in seconds)
         "INSTRUME": INSTRUMENT,   "TELESCOP": TELESCOPE,
         "FILTER":   FILTER_NAME,  "BAYERPAT": BAYER_PATTERN,
         "GAIN":     GAIN,         "FOCALLEN": FOCALLEN,
@@ -611,42 +605,20 @@ def sovereign_stamp(target: AcquisitionTarget, utc_obs: datetime,
     if moon_phase is not None: h["MOONPHASE"] = moon_phase
     if moon_alt   is not None: h["MOONALT"]   = moon_alt
     if target.auid:            h["AUID"]      = target.auid
-    h["JD"] = round(t_astropy.jd, 6)
+
     return h
 
 def write_fits(array: np.ndarray, header_dict: dict, output_path: Path) -> bool:
     output_path.parent.mkdir(parents=True, exist_ok=True)
-    array_signed = (array.astype(np.int32) - 32768).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(chr(39), chr(39)*2):<8}'".ljust(20)
-        else:                          val_str = f"'{str(value):<8}'".ljust(20)
-        return f"{key}= {val_str}{f' / {comment}' if comment else ''}"[:80].ljust(80)
-
-    priority = ["SIMPLE", "BITPIX", "NAXIS", "NAXIS1", "NAXIS2", "BZERO", "BSCALE"]
-    records  = [card(k, header_dict[k]) for k in priority if k in header_dict]
-    records += [card(k, v) for k, v in header_dict.items() if k not in priority]
-    records.append("COMMENT   SeeVar v2.0.0 -- 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")
-    data_bytes   = array_signed.tobytes()
-    remainder    = len(data_bytes) % 2880
-    if remainder:
-        data_bytes += b"\x00" * (2880 - remainder)
+
+    header = fits.Header()
+    header.update(header_dict)
+    header["COMMENT"] = "SeeVar v2.0.0 -- BZERO Signed-Integer Protected"
+
+    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
@@ -778,7 +750,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
 

core/postflight/bayer_photometry.py

@@ -7,17 +7,16 @@
 """
 
 import logging
-import math
 import numpy as np
 from pathlib import Path
 from typing import Optional, Tuple
 
-logger = logging.getLogger("seevar.bayer_photometry")
+from astropy.coordinates import SkyCoord
+from astropy.io import fits
+from astropy.wcs import WCS
+from photutils.aperture import ApertureStats, CircularAnnulus, CircularAperture
 
-# IMX585 Bayer pattern — GRBG
-# Row Even: G, R, G, R  (col 0,1,2,3...)
-# Row Odd:  B, G, B, G
-BAYER_PATTERN = "GRBG"
+logger = logging.getLogger("seevar.bayer_photometry")
 
 # Default aperture geometry (pixels) — used as fallback if PSF fit fails
 R_AP_DEFAULT      = 8
@@ -58,115 +57,72 @@ 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
+    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
 
 
-# ---------------------------------------------------------------------------
-# 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
+        self.bayer_pattern: str           = 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)
+
+        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]:
-        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]:
@@ -260,6 +216,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"):
@@ -359,7 +317,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
 
@@ -378,11 +336,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",