From a0956e216baefa4bae30794c0cdd0ae1a86ef061 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Tue, 27 Aug 2024 20:26:10 -0400 Subject: [PATCH 01/20] updating ReductionStage bitmask & propagating through science reduction steps (fixes #70) --- python/lvmdrp/functions/imageMethod.py | 73 ++++++++++++++++++-------- python/lvmdrp/functions/rssMethod.py | 7 +++ python/lvmdrp/functions/skyMethod.py | 4 ++ python/lvmdrp/utils/bitmask.py | 60 ++++++++++++++++----- 4 files changed, 107 insertions(+), 37 deletions(-) diff --git a/python/lvmdrp/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index d3482847..8acbfbae 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -28,7 +28,7 @@ from lvmdrp import log, __version__ as DRPVER from lvmdrp.core.constants import CONFIG_PATH, SPEC_CHANNELS, ARC_LAMPS, LVM_REFERENCE_COLUMN, LVM_NBLOCKS, FIDUCIAL_PLATESCALE from lvmdrp.utils.decorators import skip_on_missing_input_path, drop_missing_input_paths -from lvmdrp.utils.bitmask import QualityFlag +from lvmdrp.utils.bitmask import QualityFlag, ReductionStage from lvmdrp.core.fiberrows import FiberRows, _read_fiber_ypix from lvmdrp.core.image import ( Image, @@ -98,7 +98,7 @@ ] -def _nonlinearity_correction(ptc_params: None | numpy.ndarray, nominal_gain: float, quadrant: Image, iquad: int) -> Image: +def _nonlinearity_correction(ptc_params: None | numpy.ndarray, nominal_gain: float, quadrant: Image, iquad: int, drpstage: ReductionStage) -> Tuple[Image,ReductionStage]: """calculates non-linearity correction for input quadrant Parameters @@ -134,11 +134,12 @@ def _nonlinearity_correction(ptc_params: None | numpy.ndarray, nominal_gain: flo gain_med = numpy.nanmedian(gain_map._data) gain_min, gain_max = numpy.nanmin(gain_map._data), numpy.nanmax(gain_map._data) log.info(f"gain map stats: {gain_med = :.2f} [{gain_min = :.2f}, {gain_max = :.2f}] ({nominal_gain = :.2f} e-/ADU)") + drpstage += "LINEARITY_CORRECTED" else: log.warning("cannot apply non-linearity correction") log.info(f"using {nominal_gain = } (e-/ADU)") gain_map = Image(data=numpy.ones(quadrant._data.shape) * nominal_gain) - return gain_map + return gain_map, drpstage def _create_peaks_regions(fibermap: Table, column: int = 2000) -> None: @@ -390,6 +391,9 @@ def _fix_fiber_thermal_shifts(image, trace_cent, trace_width=None, trace_amp=Non axs_cc, axs_fb = axs # calculate thermal shifts column_shifts = image.measure_fiber_shifts(fiber_model, columns=columns, column_width=column_width, shift_range=shift_range, axs=axs_cc) + if (column_shifts!=0).any(): + image._header["DRPSTAGE"] = (ReductionStage(image._header["DRPSTAGE"]) + "FIBERS_SHIFTED").value + # shifts stats median_shift = numpy.nanmedian(column_shifts, axis=0) std_shift = numpy.nanstd(column_shifts, axis=0) @@ -441,7 +445,7 @@ def _fix_fiber_thermal_shifts(image, trace_cent, trace_width=None, trace_amp=Non # trace_cent_fixed._coeffs[ifiber] = (poly_trace + poly_deltas).coef # trace_cent_fixed.eval_coeffs() - return trace_cent_fixed, column_shifts, fiber_model + return image, trace_cent_fixed, column_shifts, fiber_model def _apply_electronic_shifts(images, out_images, drp_shifts=None, qc_shifts=None, custom_shifts=None, raw_shifts=None, @@ -2558,6 +2562,7 @@ def subtract_straylight( # load image data log.info(f"using image {os.path.basename(in_image)} for stray light subtraction") img = loadImage(in_image) + drpstage = ReductionStage(img._header["DRPSTAGE"]) unit = img._header["BUNIT"] # smooth image along dispersion axis with a median filter excluded NaN values @@ -2618,10 +2623,12 @@ def subtract_straylight( log.info("subtracting the smoothed background signal from the original image") img_out = loadImage(in_image) img_out._data = img_out._data - img_stray._data + drpstage += "STRAYLIGHT_SUBTRACTED" # include header and write out file log.info(f"writing stray light subtracted image to {os.path.basename(out_image)}") img_out.setHeader(header=img.getHeader()) + img_out._header["DRPSTAGE"] = drpstage.value img_out.writeFitsData(out_image) # plot results: polyomial fitting & smoothing, both with masked regions on @@ -3235,6 +3242,7 @@ def extract_spectra( log.info(f"extraction using aperture of {aperture} pixels") img = loadImage(in_image) + drpstage = ReductionStage(img._header["DRPSTAGE"]) mjd = img._header["SMJD"] camera = img._header["CCD"] expnum = img._header["EXPOSURE"] @@ -3274,12 +3282,12 @@ def extract_spectra( # fix centroids for thermal shifts log.info(f"measuring fiber thermal shifts @ columns: {','.join(map(str, columns))}") - trace_mask, shifts, _ = _fix_fiber_thermal_shifts(img, trace_mask, 2.5, - fiber_model=fiber_model, - trace_amp=10000, - columns=columns, - column_width=column_width, - shift_range=shift_range, axs=[axs_cc, axs_fb]) + img, trace_mask, shifts, _ = _fix_fiber_thermal_shifts(img, trace_mask, 2.5, + fiber_model=fiber_model, + trace_amp=10000, + columns=columns, + column_width=column_width, + shift_range=shift_range, axs=[axs_cc, axs_fb]) # save columns measured for thermal shifts plot_fiber_thermal_shift(columns, shifts, ax=ax_shift) save_fig(fig, product_path=out_rss, to_display=display_plots, figure_path="qa", label="fiber_thermal_shifts") @@ -3360,6 +3368,8 @@ def extract_spectra( mask |= (~(numpy.isin(slitmap_spec["orig_ifulabel"], exposed_std))&((slitmap_spec["telescope"] == "Spec")))[:, None] mask |= (slitmap_spec["fibstatus"] == 1)[:, None] + drpstage += "SPECTRA_EXTRACTED" + # propagate thermal shift to slitmap channel = img._header['CCD'][0] slitmap[f"ypix_{channel}"] = slitmap[f"ypix_{channel}"].astype("float64") @@ -3377,6 +3387,7 @@ def extract_spectra( header=img.getHeader(), slitmap=slitmap ) + rss._header["DRPSTAGE"] = drpstage.value rss.setHdrValue("NAXIS2", data.shape[0]) rss.setHdrValue("NAXIS1", data.shape[1]) rss.setHdrValue("DISPAXIS", 1) @@ -3775,6 +3786,9 @@ def preproc_raw_frame( display_plots : bool, optional whether to show plots on display or not, by default False """ + # initialize reduction status + drpstage = ReductionStage(1) + # load image log.info(f"starting preprocessing of raw image '{os.path.basename(in_image)}'") org_img = loadImage(in_image) @@ -3788,6 +3802,7 @@ def preproc_raw_frame( sjd = int(dateobs_to_sjd(org_header.get("OBSTIME"))) sjd = correct_sjd(in_image, sjd) org_header = apply_hdrfix(sjd, hdr=org_header) or org_header + drpstage += "HDRFIX_APPLIED" except ValueError as e: log.error(f"cannot apply header fix: {e}") @@ -3869,6 +3884,7 @@ def preproc_raw_frame( gain[2] *= 1.089 if camera == "z3": gain[3] /= 1.056 + gain = numpy.round(gain, 3) log.info(f"using header GAIN = {gain.tolist()} (e-/ADU)") @@ -3906,6 +3922,9 @@ def preproc_raw_frame( os_profiles.append(os_profile) os_models.append(os_model) + if subtract_overscan: + drpstage += "OVERSCAN_SUBTRACTED" + # compute overscan stats os_bias_med[i] = numpy.nanmedian(os_quad._data, axis=None) os_bias_std[i] = numpy.nanmedian(numpy.nanstd(os_quad._data, axis=1), axis=None) @@ -3995,6 +4014,7 @@ def preproc_raw_frame( if in_mask and proc_img._header["IMAGETYP"] not in {"bias", "dark", "pixflat"}: log.info(f"loading master pixel mask from {os.path.basename(in_mask)}") master_mask = loadImage(in_mask)._mask.astype(bool) + drpstage += "PIXELMASK_ADDED" else: master_mask = numpy.zeros_like(proc_img._data, dtype=bool) @@ -4018,13 +4038,12 @@ def preproc_raw_frame( drpqual = QualityFlag(0) if saturated_mask.sum() / proc_img._mask.size > 0.01: drpqual += "SATURATED" - proc_img.setHdrValue("DRPQUAL", value=drpqual.value, comment="data reduction quality flag") - - # set drp tag version - proc_img.setHdrValue("DRPVER", DRPVER, comment='data reduction pipeline software tag') # write out FITS file log.info(f"writing preprocessed image to {os.path.basename(out_image)}") + proc_img.setHdrValue("DRPVER", DRPVER, comment='data reduction pipeline version') + proc_img.setHdrValue("DRPSTAGE", drpstage.value, comment="data reduction stage") + proc_img.setHdrValue("DRPQUAL", value=drpqual.value, comment="data reduction quality flag") proc_img.writeFitsData(out_image) # plot overscan strips along X and Y axes @@ -4175,6 +4194,7 @@ def add_astrometry( # reading slitmap org_img = loadImage(in_image) + drpstage = ReductionStage(org_img._header["DRPSTAGE"]) slitmap = org_img.getSlitmap() telescope=numpy.array(slitmap['telescope'].data) x=numpy.array(slitmap['xpmm'].data) @@ -4195,7 +4215,7 @@ def copy_guider_keyword(gdrhdr, keyword, img): comment = gdrhdr.comments[keyword] if inhdr else '' img.setHdrValue(f'HIERARCH GDRCOADD {keyword}', gdrhdr.get(keyword), comment) - def getobsparam(tel): + def getobsparam(tel, drpstage): if tel!='spec': if os.path.isfile(agcfiledir[tel]): mfagc=fits.open(agcfiledir[tel]) @@ -4243,16 +4263,17 @@ def getobsparam(tel): log.warning(f"some astrometry keywords for telescope '{tel}' are missing: {RAobs = }, {DECobs = }, {PAobs = }") org_img.add_header_comment(f"no astromentry keywords '{tel}': {RAobs = }, {DECobs = }, {PAobs = }, using commanded") org_img.setHdrValue('ASTRMSRC', 'CMD position', comment='Source of astrometric solution: commanded position') + drpstage += f"{tel.upper()}_ASTROMETRY_ADDED" else: RAobs=0 DECobs=0 PAobs=0 - return RAobs, DECobs, PAobs + return RAobs, DECobs, PAobs, drpstage - RAobs_sci, DECobs_sci, PAobs_sci = getobsparam('sci') - RAobs_skye, DECobs_skye, PAobs_skye = getobsparam('skye') - RAobs_skyw, DECobs_skyw, PAobs_skyw = getobsparam('skyw') - RAobs_spec, DECobs_spec, PAobs_spec = getobsparam('spec') + RAobs_sci, DECobs_sci, PAobs_sci, drpstage = getobsparam('sci', drpstage) + RAobs_skye, DECobs_skye, PAobs_skye, drpstage = getobsparam('skye', drpstage) + RAobs_skyw, DECobs_skyw, PAobs_skyw, drpstage = getobsparam('skyw', drpstage) + RAobs_spec, DECobs_spec, PAobs_spec, drpstage = getobsparam('spec', drpstage) # Create fake IFU image WCS object for each telescope focal plane and use it to calculate RA,DEC of each fiber telcoordsdir={'sci':(RAobs_sci, DECobs_sci, PAobs_sci), 'skye':(RAobs_skye, DECobs_skye, PAobs_skye), 'skyw':(RAobs_skyw, DECobs_skyw, PAobs_skyw), 'spec':(RAobs_spec, DECobs_spec, PAobs_spec)} @@ -4290,13 +4311,13 @@ def getfibradec(tel, platescale): slitmap['ra']=RAfib slitmap['dec']=DECfib org_img._slitmap=slitmap + drpstage += "FIBERS_ASTROMETRY_ADDED" log.info(f"writing RA,DEC to slitmap in image '{os.path.basename(out_image)}'") + org_img._header["DRPSTAGE"] = drpstage.value org_img.writeFitsData(out_image) - - @skip_on_missing_input_path(["in_image"]) # @skip_if_drpqual_flags(["SATURATED"], "in_image") def detrend_frame( @@ -4351,6 +4372,7 @@ def detrend_frame( org_img = loadImage(in_image) exptime = org_img._header["EXPTIME"] img_type = org_img._header["IMAGETYP"].lower() + drpstage = ReductionStage(org_img._header["DRPSTAGE"]) log.info( "target frame parameters: " f"MJD = {org_img._header['MJD']}, " @@ -4425,7 +4447,7 @@ def detrend_frame( rdnoise = quad.getHdrValue(f"AMP{i+1} RDNOISE") # non-linearity correction - gain_map = _nonlinearity_correction(ptc_params, gain, quad, iquad=i+1) + gain_map, drpstage = _nonlinearity_correction(ptc_params, gain, quad, iquad=i+1, drpstage=drpstage) # gain-correct quadrant quad *= gain_map # propagate new NaNs to the mask @@ -4436,6 +4458,8 @@ def detrend_frame( log.info(f"median error in quadrant {i+1}: {numpy.nanmedian(quad._error):.2f} (e-)") bcorr_img.setHdrValue("BUNIT", "electron", "physical units of the image") + drpstage += "GAIN_CORRECTED" + drpstage += "POISSON_ERROR_CALCULATED" else: # convert to ADU log.info("leaving original ADU units") @@ -4450,6 +4474,7 @@ def detrend_frame( detrended_img = (bcorr_img - mdark_img.convertUnit(to=bcorr_img._header["BUNIT"])) # NOTE: this is a hack to avoid the error propagation of the division in Image detrended_img._data = detrended_img._data / numpy.nan_to_num(mflat_img._data, nan=1.0) + drpstage += "DETRENDED" # propagate pixel mask log.info("propagating pixel mask") @@ -4464,6 +4489,7 @@ def detrend_frame( rdnoise = detrended_img.getHdrValue("AMP1 RDNOISE") detrended_img.reject_cosmics(gain=1.0, rdnoise=rdnoise, rlim=1.3, iterations=5, fwhm_gauss=[2.75, 2.75], replace_box=[10,2], replace_error=1e6, verbose=True, inplace=True) + drpstage += "COSMIC_CLEANED" # replace masked pixels with NaNs if replace_with_nan: @@ -4494,6 +4520,7 @@ def detrend_frame( # save detrended image log.info(f"writing detrended image to '{os.path.basename(out_image)}'") + detrended_img._header["DRPSTAGE"] = drpstage.value detrended_img.writeFitsData(out_image) # show plots diff --git a/python/lvmdrp/functions/rssMethod.py b/python/lvmdrp/functions/rssMethod.py index 34d7c109..d7f892f9 100644 --- a/python/lvmdrp/functions/rssMethod.py +++ b/python/lvmdrp/functions/rssMethod.py @@ -23,6 +23,7 @@ from scipy import interpolate, ndimage from lvmdrp.utils.decorators import skip_on_missing_input_path, skip_if_drpqual_flags +from lvmdrp.utils.bitmask import ReductionStage from lvmdrp.core.constants import CONFIG_PATH, ARC_LAMPS from lvmdrp.core.cube import Cube from lvmdrp.core.tracemask import TraceMask @@ -741,6 +742,7 @@ def shift_wave_skylines(in_frame: str, out_frame: str, dwave: float = 8.0, skyli # write updated wobject log.info(f"writing updated wobject file '{os.path.basename(out_frame)}'") + lvmframe._header["DRPSTAGE"] = (ReductionStage(lvmframe._header["DRPSTAGE"]) + "WAVELENGTH_SHIFTED").value lvmframe.writeFitsData(out_frame) # Make QA plots showing offsets for each sky line in each spectrograph @@ -804,6 +806,7 @@ def create_pixel_table(in_rss: str, out_rss: str, in_waves: str, in_lsfs: str): lsf_traces = [TraceMask.from_file(in_lsfs) for in_lsfs in in_lsfs] lsf_trace = TraceMask.from_spectrographs(*lsf_traces) rss.set_lsf_trace(lsf_trace) + rss._header["DRPSTAGE"] = (ReductionStage(rss._header["DRPSTAGE"]) + "WAVELENGTH_CALIBRATED").value rss.writeFitsData(out_rss) return rss @@ -1630,6 +1633,7 @@ def apply_fiberflat(in_rss: str, out_frame: str, in_flat: str, clip_below: float superflat=flat._data ) lvmframe.set_header(orig_header=rss._header, flatname=os.path.basename(in_flat), ifibvar=ifibvar, ffibvar=ffibvar) + lvmframe._header["DRPSTAGE"] = (ReductionStage(lvmframe._header["DRPSTAGE"]) + "FLATFIELDED").value lvmframe.writeFitsData(out_frame) return rss, lvmframe @@ -2928,6 +2932,7 @@ def stack_spectrographs(in_rsss: List[str], out_rss: str) -> RSS: # write output log.info(f"writing stacked RSS to {os.path.basename(out_rss)}") + rss_out._header["DRPSTAGE"] = (ReductionStage(rss_out._header["DRPSTAGE"]) + "SPECTROGRAPH_STACKED").value rss_out.writeFitsData(out_rss) return rss_out @@ -2971,6 +2976,8 @@ def join_spec_channels(in_fframes: List[str], out_cframe: str, use_weights: bool sky_west=new_rss._sky_west, sky_west_error=new_rss._sky_west_error, slitmap=new_rss._slitmap) + cframe._header["DRPSTAGE"] = (ReductionStage(cframe._header["DRPSTAGE"]) + "CHANNEL_COMBINED").value + # write output RSS if out_cframe is not None: log.info(f"writing output RSS to {os.path.basename(out_cframe)}") diff --git a/python/lvmdrp/functions/skyMethod.py b/python/lvmdrp/functions/skyMethod.py index 074aeee4..124faaf1 100644 --- a/python/lvmdrp/functions/skyMethod.py +++ b/python/lvmdrp/functions/skyMethod.py @@ -29,6 +29,7 @@ SKYMODEL_CONFIG_PATH, SKYMODEL_INST_PATH, ) +from lvmdrp.utils.bitmask import ReductionStage from lvmdrp.core.plot import plt, create_subplots, save_fig from lvmdrp.core.header import Header from lvmdrp.core.passband import PassBand @@ -1455,6 +1456,7 @@ def interpolate_sky( in_frame: str, out_rss: str = None, display_plots: bool = F # write output RSS log.info(f"writing output RSS file '{os.path.basename(out_rss)}'") + new_rss._header["DRPSTAGE"] = (ReductionStage(new_rss._header["DRPSTAGE"]) + "SKY_SUPERSAMPLED").value new_rss.writeFitsData(out_rss) return new_rss, supersky, supererror, swave, ssky, svars, smask @@ -1583,6 +1585,7 @@ def combine_skies(in_rss: str, out_rss, sky_weights: Tuple[float, float] = None) rss._sky_west[std_idx] *= exptime_factors rss._sky_west_error[std_idx] *= exptime_factors + rss._header["DRPSTAGE"] = (ReductionStage(rss._header["DRPSTAGE"]) + "SKY_TELESCOPES_COMBINED").value rss.writeFitsData(out_rss) return rss, sky @@ -1658,6 +1661,7 @@ def quick_sky_subtraction(in_cframe, out_sframe, sframe = lvmSFrame(data=skysub_data, error=skysub_error, mask=cframe._mask.astype(bool), sky=skydata, sky_error=sky_error, wave=cframe._wave, lsf=cframe._lsf, header=cframe._header, slitmap=cframe._slitmap) sframe._mask |= ~np.isfinite(sframe._error) + sframe._header["DRPSTAGE"] = (ReductionStage(sframe._header["DRPSTAGE"]) + "SKY_SUBTRACTED").value sframe.writeFitsData(out_sframe) # TODO: check on expnum=7632 for halpha emission in sky fibers diff --git a/python/lvmdrp/utils/bitmask.py b/python/lvmdrp/utils/bitmask.py index eb2aa0a9..ac98f008 100644 --- a/python/lvmdrp/utils/bitmask.py +++ b/python/lvmdrp/utils/bitmask.py @@ -141,22 +141,54 @@ def __add__(self, flag): class ReductionStage(BaseBitmask): # completed reduction steps - UNREDUCED = auto() # exposure not reduced - PREPROCESSED = auto() # trimmed overscan region - CALIBRATED = auto() # bias, dark and pixelflat calibrated - COSMIC_CLEAN = auto() # cosmic ray cleaned - STRAY_CLEAN = auto() # stray light subtracted - FIBERS_FOUND = auto() # fiberflat fibers located along the column - FIBERS_TRACED = auto() # fiberflat fibers traces along the dispersion axis - SPECTRA_EXTRACTED = ( - auto() - ) # extracted the fiber spectra of any arc, flat, or science frames - WAVELENGTH_SOLVED = auto() # arc fiber wavelength solution found - WAVELENGTH_RESAMPLED = ( - auto() - ) # fiber wavelength resampled to common wavelength/LSF vector + UNREDUCED = auto() # exposure in raw stage + HDRFIX_APPLIED = auto() # header fix applied + OVERSCAN_SUBTRACTED = auto() # trimmed overscan region, overscan-subtracted + PIXELMASK_ADDED = auto() # fiducial pixel mask was added + GAIN_CORRECTED = auto() # gain correction applied + POISSON_ERROR_CALCULATED = auto() # calculated poisson error + LINEARITY_CORRECTED = auto() # linearity correction applied + DETRENDED = auto() # bias subtracted and pixel flat-fielded + COSMIC_CLEANED = auto() # cosmic ray cleaned + SCI_ASTROMETRY_ADDED = auto() # added astrometric solution for science telescope + SPEC_ASTROMETRY_ADDED = auto() # added astrometric solution for spectrophotometric telescope + SKYE_ASTROMETRY_ADDED = auto() # added astrometric solution for sky east telescope + SKYW_ASTROMETRY_ADDED = auto() # added astrometric solution for sky west telescope + FIBERS_ASTROMETRY_ADDED = auto() # added astrometric solution for all fibers in the system + STRAYLIGHT_SUBTRACTED = auto() # stray light subtracted + FIBERS_SHIFTED = auto() # fibers positions corrected for thermal shifts + SPECTRA_EXTRACTED = auto() # extracted spectra + SPECTROGRAPH_STACKED = auto() # stacked spectrograph wise + WAVELENGTH_CALIBRATED = auto() # arc fiber wavelength solution found + FLATFIELDED = auto() # fiber flat-fielded + WAVELENGTH_SHIFTED = auto() # wavelength corrected for thermal shifts + SKY_SUPERSAMPLED = auto() # extrapolated sky fibers along slit + SKY_TELESCOPES_COMBINED = auto() # telescope-combined extrapolated sky fibers + WAVELENGTH_RECTIFIED = auto() # wavelength resampled to common grid along slit + MEASURED_SENS_STD = auto() # measured sensitivity curve using standard stars + MEASURED_SENS_SCI = auto() # scaled fiducial sensitivity using science field stars + FLUX_CALIBRATED = auto() # flux calibrated all fibers in the system + CHANNEL_COMBINED = auto() # channel stitched together + SKY_SUBTRACTED = auto() # sky-subtracted all fibers in the system + def __add__(self, flag): + if isinstance(flag, self.__class__): + pass + elif isinstance(flag, str): + flag = self.__class__[flag.upper()] + elif isinstance(flag, int): + flag = self.__class__(flag) + else: + try: + return super().__add__(flag) + except: + raise + + new = copy(self) + new = self ^ self.__class__["UNREDUCED"] + return new | flag + class QualityFlag(BaseBitmask): # TODO: add flag for overscan quality OSFEATURES = auto() # Overscan region has features. From 5247dc81d4fc2b6e651c991f377b56dfee830eae Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Mon, 9 Sep 2024 14:39:22 -0300 Subject: [PATCH 02/20] updating pixel level bitmask & added comments --- python/lvmdrp/utils/bitmask.py | 85 +++++++++++++++++----------------- 1 file changed, 43 insertions(+), 42 deletions(-) diff --git a/python/lvmdrp/utils/bitmask.py b/python/lvmdrp/utils/bitmask.py index ac98f008..76345d83 100644 --- a/python/lvmdrp/utils/bitmask.py +++ b/python/lvmdrp/utils/bitmask.py @@ -98,6 +98,7 @@ def __contains__(self, flag): class RawFrameQuality(BaseBitmask): + # TODO: repurpose this to use Dmitry's QC flags GOOD = auto() # bit whether a raw frame is good for reduction TEST = auto() # bit whether a raw frame is for instrument testing purposes BAD = auto() # bit whether a raw frame is bad for reduction @@ -218,63 +219,63 @@ class QualityFlag(BaseBitmask): class PixMask(BaseBitmask): # fiber bitmasks - NOPLUG = auto() + NONEXPOSED = auto() + WEAKFIBER = auto() + DEADFIBER = auto() + INTERPOLATED = auto() + + # measure quality of tracing using polynomial fit - samples residuals BADTRACE = auto() - BADFLAT = auto() BADARC = auto() - MANYBADCOLUMNS = auto() - MANYREJECTED = auto() + + # measure offset of the fiber from the median flatfielded fiber + BADFLAT = auto() + + # offset from a preset fiber shift value LARGESHIFT = auto() + + BADSTDFIBER = auto() BADSKYFIBER = auto() - NEARWHOPPER = auto() - WHOPPER = auto() - SMEARIMAGE = auto() - SMEARHIGHSN = auto() - SMEARMEDSN = auto() - DEADFIBER = auto() + # pixel bitmasks + + # pixels with no useful information + NODATA = auto() + + # TODO: bright pixels on top and bottom edges + + # set this if X% close to saturation level SATURATION = auto() + + # from pixelmasks BADPIX = auto() - COSMIC = auto() NEARBADPIXEL = auto() + + # from CR rejection + COSMIC = auto() + + # outlying in pixelflat? possible dust spec LOWFLAT = auto() - FULLREJECT = auto() - PARTIALREJECT = auto() - SCATTEREDLIGHT = auto() + + # clipped straylight polynomial fit + STRAYLIGHT = auto() + CROSSTALK = auto() + + # missing sky lines? NOSKY = auto() + # too bright sky lines? BRIGHTSKY = auto() - NODATA = auto() - COMBINEREJ = auto() + + # pixels with sensitivities too deviant from instrumental sensitivity BADFLUXFACTOR = auto() + + # large sky residuals BADSKYCHI = auto() # define flag name constants -STATUS = list(ReductionStatus.__members__.keys()) -STAGE = list(ReductionStage.__members__.keys()) +# RAW_QUALITIES = list(RawFrameQuality.__members__.keys()) +STAGES = list(ReductionStage.__members__.keys()) FLAGS = list(QualityFlag.__members__.keys()) - -if __name__ == "__main__": - status = ReductionStatus(0) - stage = ReductionStage.UNREDUCED - print(status.get_name(), stage.get_name()) - status += "IN_PROGRESS" - print(status.get_name(), stage.get_name()) - stage += ReductionStage.PREPROCESSED - print(status.get_name(), stage.get_name()) - stage += ReductionStage.CALIBRATED - print(status.get_name(), stage.get_name()) - status += ReductionStatus.FINISHED - print(status.get_name(), stage.get_name()) - status += ReductionStatus.IN_PROGRESS - print(status.get_name(), stage.get_name()) - stage += ReductionStage.FIBERS_FOUND - print(status.get_name(), stage.get_name()) - status += ReductionStatus.FINISHED - print(status.get_name(), stage.get_name()) - print("finished" in status) - status = ReductionStatus.IN_PROGRESS - print(status.get_name(), stage.get_name()) - stage += ReductionStage.PREPROCESSED | ReductionStage.CALIBRATED - print(status.get_name(), stage.get_name()) +DRPQUALITIES = list() From 25a7a45c6cabcef373e5341b35bb9f81db47379a Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Mon, 9 Sep 2024 14:39:52 -0300 Subject: [PATCH 03/20] added first Image methods to handle pixel level bitmasks --- python/lvmdrp/core/image.py | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/python/lvmdrp/core/image.py b/python/lvmdrp/core/image.py index fb1bc1d9..b280cf8a 100644 --- a/python/lvmdrp/core/image.py +++ b/python/lvmdrp/core/image.py @@ -17,6 +17,7 @@ from scipy import interpolate from lvmdrp import log +from lvmdrp.utils.bitmask import PixMask from lvmdrp.core.constants import CON_LAMPS, ARC_LAMPS from lvmdrp.core.plot import plt from lvmdrp.core.fit_profile import gaussians, Gaussians @@ -914,7 +915,7 @@ def getData(self): """ return self._data - def getMask(self): + def get_mask(self, as_boolean=False): """ Returns the bad pixel mask of the image @@ -924,6 +925,8 @@ def getMask(self): The bad pixel mask of the image """ + if as_boolean: + return self._mask.astype(bool) return self._mask def getError(self): @@ -999,6 +1002,12 @@ def getSlice(self, slice, axis): numpy.arange(self._dim[0]), self._data[:, slice], error, mask ) + def set_mask(self, pixmask): + if isinstance(pixmask, numpy.ndarray): + assert isinstance(pixmask[0,0], PixMask) + + self._mask = pixmask + def setData( self, data=None, error=None, mask=None, header=None, select=None, inplace=True ): @@ -1303,7 +1312,7 @@ def loadFitsData( if hdu[i].header["EXTNAME"].split()[0] == "ERROR": self._error = hdu[i].data.astype("float32") elif hdu[i].header["EXTNAME"].split()[0] == "BADPIX": - self._mask = hdu[i].data.astype("bool") + self._mask = hdu[i].data.astype("int32") elif hdu[i].header["EXTNAME"].split()[0] == "FRAMES": self._individual_frames = Table(hdu[i].data) elif hdu[i].header["EXTNAME"].split()[0] == "SLITMAP": @@ -1315,7 +1324,7 @@ def loadFitsData( self._dim = self._data.shape # set dimension if extension_mask is not None: - self._mask = hdu[extension_mask].data.astype("bool") # take data + self._mask = hdu[extension_mask].data.astype("int32") # take data self._dim = self._mask.shape # set dimension if extension_error is not None: @@ -1371,7 +1380,7 @@ def writeFitsData( if self._error is not None: hdus[1] = pyfits.ImageHDU(self._error, name="ERROR") if self._mask is not None: - hdus[2] = pyfits.ImageHDU(self._mask.astype("uint8"), name="BADPIX") + hdus[2] = pyfits.ImageHDU(self._mask.astype("uint32"), name="BADPIX") if self._individual_frames is not None: hdus[3] = pyfits.BinTableHDU(self._individual_frames, name="FRAMES") if self._slitmap is not None: @@ -1384,10 +1393,10 @@ def writeFitsData( # mask hdu if extension_mask == 0: - hdu = pyfits.PrimaryHDU(self._mask.astype("uint8")) + hdu = pyfits.PrimaryHDU(self._mask.astype("uint32")) elif extension_mask > 0 and extension_mask is not None: hdus[extension_mask] = pyfits.ImageHDU( - self._mask.astype("uint8"), name="BADPIX" + self._mask.astype("uint32"), name="BADPIX" ) # error hdu @@ -1924,7 +1933,7 @@ def fitPoly(self, axis="y", order=4, plot=-1): x = x - bn.nanmean(x) # if self._mask is not None: # self._mask = numpy.logical_and(self._mask, numpy.logical_not(numpy.isnan(self._data))) - valid = ~self._mask.astype("bool") + valid = ~self._mask # iterate over the image for i in range(slices): # decide on the bad pixel mask From f4eb0ae9c14ad14e294c6c3e4511c721eaa914cd Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Mon, 9 Sep 2024 17:56:57 -0300 Subject: [PATCH 04/20] adding pixel bitmask for saturated pixels --- python/lvmdrp/utils/bitmask.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/python/lvmdrp/utils/bitmask.py b/python/lvmdrp/utils/bitmask.py index 76345d83..bb152f0d 100644 --- a/python/lvmdrp/utils/bitmask.py +++ b/python/lvmdrp/utils/bitmask.py @@ -251,6 +251,9 @@ class PixMask(BaseBitmask): BADPIX = auto() NEARBADPIXEL = auto() + # from raw data preprocessing + SATURATED = auto() + # from CR rejection COSMIC = auto() From 5e8744935018e189767dfb86a2907d79462fcfe9 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Mon, 9 Sep 2024 19:09:09 -0300 Subject: [PATCH 05/20] removing pixel mask refinement option --- python/lvmdrp/functions/imageMethod.py | 9 --------- 1 file changed, 9 deletions(-) diff --git a/python/lvmdrp/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index 574f2a84..350ecc97 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -4335,7 +4335,6 @@ def detrend_frame( calculate_error: bool = True, replace_with_nan: bool = True, reject_cr: bool = True, - median_box: list = [0, 0], display_plots: bool = False, ): """detrends input image by subtracting bias, dark and flatfielding @@ -4362,8 +4361,6 @@ def detrend_frame( whether to replace or not NaN values by zeros, by default True reject_cr : bool, optional whether to reject or not cosmic rays from detrended image, by default True - median_box : tuple, optional - size of the median box to refine pixel mask, by default [0,0] display_plots : str, optional whether to show plots on display or not, by default False """ @@ -4499,12 +4496,6 @@ def detrend_frame( log.info(f"replacing {detrended_img._mask.sum()} masked pixels with NaNs") detrended_img.apply_pixelmask() - # refine mask - if all(median_box): - log.info(f"refining pixel mask with {median_box = }") - med_img = detrended_img.medianImg(size=median_box, use_mask=True) - detrended_img.setData(mask=(detrended_img._mask | med_img._mask), inplace=True) - # normalize in case of pixel flat calibration # 'pixflat' is the imagetyp that a pixel flat can have if img_type == "pixflat": From f70cc89feb97c4ffd18ece818497dbe137c8e905 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Mon, 9 Sep 2024 19:10:35 -0300 Subject: [PATCH 06/20] implementing pixel bitmasks in 2D reduction stage --- python/lvmdrp/core/image.py | 70 ++++++++++++++------------ python/lvmdrp/functions/imageMethod.py | 22 ++++---- 2 files changed, 49 insertions(+), 43 deletions(-) diff --git a/python/lvmdrp/core/image.py b/python/lvmdrp/core/image.py index b280cf8a..0502d84e 100644 --- a/python/lvmdrp/core/image.py +++ b/python/lvmdrp/core/image.py @@ -446,7 +446,7 @@ def __add__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - new_mask = numpy.logical_or(self._mask, other._mask) + new_mask = numpy.where(self._mask != other._mask, self._mask + other._mask, self._mask) img.setData(mask=new_mask) else: img.setData(mask=self._mask) @@ -511,7 +511,7 @@ def __sub__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - new_mask = numpy.logical_or(self._mask, other._mask) + new_mask = numpy.where(self._mask != other._mask, self._mask + other._mask, self._mask) img.setData(mask=new_mask) else: img.setData(mask=self._mask) @@ -572,7 +572,7 @@ def __truediv__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - new_mask = numpy.logical_or(self._mask, other._mask) + new_mask = numpy.where(self._mask != other._mask, self._mask + other._mask, self._mask) img.setData(mask=new_mask) else: img.setData(mask=self._mask) @@ -651,7 +651,7 @@ def __mul__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - new_mask = numpy.logical_or(self._mask, other._mask) + new_mask = numpy.where(self._mask != other._mask, self._mask + other._mask, self._mask) img.setData(mask=new_mask) else: img.setData(mask=self._mask) @@ -761,7 +761,7 @@ def measure_fiber_shifts(self, ref_image, columns=[500, 1000, 1500, 2000, 2500, def apply_pixelmask(self, mask=None): """Applies the mask to the data and error arrays, setting to nan when True and leaving the same value otherwise""" if mask is None: - mask = self._mask + mask = self.get_mask(as_boolean=True) if mask is None: return self._data, self._error @@ -1486,19 +1486,22 @@ def replaceMaskMedian(self, box_x, box_y, replace_error=1e20): new_image : Image object Subsampled image """ - if self._data is None: raise RuntimeError("Image object is empty. Nothing to process.") + if self._mask is None: + return self + + mask = self.get_mask(as_boolean=True) idx = numpy.indices(self._dim) # create an index array # get x and y coordinates of bad pixels - y_cors = idx[0][self._mask] - x_cors = idx[1][self._mask] + y_cors = idx[0][mask] + x_cors = idx[1][mask] out_data = copy(self._data) msk_data = copy(self._data) - msk_data[self._mask] = numpy.nan + msk_data[mask] = numpy.nan out_error = copy(self._error) # esimate the pixel distance form the bad pixel to the filter window boundary @@ -1625,7 +1628,7 @@ def subsampleImg(self): else: new_error = None if self._mask is not None: - new_mask = numpy.zeros(new_dim, dtype="bool") + new_mask = numpy.zeros(new_dim, dtype=int) else: new_mask = None @@ -1648,10 +1651,11 @@ def subsampleImg(self): new_error[select3] = self._error.flatten() new_error[select4] = self._error.flatten() if self._mask is not None: - new_mask[select1] = self._mask.flatten() - new_mask[select2] = self._mask.flatten() - new_mask[select3] = self._mask.flatten() - new_mask[select4] = self._mask.flatten() + mask = self.get_mask(as_boolean=True) + new_mask[select1] = mask.flatten() + new_mask[select2] = mask.flatten() + new_mask[select3] = mask.flatten() + new_mask[select4] = mask.flatten() # create new Image object with the new subsample data new_image = copy(self) new_image.setData(data=new_data, error=new_error, mask=new_mask, inplace=True) @@ -1714,7 +1718,7 @@ def rebin(self, bin_x, bin_y): 1, ) # if only one bad pixel in the binning pixel exists the binned pixel will have the bad pixel status - new_mask = mask_new2 > 0 + new_mask = mask_new2# > 0 else: new_mask = None # create new Image object and return @@ -1760,7 +1764,7 @@ def convolveImg(self, kernel, mode="nearest"): new_image.setData(data=new, error=new_error, inplace=True) return new_image - def convolveGaussImg(self, sigma_x, sigma_y, mode="nearest", mask=False): + def convolveGaussImg(self, sigma_x, sigma_y, mode="nearest", use_mask=False): """ Convolves the data of the Image with a given kernel. The mask and error information will be unchanged. @@ -1781,17 +1785,18 @@ def convolveGaussImg(self, sigma_x, sigma_y, mode="nearest", mask=False): """ # convolve the data array with the 2D Gaussian convolution kernel - if self._mask is not None and mask is True: - mask_data = self._data[self._mask] - self._data[self._mask] = 0 + if self._mask is not None and use_mask is True: + mask = self.get_mask(as_boolean=True) + mask_data = self._data[mask] + self._data[mask] = 0 gauss = ndimage.filters.gaussian_filter( self._data, (sigma_y, sigma_x), mode=mode ) scale = ndimage.filters.gaussian_filter( - (~self._mask).astype('float32'), (sigma_y, sigma_x), mode=mode + (~mask).astype('float32'), (sigma_y, sigma_x), mode=mode ) new = gauss / scale - self._data[self._mask] = mask_data + self._data[mask] = mask_data else: new = ndimage.filters.gaussian_filter( self._data, (sigma_y, sigma_x), mode=mode @@ -1826,8 +1831,9 @@ def medianImg(self, size, mode="nearest", use_mask=False, propagate_error=False) median filtered image """ if self._mask is None and use_mask: + mask = self.get_mask(as_boolean=True) new_data = copy(self._data) - new_data[self._mask] = numpy.nan + new_data[mask] = numpy.nan new_data = ndimage.median_filter(new_data, size, mode=mode) new_mask = None new_error = None @@ -1840,9 +1846,10 @@ def medianImg(self, size, mode="nearest", use_mask=False, propagate_error=False) if propagate_error and self._error is not None: new_error = numpy.sqrt(ndimage.median_filter(self._error ** 2, size, mode=mode)) else: + mask = self.get_mask(as_boolean=True) # copy data and replace masked with nans new_data = copy(self._data) - new_data[self._mask] = numpy.nan + new_data[mask] = numpy.nan # perform median filter new_data = signal.medfilt2d(new_data, size) # update mask @@ -1853,7 +1860,7 @@ def medianImg(self, size, mode="nearest", use_mask=False, propagate_error=False) new_error = None if propagate_error and self._error is not None: new_error = copy(self._error) - new_error[self._mask] = numpy.nan + new_error[mask] = numpy.nan new_error = numpy.sqrt(signal.medfilt2d(new_error ** 2, size)) # reset masked errors in new array new_error[new_mask] = self._error[new_mask] @@ -3095,8 +3102,8 @@ def reject_cosmics(self, sigma_det=5, rlim=1.2, iterations=5, fwhm_gauss=[2.0,2. img_original._error = numpy.sqrt((numpy.clip(img_original._data, a_min=0.0, a_max=None) + rdnoise**2)) select = numpy.zeros(img._dim, dtype=bool) - img_original._mask = numpy.zeros(img._dim, dtype=bool) - img._mask = numpy.zeros(img._dim, dtype=bool) + img_original._mask = numpy.zeros(img._dim, dtype=int) + img._mask = numpy.zeros(img._dim, dtype=int) # start iteration out = img @@ -3119,7 +3126,7 @@ def reject_cosmics(self, sigma_det=5, rlim=1.2, iterations=5, fwhm_gauss=[2.0,2. S_prime = S-S.medianImg((5, 5)) # cleaning of the normalized Laplacian image # Perform additional clean using a 2D Gaussian smoothing kernel - fine = out.convolveGaussImg(sigma_x, sigma_y, mask=True) # convolve image with a 2D Gaussian + fine = out.convolveGaussImg(sigma_x, sigma_y, use_mask=True) # convolve image with a 2D Gaussian fine_norm = out/fine select_neg = fine_norm < 0 fine_norm.replace_subselect(select_neg, data=0) @@ -3135,7 +3142,7 @@ def reject_cosmics(self, sigma_det=5, rlim=1.2, iterations=5, fwhm_gauss=[2.0,2. log.info(f' Total number of detected cosmics: {det_pix} out of {dim[0] * dim[1]} pixels') if i == iterations-1: - img_original.replace_subselect(select, mask=True) # set the new mask + img_original.replace_subselect(select, mask=PixMask["COSMIC"]) # set the new mask if increase_radius > 0: mask_img = Image(data=img_original._mask) mask_new = mask_img.convolveImg(kernel=numpy.ones((2*increase_radius+1, 2*increase_radius+1))) @@ -3147,14 +3154,14 @@ def reject_cosmics(self, sigma_det=5, rlim=1.2, iterations=5, fwhm_gauss=[2.0,2. # leave out the dispersion direction (0 degrees), see DESI, Guy et al., ApJ, 2023, 165, 144 lse = LinearSelectionElement(11, 11, ang) bc_mask = bc_mask | ndimage.binary_closing(bmask, structure=lse.se) - img_original._mask = bc_mask + img_original._mask = bc_mask * PixMask["COSMIC"] if verbose: log.info(f' Total number after binary closing: {numpy.sum(bc_mask)} pixels') # replace possible corrput pixel with median for final output out = img_original.replaceMaskMedian(box_x, box_y, replace_error=replace_error) else: - out.replace_subselect(select, mask=True) # set the new mask + out.replace_subselect(select, mask=PixMask["COSMIC"]) # set the new mask out = out.replaceMaskMedian(box_x, box_y, replace_error=None) # replace possible corrput pixel with median if inplace: @@ -3166,10 +3173,11 @@ def reject_cosmics(self, sigma_det=5, rlim=1.2, iterations=5, fwhm_gauss=[2.0,2. if self._mask is None: self._mask = out._mask else: - self._mask |= out._mask + self._mask += out._mask else: return out + def getIndividualFrames(self): return self._individual_frames diff --git a/python/lvmdrp/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index 350ecc97..17362fea 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -28,7 +28,7 @@ from lvmdrp import log, __version__ as DRPVER from lvmdrp.core.constants import CONFIG_PATH, SPEC_CHANNELS, ARC_LAMPS, LVM_REFERENCE_COLUMN, LVM_NBLOCKS, FIDUCIAL_PLATESCALE from lvmdrp.utils.decorators import skip_on_missing_input_path, drop_missing_input_paths -from lvmdrp.utils.bitmask import QualityFlag, ReductionStage +from lvmdrp.utils.bitmask import QualityFlag, ReductionStage, PixMask from lvmdrp.core.fiberrows import FiberRows, _read_fiber_ypix from lvmdrp.core.image import ( Image, @@ -1023,7 +1023,7 @@ def detCos_drp( noise = (noise + rdnoise**2).sqrt() result = [] if cpus > 1: - fine = out.convolveGaussImg(sigma, sigma, mask=True) + fine = out.convolveGaussImg(sigma, sigma, use_mask=True) fine_norm = out / fine select_neg = fine_norm < 0 fine_norm.setData(data=0, select=select_neg) @@ -1072,7 +1072,7 @@ def detCos_drp( (5, 5) ) # cleaning of the normalized Laplacian image fine = out.convolveGaussImg( - sigma, sigma, mask=True + sigma, sigma, use_mask=True ) # convolve image with a 2D Gaussian fine_norm = out / fine @@ -4023,13 +4023,13 @@ def preproc_raw_frame( # create pixel mask on the original image log.info("building pixel mask") - proc_img._mask = master_mask + proc_img._mask = master_mask * PixMask["BADPIX"] # convert temp image to ADU for saturated pixel masking saturated_mask = proc_img._data >= 2**16 - proc_img._mask |= saturated_mask + proc_img._mask += saturated_mask * PixMask["SATURATED"] # log number of masked pixels - nmasked = proc_img._mask.sum() + nmasked = numpy.count_nonzero(proc_img._mask) log.info(f"{nmasked} ({nmasked / proc_img._mask.size * 100:.2g} %) pixels masked") # update masked pixels with NaNs if needed @@ -4451,7 +4451,7 @@ def detrend_frame( # gain-correct quadrant quad *= gain_map # propagate new NaNs to the mask - quad._mask |= numpy.isnan(quad._data) + quad._mask += numpy.isnan(quad._data) * PixMask["BADPIX"] quad.computePoissonError(rdnoise) bcorr_img.setSection(section=quad_sec, subimg=quad, inplace=True) @@ -4478,10 +4478,8 @@ def detrend_frame( # propagate pixel mask log.info("propagating pixel mask") - nanpixels = numpy.isnan(detrended_img._data) - infpixels = numpy.isinf(detrended_img._data) - detrended_img._mask = numpy.logical_or(org_img._mask, nanpixels) - detrended_img._mask = numpy.logical_or(detrended_img._mask, infpixels) + nanpixels = ~numpy.isfinite(detrended_img._data) + detrended_img._mask += numpy.where(detrended_img._mask != nanpixels * PixMask["BADPIX"], detrended_img._mask + nanpixels * PixMask["BADPIX"], detrended_img._mask) # reject cosmic rays if reject_cr: @@ -4493,7 +4491,7 @@ def detrend_frame( # replace masked pixels with NaNs if replace_with_nan: - log.info(f"replacing {detrended_img._mask.sum()} masked pixels with NaNs") + log.info(f"replacing {numpy.count_nonzero(detrended_img._mask)} masked pixels with NaNs") detrended_img.apply_pixelmask() # normalize in case of pixel flat calibration From 1ee43de02860776deaccd908c30aed2cf0fb4f15 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Wed, 11 Sep 2024 12:59:29 -0300 Subject: [PATCH 07/20] better propagation of pixel level bitmasks --- python/lvmdrp/core/image.py | 62 +++++++++++++++++++++++++++------- python/lvmdrp/utils/bitmask.py | 55 ++++++++++++++++++++++++++++++ 2 files changed, 105 insertions(+), 12 deletions(-) diff --git a/python/lvmdrp/core/image.py b/python/lvmdrp/core/image.py index 0502d84e..4436ab8e 100644 --- a/python/lvmdrp/core/image.py +++ b/python/lvmdrp/core/image.py @@ -17,7 +17,7 @@ from scipy import interpolate from lvmdrp import log -from lvmdrp.utils.bitmask import PixMask +from lvmdrp.utils.bitmask import PixMask, add_bitmask, toggle_bitmask from lvmdrp.core.constants import CON_LAMPS, ARC_LAMPS from lvmdrp.core.plot import plt from lvmdrp.core.fit_profile import gaussians, Gaussians @@ -425,7 +425,7 @@ def __add__(self, other): """ if isinstance(other, Image): # define behaviour if the other is of the same instance - img = Image(header=self._header, origin=self._origin) + img = copy(self) # add data if contained in both if self._data is not None and other._data is not None: @@ -446,8 +446,7 @@ def __add__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - new_mask = numpy.where(self._mask != other._mask, self._mask + other._mask, self._mask) - img.setData(mask=new_mask) + img.setData(mask=self._mask | other._mask) else: img.setData(mask=self._mask) return img @@ -511,8 +510,7 @@ def __sub__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - new_mask = numpy.where(self._mask != other._mask, self._mask + other._mask, self._mask) - img.setData(mask=new_mask) + img.setData(mask=self._mask | other._mask) else: img.setData(mask=self._mask) return img @@ -572,8 +570,7 @@ def __truediv__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - new_mask = numpy.where(self._mask != other._mask, self._mask + other._mask, self._mask) - img.setData(mask=new_mask) + img.setData(mask=self._mask | other._mask) else: img.setData(mask=self._mask) return img @@ -651,8 +648,7 @@ def __mul__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - new_mask = numpy.where(self._mask != other._mask, self._mask + other._mask, self._mask) - img.setData(mask=new_mask) + img.setData(mask=self._mask | other._mask) else: img.setData(mask=self._mask) return img @@ -1067,6 +1063,48 @@ def setData( return new_image + def add_bitmask(self, pixmask, where=None): + """Adds the given bitmask to the current mask + + If no mask is present in the image, a new mask will be started + with the corresponding pixels flagged with `pixmask`. + + Parameters + ---------- + pixmask : PixMask|str|int + Bitmask to add to the current mask + where : numpy.ndarray[bool], optional + Boolean selection of pixels where to add given bitmask, by default all pixels + + Returns + ------- + mask : np.ndarray[PixMask|int] + Updated pixel mask + """ + self._mask = add_bitmask(self._mask, pixmask=pixmask, where=where) + return self._mask + + def toggle_bitmask(self, pixmask, where=None): + """Toggle the given bitmask in the current mask + + If no mask is present in the image, a new mask will be started + with the corresponding pixels flagged with `pixmask`. + + Parameters + ---------- + pixmask : PixMask|str|int + Bitmask to add to the current mask + where : numpy.ndarray[bool], optional + Boolean selection of pixels where to add given bitmask, by default all pixels + + Returns + ------- + mask : np.ndarray[PixMask|int] + Updated pixel mask + """ + self._mask = toggle_bitmask(self._mask, pixmask=pixmask, where=where) + return self._mask + def convertUnit(self, to, assume="adu", gain_field="GAIN", inplace=False): """converts the unit of the image @@ -1707,11 +1745,11 @@ def rebin(self, bin_x, bin_y): if self._mask is not None: # create the new bad pixel mask - mask_new = numpy.sum( + mask_new = numpy.bitwise_or.reduce( numpy.reshape(self._mask, (self._dim[0], self._dim[1] // bin_x, bin_x)), 2, ) - mask_new2 = numpy.sum( + mask_new2 = numpy.bitwise_or.reduce( numpy.reshape( mask_new, (self._dim[0] // bin_y, bin_y, self._dim[1] // bin_x) ), diff --git a/python/lvmdrp/utils/bitmask.py b/python/lvmdrp/utils/bitmask.py index bb152f0d..ed9a51b4 100644 --- a/python/lvmdrp/utils/bitmask.py +++ b/python/lvmdrp/utils/bitmask.py @@ -8,6 +8,7 @@ from copy import copy from enum import IntFlag, auto +import numpy as np class classproperty(object): @@ -282,3 +283,57 @@ class PixMask(BaseBitmask): STAGES = list(ReductionStage.__members__.keys()) FLAGS = list(QualityFlag.__members__.keys()) DRPQUALITIES = list() + + +def _parse_bitmask(pixmask, mask_shape): + if isinstance(pixmask, PixMask): + pixmask = np.zeros(mask_shape) + elif isinstance(pixmask, str): + pixmask = PixMask[pixmask] + elif isinstance(pixmask, int): + pixmask = PixMask(pixmask) + elif isinstance(pixmask, np.ndarray[int]): + assert pixmask.shape == mask_shape, f"Wrong `pixmask` shape {pixmask.shape} not matching `mask_image` shape {mask_shape}" + else: + raise ValueError(f"Wrong type for {pixmask = }: {type(pixmask)}; expected PixMask, string or integer") + + return pixmask + + +def _parse_where(where, mask_shape): + if where is not None and isinstance(where, np.ndarray[bool]): + assert where.shape == mask_shape, f"Wrong `where` shape {where.shape} not matching `mask_image` shape {mask_shape}" + else: + where = np.ones(mask_shape, dtype=bool) + + return where + + +def add_bitmask(mask_image, pixmask, where=None): + pixmask = _parse_bitmask(pixmask, mask_shape=mask_image.shape) + where = _parse_where(where, mask_shape=mask_image.shape) + + mask = np.zeros_like(mask_image, dtype=int) + mask[where] |= pixmask + + if mask_image is None: + mask_image = mask + return mask_image + + mask_image[where] |= mask[where] + return mask_image + + +def toggle_bitmask(mask_image, pixmask, where=None): + pixmask = _parse_bitmask(pixmask, mask_shape=mask_image.shape) + where = _parse_where(where, mask_shape=mask_image.shape) + + mask = np.zeros_like(mask_image, dtype=int) + mask[where] |= pixmask + + if mask_image is None: + mask_image = mask + return mask_image + + mask_image[where] ^= mask[where] + return mask_image From e43f841e53af59d70a46c32a450467a74a694c9a Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Wed, 11 Sep 2024 19:51:41 -0300 Subject: [PATCH 08/20] implementing bitmasks to spectrum1d, rss and fiberrows modules --- python/lvmdrp/core/fiberrows.py | 77 ++++++++++------- python/lvmdrp/core/image.py | 2 +- python/lvmdrp/core/rss.py | 142 +++++++++++++++---------------- python/lvmdrp/core/spectrum1d.py | 121 +++++++++++++------------- python/lvmdrp/utils/bitmask.py | 3 + 5 files changed, 182 insertions(+), 163 deletions(-) diff --git a/python/lvmdrp/core/fiberrows.py b/python/lvmdrp/core/fiberrows.py index b5c29339..2698ac0d 100644 --- a/python/lvmdrp/core/fiberrows.py +++ b/python/lvmdrp/core/fiberrows.py @@ -7,6 +7,7 @@ from lvmdrp import log from scipy import optimize from astropy.table import Table +from lvmdrp.utils.bitmask import PixMask from lvmdrp.core.header import Header, combineHdr from lvmdrp.core.positionTable import PositionTable from lvmdrp.core.spectrum1d import Spectrum1D, _cross_match_float @@ -183,7 +184,7 @@ def __truediv__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - new_mask = numpy.logical_or(self._mask, other._mask) + new_mask = numpy.bitwise_or(self._mask, other._mask) img.setData(mask=new_mask) else: img.setData(mask=self._mask) @@ -278,7 +279,7 @@ def __add__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - new_mask = numpy.logical_or(self._mask, other._mask) + new_mask = numpy.bitwise_or(self._mask, other._mask) img.setData(mask=new_mask) else: img.setData(mask=self._mask) @@ -398,7 +399,7 @@ def __mul__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - new_mask = numpy.logical_or(self._mask, other._mask) + new_mask = numpy.bitwise_or(self._mask, other._mask) img.setData(mask=new_mask) else: img.setData(mask=self._mask) @@ -538,7 +539,7 @@ def createEmpty(self, data_dim=None, poly_deg=None, samples_columns=None): if data_dim is not None: # create empty mask all pixel assigned bad - self._mask = numpy.ones(data_dim, dtype="bool") + self._mask = numpy.ones(data_dim, dtype=int) if data_dim is not None and samples_columns is not None: self._samples = Table(data=numpy.zeros((data_dim[0], len(samples_columns))) + numpy.nan, names=samples_columns) @@ -681,6 +682,20 @@ def getData(self): mask = self._mask return data, error, mask + def get_mask(self, as_boolean=False): + """ + Returns the bad pixel mask of the image + + Returns + ----------- + _mask : numpy.ndarray + The bad pixel mask of the image + + """ + if as_boolean: + return self._mask.astype(bool) + return self._mask + def setData(self, select=None, data=None, mask=None, error=None): if select is not None: if data is not None: @@ -699,7 +714,7 @@ def setData(self, select=None, data=None, mask=None, error=None): if mask is not None: self._mask = mask nfibers, npixels = self._mask.shape - self._good_fibers = numpy.where(numpy.sum(self._mask, axis=1) != self._mask.shape[1])[0] + self._good_fibers = numpy.where(numpy.any(self._mask != 0, axis=1))[0] elif not hasattr(self, "_mask"): self._mask = None self._good_fibers = None @@ -820,8 +835,8 @@ def loadFitsData( if hdu[i].header["EXTNAME"].split()[0] == "ERROR": self._error = hdu[i].data.astype("float32") elif hdu[i].header["EXTNAME"].split()[0] == "BADPIX": - self._mask = hdu[i].data.astype("bool") - self._good_fibers = numpy.where(numpy.sum(self._mask, axis=1) != self._data.shape[1])[0] + self._mask = hdu[i].data.astype("int32") + self._good_fibers = numpy.where(numpy.any(self._mask!=0, axis=1))[0] elif hdu[i].header["EXTNAME"].split()[0] == "COEFFS": self._coeffs = hdu[i].data.astype("float32") @@ -831,8 +846,8 @@ def loadFitsData( self._fibers = self._data.shape[0] self._pixels = numpy.arange(self._data.shape[1]) if extension_mask is not None: - self._mask = hdu[extension_mask].data.astype("bool") - self._good_fibers = numpy.where(numpy.sum(self._mask, axis=1) != self._data.shape[1])[0] + self._mask = hdu[extension_mask].data.astype("int32") + self._good_fibers = numpy.where(numpy.any(self._mask!=0, axis=1))[0] if extension_error is not None: self._error = hdu[extension_error].data.astype("float32") if extension_coeffs is not None: @@ -886,7 +901,7 @@ def measureArcLines( ) for i in iterator: spec = self.getSpec(i) - if spec._mask.all(): + if (spec._mask!=0).all(): masked[i] = True continue @@ -926,7 +941,7 @@ def measureArcLines( ) for i in iterator: spec = self.getSpec(i) - if spec._mask.all(): + if (spec._mask!=0).all(): masked[i] = True continue @@ -1025,7 +1040,7 @@ def fit_spline(self, degree=3, nknots=5, knots=None, smoothing=None, weights=Non poly_table = [] poly_all_table = [] for i in range(self._fibers): - good_pix = numpy.logical_not(self._mask[i, :]) + good_pix = self._mask[i, :] == 0 if numpy.sum(good_pix) >= nknots + 1: pixels_, data_ = pixels[good_pix], self._data[i, good_pix] (t0, c0, k) = _guess_spline(pixels_, data_, k=degree, s=smoothing, w=weights) @@ -1037,8 +1052,8 @@ def fit_spline(self, degree=3, nknots=5, knots=None, smoothing=None, weights=Non poly_table.extend(numpy.column_stack([pixels_, interpolate.splev(pixels_, tck)]).tolist()) poly_all_table.extend(numpy.column_stack([pixels, interpolate.splev(pixels, tck)]).tolist()) except ValueError as e: - log.error(f'Fiber trace failure at fiber {i}: {e}') - self._mask[i, :] = True + log.error(f'Fiber spline fitting failure at fiber {i}: {e}') + self._mask[i, :] = PixMask["FAILEDSPLINE"] continue self._coeffs[i] = tck @@ -1046,9 +1061,9 @@ def fit_spline(self, degree=3, nknots=5, knots=None, smoothing=None, weights=Non if clip is not None: self._data = numpy.clip(self._data, clip[0], clip[1]) - self._mask[i, :] = False + self._mask[i, :] = 0 else: - self._mask[i, :] = True + self._mask[i, :] = PixMask["FAILEDSPLINE"] return numpy.asarray(pix_table), numpy.asarray(poly_table), numpy.asarray(poly_all_table) @@ -1074,7 +1089,7 @@ def fit_polynomial(self, deg, poly_kind="poly", clip=None): poly_table = [] poly_all_table = [] for i in range(self._fibers): - good_pix = numpy.logical_not(self._mask[i, :]) + good_pix = self._mask[i, :] == 0 if numpy.sum(good_pix) >= deg + 1: # select the polynomial class poly_cls = Spectrum1D.select_poly_class(poly_kind) @@ -1086,8 +1101,8 @@ def fit_polynomial(self, deg, poly_kind="poly", clip=None): poly_table.extend(numpy.column_stack([pixels[good_pix], poly(pixels[good_pix])]).tolist()) poly_all_table.extend(numpy.column_stack([pixels, poly(pixels)]).tolist()) except numpy.linalg.LinAlgError as e: - log.error(f'Fiber trace failure at fiber {i}: {e}') - self._mask[i, :] = True + log.error(f'Fiber polynomial fitting failure at fiber {i}: {e}') + self._mask[i, :] = PixMask["FAILEDPOLY"] continue self._coeffs[i, :] = poly.convert().coef @@ -1095,9 +1110,9 @@ def fit_polynomial(self, deg, poly_kind="poly", clip=None): if clip is not None: self._data = numpy.clip(self._data, clip[0], clip[1]) - self._mask[i, :] = False + self._mask[i, :] = 0 else: - self._mask[i, :] = True + self._mask[i, :] = PixMask["FAILEDPOLY"] return numpy.asarray(pix_table), numpy.asarray(poly_table), numpy.asarray(poly_all_table) @@ -1158,7 +1173,7 @@ def smoothTraceDist( init_dist / dist ) # compute the relative change in the fiber distance compared to the reference dispersion column change_dist[:, i] = change # store the changes into array - good_mask = numpy.logical_not(bad_mask) + good_mask = bad_mask == 0 select_good = numpy.logical_and( numpy.logical_and(change > 0.5, change < 2.0), good_mask ) # masked unrealstic changes @@ -1187,7 +1202,7 @@ def smoothTraceDist( init_dist / dist ) # compute the relative change in the fiber distance compared to the reference dispersion column change_dist[:, i] = change # store the changes into array - good_mask = numpy.logical_not(bad_mask) + good_mask = bad_mask == 0 select_good = numpy.logical_and( numpy.logical_and(change > 0.5, change < 2.0), good_mask ) # masked unrealstic changes @@ -1323,7 +1338,7 @@ def interpolate_coeffs(self): if self._coeffs is None: return self # early return if all fibers are masked - bad_fibers = self._mask.all(axis=1) + bad_fibers = (self._mask!=0).all(axis=1) if bad_fibers.sum() == self._fibers: return self @@ -1339,7 +1354,7 @@ def interpolate_coeffs(self): for ifiber in y_pixels[bad_fibers]: poly = numpy.polynomial.Polynomial(self._coeffs[ifiber, :]) self._data[ifiber, :] = poly(x_pixels) - self._mask[ifiber, :] = False + self._mask[ifiber, :] = 0 return self @@ -1372,7 +1387,7 @@ def interpolate_data(self, axis="Y", extrapolate=False, reset_mask=True): # interpolate data if axis == "Y" or axis == "y" or axis == 0: - bad_fibers = self._mask.all(axis=1) + bad_fibers = (self._mask!=0).all(axis=1) if bad_fibers.sum() == self._fibers: return self f_data = interpolate.interp1d(y_pixels[~bad_fibers], self._data[~bad_fibers, :], axis=0, bounds_error=False, fill_value="extrapolate") @@ -1382,14 +1397,14 @@ def interpolate_data(self, axis="Y", extrapolate=False, reset_mask=True): self._error = f_error(y_pixels) # unmask interpolated fibers - if self._mask is not None: - self._mask[bad_fibers, :] = False + if self._mask is not None and reset_mask: + self._mask[bad_fibers, :] = 0 elif axis == "X" or axis == "x" or axis == 1: for ifiber in y_pixels: - bad_pixels = (self._data[ifiber] <= 0) | (self._mask[ifiber, :]) + bad_pixels = (self._data[ifiber] <= 0) | (self._mask[ifiber, :] != 0) # skip fiber if all pixels are bad and set mask to True if bad_pixels.all(): - self._mask[ifiber] = True + self._mask[ifiber] = PixMask["FAILEDINTERP"] continue # skip fiber if no bad pixels are present, no need to interpolate if bad_pixels.sum() == 0: @@ -1400,7 +1415,7 @@ def interpolate_data(self, axis="Y", extrapolate=False, reset_mask=True): f_error = interpolate.interp1d(x_pixels[~bad_pixels], self._error[ifiber, ~bad_pixels], bounds_error=False, fill_value="extrapolate") self._error[ifiber, :] = f_error(x_pixels) if self._mask is not None and reset_mask: - self._mask[ifiber, bad_pixels] = False + self._mask[ifiber, bad_pixels] = 0 else: raise ValueError(f"axis {axis} not supported") diff --git a/python/lvmdrp/core/image.py b/python/lvmdrp/core/image.py index 4436ab8e..1d0b7eb0 100644 --- a/python/lvmdrp/core/image.py +++ b/python/lvmdrp/core/image.py @@ -2087,7 +2087,7 @@ def fitSpline(self, axis="y", degree=3, smoothing=0, use_weights=False, clip=Non pixels = numpy.arange(self._dim[0]) models = numpy.zeros(self._dim) for i in range(self._dim[1]): - good_pix = ~self._mask[:,i] if self._mask is not None else ~numpy.isnan(self._data[:,i]) + good_pix = self._mask[:,i] == 0 if self._mask is not None else ~numpy.isnan(self._data[:,i]) # skip column if all pixels are masked if good_pix.sum() == 0: diff --git a/python/lvmdrp/core/rss.py b/python/lvmdrp/core/rss.py index 4fca4138..d247c1a4 100644 --- a/python/lvmdrp/core/rss.py +++ b/python/lvmdrp/core/rss.py @@ -11,6 +11,7 @@ from astropy import units as u from lvmdrp import log +from lvmdrp.utils.bitmask import PixMask, _parse_bitmask from lvmdrp.core.constants import CONFIG_PATH from lvmdrp.core.apertures import Aperture from lvmdrp.core.cube import Cube @@ -142,7 +143,7 @@ def from_file(cls, in_rss): if hdu.name == "ERROR": error = hdu.data.astype("float32") if hdu.name == "BADPIX": - mask = hdu.data.astype("bool") + mask = hdu.data.astype("int32") if hdu.name == "WAVE_TRACE": wave_trace = hdu if hdu.name == "LSF_TRACE": @@ -413,8 +414,8 @@ def from_channels(cls, rss_b, rss_r, rss_z, use_weights=True): fluxes.append(f(new_wave).astype("float32")) f = interpolate.interp1d(rss._wave, rss._error, axis=1, bounds_error=False, fill_value=numpy.nan) errors.append(f(new_wave).astype("float32")) - f = interpolate.interp1d(rss._wave, rss._mask, axis=1, kind="nearest", bounds_error=False, fill_value=0) - masks.append(f(new_wave).astype("uint8")) + f = interpolate.interp1d(rss._wave, rss._mask, axis=1, kind="nearest", bounds_error=False, fill_value=PixMask["NODATA"]) + masks.append(f(new_wave).astype("int32")) f = interpolate.interp1d(rss._wave, rss._lsf, axis=1, bounds_error=False, fill_value=numpy.nan) lsfs.append(f(new_wave).astype("float32")) if rss._sky is not None: @@ -472,7 +473,7 @@ def from_channels(cls, rss_b, rss_r, rss_z, use_weights=True): new_data = bn.nansum(fluxes * weights, axis=0) new_lsf = bn.nansum(lsfs * weights, axis=0) new_error = numpy.sqrt(bn.nansum(vars, axis=0)) - new_mask = (numpy.nansum(masks, axis=0)>0) + new_mask = numpy.bitwise_or(masks, axis=0) if rss._sky is not None: new_sky = bn.nansum(skies * weights, axis=0) else: @@ -502,7 +503,7 @@ def from_channels(cls, rss_b, rss_r, rss_z, use_weights=True): new_data = bn.nanmean(fluxes, axis=0) new_lsf = bn.nanmean(lsfs, axis=0) new_error = numpy.sqrt(bn.nanmean(vars, axis=0)) - new_mask = numpy.nansum(masks, axis=0).astype("bool") + new_mask = numpy.bitwise_or(masks, axis=0) if skies.size != 0: new_sky = bn.nansum(skies, axis=0) else: @@ -582,7 +583,7 @@ def from_spectra1d( error=numpy.zeros((n_spectra, ref_spec._error.size)) if ref_spec._error is not None else None, - mask=numpy.zeros((n_spectra, ref_spec._mask.size), dtype=bool) + mask=numpy.zeros((n_spectra, ref_spec._mask.size), dtype=int) if ref_spec._mask is not None else None, sky=numpy.zeros((n_spectra, ref_spec._data.size)) @@ -748,7 +749,7 @@ def __mul__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - mask = numpy.logical_or(self._mask, other._mask) + mask = numpy.bitwise_or(self._mask, other._mask) else: mask = self._mask if data.dtype == numpy.float64: @@ -791,7 +792,7 @@ def __mul__(self, other): # combined mask of valid pixels if contained in both if self._mask is not None and other._mask is not None: - mask = numpy.logical_or(self._mask, other._mask[numpy.newaxis, :]) + mask = numpy.bitwise_or(self._mask, other._mask[numpy.newaxis, :]) else: mask = self._mask @@ -1126,10 +1127,11 @@ def match_lsf(self, target_fwhm=None, min_fwhm=0.1): return RSS.from_spectra1d(new_specs, header=self._header, slitmap=self._slitmap, good_fibers=self._good_fibers) - def maskFiber(self, fiber, replace_error=1e10): + def maskFiber(self, fiber, bitmask, replace_error=1e10): + bitmask = _parse_bitmask(bitmask) self._data[fiber, :] = 0 if self._mask is not None: - self._mask[fiber, :] = True + self._mask[fiber, :] = bitmask if self._error is not None: self._error[fiber, :] = replace_error @@ -1410,7 +1412,7 @@ def extendData(self, new_wave): else: new_error = None if self._mask is not None: - new_mask = numpy.full((self._mask.shape[0], new_wave.size), False, dtype=bool) + new_mask = numpy.full((self._mask.shape[0], new_wave.size), False, dtype=numpy.int32) new_mask[:, ipix:fpix+1] = self._mask else: new_mask = None @@ -1469,7 +1471,7 @@ def combineRSS(self, rss_in, method="mean", replace_error=1e10): dim = rss_in[0]._data.shape data = numpy.zeros((len(rss_in), dim[0], dim[1]), dtype=numpy.float32) if rss_in[0]._mask is not None: - mask = numpy.zeros((len(rss_in), dim[0], dim[1]), dtype="bool") + mask = numpy.zeros((len(rss_in), dim[0], dim[1]), dtype=numpy.int32) else: mask = None @@ -1499,20 +1501,19 @@ def combineRSS(self, rss_in, method="mean", replace_error=1e10): if method == "sum": if mask is not None: data[mask] = 0 - good_pix = bn.nansum(numpy.logical_not(mask), 0) - select_mean = good_pix > 0 - combined_data[select_mean] = bn.nansum(data, 0)[select_mean] - combined_mask = good_pix == 0 + good_pix = bn.nansum(mask == 0, 0) > 0 + combined_data[good_pix] = bn.nansum(data, 0)[good_pix] + combined_mask = numpy.bitwise_or(mask, 0) if error is not None: error[mask] = replace_error - combined_error[select_mean] = numpy.sqrt( - bn.nansum(error**2, 0)[select_mean] + combined_error[good_pix] = numpy.sqrt( + bn.nansum(error**2, 0)[good_pix] ) else: combined_error = None if sky is not None: sky[mask] = 0 - combined_sky[select_mean] = bn.nansum(sky, 0)[select_mean] + combined_sky[good_pix] = bn.nansum(sky, 0)[good_pix] else: combined_mask = None combined_data = bn.nansum(data, 0) / data.shape[0] @@ -1530,24 +1531,23 @@ def combineRSS(self, rss_in, method="mean", replace_error=1e10): elif method == "mean": if mask is not None: data[mask] = 0 - good_pix = bn.nansum(numpy.logical_not(mask), 0) - select_mean = good_pix > 0 - combined_data[select_mean] = ( - bn.nansum(data, 0)[select_mean] / good_pix[select_mean] + good_pix = bn.nansum(mask == 0, 0) > 0 + combined_data[good_pix] = ( + bn.nansum(data, 0)[good_pix] / good_pix[good_pix] ) - combined_mask = good_pix == 0 + combined_mask = numpy.bitwise_or(mask, 0) if error is not None: error[mask] = replace_error - combined_error[select_mean] = numpy.sqrt( - bn.nansum(error**2, 0)[select_mean] - / good_pix[select_mean] ** 2 + combined_error[good_pix] = numpy.sqrt( + bn.nansum(error**2, 0)[good_pix] + / good_pix[good_pix] ** 2 ) else: combined_error = None if sky is not None: sky[mask] = 0 - combined_sky[select_mean] = ( - bn.nansum(sky, 0)[select_mean] / good_pix[select_mean] + combined_sky[good_pix] = ( + bn.nansum(sky, 0)[good_pix] / good_pix[good_pix] ) else: combined_sky = None @@ -1567,15 +1567,13 @@ def combineRSS(self, rss_in, method="mean", replace_error=1e10): elif method == "weighted_mean" and error is not None: if mask is not None: - good_pix = bn.nansum(numpy.logical_not(mask), 0) - select_mean = good_pix > 0 - + good_pix = bn.nansum(mask == 0, 0) > 0 var = error**2 weights = numpy.divide(1, var, out=numpy.zeros_like(var), where=var != 0) weights /= bn.nansum(weights, 0) combined_data[good_pix] = bn.nansum(data[good_pix] * var[good_pix], 0) combined_error[good_pix] = numpy.sqrt(bn.nansum(var[good_pix], 0)) - combined_mask = ~good_pix + combined_mask = numpy.bitwise_or(mask, 0) combined_error[combined_mask] = replace_error if sky is not None: combined_sky[good_pix] = bn.nansum(sky[good_pix] * var[good_pix], 0) @@ -1595,9 +1593,9 @@ def combineRSS(self, rss_in, method="mean", replace_error=1e10): elif method == "median": if mask is not None: - good_pix = bn.nansum(numpy.logical_not(mask), 0) + good_pix = bn.nansum(mask == 0, 0) > 0 combined_data[good_pix] = bn.nanmedian(data[good_pix], 0) - combined_mask = ~good_pix + combined_mask = numpy.bitwise_or(mask, 0) if error is not None: combined_error[good_pix] = numpy.sqrt(bn.nanmedian(error[good_pix] ** 2, 0)) combined_error[combined_mask] = replace_error @@ -1669,7 +1667,7 @@ def createAperSpec(self, cent_x, cent_y, radius): def create1DSpec(self, method="mean"): if self._wave is not None and len(self._wave.shape) == 2: if self._mask is not None: - select = numpy.logical_not(self._mask) + select = self._mask == 0 else: select = numpy.ones(self._data.shape, dtype="bool") disp = self._wave[:, 1:] - self._wave[:, :-1] @@ -1698,7 +1696,7 @@ def create1DSpec(self, method="mean"): sky = self._sky[select].flatten()[idx] else: if self._mask is not None: - select = numpy.logical_not(self._mask) + select = self._mask == 0 else: select = numpy.ones(self._data.shape, dtype="bool") @@ -1734,7 +1732,7 @@ def create1DSpec(self, method="mean"): sky[i] = numpy.sum(self._sky[select[:, i], i]) if self._mask is not None: - bad = numpy.sum(self._mask, 0) + bad = numpy.sum(self._mask != 0, 0) mask = bad == self._fibers else: mask = None @@ -1764,7 +1762,7 @@ def selectSpec(self, min=0, max=0, method="median"): spec = self[i] if spec._mask is not None: - goodpix = numpy.logical_not(spec._mask) + goodpix = spec._mask == 0 else: goodpix = numpy.ones(spec._data.dim[0], dtype=numpy.float32) @@ -1857,7 +1855,7 @@ def rectify_wave(self, wave=None, wave_range=None, wave_disp=None, method="linea new_rss = RSS( data=numpy.zeros((rss._fibers, wave.size), dtype="float32"), error=numpy.zeros((rss._fibers, wave.size), dtype="float32") if rss._error is not None else None, - mask=numpy.zeros((rss._fibers, wave.size), dtype="bool"), + mask=numpy.zeros((rss._fibers, wave.size), dtype="int32"), sky=numpy.zeros((rss._fibers, wave.size), dtype="float32") if rss._sky is not None else None, sky_error=numpy.zeros((rss._fibers, wave.size), dtype="float32") if rss._sky_error is not None else None, sky_east=numpy.zeros((rss._fibers, wave.size), dtype="float32") if rss._sky_east is not None else None, @@ -1887,10 +1885,10 @@ def rectify_wave(self, wave=None, wave_range=None, wave_disp=None, method="linea if rss._error is not None: f = interpolate.interp1d(rss._wave[ifiber], rss._error[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan) new_rss._error[ifiber] = f(wave).astype("float32") - f = interpolate.interp1d(rss._wave[ifiber], rss._mask[ifiber], kind="nearest", bounds_error=False, fill_value=1) - new_rss._mask[ifiber] = f(wave).astype("bool") - new_rss._mask[ifiber] |= numpy.isnan(new_rss._data[ifiber])|(new_rss._data[ifiber]==0) - new_rss._mask[ifiber] |= ~numpy.isfinite(new_rss._error[ifiber]) + f = interpolate.interp1d(rss._wave[ifiber], rss._mask[ifiber], kind="nearest", bounds_error=False, fill_value=PixMask["NODATA"]) + new_rss._mask[ifiber] = f(wave).astype("int32") + new_rss._mask[ifiber] |= numpy.isnan(new_rss._data[ifiber]) * PixMask["NODATA"] + new_rss._mask[ifiber] |= ~numpy.isfinite(new_rss._error[ifiber]) * PixMask["BADPIX"] if rss._lsf is not None: f = interpolate.interp1d(rss._wave[ifiber], rss._lsf[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan) new_rss._lsf[ifiber] = f(wave).astype("float32") @@ -2006,7 +2004,7 @@ def to_native_wave(self, method="linear", interp_density=True, return_density=Fa new_rss = RSS( data=numpy.zeros((rss._fibers, wave.shape[1]), dtype="float32"), error=numpy.zeros((rss._fibers, wave.shape[1]), dtype="float32"), - mask=numpy.zeros((rss._fibers, wave.shape[1]), dtype="bool"), + mask=numpy.zeros((rss._fibers, wave.shape[1]), dtype="int32"), sky=numpy.zeros((rss._fibers, wave.shape[1]), dtype="float32") if rss._sky is not None else None, sky_error=numpy.zeros((rss._fibers, wave.shape[1]), dtype="float32") if rss._sky_error is not None else None, cent_trace=rss._cent_trace, @@ -2023,8 +2021,8 @@ def to_native_wave(self, method="linear", interp_density=True, return_density=Fa new_rss._data[ifiber] = f(wave[ifiber]).astype("float32") f = interpolate.interp1d(rss._wave, rss._error[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan) new_rss._error[ifiber] = f(wave[ifiber]).astype("float32") - f = interpolate.interp1d(rss._wave, rss._mask[ifiber], kind="nearest", bounds_error=False, fill_value=1) - new_rss._mask[ifiber] = f(wave[ifiber]).astype("bool") + f = interpolate.interp1d(rss._wave, rss._mask[ifiber], kind="nearest", bounds_error=False, fill_value=PixMask["NODATA"]) + new_rss._mask[ifiber] = f(wave[ifiber]).astype("int32") if rss._sky is not None: f = interpolate.interp1d(rss._wave, rss._sky[ifiber], kind=method, bounds_error=False, fill_value=numpy.nan) new_rss._sky[ifiber] = f(wave[ifiber]).astype("float32") @@ -2159,7 +2157,7 @@ def createCubeInterpolation( if self._error is not None: var = self._error**2 inv_var = numpy.zeros_like(var) - good_pix = numpy.logical_not(self._mask) + good_pix = self._mask == 0 inv_var[good_pix] = 1.0 / var[good_pix] else: inv_var = numpy.ones_like(self._data) @@ -2201,12 +2199,12 @@ def createCubeInterpolation( * good_pix[i, :][:, numpy.newaxis] ) mask[:, select_bad] = numpy.logical_or( - mask[:, select_bad], self._mask[i, :][:, numpy.newaxis] + mask[:, select_bad], self._mask[i, :][:, numpy.newaxis] != 0 ) mask2[:, select_bad] = numpy.logical_or( mask2[:, select_bad], numpy.logical_and( - self._mask[i, :][:, numpy.newaxis], + self._mask[i, :][:, numpy.newaxis] != 0, self._data[i, :][:, numpy.newaxis] == 0, ), ) @@ -2255,12 +2253,12 @@ def createCubeInterpolation( # fiber_mask = (self._mask[i, :])[:, numpy.newaxis, numpy.newaxis] mask[:, select] = numpy.logical_or( - mask[:, select], self._mask[i, :][:, numpy.newaxis] + mask[:, select], self._mask[i, :][:, numpy.newaxis] != 0 ) mask2[:, select] = numpy.logical_or( mask2[:, select], numpy.logical_and( - self._mask[i, :][:, numpy.newaxis], + self._mask[i, :][:, numpy.newaxis] != 0, self._data[i, :][:, numpy.newaxis] == 0, ), ) @@ -2320,7 +2318,7 @@ def createCubeInterpolation( if self._error is not None: var = self._error**2 inv_var = numpy.zeros_like(var) - good_pix = numpy.logical_not(self._mask) + good_pix = self._mask == 0 inv_var[good_pix] = 1.0 / var[good_pix] else: inv_var = numpy.ones_like(self._data) @@ -2348,12 +2346,12 @@ def createCubeInterpolation( ) cover_fraction[select] = area mask[:, select] = numpy.logical_or( - mask[:, select], self._mask[i, :][:, numpy.newaxis] + mask[:, select], self._mask[i, :][:, numpy.newaxis] != 0 ) mask2[:, select] = numpy.logical_or( mask2[:, select], numpy.logical_and( - self._mask[i, :][:, numpy.newaxis], + self._mask[i, :][:, numpy.newaxis] != 0, self._data[i, :][:, numpy.newaxis] == 0, ), ) @@ -2442,7 +2440,7 @@ def createCubeInterDAR( if self._error is not None: var = self._error**2 inv_var = numpy.zeros_like(var) - good_pix = numpy.logical_not(self._mask) + good_pix = self._mask == 0 inv_var[good_pix] = 1.0 / var[good_pix] else: inv_var = numpy.ones_like(self._data) @@ -2473,13 +2471,13 @@ def createCubeInterDAR( mask[i, select_bad] = numpy.logical_or( mask[i, select_bad], numpy.logical_and( - self._mask[j, i], weights_0[i, select_bad] > 0 + self._mask[j, i] != 0, weights_0[i, select_bad] > 0 ), ) mask2[i, select_bad] = numpy.logical_or( mask2[i, select_bad], numpy.logical_and( - numpy.logical_and(self._mask[j, i], self._data[j, i] == 0), + numpy.logical_and(self._mask[j, i] != 0, self._data[j, i] == 0), weights_0[i, select_bad] > 0, ), ) @@ -2545,7 +2543,7 @@ def createCubeInterDAR( mask2, numpy.logical_and( cover_fraction > 0, - numpy.logical_and(self._mask[j, :], (self._data[j, :] == 0))[ + numpy.logical_and(self._mask[j, :] != 0, (self._data[j, :] == 0))[ :, numpy.newaxis, numpy.newaxis ], ), @@ -2811,7 +2809,7 @@ def createCubeInterDAR_new( mask, numpy.logical_and( cover_fraction > 0, - self._mask[j, :][:, numpy.newaxis, numpy.newaxis], + self._mask[j, :][:, numpy.newaxis, numpy.newaxis] != 0, ), ) mask2 = numpy.logical_or( @@ -2819,7 +2817,7 @@ def createCubeInterDAR_new( numpy.logical_and( cover_fraction > 0, numpy.logical_and( - self._mask[j, :], (self._data[j, :] == 0) + self._mask[j, :] != 0, (self._data[j, :] == 0) )[:, numpy.newaxis, numpy.newaxis], ), ) @@ -2864,7 +2862,7 @@ def createCubeInterDAR_new( resolution = float(resolution) # points = self._res_elements # fibers = self._fibers - good_pix = numpy.logical_not(self._mask) + good_pix = self._mask == 0 # arc_position_x = self._arc_position_x.astype(numpy.float32) # arc_position_y = self._arc_position_y.astype(numpy.float32) # size_x = self._size[0] @@ -3325,7 +3323,7 @@ def setSlitmap(self, slitmap): def apply_pixelmask(self, mask=None): """Replaces masked pixels in RSS by NaN values""" if mask is None: - mask = self._mask + mask = self._mask != 0 if mask is None: return self._data, self._error @@ -3459,7 +3457,7 @@ def writeFitsData(self, out_rss, replace_masked=True, include_wave=False): if self._error is not None: hdus.append(pyfits.ImageHDU(self._error.astype("float32"), name="ERROR")) if self._mask is not None: - hdus.append(pyfits.ImageHDU(self._mask.astype("uint8"), name="BADPIX")) + hdus.append(pyfits.ImageHDU(self._mask.astype("int32"), name="BADPIX")) # include wavelength extension for rectified RSSs if include_wave and self._wave and len(self._wave.shape) == 1: @@ -3591,7 +3589,7 @@ def from_hdulist(cls, hdulist): data = hdulist["FLUX"].data error = numpy.divide(1, hdulist["IVAR"].data, where=hdulist["IVAR"].data != 0, out=numpy.zeros_like(hdulist["IVAR"].data)) error = numpy.sqrt(error) - mask = hdulist["MASK"].data.astype("bool") + mask = hdulist["MASK"].data.astype("int32") cent_trace = Table(hdulist["CENT_TRACE"].data) width_trace = Table(hdulist["WIDTH_TRACE"].data) wave_trace = Table(hdulist["WAVE_TRACE"].data) @@ -3640,7 +3638,7 @@ def writeFitsData(self, out_file, replace_masked=True): # fill in rest of the template self._template["FLUX"].data = self._data self._template["IVAR"].data = numpy.divide(1, self._error**2, where=self._error != 0, out=numpy.zeros_like(self._error)) - self._template["MASK"].data = self._mask.astype("uint8") + self._template["MASK"].data = self._mask.astype("int32") self._template["WAVE_TRACE"] = pyfits.BinTableHDU(data=self._wave_trace, name="WAVE_TRACE") self._template["LSF_TRACE"] = pyfits.BinTableHDU(data=self._lsf_trace, name="LSF_TRACE") self._template["CENT_TRACE"] = pyfits.BinTableHDU(data=self._cent_trace, name="CENT_TRACE") @@ -3665,7 +3663,7 @@ def from_hdulist(cls, hdulist): data = hdulist["FLUX"].data error = numpy.divide(1, hdulist["IVAR"].data, where=hdulist["IVAR"].data != 0, out=numpy.zeros_like(hdulist["IVAR"].data)) error = numpy.sqrt(error) - mask = hdulist["MASK"].data.astype("bool") + mask = hdulist["MASK"].data.astype("int32") wave = hdulist["WAVE"].data lsf = hdulist["LSF"].data sky_east = hdulist["SKY_EAST"].data @@ -3715,7 +3713,7 @@ def writeFitsData(self, out_file, replace_masked=True): # fill in rest of the template self._template["FLUX"].data = self._data self._template["IVAR"].data = numpy.divide(1, self._error**2, where=self._error != 0, out=numpy.zeros_like(self._error)) - self._template["MASK"].data = self._mask.astype("uint8") + self._template["MASK"].data = self._mask.astype("int32") self._template["WAVE"].data = self._wave self._template["LSF"].data = self._lsf self._template["SKY_EAST"].data = self._sky_east @@ -3743,7 +3741,7 @@ def from_hdulist(cls, hdulist): data = hdulist["FLUX"].data error = numpy.divide(1, hdulist["IVAR"].data, where=hdulist["IVAR"].data != 0, out=numpy.zeros_like(hdulist["IVAR"].data)) error = numpy.sqrt(error) - mask = hdulist["MASK"].data.astype("bool") + mask = hdulist["MASK"].data.astype("int32") wave = hdulist["WAVE"].data lsf = hdulist["LSF"].data sky_east = hdulist["SKY_EAST"].data @@ -3789,7 +3787,7 @@ def writeFitsData(self, out_file, replace_masked=True): # fill in rest of the template self._template["FLUX"].data = self._data self._template["IVAR"].data = numpy.divide(1, self._error**2, where=self._error != 0, out=numpy.zeros_like(self._error)) - self._template["MASK"].data = self._mask.astype("uint8") + self._template["MASK"].data = self._mask.astype("int32") self._template["WAVE"].data = self._wave self._template["LSF"].data = self._lsf self._template["SKY_EAST"].data = self._sky_east @@ -3815,7 +3813,7 @@ def from_hdulist(cls, hdulist): data = hdulist["FLUX"].data error = numpy.divide(1, hdulist["IVAR"].data, where=hdulist["IVAR"].data != 0, out=numpy.zeros_like(hdulist["IVAR"].data)) error = numpy.sqrt(error) - mask = hdulist["MASK"].data.astype("bool") + mask = hdulist["MASK"].data.astype("int32") wave = hdulist["WAVE"].data lsf = hdulist["LSF"].data sky = hdulist["SKY"].data @@ -3852,7 +3850,7 @@ def writeFitsData(self, out_file, replace_masked=True): # fill in rest of the template self._template["FLUX"].data = self._data self._template["IVAR"].data = numpy.divide(1, self._error**2, where=self._error != 0, out=numpy.zeros_like(self._error)) - self._template["MASK"].data = self._mask.astype("uint8") + self._template["MASK"].data = self._mask.astype("int32") self._template["WAVE"].data = self._wave self._template["LSF"].data = self._lsf self._template["SKY"].data = self._sky.astype("float32") diff --git a/python/lvmdrp/core/spectrum1d.py b/python/lvmdrp/core/spectrum1d.py index f0c60742..47fac3d5 100644 --- a/python/lvmdrp/core/spectrum1d.py +++ b/python/lvmdrp/core/spectrum1d.py @@ -10,6 +10,7 @@ from typing import List, Tuple from lvmdrp import log +from lvmdrp.utils.bitmask import PixMask from lvmdrp.utils import gaussian from lvmdrp.core import fit_profile from lvmdrp.core import plot @@ -593,14 +594,14 @@ def __sub__(self, other): data = self._data - other._data if self._mask is not None and other._mask is not None: - mask = numpy.logical_or(self._mask, other._mask) - select_zero = numpy.logical_and(select_zero, mask) + mask = numpy.bitwise_or(self._mask, other._mask) + select_zero = numpy.logical_and(select_zero, mask!=0) data[select_zero] = 0 elif self._mask is None and other._mask is not None: mask = other._mask elif self._mask is not None and other._mask is None: mask = self._mask - select_zero = numpy.logical_and(select_zero, mask) + select_zero = numpy.logical_and(select_zero, mask!=0) data[select_zero] = 0 else: mask = None @@ -724,14 +725,14 @@ def __add__(self, other): data = self._data + other._data if self._mask is not None and other._mask is not None: - mask = numpy.logical_or(self._mask, other._mask) - select_zero = numpy.logical_and(select_zero, mask) + mask = numpy.bitwise_or(self._mask, other._mask) + select_zero = numpy.logical_and(select_zero, mask!=0) data[select_zero] = 0 elif self._mask is None and other._mask is not None: mask = other._mask elif self._mask is not None and other._mask is None: mask = self._mask - select_zero = numpy.logical_and(select_zero, mask) + select_zero = numpy.logical_and(select_zero, mask!=0) data[select_zero] = 0 else: mask = None @@ -888,14 +889,14 @@ def __truediv__(self, other): ) if self._mask is not None and other._mask is not None: - mask = numpy.logical_or(self._mask, other._mask) - mask[~select] = True + mask = numpy.bitwise_or(self._mask, other._mask) + mask[~select] = PixMask["BADPIX"] elif other._mask is not None: mask = other._mask - mask[~select] = True + mask[~select] = PixMask["BADPIX"] elif self._mask is not None: mask = self._mask - mask[~select] = True + mask[~select] = PixMask["BADPIX"] else: mask = None @@ -1011,7 +1012,7 @@ def __truediv__(self, other): if self._mask is not None: mask = self._mask - mask[~select] = True + mask[~select] = PixMask["BADPIX"] else: mask = None @@ -1069,7 +1070,7 @@ def __truediv__(self, other): if self._mask is not None: mask = self._mask - mask[~select] = True + mask[~select] = PixMask["BADPIX"] if self._sky is not None: sky = numpy.divide( @@ -1130,14 +1131,14 @@ def __rtruediv__(self, other): ) if self._mask is not None and other._mask is not None: - mask = numpy.logical_or(self._mask, other._mask) - mask[~select] = True + mask = numpy.bitwise_or(self._mask, other._mask) + mask[~select] = PixMask["BADPIX"] elif other._mask is not None: mask = other._mask - mask[~select] = True + mask[~select] = PixMask["BADPIX"] elif self._mask is not None: mask = self._mask - mask[~select] = True + mask[~select] = PixMask["BADPIX"] else: mask = None @@ -1256,7 +1257,7 @@ def __rtruediv__(self, other): if self._mask is not None: mask = self._mask - mask[~select] = True + mask[~select] = PixMask["BADPIX"] else: mask = None @@ -1322,7 +1323,7 @@ def __rtruediv__(self, other): if self._mask is not None: mask = self._mask - mask[~select] = True + mask[~select] = PixMask["BADPIX"] else: mask = None @@ -1783,7 +1784,7 @@ def loadFitsData( if hdu[i].header["EXTNAME"].split()[0] == "ERROR": self._error = hdu[i].data elif hdu[i].header["EXTNAME"].split()[0] == "BADPIX": - self._mask = hdu[i].data.astype("bool") + self._mask = hdu[i].data.astype("int32") elif hdu[i].header["EXTNAME"].split()[0] == "WAVE": self._wave = hdu[i].data elif hdu[i].header["EXTNAME"].split()[0] == "INSTFWHM": @@ -1801,7 +1802,7 @@ def loadFitsData( self._data = hdu[extension_data].data if extension_mask is not None: - self._mask = hdu[extension_mask].data + self._mask = hdu[extension_mask].data.astype("int32") if extension_error is not None: self._error = hdu[extension_error].data if extension_wave is not None: @@ -1891,7 +1892,7 @@ def writeFitsData( if self._error is not None: hdus[3] = pyfits.ImageHDU(self._error, name="ERROR") if self._mask is not None: - hdus[4] = pyfits.ImageHDU(self._mask.astype("uint8"), name="BADPIX") + hdus[4] = pyfits.ImageHDU(self._mask.astype("int32"), name="BADPIX") if self._sky is not None: hdus[5] = pyfits.ImageHDU(self._sky, name="SKY") if self._sky_error is not None: @@ -1916,10 +1917,10 @@ def writeFitsData( # mask hdu if extension_mask == 0: - hdu = pyfits.PrimaryHDU(self._mask.astype("uint8")) + hdu = pyfits.PrimaryHDU(self._mask.astype("int32")) elif extension_mask > 0 and extension_mask is not None: hdus[extension_mask] = pyfits.ImageHDU( - self._mask.astype("uint8"), name="BADPIX" + self._mask.astype("int32"), name="BADPIX" ) # error hdu @@ -2074,7 +2075,7 @@ def getData(self): Array of the data value error : numpy.ndarray (float) Array of the corresponding errors - mask : numpy.ndarray (bool) + mask : numpy.ndarray (int) Array of the bad pixel mask sky : numpy.ndarray (float) Array of the sky spectrum @@ -2105,10 +2106,10 @@ def resampleSpec( # case where input spectrum has more than half the pixels masked if numpy.nansum(self._data) == 0.0 or ( - self._mask is not None and numpy.sum(self._mask) > self._dim / 2 + self._mask is not None and numpy.sum(self._mask!=0) > self._dim / 2 ): # all pixels masked - new_mask = numpy.ones(len(ref_wave), dtype=bool) + new_mask = numpy.ones(len(ref_wave), dtype=numpy.int32) * PixMask["FAILEDINTERP"] # all data points to zero new_data = numpy.zeros(len(ref_wave), numpy.float32) # all LSF pixels zero (if present) @@ -2147,8 +2148,8 @@ def resampleSpec( else: # good pixels selection if self._mask is not None: - select_badpix = self._mask - select_goodpix = numpy.logical_not(self._mask) + select_badpix = self._mask != 0 + select_goodpix = self._mask == 0 else: select_badpix = numpy.zeros(self._dim, dtype=bool) select_goodpix = numpy.ones(self._dim, dtype=bool) @@ -2266,11 +2267,11 @@ def resampleSpec( # interpolate mask -------------------------------------------------------------------------------------------------------------------------------- if self._mask is not None: - badpix = numpy.zeros(ref_wave.shape[0], dtype=bool) + badpix = numpy.zeros(ref_wave.shape[0], dtype=numpy.float32) indices = numpy.arange(self._wave.shape[0]) - nbadpix = numpy.sum(self._mask) + nbadpix = numpy.sum(self._mask!=0) if nbadpix > 0: - badpix_id = indices[self._mask] + badpix_id = indices[self._mask!=0] for i in range(len(badpix_id)): badpix_min = badpix_id[i] - 2 badpix_max = badpix_id[i] + 2 @@ -2281,17 +2282,17 @@ def resampleSpec( ref_wave >= self._wave[bound[0]], ref_wave <= self._wave[bound[1]], ) - badpix = numpy.logical_or(badpix, select_bad) - new_mask = numpy.logical_or( + badpix = numpy.bitwise_or(badpix, select_bad * PixMask["FAILEDINTERP"]) + new_mask = numpy.bitwise_or( badpix, numpy.logical_or( ref_wave < self._wave[0], ref_wave > self._wave[-1] - ), + ) * PixMask["FAILEDINTERP"], ) else: new_mask = numpy.logical_or( ref_wave < self._wave[0], ref_wave > self._wave[-1] - ) + ) * PixMask["FAILEDINTERP"] # replace error values in masked pixels if new_error is not None: new_error[new_mask] = replace_error @@ -2442,7 +2443,9 @@ def resampleSpec_flux_conserving( def apply_pixelmask(self, mask=None, inplace=False): if mask is None: - mask = self._mask + mask = self._mask != 0 + else: + mask = mask != 0 if mask is None: return self @@ -2463,7 +2466,7 @@ def apply_pixelmask(self, mask=None, inplace=False): def interpolate_masked(self, mask=None, inplace=False): mask = mask if mask is not None else self._mask - if mask is None or mask.all(): + if mask is None or (mask!=0).all(): return self if inplace: @@ -2471,7 +2474,7 @@ def interpolate_masked(self, mask=None, inplace=False): else: new_spec = deepcopy(self) - good_pix = ~mask + good_pix = mask==0 new_spec._data = numpy.interp(new_spec._wave, new_spec._wave[good_pix], new_spec._data[good_pix], left=new_spec._data[good_pix][0], right=new_spec._data[good_pix][-1]) if new_spec._error is not None: new_spec._error = numpy.interp(new_spec._wave, new_spec._wave[good_pix], new_spec._error[good_pix], left=new_spec._error[good_pix][0], right=new_spec._error[good_pix][-1]) @@ -2582,12 +2585,12 @@ def binSpec(self, new_wave): new_disp = new_wave[1:] - new_wave[:-1] new_disp = numpy.insert(new_disp, 0, new_disp[0]) data_out = numpy.zeros(len(new_wave), dtype=numpy.float32) - mask_out = numpy.zeros(len(new_wave), dtype="bool") + mask_out = numpy.zeros(len(new_wave), dtype=numpy.int32) sky_out = numpy.zeros(len(new_wave), dtype=numpy.float32) sky_error_out = numpy.zeros(len(new_wave), dtype=numpy.float32) if self._mask is not None: - mask_in = numpy.logical_and(self._mask) + mask_in = self._mask == 0 else: mask_in = numpy.ones(len(self._wave), dtype="bool") # masked_data = self._wave[mask_in] @@ -2630,7 +2633,7 @@ def binSpec(self, new_wave): / numpy.sum(select) ** 2 ) else: - mask_out[i] = True + mask_out[i] = PixMask["FAILEDINTERP"] data_out = numpy.interp(new_wave, masked_wave, self._data[mask_in]) if self._sky is not None: sky_out = numpy.interp(new_wave, masked_wave, self._sky[mask_in]) @@ -2662,9 +2665,9 @@ def smoothSpec(self, size, method="gauss", mode="nearest"): self._data = ndimage.filters.median_filter(self._data, size, mode=mode) elif method == "BSpline": smooth = interpolate.splrep( - self._wave[~self._mask], - self._data[~self._mask], - w=1.0 / numpy.sqrt(numpy.fabs(self._data[~self._mask])), + self._wave[self._mask==0], + self._data[self._mask==0], + w=1.0 / numpy.sqrt(numpy.fabs(self._data[self._mask==0])), s=size, ) self._data = interpolate.splev(self._wave, smooth, der=0) @@ -2672,7 +2675,7 @@ def smoothSpec(self, size, method="gauss", mode="nearest"): def smoothGaussVariable(self, diff_fwhm): fact = numpy.sqrt(2.0 * numpy.pi) if self._mask is not None: - mask = numpy.logical_not(self._mask) + mask = self._mask == 0 else: mask = numpy.ones(self._data.shape[0], dtype="bool") @@ -2747,7 +2750,7 @@ def smoothPoly( self, deg=5, poly_kind="legendre", start_wave=None, end_wave=None, ref_base=None ): if self._mask is not None: - mask = numpy.logical_not(self._mask) + mask = self._mask == 0 else: mask = numpy.ones(self._dim, dtype="bool") if start_wave is not None: @@ -2778,7 +2781,7 @@ def smoothPoly( else: self._data[:] = 0 if self._mask is not None: - self._mask[:] = True + self._mask[:] = PixMask["FAILEDPOLY"] out_par = 0 return out_par @@ -2830,7 +2833,7 @@ def findPeaks( wave = self._wave[s] pixels = self._pixels[s] if self._mask is not None: - mask = self._mask[s] + mask = self._mask[s] != 0 else: mask = numpy.zeros_like(data, dtype=bool) @@ -3345,7 +3348,7 @@ def fitSepGauss( warning=False, ): error = self._error if self._error is not None else numpy.ones(self._dim, dtype=numpy.float32) - mask = self._mask if self._mask is not None else numpy.zeros(self._dim, dtype=bool) + mask = self._mask != 0 if self._mask is not None else numpy.zeros(self._dim, dtype=bool) error[mask] = numpy.inf data = self._data.copy() data[mask] = 0.0 @@ -3398,7 +3401,7 @@ def fitSepGauss( fwhm[i] *= 2.354 if axs is not None: - axs[i] = gauss.plot(self._wave[select], self._data[select], mask=self._mask[select], ax=axs[i]) + axs[i] = gauss.plot(self._wave[select], self._data[select], mask=self._mask[select]!=0, ax=axs[i]) axs[i].axvline(cent_guess[i], ls="--", lw=1, color="tab:red", label="cent. guess") axs[i].set_title(f"{axs[i].get_title()} @ {cent[i]:.1f} (pixel)") axs[i].text(0.05, 0.9, f"flux = {flux[i]:.2f}", va="bottom", ha="left", transform=axs[i].transAxes, fontsize=11) @@ -3453,7 +3456,7 @@ def obtainGaussFluxPeaks(self, pos, sigma, indices, replace_error=1e10, plot=Fal bad_pix[select] = True # masking fibers if all pixels are bad within aperture - bad_pix[nselect] = bn.nansum(self._mask[pixels[nselect, :]], 1) == aperture + bad_pix[nselect] = bn.nansum(self._mask[pixels[nselect, :]] != 0, 1) == aperture else: bad_pix = None if self._error is None: @@ -3508,7 +3511,7 @@ def collapseSpec(self, method="mean", start=None, end=None, transmission_func=No select_end = numpy.ones(self._dim, dtype="bool") select = numpy.logical_and(select_start, select_end) if self._mask is not None: - select = numpy.logical_and(select, numpy.logical_not(self._mask)) + select = numpy.logical_and(select, self._mask == 0) if method != "mean" and method != "median" and method != "sum": raise ValueError("method must be either 'mean', 'median' or 'sum'") @@ -3552,7 +3555,7 @@ def coaddSpec(self, other, wave=None): fluxes = numpy.ma.zeros((2, len(wave))) errors = numpy.zeros_like(fluxes) fwhms = numpy.zeros_like(fluxes) - masks = numpy.zeros_like(fluxes, dtype=bool) + masks = numpy.zeros_like(fluxes, dtype=numpy.int32) skies = numpy.zeros_like(fluxes) sky_errors = numpy.zeros_like(fluxes) for i, s in enumerate(spectra): @@ -3594,11 +3597,11 @@ def coaddSpec(self, other, wave=None): skies = bn.nansum(skies * weights, axis=0) sky_errors = numpy.sqrt(bn.nansum(sky_errors**2 * weights**2), axis=0) - masks = numpy.logical_and(masks[0], masks[1]) - masks = numpy.logical_or(masks, numpy.isnan(fluxes)) - masks = numpy.logical_or(masks, numpy.isnan(fwhms)) - masks = numpy.logical_or(masks, numpy.isnan(errors)) - masks = numpy.logical_or(masks, numpy.isnan(skies)) - masks = numpy.logical_or(masks, numpy.isnan(sky_errors)) + masks = numpy.bitwise_and(masks[0], masks[1]) + masks = numpy.bitwise_or(masks, numpy.isnan(fluxes)) + masks = numpy.bitwise_or(masks, numpy.isnan(fwhms)) + masks = numpy.bitwise_or(masks, numpy.isnan(errors)) + masks = numpy.bitwise_or(masks, numpy.isnan(skies)) + masks = numpy.bitwise_or(masks, numpy.isnan(sky_errors)) return Spectrum1D(wave=wave, data=fluxes, error=errors, lsf=fwhms, mask=masks, sky=skies, sky_error=sky_errors) diff --git a/python/lvmdrp/utils/bitmask.py b/python/lvmdrp/utils/bitmask.py index ed9a51b4..01839ac5 100644 --- a/python/lvmdrp/utils/bitmask.py +++ b/python/lvmdrp/utils/bitmask.py @@ -226,6 +226,9 @@ class PixMask(BaseBitmask): INTERPOLATED = auto() # measure quality of tracing using polynomial fit - samples residuals + FAILEDPOLY = auto() + FAILEDSPLINE = auto() + FAILEDINTERP = auto() BADTRACE = auto() BADARC = auto() From e3675225fd2b30e7185c5549ca6321f99f0775cf Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Thu, 12 Sep 2024 14:25:58 -0300 Subject: [PATCH 09/20] tweaking bitmask classes & routines --- python/lvmdrp/core/image.py | 2 + python/lvmdrp/utils/bitmask.py | 85 ++++++++++++++-------------------- 2 files changed, 38 insertions(+), 49 deletions(-) diff --git a/python/lvmdrp/core/image.py b/python/lvmdrp/core/image.py index 0969c6e8..cd6ce261 100644 --- a/python/lvmdrp/core/image.py +++ b/python/lvmdrp/core/image.py @@ -1081,6 +1081,8 @@ def add_bitmask(self, pixmask, where=None): mask : np.ndarray[PixMask|int] Updated pixel mask """ + if self._mask is None and self._dim is not None: + self._mask = numpy.zeros(self._dim) self._mask = add_bitmask(self._mask, pixmask=pixmask, where=where) return self._mask diff --git a/python/lvmdrp/utils/bitmask.py b/python/lvmdrp/utils/bitmask.py index 01839ac5..edca5ee0 100644 --- a/python/lvmdrp/utils/bitmask.py +++ b/python/lvmdrp/utils/bitmask.py @@ -219,38 +219,7 @@ class QualityFlag(BaseBitmask): class PixMask(BaseBitmask): - # fiber bitmasks - NONEXPOSED = auto() - WEAKFIBER = auto() - DEADFIBER = auto() - INTERPOLATED = auto() - - # measure quality of tracing using polynomial fit - samples residuals - FAILEDPOLY = auto() - FAILEDSPLINE = auto() - FAILEDINTERP = auto() - BADTRACE = auto() - BADARC = auto() - - # measure offset of the fiber from the median flatfielded fiber - BADFLAT = auto() - - # offset from a preset fiber shift value - LARGESHIFT = auto() - - BADSTDFIBER = auto() - BADSKYFIBER = auto() - - # pixel bitmasks - - # pixels with no useful information - NODATA = auto() - - # TODO: bright pixels on top and bottom edges - - # set this if X% close to saturation level - SATURATION = auto() - + # pixel bitmasks ------------------------------------------------------ # from pixelmasks BADPIX = auto() NEARBADPIXEL = auto() @@ -280,6 +249,36 @@ class PixMask(BaseBitmask): # large sky residuals BADSKYCHI = auto() + # fiber bitmasks ------------------------------------------------------ + NONEXPOSED = auto() + WEAKFIBER = auto() + DEADFIBER = auto() + INTERPOLATED = auto() + + # measure quality of tracing using polynomial fit - samples residuals + FAILEDPOLY = auto() + FAILEDSPLINE = auto() + FAILEDINTERP = auto() + BADTRACE = auto() + BADARC = auto() + + # measure offset of the fiber from the median flatfielded fiber + BADFLAT = auto() + + # offset from a preset fiber shift value + LARGESHIFT = auto() + + BADSTDFIBER = auto() + BADSKYFIBER = auto() + + # pixels with no useful information + NODATA = auto() + + # TODO: bright pixels on top and bottom edges + + # set this if X% close to saturation level + SATURATION = auto() + # define flag name constants # RAW_QUALITIES = list(RawFrameQuality.__members__.keys()) @@ -304,8 +303,10 @@ def _parse_bitmask(pixmask, mask_shape): def _parse_where(where, mask_shape): - if where is not None and isinstance(where, np.ndarray[bool]): + if where is not None: + assert isinstance(where, np.ndarray), f"Wrong type for `where` {type(where)}, expected `numpy.ndarray`" assert where.shape == mask_shape, f"Wrong `where` shape {where.shape} not matching `mask_image` shape {mask_shape}" + assert isinstance(where[0,0], np.bool_), f"Wrong `where` Numpy dtype {type(where[0,0])}, expected a boolean array" else: where = np.ones(mask_shape, dtype=bool) @@ -316,14 +317,7 @@ def add_bitmask(mask_image, pixmask, where=None): pixmask = _parse_bitmask(pixmask, mask_shape=mask_image.shape) where = _parse_where(where, mask_shape=mask_image.shape) - mask = np.zeros_like(mask_image, dtype=int) - mask[where] |= pixmask - - if mask_image is None: - mask_image = mask - return mask_image - - mask_image[where] |= mask[where] + mask_image[where] |= pixmask return mask_image @@ -331,12 +325,5 @@ def toggle_bitmask(mask_image, pixmask, where=None): pixmask = _parse_bitmask(pixmask, mask_shape=mask_image.shape) where = _parse_where(where, mask_shape=mask_image.shape) - mask = np.zeros_like(mask_image, dtype=int) - mask[where] |= pixmask - - if mask_image is None: - mask_image = mask - return mask_image - - mask_image[where] ^= mask[where] + mask_image[where] ^= pixmask return mask_image From dd4afbe575312d1f5d8453ccdd68cbf3d722c858 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Thu, 12 Sep 2024 16:43:01 -0300 Subject: [PATCH 10/20] adding bitmask handling methods to Image and RSS classes --- python/lvmdrp/core/image.py | 19 ++++++++++- python/lvmdrp/core/rss.py | 63 ++++++++++++++++++++++++++++++++++++- 2 files changed, 80 insertions(+), 2 deletions(-) diff --git a/python/lvmdrp/core/image.py b/python/lvmdrp/core/image.py index cd6ce261..126fbe04 100644 --- a/python/lvmdrp/core/image.py +++ b/python/lvmdrp/core/image.py @@ -1082,7 +1082,7 @@ def add_bitmask(self, pixmask, where=None): Updated pixel mask """ if self._mask is None and self._dim is not None: - self._mask = numpy.zeros(self._dim) + self._mask = numpy.zeros(self._dim, dtype=int) self._mask = add_bitmask(self._mask, pixmask=pixmask, where=where) return self._mask @@ -1107,6 +1107,23 @@ def toggle_bitmask(self, pixmask, where=None): self._mask = toggle_bitmask(self._mask, pixmask=pixmask, where=where) return self._mask + def print_bitmasks(self, logger=None): + """Prints or logs bitmask value counts + + Parameters + ---------- + logger : logger, optional + logs with level `info` if given, by default None (uses print) + """ + if self._mask is None: + return + uniques, counts = numpy.unique(self._mask, return_counts=True) + bitmasks = dict(zip(map(lambda p: PixMask(p).name if p>0 else "GOODPIX", uniques), counts)) + if logger: + logger.info(f"{bitmasks}") + return + print(bitmasks) + def convertUnit(self, to, assume="adu", gain_field="GAIN", inplace=False): """converts the unit of the image diff --git a/python/lvmdrp/core/rss.py b/python/lvmdrp/core/rss.py index 8315bf1d..96101804 100644 --- a/python/lvmdrp/core/rss.py +++ b/python/lvmdrp/core/rss.py @@ -14,7 +14,7 @@ from astropy import units as u from lvmdrp import log -from lvmdrp.utils.bitmask import PixMask, _parse_bitmask +from lvmdrp.utils.bitmask import PixMask, _parse_bitmask, add_bitmask, toggle_bitmask from lvmdrp.core.constants import CONFIG_PATH from lvmdrp.core.apertures import Aperture from lvmdrp.core.cube import Cube @@ -957,6 +957,67 @@ def __setitem__(self, fiber, spec): if self._sky_error is not None and spec._sky_error is not None: self._sky_error[fiber] = spec._sky_error + def add_bitmask(self, pixmask, where=None): + """Adds the given bitmask to the current mask + + If no mask is present in the image, a new mask will be started + with the corresponding pixels flagged with `pixmask`. + + Parameters + ---------- + pixmask : PixMask|str|int + Bitmask to add to the current mask + where : numpy.ndarray[bool], optional + Boolean selection of pixels where to add given bitmask, by default all pixels + + Returns + ------- + mask : np.ndarray[PixMask|int] + Updated pixel mask + """ + if self._mask is None and self._dim is not None: + self._mask = numpy.zeros(self._dim, dtype=int) + self._mask = add_bitmask(self._mask, pixmask=pixmask, where=where) + return self._mask + + def toggle_bitmask(self, pixmask, where=None): + """Toggle the given bitmask in the current mask + + If no mask is present in the image, a new mask will be started + with the corresponding pixels flagged with `pixmask`. + + Parameters + ---------- + pixmask : PixMask|str|int + Bitmask to add to the current mask + where : numpy.ndarray[bool], optional + Boolean selection of pixels where to add given bitmask, by default all pixels + + Returns + ------- + mask : np.ndarray[PixMask|int] + Updated pixel mask + """ + self._mask = toggle_bitmask(self._mask, pixmask=pixmask, where=where) + return self._mask + + def print_bitmasks(self, logger=None): + """Prints or logs bitmask value counts + + Parameters + ---------- + logger : logger, optional + logs with level `info` if given, by default None (uses print) + """ + if self._mask is None: + return + uniques, counts = numpy.unique(self._mask, return_counts=True) + bitmasks = dict(zip(map(lambda p: PixMask(p).name if p>0 else "GOODPIX", uniques), counts)) + if logger: + logger.info(f"{bitmasks}") + return + print(bitmasks) + def add_header_comment(self, comstr): ''' Append a COMMENT card at the end of the FITS header. From 8030e82035bc491f4b12ac63fe96af9d24aa0e28 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Thu, 12 Sep 2024 20:30:43 -0300 Subject: [PATCH 11/20] consolidating bitmask print methods --- python/lvmdrp/core/image.py | 7 +------ python/lvmdrp/core/rss.py | 9 ++------- python/lvmdrp/utils/bitmask.py | 9 +++++++++ 3 files changed, 12 insertions(+), 13 deletions(-) diff --git a/python/lvmdrp/core/image.py b/python/lvmdrp/core/image.py index 37564d6d..f1c5f56b 100644 --- a/python/lvmdrp/core/image.py +++ b/python/lvmdrp/core/image.py @@ -1117,12 +1117,7 @@ def print_bitmasks(self, logger=None): """ if self._mask is None: return - uniques, counts = numpy.unique(self._mask, return_counts=True) - bitmasks = dict(zip(map(lambda p: PixMask(p).name if p>0 else "GOODPIX", uniques), counts)) - if logger: - logger.info(f"{bitmasks}") - return - print(bitmasks) + print_bitmasks(self._mask, logger=log) def convertUnit(self, to, assume="adu", gain_field="GAIN", inplace=False): """converts the unit of the image diff --git a/python/lvmdrp/core/rss.py b/python/lvmdrp/core/rss.py index 96101804..c4b0d47a 100644 --- a/python/lvmdrp/core/rss.py +++ b/python/lvmdrp/core/rss.py @@ -14,7 +14,7 @@ from astropy import units as u from lvmdrp import log -from lvmdrp.utils.bitmask import PixMask, _parse_bitmask, add_bitmask, toggle_bitmask +from lvmdrp.utils.bitmask import PixMask, _parse_bitmask, add_bitmask, toggle_bitmask, print_bitmasks from lvmdrp.core.constants import CONFIG_PATH from lvmdrp.core.apertures import Aperture from lvmdrp.core.cube import Cube @@ -1011,12 +1011,7 @@ def print_bitmasks(self, logger=None): """ if self._mask is None: return - uniques, counts = numpy.unique(self._mask, return_counts=True) - bitmasks = dict(zip(map(lambda p: PixMask(p).name if p>0 else "GOODPIX", uniques), counts)) - if logger: - logger.info(f"{bitmasks}") - return - print(bitmasks) + print_bitmasks(self._mask, logger=log) def add_header_comment(self, comstr): ''' diff --git a/python/lvmdrp/utils/bitmask.py b/python/lvmdrp/utils/bitmask.py index edca5ee0..d49068be 100644 --- a/python/lvmdrp/utils/bitmask.py +++ b/python/lvmdrp/utils/bitmask.py @@ -327,3 +327,12 @@ def toggle_bitmask(mask_image, pixmask, where=None): mask_image[where] ^= pixmask return mask_image + + +def print_bitmasks(mask_array, logger=None): + uniques, counts = np.unique(mask_array, return_counts=True) + bitmasks = dict(zip(map(lambda p: PixMask(int(p)).get_name() if p>0 else "GOODPIX", uniques), counts)) + if logger: + logger.info(f"{bitmasks}") + return + print(bitmasks) From 7979be20ccb38f56cc8b10649bf27d9981188092 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Thu, 12 Sep 2024 21:13:36 -0300 Subject: [PATCH 12/20] fixing apply_mask method in RSS class --- python/lvmdrp/core/rss.py | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/python/lvmdrp/core/rss.py b/python/lvmdrp/core/rss.py index c4b0d47a..5ab05a8a 100644 --- a/python/lvmdrp/core/rss.py +++ b/python/lvmdrp/core/rss.py @@ -3382,25 +3382,26 @@ def apply_pixelmask(self, mask=None): """Replaces masked pixels in RSS by NaN values""" if mask is None: mask = self._mask != 0 + else: + mask = mask != 0 if mask is None: return self._data, self._error - if self._mask is not None: - self._data[self._mask] = numpy.nan - if self._error is not None: - self._error[self._mask] = numpy.nan - if self._sky is not None: - self._sky[self._mask] = numpy.nan - if self._sky_error is not None: - self._sky_error[self._mask] = numpy.nan - if self._sky_east is not None: - self._sky_east[self._mask] = numpy.nan - if self._sky_east_error is not None: - self._sky_east_error[self._mask] = numpy.nan - if self._sky_west is not None: - self._sky_west[self._mask] = numpy.nan - if self._sky_west_error is not None: - self._sky_west_error[self._mask] = numpy.nan + self._data[mask] = numpy.nan + if self._error is not None: + self._error[mask] = numpy.nan + if self._sky is not None: + self._sky[mask] = numpy.nan + if self._sky_error is not None: + self._sky_error[mask] = numpy.nan + if self._sky_east is not None: + self._sky_east[mask] = numpy.nan + if self._sky_east_error is not None: + self._sky_east_error[mask] = numpy.nan + if self._sky_west is not None: + self._sky_west[mask] = numpy.nan + if self._sky_west_error is not None: + self._sky_west_error[mask] = numpy.nan return self._data, self._error From dcea1a6fa98a58e063386254931c69bc48c8df78 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Thu, 12 Sep 2024 21:54:11 -0300 Subject: [PATCH 13/20] fixing propagation of bitmasks during extraction --- python/lvmdrp/core/image.py | 39 ++++++++++---------------- python/lvmdrp/core/spectrum1d.py | 15 ++++------ python/lvmdrp/functions/imageMethod.py | 10 +++---- 3 files changed, 26 insertions(+), 38 deletions(-) diff --git a/python/lvmdrp/core/image.py b/python/lvmdrp/core/image.py index f1c5f56b..d8fb1628 100644 --- a/python/lvmdrp/core/image.py +++ b/python/lvmdrp/core/image.py @@ -17,7 +17,7 @@ from scipy import interpolate from lvmdrp import log -from lvmdrp.utils.bitmask import PixMask, add_bitmask, toggle_bitmask +from lvmdrp.utils.bitmask import PixMask, add_bitmask, toggle_bitmask, print_bitmasks from lvmdrp.core.constants import CON_LAMPS, ARC_LAMPS from lvmdrp.core.plot import plt from lvmdrp.core.fit_profile import gaussians, Gaussians @@ -2602,7 +2602,7 @@ def extractSpecAperture(self, TraceMask, aperture): error = numpy.zeros((TraceMask._fibers, self._dim[1]), dtype=numpy.float32) else: error = None - mask = numpy.zeros((TraceMask._fibers, self._dim[1]), dtype="bool") + mask = numpy.zeros((TraceMask._fibers, self._dim[1]), dtype=numpy.int32) for i in range(self._dim[1]): pixels = numpy.round( pos[:, i][:, numpy.newaxis] @@ -2630,10 +2630,10 @@ def extractSpecAperture(self, TraceMask, aperture): numpy.sum(self._error[:, i][pixels] ** 2, 1) ) if self._mask is not None: - mask[good_pix[:, i], i] = numpy.sum(self._mask[:, i][pixels], 1) > 0 + mask[good_pix[:, i], i] = numpy.bitwise_or.reduce(self._mask[:, i][pixels], 1) # update mask with trace mask - mask |= bad_pix + # mask |= TraceMask._mask.astype(int) return data, error, mask def extractSpecOptimal(self, cent_trace, trace_fwhm, plot_fig=False): @@ -2643,7 +2643,7 @@ def extractSpecOptimal(self, cent_trace, trace_fwhm, plot_fig=False): error = numpy.zeros((cent_trace._fibers, self._dim[1]), dtype=numpy.float32) else: error = None - mask = numpy.zeros((cent_trace._fibers, self._dim[1]), dtype="bool") + mask = numpy.zeros((cent_trace._fibers, self._dim[1]), dtype=numpy.int32) self._data = numpy.nan_to_num(self._data) self._error = numpy.nan_to_num(self._error, nan=numpy.inf) @@ -2654,19 +2654,10 @@ def extractSpecOptimal(self, cent_trace, trace_fwhm, plot_fig=False): for i in range(self._dim[1]): # get i-column from image and trace slice_img = self.getSlice(i, axis="y") - slice_cent = cent_trace.getSlice(i, axis="y") - cent = slice_cent[0] + cent, cent_err, cent_mask = cent_trace.getSlice(i, axis="y") # define fiber mask - bad_fiber = (slice_cent[2] == 1) | ( - (slice_cent[0] < 0) | (slice_cent[0] > len(slice_img._data) - 1) - ) - # bad_fiber = numpy.logical_or( - # (slice_cent[2] == 1), - # numpy.logical_or( - # slice_cent[0] < 0, slice_cent[0] > len(slice_img._data) - 1 - # ), - # ) + bad_fiber = (cent_mask != 0) | ((cent < 0) | (cent > len(slice_img._data) - 1)) good_fiber = ~bad_fiber # get i-column from sigma trace @@ -2677,13 +2668,13 @@ def extractSpecOptimal(self, cent_trace, trace_fwhm, plot_fig=False): slice_img._data[select_nan] = 0 # measure flux along the given columns - result = slice_img.obtainGaussFluxPeaks(cent[good_fiber], sigma[good_fiber], plot=plot_fig) - data[good_fiber, i] = result[0] + edata, eerror, emask = slice_img.obtainGaussFluxPeaks(cent[good_fiber], sigma[good_fiber], plot=plot_fig) + data[good_fiber, i] = edata if self._error is not None: - error[good_fiber, i] = result[1] + error[good_fiber, i] = eerror if self._mask is not None: - mask[good_fiber, i] = result[2] - mask[bad_fiber, i] = True + mask[good_fiber, i] |= emask + mask[bad_fiber, i] |= cent_mask[bad_fiber] return data, error, mask def maskFiberTraces(self, TraceMask, aperture=3, parallel="auto"): @@ -2691,7 +2682,7 @@ def maskFiberTraces(self, TraceMask, aperture=3, parallel="auto"): dx = numpy.arange(-n0, n0 + 1, 0.5) trace = TraceMask.getData()[0] if self._mask is None: - self._mask = numpy.zeros(self._dim, dtype="bool") + self._mask = numpy.zeros(self._dim, dtype="int32") if parallel == "auto": cpus = cpu_count() else: @@ -2722,7 +2713,7 @@ def maskFiberTraces(self, TraceMask, aperture=3, parallel="auto"): for i in range(self._dim[1]): if cpus > 1: - self._mask[mask_pixels[i].get(), i] = True + self._mask[mask_pixels[i].get(), i] |= PixMask["NODATA"] else: select = numpy.unique( numpy.clip( @@ -2735,7 +2726,7 @@ def maskFiberTraces(self, TraceMask, aperture=3, parallel="auto"): .astype("int16") .flatten() ) - self._mask[select, i] = True + self._mask[select, i] |= PixMask["NODATA"] def peakPosition(self, guess_x=None, guess_y=None, box_x=None, box_y=None): image = self._data * numpy.logical_not(self._mask) diff --git a/python/lvmdrp/core/spectrum1d.py b/python/lvmdrp/core/spectrum1d.py index e66af5ce..0a39fc4c 100644 --- a/python/lvmdrp/core/spectrum1d.py +++ b/python/lvmdrp/core/spectrum1d.py @@ -3440,21 +3440,20 @@ def obtainGaussFluxPeaks(self, pos, sigma, replace_error=1e10, plot=False): """ nfibers = len(pos) - aperture = 3 + nsigmas = 1 # round up fiber locations pixels = numpy.round( - pos[:, None] + numpy.arange(-aperture / 2.0, aperture / 2.0, 1.0)[None, :] + pos[:, None] + numpy.arange(-nsigmas, nsigmas+1, 1)[None, :] * bn.nanmedian(sigma) ).astype("int") # defining bad pixels for each fiber if needed if self._mask is not None: # select: fibers in the boundary of the chip - bad_pix = numpy.zeros(nfibers, dtype="bool") - select = bn.nansum(pixels >= self._mask.shape[0], 1) - nselect = numpy.logical_not(select) - bad_pix[select] = True + bad_pix = numpy.zeros(nfibers, dtype=numpy.int32) + select = numpy.logical_or.reduce(pixels >= self._mask.shape[0], 1) + bad_pix[select] |= numpy.bitwise_or.reduce(self._mask[pixels[select, :]], 1) # masking fibers if all pixels are bad within aperture - bad_pix[nselect] = bn.nansum(self._mask[pixels[nselect, :]] != 0, 1) == aperture + bad_pix[~select] |= numpy.bitwise_or.reduce(self._mask[pixels[~select, :]], 1) else: bad_pix = None if self._error is None: @@ -3483,8 +3482,6 @@ def obtainGaussFluxPeaks(self, pos, sigma, replace_error=1e10, plot=False): error = numpy.sqrt(1 / ((B.multiply(B)).sum(axis=0))).A error = error[0,:] - if bad_pix is not None and bn.nansum(bad_pix) > 0: - error[bad_pix] = replace_error # pyfits.writeto('B1.fits', B1.toarray(), overwrite=True) # if plot: diff --git a/python/lvmdrp/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index 989dc502..d13f01ea 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -2743,13 +2743,15 @@ def extract_spectra( select_spec = slitmap["spectrographid"] == int(img._header["CCD"][1]) slitmap_spec = slitmap[select_spec] exposed_selection = numpy.array(list(img._header["STD*ACQ"].values())) - # mask fibers that are not exposed # TODO: use the more reliable routine get_exposed_std_fibers once is merged from addqa branch if len(exposed_selection) != 0: exposed_std = numpy.array(list(img._header["STD*FIB"].values()))[exposed_selection] - mask |= (~(numpy.isin(slitmap_spec["orig_ifulabel"], exposed_std))&((slitmap_spec["telescope"] == "Spec")))[:, None] - mask |= (slitmap_spec["fibstatus"] == 1)[:, None] + mask |= (~(numpy.isin(slitmap_spec["orig_ifulabel"], exposed_std))&((slitmap_spec["telescope"] == "Spec")))[:, None] * PixMask["NONEXPOSED"] + # mask dead fibers + mask |= (slitmap_spec["fibstatus"] == 1)[:, None] * PixMask["DEADFIBER"] + # print(numpy.unique(mask)) + # exit() drpstage += "SPECTRA_EXTRACTED" # propagate thermal shift to slitmap @@ -2757,8 +2759,6 @@ def extract_spectra( slitmap[f"ypix_{channel}"] = slitmap[f"ypix_{channel}"].astype("float64") slitmap[f"ypix_{channel}"][select_spec] += bn.nanmedian(shifts, axis=0) - if error is not None: - error[mask] = replace_error rss = RSS( data=data, mask=mask, From 8fa05b1c4abce80f4612256698b007f1e46a7977 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Fri, 13 Sep 2024 12:36:59 -0300 Subject: [PATCH 14/20] fixing bitmask propagation in channel combine method --- python/lvmdrp/core/rss.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/lvmdrp/core/rss.py b/python/lvmdrp/core/rss.py index 5ab05a8a..9c1c3e3b 100644 --- a/python/lvmdrp/core/rss.py +++ b/python/lvmdrp/core/rss.py @@ -476,7 +476,7 @@ def from_channels(cls, rss_b, rss_r, rss_z, use_weights=True): new_data = bn.nansum(fluxes * weights, axis=0) new_lsf = bn.nansum(lsfs * weights, axis=0) new_error = numpy.sqrt(bn.nansum(vars, axis=0)) - new_mask = numpy.bitwise_or(masks, axis=0) + new_mask = numpy.bitwise_or.reduce(masks, axis=0) if rss._sky is not None: new_sky = bn.nansum(skies * weights, axis=0) else: @@ -506,7 +506,7 @@ def from_channels(cls, rss_b, rss_r, rss_z, use_weights=True): new_data = bn.nanmean(fluxes, axis=0) new_lsf = bn.nanmean(lsfs, axis=0) new_error = numpy.sqrt(bn.nanmean(vars, axis=0)) - new_mask = numpy.bitwise_or(masks, axis=0) + new_mask = numpy.bitwise_or.reduce(masks, axis=0) if skies.size != 0: new_sky = bn.nansum(skies, axis=0) else: From 5175d5e269bd43d4257a366db3a7e2d74b3a4cd7 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Fri, 13 Sep 2024 12:38:24 -0300 Subject: [PATCH 15/20] implementing bitmask propagation in high-level Image and RSS routines --- python/lvmdrp/functions/imageMethod.py | 88 +++++++++++++------------- python/lvmdrp/functions/rssMethod.py | 24 +++---- 2 files changed, 57 insertions(+), 55 deletions(-) diff --git a/python/lvmdrp/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index d13f01ea..1fba2ab6 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -256,13 +256,13 @@ def _channel_combine_fiber_params(in_cent_traces, in_waves, add_overscan_columns if add_overscan_columns: os_region = numpy.zeros((mtraces[i]._data.shape[0], 34), dtype=int) trace_data = numpy.split(mtraces[i]._data, 2, axis=1) - trace_mask = numpy.split(mtraces[i]._mask, 2, axis=1) + trace_mask = numpy.split(mtraces[i]._mask.astype(bool), 2, axis=1) mtraces[i]._data = numpy.concatenate([trace_data[0], os_region, trace_data[1]], axis=1) - mtraces[i]._mask = numpy.concatenate([trace_mask[0], ~os_region.astype(bool), trace_mask[1]], axis=1) + mtraces[i]._mask = numpy.concatenate([trace_mask[0], ~os_region.astype("bool"), trace_mask[1]], axis=1) wave_data = numpy.split(mwaves[i]._data, 2, axis=1) wave_mask = numpy.split(mwaves[i]._mask, 2, axis=1) mwaves[i]._data = numpy.concatenate([wave_data[0], os_region, wave_data[1]], axis=1) - mwaves[i]._mask = numpy.concatenate([wave_mask[0], ~os_region.astype(bool), wave_mask[1]], axis=1) + mwaves[i]._mask = numpy.concatenate([wave_mask[0], ~os_region.astype("bool"), wave_mask[1]], axis=1) # get pixels where spectropgraph is not blocked channel_mask = (SPEC_CHANNELS[channels[i]][0]>=mwaves[i]._data)|(mwaves[i]._data>=SPEC_CHANNELS[channels[i]][1]) @@ -298,13 +298,13 @@ def _get_fiber_selection(traces, image_shape=(4080, 4120), y_widths=3): numpy.ndarray fiber selection mask """ - images = [Image(data=numpy.zeros(image_shape), mask=numpy.zeros(image_shape, dtype=bool)) for _ in range(len(traces))] + images = [Image(data=numpy.zeros(image_shape), mask=numpy.zeros(image_shape, dtype=numpy.int32)) for _ in range(len(traces))] for i in range(len(traces)): images[i].maskFiberTraces(traces[i], aperture=y_widths, parallel=1) image = Image() image.unsplit(images) - return image._mask + return image._mask != 0 def _get_wave_selection(waves, lines_list, window): @@ -635,7 +635,7 @@ def select_lines_2d(in_images, out_mask, in_cent_traces, in_waves, lines_list=No # get lines selection mask log.info(f"selecting lines with {wave_widths = } angstrom") lines_mask = _get_wave_selection(mwave._data, lines_list=lines_list, window=wave_widths) - lines_mask &= ~mwave._mask + lines_mask &= (mwave._mask == 0) # make sky mask 2d log.info("interpolating mask into 2D frame") @@ -892,7 +892,7 @@ def addCCDMask_drp(image, mask, replaceError="1e10"): img = loadImage(image) bad_pixel = loadImage(mask, extension_mask=0) if img._mask is not None: - mask_comb = numpy.logical_or(img._mask, bad_pixel._mask) + mask_comb = numpy.bitwise_or(img._mask, bad_pixel._mask) else: mask_comb = bad_pixel._mask if img._error is not None: @@ -1650,12 +1650,12 @@ def trace_peaks( trace_data = copy(trace) # mask zeros and data outside threshold and max_diff - trace._mask |= (trace._data <= 0) + trace._mask |= (trace._data <= 0) * PixMask["BADPIX"] # smooth all trace by a polynomial log.info(f"fitting trace with {numpy.abs(poly_disp)}-deg polynomial") table, table_poly, table_poly_all = trace.fit_polynomial(poly_disp, poly_kind="poly") # set bad fibers in trace mask - trace._mask[bad_fibers] = True + trace._mask[bad_fibers] = PixMask["BADTRACE"] | PixMask["WEAKFIBER"] | PixMask["NONEXPOSED"] if write_trace_data: _create_trace_regions(out_trace, table, table_poly, table_poly_all, display_plots=display_plots) @@ -1665,7 +1665,7 @@ def trace_peaks( x_pixels = numpy.arange(trace._data.shape[1]) y_pixels = numpy.arange(trace._fibers) for column in range(trace._data.shape[1]): - mask = trace._mask[:, column] + mask = trace._mask[:, column] != 0 for order in range(trace._coeffs.shape[1]): trace._coeffs[mask, order] = numpy.interp(y_pixels[mask], y_pixels[~mask], trace._coeffs[~mask, order]) # evaluate trace at interpolated fibers @@ -1961,8 +1961,8 @@ def subtract_straylight( # update mask if img_median._mask is None: - img_median._mask = numpy.zeros(img_median._data.shape, dtype=bool) - img_median._mask = img_median._mask | numpy.isnan(img_median._data) | numpy.isinf(img_median._data) | (img_median._data == 0) + img_median._mask = numpy.zeros(img_median._data.shape, dtype=numpy.int32) + img_median._mask |= ((numpy.isnan(img_median._data) | numpy.isinf(img_median._data) | (img_median._data == 0)) * PixMask["BADPIX"]) # mask regions around each fiber within a given cross-dispersion aperture log.info(f"masking fibers with an aperture of {aperture} pixels") @@ -1981,11 +1981,11 @@ def subtract_straylight( for icol in range(img_median._dim[1]): # mask top/bottom rows before/after first/last fiber - img_median._mask[tfiber[icol]:, icol] = True - img_median._mask[:bfiber[icol], icol] = True + img_median._mask[tfiber[icol]:, icol] = PixMask["NODATA"] + img_median._mask[:bfiber[icol], icol] = PixMask["NODATA"] # unmask select_nrows around each region - img_median._mask[(tfiber[icol]+aperture//2):(tfiber[icol]+aperture//2+select_tnrows), icol] = False - img_median._mask[(bfiber[icol]-aperture//2-select_bnrows):(bfiber[icol]-aperture//2), icol] = False + img_median._mask[(tfiber[icol]+aperture//2):(tfiber[icol]+aperture//2+select_tnrows), icol] = 0 + img_median._mask[(bfiber[icol]-aperture//2-select_bnrows):(bfiber[icol]-aperture//2), icol] = 0 # fit the signal in unmaksed areas along cross-dispersion axis by a polynomial log.info(f"fitting spline with {smoothing = } to the background signal along cross-dispersion axis") @@ -2155,7 +2155,7 @@ def traceFWHM_drp( median_cross = max(median_cross, 1) img = img.medianImg((median_cross, median_box)) - img._mask[...] = False + img._mask[...] = 0 # plt.figure(figsize=(20,10)) # plt.plot(img.getSlice(1300, axis="y")._error.tolist(), lw=0.6, color="0.7") @@ -2185,7 +2185,7 @@ def traceFWHM_drp( trace_mask.loadFitsData(in_trace) orig_trace = copy(trace_mask) - trace_mask._mask[...] = False + trace_mask._mask[...] = 0 # create a trace mask for the image traceFWHM = TraceMask() @@ -2234,13 +2234,13 @@ def traceFWHM_drp( fwhm, mask = img.traceFWHM(select_steps, trace_mask, blocks, init_fwhm, threshold_flux, max_pix=dim[0]) for ifiber in range(orig_trace._fibers): - if orig_trace._mask[ifiber].all(): + if (orig_trace._mask[ifiber] != 0).all(): continue good_pix = (~mask[ifiber]) & (~numpy.isnan(fwhm[ifiber])) & (fwhm[ifiber] != 0.0) & ((clip[0]= 2**16 - proc_img._mask += saturated_mask * PixMask["SATURATED"] + proc_img.add_bitmask("SATURATED", where=saturated_mask) # log number of masked pixels nmasked = numpy.count_nonzero(proc_img._mask) @@ -3830,7 +3831,7 @@ def detrend_frame( # gain-correct quadrant quad *= gain_map # propagate new NaNs to the mask - quad._mask += numpy.isnan(quad._data) * PixMask["BADPIX"] + quad.add_bitmask("BADPIX", where=numpy.isnan(quad._data)) quad.computePoissonError(rdnoise) bcorr_img.setSection(section=quad_sec, subimg=quad, inplace=True) @@ -3858,7 +3859,7 @@ def detrend_frame( # propagate pixel mask log.info("propagating pixel mask") nanpixels = ~numpy.isfinite(detrended_img._data) - detrended_img._mask += numpy.where(detrended_img._mask != nanpixels * PixMask["BADPIX"], detrended_img._mask + nanpixels * PixMask["BADPIX"], detrended_img._mask) + detrended_img.add_bitmask("BADPIX", where=nanpixels) # reject cosmic rays if reject_cr: @@ -3877,7 +3878,7 @@ def detrend_frame( # 'pixflat' is the imagetyp that a pixel flat can have if img_type == "pixflat": flat_array = numpy.ma.masked_array( - detrended_img._data, mask=detrended_img._mask + detrended_img._data, mask=detrended_img._mask != 0 ) detrended_img = detrended_img / numpy.ma.median(flat_array) @@ -4216,7 +4217,7 @@ def create_pixelmask(in_short_dark, in_long_dark, out_pixmask, in_flat_a=None, i sections = short_dark.getHdrValue("AMP? TRIMSEC") # define pixelmask image - pixmask = Image(data=numpy.zeros_like(short_dark._data), mask=numpy.zeros_like(short_dark._data, dtype="bool")) + pixmask = Image(data=numpy.zeros_like(short_dark._data), mask=numpy.zeros_like(short_dark._data, dtype="int32")) # create histogram figure fig_dis, axs_dis = create_subplots(to_display=display_plots, nrows=2, ncols=2, figsize=(15,10), sharex=True, sharey=True) @@ -4315,11 +4316,12 @@ def create_pixelmask(in_short_dark, in_long_dark, out_pixmask, in_flat_a=None, i log.info(f"masking {flat_mask.sum()} pixels with flats ratio < {flat_low} or > {flat_high}") # update pixel mask - pixmask.setData(mask=(pixmask._mask | flat_mask), inplace=True) + pixmask.add_bitmask("BADPIX", where=flat_mask) # set masked pixels to NaN - log.info(f"masked {pixmask._mask.sum()} pixels in total ({pixmask._mask.sum()/pixmask._mask.size*100:.2f}%)") - pixmask.setData(data=numpy.nan, select=pixmask._mask) + nmasked = numpy.count_nonzero(pixmask._mask) + log.info(f"masked {nmasked} pixels in total ({nmasked/pixmask._mask.size*100:.2f}%)") + pixmask.setData(data=numpy.nan) # write output mask log.info(f"writing pixel mask to '{os.path.basename(out_pixmask)}'") @@ -4373,7 +4375,7 @@ def trace_centroids(in_image: str, counts_threshold = counts_threshold * coadd # handle invalid error values - img._error[img._mask|(img._error<=0)] = numpy.inf + img._error[(img._mask!=0)|(img._error<=0)] = numpy.inf # calculate centroids for reference column if correct_ref: @@ -4394,7 +4396,7 @@ def trace_centroids(in_image: str, _create_trace_regions(out_trace_cent, table_data, table_poly, table_poly_all, display_plots=display_plots) # set bad fibers in trace mask - centroids._mask[bad_fibers] = True + centroids._mask[bad_fibers] = PixMask["DEADFIBER"] # linearly interpolate coefficients at masked fibers log.info(f"interpolating coefficients at {bad_fibers.sum()} masked fibers") centroids.interpolate_coeffs() @@ -4403,7 +4405,7 @@ def trace_centroids(in_image: str, centroids.interpolate_data(axis="X") # set bad fibers in trace mask - centroids._mask[bad_fibers] = True + centroids._mask[bad_fibers] = PixMask["DEADFIBER"] log.info(f"interpolating data at {bad_fibers.sum()} masked fibers") centroids.interpolate_data(axis="Y") @@ -4569,7 +4571,7 @@ def trace_fibers( counts_threshold = counts_threshold * coadd # handle invalid error values - img._error[img._mask|(img._error<=0)] = numpy.inf + img._error[(img._mask != 0)|(img._error<=0)] = numpy.inf # trace centroids in each column log.info(f"loading guess fiber centroids from '{os.path.basename(in_trace_cent_guess)}'") @@ -4614,9 +4616,9 @@ def trace_fibers( _create_trace_regions(out_trace_fwhm, table_data, table_poly, table_poly_all, display_plots=display_plots) # set bad fibers in trace mask - trace_amp._mask[bad_fibers] = True - trace_cent._mask[bad_fibers] = True - trace_fwhm._mask[bad_fibers] = True + trace_amp._mask[bad_fibers] = PixMask["DEADFIBER"] + trace_cent._mask[bad_fibers] = PixMask["DEADFIBER"] + trace_fwhm._mask[bad_fibers] = PixMask["DEADFIBER"] # linearly interpolate coefficients at masked fibers if interpolate_missing: @@ -4631,9 +4633,9 @@ def trace_fibers( trace_cent.interpolate_data(axis="X") trace_fwhm.interpolate_data(axis="X") # set bad fibers in trace mask - trace_amp._mask[bad_fibers] = True - trace_cent._mask[bad_fibers] = True - trace_fwhm._mask[bad_fibers] = True + trace_amp._mask[bad_fibers] = PixMask["DEADFIBER"] + trace_cent._mask[bad_fibers] = PixMask["DEADFIBER"] + trace_fwhm._mask[bad_fibers] = PixMask["DEADFIBER"] if interpolate_missing: log.info(f"interpolating data at {bad_fibers.sum()} masked fibers") @@ -4691,7 +4693,7 @@ def trace_fibers( img_slice = img_.getSlice(icolumn, axis="y") joint_mod = mod_columns[i](img_slice._pixels) - img_slice._data[(img_slice._mask)|(joint_mod<=0)] = numpy.nan + img_slice._data[(img_slice._mask!=0)|(joint_mod<=0)] = numpy.nan weights = img_slice._data / bn.nansum(img_slice._data) * 500 residuals = (joint_mod - img_slice._data) / img_slice._data * 100 diff --git a/python/lvmdrp/functions/rssMethod.py b/python/lvmdrp/functions/rssMethod.py index c5ff92a5..39cae775 100644 --- a/python/lvmdrp/functions/rssMethod.py +++ b/python/lvmdrp/functions/rssMethod.py @@ -22,7 +22,7 @@ from scipy import interpolate, ndimage from lvmdrp.utils.decorators import skip_on_missing_input_path, skip_if_drpqual_flags -from lvmdrp.utils.bitmask import ReductionStage +from lvmdrp.utils.bitmask import ReductionStage, PixMask from lvmdrp.core.constants import CONFIG_PATH, ARC_LAMPS from lvmdrp.core.cube import Cube from lvmdrp.core.tracemask import TraceMask @@ -95,7 +95,7 @@ def _illumination_correction(fiberflat, apply_correction=True): # define data and set to NaN bad pixels data = copy(fiberflat._data) - data[(fiberflat._mask)|(data <= 0)] = numpy.nan + data[(fiberflat._mask != 0)|(data <= 0)] = numpy.nan # compute median factors sci_factor = bn.nanmedian(data[sci_fibers, 1000:3000]) @@ -318,7 +318,7 @@ def determine_wavelength_solution(in_arcs: List[str]|str, out_wave: str, out_lsf arc, _, _ = arc.subtract_continuum(niter=cont_niter, thresh=cont_thresh, median_box_range=cont_box_range) # replace NaNs - mask = arc._mask | numpy.isnan(arc._data) | numpy.isnan(arc._error) + mask = (arc._mask != 0) | numpy.isnan(arc._data) | numpy.isnan(arc._error) mask |= (arc._data < 0.0) | (arc._error <= 0.0) arc._data[mask] = 0.0 arc._error[mask] = numpy.inf @@ -697,7 +697,7 @@ def shift_wave_skylines(in_frame: str, out_frame: str, dwave: float = 8.0, skyli iterator = tqdm(range(lvmframe._fibers), total=lvmframe._fibers, desc=f"measuring offsets using {len(skylines)} sky line(s)", ascii=True, unit="fiber") for ifiber in iterator: spec = lvmframe.getSpec(ifiber) - if spec._mask.all() or lvmframe._slitmap[ifiber]["telescope"] == "Spec" or lvmframe._slitmap[ifiber]["fibstatus"] in [1, 2]: + if (spec._mask!=0).all() or lvmframe._slitmap[ifiber]["telescope"] == "Spec" or lvmframe._slitmap[ifiber]["fibstatus"] in [1, 2]: continue fwhm_guess = numpy.nanmean(numpy.interp(skylines, lvmframe._wave[ifiber], lvmframe._lsf[ifiber])) @@ -1324,7 +1324,7 @@ def create_fiberflat(in_rsss: List[str], out_rsss: List[str], median_box: int = log.info(f"computing fiberflat across {fiberflat._fibers} fibers and {(~numpy.isnan(norm)).sum()} wavelength bins") normalized = fiberflat._data / norm[None, :] fiberflat._data = normalized - fiberflat._mask |= normalized <= 0 + fiberflat._mask |= PixMask["BADFLAT"] # apply clipping if clip_range is not None: @@ -1350,19 +1350,19 @@ def create_fiberflat(in_rsss: List[str], out_rsss: List[str], median_box: int = # interpolate masked pixels in fiberflat for ifiber in range(fiberflat._fibers): - wave, data, mask = fiberflat._wave[ifiber], fiberflat._data[ifiber], fiberflat._mask[ifiber] + wave, data, mask = fiberflat._wave[ifiber], fiberflat._data[ifiber], fiberflat._mask[ifiber] != 0 mask |= ~numpy.isfinite(data) if numpy.sum(~mask) == 0: continue fiberflat._data[ifiber, mask] = interpolate.interp1d(wave[~mask], data[~mask], bounds_error=False, assume_sorted=True)(wave[mask]) - fiberflat._mask[ifiber, mask] = False + fiberflat._mask[ifiber, mask] = 0 # create diagnostic plots log.info("creating diagnostic plots for fiberflat") fig, axs = create_subplots(to_display=display_plots, nrows=3, ncols=1, figsize=(12, 15), sharex=True) # plot original continuum exposure, fiberflat and corrected fiberflat per fiber colors = plt.cm.Spectral(numpy.linspace(0, 1, fiberflat._fibers)) - rss._data[rss._mask] = numpy.nan + rss._data[rss._mask != 0] = numpy.nan stdev_ori = biweight_scale(rss._data, axis=0, ignore_nan=True)[1950:2050].mean() stdev_new = biweight_scale(rss._data/fiberflat._data, axis=0, ignore_nan=True)[1950:2050].mean() log.info(f"flatfield statistics: {stdev_ori = :.2f}, {stdev_new = :.2f} ({stdev_new/stdev_ori*100:.2f}%)") @@ -1414,7 +1414,7 @@ def create_fiberflat(in_rsss: List[str], out_rsss: List[str], median_box: int = # perform some statistic about the fiberflat if fiberflat_cam._mask is not None: - select = numpy.logical_not(fiberflat_cam._mask) + select = fiberflat_cam._mask == 0 else: select = fiberflat_cam._data == fiberflat_cam._data min = bn.nanmin(fiberflat_cam._data[select]) @@ -2450,7 +2450,7 @@ def maskFibers_drp(in_rss, out_rss, fibers, replace_error="1e10"): rss.writeFitsData(out_rss) -def maskNAN_drp(in_rss, replace_error="1e12"): +def maskNAN_drp(in_rss, bitmask, replace_error="1e12"): """mask NaN values in given RSS Parameters @@ -2469,7 +2469,7 @@ def maskNAN_drp(in_rss, replace_error="1e12"): if rss._error is not None: rss._error[i, :] = float(replace_error) if rss._mask is not None: - rss._mask[i, :] = True + rss._mask[i, :] = bitmask rss.writeFitsData(in_rss) @@ -3054,7 +3054,7 @@ def createMasterFiberFlat_drp( for ifiber in range(fiberflat._fibers): print(ifiber) fiber = fiberflat.getSpec(ifiber) - select = numpy.logical_not(fiber._mask) + select = fiber._mask == 0 wave_ori = fiber._wave norm = superflat_func(fiber._wave[select]) * (1 / fiber[select]) From af9f169b6319760c296883e32e8c4998c7917ed9 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Fri, 13 Sep 2024 12:38:59 -0300 Subject: [PATCH 16/20] implementing bitmask propagation to sky routines --- python/lvmdrp/core/sky.py | 8 +++++--- python/lvmdrp/functions/skyMethod.py | 23 +++++++++++------------ python/lvmdrp/utils/bitmask.py | 2 +- 3 files changed, 17 insertions(+), 16 deletions(-) diff --git a/python/lvmdrp/core/sky.py b/python/lvmdrp/core/sky.py index af394bf3..9fae7cf5 100644 --- a/python/lvmdrp/core/sky.py +++ b/python/lvmdrp/core/sky.py @@ -43,6 +43,7 @@ SKYMODEL_INST_PATH, SKYMODEL_MODEL_CONFIG_PATH, ) +from lvmdrp.utils.bitmask import PixMask from lvmdrp.external.skycorr import createParFile, fitstabSkyCorrWrapper, runSkyCorr from lvmdrp import log @@ -713,11 +714,12 @@ def fit_supersky(sky_wave, sky_data, sky_vars, sky_mask, sci_wave, sci_data): weights = 1 / svars # define interpolation functions + good_pix = smask == 0 # NOTE: store a super sampled version of the splines as an extension of the sky RSS - f_data = interpolate.make_smoothing_spline(swave[~smask], ssky[~smask], w=weights[~smask], lam=1e-6) + f_data = interpolate.make_smoothing_spline(swave[good_pix], ssky[good_pix], w=weights[good_pix], lam=1e-6) # NOTE: verify that the evaluated errors are not in variance at this stage - f_error = interpolate.make_smoothing_spline(swave[~smask], svars[~smask] / nsky_fibers, w=weights[~smask], lam=1e-6) - f_mask = interpolate.interp1d(swave, smask, kind="nearest", bounds_error=False, fill_value=0) + f_error = interpolate.make_smoothing_spline(swave[good_pix], svars[good_pix] / nsky_fibers, w=weights[good_pix], lam=1e-6) + f_mask = interpolate.interp1d(swave, smask, kind="nearest", bounds_error=False, fill_value=PixMask["NODATA"]) return f_data, f_error, f_mask, swave, ssky, svars, smask diff --git a/python/lvmdrp/functions/skyMethod.py b/python/lvmdrp/functions/skyMethod.py index 124faaf1..414df022 100644 --- a/python/lvmdrp/functions/skyMethod.py +++ b/python/lvmdrp/functions/skyMethod.py @@ -29,7 +29,7 @@ SKYMODEL_CONFIG_PATH, SKYMODEL_INST_PATH, ) -from lvmdrp.utils.bitmask import ReductionStage +from lvmdrp.utils.bitmask import ReductionStage, PixMask from lvmdrp.core.plot import plt, create_subplots, save_fig from lvmdrp.core.header import Header from lvmdrp.core.passband import PassBand @@ -1396,10 +1396,10 @@ def interpolate_sky( in_frame: str, out_rss: str = None, display_plots: bool = F dlambda = np.column_stack((dlambda, dlambda[:, -1])) new_sky = s_ssky(frame._wave).astype("float32") * dlambda new_error = np.sqrt(s_error(frame._wave).astype("float32")) * dlambda - new_mask = s_mask(frame._wave).astype(bool) + new_mask = s_mask(frame._wave).astype("int32") # update mask with new bad pixels - new_mask = (new_sky < 0) | np.isnan(new_sky) - new_mask |= (new_error < 0) | np.isnan(new_error) + new_mask |= ((new_sky < 0) | np.isnan(new_sky)) * PixMask["BADSKYFIT"] + new_mask |= ((new_error < 0) | np.isnan(new_error)) * PixMask["BADSKYFIT"] # store outputs supersky[telescope] = s_ssky @@ -1658,9 +1658,8 @@ def quick_sky_subtraction(in_cframe, out_sframe, # print("writing lvmSFrame") log.info(f"writing lvmSframe to {out_sframe}") - sframe = lvmSFrame(data=skysub_data, error=skysub_error, mask=cframe._mask.astype(bool), sky=skydata, sky_error=sky_error, + sframe = lvmSFrame(data=skysub_data, error=skysub_error, mask=cframe._mask, sky=skydata, sky_error=sky_error, wave=cframe._wave, lsf=cframe._lsf, header=cframe._header, slitmap=cframe._slitmap) - sframe._mask |= ~np.isfinite(sframe._error) sframe._header["DRPSTAGE"] = (ReductionStage(sframe._header["DRPSTAGE"]) + "SKY_SUBTRACTED").value sframe.writeFitsData(out_sframe) # TODO: check on expnum=7632 for halpha emission in sky fibers @@ -1922,22 +1921,22 @@ def create_skysub_spectrum(hdu: fits.HDUList, tel: str, # do continuum vs line separation log.info("separating line and continuum in Sci, SkyE ,and SkyW spectra") - specsci=Spectrum1D(wave=Table(hdusci.data)['WAVE'].data, data=Table(hdusci.data)['FLUX'].data, mask=np.zeros(len(Table(hdusci.data)), dtype=bool)) - specskye=Spectrum1D(wave=Table(hduskye.data)['WAVE'].data, data=Table(hduskye.data)['FLUX'].data, mask=np.zeros(len(Table(hduskye.data)), dtype=bool)) - specskyw=Spectrum1D(wave=Table(hduskyw.data)['WAVE'].data, data=Table(hduskyw.data)['FLUX'].data, mask=np.zeros(len(Table(hduskyw.data)), dtype=bool)) + specsci=Spectrum1D(wave=Table(hdusci.data)['WAVE'].data, data=Table(hdusci.data)['FLUX'].data, mask=np.zeros(len(Table(hdusci.data)), dtype=np.int32)) + specskye=Spectrum1D(wave=Table(hduskye.data)['WAVE'].data, data=Table(hduskye.data)['FLUX'].data, mask=np.zeros(len(Table(hduskye.data)), dtype=np.int32)) + specskyw=Spectrum1D(wave=Table(hduskyw.data)['WAVE'].data, data=Table(hduskyw.data)['FLUX'].data, mask=np.zeros(len(Table(hduskyw.data)), dtype=np.int32)) # specsci_error = Spectrum1D( # wave=Table(hdusci.data)["WAVE"].data, # data=Table(hdusci.data)['ERROR'].data, - # mask=np.zeros(len(Table(hdusci.data)), dtype=bool)) + # mask=np.zeros(len(Table(hdusci.data)), dtype=np.int32)) specskye_error = Spectrum1D( wave=Table(hduskye.data)["WAVE"].data, data=Table(hduskye.data)['ERROR'].data, - mask=np.zeros(len(Table(hduskye.data)), dtype=bool)) + mask=np.zeros(len(Table(hduskye.data)), dtype=np.int32)) specskyw_error = Spectrum1D( wave=Table(hduskyw.data)["WAVE"].data, data=Table(hduskyw.data)['ERROR'].data, - mask=np.zeros(len(Table(hduskyw.data)), dtype=bool)) + mask=np.zeros(len(Table(hduskyw.data)), dtype=np.int32)) csci, cscimask = find_continuum(specsci._data) cskye, cskyemask = find_continuum(specskye._data) diff --git a/python/lvmdrp/utils/bitmask.py b/python/lvmdrp/utils/bitmask.py index d49068be..8b6cf92c 100644 --- a/python/lvmdrp/utils/bitmask.py +++ b/python/lvmdrp/utils/bitmask.py @@ -247,7 +247,7 @@ class PixMask(BaseBitmask): BADFLUXFACTOR = auto() # large sky residuals - BADSKYCHI = auto() + BADSKYFIT = auto() # fiber bitmasks ------------------------------------------------------ NONEXPOSED = auto() From a4fb724a6aad0740bdc2f5e25b1410718755718d Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Fri, 13 Sep 2024 12:58:37 -0300 Subject: [PATCH 17/20] cleaning lingering debug print --- python/lvmdrp/functions/imageMethod.py | 1 - 1 file changed, 1 deletion(-) diff --git a/python/lvmdrp/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index 1fba2ab6..9592bbf7 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -3402,7 +3402,6 @@ def preproc_raw_frame( # create pixel mask on the original image log.info("building pixel mask") - print(proc_img._mask) proc_img.add_bitmask("BADPIX", where=master_mask) # convert temp image to ADU for saturated pixel masking saturated_mask = proc_img._data >= 2**16 From 857f779b2b5c30e1a3c2eb1ef17b16b2f0d45c1a Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Fri, 13 Sep 2024 15:57:45 -0300 Subject: [PATCH 18/20] adding missing reduction stage & improving/adding bitmask handling routines --- python/lvmdrp/utils/bitmask.py | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/python/lvmdrp/utils/bitmask.py b/python/lvmdrp/utils/bitmask.py index 8b6cf92c..8a938b83 100644 --- a/python/lvmdrp/utils/bitmask.py +++ b/python/lvmdrp/utils/bitmask.py @@ -147,6 +147,7 @@ class ReductionStage(BaseBitmask): HDRFIX_APPLIED = auto() # header fix applied OVERSCAN_SUBTRACTED = auto() # trimmed overscan region, overscan-subtracted PIXELMASK_ADDED = auto() # fiducial pixel mask was added + PREPROCESSED = auto() GAIN_CORRECTED = auto() # gain correction applied POISSON_ERROR_CALCULATED = auto() # calculated poisson error LINEARITY_CORRECTED = auto() # linearity correction applied @@ -287,19 +288,17 @@ class PixMask(BaseBitmask): DRPQUALITIES = list() -def _parse_bitmask(pixmask, mask_shape): - if isinstance(pixmask, PixMask): - pixmask = np.zeros(mask_shape) - elif isinstance(pixmask, str): - pixmask = PixMask[pixmask] - elif isinstance(pixmask, int): - pixmask = PixMask(pixmask) - elif isinstance(pixmask, np.ndarray[int]): - assert pixmask.shape == mask_shape, f"Wrong `pixmask` shape {pixmask.shape} not matching `mask_image` shape {mask_shape}" +def _parse_bitmask(bitmask, kind): + if isinstance(bitmask, kind): + return bitmask + elif isinstance(bitmask, str): + bitmask = kind[bitmask] + elif isinstance(bitmask, int): + bitmask = kind(bitmask) else: - raise ValueError(f"Wrong type for {pixmask = }: {type(pixmask)}; expected PixMask, string or integer") + raise ValueError(f"Wrong type for {bitmask = }: {type(bitmask)}; expected {kind}, string or integer") - return pixmask + return bitmask def _parse_where(where, mask_shape): @@ -314,7 +313,7 @@ def _parse_where(where, mask_shape): def add_bitmask(mask_image, pixmask, where=None): - pixmask = _parse_bitmask(pixmask, mask_shape=mask_image.shape) + pixmask = _parse_bitmask(pixmask, kind=PixMask) where = _parse_where(where, mask_shape=mask_image.shape) mask_image[where] |= pixmask @@ -322,13 +321,22 @@ def add_bitmask(mask_image, pixmask, where=None): def toggle_bitmask(mask_image, pixmask, where=None): - pixmask = _parse_bitmask(pixmask, mask_shape=mask_image.shape) + pixmask = _parse_bitmask(pixmask, kind=PixMask) where = _parse_where(where, mask_shape=mask_image.shape) mask_image[where] ^= pixmask return mask_image +def update_header_bitmask(header, kind, bitmask, key, comment): + bitmask = _parse_bitmask(bitmask, kind=kind) + if key not in header: + header[key] = (bitmask.value, comment) + header[key] |= bitmask.value + + return header + + def print_bitmasks(mask_array, logger=None): uniques, counts = np.unique(mask_array, return_counts=True) bitmasks = dict(zip(map(lambda p: PixMask(int(p)).get_name() if p>0 else "GOODPIX", uniques), counts)) From b89069d29336d17a385badf4156455f98a85a4a8 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Fri, 13 Sep 2024 15:59:13 -0300 Subject: [PATCH 19/20] improving propagation of header bitmasks --- python/lvmdrp/core/image.py | 24 +++++++--- python/lvmdrp/core/rss.py | 12 ++++- python/lvmdrp/functions/imageMethod.py | 66 ++++++++++---------------- python/lvmdrp/functions/rssMethod.py | 12 ++--- 4 files changed, 59 insertions(+), 55 deletions(-) diff --git a/python/lvmdrp/core/image.py b/python/lvmdrp/core/image.py index 91ad857b..77c977d4 100644 --- a/python/lvmdrp/core/image.py +++ b/python/lvmdrp/core/image.py @@ -17,7 +17,7 @@ from scipy import interpolate from lvmdrp import log -from lvmdrp.utils.bitmask import PixMask, add_bitmask, toggle_bitmask, print_bitmasks +from lvmdrp.utils.bitmask import ReductionStage, PixMask, QualityFlag, add_bitmask, toggle_bitmask, update_header_bitmask, print_bitmasks from lvmdrp.core.constants import CON_LAMPS, ARC_LAMPS from lvmdrp.core.plot import plt from lvmdrp.core.fit_profile import gaussians, Gaussians @@ -998,12 +998,6 @@ def getSlice(self, slice, axis): numpy.arange(self._dim[0]), self._data[:, slice], error, mask ) - def set_mask(self, pixmask): - if isinstance(pixmask, numpy.ndarray): - assert isinstance(pixmask[0,0], PixMask) - - self._mask = pixmask - def setData( self, data=None, error=None, mask=None, header=None, select=None, inplace=True ): @@ -1119,6 +1113,22 @@ def print_bitmasks(self, logger=None): return print_bitmasks(self._mask, logger=log) + def update_drpstage(self, stage): + if self._header is None: + return + + self._header = update_header_bitmask(header=self._header, + kind=ReductionStage, bitmask=stage, + key="DRPSTAGE", comment="data reduction stage") + + def update_drpqual(self, quality): + if self._header is None: + return + + self._header = update_header_bitmask(header=self._header, + kind=QualityFlag, bitmask=quality, + key="DRPQUAL", comment="data reduction quality") + def convertUnit(self, to, assume="adu", gain_field="GAIN", inplace=False): """converts the unit of the image diff --git a/python/lvmdrp/core/rss.py b/python/lvmdrp/core/rss.py index 9c1c3e3b..00b7c9b1 100644 --- a/python/lvmdrp/core/rss.py +++ b/python/lvmdrp/core/rss.py @@ -14,7 +14,7 @@ from astropy import units as u from lvmdrp import log -from lvmdrp.utils.bitmask import PixMask, _parse_bitmask, add_bitmask, toggle_bitmask, print_bitmasks +from lvmdrp.utils.bitmask import ReductionStage, PixMask, _parse_bitmask, add_bitmask, toggle_bitmask, update_header_bitmask, print_bitmasks from lvmdrp.core.constants import CONFIG_PATH from lvmdrp.core.apertures import Aperture from lvmdrp.core.cube import Cube @@ -1013,6 +1013,14 @@ def print_bitmasks(self, logger=None): return print_bitmasks(self._mask, logger=log) + def update_drpstage(self, stage): + if self._header is None: + return + + self._header = update_header_bitmask(header=self._header, + kind=ReductionStage, bitmask=stage, + key="DRPSTAGE", comment="data reduction stage") + def add_header_comment(self, comstr): ''' Append a COMMENT card at the end of the FITS header. @@ -1186,7 +1194,7 @@ def match_lsf(self, target_fwhm=None, min_fwhm=0.1): return RSS.from_spectra1d(new_specs, header=self._header, slitmap=self._slitmap, good_fibers=self._good_fibers) def maskFiber(self, fiber, bitmask, replace_error=1e10): - bitmask = _parse_bitmask(bitmask) + bitmask = _parse_bitmask(bitmask, kind=PixMask) self._data[fiber, :] = 0 if self._mask is not None: self._mask[fiber, :] = bitmask diff --git a/python/lvmdrp/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index 9592bbf7..23e346b7 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -28,7 +28,7 @@ from lvmdrp import log, __version__ as DRPVER from lvmdrp.core.constants import CONFIG_PATH, SPEC_CHANNELS, ARC_LAMPS, LVM_REFERENCE_COLUMN, LVM_NBLOCKS, FIDUCIAL_PLATESCALE from lvmdrp.utils.decorators import skip_on_missing_input_path, drop_missing_input_paths -from lvmdrp.utils.bitmask import QualityFlag, ReductionStage, PixMask +from lvmdrp.utils.bitmask import PixMask from lvmdrp.core.fiberrows import FiberRows, _read_fiber_ypix from lvmdrp.core.image import ( Image, @@ -97,7 +97,7 @@ ] -def _nonlinearity_correction(ptc_params: None | numpy.ndarray, nominal_gain: float, quadrant: Image, iquad: int, drpstage: ReductionStage) -> Tuple[Image,ReductionStage]: +def _nonlinearity_correction(ptc_params: None | numpy.ndarray, nominal_gain: float, quadrant: Image, iquad: int) -> Tuple[Image]: """calculates non-linearity correction for input quadrant Parameters @@ -133,12 +133,12 @@ def _nonlinearity_correction(ptc_params: None | numpy.ndarray, nominal_gain: flo gain_med = bn.nanmedian(gain_map._data) gain_min, gain_max = bn.nanmin(gain_map._data), bn.nanmax(gain_map._data) log.info(f"gain map stats: {gain_med = :.2f} [{gain_min = :.2f}, {gain_max = :.2f}] ({nominal_gain = :.2f} e-/ADU)") - drpstage += "LINEARITY_CORRECTED" + quadrant.update_drpstage("LINEARITY_CORRECTED") else: log.warning("cannot apply non-linearity correction") log.info(f"using {nominal_gain = } (e-/ADU)") gain_map = Image(data=numpy.ones(quadrant._data.shape) * nominal_gain) - return gain_map, drpstage + return gain_map def _create_peaks_regions(fibermap: Table, column: int = 2000) -> None: @@ -391,7 +391,7 @@ def _fix_fiber_thermal_shifts(image, trace_cent, trace_width=None, trace_amp=Non # calculate thermal shifts column_shifts = image.measure_fiber_shifts(fiber_model, columns=columns, column_width=column_width, shift_range=shift_range, axs=axs_cc) if (column_shifts!=0).any(): - image._header["DRPSTAGE"] = (ReductionStage(image._header["DRPSTAGE"]) + "FIBERS_SHIFTED").value + image.update_drpstage("FIBERS_SHIFTED") # shifts stats median_shift = bn.nanmedian(column_shifts, axis=0) @@ -1941,7 +1941,6 @@ def subtract_straylight( # load image data log.info(f"using image {os.path.basename(in_image)} for stray light subtraction") img = loadImage(in_image) - drpstage = ReductionStage(img._header["DRPSTAGE"]) unit = img._header["BUNIT"] # smooth image along dispersion axis with a median filter excluded NaN values @@ -2002,12 +2001,11 @@ def subtract_straylight( log.info("subtracting the smoothed background signal from the original image") img_out = loadImage(in_image) img_out._data = img_out._data - img_stray._data - drpstage += "STRAYLIGHT_SUBTRACTED" # include header and write out file log.info(f"writing stray light subtracted image to {os.path.basename(out_image)}") img_out.setHeader(header=img.getHeader()) - img_out._header["DRPSTAGE"] = drpstage.value + img_out.update_drpstage("STRAYLIGHT_SUBTRACTED") img_out.writeFitsData(out_image) # plot results: polyomial fitting & smoothing, both with masked regions on @@ -2621,7 +2619,6 @@ def extract_spectra( log.info(f"extraction using aperture of {aperture} pixels") img = loadImage(in_image) - drpstage = ReductionStage(img._header["DRPSTAGE"]) mjd = img._header["SMJD"] camera = img._header["CCD"] expnum = img._header["EXPOSURE"] @@ -2752,7 +2749,6 @@ def extract_spectra( # print(numpy.unique(mask)) # exit() - drpstage += "SPECTRA_EXTRACTED" # propagate thermal shift to slitmap channel = img._header['CCD'][0] @@ -2769,7 +2765,7 @@ def extract_spectra( header=img.getHeader(), slitmap=slitmap ) - rss._header["DRPSTAGE"] = drpstage.value + rss.update_drpstage("SPECTRA_EXTRACTED") rss.setHdrValue("NAXIS2", data.shape[0]) rss.setHdrValue("NAXIS1", data.shape[1]) rss.setHdrValue("DISPAXIS", 1) @@ -3168,12 +3164,10 @@ def preproc_raw_frame( display_plots : bool, optional whether to show plots on display or not, by default False """ - # initialize reduction status - drpstage = ReductionStage(1) - # load image log.info(f"starting preprocessing of raw image '{os.path.basename(in_image)}'") org_img = loadImage(in_image) + org_img.update_drpqual(0) org_header = org_img.getHeader() camera = org_header["CCD"] @@ -3184,7 +3178,7 @@ def preproc_raw_frame( sjd = int(dateobs_to_sjd(org_header.get("OBSTIME"))) sjd = correct_sjd(in_image, sjd) org_header = apply_hdrfix(sjd, hdr=org_header) or org_header - drpstage += "HDRFIX_APPLIED" + org_img.update_drpstage("HDRFIX_APPLIED") except ValueError as e: log.error(f"cannot apply header fix: {e}") @@ -3300,13 +3294,11 @@ def preproc_raw_frame( log.info(f"masked {os_nmask} ({os_nmask/os_data.size*100:.2f}%) pixels in overscan above {overscan_threshold} standard deviations") sc_quad = sc_quad - os_model + sc_quad.update_drpstage("OVERSCAN_SUBTRACTED") os_profiles.append(os_profile) os_models.append(os_model) - if subtract_overscan: - drpstage += "OVERSCAN_SUBTRACTED" - # compute overscan stats os_bias_med[i] = bn.nanmedian(os_quad._data, axis=None) os_bias_std[i] = bn.nanmedian(bn.nanstd(os_quad._data, axis=1), axis=None) @@ -3396,13 +3388,13 @@ def preproc_raw_frame( if in_mask and proc_img._header["IMAGETYP"] not in {"bias", "dark", "pixflat"}: log.info(f"loading master pixel mask from {os.path.basename(in_mask)}") master_mask = loadImage(in_mask)._mask.astype(bool) - drpstage += "PIXELMASK_ADDED" else: master_mask = numpy.zeros_like(proc_img._data, dtype=bool) # create pixel mask on the original image log.info("building pixel mask") proc_img.add_bitmask("BADPIX", where=master_mask) + proc_img.update_drpstage("PIXELMASK_ADDED") # convert temp image to ADU for saturated pixel masking saturated_mask = proc_img._data >= 2**16 proc_img.add_bitmask("SATURATED", where=saturated_mask) @@ -3417,15 +3409,13 @@ def preproc_raw_frame( proc_img.apply_pixelmask() # update data reduction quality flag - drpqual = QualityFlag(0) if saturated_mask.sum() / proc_img._mask.size > 0.01: - drpqual += "SATURATED" + proc_img.update_drpqual("SATURATED") # write out FITS file log.info(f"writing preprocessed image to {os.path.basename(out_image)}") + proc_img.update_drpstage("PREPROCESSED") proc_img.setHdrValue("DRPVER", DRPVER, comment='data reduction pipeline version') - proc_img.setHdrValue("DRPSTAGE", drpstage.value, comment="data reduction stage") - proc_img.setHdrValue("DRPQUAL", value=drpqual.value, comment="data reduction quality flag") proc_img.writeFitsData(out_image) # plot overscan strips along X and Y axes @@ -3576,7 +3566,6 @@ def add_astrometry( # reading slitmap org_img = loadImage(in_image) - drpstage = ReductionStage(org_img._header["DRPSTAGE"]) slitmap = org_img.getSlitmap() telescope=numpy.array(slitmap['telescope'].data) x=numpy.array(slitmap['xpmm'].data) @@ -3597,7 +3586,7 @@ def copy_guider_keyword(gdrhdr, keyword, img): comment = gdrhdr.comments[keyword] if inhdr else '' img.setHdrValue(f'HIERARCH GDRCOADD {keyword}', gdrhdr.get(keyword), comment) - def getobsparam(tel, drpstage): + def getobsparam(tel): if tel!='spec': if os.path.isfile(agcfiledir[tel]): mfagc=fits.open(agcfiledir[tel]) @@ -3645,17 +3634,17 @@ def getobsparam(tel, drpstage): log.warning(f"some astrometry keywords for telescope '{tel}' are missing: {RAobs = }, {DECobs = }, {PAobs = }") org_img.add_header_comment(f"no astromentry keywords '{tel}': {RAobs = }, {DECobs = }, {PAobs = }, using commanded") org_img.setHdrValue('ASTRMSRC', 'CMD position', comment='Source of astrometric solution: commanded position') - drpstage += f"{tel.upper()}_ASTROMETRY_ADDED" + org_img.update_drpstage(f"{tel.upper()}_ASTROMETRY_ADDED") else: RAobs=0 DECobs=0 PAobs=0 - return RAobs, DECobs, PAobs, drpstage + return RAobs, DECobs, PAobs - RAobs_sci, DECobs_sci, PAobs_sci, drpstage = getobsparam('sci', drpstage) - RAobs_skye, DECobs_skye, PAobs_skye, drpstage = getobsparam('skye', drpstage) - RAobs_skyw, DECobs_skyw, PAobs_skyw, drpstage = getobsparam('skyw', drpstage) - RAobs_spec, DECobs_spec, PAobs_spec, drpstage = getobsparam('spec', drpstage) + RAobs_sci, DECobs_sci, PAobs_sci = getobsparam('sci') + RAobs_skye, DECobs_skye, PAobs_skye = getobsparam('skye') + RAobs_skyw, DECobs_skyw, PAobs_skyw = getobsparam('skyw') + RAobs_spec, DECobs_spec, PAobs_spec = getobsparam('spec') # Create fake IFU image WCS object for each telescope focal plane and use it to calculate RA,DEC of each fiber telcoordsdir={'sci':(RAobs_sci, DECobs_sci, PAobs_sci), 'skye':(RAobs_skye, DECobs_skye, PAobs_skye), 'skyw':(RAobs_skyw, DECobs_skyw, PAobs_skyw), 'spec':(RAobs_spec, DECobs_spec, PAobs_spec)} @@ -3693,10 +3682,9 @@ def getfibradec(tel, platescale): slitmap['ra']=RAfib slitmap['dec']=DECfib org_img._slitmap=slitmap - drpstage += "FIBERS_ASTROMETRY_ADDED" + org_img.update_drpstage("FIBERS_ASTROMETRY_ADDED") log.info(f"writing RA,DEC to slitmap in image '{os.path.basename(out_image)}'") - org_img._header["DRPSTAGE"] = drpstage.value org_img.writeFitsData(out_image) @@ -3751,7 +3739,6 @@ def detrend_frame( org_img = loadImage(in_image) exptime = org_img._header["EXPTIME"] img_type = org_img._header["IMAGETYP"].lower() - drpstage = ReductionStage(org_img._header["DRPSTAGE"]) log.info( "target frame parameters: " f"MJD = {org_img._header['MJD']}, " @@ -3826,7 +3813,7 @@ def detrend_frame( rdnoise = quad.getHdrValue(f"AMP{i+1} RDNOISE") # non-linearity correction - gain_map, drpstage = _nonlinearity_correction(ptc_params, gain, quad, iquad=i+1, drpstage=drpstage) + gain_map = _nonlinearity_correction(ptc_params, gain, quad, iquad=i+1) # gain-correct quadrant quad *= gain_map # propagate new NaNs to the mask @@ -3837,8 +3824,8 @@ def detrend_frame( log.info(f"median error in quadrant {i+1}: {bn.nanmedian(quad._error):.2f} (e-)") bcorr_img.setHdrValue("BUNIT", "electron", "physical units of the image") - drpstage += "GAIN_CORRECTED" - drpstage += "POISSON_ERROR_CALCULATED" + bcorr_img.update_drpstage("GAIN_CORRECTED") + bcorr_img.update_drpstage("POISSON_ERROR_CALCULATED") else: # convert to ADU log.info("leaving original ADU units") @@ -3853,7 +3840,7 @@ def detrend_frame( detrended_img = (bcorr_img - mdark_img.convertUnit(to=bcorr_img._header["BUNIT"])) # NOTE: this is a hack to avoid the error propagation of the division in Image detrended_img._data = detrended_img._data / numpy.nan_to_num(mflat_img._data, nan=1.0) - drpstage += "DETRENDED" + detrended_img.update_drpstage("DETRENDED") # propagate pixel mask log.info("propagating pixel mask") @@ -3866,7 +3853,7 @@ def detrend_frame( rdnoise = detrended_img.getHdrValue("AMP1 RDNOISE") detrended_img.reject_cosmics(gain=1.0, rdnoise=rdnoise, rlim=1.3, iterations=5, fwhm_gauss=[2.75, 2.75], replace_box=[10,2], replace_error=1e6, verbose=True, inplace=True) - drpstage += "COSMIC_CLEANED" + detrended_img.update_drpstage("COSMIC_CLEANED") # replace masked pixels with NaNs if replace_with_nan: @@ -3891,7 +3878,6 @@ def detrend_frame( # save detrended image log.info(f"writing detrended image to '{os.path.basename(out_image)}'") - detrended_img._header["DRPSTAGE"] = drpstage.value detrended_img.writeFitsData(out_image) # show plots diff --git a/python/lvmdrp/functions/rssMethod.py b/python/lvmdrp/functions/rssMethod.py index 39cae775..72c85c2d 100644 --- a/python/lvmdrp/functions/rssMethod.py +++ b/python/lvmdrp/functions/rssMethod.py @@ -22,7 +22,7 @@ from scipy import interpolate, ndimage from lvmdrp.utils.decorators import skip_on_missing_input_path, skip_if_drpqual_flags -from lvmdrp.utils.bitmask import ReductionStage, PixMask +from lvmdrp.utils.bitmask import PixMask from lvmdrp.core.constants import CONFIG_PATH, ARC_LAMPS from lvmdrp.core.cube import Cube from lvmdrp.core.tracemask import TraceMask @@ -743,7 +743,7 @@ def shift_wave_skylines(in_frame: str, out_frame: str, dwave: float = 8.0, skyli # write updated wobject log.info(f"writing updated wobject file '{os.path.basename(out_frame)}'") - lvmframe._header["DRPSTAGE"] = (ReductionStage(lvmframe._header["DRPSTAGE"]) + "WAVELENGTH_SHIFTED").value + lvmframe.update_drpstage("WAVELENGTH_SHIFTED") lvmframe.writeFitsData(out_frame) # Make QA plots showing offsets for each sky line in each spectrograph @@ -826,7 +826,7 @@ def create_pixel_table(in_rss: str, out_rss: str, in_waves: str, in_lsfs: str, c log.info(f"heliocentric velocities [km/s]: {helio_rvs}") log.info(f"writing output RSS to {out_rss}") - rss._header["DRPSTAGE"] = (ReductionStage(rss._header["DRPSTAGE"]) + "WAVELENGTH_CALIBRATED").value + rss.update_drpstage("WAVELENGTH_CALIBRATED") rss.writeFitsData(out_rss) return rss @@ -1631,7 +1631,7 @@ def apply_fiberflat(in_rss: str, out_frame: str, in_flat: str, clip_below: float superflat=flat._data ) lvmframe.set_header(orig_header=rss._header, flatname=os.path.basename(in_flat), ifibvar=ifibvar, ffibvar=ffibvar) - lvmframe._header["DRPSTAGE"] = (ReductionStage(lvmframe._header["DRPSTAGE"]) + "FLATFIELDED").value + lvmframe.update_drpstage("FLATFIELDED") lvmframe.writeFitsData(out_frame) return rss, lvmframe @@ -2930,7 +2930,7 @@ def stack_spectrographs(in_rsss: List[str], out_rss: str) -> RSS: # write output log.info(f"writing stacked RSS to {os.path.basename(out_rss)}") - rss_out._header["DRPSTAGE"] = (ReductionStage(rss_out._header["DRPSTAGE"]) + "SPECTROGRAPH_STACKED").value + rss_out.update_drpstage("SPECTROGRAPH_STACKED") rss_out.writeFitsData(out_rss) return rss_out @@ -2974,7 +2974,7 @@ def join_spec_channels(in_fframes: List[str], out_cframe: str, use_weights: bool sky_west=new_rss._sky_west, sky_west_error=new_rss._sky_west_error, slitmap=new_rss._slitmap) - cframe._header["DRPSTAGE"] = (ReductionStage(cframe._header["DRPSTAGE"]) + "CHANNEL_COMBINED").value + cframe.update_drpstage("CHANNEL_COMBINED") # write output RSS if out_cframe is not None: From 56ed8e99bd4c6d5af424bb12b4bb9f964560e046 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Sun, 15 Sep 2024 17:49:46 -0300 Subject: [PATCH 20/20] updating overall reduction quality bitmasks --- python/lvmdrp/utils/bitmask.py | 39 +++++++++++++++++++++++++++++++--- 1 file changed, 36 insertions(+), 3 deletions(-) diff --git a/python/lvmdrp/utils/bitmask.py b/python/lvmdrp/utils/bitmask.py index 8a938b83..b1108ed9 100644 --- a/python/lvmdrp/utils/bitmask.py +++ b/python/lvmdrp/utils/bitmask.py @@ -193,8 +193,26 @@ def __add__(self, flag): return new | flag class QualityFlag(BaseBitmask): - # TODO: add flag for overscan quality - OSFEATURES = auto() # Overscan region has features. + # overall quality of the reduction + + # raw frame quality + SATURATED10P = auto() # 10% of saturated pixels in this frame. + SATURATED20P = auto() # 20% of saturated pixels in this frame. + SATURATED50P = auto() # 50% of saturated pixels in this frame. + + # OS quality + BADOVERSCANQ1 = auto() # bad overscan (too high/low values) + BADOVERSCANQ2 = auto() + BADOVERSCANQ3 = auto() + BADOVERSCANQ4 = auto() + BADOSMODELQ1 = auto() # bad overscan model + BADOSMODELQ2 = auto() + BADOSMODELQ3 = auto() + BADOSMODELQ4 = auto() + + # detrending quality + + EXTRACTBAD = auto() # Many bad values in extracted frame. EXTRACTBRIGHT = auto() # Extracted spectra abnormally bright. LOWEXPTIME = auto() # Exposure time less than 10 minutes. @@ -204,9 +222,25 @@ class QualityFlag(BaseBitmask): BADDITHER = auto() # Bad dither location information. ARCFOCUS = auto() # Bad focus on arc frames. RAMPAGINGBUNNY = auto() # Rampaging dust bunnies in IFU flats. + + + # flux calibration quality + FEWSTDSTARS = auto() + FEWSCISTARS = auto() + NOSTDSTARS = auto() + NOSCISTARS = auto() + + BADSENSITIVITY = auto() # measured sensitivity curve too different from instrument one + SCIZEROPOINT = auto() # too low/high zero point in science stars + + SKYSUBBAD = auto() # Bad sky subtraction. SKYSUBFAIL = auto() # Failed sky subtraction. + + FULLCLOUD = auto() # Completely cloudy exposure. + + # thermal shifts BADFLEXURE = auto() # Abnormally high flexure LSF correction. BGROOVEFAIL = auto() # Possible B-groove glue failure. RGROOVEFAIL = auto() # Possible R-groove glue failure. @@ -216,7 +250,6 @@ class QualityFlag(BaseBitmask): NOSPEC3 = auto() # No data from spec3. BLOWTORCH = auto() # Blowtorch artifact detected. SEVEREBT = auto() # Severe blowtorch artifact. - SATURATED = auto() # X% of saturated pixels in this frame. class PixMask(BaseBitmask):