From eca47a2904a43fbcf6f54a3d5f7f869afe120dfd Mon Sep 17 00:00:00 2001 From: charliekilpatrick Date: Thu, 5 Mar 2026 18:45:31 -0600 Subject: [PATCH] Enhance create_sky method with iterative sigma clipping and robust masking - Updated the create_sky method to include configurable parameters for iterative sigma clipping, allowing for better estimation of sky background and masking of bright sources. - Improved handling of NaN values during normalization and masking processes. - Added unit tests to validate the new functionality, ensuring robust performance with both default and custom parameters. --- potpyri/instruments/instrument.py | 123 +++++++++++++++++++++--------- tests/test_instruments.py | 68 +++++++++++++++++ 2 files changed, 157 insertions(+), 34 deletions(-) diff --git a/potpyri/instruments/instrument.py b/potpyri/instruments/instrument.py index 3d7c0e9..5a644ba 100755 --- a/potpyri/instruments/instrument.py +++ b/potpyri/instruments/instrument.py @@ -845,9 +845,18 @@ def create_flat(self, flat_list, fil, amp, binn, paths, mbias=None, return - def create_sky(self, sky_list, fil, amp, binn, paths, log=None, **kwargs): + def create_sky(self, sky_list, fil, amp, binn, paths, log=None, + sky_sigma_upper=3.0, sky_sigma_lower=4.0, sky_maxiters=5, + sky_n_sigma_high=3.0, sky_n_sigma_low=4.0, + msky_n_sigma_high=5.0, msky_n_sigma_low=5.0, msky_maxiters=5, + **kwargs): """Combine sky frames (normalize by median) and write master sky. + Uses iterative sigma clipping to estimate the sky background and mask + bright sources (and optionally low outliers) so the median is not + biased. Per-frame and combined master sky masking use configurable + sigma and iteration limits. + Parameters ---------- sky_list : list of str @@ -862,6 +871,24 @@ def create_sky(self, sky_list, fil, amp, binn, paths, log=None, **kwargs): Paths dict (cal key). log : ColoredLogger, optional Logger for progress. + sky_sigma_upper : float, optional + Upper sigma for iterative clip when estimating per-frame sky. + Default 3.0. + sky_sigma_lower : float, optional + Lower sigma for iterative clip when estimating per-frame sky. + Default 4.0. + sky_maxiters : int, optional + Max iterations for sigma clipping per sky frame. Default 5. + sky_n_sigma_high : float, optional + Mask pixels above med + this many * std (bright sources). Default 3.0. + sky_n_sigma_low : float, optional + Mask pixels below med - this many * std (defects). Default 4.0. + msky_n_sigma_high : float, optional + Mask master sky pixels above median + this * std. Default 5.0. + msky_n_sigma_low : float, optional + Mask master sky pixels below median - this * std. Default 5.0. + msky_maxiters : int, optional + Max iterations for sigma clipping on combined master sky. Default 5. **kwargs Unused; for subclass compatibility. @@ -874,59 +901,82 @@ def create_sky(self, sky_list, fil, amp, binn, paths, log=None, **kwargs): log.info(f'Processing files for filter: {fil}') log.info(f'{len(sky_list)} files found.') - scale = [] skys = [] for i, sky in enumerate(sky_list): if log: log.info(f'Importing {sky}') sky_full = self.import_sci_image(sky, log=log) - mean, med, stddev = sigma_clipped_stats(sky_full.data, - sigma_upper=2.5, sigma_lower=3.5) - - # Mask outliers - mask = sky_full.data > med + 2.5 * stddev - sky_full.data[mask]=np.nan + # Iterative sigma-clipped stats for robust sky level (excludes + # bright sources and defects so median is unbiased) + mean, med, stddev = sigma_clipped_stats( + sky_full.data, + sigma_upper=sky_sigma_upper, + sigma_lower=sky_sigma_lower, + maxiters=sky_maxiters, + ) - # Normalize by median sky background - mean, med, stddev = sigma_clipped_stats(sky_full.data, - sigma_upper=2.5, sigma_lower=3.5) - norm = 1./med + # Mask bright sources and low outliers so they don't affect norm + mask_high = sky_full.data > med + sky_n_sigma_high * stddev + mask_low = sky_full.data < med - sky_n_sigma_low * stddev + sky_full.data[mask_high] = np.nan + sky_full.data[mask_low] = np.nan + + # Recompute stats on masked data (NaNs excluded) for normalization + mean, med, stddev = sigma_clipped_stats( + sky_full.data, + sigma_upper=sky_sigma_upper, + sigma_lower=sky_sigma_lower, + maxiters=sky_maxiters, + ) + norm = 1.0 / med sky_full = sky_full.multiply(norm) - + # Vet the sky normalization - it should not be negative if norm > 0.: - log.info(f'Sky normalization: {norm}') + if log: log.info(f'Sky normalization: {norm}') else: - # Skip this file - log.error(f'Sky normalization: {norm}') + if log: log.error(f'Sky normalization: {norm}') continue - sky_full.mask[np.isnan(sky_full.data)]=True - sky_full.data[np.isnan(sky_full.data)]=1.0 - + nan_mask = np.isnan(sky_full.data) + if sky_full.mask is None: + sky_full.mask = np.zeros(sky_full.data.shape, dtype=bool) + sky_full.mask[nan_mask] = True + sky_full.data[nan_mask] = 1.0 + skys.append(sky_full) - msky = ccdproc.combine(skys, method='median', sigma_clip=True, + msky = ccdproc.combine(skys, method='median', sigma_clip=True, clip_extrema=True) - # Mask sky image - msky.data[np.isinf(msky.data)]=1.0 - msky.data[msky.data==0.0]=1.0 - mean, median, stddev = sigma_clipped_stats(msky.data) - mask = msky.data > median + 10 * stddev - msky.data[mask]=1.0 - - if log: + # Robust masking of combined master sky: iterative sigma clipping + # for accurate median/std, then mask excess flux and defects + msky.data[np.isinf(msky.data)] = np.nan + msky.data[msky.data == 0.0] = np.nan + mean, median, stddev = sigma_clipped_stats( + msky.data, + sigma_upper=3.0, + sigma_lower=3.0, + maxiters=msky_maxiters, + ) + mask_high = msky.data > median + msky_n_sigma_high * stddev + mask_low = msky.data < median - msky_n_sigma_low * stddev + msky.data[mask_high] = 1.0 + msky.data[mask_low] = 1.0 + # Restore valid fill for normalized sky (1.0 = no correction) + msky.data[np.isnan(msky.data)] = 1.0 + + if log: log.info(f'Made sky for filter: {fil}, amp: {amp}, bin: {binn}.') - - msky.header['VER'] = (__version__, + + msky.header['VER'] = (__version__, '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') if log: log.info(f'Master sky written to {sky_filename}') - + return def load_staticmask(self, hdr, paths): @@ -1161,9 +1211,14 @@ def process_science(self, sci_list, fil, amp, binn, paths, mbias=None, sky_frame = self.load_sky(paths, fil, amp, binn) - for i,frame in enumerate(processed): - - mean, med, stddev = sigma_clipped_stats(frame.data) + for i, frame in enumerate(processed): + # Robust scale: sigma-clipped median so bright sources don't bias sky level + mean, med, stddev = sigma_clipped_stats( + frame.data, + sigma_upper=3.0, + sigma_lower=3.0, + maxiters=5, + ) # 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, diff --git a/tests/test_instruments.py b/tests/test_instruments.py index c2f7b7c..9746dac 100644 --- a/tests/test_instruments.py +++ b/tests/test_instruments.py @@ -330,3 +330,71 @@ def test_format_datasec_binning_one(): base = Instrument() out = base.format_datasec("[100:200,50:150]", binning=1) assert "[100:200,50:150]" == out + + +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() + sky_dir = tmp_path / "sky" + sky_dir.mkdir() + paths = {"cal": str(cal_dir)} + + # Two small sky frames: constant background + one frame with a bright source + sky_val = 1000.0 + shape = (12, 12) + for i in range(2): + data = np.full(shape, sky_val, dtype=np.float32) + if i == 1: + data[5, 5] = 8000.0 # bright source + data[6, 5] = 6000.0 + ccd = CCDData(data, unit=u.electron) + path = sky_dir / f"sky{i}.fits" + ccd.write(str(path), overwrite=True) + + sky_list = [str(sky_dir / "sky0.fits"), str(sky_dir / "sky1.fits")] + tel.create_sky(sky_list, "r", "1", "22", paths, log=None, + sky_sigma_upper=3.0, sky_sigma_lower=4.0, sky_maxiters=5, + sky_n_sigma_high=3.0, sky_n_sigma_low=4.0, + msky_n_sigma_high=5.0, msky_n_sigma_low=5.0, msky_maxiters=5) + + msky_path = cal_dir / "msky_r_1_22.fits" + assert msky_path.exists() + with fits.open(msky_path) as hdu: + msky_data = hdu[0].data + # Normalized sky: masked (bright) pixels set to 1.0, rest ~1.0 + assert msky_data.shape == shape + # At (5,5) and (6,5) we had bright flux; after masking they should be 1.0 + assert msky_data[5, 5] == 1.0 + assert msky_data[6, 5] == 1.0 + # Median of combined sky should be close to 1.0 (normalized) + med = np.nanmedian(msky_data) + assert 0.95 < med < 1.05 + + +def test_create_sky_accepts_default_and_custom_sigma_params(tmp_path): + """create_sky runs with default params and with custom sigma/maxiters.""" + tel = GMOS() + cal_dir = tmp_path / "cal" + cal_dir.mkdir() + sky_dir = tmp_path / "sky" + sky_dir.mkdir() + paths = {"cal": str(cal_dir)} + data = np.full((8, 8), 500.0, dtype=np.float32) + ccd = CCDData(data, unit=u.electron) + sky_path = sky_dir / "sky.fits" + ccd.write(str(sky_path), overwrite=True) + sky_list = [str(sky_path)] + + # Default params + tel.create_sky(sky_list, "r", "1", "22", paths, log=None) + assert (cal_dir / "msky_r_1_22.fits").exists() + + # Custom params (tight sigma) + (cal_dir / "msky_r_1_22.fits").unlink() + tel.create_sky(sky_list, "r", "1", "22", paths, log=None, + sky_n_sigma_high=2.0, msky_n_sigma_high=4.0) + assert (cal_dir / "msky_r_1_22.fits").exists()