diff --git a/.gitignore b/.gitignore index 2504024..4dc879b 100644 --- a/.gitignore +++ b/.gitignore @@ -5,6 +5,8 @@ examples/ *.pyc __pycache__/ .coverage +# Generated at build/install time by setuptools_scm (see pyproject.toml) +potpyri/_version.py *.pyo *.pyd *.pyz diff --git a/potpyri/_version.py b/potpyri/_version.py deleted file mode 100644 index e9a3f8a..0000000 --- a/potpyri/_version.py +++ /dev/null @@ -1,24 +0,0 @@ -# file generated by vcs-versioning -# don't change, don't track in version control -from __future__ import annotations - -__all__ = [ - "__version__", - "__version_tuple__", - "version", - "version_tuple", - "__commit_id__", - "commit_id", -] - -version: str -__version__: str -__version_tuple__: tuple[int | str, ...] -version_tuple: tuple[int | str, ...] -commit_id: str | None -__commit_id__: str | None - -__version__ = version = '1.0.17.dev18' -__version_tuple__ = version_tuple = (1, 0, 17, 'dev18') - -__commit_id__ = commit_id = 'g185939e0b' diff --git a/potpyri/instruments/BINOSPEC.py b/potpyri/instruments/BINOSPEC.py index d90cdaa..e4ee84e 100755 --- a/potpyri/instruments/BINOSPEC.py +++ b/potpyri/instruments/BINOSPEC.py @@ -91,11 +91,10 @@ def __init__(self): self.out_size = 5000 - def raw_format(self, proc): + def _default_raw_format(self, proc): if proc: - return('sci_img_*proc.fits*') - else: - return('sci_img*[!proc].fits*') + return 'sci_img_*proc.fits*' + return 'sci_img*[!proc].fits*' # Get a unique image number that can be derived only from the file header def get_number(self, hdr): diff --git a/potpyri/instruments/DEIMOS.py b/potpyri/instruments/DEIMOS.py index 3db0835..98b02d3 100755 --- a/potpyri/instruments/DEIMOS.py +++ b/potpyri/instruments/DEIMOS.py @@ -89,11 +89,10 @@ def __init__(self): self.out_size = 9000 - def raw_format(self, proc): - if proc and str(proc)=='raw': - return('d*.fits') - else: - return('DE*.fits.gz') + def _default_raw_format(self, proc): + if proc and str(proc) == 'raw': + return 'd*.fits' + return 'DE*.fits.gz' def get_number(self, hdr): elap = Time(hdr['MJD'], format='mjd')-Time('1980-01-01') diff --git a/potpyri/instruments/F2.py b/potpyri/instruments/F2.py index b8529ee..9925540 100755 --- a/potpyri/instruments/F2.py +++ b/potpyri/instruments/F2.py @@ -92,8 +92,9 @@ def get_saturation(self, hdr): # Gives the saturation level in e- return(self.saturation*hdr['NREADS']*hdr['GAIN']) - def raw_format(self, proc): - return('*.fits.bz2') + def _default_raw_format(self, proc): + """Gemini archive default: ``*.fits.bz2`` (use ``--proc fits`` for ``*.fits``).""" + return '*.fits.bz2' def get_filter(self, hdr): filt = hdr['FILTER'].replace(' ','').split('_')[0] diff --git a/potpyri/instruments/FOURSTAR.py b/potpyri/instruments/FOURSTAR.py index 6c7c9d7..73474c4 100755 --- a/potpyri/instruments/FOURSTAR.py +++ b/potpyri/instruments/FOURSTAR.py @@ -96,8 +96,8 @@ def get_number(self, hdr): elap = int(np.round(elap.to(u.second).value)) return(elap) - def raw_format(self, proc): - return('*.fits') + def _default_raw_format(self, proc): + return '*.fits' def get_rdnoise(self, hdr): return(hdr['RDNOISE']) diff --git a/potpyri/instruments/GMOS.py b/potpyri/instruments/GMOS.py index 298484b..2186a32 100755 --- a/potpyri/instruments/GMOS.py +++ b/potpyri/instruments/GMOS.py @@ -110,11 +110,10 @@ def get_gain(self, hdr): except KeyError: return(self.gain) - def raw_format(self, proc): - if str(proc).lower()=='dragons': - return('*.fits') - else: - return('*.fits.bz2') + def _default_raw_format(self, proc): + if str(proc).lower() == 'dragons': + return '*.fits' + return '*.fits.bz2' def get_ampl(self, hdr): nccd = hdr['NCCDS'] diff --git a/potpyri/instruments/IMACS.py b/potpyri/instruments/IMACS.py index 1865e3f..83832eb 100755 --- a/potpyri/instruments/IMACS.py +++ b/potpyri/instruments/IMACS.py @@ -117,8 +117,8 @@ def get_ampl(self, hdr): return(amp) # Raw image format for ingestion - def raw_format(self, proc): - return('iff*.fits*') + def _default_raw_format(self, proc): + return 'iff*.fits*' def import_image(self, filename, amp, log=None): filename = os.path.abspath(filename) diff --git a/potpyri/instruments/LRIS.py b/potpyri/instruments/LRIS.py index 1e19129..cc87bf0 100755 --- a/potpyri/instruments/LRIS.py +++ b/potpyri/instruments/LRIS.py @@ -97,14 +97,12 @@ def get_number(self, header): number = str(header['FRAMENO']).zfill(5) return(number) - def raw_format(self, proc): - - if str(proc)=='archive': - return('*.fits*') - elif str(proc)=='raw': - return('*[b,r]*.fits') - else: - return('*.fits*') + def _default_raw_format(self, proc): + if str(proc) == 'archive': + return '*.fits*' + if str(proc) == 'raw': + return '*[b,r]*.fits' + return '*.fits*' def get_instrument_name(self, hdr): instrument = hdr['INSTRUME'] diff --git a/potpyri/instruments/MMIRS.py b/potpyri/instruments/MMIRS.py index 968d21c..58562d1 100755 --- a/potpyri/instruments/MMIRS.py +++ b/potpyri/instruments/MMIRS.py @@ -95,8 +95,8 @@ def get_number(self, hdr): elap = int(np.round(elap.to(u.second).value)) return(elap) - def raw_format(self, proc): - return('*.fits') + def _default_raw_format(self, proc): + return '*.fits' def get_rdnoise(self, hdr): return(hdr['RDNOISE']) diff --git a/potpyri/instruments/MOSFIRE.py b/potpyri/instruments/MOSFIRE.py index e4250da..83dabc5 100755 --- a/potpyri/instruments/MOSFIRE.py +++ b/potpyri/instruments/MOSFIRE.py @@ -91,8 +91,8 @@ def __init__(self): def get_saturation(self, hdr): return(hdr['SATURATE']*hdr['SYSGAIN']) - def raw_format(self, proc): - return('*.fits.gz') + def _default_raw_format(self, proc): + return '*.fits.gz' def get_filter(self, hdr): filt = hdr['FILTER'].replace(' ','').split('_')[0] diff --git a/potpyri/instruments/__init__.py b/potpyri/instruments/__init__.py index d3a721b..13053e4 100755 --- a/potpyri/instruments/__init__.py +++ b/potpyri/instruments/__init__.py @@ -1,71 +1,129 @@ """Instrument implementations and factory for POTPyRI. -Each instrument module defines a subclass of Instrument with detector keywords, -calibration behavior, and file-sorting rules. New instruments must be added -here and to __all__. +Each instrument module defines a subclass of :class:`~potpyri.instruments.instrument.Instrument` +with detector keywords, calibration behavior, and file-sorting rules. To add a +new instrument, create its module and add the canonical name to ``__all__`` below. """ -from . import BINOSPEC -from . import DEIMOS -from . import F2 -from . import FOURSTAR -from . import GMOS -from . import IMACS -from . import LRIS -from . import MMIRS -from . import MOSFIRE +from __future__ import annotations -# New instruments need to be added here and to the function below -__all__ = ["BINOSPEC","DEIMOS","F2","FOURSTAR","GMOS","IMACS","LRIS","MMIRS","MOSFIRE"] +from importlib import import_module +from typing import TYPE_CHECKING +if TYPE_CHECKING: + from .instrument import Instrument -def instrument_getter(instname, log=None): - """Return an Instrument instance for the given instrument name. +# Canonical instrument names (single source of truth). Each name must match a +# submodule ``potpyri.instruments.`` with class ````. +__all__ = [ + 'BINOSPEC', + 'DEIMOS', + 'F2', + 'FOURSTAR', + 'GMOS', + 'IMACS', + 'LRIS', + 'MMIRS', + 'MOSFIRE', +] + +# Optional CLI / shorthand aliases (not duplicated in __all__). +_INSTRUMENT_ALIASES: dict[str, str] = { + 'BINO': 'BINOSPEC', + 'MMIR': 'MMIRS', +} + +_MODULES = {name: import_module(f'.{name}', __name__) for name in __all__} + +# Re-export instrument submodules: ``from potpyri.instruments import GMOS`` +for _name, _mod in _MODULES.items(): + globals()[_name] = _mod + + +class UnknownInstrumentError(ValueError): + """Raised when an instrument name is not supported.""" + + def __init__(self, name: str, *, resolved: str | None = None) -> None: + self.name = name + self.resolved = resolved + supported = ', '.join(__all__) + aliases = ', '.join(f'{k}->{v}' for k, v in sorted(_INSTRUMENT_ALIASES.items())) + hint = f' Supported instruments: {supported}.' + if aliases: + hint += f' Aliases: {aliases}.' + if resolved and resolved != name.strip().upper(): + detail = f'{name!r} (resolved to {resolved!r})' + else: + detail = repr(name) + super().__init__(f'Instrument {detail} is not supported by POTPyRI.{hint}') + + +def supported_instruments() -> tuple[str, ...]: + """Return supported canonical instrument names.""" + return tuple(__all__) + + +def resolve_instrument_name(name: str) -> str: + """Normalize *name* to a canonical instrument key (uppercase, aliases applied). + + Parameters + ---------- + name : str + User-supplied instrument name or alias (e.g. ``'gmos'``, ``'BINO'``). + + Returns + ------- + str + Canonical name from :data:`__all__` when recognized, otherwise the + uppercased input (may still be unknown to :func:`instrument_getter`). + """ + key = name.strip().upper() + if key in __all__: + return key + if key in _INSTRUMENT_ALIASES: + return _INSTRUMENT_ALIASES[key] + # Legacy substring aliases used by older scripts / partial names. + if 'BINO' in key: + return 'BINOSPEC' + if 'MMIR' in key: + return 'MMIRS' + return key + + +def instrument_getter(instname: str, log=None) -> Instrument | None: + """Return an :class:`~potpyri.instruments.instrument.Instrument` for *instname*. Parameters ---------- instname : str - Instrument name (e.g. 'GMOS', 'LRIS', 'BINOSPEC'). + Instrument name or alias (e.g. ``'GMOS'``, ``'BINO'``, ``'MMIR'``). log : ColoredLogger, optional - Logger for error messages. If None, raises Exception on unsupported - instrument. + If provided and the instrument is unknown, log an error and return + ``None`` instead of raising. Returns ------- - Instrument - Subclass instance for the requested instrument. + Instrument or None + Configured instrument instance, or ``None`` if unknown and *log* is set. Raises ------ - Exception - If instname is not in __all__ and log is None. + UnknownInstrumentError + If *instname* is not supported and *log* is ``None``. + TypeError + If *instname* is missing or not a string. """ - instname = instname.upper() + if not instname or not isinstance(instname, str): + raise TypeError( + f'instrument name must be a non-empty string, got {instname!r}' + ) - if instname not in __all__: - if log: - log.error(f'Instrument {instname} not supported by POTPyRI') - else: - raise Exception(f'Instrument {instname} not supported by POTPyRI') - - tel = None - - if instname=="BINOSPEC": - tel = BINOSPEC.BINOSPEC() - elif instname=="DEIMOS": - tel = DEIMOS.DEIMOS() - elif instname=="F2": - tel = F2.F2() - elif instname=="FOURSTAR": - tel = FOURSTAR.FOURSTAR() - elif instname=="GMOS": - tel = GMOS.GMOS() - elif instname=="IMACS": - tel = IMACS.IMACS() - elif instname=="LRIS": - tel = LRIS.LRIS() - elif instname=="MMIRS": - tel = MMIRS.MMIRS() - elif instname=="MOSFIRE": - tel = MOSFIRE.MOSFIRE() - - return(tel) + canonical = resolve_instrument_name(instname) + if canonical not in __all__: + exc = UnknownInstrumentError(instname, resolved=canonical) + if log is not None: + log.error(str(exc)) + return None + raise exc + + cls = getattr(_MODULES[canonical], canonical) + return cls() diff --git a/potpyri/instruments/instrument.py b/potpyri/instruments/instrument.py index 9340dc1..cc6b611 100755 --- a/potpyri/instruments/instrument.py +++ b/potpyri/instruments/instrument.py @@ -7,6 +7,7 @@ """ from potpyri._version import __version__ +import glob import os import astropy import datetime @@ -25,6 +26,12 @@ from astropy.stats import SigmaClip from astropy.time import Time +# Per-frame optical background subtraction modes (see apply_optical_background_subtraction). +BKG_SUB_LOCAL = 'local' +BKG_SUB_CONSTANT = 'constant' +BKG_SUB_NONE = 'none' +VALID_BKG_SUB_MODES = (BKG_SUB_LOCAL, BKG_SUB_CONSTANT, BKG_SUB_NONE) + # WCS-related keywords to remove from calibration headers so saved files never # trigger InvalidTransformError (e.g. DEC=-100, ill-conditioned CD matrix). _CAL_HEADER_WCS_KEYS = ( @@ -124,6 +131,10 @@ class Instrument(object): science processing. Subclasses override attributes as needed. """ + # ``--proc fits`` (and aliases) force this glob on every instrument. + PROC_FITS_GLOB = '*.fits' + PROC_FITS_ALIASES = frozenset({'fits', '.fits', 'uncompressed'}) + def __init__(self): # Version @@ -388,12 +399,92 @@ def format_datasec(self, sec_string, binning=1): return(sec_string) + def _proc_requests_fits(self, proc): + """True when ``--proc`` requests uncompressed ``*.fits`` (all instruments).""" + return proc is not None and str(proc).lower() in self.PROC_FITS_ALIASES + def raw_format(self, proc): - """Return glob pattern for raw files (e.g. 'sci_img_*.fits' or 'sci_img*[!proc].fits').""" + """Return glob pattern for raw file discovery under ``data/raw``, ``data``, ``bad``. + + If ``proc`` is ``fits`` / ``.fits`` / ``uncompressed``, returns ``*.fits`` + regardless of instrument. Otherwise delegates to :meth:`_default_raw_format`. + """ + if self._proc_requests_fits(proc): + return self.PROC_FITS_GLOB + return self._default_raw_format(proc) + + def _default_raw_format(self, proc): + """Instrument-specific glob; override in subclasses.""" if proc: - return('sci_img_*.fits') + return 'sci_img_*.fits' + return 'sci_img*[!proc].fits' + + def discover_raw_files(self, paths, proc=None): + """Glob raw inputs under ``paths['raw']``, ``paths['data']``, and ``paths['bad']``. + + Returns + ------- + dict + ``pattern``, ``instrument_pattern``, ``proc``, ``proc_fits_override``, + ``per_dir`` (list of label, directory, glob_path, matched paths), + and ``files`` (combined matches). + """ + pattern = self.raw_format(proc) + instrument_pattern = self._default_raw_format(proc) + per_dir = [] + files = [] + for label, key in (('raw', 'raw'), ('data', 'data'), ('bad', 'bad')): + directory = paths[key] + glob_path = os.path.join(directory, pattern) + matched = sorted(glob.glob(glob_path)) + per_dir.append((label, directory, glob_path, matched)) + files.extend(matched) + return { + 'pattern': pattern, + 'instrument_pattern': instrument_pattern, + 'proc': proc, + 'proc_fits_override': self._proc_requests_fits(proc), + 'per_dir': per_dir, + 'files': files, + } + + def format_raw_discovery_message(self, paths, proc=None): + """Human-readable summary of where raw files are searched and what matched.""" + info = self.discover_raw_files(paths, proc) + n_total = len(info['files']) + lines = [ + f'Raw input discovery ({self.name}, --proc={info["proc"]!r}):', + ( + f' Glob pattern: {info["pattern"]!r} ' + '(Python glob: * = any substring; ? = one character; ' + '[seq] = one character in seq)' + ), + ] + if info['proc_fits_override']: + lines.append( + f' --proc fits: using {self.PROC_FITS_GLOB!r} on all instruments ' + f'(instrument-specific pattern for this --proc would be ' + f'{info["instrument_pattern"]!r})' + ) else: - return('sci_img*[!proc].fits') + lines.append( + f' Instrument pattern for --proc={info["proc"]!r}: ' + f'{info["instrument_pattern"]!r}' + ) + lines.append(' Search directories (pattern applied in each):') + for label, directory, glob_path, matched in info['per_dir']: + lines.append( + f' [{label}] {glob_path} -> {len(matched)} file(s)' + ) + if n_total == 0: + lines.append( + ' No files matched. Place raw data under data/raw/ or data/ ' + f'with basenames matching {info["pattern"]!r}. ' + 'For uncompressed .fits on any instrument, pass --proc fits.' + ) + else: + lines.append(f' Total: {n_total} raw file(s) to classify.') + return '\n'.join(lines) def get_stk_name(self, hdr, red_path): """Return full path for stacked output FITS (target.filter.utYYMMDD.amp.binn.stk.fits).""" @@ -1101,8 +1192,118 @@ def expand_mask(self, input_data, input_mask=None): return(input_mask) + def apply_optical_background_subtraction(self, processed_data, bkg_sub='local', + save_bkg=False, paths=None, log=None): + """Subtract per-frame background for optical data (not NIR master sky). + + Parameters + ---------- + processed_data : ccdproc.CCDData + Bias/dark/flat-corrected science frame with mask set. + bkg_sub : str, optional + ``'local'`` (default): photutils ``Background2D`` mesh model. + ``'constant'``: single sigma-clipped median subtracted from every pixel. + ``'none'``: no subtraction; ``SKYBKG`` set to 0. + save_bkg : bool, optional + If True and ``bkg_sub='local'``, write the 2D background model. + paths : dict, optional + Paths dict (``work`` key) required when ``save_bkg`` is True. + log : ColoredLogger, optional + Logger for progress. + + Returns + ------- + ccdproc.CCDData + Background-subtracted frame with ``SKYBKG`` and ``BKGSUB`` header keys. + """ + mode = (bkg_sub or BKG_SUB_LOCAL).lower() + if mode not in VALID_BKG_SUB_MODES: + raise ValueError( + f'bkg_sub must be one of {VALID_BKG_SUB_MODES}, got {bkg_sub!r}' + ) + + if mode == BKG_SUB_NONE: + final_img = processed_data + final_img.header['SKYBKG'] = 0.0 + final_img.header['BKGSUB'] = BKG_SUB_NONE + return final_img + + if mode == BKG_SUB_CONSTANT: + if log: + log.info('Calculating constant frame background.') + _, med, _ = sigma_clipped_stats( + processed_data.data, + mask=processed_data.mask, + sigma=3.0, + maxiters=5, + ) + if not np.isfinite(med): + med = np.nanmedian(processed_data.data) + background = med * u.electron + if log: + log.info(f'Constant background value: {background}') + final_img = processed_data.subtract(background) + final_img.header = processed_data.header + if 'SATURATE' in final_img.header: + final_img.header['SATURATE'] -= background.value + final_img.header['SKYBKG'] = background.value + final_img.header['BKGSUB'] = BKG_SUB_CONSTANT + return final_img + + if log: + log.info('Calculating 2D background.') + approx_background = np.nanmedian(processed_data.data) * u.electron + if log: + log.info(f'Approximate background value: {approx_background}') + bkg = Background2D( + processed_data, + (64, 64), + filter_size=(3, 3), + sigma_clip=SigmaClip(sigma=3), + exclude_percentile=80, + bkg_estimator=MeanBackground(), + mask=processed_data.mask, + fill_value=np.nanmedian(processed_data.data), + ) + + med_background = np.nanmedian(bkg.background) + if log: + log.info(f'Median background: {med_background}') + + if np.isnan(med_background.value): + final_img = processed_data.subtract(approx_background) + final_img.header = processed_data.header + if 'SATURATE' in final_img.header: + final_img.header['SATURATE'] -= approx_background.value + final_img.header['SKYBKG'] = approx_background.value + final_img.header['BKGSUB'] = BKG_SUB_CONSTANT + else: + if save_bkg: + if paths is None: + raise ValueError('paths is required when save_bkg=True') + bkg_filename = self.get_bkg_name(processed_data.header, paths['work']) + if log: + log.info(f'Writing background file: {bkg_filename}') + bkg_hdu = fits.PrimaryHDU(bkg.background.value) + bkg_hdu.header = processed_data.header + bkg_hdu.writeto(bkg_filename, overwrite=True, output_verify='silentfix') + + final_img = processed_data.subtract( + CCDData(bkg.background, unit=u.electron), + propagate_uncertainties=True, + handle_meta='first_found', + ) + if 'SATURATE' in final_img.header: + final_img.header['SATURATE'] -= med_background.value + final_img.header['SKYBKG'] = med_background.value + final_img.header['BKGSUB'] = BKG_SUB_LOCAL + if log: + log.info('Updating background and saturation values.') + return final_img + def process_science(self, sci_list, fil, amp, binn, paths, mbias=None, - mflat=None, mdark=None, skip_skysub=False, save_bkg=False, log=None): + mflat=None, mdark=None, skip_skysub=False, bkg_sub='local', + save_bkg=False, log=None): """Reduce science frames: bias/dark/flat, optional sky subtraction; return CCDData list. Parameters @@ -1124,7 +1325,13 @@ def process_science(self, sci_list, fil, amp, binn, paths, mbias=None, mdark : ccdproc.CCDData, optional Master dark. skip_skysub : bool, optional - If True, skip 2D background subtraction. Default is False. + If True, skip 2D background subtraction. Deprecated; use + ``bkg_sub='none'``. Default is False. + bkg_sub : str, optional + Per-frame background subtraction for optical data: ``'local'`` + (2D mesh, default), ``'constant'`` (single value per frame), or + ``'none'``. Ignored when ``needs_sky_subtraction(fil)`` is True + (NIR / GMOS z-band use master sky instead). save_bkg : bool, optional If True, save background model. Default is False. log : ColoredLogger, optional @@ -1137,6 +1344,8 @@ def process_science(self, sci_list, fil, amp, binn, paths, mbias=None, """ processed = [] processed_names = [] + if skip_skysub: + bkg_sub = BKG_SUB_NONE for sci in sorted(sci_list): if log: log.info(f'Importing {sci}') sci_full = self.import_image(sci, amp, log=log) @@ -1176,41 +1385,18 @@ def process_science(self, sci_list, fil, amp, binn, paths, mbias=None, processed_data.data[processed_data.mask]=np.nan if log: log.info(f'Wavelength is {self.wavelength}') - if not skip_skysub and not self.needs_sky_subtraction(fil): - if log: log.info('Calculating 2D background.') - approx_background=np.nanmedian(processed_data.data) * u.electron - if log: log.info(f'Approximate background value: {approx_background}') - bkg = Background2D(processed_data, (64, 64), filter_size=(3, 3), - sigma_clip=SigmaClip(sigma=3), exclude_percentile=80, - bkg_estimator=MeanBackground(), mask=processed_data.mask, - fill_value=np.nanmedian(processed_data.data)) - - med_background = np.nanmedian(bkg.background) - if log: log.info(f'Median background: {med_background}') - - if np.isnan(med_background.value): - final_img = processed_data.subtract(approx_background) - final_img.header = processed_data.header - final_img.header['SATURATE'] -= approx_background.value - final_img.header['SKYBKG'] = approx_background.value - else: - if save_bkg: - bkg_filename = self.get_bkg_name(processed_data.header, paths['work']) - if log: log.info(f'Writing background file: {bkg_filename}') - bkg_hdu = fits.PrimaryHDU(bkg.background.value) - bkg_hdu.header = processed_data.header - bkg_hdu.writeto(bkg_filename, overwrite=True, - output_verify='silentfix') - - final_img = processed_data.subtract(CCDData(bkg.background, - unit=u.electron), propagate_uncertainties=True, - handle_meta='first_found') - final_img.header['SATURATE'] -= med_background.value - final_img.header['SKYBKG'] = med_background.value - if log: log.info('Updating background and saturation values.') + if not self.needs_sky_subtraction(fil): + final_img = self.apply_optical_background_subtraction( + processed_data, + bkg_sub=bkg_sub, + save_bkg=save_bkg, + paths=paths, + log=log, + ) else: final_img = processed_data final_img.header['SKYBKG'] = 0.0 + final_img.header['BKGSUB'] = BKG_SUB_NONE # Apply final masking based on excessively negative values finite = np.isfinite(final_img.data) diff --git a/potpyri/primitives/image_procs.py b/potpyri/primitives/image_procs.py index e6b577b..81bb915 100755 --- a/potpyri/primitives/image_procs.py +++ b/potpyri/primitives/image_procs.py @@ -298,7 +298,7 @@ def align_images(reduced_files, paths, tel, binn, use_wcs=None, fieldcenter=None return(aligned_images, aligned_data) -def image_proc(image_data, tel, paths, skip_skysub=False, +def image_proc(image_data, tel, paths, skip_skysub=False, bkg_sub='local', fieldcenter=None, out_size=None, satellites=True, cosmic_ray=True, skip_fine_align=False, skip_gaia=None, fine_align_catalog='gaia', skip_external_astrometry=False, keep_all_astro=False, @@ -319,7 +319,12 @@ def image_proc(image_data, tel, paths, skip_skysub=False, paths : dict Paths dict from options.add_paths. skip_skysub : bool, optional - If True, skip sky subtraction. Default is False. + If True, skip sky subtraction. Deprecated; use ``bkg_sub='none'``. + Default is False. + bkg_sub : str, optional + Per-frame background mode for optical data: ``'local'`` (2D mesh), + ``'constant'`` (single median per frame), or ``'none'``. Default + is ``'local'``. fieldcenter : sequence, optional [ra, dec] for alignment. out_size : int, optional @@ -415,8 +420,10 @@ def image_proc(image_data, tel, paths, skip_skysub=False, # Bias subtraction, gain correction, flat correction, and flat fielding files = image_data['File'] + if skip_skysub: + bkg_sub = 'none' processed = tel.process_science(files, fil, amp, binn, paths, - mbias=mbias, mflat=mflat, mdark=mdark, skip_skysub=skip_skysub, log=log) + mbias=mbias, mflat=mflat, mdark=mdark, bkg_sub=bkg_sub, log=log) # Get filenames for output processed data reduced_files = [p.header['FILENAME'] for p in processed] diff --git a/potpyri/primitives/photometry.py b/potpyri/primitives/photometry.py index 154a04a..37ca044 100755 --- a/potpyri/primitives/photometry.py +++ b/potpyri/primitives/photometry.py @@ -247,7 +247,7 @@ def run_sextractor(img_file, log=None, sextractor_path=None): return(table) def extract_aperture_stats(img_data, img_mask, img_error, stars, - aperture_radius=10.0, log=None): + aperture_radius=10.0, fwhm_measure_radius=None, log=None): """Compute aperture flux and error for a star table; return table with added columns. Parameters @@ -261,7 +261,10 @@ def extract_aperture_stats(img_data, img_mask, img_error, stars, stars : astropy.table.Table Table with xcentroid, ycentroid (modified in place with refined centroids). aperture_radius : float, optional - Aperture radius in pixels. Default is 10.0. + Aperture radius in pixels for flux and centroid statistics. Default 10.0. + fwhm_measure_radius : float, optional + Smaller aperture for ``ApertureStats.fwhm`` only. If omitted, uses + ``aperture_radius``. FWHM from large apertures is biased high. log : ColoredLogger, optional Logger for progress. @@ -285,39 +288,49 @@ def extract_aperture_stats(img_data, img_mask, img_error, stars, # get_star_catalog() (or DAOStarFinder output with new column names). stars = _normalize_daofind_catalog(stars) - # Estimate a reasonable aperture radius and centroid for sources - fwhms=[] - for i,star in enumerate(stars): + # Refine centroids at a fixed aperture. Do not grow the aperture from + # ApertureStats FWHM values: that statistic scales with aperture radius + # and can diverge (e.g. 12 px -> 30 px -> ~34 px FWHM estimates). + # FWHM is measured in a compact aperture; flux uses a larger aperture. + fwhm_radius = fwhm_measure_radius if fwhm_measure_radius is not None else aperture_radius + fwhms = [] + for i, star in enumerate(stars): aper = CircularAperture((star['xcentroid'], star['ycentroid']), - aperture_radius) - aperstats = ApertureStats(img_data, aper, mask=img_mask, + fwhm_radius) + aperstats = ApertureStats(img_data, aper, mask=img_mask, error=img_error) fwhms.append(_apstats_float(aperstats, 'fwhm')) stars[i]['xcentroid'] = _apstats_float(aperstats, 'x_centroid', 'xcentroid') stars[i]['ycentroid'] = _apstats_float(aperstats, 'y_centroid', 'ycentroid') - if aperture_radius<2.5*np.nanmean(fwhms): - aperture_radius=2.5*np.nanmean(fwhms) - + flux_radius = max(aperture_radius, 2.5 * np.nanmedian(fwhms)) if log: - log.info(f'New aperture radius={aperture_radius}') + log.info( + f'Aperture radii: fwhm={fwhm_radius:.2f} px, flux={flux_radius:.2f} px' + ) else: - print(f'New aperture radius={aperture_radius}') + print( + f'Aperture radii: fwhm={fwhm_radius:.2f} px, flux={flux_radius:.2f} px' + ) for star in stars: aper = CircularAperture((star['xcentroid'], star['ycentroid']), - aperture_radius) + flux_radius) - aperstats = ApertureStats(img_data, aper, mask=img_mask, + aperstats = ApertureStats(img_data, aper, mask=img_mask, error=img_error) covx = np.maximum( _apstats_float(aperstats, 'covariance_xx', 'covar_sigx2'), 0.0) covy = np.maximum( _apstats_float(aperstats, 'covariance_yy', 'covar_sigy2'), 0.0) + aper_fwhm = CircularAperture((star['xcentroid'], star['ycentroid']), + fwhm_radius) + fwhm_stats = ApertureStats(img_data, aper_fwhm, mask=img_mask, + error=img_error) apertable.add_row([ - _apstats_float(aperstats, 'fwhm'), + _apstats_float(fwhm_stats, 'fwhm'), _apstats_float(aperstats, 'semimajor_axis', 'semimajor_sigma'), _apstats_float(aperstats, 'semiminor_axis', 'semiminor_sigma'), _apstats_float(aperstats, 'orientation'), @@ -332,6 +345,49 @@ def extract_aperture_stats(img_data, img_mask, img_error, stars, return(apertable) +def _radial_fwhm_from_array(data, center=None, max_radius=None): + """Estimate FWHM (pixels) from a 2D profile via a circular radial average.""" + data = np.asarray(data, dtype=float) + cy, cx = center if center is not None else (data.shape[0] // 2, data.shape[1] // 2) + if max_radius is None: + max_radius = min(cy, cx, data.shape[0] - cy - 1, data.shape[1] - cx - 1) + yy, xx = np.indices(data.shape) + r = np.sqrt((xx - cx) ** 2 + (yy - cy) ** 2) + peak = data[cy, cx] + if not np.isfinite(peak) or peak <= 0: + peak = np.nanmax(data) + if not np.isfinite(peak) or peak <= 0: + return np.nan + + rs, vs = [], [] + step = 0.1 + rb = 0.0 + while rb < max_radius: + m = (r >= rb) & (r < rb + step) + if m.any(): + rs.append(rb + step / 2.0) + vs.append(float(np.nanmean(data[m]))) + rb += step + if len(vs) < 3: + return np.nan + + rs = np.asarray(rs) + vs = np.asarray(vs) + half_level = peak / 2.0 + below = np.where(vs < half_level)[0] + if len(below) == 0: + return np.nan + i = int(below[0]) + if i == 0: + return 2.0 * rs[0] + r1, r2 = rs[i - 1], rs[i] + v1, v2 = vs[i - 1], vs[i] + if v1 == v2: + return 2.0 * r1 + r_half = r1 + (half_level - v1) * (r2 - r1) / (v2 - v1) + return 2.0 * r_half + + def generate_epsf(img_file, x, y, size=11, oversampling=2, maxiters=11, log=None): """Build an ePSF from cutouts at (x, y) positions in the image. @@ -381,7 +437,7 @@ def generate_epsf(img_file, x, y, size=11, oversampling=2, maxiters=11, return(epsf) def extract_fwhm_from_epsf(epsf, fwhm_init): - """Estimate FWHM from ePSF model (Moffat2D fit). + """Estimate FWHM from ePSF model (Moffat2D fit on the PSF core). Parameters ---------- @@ -395,27 +451,36 @@ def extract_fwhm_from_epsf(epsf, fwhm_init): float FWHM in pixels from fitted Moffat2D. """ - # Get the raw data for the FWHM and size in x and y - data = epsf.data - x = np.arange(data.shape[0]) - y = np.arange(data.shape[1]) - xx, yy = np.meshgrid(x, y) - - # Fit to Moffat2D model in astropy, initial guess is amplitude, - # centroid in x and y, core width of - # Moffat model and power index scaling of model - p_init = functional_models.Moffat2D(amplitude=0.5, x_0=data.shape[0]/2., - y_0=data.shape[1]/2., gamma=fwhm_init, alpha=1.) - - fit_p = fitting.LevMarLSQFitter() - - # Fit functional model to the data - p = fit_p(p_init, xx, yy, data) + data = np.asarray(epsf.data, dtype=float) + cy, cx = data.shape[0] // 2, data.shape[1] // 2 + # Fit only the core: a Moffat2D on the full ePSF array matches extended + # wings in large cutouts and can return FWHM >> the true PSF width. + half = min(cy, cx, max(5, min(int(2.5 * float(fwhm_init)), 12))) + sub = data[cy - half:cy + half + 1, cx - half:cx + half + 1] + xs = np.arange(sub.shape[0]) + ys = np.arange(sub.shape[1]) + xxs, yys = np.meshgrid(xs, ys) + + peak = float(np.nanmax(sub)) + p_init = functional_models.Moffat2D( + amplitude=peak if peak > 0 else 0.5, + x_0=half, + y_0=half, + gamma=max(float(fwhm_init), 1.0), + alpha=2.0, + ) - # Extract and round FWHM - fwhm = float('%.4f'%p.fwhm) + fit_p = fitting.LevMarLSQFitter() + with warnings.catch_warnings(): + warnings.simplefilter('ignore') + p = fit_p(p_init, xxs, yys, sub) - return(p.fwhm) + fwhm = float(p.fwhm) + if not np.isfinite(fwhm) or fwhm <= 0: + fwhm = _radial_fwhm_from_array(sub, center=(half, half), max_radius=half) + if not np.isfinite(fwhm) or fwhm <= 0: + fwhm = float(fwhm_init) + return float('%.4f' % fwhm) def run_photometry(img_file, epsf, fwhm, threshold, shape, stars): """Run PSF photometry and aperture photometry; append result tables to FITS. @@ -585,7 +650,9 @@ def get_star_catalog(img_data, img_mask, img_error, fwhm_init=5.0, # Extract the aperture stats from each star and append to the output catalog stats = extract_aperture_stats(img_data, img_mask, img_error, stars, - aperture_radius=2.5*fwhm_init, log=log) + aperture_radius=2.5*fwhm_init, + fwhm_measure_radius=max(4.0, float(fwhm_init)), + log=log) stars = hstack([stars, stats]) return(stars) @@ -714,9 +781,14 @@ def do_phot(img_file, metadata['NPSFSTAR']=len(bright) + fwhm_catalog = fwhm + # Instantiate EPSF - size=int(fwhm*fwhm_scale_psf) - if size%2==0: size=size+1 + size = int(max(fwhm_catalog, star_param['fwhm_init']) * fwhm_scale_psf) + size_cap = max(25, int(star_param['fwhm_init'] * fwhm_scale_psf * 2)) + size = min(size, size_cap) + if size % 2 == 0: + size = size + 1 if log: log.info(f'EPSF size will be {size} pixels') @@ -726,10 +798,19 @@ def do_phot(img_file, epsf = generate_epsf(img_file, bright['xcentroid'], bright['ycentroid'], size=size, oversampling=oversampling, maxiters=11, log=log) - fwhm = extract_fwhm_from_epsf(epsf, fwhm*oversampling) - # Scale by oversampling - fwhm = fwhm/oversampling - fwhm = float('%.4f'%fwhm) + fwhm_epsf = extract_fwhm_from_epsf(epsf, fwhm_catalog * oversampling) + fwhm_epsf = fwhm_epsf / oversampling + if (np.isfinite(fwhm_epsf) and fwhm_epsf > 0 + and 0.5 * fwhm_catalog <= fwhm_epsf <= 3.0 * fwhm_catalog): + fwhm = fwhm_epsf + else: + fwhm = fwhm_catalog + if log: + log.info( + f'Using catalog FWHM={fwhm:.4f} px ' + f'(ePSF fit {fwhm_epsf:.4f} px outside 0.5-3x range)' + ) + fwhm = float('%.4f' % fwhm) if log: log.info(f'FWHM={fwhm} pixels') diff --git a/potpyri/primitives/sort_files.py b/potpyri/primitives/sort_files.py index 754fc59..c9d7644 100755 --- a/potpyri/primitives/sort_files.py +++ b/potpyri/primitives/sort_files.py @@ -260,18 +260,32 @@ def handle_files(file_list, paths, tel, incl_bad=False, proc=None, elif os.path.exists(file_list): os.remove(file_list) - files = glob.glob(os.path.join(paths['raw'], tel.raw_format(proc)))+\ - glob.glob(os.path.join(paths['data'], tel.raw_format(proc)))+\ - glob.glob(os.path.join(paths['bad'], tel.raw_format(proc))) + discovery_msg = tel.format_raw_discovery_message(paths, proc) + if log: + log.info(discovery_msg) + else: + print(discovery_msg, flush=True) + files = tel.discover_raw_files(paths, proc)['files'] - if log: log.info('Sorting files and creating file lists.') - if len(files)!=0: - if log: log.info(f'{len(files)} files found.') - file_table = sort_files(files, file_list, tel, paths, incl_bad=incl_bad, + if log: + log.info('Sorting files and creating file lists.') + if len(files) != 0: + if log: + log.info(f'{len(files)} files found.') + file_table = sort_files(files, file_list, tel, paths, incl_bad=incl_bad, log=log) else: - if log: log.critical('No files found, please check data path and rerun.') + if log: + log.critical( + 'No raw input files found. See the discovery summary above for ' + 'directories searched and the glob pattern required.' + ) + else: + print( + 'No raw input files found. See discovery summary above.', + flush=True, + ) logging.shutdown() sys.exit(-1) diff --git a/potpyri/scripts/main_pipeline.py b/potpyri/scripts/main_pipeline.py index abb3111..e96a02b 100755 --- a/potpyri/scripts/main_pipeline.py +++ b/potpyri/scripts/main_pipeline.py @@ -40,6 +40,7 @@ def main_pipeline(instrument: str, phot_sn_max: float = None, fwhm_init: float = None, skip_skysub: bool = None, + bkg_sub: str = None, file_list_name: str = None, fieldcenter: list = None, out_size: int = None, @@ -94,7 +95,8 @@ def main_pipeline(instrument: str, ####################################################### log.info(f'Generating stack for {tar}') stack = image_procs.image_proc(file_table[file_table['TargType']==tar], tel, paths, - skip_skysub=skip_skysub, fieldcenter=fieldcenter, out_size=out_size, + skip_skysub=skip_skysub, bkg_sub=bkg_sub or 'local', + fieldcenter=fieldcenter, out_size=out_size, cosmic_ray=not skip_cr, skip_fine_align=skip_fa, fine_align_catalog=fac, skip_external_astrometry=skip_external_astrometry or False, diff --git a/potpyri/utils/options.py b/potpyri/utils/options.py index f7ba008..4bb81f4 100755 --- a/potpyri/utils/options.py +++ b/potpyri/utils/options.py @@ -40,10 +40,13 @@ def init_options(): help='''Option to only reduce a specific target. String used here must be contained within the target name in file headers. Optional parameter.''') - params.add_argument('--proc', - type=str, - default=True, - help='''Option to specify file processing for data ingestion.''') + params.add_argument('--proc', + type=str, + default=True, + help='''File-ingestion mode for raw discovery (instrument-specific glob). + Use ``fits`` (or ``.fits``) on any instrument for uncompressed ``*.fits``. + Other values depend on the instrument (e.g. GMOS ``dragons``, LRIS + ``archive``/``raw``, BINOSPEC ``proc``).''') params.add_argument('--include-bad','--incl-bad', default=False, action='store_true', @@ -71,7 +74,15 @@ def init_options(): params.add_argument('--skip-skysub', default=False, action='store_true', - help='''Do not perform sky subtraction during image processing.''') + help='''Do not perform sky subtraction during image processing. Equivalent + to ``--bkg-sub none``.''') + params.add_argument('--bkg-sub', + default='local', + choices=['local', 'constant', 'none'], + dest='bkg_sub', + help='''Per-frame background subtraction for optical data: ``local`` uses a + spatially varying 2D mesh (default); ``constant`` subtracts one + sigma-clipped median per frame; ``none`` skips subtraction.''') params.add_argument('--fieldcenter', default=None, nargs=2, @@ -130,7 +141,7 @@ def init_options(): def add_options(): """Parse command-line options and return normalized args. - Handles instrument aliases (e.g. BINO -> BINOSPEC, MMIR -> MMIRS). + Normalizes instrument aliases via :func:`potpyri.instruments.resolve_instrument_name`. Returns ------- @@ -140,11 +151,11 @@ def add_options(): params = init_options() args = params.parse_args() - # Handle/parse options - if 'BINO' in args.instrument: - args.instrument = 'BINOSPEC' - if 'MMIR' in args.instrument: - args.instrument = 'MMIRS' + if args.instrument is not None: + args.instrument = instruments.resolve_instrument_name(args.instrument) + + if args.skip_skysub: + args.bkg_sub = 'none' return(args) diff --git a/tests/test_background_subtraction.py b/tests/test_background_subtraction.py new file mode 100644 index 0000000..2a0b809 --- /dev/null +++ b/tests/test_background_subtraction.py @@ -0,0 +1,71 @@ +"""Tests for per-frame optical background subtraction modes.""" +import numpy as np +import pytest +import astropy.units as u +from astropy.nddata import CCDData + +from potpyri.instruments import instrument_getter +from potpyri.instruments.instrument import ( + BKG_SUB_CONSTANT, + BKG_SUB_LOCAL, + BKG_SUB_NONE, +) + + +def _flat_frame(level=1000.0, shape=(128, 128)): + data = np.full(shape, level, dtype=np.float64) + mask = np.zeros(shape, dtype=bool) + return CCDData(data, unit=u.electron, mask=mask) + + +def _gradient_frame(shape=(128, 128)): + y, x = np.indices(shape) + data = 500.0 + 2.0 * x + 3.0 * y + mask = np.zeros(shape, dtype=bool) + return CCDData(data, unit=u.electron, mask=mask) + + +def test_constant_background_subtracts_uniform_level(): + tel = instrument_getter('GMOS') + frame = _flat_frame(level=1000.0) + result = tel.apply_optical_background_subtraction( + frame, bkg_sub=BKG_SUB_CONSTANT) + assert result.header['BKGSUB'] == BKG_SUB_CONSTANT + assert result.header['SKYBKG'] == pytest.approx(1000.0) + finite = result.data[np.isfinite(result.data)] + assert np.nanmedian(finite) == pytest.approx(0.0, abs=1e-5) + + +def test_constant_background_preserves_spatial_gradient(): + tel = instrument_getter('GMOS') + frame = _gradient_frame() + result = tel.apply_optical_background_subtraction( + frame, bkg_sub=BKG_SUB_CONSTANT) + # Constant subtraction removes only the median level, not the gradient. + assert np.std(result.data) > np.std(frame.data) * 0.5 + + +def test_none_background_leaves_data_unchanged(): + tel = instrument_getter('GMOS') + frame = _flat_frame(level=1000.0) + result = tel.apply_optical_background_subtraction(frame, bkg_sub=BKG_SUB_NONE) + assert result.header['BKGSUB'] == BKG_SUB_NONE + assert result.header['SKYBKG'] == 0.0 + np.testing.assert_array_equal(result.data, frame.data) + + +def test_local_background_reduces_spatial_gradient_more_than_constant(): + tel = instrument_getter('GMOS') + frame = _gradient_frame(shape=(256, 256)) + constant = tel.apply_optical_background_subtraction( + frame, bkg_sub=BKG_SUB_CONSTANT) + local = tel.apply_optical_background_subtraction(frame, bkg_sub=BKG_SUB_LOCAL) + assert local.header['BKGSUB'] == BKG_SUB_LOCAL + assert np.nanstd(local.data) < np.nanstd(constant.data) + + +def test_invalid_bkg_sub_raises(): + tel = instrument_getter('GMOS') + frame = _flat_frame() + with pytest.raises(ValueError, match='bkg_sub must be one of'): + tel.apply_optical_background_subtraction(frame, bkg_sub='polynomial') diff --git a/tests/test_binspec_frb_stack_photometry.py b/tests/test_binspec_frb_stack_photometry.py deleted file mode 100644 index 0af2f42..0000000 --- a/tests/test_binspec_frb_stack_photometry.py +++ /dev/null @@ -1,113 +0,0 @@ -"""Integration test: photometry + zeropoint on FRB BINOSPEC stack. - -Motivation: photutils 3.x ``DAOStarFinder`` uses ``x_centroid``/``y_centroid`` -column names; POTPyRI historically expected ``xcentroid``/``ycentroid``, which -raised ``KeyError: 'xcentroid'`` inside ``extract_aperture_stats``. - -The stack must be provided locally (large file, not in git). Set:: - - export POTPYRI_FRB_STACK_TEST=/absolute/path/to/FRB20260213A.r.ut260322.2.11.stk.fits - -or use the default ``$HOME/FRB20260213A.r.ut260322.2.11.stk.fits``. The test skips -if that path is missing. - -Network (Vizier / Pan-STARRS1) is required for the zeropoint step; unreachable -catalog service skips the test similarly to ``test_absphot``. -""" -from __future__ import annotations - -import os -import shutil - -import numpy as np -import pytest -import requests -from astropy.io import fits - -from potpyri.instruments import instrument_getter -from potpyri.primitives import absphot, photometry -from potpyri.utils import logger, options - -_DEFAULT_FRB_STACK = os.path.join( - os.path.expanduser('~'), 'FRB20260213A.r.ut260322.2.11.stk.fits' -) - - -def _frb_stack_path() -> str: - return os.environ.get('POTPYRI_FRB_STACK_TEST', _DEFAULT_FRB_STACK) - - -def _strip_optional_hdu(hdul): - for key in ('APPPHOT', 'PSFPHOT', 'PSFSTARS', 'RESIDUAL', 'PSF'): - if key in hdul: - del hdul[key] - - -def _sanitize_error(hdul): - err = hdul['ERROR'].data - med = np.nanmedian(err) - err[np.isnan(err)] = med - err[err < 0.0] = med - err[np.isinf(err)] = np.max(hdul['SCI'].data) - - -def _prepare_stack_copy(src: str, dest: str) -> None: - shutil.copy2(src, dest) - with fits.open(dest, mode='update') as hdul: - _strip_optional_hdu(hdul) - _sanitize_error(hdul) - hdul.flush() - - -@pytest.mark.integration -def test_binspec_frb_stack_photometry_and_absphot(tmp_path): - """BINOSPEC FRB stack completes photloop and absphot.find_zeropoint.""" - stack_src = _frb_stack_path() - if not os.path.isfile(stack_src): - pytest.skip( - f'Missing FRB stack at {stack_src!r}; copy the FITS locally or set ' - 'POTPYRI_FRB_STACK_TEST to its absolute path.' - ) - - stack = str(tmp_path / 'frb_stack.fits') - _prepare_stack_copy(stack_src, stack) - - tel = instrument_getter('BINOSPEC') - paths = options.add_paths(str(tmp_path), 'files.txt', tel) - log = logger.get_log(paths['log']) - - try: - try: - photometry.photloop( - stack, - phot_sn_min=3.0, - phot_sn_max=20.0, - fwhm_init=5.0, - log=log, - ) - except photometry.PhotometryError as exc: - pytest.fail(f'photometry.photloop failed: {exc}') - - with fits.open(stack) as hdul: - names = [h.name for h in hdul] - assert 'APPPHOT' in names, f'expected APPPHOT HDU, got {names!r}' - assert 'PSFPHOT' in names - assert 'PSFSTARS' in names - assert 'PSF' in names - n_app = len(hdul['APPPHOT'].data) - assert n_app > 100, f'expected substantial catalog, got NOBJECT={n_app}' - - try: - ok = absphot.find_zeropoint(stack, tel, log=log) - except (requests.exceptions.ConnectionError, OSError) as exc: - pytest.skip(f'Catalog / network unavailable for zeropoint: {exc}') - - assert ok is True, 'absphot.find_zeropoint returned False (see log)' - - with fits.open(stack) as hdul: - for hdr in (hdul['PRIMARY'].header, hdul['SCI'].header): - assert 'ZPTMAG' in hdr, 'ZPTMAG not written to PRIMARY/SCI' - assert np.isfinite(hdr['ZPTMAG']) - assert hdr.get('ZPTNSTAR', 0) >= 5 - finally: - log.close() diff --git a/tests/test_instruments.py b/tests/test_instruments.py index ca0a772..6a7f0b0 100644 --- a/tests/test_instruments.py +++ b/tests/test_instruments.py @@ -8,7 +8,13 @@ from astropy.nddata import CCDData from astropy import units as u -from potpyri.instruments import instrument_getter +from potpyri.instruments import ( + UnknownInstrumentError, + __all__ as INSTRUMENT_NAMES, + instrument_getter, + resolve_instrument_name, + supported_instruments, +) from potpyri.instruments.GMOS import GMOS from potpyri.instruments.instrument import ( Instrument, @@ -32,15 +38,31 @@ def test_fix_deprecated_wcs_header_cards_radecsys_and_mjd(): assert 'MJD-OBS' not in h +def test_resolve_instrument_name_aliases(): + """resolve_instrument_name maps aliases and canonical names.""" + assert resolve_instrument_name('gmos') == 'GMOS' + assert resolve_instrument_name('BINO') == 'BINOSPEC' + assert resolve_instrument_name('MMIR') == 'MMIRS' + assert resolve_instrument_name('BINOSPEC') == 'BINOSPEC' + + +def test_supported_instruments_matches_all(): + assert set(supported_instruments()) == set(INSTRUMENT_NAMES) + + def test_instrument_getter_unsupported_raises(): - """instrument_getter raises when instrument is not supported and log is None.""" - with pytest.raises(Exception, match="not supported"): + """instrument_getter raises UnknownInstrumentError when log is None.""" + with pytest.raises(UnknownInstrumentError, match="not supported"): instrument_getter("UNKNOWN_INSTRUMENT", log=None) +def test_instrument_getter_invalid_name_type(): + with pytest.raises(TypeError): + instrument_getter("", log=None) + + def test_instrument_getter_unsupported_with_log(tmp_path): """instrument_getter with log calls log.error and returns None for unsupported name.""" - from potpyri.instruments import __init__ as instruments_init tel = instrument_getter("GMOS") paths = options.add_paths(str(tmp_path), "files.txt", tel) log = logger.get_log(paths["log"]) @@ -154,13 +176,56 @@ def test_format_datasec(): def test_raw_format(): - """Base raw_format and GMOS raw_format.""" + """raw_format and --proc fits override across instruments.""" base = Instrument() assert base.raw_format(True) == "sci_img_*.fits" assert base.raw_format(False) == "sci_img*[!proc].fits" + assert base.raw_format("fits") == "*.fits" gmos = GMOS() assert gmos.raw_format("dragons") == "*.fits" assert gmos.raw_format("other") == "*.fits.bz2" + assert gmos.raw_format("fits") == "*.fits" + f2 = instrument_getter("F2") + assert f2.raw_format(True) == "*.fits.bz2" + assert f2.raw_format("fits") == "*.fits" + assert f2.raw_format(".fits") == "*.fits" + assert f2.raw_format("uncompressed") == "*.fits" + mosfire = instrument_getter("MOSFIRE") + assert mosfire.raw_format(None) == "*.fits.gz" + assert mosfire.raw_format("fits") == "*.fits" + lris = instrument_getter("LRIS") + assert lris.raw_format("archive") == "*.fits*" + assert lris.raw_format("fits") == "*.fits" + + +def test_discover_raw_files_and_message(tmp_path): + """discover_raw_files reports search paths and glob matches.""" + raw_dir = tmp_path / 'raw' + data_dir = tmp_path / 'data' + bad_dir = tmp_path / 'bad' + for d in (raw_dir, data_dir, bad_dir): + d.mkdir() + (raw_dir / 'science.fits').write_bytes(b'\x00') + (data_dir / 'other.fits.bz2').write_bytes(b'\x00') + + tel = instrument_getter("F2") + paths = {'raw': str(raw_dir), 'data': str(data_dir), 'bad': str(bad_dir)} + + info = tel.discover_raw_files(paths, proc='fits') + assert info['pattern'] == '*.fits' + assert len(info['files']) == 1 + assert info['files'][0].endswith('science.fits') + + msg = tel.format_raw_discovery_message(paths, proc='fits') + assert 'Glob pattern' in msg + assert str(raw_dir) in msg + assert 'science.fits' not in msg or '1 file(s)' in msg + assert 'Search directories' in msg + + info_bz2 = tel.discover_raw_files(paths, proc=True) + assert info_bz2['pattern'] == '*.fits.bz2' + assert len(info_bz2['files']) == 1 + assert info_bz2['files'][0].endswith('other.fits.bz2') def test_get_stk_name_get_sci_name_get_bkg_name(): diff --git a/tests/test_options.py b/tests/test_options.py index 540dd31..65ae10e 100644 --- a/tests/test_options.py +++ b/tests/test_options.py @@ -43,3 +43,13 @@ def test_add_paths_creates_dirs(tmp_path): assert os.path.exists(paths["raw"]) assert os.path.exists(paths["red"]) assert os.path.exists(paths["log"]) + + +def test_skip_skysub_sets_bkg_sub_none(monkeypatch): + """--skip-skysub overrides --bkg-sub to none.""" + monkeypatch.setattr( + 'sys.argv', + ['main_pipeline', 'GMOS', '/tmp/data', '--skip-skysub', '--bkg-sub', 'constant'], + ) + args = options.add_options() + assert args.bkg_sub == 'none' diff --git a/tests/test_photometry.py b/tests/test_photometry.py index c858b72..ee7480d 100644 --- a/tests/test_photometry.py +++ b/tests/test_photometry.py @@ -221,3 +221,42 @@ def test_extract_fwhm_from_epsf(): fwhm = photometry.extract_fwhm_from_epsf(epsf, fwhm_init=3.0) assert np.isfinite(fwhm) assert fwhm > 0 + + +def test_extract_fwhm_from_epsf_wide_array_uses_core(): + """A broad flat ePSF array must not return FWHM ~ array size.""" + size = 45 + y, x = np.ogrid[-size // 2:size // 2 + 1, -size // 2:size // 2 + 1] + data = np.exp(-(x * x + y * y) / (2 * 2.0**2)).astype(float) + data += 0.02 # extended pedestal like stacked ePSF wings + epsf = SimpleNamespace(data=data) + fwhm = photometry.extract_fwhm_from_epsf(epsf, fwhm_init=5.0) + assert np.isfinite(fwhm) + assert fwhm < 15.0 + + +def test_extract_aperture_stats_fwhm_uses_compact_radius(): + """FWHM must be measured at fwhm_measure_radius, not the flux aperture.""" + stars = Table( + { + 'xcentroid': Column([32.0]), + 'ycentroid': Column([32.0]), + 'peak': Column([100.0]), + 'flux': Column([50.0]), + }, + ) + img = np.ones((64, 64), dtype=float) * 50.0 + y, x = np.ogrid[:64, :64] + img += 500.0 * np.exp(-((x - 32) ** 2 + (y - 32) ** 2) / (2 * 2.0**2)) + mask = np.zeros((64, 64), dtype=bool) + err = np.ones((64, 64), dtype=float) * 5.0 + + stats_small = photometry.extract_aperture_stats( + img, mask, err, stars, aperture_radius=20.0, + fwhm_measure_radius=5.0, log=None, + ) + stats_large = photometry.extract_aperture_stats( + img, mask, err, stars, aperture_radius=20.0, + fwhm_measure_radius=20.0, log=None, + ) + assert stats_small['fwhm'][0] < stats_large['fwhm'][0] diff --git a/tests/test_sort.py b/tests/test_sort.py index 21d471d..e078751 100644 --- a/tests/test_sort.py +++ b/tests/test_sort.py @@ -134,3 +134,28 @@ def test_is_bias_gmos(): assert sort_files.is_bias(hdr, tel) == True hdr['OBSTYPE'] = 'object' assert sort_files.is_bias(hdr, tel) == False + + +def test_handle_files_logs_discovery_when_no_files(tmp_path, capsys): + """handle_files prints search paths and glob pattern before exiting on empty input.""" + tel = instrument_getter('F2') + raw_dir = tmp_path / 'raw' + data_dir = tmp_path + bad_dir = tmp_path / 'bad' + for d in (raw_dir, bad_dir): + d.mkdir() + paths = { + 'raw': str(raw_dir), + 'data': str(data_dir), + 'bad': str(bad_dir), + 'filelist': str(tmp_path / 'file_list.txt'), + } + with pytest.raises(SystemExit): + sort_files.handle_files( + paths['filelist'], paths, tel, proc='fits', log=None, + ) + out = capsys.readouterr().out + assert 'Glob pattern' in out + assert '*.fits' in out + assert str(raw_dir) in out + assert 'No files matched' in out