22 """Extract small cutouts around bright stars, normalize and warp them to the
23 same arbitrary pixel grid.
26 __all__ = [
"ProcessBrightStarsTask"]
29 import astropy.units
as u
33 from lsst.afw import image
as afwImage
34 from lsst.afw import detection
as afwDetect
35 from lsst.afw import cameraGeom
as cg
44 from lsst.utils.timer
import timeMethod
48 dimensions=(
"instrument",
"visit",
"detector")):
49 inputExposure = cT.Input(
50 doc=
"Input exposure from which to extract bright star stamps",
52 storageClass=
"ExposureF",
53 dimensions=(
"visit",
"detector")
56 doc=
"Input Sky Correction to be subtracted from the calexp if doApplySkyCorr=True",
58 storageClass=
"Background",
59 dimensions=(
"instrument",
"visit",
"detector")
61 refCat = cT.PrerequisiteInput(
62 doc=
"Reference catalog that contains bright star positions",
63 name=
"gaia_dr2_20200414",
64 storageClass=
"SimpleCatalog",
65 dimensions=(
"skypix",),
69 brightStarStamps = cT.Output(
70 doc=
"Set of preprocessed postage stamps, each centered on a single bright star.",
71 name=
"brightStarStamps",
72 storageClass=
"BrightStarStamps",
73 dimensions=(
"visit",
"detector")
78 if not config.doApplySkyCorr:
79 self.inputs.remove(
"skyCorr")
83 pipelineConnections=ProcessBrightStarsConnections):
84 """Configuration parameters for ProcessBrightStarsTask
86 magLimit = pexConfig.Field(
88 doc=
"Magnitude limit, in Gaia G; all stars brighter than this value will be processed",
91 stampSize = pexConfig.ListField(
93 doc=
"Size of the stamps to be extracted, in pixels",
96 modelStampBuffer = pexConfig.Field(
98 doc=
"'Buffer' factor to be applied to determine the size of the stamp the processed stars will "
99 "be saved in. This will also be the size of the extended PSF model.",
102 doRemoveDetected = pexConfig.Field(
104 doc=
"Whether DETECTION footprints, other than that for the central object, should be changed to "
108 doApplyTransform = pexConfig.Field(
110 doc=
"Apply transform to bright star stamps to correct for optical distortions?",
113 warpingKernelName = pexConfig.ChoiceField(
115 doc=
"Warping kernel",
118 "bilinear":
"bilinear interpolation",
119 "lanczos3":
"Lanczos kernel of order 3",
120 "lanczos4":
"Lanczos kernel of order 4",
121 "lanczos5":
"Lanczos kernel of order 5",
124 annularFluxRadii = pexConfig.ListField(
126 doc=
"Inner and outer radii of the annulus used to compute the AnnularFlux for normalization, "
130 annularFluxStatistic = pexConfig.ChoiceField(
132 doc=
"Type of statistic to use to compute annular flux.",
137 "MEANCLIP":
"clipped mean",
140 numSigmaClip = pexConfig.Field(
142 doc=
"Sigma for outlier rejection; ignored if annularFluxStatistic != 'MEANCLIP'.",
145 numIter = pexConfig.Field(
147 doc=
"Number of iterations of outlier rejection; ignored if annularFluxStatistic != 'MEANCLIP'.",
150 badMaskPlanes = pexConfig.ListField(
152 doc=
"Mask planes that, if set, lead to associated pixels not being included in the computation of the"
154 default=(
'BAD',
'CR',
'CROSSTALK',
'EDGE',
'NO_DATA',
'SAT',
'SUSPECT',
'UNMASKEDNAN')
156 minPixelsWithinFrame = pexConfig.Field(
158 doc=
"Minimum number of pixels that must fall within the stamp boundary for the bright star to be"
159 " saved when its center is beyond the exposure boundary.",
162 doApplySkyCorr = pexConfig.Field(
164 doc=
"Apply full focal plane sky correction before extracting stars?",
167 discardNanFluxStars = pexConfig.Field(
169 doc=
"Should stars with NaN annular flux be discarded?",
172 refObjLoader = pexConfig.ConfigurableField(
173 target=LoadIndexedReferenceObjectsTask,
174 doc=
"Reference object loader for astrometric calibration.",
178 self.
refObjLoaderrefObjLoader.ref_dataset_name =
"gaia_dr2_20200414"
182 """The description of the parameters for this Task are detailed in
183 :lsst-task:`~lsst.pipe.base.PipelineTask`.
187 `ProcessBrightStarsTask` is used to extract, process, and store small
188 image cut-outs (or "postage stamps") around bright stars. It relies on
189 three methods, called in succession:
192 Find bright stars within the exposure using a reference catalog and
193 extract a stamp centered on each.
195 Shift and warp each stamp to remove optical distortions and sample all
196 stars on the same pixel grid.
197 `measureAndNormalize`
198 Compute the flux of an object in an annulus and normalize it. This is
199 required to normalize each bright star stamp as their central pixels
200 are likely saturated and/or contain ghosts, and cannot be used.
202 ConfigClass = ProcessBrightStarsConfig
203 _DefaultName =
"processBrightStars"
204 RunnerClass = pipeBase.ButlerInitializedTaskRunner
206 def __init__(self, butler=None, initInputs=None, *args, **kwargs):
209 self.
modelStampSizemodelStampSize = [int(self.config.stampSize[0]*self.config.modelStampBuffer),
210 int(self.config.stampSize[1]*self.config.modelStampBuffer)]
219 if butler
is not None:
220 self.makeSubtask(
'refObjLoader', butler=butler)
223 """Apply correction to the sky background level.
225 Sky corrections can be generated with the 'skyCorrection.py'
226 executable in pipe_drivers. Because the sky model used by that
227 code extends over the entire focal plane, this can produce
228 better sky subtraction.
229 The calexp is updated in-place.
233 calexp : `lsst.afw.image.Exposure` or `lsst.afw.image.MaskedImage`
235 skyCorr : `lsst.afw.math.backgroundList.BackgroundList` or None,
237 Full focal plane sky correction, obtained by running
238 `lsst.pipe.drivers.skyCorrection.SkyCorrectionTask`.
241 calexp = calexp.getMaskedImage()
242 calexp -= skyCorr.getImage()
245 """ Read position of bright stars within `inputExposure` from refCat
250 inputExposure : `afwImage.exposure.exposure.ExposureF`
251 The image from which bright star stamps should be extracted.
252 refObjLoader : `LoadIndexedReferenceObjectsTask`, optional
253 Loader to find objects within a reference catalog.
257 result : `lsst.pipe.base.Struct`
258 Result struct with components:
260 - ``starIms``: `list` of stamps
261 - ``pixCenters``: `list` of corresponding coordinates to each
262 star's center, in pixels.
263 - ``GMags``: `list` of corresponding (Gaia) G magnitudes.
264 - ``gaiaIds``: `np.ndarray` of corresponding unique Gaia
267 if refObjLoader
is None:
268 refObjLoader = self.refObjLoader
273 wcs = inputExposure.getWcs()
275 inputIm = inputExposure.maskedImage
276 inputExpBBox = inputExposure.getBBox()
277 dilatationExtent =
geom.Extent2I(np.array(self.config.stampSize) - self.config.minPixelsWithinFrame)
279 withinCalexp = refObjLoader.loadPixelBox(inputExpBBox.dilatedBy(dilatationExtent), wcs,
280 filterName=
"phot_g_mean")
281 refCat = withinCalexp.refCat
283 fluxLimit = ((self.config.magLimit*u.ABmag).
to(u.nJy)).to_value()
284 GFluxes = np.array(refCat[
'phot_g_mean_flux'])
285 bright = GFluxes > fluxLimit
287 allGMags = [((gFlux*u.nJy).
to(u.ABmag)).to_value()
for gFlux
in GFluxes[bright]]
288 allIds = refCat.columns.extract(
"id", where=bright)[
"id"]
289 selectedColumns = refCat.columns.extract(
'coord_ra',
'coord_dec', where=bright)
290 for j, (ra, dec)
in enumerate(zip(selectedColumns[
"coord_ra"], selectedColumns[
"coord_dec"])):
292 cpix = wcs.skyToPixel(sp)
294 starIm = inputExposure.getCutout(sp,
geom.Extent2I(self.config.stampSize))
295 except InvalidParameterError:
297 bboxCorner = np.array(cpix) - np.array(self.config.stampSize)/2
301 clippedStarBBox.clip(inputExpBBox)
302 if clippedStarBBox.getArea() > 0:
305 starIm = afwImage.ExposureF(bbox=idealBBox)
306 starIm.image[:] = np.nan
307 starIm.mask.set(inputExposure.mask.getPlaneBitMask(
"NO_DATA"))
309 clippedIm = inputIm.Factory(inputIm, clippedStarBBox)
310 starIm.maskedImage[clippedStarBBox] = clippedIm
312 starIm.setDetector(inputExposure.getDetector())
313 starIm.setWcs(inputExposure.getWcs())
316 if self.config.doRemoveDetected:
319 afwDetect.Threshold.BITMASK)
321 allFootprints = omask.getFootprints()
323 for fs
in allFootprints:
325 otherFootprints.append(fs)
326 nbMatchingFootprints = len(allFootprints) - len(otherFootprints)
327 if not nbMatchingFootprints == 1:
328 self.log.
warning(
"Failed to uniquely identify central DETECTION footprint for star "
329 "%s; found %d footprints instead.",
330 allIds[j], nbMatchingFootprints)
331 omask.setFootprints(otherFootprints)
332 omask.setMask(starIm.mask,
"BAD")
333 starIms.append(starIm)
334 pixCenters.append(cpix)
335 GMags.append(allGMags[j])
336 ids.append(allIds[j])
337 return pipeBase.Struct(starIms=starIms,
338 pixCenters=pixCenters,
343 """Warps and shifts all given stamps so they are sampled on the same
344 pixel grid and centered on the central pixel. This includes rotating
345 the stamp depending on detector orientation.
349 stamps : `collections.abc.Sequence`
350 [`afwImage.exposure.exposure.ExposureF`]
351 Image cutouts centered on a single object.
352 pixCenters : `collections.abc.Sequence` [`geom.Point2D`]
353 Positions of each object's center (as obtained from the refCat),
358 result : `lsst.pipe.base.Struct`
359 Result struct with components:
362 `list` [`afwImage.maskedImage.maskedImage.MaskedImage`] of
363 stamps of warped stars
364 - ``warpTransforms``:
365 `list` [`afwGeom.TransformPoint2ToPoint2`] of
366 the corresponding Transform from the initial star stamp to
367 the common model grid
369 `list` [`geom.Point2I`] of coordinates of the bottom-left
370 pixels of each stamp, before rotation
371 - ``nb90Rots``: `int`, the number of 90 degrees rotations required
372 to compensate for detector orientation
377 bufferPix = (self.
modelStampSizemodelStampSize[0] - self.config.stampSize[0],
381 det = stamps[0].getDetector()
383 if self.config.doApplyTransform:
384 pixToTan = det.getTransform(cg.PIXELS, cg.TAN_PIXELS)
386 pixToTan = tFactory.makeIdentityTransform()
388 possibleRots = np.array([k*np.pi/2
for k
in range(4)])
390 yaw = det.getOrientation().getYaw()
391 nb90Rots = np.argmin(np.abs(possibleRots - float(yaw)))
394 warpedStars, warpTransforms, xy0s = [], [], []
395 for star, cent
in zip(stamps, pixCenters):
397 destImage = afwImage.MaskedImageF(*self.
modelStampSizemodelStampSize)
399 newBottomLeft = pixToTan.applyForward(bottomLeft)
400 newBottomLeft.setX(newBottomLeft.getX() - bufferPix[0]/2)
401 newBottomLeft.setY(newBottomLeft.getY() - bufferPix[1]/2)
405 destImage.setXY0(newBottomLeft)
406 xy0s.append(newBottomLeft)
409 newCenter = pixToTan.applyForward(cent)
410 shift = self.
modelCentermodelCenter[0] + newBottomLeft[0] - newCenter[0],\
411 self.
modelCentermodelCenter[1] + newBottomLeft[1] - newCenter[1]
413 shiftTransform = tFactory.makeTransform(affineShift)
416 starWarper = pixToTan.then(shiftTransform)
420 starWarper, warpCont)
422 self.log.
debug(
"Warping of a star failed: no good pixel in output")
425 destImage.setXY0(0, 0)
430 warpedStars.append(destImage.clone())
431 warpTransforms.append(starWarper)
432 return pipeBase.Struct(warpedStars=warpedStars, warpTransforms=warpTransforms, xy0s=xy0s,
436 def run(self, inputExposure, refObjLoader=None, dataId=None, skyCorr=None):
437 """Identify bright stars within an exposure using a reference catalog,
438 extract stamps around each, then preprocess them. The preprocessing
439 steps are: shifting, warping and potentially rotating them to the same
440 pixel grid; computing their annular flux and normalizing them.
444 inputExposure : `afwImage.exposure.exposure.ExposureF`
445 The image from which bright star stamps should be extracted.
446 refObjLoader : `LoadIndexedReferenceObjectsTask`, optional
447 Loader to find objects within a reference catalog.
448 dataId : `dict` or `lsst.daf.butler.DataCoordinate`
449 The dataId of the exposure (and detector) bright stars should be
451 skyCorr : `lsst.afw.math.backgroundList.BackgroundList` or ``None``,
453 Full focal plane sky correction, obtained by running
454 `lsst.pipe.drivers.skyCorrection.SkyCorrectionTask`.
458 result : `lsst.pipe.base.Struct`
459 Result struct with component:
461 - ``brightStarStamps``: ``bSS.BrightStarStamps``
463 if self.config.doApplySkyCorr:
464 self.log.
info(
"Applying sky correction to exposure %s (exposure will be modified in-place).",
467 self.log.
info(
"Extracting bright stars from exposure %s", dataId)
469 extractedStamps = self.
extractStampsextractStamps(inputExposure, refObjLoader=refObjLoader)
470 if not extractedStamps.starIms:
471 self.log.
info(
"No suitable bright star found.")
474 self.log.
info(
"Applying warp and/or shift to %i star stamps from exposure %s",
475 len(extractedStamps.starIms), dataId)
476 warpOutputs = self.
warpStampswarpStamps(extractedStamps.starIms, extractedStamps.pixCenters)
477 warpedStars = warpOutputs.warpedStars
478 xy0s = warpOutputs.xy0s
479 brightStarList = [bSS.BrightStarStamp(stamp_im=warp,
480 archive_element=transform,
482 gaiaGMag=extractedStamps.GMags[j],
483 gaiaId=extractedStamps.gaiaIds[j])
484 for j, (warp, transform)
in
485 enumerate(zip(warpedStars, warpOutputs.warpTransforms))]
487 self.log.
info(
"Computing annular flux and normalizing %i bright stars from exposure %s",
488 len(warpedStars), dataId)
491 statsControl.setNumSigmaClip(self.config.numSigmaClip)
492 statsControl.setNumIter(self.config.numIter)
493 innerRadius, outerRadius = self.config.annularFluxRadii
495 brightStarStamps = bSS.BrightStarStamps.initAndNormalize(brightStarList,
496 innerRadius=innerRadius,
497 outerRadius=outerRadius,
498 nb90Rots=warpOutputs.nb90Rots,
501 statsControl=statsControl,
503 badMaskPlanes=self.config.badMaskPlanes,
504 discardNanFluxObjects=(
505 self.config.discardNanFluxStars))
506 return pipeBase.Struct(brightStarStamps=brightStarStamps)
509 """Read in required calexp, extract and process stamps around bright
510 stars and write them to disk.
514 dataRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
515 Data reference to the calexp to extract bright stars from.
517 calexp = dataRef.get(
"calexp")
518 skyCorr = dataRef.get(
"skyCorr")
if self.config.doApplySkyCorr
else None
519 output = self.
runrun(calexp, dataId=dataRef.dataId, skyCorr=skyCorr)
521 dataRef.put(output.brightStarStamps,
"brightStarStamps")
522 return pipeBase.Struct(brightStarStamps=output.brightStarStamps)
525 inputs = butlerQC.get(inputRefs)
526 inputs[
'dataId'] = str(butlerQC.quantum.dataId)
528 for ref
in inputRefs.refCat],
529 refCats=inputs.pop(
"refCat"),
530 config=self.config.refObjLoader)
531 output = self.
runrun(**inputs, refObjLoader=refObjLoader)
533 butlerQC.put(output, outputRefs)
A Threshold is used to pass a threshold value to detection algorithms.
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Pass parameters to a Statistics object.
Parameters to control convolution.
An integer coordinate rectangle.
Point in an unspecified spherical coordinate system.
def __init__(self, *config=None)
def runDataRef(self, dataRef)
def runQuantum(self, butlerQC, inputRefs, outputRefs)
def warpStamps(self, stamps, pixCenters)
def run(self, inputExposure, refObjLoader=None, dataId=None, skyCorr=None)
def applySkyCorr(self, calexp, skyCorr)
def extractStamps(self, inputExposure, refObjLoader=None)
def __init__(self, butler=None, initInputs=None, *args, **kwargs)
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
std::shared_ptr< ImageT > rotateImageBy90(ImageT const &image, int nQuarter)
Rotate an image by an integral number of quarter turns.
int warpImage(DestImageT &destImage, geom::SkyWcs const &destWcs, SrcImageT const &srcImage, geom::SkyWcs const &srcWcs, WarpingControl const &control, typename DestImageT::SinglePixel padValue=lsst::afw::math::edgePixel< DestImageT >(typename lsst::afw::image::detail::image_traits< DestImageT >::image_category()))
Warp an Image or MaskedImage to a new Wcs.
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.