From fa4c7c50300f7ab5fdfdef419dbd6338a96a241d Mon Sep 17 00:00:00 2001 From: "Amir E. Bazkiaei" Date: Sun, 25 Jun 2023 16:23:37 +1000 Subject: [PATCH 1/4] DM-38703: Adding stars that are missing from bright stars stamps to Bright Star Subtracting pipeline --- .../lsst/meas/algorithms/brightStarStamps.py | 56 ++++++++++++++++++- 1 file changed, 54 insertions(+), 2 deletions(-) diff --git a/python/lsst/meas/algorithms/brightStarStamps.py b/python/lsst/meas/algorithms/brightStarStamps.py index ba4be52b2..feb740596 100644 --- a/python/lsst/meas/algorithms/brightStarStamps.py +++ b/python/lsst/meas/algorithms/brightStarStamps.py @@ -68,6 +68,8 @@ class BrightStarStamp(AbstractStamp): archive_element: Persistable | None = None annularFlux: float | None = None minValidAnnulusFraction: float = 0.0 + optimalInnerRadius: int | None = None + optimalOuterRadius: int | None = None @classmethod def factory(cls, stamp_im, metadata, idx, archive_element=None, minValidAnnulusFraction=0.0): @@ -267,6 +269,7 @@ def initAndNormalize( use_archive=False, imCenter=None, discardNanFluxObjects=True, + forceFindFlux=False, statsControl=StatisticsControl(), statsFlag=stringToStatisticsProperty("MEAN"), badMaskPlanes=("BAD", "SAT", "NO_DATA"), @@ -331,12 +334,13 @@ def initAndNormalize( stamps, stamps normalized with different annulus definitions, or if stamps are to be normalized but annular radii were not provided. """ + stampSize = starStamps[0].stamp_im.getDimensions() if imCenter is None: - stampSize = starStamps[0].stamp_im.getDimensions() imCenter = stampSize[0] // 2, stampSize[1] // 2 # Create SpanSet of annulus outerCircle = SpanSet.fromShape(outerRadius, Stencil.CIRCLE, offset=imCenter) innerCircle = SpanSet.fromShape(innerRadius, Stencil.CIRCLE, offset=imCenter) + annulusWidth = outerRadius - innerRadius annulus = outerCircle.intersectNot(innerCircle) # Initialize (unnormalized) brightStarStamps instance bss = cls( @@ -354,6 +358,7 @@ def initAndNormalize( bss._innerRadius, bss._outerRadius = innerRadius, outerRadius # Create a list to contain rejected stamps. rejects = [] + badStamps = [] # Apply normalization for stamp in bss._stamps: try: @@ -367,13 +372,60 @@ def initAndNormalize( # steps needed before bright stars can be subtracted. if discardNanFluxObjects: rejects.append(stamp) + elif forceFindFlux: + newInnerRadius = innerRadius + newOuterRadius = outerRadius + while True: + newOuterRadius += annulusWidth + newInnerRadius += annulusWidth + if newOuterRadius > min(imCenter): + logger.info(f"No flux found for the star with Gaia ID of {stamp.gaiaId}") + stamp.annularFlux = None + badStamps.append(stamp) + break + newOuterCircle = SpanSet.fromShape(newOuterRadius, Stencil.CIRCLE, offset=imCenter) + newInnerCircle = SpanSet.fromShape(newInnerRadius, Stencil.CIRCLE, offset=imCenter) + newAnnulus = newOuterCircle.intersectNot(newInnerCircle) + try: + stamp.measureAndNormalize( + newAnnulus, + statsControl=statsControl, + statsFlag=statsFlag, + badMaskPlanes=badMaskPlanes, + ) + + except RuntimeError as err: + stamp.annularFlux = np.nan + logger.error( + "The annulux flux was not found for radii {} and {}".format( + newInnerRadius, newOuterRadius + ) + ) + if stamp.annularFlux and stamp.annularFlux > 0: + logger.info("The flux is found within an optimized annulus.") + logger.info( + "The optimized annulus raddi are {} and {} and the flux is {}".format( + newInnerRadius, newOuterRadius, stamp.annularFlux + ) + ) + stamp.optimalOuterRadius = newOuterRadius + stamp.optimalInnerRadius = newInnerRadius + break + # elif newInnerRadius > maxInnerRadius: + # logger.info("The star with gaiaId of {} is impossible to normalize!".format(stamp.gaiaId)) + # break else: stamp.annularFlux = np.nan # Remove rejected stamps. + bss.normalized = True if discardNanFluxObjects: for reject in rejects: bss._stamps.remove(reject) - bss.normalized = True + elif forceFindFlux: + for badStamp in badStamps: + bss._stamps.remove(badStamp) + bss._innerRadius, bss._outerRadius = None, None + return bss, badStamps return bss def _refresh_metadata(self): From 5473c6ef4b57040b418ff8ab6b473b533489f54b Mon Sep 17 00:00:00 2001 From: "Amir E. Bazkiaei" Date: Thu, 6 Jul 2023 11:59:11 +1000 Subject: [PATCH 2/4] minor correction to pass testes --- python/lsst/meas/algorithms/brightStarStamps.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/python/lsst/meas/algorithms/brightStarStamps.py b/python/lsst/meas/algorithms/brightStarStamps.py index feb740596..4371f29a0 100644 --- a/python/lsst/meas/algorithms/brightStarStamps.py +++ b/python/lsst/meas/algorithms/brightStarStamps.py @@ -394,7 +394,7 @@ def initAndNormalize( badMaskPlanes=badMaskPlanes, ) - except RuntimeError as err: + except RuntimeError: stamp.annularFlux = np.nan logger.error( "The annulux flux was not found for radii {} and {}".format( @@ -411,9 +411,6 @@ def initAndNormalize( stamp.optimalOuterRadius = newOuterRadius stamp.optimalInnerRadius = newInnerRadius break - # elif newInnerRadius > maxInnerRadius: - # logger.info("The star with gaiaId of {} is impossible to normalize!".format(stamp.gaiaId)) - # break else: stamp.annularFlux = np.nan # Remove rejected stamps. @@ -421,7 +418,7 @@ def initAndNormalize( if discardNanFluxObjects: for reject in rejects: bss._stamps.remove(reject) - elif forceFindFlux: + elif forceFindFlux: for badStamp in badStamps: bss._stamps.remove(badStamp) bss._innerRadius, bss._outerRadius = None, None From 5404897a9c8db350c99c4c664ba21fb37ba82758 Mon Sep 17 00:00:00 2001 From: "Amir E. Bazkiaei" Date: Thu, 20 Jul 2023 16:24:34 +1000 Subject: [PATCH 3/4] some typos corrected --- python/lsst/meas/algorithms/brightStarStamps.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/lsst/meas/algorithms/brightStarStamps.py b/python/lsst/meas/algorithms/brightStarStamps.py index 4371f29a0..f114934f9 100644 --- a/python/lsst/meas/algorithms/brightStarStamps.py +++ b/python/lsst/meas/algorithms/brightStarStamps.py @@ -397,14 +397,14 @@ def initAndNormalize( except RuntimeError: stamp.annularFlux = np.nan logger.error( - "The annulux flux was not found for radii {} and {}".format( + "The annular flux was not found for radii {} and {}".format( newInnerRadius, newOuterRadius ) ) if stamp.annularFlux and stamp.annularFlux > 0: logger.info("The flux is found within an optimized annulus.") logger.info( - "The optimized annulus raddi are {} and {} and the flux is {}".format( + "The optimized annulus radii are {} and {} and the flux is {}".format( newInnerRadius, newOuterRadius, stamp.annularFlux ) ) From 0dcb3da40c4346e5bcdb91716919d0da0f29027d Mon Sep 17 00:00:00 2001 From: "Amir E. Bazkiaei" Date: Mon, 4 Sep 2023 22:50:26 +0800 Subject: [PATCH 4/4] some modifications that are due to some changes in pipe_tasks are applied' --- .../lsst/meas/algorithms/brightStarStamps.py | 126 ++++++++++-------- 1 file changed, 72 insertions(+), 54 deletions(-) diff --git a/python/lsst/meas/algorithms/brightStarStamps.py b/python/lsst/meas/algorithms/brightStarStamps.py index f114934f9..2b889f3aa 100644 --- a/python/lsst/meas/algorithms/brightStarStamps.py +++ b/python/lsst/meas/algorithms/brightStarStamps.py @@ -145,27 +145,26 @@ def measureAndNormalize( Collection of mask planes to ignore when computing annularFlux. """ stampSize = self.stamp_im.getDimensions() - # create image: same pixel values within annulus, NO_DATA elsewhere + # Create image: science pixel values within annulus, NO_DATA elsewhere maskPlaneDict = self.stamp_im.mask.getMaskPlaneDict() annulusImage = MaskedImageF(stampSize, planeDict=maskPlaneDict) annulusMask = annulusImage.mask annulusMask.array[:] = 2 ** maskPlaneDict["NO_DATA"] annulus.copyMaskedImage(self.stamp_im, annulusImage) - # set mask planes to be ignored + # Set mask planes to be ignored. andMask = reduce(ior, (annulusMask.getPlaneBitMask(bm) for bm in badMaskPlanes)) statsControl.setAndMask(andMask) annulusStat = makeStatistics(annulusImage, statsFlag, statsControl) - # determining the number of valid (unmasked) pixels within the annulus + # Determine the number of valid (unmasked) pixels within the annulus. unMasked = annulusMask.array.size - np.count_nonzero(annulusMask.array) - # should we have some messages printed here? logger.info( "The Star's annulus contains %s valid pixels and the annulus itself contains %s pixels.", unMasked, annulus.getArea(), ) if unMasked > (annulus.getArea() * self.minValidAnnulusFraction): - # compute annularFlux + # Compute annularFlux. self.annularFlux = annulusStat.getValue() logger.info("Annular flux is: %s", self.annularFlux) else: @@ -176,7 +175,7 @@ def measureAndNormalize( raise RuntimeError("Annular flux computation failed, likely because there are no valid pixels.") if self.annularFlux < 0: raise RuntimeError("The annular flux is negative. The stamp can not be normalized!") - # normalize stamps + # Normalize stamps. self.stamp_im.image.array /= self.annularFlux return None @@ -246,8 +245,7 @@ def __init__( use_archive=False, ): super().__init__(starStamps, metadata, use_mask, use_variance, use_archive) - # Ensure stamps contain a flux measurement if and only if they are - # already expected to be normalized + # Ensure stamps contain a flux measure if expected to be normalized. self._checkNormalization(False, innerRadius, outerRadius) self._innerRadius, self._outerRadius = innerRadius, outerRadius if innerRadius is not None and outerRadius is not None: @@ -337,12 +335,14 @@ def initAndNormalize( stampSize = starStamps[0].stamp_im.getDimensions() if imCenter is None: imCenter = stampSize[0] // 2, stampSize[1] // 2 - # Create SpanSet of annulus + + # Create SpanSet of annulus. outerCircle = SpanSet.fromShape(outerRadius, Stencil.CIRCLE, offset=imCenter) innerCircle = SpanSet.fromShape(innerRadius, Stencil.CIRCLE, offset=imCenter) annulusWidth = outerRadius - innerRadius annulus = outerCircle.intersectNot(innerCircle) - # Initialize (unnormalized) brightStarStamps instance + + # Initialize (unnormalized) brightStarStamps instance. bss = cls( starStamps, innerRadius=None, @@ -353,18 +353,31 @@ def initAndNormalize( use_variance=use_variance, use_archive=use_archive, ) - # Ensure no stamps had already been normalized. + + # Ensure that no stamps have already been normalized. bss._checkNormalization(True, innerRadius, outerRadius) bss._innerRadius, bss._outerRadius = innerRadius, outerRadius - # Create a list to contain rejected stamps. + + # Apply normalization. rejects = [] badStamps = [] - # Apply normalization for stamp in bss._stamps: try: stamp.measureAndNormalize( annulus, statsControl=statsControl, statsFlag=statsFlag, badMaskPlanes=badMaskPlanes ) + # Stars that are missing from the first steps results might + # still have a flux within the normalization annulus. The + # following two lines make sure that these stars are included + # in the subtraction process. Failing to assign the optimal + # radii values may result in an error in `createAnnulus` + # method of `SubtractBrightStarsTask` class. An alternative to + # handle this, is to create to types of stams that are missing + # from the input brightStarStamps. One for those that have + # flux within the normalization annulus and another for those + # that do not have a flux within the normalization annulus. + stamp.optimalOuterRadius = outerRadius + stamp.optimalInnerRadius = innerRadius except RuntimeError as err: logger.error(err) # Optionally keep NaN flux objects, for bookkeeping purposes, @@ -379,7 +392,7 @@ def initAndNormalize( newOuterRadius += annulusWidth newInnerRadius += annulusWidth if newOuterRadius > min(imCenter): - logger.info(f"No flux found for the star with Gaia ID of {stamp.gaiaId}") + logger.info("No flux found for the star with Gaia ID of %s", stamp.gaiaId) stamp.annularFlux = None badStamps.append(stamp) break @@ -397,22 +410,24 @@ def initAndNormalize( except RuntimeError: stamp.annularFlux = np.nan logger.error( - "The annular flux was not found for radii {} and {}".format( - newInnerRadius, newOuterRadius - ) + "The annular flux was not found for radii %d and %d", + newInnerRadius, + newOuterRadius, ) if stamp.annularFlux and stamp.annularFlux > 0: logger.info("The flux is found within an optimized annulus.") logger.info( - "The optimized annulus radii are {} and {} and the flux is {}".format( - newInnerRadius, newOuterRadius, stamp.annularFlux - ) + "The optimized annulus radii are %d and %d and the flux is %f", + newInnerRadius, + newOuterRadius, + stamp.annularFlux, ) stamp.optimalOuterRadius = newOuterRadius stamp.optimalInnerRadius = newInnerRadius break else: stamp.annularFlux = np.nan + # Remove rejected stamps. bss.normalized = True if discardNanFluxObjects: @@ -426,9 +441,11 @@ def initAndNormalize( return bss def _refresh_metadata(self): - """Refresh metadata. Should be called before writing the object out.""" - # add full list of positions, Gaia magnitudes, IDs and annularFlxes to - # shared metadata + """Refresh metadata. Should be called before writing the object out. + + This method adds full lists of positions, Gaia magnitudes, IDs and + annular fluxes to the shared metadata. + """ self._metadata["G_MAGS"] = self.getMagnitudes() self._metadata["GAIA_IDS"] = self.getGaiaIds() positions = self.getPositions() @@ -449,7 +466,7 @@ def readFits(cls, filename): Parameters ---------- filename : `str` - Name of the file to read + Name of the file to read. """ return cls.readFitsWithOptions(filename, None) @@ -460,9 +477,9 @@ def readFitsWithOptions(cls, filename, options): Parameters ---------- filename : `str` - Name of the file to read + Name of the file to read. options : `PropertyList` - Collection of metadata parameters + Collection of metadata parameters. """ stamps, metadata = readFitsWithOptions(filename, BrightStarStamp.factory, options) nb90Rots = metadata["NB_90_ROTS"] if "NB_90_ROTS" in metadata else None @@ -530,11 +547,12 @@ def extend(self, bss): self._stamps += bss._stamps def getMagnitudes(self): - """Retrieve Gaia G magnitudes for each star. + """Retrieve Gaia G-band magnitudes for each star. Returns ------- gaiaGMags : `list` [`float`] + Gaia G-band magnitudes for each star. """ return [stamp.gaiaGMag for stamp in self._stamps] @@ -544,20 +562,23 @@ def getGaiaIds(self): Returns ------- gaiaIds : `list` [`int`] + Gaia IDs for each star. """ return [stamp.gaiaId for stamp in self._stamps] def getAnnularFluxes(self): - """Retrieve normalization factors for each star. + """Retrieve normalization factor for each star. These are computed by integrating the flux in annulus centered on the bright star, far enough from center to be beyond most severe ghosts and - saturation. The inner and outer radii that define the annulus can be - recovered from the metadata. + saturation. + The inner and outer radii that define the annulus can be recovered from + the metadata. Returns ------- annularFluxes : `list` [`float`] + Annular fluxes which give the normalization factor for each star. """ return [stamp.annularFlux for stamp in self._stamps] @@ -577,8 +598,7 @@ def selectByMag(self, magMin=None, magMax=None): for stamp in self._stamps if (magMin is None or stamp.gaiaGMag > magMin) and (magMax is None or stamp.gaiaGMag < magMax) ] - # This is an optimization to save looping over the init argument when - # it is already guaranteed to be the correct type + # This saves looping over init when guaranteed to be the correct type. instance = BrightStarStamps( (), innerRadius=self._innerRadius, outerRadius=self._outerRadius, metadata=self._metadata ) @@ -591,8 +611,8 @@ def _checkRadius(self, innerRadius, outerRadius): """ if innerRadius != self._innerRadius or outerRadius != self._outerRadius: raise AttributeError( - "Trying to mix stamps normalized with annulus radii " - f"{innerRadius, outerRadius} with those of BrightStarStamp instance\n" + f"Trying to mix stamps normalized with annulus radii {innerRadius, outerRadius} with those " + "of BrightStarStamp instance\n" f"(computed with annular radii {self._innerRadius, self._outerRadius})." ) @@ -604,42 +624,40 @@ def _checkNormalization(self, normalize, innerRadius, outerRadius): nStamps = len(self) nFluxVals = nStamps - noneFluxCount if noneFluxCount and noneFluxCount < nStamps: - # at least one stamp contains an annularFlux value (i.e. has been - # normalized), but not all of them do + # At least one stamp contains an annularFlux value (i.e. has been + # normalized), but not all of them do. raise AttributeError( - f"Only {nFluxVals} stamps contain an annularFlux value.\nAll stamps in a " - "BrightStarStamps instance must either be normalized with the same annulus " - "definition, or none of them can contain an annularFlux value." + f"Only {nFluxVals} stamps contain an annularFlux value.\nAll stamps in a BrightStarStamps " + "instance must either be normalized with the same annulus definition, or none of them can " + "contain an annularFlux value." ) elif normalize: - # stamps are to be normalized; ensure annular radii are specified - # and they have no annularFlux + # Stamps are to be normalized; ensure annular radii are specified + # and they have no annularFlux. if innerRadius is None or outerRadius is None: raise AttributeError( - "For stamps to be normalized (normalize=True), please provide a valid " - "value (in pixels) for both innerRadius and outerRadius." + "For stamps to be normalized (normalize=True), please provide a valid value (in pixels) " + "for both innerRadius and outerRadius." ) elif noneFluxCount < nStamps: raise AttributeError( - f"{nFluxVals} stamps already contain an annularFlux value. For stamps to " - "be normalized, all their annularFlux must be None." + f"{nFluxVals} stamps already contain an annularFlux value. For stamps to be normalized, " + "all their annularFlux must be None." ) elif innerRadius is not None and outerRadius is not None: - # Radii provided, but normalize=False; check that stamps - # already contain annularFluxes + # Radii provided, but normalize=False; check that stamps already + # contain annularFluxes. if noneFluxCount: raise AttributeError( - f"{noneFluxCount} stamps contain no annularFlux, but annular radius " - "values were provided and normalize=False.\nTo normalize stamps, set " - "normalize to True." + f"{noneFluxCount} stamps contain no annularFlux, but annular radius values were provided " + "and normalize=False.\nTo normalize stamps, set normalize to True." ) else: # At least one radius value is missing; ensure no stamps have - # already been normalized + # already been normalized. if nFluxVals: raise AttributeError( - f"{nFluxVals} stamps contain an annularFlux value. If stamps have " - "been normalized, the innerRadius and outerRadius values used must " - "be provided." + f"{nFluxVals} stamps contain an annularFlux value. If stamps have been normalized, the " + "innerRadius and outerRadius values used must be provided." ) return None