diff --git a/potpyri/primitives/image_procs.py b/potpyri/primitives/image_procs.py index d0b98d87..b4b1df1b 100755 --- a/potpyri/primitives/image_procs.py +++ b/potpyri/primitives/image_procs.py @@ -18,6 +18,7 @@ import acstools import astropy.units as u +from astropy.coordinates import SkyCoord from astropy.io import fits, ascii from astropy.stats import sigma_clipped_stats from astropy.wcs import WCS @@ -266,11 +267,13 @@ def align_images(reduced_files, paths, tel, binn, use_wcs=None, fieldcenter=None def image_proc(image_data, tel, paths, skip_skysub=False, fieldcenter=None, out_size=None, satellites=True, cosmic_ray=True, - skip_gaia=False, keep_all_astro=False, log=None): + skip_gaia=False, keep_all_astro=False, relative_calibration=False, + log=None): """Full image processing: align, mask, stack, and optionally detrend. Orchestrates WCS solving, alignment, satellite/cosmic-ray masking, - stacking, and optional sky subtraction. + stacking, and optional sky subtraction. Optionally calibrates frames + to each other using SExtractor catalogs and RA/Dec cross-matching before stacking. Parameters ---------- @@ -294,6 +297,10 @@ def image_proc(image_data, tel, paths, skip_skysub=False, If True, skip Gaia alignment. Default is False. keep_all_astro : bool, optional If True, keep all images regardless of astrometric dispersion. + relative_calibration : bool, optional + If True, run SExtractor on each aligned frame, cross-match sources in + RA/Dec, and use relative fluxes to set scale factors for the stack + combine. Default is False. log : ColoredLogger, optional Logger for progress. @@ -460,7 +467,12 @@ def image_proc(image_data, tel, paths, skip_skysub=False, if log: log.info('Creating median stack.') if len(aligned_data)>1: - sci_med = stack_data(aligned_data, tel, masks, errors, log=log) + relative_scales = None + if relative_calibration: + exptimes = np.array([tel.get_exptime(d.header) for d in aligned_data]) + relative_scales = compute_relative_scales(data_images, paths, exptimes, log=log) + sci_med = stack_data(aligned_data, tel, masks, errors, log=log, + relative_scales=relative_scales) sci_med = add_stack_mask(sci_med, aligned_data) if tel.detrend: @@ -567,10 +579,139 @@ def detrend_stack(stack): return(stack) -def stack_data(stacking_data, tel, masks, errors, mem_limit=8.0e9, log=None): +def compute_relative_scales(data_images, paths, exptimes, log=None, + match_radius_arcsec=2.0, min_sources=5, min_rel_err=0.01): + """Compute relative flux scale factors from SExtractor catalogs and RA/Dec cross-matching. + + Runs Source Extractor on each aligned science image, cross-matches sources + between frames in RA/Dec, and uses the inverse-variance-weighted mean of + flux ratios (using FLUXERR_AUTO) to define scale factors so that combine() + can stack on a common flux scale. For frames where matching fails (too few + matches), the scale is set to 1.0/exptime for that frame. + + Parameters + ---------- + data_images : list of str + Paths to FITS with SCI extension (e.g. *_data.fits from image_proc). + paths : dict + Paths dict; paths.get('source_extractor') used for SExtractor binary. + exptimes : array-like + Exposure time in seconds for each image (same order as data_images). + log : ColoredLogger, optional + Logger for progress. + match_radius_arcsec : float, optional + Match radius in arcsec for cross-matching sources. Default 2.0. + min_sources : int, optional + Minimum matched sources per frame to compute scale. Default 5. + min_rel_err : float, optional + Minimum relative flux error (err/flux) used in variance of ratio to + avoid zero variance. Default 0.01. + + Returns + ------- + np.ndarray or None + Relative scale factors, one per image (reference frame scale = 1.0). + Frames that fail get scale = 1.0/exptime. None if catalogs cannot be + built at all (caller should use exposure-time-only scaling). + """ + nimg = len(data_images) + if nimg < 2: + return None + exptimes = np.atleast_1d(np.asarray(exptimes, dtype=float)) + if len(exptimes) != nimg: + if log: + log.warning('Relative calibration: exptimes length does not match data_images.') + return None + + sextractor_path = paths.get('source_extractor') + def _col(cat, name): + """Get column by case-insensitive name.""" + for c in cat.colnames: + if c.upper() == name.upper(): + return c + return None + + catalogs = [] + for p in data_images: + cat = photometry.run_sextractor(p, log=log, sextractor_path=sextractor_path) + if cat is None or len(cat) < min_sources: + if log: + log.warning(f'Relative calibration: insufficient catalog for {p}, skipping.') + return None + ra_col = _col(cat, 'ALPHA_J2000') + dec_col = _col(cat, 'DELTA_J2000') + flux_col = _col(cat, 'FLUX_AUTO') + fluxerr_col = _col(cat, 'FLUXERR_AUTO') + if not (ra_col and dec_col and flux_col): + if log: + log.warning('Relative calibration: catalog missing RA/Dec or FLUX_AUTO.') + return None + catalogs.append((cat, ra_col, dec_col, flux_col, fluxerr_col)) + + # Reference = frame with most sources + ref_idx = int(np.argmax([len(c[0]) for c in catalogs])) + ref, ref_ra, ref_dec, ref_flux, ref_fluxerr = catalogs[ref_idx] + ref_coords = SkyCoord(ref[ref_ra], ref[ref_dec], unit='deg') + + match_radius = match_radius_arcsec * u.arcsec + relative_scales = np.ones(nimg, dtype=float) + # Reference frame scale = 1.0; failed frames will get 1.0/exptime + for i in range(nimg): + if i == ref_idx: + continue + cat, ra_col, dec_col, flux_col, fluxerr_col = catalogs[i] + coords = SkyCoord(cat[ra_col], cat[dec_col], unit='deg') + # For each source in this frame, find nearest in reference + idx_ref, d2d, _ = coords.match_to_catalog_sky(ref_coords) + keep = d2d < match_radius + n_match = np.sum(keep) + if n_match < min_sources: + if log: + log.warning(f'Relative calibration: frame {i} has {n_match} matches (< {min_sources}), using 1/exptime.') + relative_scales[i] = 1.0 / exptimes[i] + continue + flux_ref = np.array(ref[ref_flux][idx_ref[keep]], dtype=float) + flux_i = np.array(cat[flux_col][keep], dtype=float) + # Flux errors: use FLUXERR_AUTO if present, else fall back to unweighted + if ref_fluxerr and fluxerr_col: + err_ref = np.array(ref[ref_fluxerr][idx_ref[keep]], dtype=float) + err_i = np.array(cat[fluxerr_col][keep], dtype=float) + else: + err_ref = err_i = None + # Avoid zeros and negatives in flux + valid = (flux_i > 0) & (flux_ref > 0) + if np.sum(valid) < min_sources: + if log: + log.warning(f'Relative calibration: frame {i} has too few valid flux pairs, using 1/exptime.') + relative_scales[i] = 1.0 / exptimes[i] + continue + ratio = flux_ref[valid] / flux_i[valid] + if err_ref is not None and err_i is not None: + err_ref = err_ref[valid] + err_i = err_i[valid] + # Relative errors with floor to avoid zero variance + rel_err_ref = np.maximum(np.abs(err_ref) / np.maximum(flux_ref[valid], 1e-30), min_rel_err) + rel_err_i = np.maximum(np.abs(err_i) / np.maximum(flux_i[valid], 1e-30), min_rel_err) + # Variance of ratio R = A/B: var(R) ≈ R^2 * ( (σ_A/A)^2 + (σ_B/B)^2 ) + var_ratio = ratio ** 2 * (rel_err_ref ** 2 + rel_err_i ** 2) + var_ratio = np.maximum(var_ratio, 1e-30) + weight = 1.0 / var_ratio + relative_scales[i] = float(np.average(ratio, weights=weight)) + else: + relative_scales[i] = float(np.median(ratio)) + + if log: + log.info(f'Relative calibration: scales (ref frame {ref_idx} = 1.0): {relative_scales}') + return relative_scales + + +def stack_data(stacking_data, tel, masks, errors, mem_limit=8.0e9, log=None, + relative_scales=None): """Combine aligned CCDData list with exposure scaling (median or average). Masks are applied (masked pixels set to nan) before combination. + Optionally uses relative flux scales from source cross-matching so that + frames are combined on a common flux scale. Parameters ---------- @@ -586,6 +727,9 @@ def stack_data(stacking_data, tel, masks, errors, mem_limit=8.0e9, log=None): Memory limit in bytes for ccdproc.combine. Default is 8e9. log : ColoredLogger, optional Logger for progress. + relative_scales : np.ndarray, optional + Relative flux scale factors (one per image; reference = 1.0). If + provided, scale = relative_scales (exposure-time scaling is not applied). Returns ------- @@ -608,7 +752,9 @@ def stack_data(stacking_data, tel, masks, errors, mem_limit=8.0e9, log=None): exptimes.append(float(tel.get_exptime(stk.header))) exptimes = np.array(exptimes) - scale = 1./exptimes + scale = 1.0 / exptimes + if relative_scales is not None: + scale = np.asarray(relative_scales, dtype=float) if len(scale)==1: stack_method='average' diff --git a/potpyri/primitives/photometry.py b/potpyri/primitives/photometry.py index 688211af..cd1796c2 100755 --- a/potpyri/primitives/photometry.py +++ b/potpyri/primitives/photometry.py @@ -79,7 +79,7 @@ def create_params(outfile): f.write('FLUXERR_AUTO \n') f.write('FWHM_IMAGE \n') -def run_sextractor(img_file, log=None): +def run_sextractor(img_file, log=None, sextractor_path=None): """Run Source Extractor on SCI extension; return catalog table or None. Parameters @@ -88,6 +88,8 @@ def run_sextractor(img_file, log=None): Path to FITS with SCI extension (and SATURATE in header). log : ColoredLogger, optional Logger for progress. + sextractor_path : str, optional + Path to Source Extractor binary. If None, uses 'sex' from PATH. Returns ------- @@ -113,12 +115,13 @@ def run_sextractor(img_file, log=None): datahdu.header = hdu['SCI'].header datahdu.writeto(tmpfile, overwrite=True) + sex_cmd = (sextractor_path if sextractor_path else 'sex').strip() if log: log.info(f'Running source extractor on {img_file}') else: print(f'Running source extractor on {img_file}') - cmd = f'sex {tmpfile} -CATALOG_NAME {catfile} -CATALOG_TYPE ASCII_HEAD ' + cmd = f'{sex_cmd} {tmpfile} -CATALOG_NAME {catfile} -CATALOG_TYPE ASCII_HEAD ' cmd += f'-PARAMETERS_NAME {paramfile} -FILTER_NAME {convfile} ' cmd += f'-SATUR_LEVEL {saturate} > /dev/null 2> /dev/null' diff --git a/potpyri/scripts/main_pipeline.py b/potpyri/scripts/main_pipeline.py index 18fafacf..5583e672 100755 --- a/potpyri/scripts/main_pipeline.py +++ b/potpyri/scripts/main_pipeline.py @@ -47,6 +47,7 @@ def main_pipeline(instrument: str, skip_cr: bool = None, skip_gaia: bool = None, keep_all_astro: bool = None, + relative_calibration: bool = None, **kwargs) -> None: """Run the full reduction pipeline: sort files, calibrate, process, stack, photometry, zeropoint. @@ -89,7 +90,7 @@ def main_pipeline(instrument: str, stack = image_procs.image_proc(file_table[file_table['TargType']==tar], tel, paths, skip_skysub=skip_skysub, fieldcenter=fieldcenter, out_size=out_size, cosmic_ray=not skip_cr, skip_gaia=skip_gaia, keep_all_astro=keep_all_astro, - log=log) + relative_calibration=relative_calibration or False, log=log) # Photometry step ################# diff --git a/potpyri/utils/options.py b/potpyri/utils/options.py index abcecd09..b90a59c2 100755 --- a/potpyri/utils/options.py +++ b/potpyri/utils/options.py @@ -98,6 +98,11 @@ def init_options(): default=False, action='store_true', help='Tell the pipeline to keep all images regardless of astrometric dispersion.') + params.add_argument('--relative-calibration', + default=False, + action='store_true', + help='Before stacking, calibrate frames to each other using SExtractor catalogs ' + 'and RA/Dec cross-matching; use relative source fluxes to set combine scales.') return(params) diff --git a/tests/test_image_procs.py b/tests/test_image_procs.py index ffbe1633..c04f313a 100644 --- a/tests/test_image_procs.py +++ b/tests/test_image_procs.py @@ -1,11 +1,13 @@ -"""Tests for image_procs: remove_pv_distortion, get_fieldcenter, generate_wcs, detrend_stack, add_stack_mask.""" +"""Tests for image_procs: remove_pv_distortion, get_fieldcenter, generate_wcs, detrend_stack, add_stack_mask, compute_relative_scales.""" from potpyri.primitives import image_procs from potpyri.instruments import instrument_getter import numpy as np import os +from unittest.mock import patch from astropy.io import fits +from astropy.table import Table from astropy.wcs import WCS from ccdproc import CCDData from astropy import units as u @@ -122,3 +124,167 @@ def test_create_error(tmp_path): assert err_hdu.header.get("BUNIT") == "ELECTRONS" assert np.all(err_hdu.data > 0) assert np.all(np.isfinite(err_hdu.data)) + + +def _make_catalog(n, ra=30.0, dec=-30.0, flux=100.0, fluxerr=5.0, dra=0.001, ddec=0.001): + """Build a minimal SExtractor-like table with ALPHA_J2000, DELTA_J2000, FLUX_AUTO, FLUXERR_AUTO.""" + rng = np.random.default_rng(42) + ra_arr = ra + (rng.random(n) - 0.5) * dra + dec_arr = dec + (rng.random(n) - 0.5) * ddec + flux_arr = np.full(n, flux, dtype=float) + (rng.random(n) - 0.5) * 0.1 * flux + err_arr = np.full(n, fluxerr, dtype=float) + t = Table() + t['ALPHA_J2000'] = ra_arr + t['DELTA_J2000'] = dec_arr + t['FLUX_AUTO'] = np.maximum(flux_arr, 1.0) + t['FLUXERR_AUTO'] = err_arr + return t + + +def test_compute_relative_scales_returns_none_for_single_image(): + """compute_relative_scales returns None when fewer than 2 images.""" + paths = {} + data_images = ['/fake/path0.fits'] + exptimes = [60.0] + out = image_procs.compute_relative_scales(data_images, paths, exptimes) + assert out is None + + +def test_compute_relative_scales_returns_none_for_exptimes_length_mismatch(): + """compute_relative_scales returns None when exptimes length != number of images.""" + paths = {} + data_images = ['/fake/a.fits', '/fake/b.fits'] + exptimes = [60.0, 30.0, 90.0] + with patch.object(image_procs.photometry, 'run_sextractor') as m_sextractor: + m_sextractor.return_value = _make_catalog(10) + out = image_procs.compute_relative_scales(data_images, paths, exptimes) + assert out is None + + +def test_compute_relative_scales_returns_none_when_sextractor_fails(): + """compute_relative_scales returns None when run_sextractor returns None for one image.""" + paths = {} + data_images = ['/fake/a.fits', '/fake/b.fits'] + exptimes = [60.0, 60.0] + + def side_effect(path, log=None, sextractor_path=None): + if 'a.fits' in path: + return _make_catalog(10) + return None + + with patch.object(image_procs.photometry, 'run_sextractor', side_effect=side_effect): + out = image_procs.compute_relative_scales(data_images, paths, exptimes) + assert out is None + + +def test_compute_relative_scales_returns_none_when_catalog_too_small(): + """compute_relative_scales returns None when a catalog has fewer than min_sources.""" + paths = {} + data_images = ['/fake/a.fits', '/fake/b.fits'] + exptimes = [60.0, 60.0] + small = _make_catalog(2) + big = _make_catalog(10) + + def side_effect(path, log=None, sextractor_path=None): + return small if 'a.fits' in path else big + + with patch.object(image_procs.photometry, 'run_sextractor', side_effect=side_effect): + out = image_procs.compute_relative_scales(data_images, paths, exptimes, min_sources=5) + assert out is None + + +def test_compute_relative_scales_two_frames_flux_ratio(): + """compute_relative_scales returns [1, scale] when frame 1 has half the flux of reference.""" + paths = {} + data_images = ['/fake/ref.fits', '/fake/f1.fits'] + exptimes = [60.0, 60.0] + # Same positions so they match; ref flux 100, frame1 flux 50 -> ratio 2.0 + ref_cat = _make_catalog(10, flux=100.0, fluxerr=5.0) + f1_cat = _make_catalog(10, flux=50.0, fluxerr=2.5) + # Use identical RA/Dec so match_to_catalog_sky finds pairs (same seed and pattern) + f1_cat['ALPHA_J2000'] = ref_cat['ALPHA_J2000'].copy() + f1_cat['DELTA_J2000'] = ref_cat['DELTA_J2000'].copy() + + def side_effect(path, log=None, sextractor_path=None): + return ref_cat if 'ref.fits' in path else f1_cat + + with patch.object(image_procs.photometry, 'run_sextractor', side_effect=side_effect): + out = image_procs.compute_relative_scales(data_images, paths, exptimes, min_sources=5) + assert out is not None + assert len(out) == 2 + assert out[0] == 1.0 + np.testing.assert_allclose(out[1], 2.0, rtol=0.15) + + +def test_compute_relative_scales_frame_fails_gets_exptime_scale(): + """When one frame has too few matches, that frame gets scale = 1.0/exptime.""" + paths = {} + data_images = ['/fake/ref.fits', '/fake/f1.fits', '/fake/f2.fits'] + exptimes = [60.0, 30.0, 90.0] + ref_cat = _make_catalog(10, ra=30.0, dec=-30.0, flux=100.0) + f1_cat = _make_catalog(10, ra=30.0, dec=-30.0, flux=50.0) + f1_cat['ALPHA_J2000'] = ref_cat['ALPHA_J2000'].copy() + f1_cat['DELTA_J2000'] = ref_cat['DELTA_J2000'].copy() + # Frame 2: completely different sky position -> no matches + f2_cat = _make_catalog(10, ra=100.0, dec=50.0, flux=100.0) + + def side_effect(path, log=None, sextractor_path=None): + if 'ref.fits' in path: + return ref_cat + if 'f1.fits' in path: + return f1_cat + return f2_cat + + with patch.object(image_procs.photometry, 'run_sextractor', side_effect=side_effect): + out = image_procs.compute_relative_scales(data_images, paths, exptimes, min_sources=5, match_radius_arcsec=2.0) + assert out is not None + assert len(out) == 3 + # Ref (most sources) is frame 0; scale[0]=1 + assert out[0] == 1.0 + # Frame 1 matches ref -> computed ratio ~2.0 + np.testing.assert_allclose(out[1], 2.0, rtol=0.2) + # Frame 2 has no matches within 2" -> 1/exptime = 1/90 + np.testing.assert_allclose(out[2], 1.0 / 90.0, rtol=1e-5) + + +def test_compute_relative_scales_uses_fluxerr_weighting(): + """compute_relative_scales uses FLUXERR_AUTO when present (weighted mean).""" + paths = {} + data_images = ['/fake/ref.fits', '/fake/f1.fits'] + exptimes = [60.0, 60.0] + ref_cat = _make_catalog(10, flux=100.0, fluxerr=2.0) + f1_cat = _make_catalog(10, flux=50.0, fluxerr=1.0) + f1_cat['ALPHA_J2000'] = ref_cat['ALPHA_J2000'].copy() + f1_cat['DELTA_J2000'] = ref_cat['DELTA_J2000'].copy() + + def side_effect(path, log=None, sextractor_path=None): + return ref_cat if 'ref.fits' in path else f1_cat + + with patch.object(image_procs.photometry, 'run_sextractor', side_effect=side_effect): + out = image_procs.compute_relative_scales(data_images, paths, exptimes, min_sources=5) + assert out is not None + assert out[0] == 1.0 + # Weighted mean of ratio ~2.0 (with small errors) should be close to 2 + np.testing.assert_allclose(out[1], 2.0, rtol=0.15) + + +def test_compute_relative_scales_without_fluxerr_uses_median(): + """When FLUXERR_AUTO is missing, compute_relative_scales uses median of ratios.""" + paths = {} + data_images = ['/fake/ref.fits', '/fake/f1.fits'] + exptimes = [60.0, 60.0] + ref_cat = _make_catalog(10, flux=100.0, fluxerr=5.0) + ref_cat.remove_column('FLUXERR_AUTO') + f1_cat = _make_catalog(10, flux=50.0, fluxerr=2.5) + f1_cat.remove_column('FLUXERR_AUTO') + f1_cat['ALPHA_J2000'] = ref_cat['ALPHA_J2000'].copy() + f1_cat['DELTA_J2000'] = ref_cat['DELTA_J2000'].copy() + + def side_effect(path, log=None, sextractor_path=None): + return ref_cat if 'ref.fits' in path else f1_cat + + with patch.object(image_procs.photometry, 'run_sextractor', side_effect=side_effect): + out = image_procs.compute_relative_scales(data_images, paths, exptimes, min_sources=5) + assert out is not None + assert out[0] == 1.0 + np.testing.assert_allclose(out[1], 2.0, rtol=0.2) diff --git a/tests/test_instruments.py b/tests/test_instruments.py index 9746dac6..b67c1607 100644 --- a/tests/test_instruments.py +++ b/tests/test_instruments.py @@ -334,8 +334,6 @@ def test_format_datasec_binning_one(): def test_create_sky_masks_bright_sources(tmp_path): """create_sky uses iterative sigma clipping and masks bright sources so sky estimate is robust.""" - from astropy.stats import sigma_clipped_stats - tel = GMOS() cal_dir = tmp_path / "cal" cal_dir.mkdir()