WCS/catalog refactor, photometry and zeropoint reliability, and FITS header fixes

.gitignore

@@ -1,6 +1,10 @@
+# Examples: entire directory local-only (scripts, data, outputs); not versioned
+examples/
+
 # Python
 *.pyc
 __pycache__/
+.coverage
 *.pyo
 *.pyd
 *.pyz

docs/source/api/generated/potpyri.primitives.solve_wcs.rst

@@ -7,9 +7,10 @@
    .. rubric:: Functions
 
    .. autosummary::
-
-      align_to_gaia
+
       clean_up_astrometry
-      get_gaia_catalog
+      fine_align_wcs
+      get_fine_align_reference_catalog
       solve_astrometry
+      validate_existing_wcs_header
 

potpyri/_version.py

@@ -1,5 +1,6 @@
-# file generated by setuptools-scm
+# file generated by vcs-versioning
 # don't change, don't track in version control
+from __future__ import annotations
 
 __all__ = [
     "__version__",
@@ -10,25 +11,14 @@
     "commit_id",
 ]
 
-TYPE_CHECKING = False
-if TYPE_CHECKING:
-    from typing import Tuple
-    from typing import Union
-
-    VERSION_TUPLE = Tuple[Union[int, str], ...]
-    COMMIT_ID = Union[str, None]
-else:
-    VERSION_TUPLE = object
-    COMMIT_ID = object
-
 version: str
 __version__: str
-__version_tuple__: VERSION_TUPLE
-version_tuple: VERSION_TUPLE
-commit_id: COMMIT_ID
-__commit_id__: COMMIT_ID
+__version_tuple__: tuple[int | str, ...]
+version_tuple: tuple[int | str, ...]
+commit_id: str | None
+__commit_id__: str | None
 
-__version__ = version = '1.0.17.dev6'
-__version_tuple__ = version_tuple = (1, 0, 17, 'dev6')
+__version__ = version = '1.0.17.dev18'
+__version_tuple__ = version_tuple = (1, 0, 17, 'dev18')
 
-__commit_id__ = commit_id = 'gac532c0e6'
+__commit_id__ = commit_id = 'g185939e0b'

potpyri/instruments/GMOS.py

@@ -150,12 +150,20 @@ def get_instrument_name(self, hdr):
 
     def import_image(self, filename, amp, log=None):
 
-        with fits.open(filename) as hdr:
-            header = hdr[0].header
+        # Load extensions without CCDData.read(...): that parses WCS from each
+        # extension header and raises FITSFixedWarning for deprecated RADECSYS /
+        # MJD-OBS. Use fixed headers and wcs=None until a reduced WCS exists.
+        with fits.open(filename, memmap=False) as hdr:
+            header = hdr[0].header.copy()
+            instrument.fix_deprecated_wcs_header_cards(header)
             binn = self.get_binning(hdr[1].header)
-
-        raw = [CCDData.read(filename, hdu=x+1, unit='adu') 
-            for x in range(int(amp))]
+            raw = []
+            for x in range(int(amp)):
+                hdu = hdr[x + 1]
+                meta = hdu.header.copy()
+                instrument.fix_deprecated_wcs_header_cards(meta)
+                raw.append(
+                    CCDData(np.asarray(hdu.data), meta=meta, unit=u.adu, wcs=None))
         red = [ccdproc.ccd_process(x, oscan=x.header['BIASSEC'], 
             oscan_model=models.Chebyshev1D(3), 
             trim=x.header['DATASEC'], 

potpyri/instruments/instrument.py

@@ -37,6 +37,52 @@
 )
 
 
+def fix_deprecated_wcs_header_cards(header):
+    """Normalize WCS-related FITS cards before :class:`~astropy.wcs.WCS` parses the header.
+
+    Archive and Gemini raw frames often use deprecated ``RADECSYS`` instead of
+    ``RADESYS``, and leave ``MJD-OBS`` in a form that makes WCSLib run ``datfix``
+    and emit ``FITSFixedWarning``. Updating the header in place avoids those
+    warnings without globally filtering them.
+
+    Parameters
+    ----------
+    header : astropy.io.fits.Header
+        Header to modify in place (pass a copy if the on-disk HDU must be
+        unchanged).
+    """
+    if header is None:
+        return
+
+    if 'RADECSYS' in header:
+        if 'RADESYS' not in header:
+            header['RADESYS'] = str(header['RADECSYS']).strip()
+        try:
+            del header['RADECSYS']
+        except KeyError:
+            pass
+
+    if 'RADESYS' in header:
+        v = header['RADESYS']
+        if isinstance(v, str):
+            header['RADESYS'] = v.strip()
+
+    if 'MJD-OBS' in header:
+        try:
+            mjd = float(header['MJD-OBS'])
+            t = Time(mjd, format='mjd', scale='utc')
+            isot = t.isot
+            if 'T' in isot:
+                date_part, time_part = isot.split('T', 1)
+                header['DATE-OBS'] = date_part
+                header['TIME-OBS'] = time_part
+            # Having both MJD-OBS and DATE-OBS can still trigger WCSLib datfix;
+            # remove MJD-OBS from this in-memory header once UTC date/time are set.
+            del header['MJD-OBS']
+        except (ValueError, TypeError, KeyError):
+            pass
+
+
 def _sanitize_calibration_header(header):
     """Remove WCS and coordinate keywords from a header so it can be written safely.
 
@@ -65,7 +111,9 @@ def _read_calibration_ccd(path, unit, hdu_index=0):
     """
     with fits.open(path) as hdu_list:
         hdu = hdu_list[hdu_index]
-        return CCDData(hdu.data, meta=hdu.header.copy(), unit=unit, wcs=None)
+        meta = hdu.header.copy()
+        fix_deprecated_wcs_header_cards(meta)
+        return CCDData(hdu.data, meta=meta, unit=unit, wcs=None)
 
 
 class Instrument(object):
@@ -946,8 +994,23 @@ def create_sky(self, sky_list, fil, amp, binn, paths, log=None,
 
             skys.append(sky_full)
 
-        msky = ccdproc.combine(skys, method='median', sigma_clip=True,
-            clip_extrema=True)
+        n_sky = len(skys)
+        if n_sky == 0:
+            if log:
+                log.error('create_sky: no valid sky frames after masking; skipping master sky.')
+            return
+
+        # Per-frame masking already sigma-clips; stacking sigma_clip needs >=3
+        # samples per pixel to be stable. With 1–2 frames it produces all-NaN
+        # slices in ccdproc (RuntimeWarning) and no real outlier rejection gain.
+        if n_sky == 1:
+            msky = skys[0].copy()
+        elif n_sky == 2:
+            msky = ccdproc.combine(
+                skys, method='median', sigma_clip=False, clip_extrema=False)
+        else:
+            msky = ccdproc.combine(
+                skys, method='median', sigma_clip=True, clip_extrema=True)
 
         # Robust masking of combined master sky: iterative sigma clipping
         # for accurate median/std, then mask excess flux and defects
@@ -959,6 +1022,10 @@ def create_sky(self, sky_list, fil, amp, binn, paths, log=None,
             sigma_lower=3.0,
             maxiters=msky_maxiters,
         )
+        if not np.isfinite(stddev) or stddev <= 0:
+            stddev = float(np.nanstd(msky.data, dtype=float))
+        if not np.isfinite(stddev) or stddev <= 0:
+            stddev = 1e-12
         mask_high = msky.data > median + msky_n_sigma_high * stddev
         mask_low = msky.data < median - msky_n_sigma_low * stddev
         msky.data[mask_high] = 1.0
@@ -1022,10 +1089,6 @@ def expand_mask(self, input_data, input_mask=None):
         np.ndarray
             Boolean mask (True = bad).
         """
-        # Get image statistics
-        mean, median, stddev = sigma_clipped_stats(input_data.data)
-
-        # Add to input mask
         add_mask = np.isnan(input_data.data)
         add_mask = add_mask  | np.isinf(input_data.data)
         add_mask = add_mask  | (input_data.data==0.0) 
@@ -1150,7 +1213,12 @@ def process_science(self, sci_list, fil, amp, binn, paths, mbias=None,
                 final_img.header['SKYBKG'] = 0.0
 
             # Apply final masking based on excessively negative values
-            mean, median, stddev = sigma_clipped_stats(final_img.data)
+            finite = np.isfinite(final_img.data)
+            if np.any(finite):
+                mean, median, stddev = sigma_clipped_stats(
+                    final_img.data, mask=~finite)
+            else:
+                mean = median = stddev = np.nan
             mask = final_img.data < median - 10 * stddev
             final_img.data[mask] = np.nan
             final_img.mask = final_img.mask | mask