diff --git a/.coverage b/.coverage new file mode 100644 index 00000000..397322bf Binary files /dev/null and b/.coverage differ diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index afa1e085..4b1e469a 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -11,9 +11,9 @@ If you encounter a bug, have a question about usage, or have a feature request: ## Contributing code or documentation -1. **Fork the repository** and create a branch from the default branch. -2. **Make your changes** and add or update tests if applicable. Run the test suite with `pytest tests` (use `pytest tests -m "not integration"` for offline runs). -3. **Open a pull request** against the main repository. Describe your changes clearly and reference any related issues. +1. **Fork the repository** and create a branch (often from `dev` for ongoing development; the default branch may be `main`). +2. **Make your changes** and add or update tests if applicable. Run the test suite with `pytest tests` (use `pytest tests -m "not integration"` for offline runs without fixture downloads). +3. **Open a pull request** against the main repository (or against `dev` when that branch is used for integration). Describe your changes clearly and reference any related issues. 4. **Code style:** Follow the existing style in the codebase. The project uses standard Python packaging and type hints where appropriate. Instrument-specific changes (e.g. new instruments or header/sorting logic) should include a brief justification and, if possible, a note in the PR description on how the change was tested. diff --git a/README.md b/README.md index 8b4d5b50..336591c8 100755 --- a/README.md +++ b/README.md @@ -76,8 +76,12 @@ From the project root (with test deps installed: `pip install -e ".[test]"`): pytest tests ``` -- **Offline (no network):** - `pytest tests -m "not integration"` +- **Unit tests only (no network):** + `pytest tests -m "not integration"` + This skips the five integration tests that download fixtures from external storage. CI runs this by default. +- **All tests (including integration):** + `pytest tests` + Requires network access for fixture downloads. - **Verbose:** `pytest tests -v` @@ -227,7 +231,9 @@ POTPyRI can stack data from multiple nights, although it will not use unique cal ## Documentation -Full API documentation (generated from docstrings) is available at [https://CIERA-Transients.github.io/POTPyRI/](https://CIERA-Transients.github.io/POTPyRI/). To build the docs locally, install with `pip install -e ".[docs]"` and run `cd docs && make html` (or use the Sphinx `make.bat` on Windows). +- **Online:** Full API documentation (generated from docstrings) is at [https://ciera-transients.github.io/POTPyRI/](https://ciera-transients.github.io/POTPyRI/). +- **Local build:** Install with `pip install -e ".[docs]"` and run `cd docs && make html` (on Windows use `make.bat html`). Output is in `docs/build/html/`. +- **Deployment:** Docs are built and deployed to GitHub Pages on pushes to `main`; see `docs/DEPLOY.md` for setup details. ## Contributing and support @@ -237,7 +243,7 @@ We welcome contributions and feedback. See [CONTRIBUTING.md](CONTRIBUTING.md) fo - How to **contribute** code or documentation (pull requests) - Where to **seek support** (issues or developer contact) -If you encounter an error in the pipeline, have a special data setup that does not run through the pipeline, or wish to add an instrument, please open an issue on [GitHub](https://github.com/CIERA-Transients/POTPyRI/issues) or contact the developers at `ckilpatrick@northwestern.edu`. +Development often uses the `dev` branch; pull requests may target `main` or `dev` as noted in the repository. If you encounter an error in the pipeline, have a special data setup that does not run through the pipeline, or wish to add an instrument, please open an issue on [GitHub](https://github.com/CIERA-Transients/POTPyRI/issues) or contact the developers at `ckilpatrick@northwestern.edu`. ## License diff --git a/docs/source/index.rst b/docs/source/index.rst index 72c22363..69106792 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -7,6 +7,8 @@ It provides instrument-specific reduction workflows for bias, dark, and flat calibration; image alignment and stacking; WCS solving with astrometry.net; photometry with photutils; and flux calibration with astroquery. +The package includes nine instrument modules (BINOSPEC, DEIMOS, F2, FOURSTAR, GMOS, IMACS, LRIS, MMIRS, MOSFIRE), primitives for calibration, image processing, file sorting, WCS solving, photometry, and absolute photometry, plus CLI scripts for the main pipeline and archive data download. For installation, usage, and contributing, see the repository `README `_ and `CONTRIBUTING.md `_. + Contents ======== diff --git a/docs/source/installation.rst b/docs/source/installation.rst index f2da5374..fe2ece20 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -1,7 +1,7 @@ Installation ============ -POTPyRI requires **Python 3.11 or later**. +POTPyRI requires **Python 3.11 or later**. For full platform notes (conda, Apple Silicon, etc.) and the ``environment.yml``-based workflow, see the main **README** in the repository root. From PyPI --------- @@ -43,7 +43,15 @@ The ``-e`` (editable) flag uses the current directory as the live source, so cha Non-Python dependencies ------------------------ -The pipeline uses **astrometry.net** and **Source Extractor** (SExtractor). Install them via your system package manager or conda (see the main README). Index files for astrometry.net can be installed with the ``download_anet_index`` script after astrometry.net is on your path. +The pipeline uses **astrometry.net** and **Source Extractor** (SExtractor). Install them via your system package manager or conda (see the main README for a platform table). Index files for astrometry.net (~59 GB) can be installed with the ``download_anet_index`` script after astrometry.net is on your path. + +Testing +------- + +From the project root with ``.[test]`` installed: + +- **Unit tests only (no network):** ``pytest tests -m "not integration"`` — skips the five integration tests that download fixtures. This is what CI runs. +- **All tests:** ``pytest tests`` — requires network access for fixture downloads. Building the documentation --------------------------- @@ -54,4 +62,4 @@ With ``.[docs]`` installed, from the project root: cd docs && make html -Output is in ``docs/build/html/``. On Windows use ``make.bat html``. +Output is in ``docs/build/html/``. On Windows use ``make.bat html``. The live API docs are published at https://ciera-transients.github.io/POTPyRI/; see ``docs/DEPLOY.md`` for deployment setup. diff --git a/potpyri/instruments/instrument.py b/potpyri/instruments/instrument.py index 074ced65..3d7c0e97 100755 --- a/potpyri/instruments/instrument.py +++ b/potpyri/instruments/instrument.py @@ -25,6 +25,49 @@ from astropy.stats import SigmaClip from astropy.time import Time +# 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 = ( + 'CTYPE1', 'CTYPE2', 'CRVAL1', 'CRVAL2', 'CRPIX1', 'CRPIX2', + 'CD1_1', 'CD1_2', 'CD2_1', 'CD2_2', 'CDELT1', 'CDELT2', 'CROTA2', + 'RADECSYS', 'RADESYS', 'EQUINOX', 'RA', 'DEC', 'LONPOLE', 'LATPOLE', + 'CUNIT1', 'CUNIT2', 'MJD-OBS', + 'LTM1_1', 'LTM1_2', 'LTM2_1', 'LTM2_2', 'LTV1', 'LTV2', + 'DTV1', 'DTV2', 'DTM1_1', 'DTM2_2', +) + + +def _sanitize_calibration_header(header): + """Remove WCS and coordinate keywords from a header so it can be written safely. + + Calibration frames (bias, dark, flat, sky) inherit headers from input frames + that may contain invalid or ill-conditioned WCS (e.g. DEC=-100). Removing + these keywords ensures saved calibration files do not cause + InvalidTransformError when later loaded by code that parses WCS. + Modifies header in place. + """ + for key in _CAL_HEADER_WCS_KEYS: + if key in header: + del header[key] + # Remove PV* distortion keywords (can be ill-conditioned) + to_del = [k for k in header.keys() if k.startswith('PV')] + for k in to_del: + del header[k] + + +def _read_calibration_ccd(path, unit, hdu_index=0): + """Read a calibration FITS into CCDData without parsing WCS. + + Master bias/dark/flat frames often have WCS keywords that are ill-conditioned + or copied from templates. Parsing them can raise InvalidTransformError and + log 'Ill-conditioned coordinate transformation parameter'. Calibration + frames do not need WCS, so we load data and header explicitly with wcs=None. + """ + with fits.open(path) as hdu_list: + hdu = hdu_list[hdu_index] + return CCDData(hdu.data, meta=hdu.header.copy(), unit=unit, wcs=None) + + class Instrument(object): """Base class for all POTPyRI instruments. @@ -384,10 +427,9 @@ def load_bias(self, paths, amp, binn): """ bias = self.get_mbias_name(paths, amp, binn) if os.path.exists(bias): - mbias = CCDData.read(bias) + mbias = _read_calibration_ccd(bias, u.electron, hdu_index=0) elif os.path.exists(bias+'.fz'): - hdu = fits.open(bias+'.fz') - mbias = CCDData(hdu[1].data, header=hdu[1].header, unit=u.electron) + mbias = _read_calibration_ccd(bias+'.fz', u.electron, hdu_index=1) else: raise Exception(f'Could not find bias: {bias}') return(mbias) @@ -416,10 +458,9 @@ def load_dark(self, paths, amp, binn): """ dark = self.get_mdark_name(paths, amp, binn) if os.path.exists(dark): - mdark = CCDData.read(dark) + mdark = _read_calibration_ccd(dark, u.electron, hdu_index=0) elif os.path.exists(dark+'.fz'): - hdu = fits.open(dark+'.fz') - mdark = CCDData(hdu[1].data, header=hdu[1].header, unit=u.electron) + mdark = _read_calibration_ccd(dark+'.fz', u.electron, hdu_index=1) else: raise Exception(f'Could not find dark: {dark}') return(mdark) @@ -450,10 +491,9 @@ def load_flat(self, paths, fil, amp, binn): """ flat = self.get_mflat_name(paths, fil, amp, binn) if os.path.exists(flat): - mflat = CCDData.read(flat) + mflat = _read_calibration_ccd(flat, u.dimensionless_unscaled, hdu_index=0) elif os.path.exists(flat+'.fz'): - hdu = fits.open(flat+'.fz') - mflat = CCDData(hdu[1].data, header=hdu[1].header, unit=u.electron) + mflat = _read_calibration_ccd(flat+'.fz', u.dimensionless_unscaled, hdu_index=1) else: raise Exception(f'Could not find flat: {flat}') return(mflat) @@ -484,10 +524,9 @@ def load_sky(self, paths, fil, amp, binn): """ sky = self.get_msky_name(paths, fil, amp, binn) if os.path.exists(sky): - msky = CCDData.read(sky) + msky = _read_calibration_ccd(sky, u.dimensionless_unscaled, hdu_index=0) elif os.path.exists(sky+'.fz'): - hdu = fits.open(sky+'.fz') - msky = CCDData(hdu[1].data, header=hdu[1].header, unit=u.electron) + msky = _read_calibration_ccd(sky+'.fz', u.dimensionless_unscaled, hdu_index=1) else: raise Exception(f'Could not find sky frame: {sky}') return(msky) @@ -615,8 +654,9 @@ def create_bias(self, bias_list, amp, binn, paths, 'Version of telescope parameter file used.') bias_filename = self.get_mbias_name(paths, amp, binn) + _sanitize_calibration_header(mbias.header) mbias.write(bias_filename, overwrite=True, output_verify='silentfix') - log.info(f'Master bias written to {bias_filename}') + if log: log.info(f'Master bias written to {bias_filename}') return @@ -686,6 +726,7 @@ def create_dark(self, dark_list, amp, binn, paths, mbias=None, log=None): darkname = self.get_mdark_name(paths, amp, binn) if log: log.info(f'Writing master dark to {darkname}') + _sanitize_calibration_header(mdark.header) mdark.write(darkname, overwrite=True, output_verify='silentfix') return @@ -798,8 +839,9 @@ def create_flat(self, flat_list, fil, amp, binn, paths, mbias=None, 'Version of telescope parameter file used.') flat_filename = self.get_mflat_name(paths, fil, amp, binn) + _sanitize_calibration_header(mflat.header) mflat.write(flat_filename, overwrite=True, output_verify='silentfix') - log.info(f'Master flat written to {flat_filename}') + if log: log.info(f'Master flat written to {flat_filename}') return @@ -881,8 +923,9 @@ def create_sky(self, sky_list, fil, amp, binn, paths, log=None, **kwargs): 'Version of telescope parameter file used.') sky_filename = self.get_msky_name(paths, fil, amp, binn) + _sanitize_calibration_header(msky.header) msky.write(sky_filename, overwrite=True, output_verify='silentfix') - log.info(f'Master sky written to {sky_filename}') + if log: log.info(f'Master sky written to {sky_filename}') return @@ -1121,10 +1164,12 @@ def process_science(self, sci_list, fil, amp, binn, paths, mbias=None, for i,frame in enumerate(processed): mean, med, stddev = sigma_clipped_stats(frame.data) - frame_sky = sky_frame.multiply(med, + # Scale normalized sky to same units as science (electrons) so subtract is valid + science_unit = frame.unit if frame.unit is not None else u.electron + frame_sky = sky_frame.multiply(med * science_unit, propagate_uncertainties=True, handle_meta='first_found') - processed[i] = frame.subtract(frame_sky, + processed[i] = frame.subtract(frame_sky, propagate_uncertainties=True, handle_meta='first_found') processed[i].header['SKYBKG']=med processed[i].header['SATURATE']-=med diff --git a/potpyri/primitives/absphot.py b/potpyri/primitives/absphot.py index fdf55801..6d5af552 100755 --- a/potpyri/primitives/absphot.py +++ b/potpyri/primitives/absphot.py @@ -463,6 +463,9 @@ def convert_filter_name(self, filt): return 'z' if filt=='Y': return 'J' + # K, Ks, Kspec all use 2MASS K-band for calibration + if filt in ('K', 'Ks', 'Kspec'): + return 'K' else: return filt diff --git a/potpyri/utils/utilities.py b/potpyri/utils/utilities.py index 7bf6af9e..2602cadc 100755 --- a/potpyri/utils/utilities.py +++ b/potpyri/utils/utilities.py @@ -71,19 +71,24 @@ def find_catalog(catalog, fil, coord_ra, coord_dec): # If these catalogs are to be updated in the future, select mag columns that correspond to the # PSF mags. if catalog.upper() == 'SDSS': - if fil.lower() not in ['u','g','r','i','z']: return(catalog_ID, ra, dec, mag, err) + if fil.lower() not in ['u','g','r','i','z']: return(catalog, catalog_ID, ra, dec, mag, err) catalog_ID, ra, dec, mag, err = 'V/154', 'RA_ICRS', 'DE_ICRS', fil.lower()+'mag', 'e_'+fil.lower()+'mag' elif catalog.upper() == '2MASS': - if fil.upper() not in ['J','H','K']: return(catalog_ID, ra, dec, mag, err) - catalog_ID, ra, dec, mag, err = 'II/246', 'RAJ2000', 'DEJ2000', fil.upper()+'mag', 'e_'+fil.upper()+'mag' + # K, Ks, Kspec all use 2MASS K-band columns (Kmag, e_Kmag) + fil_2mass = fil.upper() + if fil_2mass in ('KS', 'KSPEC'): + fil_2mass = 'K' + if fil_2mass not in ['J', 'H', 'K']: + return(catalog, catalog_ID, ra, dec, mag, err) + catalog_ID, ra, dec, mag, err = 'II/246', 'RAJ2000', 'DEJ2000', fil_2mass+'mag', 'e_'+fil_2mass+'mag' elif catalog.upper() == 'UKIRT': - if fil.upper() not in ['Y','J','H','K']: return(catalog_ID, ra, dec, mag, err) + if fil.upper() not in ['Y','J','H','K']: return(catalog, catalog_ID, ra, dec, mag, err) catalog_ID, ra, dec, mag, err = 'II/319', 'ra', 'dec', fil.upper() + 'mag', 'e_'+fil.upper()+'mag' elif catalog.upper() == 'PS1': - if fil.lower() not in ['g','r','i','z','y']: return(catalog_ID, ra, dec, mag, err) + if fil.lower() not in ['g','r','i','z','y']: return(catalog, catalog_ID, ra, dec, mag, err) catalog_ID, ra, dec, mag, err = 'II/349', 'RAJ2000', 'DEJ2000', fil.lower()+'mag', 'e_'+fil.lower()+'mag' elif catalog.upper() == 'SKYMAPPER': - if fil.lower() not in ['u','v','g','r','i','z']: return(catalog_ID, ra, dec, mag, err) + if fil.lower() not in ['u','v','g','r','i','z']: return(catalog, catalog_ID, ra, dec, mag, err) catalog_ID, ra, dec, mag, err = 'II/379/smssdr4', 'RAICRS', 'DEICRS', fil.lower()+'PSF', 'e_'+fil.lower()+'PSF' return(catalog, catalog_ID, ra, dec, mag, err) diff --git a/tests/test_absphot.py b/tests/test_absphot.py index a7f6905c..81323d3e 100644 --- a/tests/test_absphot.py +++ b/tests/test_absphot.py @@ -91,6 +91,10 @@ def test_convert_filter_name(): assert cal.convert_filter_name('Y') == 'J' assert cal.convert_filter_name('V') == 'g' assert cal.convert_filter_name('I') == 'i' + # K, Ks, Kspec all map to 2MASS K-band + assert cal.convert_filter_name('K') == 'K' + assert cal.convert_filter_name('Ks') == 'K' + assert cal.convert_filter_name('Kspec') == 'K' def test_get_minmag(): diff --git a/tests/test_image_procs.py b/tests/test_image_procs.py index 40f0475c..ffbe1633 100644 --- a/tests/test_image_procs.py +++ b/tests/test_image_procs.py @@ -103,3 +103,22 @@ def test_add_stack_mask(tmp_path): assert np.any(stack[1].data >= 0) finally: log.close() + + +def test_create_error(tmp_path): + """create_error returns error HDU from science path, mask HDU, and rdnoise.""" + science_path = os.path.join(tmp_path, "sci.fits") + data = np.ones((16, 16), dtype=np.float32) * 100.0 + hdu_sci = fits.PrimaryHDU(data) + hdu_sci.header["SATURATE"] = 50000.0 + hdu_sci.writeto(science_path, overwrite=True) + + mask_data = fits.ImageHDU(np.zeros((16, 16), dtype=np.uint8)) + rdnoise = 4.0 + + err_hdu = image_procs.create_error(science_path, mask_data, rdnoise) + assert err_hdu is not None + assert err_hdu.data.shape == (16, 16) + assert err_hdu.header.get("BUNIT") == "ELECTRONS" + assert np.all(err_hdu.data > 0) + assert np.all(np.isfinite(err_hdu.data)) diff --git a/tests/test_instruments.py b/tests/test_instruments.py new file mode 100644 index 00000000..c2f7b7c9 --- /dev/null +++ b/tests/test_instruments.py @@ -0,0 +1,332 @@ +"""Unit tests for instrument base class methods and instrument_getter.""" +import os + +import numpy as np +import pytest +from astropy.io import fits +from astropy.table import Table +from astropy.nddata import CCDData +from astropy import units as u + +from potpyri.instruments import instrument_getter +from potpyri.instruments.GMOS import GMOS +from potpyri.instruments.instrument import ( + Instrument, + _sanitize_calibration_header, + _read_calibration_ccd, +) +from potpyri.utils import options +from potpyri.utils import logger + + +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("UNKNOWN_INSTRUMENT", 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"]) + try: + # When log is provided, getter still raises (code path: log.error then no return) + # Actually re-reading the code: if log, it calls log.error and then falls through + # and tel stays None, so it returns None. So we get None, not an exception. + result = instrument_getter("UNSUPPORTED", log=log) + assert result is None + finally: + log.close() + + +def test_match_type_keywords(): + """Instrument.match_type_keywords returns mask for Type column.""" + tel = GMOS() + file_table = Table({"Type": ["SCIENCE", "BIAS", "SCIENCE", "FLAT"]}) + mask = tel.match_type_keywords("SCIENCE,FLAT", file_table) + assert np.array_equal(mask, [True, False, True, True]) + + +def test_needs_sky_subtraction(): + """GMOS needs_sky_subtraction True for z-band, False otherwise.""" + tel = GMOS() + assert tel.needs_sky_subtraction("z") is True + assert tel.needs_sky_subtraction("r") is False + assert tel.needs_sky_subtraction("Z") is True + + +def test_get_pixscale(): + """get_pixscale returns instrument pixel scale.""" + tel = GMOS() + assert tel.get_pixscale() == 0.0803 + + +def test_get_rdnoise_get_gain(): + """get_rdnoise and get_gain use header if present else default.""" + tel = GMOS() + hdr = fits.Header() + assert tel.get_rdnoise(hdr) == 4.14 + assert tel.get_gain(hdr) == 1.63 + hdr["RDNOISE"] = 5.0 + hdr["GAIN"] = 2.0 + assert tel.get_rdnoise(hdr) == 5.0 + assert tel.get_gain(hdr) == 2.0 + + +def test_get_target_get_filter_get_exptime(): + """Header getters for target, filter, exptime.""" + tel = GMOS() + hdr = fits.Header({"OBJECT": " NGC1234 ", "FILTER2": " r_G0326 ", "EXPTIME": 60.0}) + assert tel.get_target(hdr) == "NGC1234" + assert tel.get_filter(hdr) == "r" + assert tel.get_exptime(hdr) == 60.0 + + +def test_get_ampl_get_binning(): + """GMOS get_ampl and get_binning from header.""" + tel = GMOS() + hdr = fits.Header({"NCCDS": "1", "CCDSUM": "2 2"}) + assert tel.get_ampl(hdr) == "4" + hdr["NCCDS"] = "2" + assert tel.get_ampl(hdr) == "12" + assert tel.get_binning(hdr) == "22" + + +def test_get_out_size(): + """get_out_size scales out_size by binning.""" + tel = GMOS() + hdr = fits.Header({"CCDSUM": "2 2"}) + # GMOS out_size 3200, binn 2 -> 1600 + assert tel.get_out_size(hdr) == 1600 + + +def test_get_time_get_number(): + """get_time and get_number from DATE-OBS and TIME-OBS.""" + tel = GMOS() + hdr = fits.Header({ + "DATE-OBS": "2024-06-18", + "TIME-OBS": "12:00:00", + "NCCDS": "1", + "CCDSUM": "2 2", + }) + t = tel.get_time(hdr) + assert t > 0 and np.isfinite(t) + n = tel.get_number(hdr) + assert isinstance(n, (int, np.integer)) + + +def test_get_instrument_name(): + """get_instrument_name returns lowercase name.""" + tel = GMOS() + hdr = fits.Header() + assert tel.get_instrument_name(hdr) == "gmos" + + +def test_get_catalog(): + """GMOS get_catalog returns SkyMapper for dec < -30, PS1 otherwise.""" + tel = GMOS() + hdr_n = fits.Header({"RA": 180.0, "DEC": 0.0}) + hdr_s = fits.Header({"RA": 180.0, "DEC": -35.0}) + assert tel.get_catalog(hdr_n) == "PS1" + assert tel.get_catalog(hdr_s) == "SkyMapper" + + +def test_format_datasec(): + """format_datasec converts section string with binning.""" + tel = Instrument() + out = tel.format_datasec("[1055:3024,217:3911]", binning=2) + assert "[527:1512,108:1955]" == out or out == "[527:1512,109:1956]" + + +def test_raw_format(): + """Base raw_format and GMOS raw_format.""" + base = Instrument() + assert base.raw_format(True) == "sci_img_*.fits" + assert base.raw_format(False) == "sci_img*[!proc].fits" + gmos = GMOS() + assert gmos.raw_format("dragons") == "*.fits" + assert gmos.raw_format("other") == "*.fits.bz2" + + +def test_get_stk_name_get_sci_name_get_bkg_name(): + """Naming helpers return paths under red_path.""" + tel = GMOS() + hdr = fits.Header({ + "OBJECT": "Target", + "FILTER2": "r", + "DATE-OBS": "2024-06-18", + "TIME-OBS": "12:00:00", + "NCCDS": "1", + "CCDSUM": "2 2", + }) + red_path = "/data/red" + stk = tel.get_stk_name(hdr, red_path) + assert "Target" in stk and "r" in stk and "stk.fits" in stk and red_path in stk + sci = tel.get_sci_name(hdr, red_path) + assert "Target" in sci and ".fits" in sci and red_path in sci + bkg = tel.get_bkg_name(hdr, red_path) + assert "_bkg.fits" in bkg and red_path in bkg + + +def test_get_mbias_name_get_mdark_name_get_mflat_name_get_msky_name(): + """Calibration filename helpers.""" + tel = GMOS() + paths = {"cal": "/data/red/cals"} + assert "mbias_1_22.fits" in tel.get_mbias_name(paths, "1", "22") + assert "mdark_1_22.fits" in tel.get_mdark_name(paths, "1", "22") + assert "mflat_r_1_22.fits" in tel.get_mflat_name(paths, "r", "1", "22") + assert "msky_r_1_22.fits" in tel.get_msky_name(paths, "r", "1", "22") + + +def test_base_instrument_get_ampl_get_binning_missing_keyword(): + """Base Instrument get_ampl/get_binning return default when keyword missing.""" + base = Instrument() + hdr = fits.Header() + assert base.get_ampl(hdr) == "2" + assert base.get_binning(hdr) == "CCDSUM" + + +def test_sanitize_calibration_header(): + """_sanitize_calibration_header removes WCS/coord keywords so saved cals don't trigger InvalidTransformError.""" + h = fits.Header() + h["CTYPE1"] = "RA---TAN" + h["CTYPE2"] = "DEC--TAN" + h["CRVAL1"] = 0.0 + h["CRVAL2"] = -100.0 + h["RA"] = "00:00:00" + h["DEC"] = "-100:00:00" + h["PV1_1"] = 1.0 + h["EXPTIME"] = 60.0 + h["VER"] = "1.0" + _sanitize_calibration_header(h) + assert "CTYPE1" not in h + assert "CRVAL2" not in h + assert "RA" not in h + assert "DEC" not in h + assert "PV1_1" not in h + assert h["EXPTIME"] == 60.0 + assert h["VER"] == "1.0" + + +def test_read_calibration_ccd(tmp_path): + """_read_calibration_ccd loads calibration FITS without parsing WCS (avoids ill-conditioned header errors).""" + path = tmp_path / "cal.fits" + hdu = fits.PrimaryHDU(np.zeros((10, 10), dtype=np.float32)) + hdu.header["CRVAL2"] = -100.0 # invalid; would raise if WCS were parsed + hdu.header["CTYPE1"] = "RA---TAN" + hdu.writeto(path, overwrite=True) + ccd = _read_calibration_ccd(str(path), u.electron, hdu_index=0) + assert ccd.wcs is None + assert ccd.unit == u.electron + assert ccd.data.shape == (10, 10) + + +def test_sky_subtraction_units(): + """Scaled sky (normalized * med * electron) has same unit as science so subtract is valid.""" + # Normalized master sky (dimensionless) * (med * u.electron) -> electron; then science - sky is valid + sky = CCDData(np.ones((5, 5)), unit=u.dimensionless_unscaled) + frame = CCDData(np.ones((5, 5)) * 100.0, unit=u.electron) + med = 50.0 + science_unit = frame.unit if frame.unit is not None else u.electron + frame_sky = sky.multiply(med * science_unit, propagate_uncertainties=True, handle_meta="first_found") + result = frame.subtract(frame_sky, propagate_uncertainties=True, handle_meta="first_found") + assert result.unit == u.electron + np.testing.assert_allclose(result.data, 50.0) + + +def test_get_staticmask_filename(tmp_path): + """get_staticmask_filename returns [path] when mask exists, [None] otherwise.""" + tel = GMOS() + hdr = fits.Header({"CCDSUM": "2 2"}) + # Path is paths['code']/../data/staticmasks/{instname}.{binn}.staticmask.fits.fz + code_dir = tmp_path / "code" + code_dir.mkdir() + data_dir = tmp_path / "data" / "staticmasks" + data_dir.mkdir(parents=True) + paths = {"code": str(code_dir)} + # When file does not exist + out = tel.get_staticmask_filename(hdr, paths) + assert out == [None] + # When file exists + mask_file = data_dir / "gmos.22.staticmask.fits.fz" + mask_file.touch() + out = tel.get_staticmask_filename(hdr, paths) + assert out[0] is not None + assert "gmos.22.staticmask.fits.fz" in out[0] + + +def test_load_bias_load_dark_load_flat_load_sky(tmp_path): + """load_bias, load_dark, load_flat, load_sky load from temp FITS via _read_calibration_ccd.""" + tel = GMOS() + cal_dir = tmp_path / "cals" + cal_dir.mkdir() + paths = {"cal": str(cal_dir)} + + # Bias: mbias_1_22.fits + bias_path = cal_dir / "mbias_1_22.fits" + fits.PrimaryHDU(np.zeros((10, 10), dtype=np.float32)).writeto(bias_path, overwrite=True) + mbias = tel.load_bias(paths, "1", "22") + assert mbias.unit == u.electron + assert mbias.data.shape == (10, 10) + + # Dark: mdark_1_22.fits + dark_path = cal_dir / "mdark_1_22.fits" + fits.PrimaryHDU(np.zeros((10, 10), dtype=np.float32)).writeto(dark_path, overwrite=True) + mdark = tel.load_dark(paths, "1", "22") + assert mdark.unit == u.electron + + # Flat: mflat_r_1_22.fits + flat_path = cal_dir / "mflat_r_1_22.fits" + fits.PrimaryHDU(np.ones((10, 10), dtype=np.float32)).writeto(flat_path, overwrite=True) + mflat = tel.load_flat(paths, "r", "1", "22") + assert mflat.unit == u.dimensionless_unscaled + + # Sky: msky_r_1_22.fits + sky_path = cal_dir / "msky_r_1_22.fits" + fits.PrimaryHDU(np.ones((10, 10), dtype=np.float32)).writeto(sky_path, overwrite=True) + msky = tel.load_sky(paths, "r", "1", "22") + assert msky.unit == u.dimensionless_unscaled + + +def test_load_bias_raises_when_missing(tmp_path): + """load_bias raises when no bias file exists at expected path.""" + tel = GMOS() + paths = {"cal": str(tmp_path)} + with pytest.raises(Exception, match="Could not find bias"): + tel.load_bias(paths, "1", "22") + + +def test_expand_mask(): + """expand_mask combines NaN, inf, zero pixels with optional input mask.""" + tel = Instrument() + data = np.ones((8, 8), dtype=float) + data[0, 0] = np.nan + data[1, 1] = np.inf + data[2, 2] = 0.0 + ccd = CCDData(data, unit=u.electron) + out = tel.expand_mask(ccd, input_mask=None) + assert out[0, 0] + assert out[1, 1] + assert out[2, 2] + assert not out[3, 3] + # With input mask + extra = np.zeros((8, 8), dtype=bool) + extra[4, 4] = True + out2 = tel.expand_mask(ccd, input_mask=extra) + assert out2[4, 4] + + +def test_base_get_catalog(): + """Base Instrument get_catalog returns catalog_zp.""" + base = Instrument() + hdr = fits.Header() + assert base.get_catalog(hdr) == "PS1" + + +def test_format_datasec_binning_one(): + """format_datasec with binning=1 preserves integer bounds.""" + base = Instrument() + out = base.format_datasec("[100:200,50:150]", binning=1) + assert "[100:200,50:150]" == out diff --git a/tests/test_options.py b/tests/test_options.py new file mode 100644 index 00000000..540dd311 --- /dev/null +++ b/tests/test_options.py @@ -0,0 +1,45 @@ +"""Unit tests for potpyri.utils.options.""" +import os + +import pytest + +from potpyri.utils import options +from potpyri.instruments import instrument_getter + + +def test_init_options(): + """init_options returns parser with instrument choices.""" + parser = options.init_options() + assert parser is not None + for action in parser._get_positional_actions(): + if getattr(action, "dest", None) == "instrument": + assert hasattr(action, "choices") + assert "GMOS" in action.choices + assert "LRIS" in action.choices + return + pytest.fail("instrument positional not found") + + +def test_add_paths_raises_when_data_path_missing(): + """add_paths raises when data_path does not exist.""" + tel = instrument_getter("GMOS") + with pytest.raises(Exception, match="does not exist"): + options.add_paths("/nonexistent/path/12345", "files.txt", tel) + + +def test_add_paths_creates_dirs(tmp_path): + """add_paths creates raw, bad, red, log, cal, work and returns paths dict.""" + tel = instrument_getter("GMOS") + paths = options.add_paths(str(tmp_path), "files.txt", tel) + assert "data" in paths + assert paths["data"] == os.path.abspath(tmp_path) + assert "raw" in paths + assert "red" in paths + assert "log" in paths + assert "cal" in paths + assert "work" in paths + assert "filelist" in paths + assert "code" in paths + assert os.path.exists(paths["raw"]) + assert os.path.exists(paths["red"]) + assert os.path.exists(paths["log"]) diff --git a/tests/test_utilities.py b/tests/test_utilities.py new file mode 100644 index 00000000..d26cf19b --- /dev/null +++ b/tests/test_utilities.py @@ -0,0 +1,99 @@ +"""Unit tests for potpyri.utils.utilities (find_catalog, is_number, parse_coord).""" +import pytest + +from potpyri.utils import utilities + + +def test_find_catalog_ps1(): + """find_catalog returns PS1 columns for g,r,i,z,y.""" + cat, cid, ra, dec, mag, err = utilities.find_catalog("PS1", "r", 180.0, 0.0) + assert cat == "PS1" + assert cid == "II/349" + assert ra == "RAJ2000" + assert dec == "DEJ2000" + assert mag == "rmag" + assert err == "e_rmag" + + +def test_find_catalog_ps1_unsupported_filter(): + """find_catalog returns None columns for unsupported filter.""" + cat, cid, ra, dec, mag, err = utilities.find_catalog("PS1", "U", 180.0, 0.0) + assert cid is None + assert mag is None + + +def test_find_catalog_2mass(): + """find_catalog returns 2MASS columns for J,H,K.""" + cat, cid, ra, dec, mag, err = utilities.find_catalog("2MASS", "J", 0.0, 0.0) + assert cid == "II/246" + assert mag == "Jmag" + assert err == "e_Jmag" + + +def test_find_catalog_2mass_k_band(): + """find_catalog returns 2MASS K-band (Kmag, e_Kmag) for K, Ks, and Kspec.""" + for filt in ("K", "Ks", "Kspec"): + cat, cid, ra, dec, mag, err = utilities.find_catalog("2MASS", filt, 0.0, 0.0) + assert cid == "II/246", f"2MASS catalog ID for filter {filt!r}" + assert mag == "Kmag", f"2MASS K-band mag column for filter {filt!r}" + assert err == "e_Kmag", f"2MASS K-band err column for filter {filt!r}" + + +def test_find_catalog_sdss(): + """find_catalog returns SDSS columns for u,g,r,i,z.""" + cat, cid, ra, dec, mag, err = utilities.find_catalog("SDSS", "g", 0.0, 0.0) + assert cid == "V/154" + assert mag == "gmag" + + +def test_find_catalog_skymapper_u_south(): + """find_catalog uses SkyMapper for u-band when dec < 0.""" + cat, cid, ra, dec, mag, err = utilities.find_catalog("PS1", "u", 180.0, -40.0) + assert cat == "skymapper" + assert cid == "II/379/smssdr4" + + +def test_is_number(): + """is_number returns True for numeric strings and numbers.""" + assert utilities.is_number("1.5") is True + assert utilities.is_number(1.5) is True + assert utilities.is_number("0") is True + assert utilities.is_number("abc") is False + assert utilities.is_number("") is False + + +def test_parse_coord_degrees(): + """parse_coord accepts decimal degrees.""" + coord = utilities.parse_coord(180.0, -30.0) + assert coord is not None + assert abs(coord.ra.deg - 180.0) < 0.01 + assert abs(coord.dec.deg - (-30.0)) < 0.01 + + +def test_parse_coord_sexagesimal(): + """parse_coord accepts sexagesimal strings.""" + coord = utilities.parse_coord("12:00:00", "-30:00:00") + assert coord is not None + assert abs(coord.ra.deg - 180.0) < 1.0 + assert abs(coord.dec.deg - (-30.0)) < 1.0 + + +def test_parse_coord_invalid(): + """parse_coord returns None for invalid input.""" + assert utilities.parse_coord("not", "numbers") is None + + +def test_find_catalog_ukirt(): + """find_catalog returns UKIRT columns for Y,J,H,K.""" + cat, cid, ra, dec, mag, err = utilities.find_catalog("UKIRT", "K", 0.0, 0.0) + assert cid == "II/319" + assert mag == "Kmag" + assert err == "e_Kmag" + + +def test_parse_coord_value_error(): + """parse_coord returns None when SkyCoord raises ValueError.""" + # Valid format but invalid values can trigger ValueError + result = utilities.parse_coord("99:99:99", "99:99:99") + # May return None or raise; doc says returns None on parse failure + assert result is None or hasattr(result, "ra")