diff --git a/python/lvmdrp/core/fiberrows.py b/python/lvmdrp/core/fiberrows.py index 9be38c2f..fc1080fa 100644 --- a/python/lvmdrp/core/fiberrows.py +++ b/python/lvmdrp/core/fiberrows.py @@ -8,6 +8,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 @@ -200,7 +201,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) @@ -295,7 +296,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) @@ -415,7 +416,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) @@ -555,7 +556,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) @@ -698,6 +699,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: @@ -716,7 +731,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 @@ -837,8 +852,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") @@ -848,8 +863,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: @@ -908,7 +923,7 @@ def measureArcLines( ) for i in iterator: spec = self.getSpec(i) - if spec._mask.all(): + if (spec._mask!=0).all(): masked[i] = True continue @@ -952,7 +967,7 @@ def measureArcLines( ) for i in iterator: spec = self.getSpec(i) - if spec._mask.all(): + if (spec._mask!=0).all(): masked[i] = True continue @@ -1061,7 +1076,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) @@ -1073,8 +1088,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 @@ -1082,9 +1097,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) @@ -1111,7 +1126,7 @@ def fit_polynomial(self, deg, poly_kind="poly", clip=None, min_samples_frac=0.0) 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 good_sam = numpy.isfinite(samples[i, :]) if numpy.sum(good_pix) >= deg + 1 and good_sam.sum() / good_sam.size > min_samples_frac: # select the polynomial class @@ -1124,8 +1139,8 @@ def fit_polynomial(self, deg, poly_kind="poly", clip=None, min_samples_frac=0.0) 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 @@ -1133,10 +1148,10 @@ def fit_polynomial(self, deg, poly_kind="poly", clip=None, min_samples_frac=0.0) if clip is not None: self._data = numpy.clip(self._data, clip[0], clip[1]) - self._mask[i, :] = False + self._mask[i, :] = 0 else: log.warning(f"fiber {i} does not meet criteria: {good_pix.sum() = } >= {deg + 1 = } or {good_sam.sum() / good_sam.size = } > {min_samples_frac = }") - self._mask[i, :] = True + self._mask[i, :] = PixMask["FAILEDPOLY"] return numpy.asarray(pix_table), numpy.asarray(poly_table), numpy.asarray(poly_all_table) @@ -1197,7 +1212,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 @@ -1226,7 +1241,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 @@ -1362,7 +1377,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 @@ -1378,7 +1393,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 @@ -1411,7 +1426,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") @@ -1421,14 +1436,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: @@ -1439,7 +1454,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 09074036..03a27758 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 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 @@ -426,7 +427,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: @@ -447,8 +448,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) - img.setData(mask=new_mask) + img.setData(mask=self._mask | other._mask) else: img.setData(mask=self._mask) return img @@ -512,8 +512,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) - img.setData(mask=new_mask) + img.setData(mask=self._mask | other._mask) else: img.setData(mask=self._mask) return img @@ -573,8 +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) - img.setData(mask=new_mask) + img.setData(mask=self._mask | other._mask) else: img.setData(mask=self._mask) return img @@ -652,8 +650,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) - img.setData(mask=new_mask) + img.setData(mask=self._mask | other._mask) else: img.setData(mask=self._mask) return img @@ -796,7 +793,7 @@ def measure_fiber_shifts(self, ref_image, trace_cent, columns=[500, 1000, 1500, 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 @@ -950,7 +947,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 @@ -960,6 +957,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): @@ -1094,6 +1093,78 @@ 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 + """ + 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 + 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 @@ -1340,7 +1411,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": @@ -1352,7 +1423,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: @@ -1408,7 +1479,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: @@ -1421,10 +1492,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 @@ -1511,19 +1582,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 @@ -1650,7 +1724,7 @@ def subsampleImg(self): else: new_error = None if self._mask is not None: - new_mask = numpy.zeros(new_data.shape, dtype="bool") + new_mask = numpy.zeros(new_data.shape, dtype=int) else: new_mask = None @@ -1705,18 +1779,18 @@ 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) ), 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 @@ -1762,7 +1836,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. @@ -1783,17 +1857,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 @@ -1900,7 +1975,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 @@ -2009,7 +2084,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: @@ -2571,7 +2646,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] @@ -2599,10 +2674,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): @@ -2612,7 +2687,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) @@ -2623,19 +2698,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 @@ -2646,13 +2712,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"): @@ -2660,7 +2726,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: @@ -2691,7 +2757,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( @@ -2704,7 +2770,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) @@ -3057,8 +3123,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 @@ -3081,7 +3147,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) @@ -3097,7 +3163,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))) @@ -3109,14 +3175,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: @@ -3128,10 +3194,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/core/rss.py b/python/lvmdrp/core/rss.py index 2f2d2634..46667442 100644 --- a/python/lvmdrp/core/rss.py +++ b/python/lvmdrp/core/rss.py @@ -13,6 +13,7 @@ from astropy import units as u from lvmdrp import log +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, CON_LAMPS, ARC_LAMPS from lvmdrp.core.apertures import Aperture from lvmdrp.core.cube import Cube @@ -130,7 +131,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": @@ -419,7 +420,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 = (bn.nansum(masks, axis=0)>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: @@ -449,7 +450,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 = bn.nansum(masks, axis=0).astype("bool") + new_mask = numpy.bitwise_or.reduce(masks, axis=0) if skies.size != 0: new_sky = bn.nansum(skies, axis=0) else: @@ -529,7 +530,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)) @@ -695,7 +696,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: @@ -738,7 +739,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 @@ -900,6 +901,70 @@ 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 + 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. @@ -1080,10 +1145,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, kind=PixMask) 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 @@ -1364,7 +1430,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 @@ -1423,7 +1489,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 @@ -1453,20 +1519,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_sky = None else: @@ -1486,24 +1551,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 @@ -1523,15 +1587,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) @@ -1551,9 +1613,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 @@ -1643,7 +1705,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] @@ -1672,7 +1734,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") @@ -1708,7 +1770,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 @@ -1738,7 +1800,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) @@ -1826,7 +1888,7 @@ def rectify_wave(self, wave=None, wave_range=None, wave_disp=None): 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, @@ -1851,9 +1913,9 @@ def rectify_wave(self, wave=None, wave_range=None, wave_disp=None): f = resample_flux_density(wave, rss._wave[ifiber], rss._data[ifiber]) new_rss._data[ifiber] = f.astype("float32") f = resample_flux_density(wave, rss._wave[ifiber], rss._mask[ifiber]) - new_rss._mask[ifiber] = f.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]) + new_rss._mask[ifiber] = f.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 = numpy.interp(wave, rss._wave[ifiber], rss._lsf[ifiber]) new_rss._lsf[ifiber] = f.astype("float32") @@ -1955,7 +2017,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) @@ -1997,12 +2059,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, ), ) @@ -2051,12 +2113,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, ), ) @@ -2116,7 +2178,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) @@ -2144,12 +2206,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, ), ) @@ -2238,7 +2300,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) @@ -2269,13 +2331,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, ), ) @@ -2341,7 +2403,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 ], ), @@ -2607,7 +2669,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( @@ -2615,7 +2677,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], ), ) @@ -2660,7 +2722,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] @@ -3121,26 +3183,27 @@ 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 + 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 @@ -3323,7 +3386,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: @@ -3450,7 +3513,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) @@ -3506,7 +3569,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") @@ -3532,7 +3595,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 @@ -3595,7 +3658,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 @@ -3624,7 +3687,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 @@ -3683,7 +3746,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 @@ -3710,7 +3773,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 @@ -3758,7 +3821,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/sky.py b/python/lvmdrp/core/sky.py index 031db3e4..cf165655 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 @@ -786,11 +787,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/core/spectrum1d.py b/python/lvmdrp/core/spectrum1d.py index 90e0a9f1..44a6353d 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 @@ -783,14 +784,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 @@ -914,14 +915,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 @@ -1078,14 +1079,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 @@ -1201,7 +1202,7 @@ def __truediv__(self, other): if self._mask is not None: mask = self._mask - mask[~select] = True + mask[~select] = PixMask["BADPIX"] else: mask = None @@ -1259,7 +1260,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( @@ -1320,14 +1321,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 @@ -1446,7 +1447,7 @@ def __rtruediv__(self, other): if self._mask is not None: mask = self._mask - mask[~select] = True + mask[~select] = PixMask["BADPIX"] else: mask = None @@ -1512,7 +1513,7 @@ def __rtruediv__(self, other): if self._mask is not None: mask = self._mask - mask[~select] = True + mask[~select] = PixMask["BADPIX"] else: mask = None @@ -1973,7 +1974,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": @@ -1991,7 +1992,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: @@ -2081,7 +2082,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: @@ -2106,10 +2107,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 @@ -2264,7 +2265,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 @@ -2295,10 +2296,10 @@ def resampleSpec( # case where input spectrum has more than half the pixels masked if bn.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) @@ -2337,8 +2338,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) @@ -2456,11 +2457,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 @@ -2471,17 +2472,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 @@ -2632,7 +2633,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 @@ -2653,7 +2656,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: @@ -2661,7 +2664,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]) @@ -2772,12 +2775,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] @@ -2820,7 +2823,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]) @@ -2852,9 +2855,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) @@ -2862,7 +2865,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") @@ -2937,7 +2940,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: @@ -2968,7 +2971,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 @@ -3020,7 +3023,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) @@ -3543,6 +3546,7 @@ def fitSepGauss( # mask |= (~numpy.isfinite(data) | ~numpy.isfinite(error)) # reset bad pixels in data and error + mask = mask != 0 error[mask] = numpy.inf data[mask] = 0.0 @@ -3605,7 +3609,7 @@ def fitSepGauss( if axs is not None: axs[i].axvspan(x[0], x[-1], alpha=0.1, fc="0.5", label="reg. of masking") - 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) @@ -3644,21 +3648,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, :]], 1) == aperture + bad_pix[~select] |= numpy.bitwise_or.reduce(self._mask[pixels[~select, :]], 1) else: bad_pix = None if self._error is None: @@ -3698,10 +3701,8 @@ 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('B.fits', B.toarray(), overwrite=True) + # pyfits.writeto('B1.fits', B1.toarray(), overwrite=True) # if plot: # plt.plot(self._data, "ok") # plt.plot(numpy.dot(A * self._error[:, None], out[0]), "-r") @@ -3720,7 +3721,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'") @@ -3764,7 +3765,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): @@ -3806,11 +3807,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/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index 3aa8e9c2..53f145c3 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -29,7 +29,7 @@ from lvmdrp import log, DRP_COMMIT, __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 PixMask from lvmdrp.core.fiberrows import FiberRows, _read_fiber_ypix from lvmdrp.core.image import ( Image, @@ -99,7 +99,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) -> Tuple[Image]: """calculates non-linearity correction for input quadrant Parameters @@ -135,6 +135,7 @@ 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)") + quadrant.update_drpstage("LINEARITY_CORRECTED") else: log.warning("cannot apply non-linearity correction") log.info(f"using {nominal_gain = } (e-/ADU)") @@ -257,13 +258,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]) @@ -299,13 +300,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): @@ -388,6 +389,9 @@ 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, trace_cent, columns=columns, column_width=column_width, shift_range=shift_range, axs=axs) + if (column_shifts!=0).any(): + image.update_drpstage("FIBERS_SHIFTED") + # shifts stats median_shift = numpy.nan_to_num(bn.nanmedian(column_shifts, axis=0)) std_shift = numpy.nan_to_num(bn.nanstd(column_shifts, axis=0)) @@ -403,7 +407,7 @@ def _fix_fiber_thermal_shifts(image, trace_cent, trace_width=None, trace_amp=Non trace_cent_fixed._coeffs[:, 0] += median_shift trace_cent_fixed.eval_coeffs() - return trace_cent_fixed, column_shifts, median_shift, std_shift, fiber_model + return image, trace_cent_fixed, column_shifts, median_shift, std_shift, fiber_model def _apply_electronic_shifts(images, out_images, drp_shifts=None, qc_shifts=None, custom_shifts=None, raw_shifts=None, @@ -594,7 +598,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") @@ -859,7 +863,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: @@ -1617,12 +1621,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) @@ -1632,7 +1636,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 @@ -1927,8 +1931,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") @@ -1947,11 +1951,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") @@ -1972,6 +1976,7 @@ def subtract_straylight( # 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.update_drpstage("STRAYLIGHT_SUBTRACTED") img_out.writeFitsData(out_image) # plot results: polyomial fitting & smoothing, both with masked regions on @@ -2119,7 +2124,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") @@ -2149,7 +2154,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() @@ -2198,13 +2203,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 + proc_img.add_bitmask("SATURATED", where=saturated_mask) # NOTE: this is a patch to mask first/last 3 columns as they don't have any useful data - proc_img._mask[:, :3] |= True - proc_img._mask[:, -3:] |= True + proc_img._mask[:, :3] |= PixMask["NODATA"] + proc_img._mask[:, -3:] |= PixMask["NODATA"] # 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 @@ -3342,13 +3353,8 @@ 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.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') + proc_img.update_drpqual("SATURATED") # set drp commit SHA proc_img.setHdrValue("COMMIT", DRP_COMMIT, comment="data reduction pipeline commit HASH") # add calibrations used to header @@ -3356,6 +3362,8 @@ def preproc_raw_frame( # 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.writeFitsData(out_image) # plot overscan strips along X and Y axes @@ -3579,6 +3587,7 @@ 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 astrometry: commanded position') + org_img.update_drpstage(f"{tel.upper()}_ASTROMETRY_ADDED") else: RAobs=0 DECobs=0 @@ -3626,11 +3635,11 @@ def getfibradec(tel, platescale): slitmap['ra']=RAfib * u.deg slitmap['dec']=DECfib * u.deg org_img._slitmap=slitmap - # set header keyword with best knowledge of IFU center org_img.setHdrValue('IFUCENRA', RAobs_sci, 'best SCI IFU RA (ASTRMSRC)') org_img.setHdrValue('IFUCENDE', DECobs_sci, 'best SCI IFU DEC (ASTRMSRC)') org_img.setHdrValue('IFUCENPA', PAobs_sci, 'best SCI IFU PA (ASTRMSRC)') + org_img.update_drpstage("FIBERS_ASTROMETRY_ADDED") log.info(f"writing RA,DEC to slitmap in image '{os.path.basename(out_image)}'") org_img.writeFitsData(out_image) @@ -3650,7 +3659,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 @@ -3677,8 +3685,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 """ @@ -3769,13 +3775,15 @@ def detrend_frame( # gain-correct quadrant quad *= gain_map # propagate new NaNs to the mask - quad._mask |= numpy.isnan(quad._data) + quad.add_bitmask("BADPIX", where=numpy.isnan(quad._data)) quad.computePoissonError(rdnoise) bcorr_img.setSection(section=quad_sec, subimg=quad, inplace=True) log.info(f"median error in quadrant {i+1}: {bn.nanmedian(quad._error):.2f} (e-)") bcorr_img.setHdrValue("BUNIT", "electron", "physical units of the array values") + bcorr_img.update_drpstage("GAIN_CORRECTED") + bcorr_img.update_drpstage("POISSON_ERROR_CALCULATED") else: # convert to ADU log.info("leaving original ADU units") @@ -3790,13 +3798,12 @@ 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) + detrended_img.update_drpstage("DETRENDED") # 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.add_bitmask("BADPIX", where=nanpixels) # reject cosmic rays if reject_cr: @@ -3804,17 +3811,18 @@ def detrend_frame( rdnoise = detrended_img.getHdrValue(f"{camera.upper()} 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) + detrended_img.update_drpstage("COSMIC_CLEANED") # 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 # '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) @@ -4159,7 +4167,7 @@ def create_pixelmask(in_short_dark, in_long_dark, out_pixmask, in_flat_a=None, i sections = short_dark.getHdrValue(f"{camera.upper()} 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) @@ -4258,11 +4266,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)}'") @@ -4316,7 +4325,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: @@ -4337,7 +4346,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() @@ -4346,7 +4355,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") @@ -4512,7 +4521,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)}'") @@ -4557,9 +4566,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: @@ -4574,9 +4583,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") @@ -4634,7 +4643,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 fdeb8a0b..2ecd8b17 100644 --- a/python/lvmdrp/functions/rssMethod.py +++ b/python/lvmdrp/functions/rssMethod.py @@ -22,6 +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 PixMask from lvmdrp.core.constants import CONFIG_PATH, ARC_LAMPS from lvmdrp.core.cube import Cube from lvmdrp.core.tracemask import TraceMask @@ -94,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]) @@ -382,7 +383,7 @@ def determine_wavelength_solution(in_arcs: List[str]|str, out_wave: str, out_lsf log.info(f"found {len(exposed)} exposed standard fibers: {fibermap['orig_ifulabel'][exposed]}") # 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 @@ -752,7 +753,7 @@ def shift_wave_skylines(in_frame: str, out_frame: str, dwave: float = 8.0, skyli for ifiber in iterator: # skip dead/non-exposed fibers 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 # skip fibers with low S/N @@ -807,6 +808,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.update_drpstage("WAVELENGTH_SHIFTED") lvmframe.writeFitsData(out_frame) # Make QA plots showing offsets for each sky line in each spectrograph @@ -897,6 +899,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.update_drpstage("WAVELENGTH_CALIBRATED") rss.writeFitsData(out_rss) return rss @@ -1394,7 +1397,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: @@ -1420,19 +1423,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}%)") @@ -1484,7 +1487,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]) @@ -1693,6 +1696,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)) + lvmframe.update_drpstage("FLATFIELDED") lvmframe.writeFitsData(out_frame) return rss, lvmframe @@ -2511,7 +2515,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 @@ -2530,7 +2534,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) @@ -2991,6 +2995,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.update_drpstage("SPECTROGRAPH_STACKED") rss_out.writeFitsData(out_rss) return rss_out @@ -3034,6 +3039,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.update_drpstage("CHANNEL_COMBINED") + # write output RSS if out_cframe is not None: log.info(f"writing output RSS to {os.path.basename(out_cframe)}") @@ -3112,7 +3119,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]) diff --git a/python/lvmdrp/functions/skyMethod.py b/python/lvmdrp/functions/skyMethod.py index 374cc634..dc555036 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, 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 +1397,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 @@ -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 @@ -1581,6 +1583,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 @@ -1653,9 +1656,9 @@ def quick_sky_subtraction(in_cframe, out_sframe, skymethod: str = 'farlines_near # 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 @@ -1921,22 +1924,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 1325ef2d..02ded2ca 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): @@ -98,6 +99,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 BAD = auto() # bit whether a raw frame is bad for reduction TEST = auto() # bit whether a raw frame is for instrument testing purposes @@ -141,25 +143,76 @@ 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 + PREPROCESSED = auto() + 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. + # 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. @@ -169,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. @@ -181,68 +250,130 @@ 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): - # fiber bitmasks - NOPLUG = auto() - BADTRACE = auto() - BADFLAT = auto() - BADARC = auto() - MANYBADCOLUMNS = auto() - MANYREJECTED = auto() - LARGESHIFT = auto() - BADSKYFIBER = auto() - NEARWHOPPER = auto() - WHOPPER = auto() - SMEARIMAGE = auto() - SMEARHIGHSN = auto() - SMEARMEDSN = auto() - DEADFIBER = auto() - # pixel bitmasks - SATURATION = auto() + # pixel bitmasks ------------------------------------------------------ + # from pixelmasks BADPIX = auto() - COSMIC = auto() NEARBADPIXEL = auto() + + # from raw data preprocessing + SATURATED = 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() - BADSKYCHI = auto() + + # large sky residuals + BADSKYFIT = 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 -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()) +DRPQUALITIES = list() + + +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 {bitmask = }: {type(bitmask)}; expected {kind}, string or integer") + + return bitmask + + +def _parse_where(where, mask_shape): + 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) + + return where + + +def add_bitmask(mask_image, pixmask, where=None): + pixmask = _parse_bitmask(pixmask, kind=PixMask) + where = _parse_where(where, mask_shape=mask_image.shape) + + mask_image[where] |= pixmask + return mask_image + + +def toggle_bitmask(mask_image, pixmask, where=None): + 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 + -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()) +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)