27import lsst.afw.image
as afwImage
33 populate_sattle_visit_cache)
34from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask
36from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig
39import lsst.meas.extensions.shapeHSM
44from lsst.utils.timer
import timeMethod
46from .
import DipoleFitTask
48__all__ = [
"DetectAndMeasureConfig",
"DetectAndMeasureTask",
49 "DetectAndMeasureScoreConfig",
"DetectAndMeasureScoreTask"]
53 """Raised when the residuals in footprints of stars used to compute the
54 psf-matching kernel exceeds the configured maximum.
57 msg = (
"The ratio of residual power in source footprints on the"
58 " difference image to the power in the footprints on the"
59 f
" science image was {ratio}, which exceeds the maximum"
60 f
" threshold of {threshold}")
67 return {
"ratio": self.
ratio,
73 """Raised when there are no diaSources detected on an image difference.
76 msg = (
"No diaSources detected!")
85 dimensions=(
"instrument",
"visit",
"detector"),
86 defaultTemplates={
"coaddName":
"deep",
89 science = pipeBase.connectionTypes.Input(
90 doc=
"Input science exposure.",
91 dimensions=(
"instrument",
"visit",
"detector"),
92 storageClass=
"ExposureF",
93 name=
"{fakesType}calexp"
95 matchedTemplate = pipeBase.connectionTypes.Input(
96 doc=
"Warped and PSF-matched template used to create the difference image.",
97 dimensions=(
"instrument",
"visit",
"detector"),
98 storageClass=
"ExposureF",
99 name=
"{fakesType}{coaddName}Diff_matchedExp",
101 difference = pipeBase.connectionTypes.Input(
102 doc=
"Result of subtracting template from science.",
103 dimensions=(
"instrument",
"visit",
"detector"),
104 storageClass=
"ExposureF",
105 name=
"{fakesType}{coaddName}Diff_differenceTempExp",
107 kernelSources = pipeBase.connectionTypes.Input(
108 doc=
"Final selection of sources used for psf matching.",
109 dimensions=(
"instrument",
"visit",
"detector"),
110 storageClass=
"SourceCatalog",
111 name=
"{fakesType}{coaddName}Diff_psfMatchSources"
113 outputSchema = pipeBase.connectionTypes.InitOutput(
114 doc=
"Schema (as an example catalog) for output DIASource catalog.",
115 storageClass=
"SourceCatalog",
116 name=
"{fakesType}{coaddName}Diff_diaSrc_schema",
118 diaSources = pipeBase.connectionTypes.Output(
119 doc=
"Detected diaSources on the difference image.",
120 dimensions=(
"instrument",
"visit",
"detector"),
121 storageClass=
"SourceCatalog",
122 name=
"{fakesType}{coaddName}Diff_diaSrc",
124 subtractedMeasuredExposure = pipeBase.connectionTypes.Output(
125 doc=
"Difference image with detection mask plane filled in.",
126 dimensions=(
"instrument",
"visit",
"detector"),
127 storageClass=
"ExposureF",
128 name=
"{fakesType}{coaddName}Diff_differenceExp",
130 differenceBackground = pipeBase.connectionTypes.Output(
131 doc=
"Background model that was subtracted from the difference image.",
132 dimensions=(
"instrument",
"visit",
"detector"),
133 storageClass=
"Background",
134 name=
"difference_background",
136 maskedStreaks = pipeBase.connectionTypes.Output(
137 doc=
'Catalog of streak fit parameters for the difference image.',
138 storageClass=
"ArrowNumpyDict",
139 dimensions=(
"instrument",
"visit",
"detector"),
140 name=
"{fakesType}{coaddName}Diff_streaks",
142 glintTrailInfo = pipeBase.connectionTypes.Output(
143 doc=
'Dict of fit parameters for glint trails in the catalog.',
144 storageClass=
"ArrowNumpyDict",
145 dimensions=(
"instrument",
"visit",
"detector"),
146 name=
"trailed_glints",
149 def __init__(self, *, config):
150 super().__init__(config=config)
151 if not (self.config.writeStreakInfo
and self.config.doMaskStreaks):
152 self.outputs.remove(
"maskedStreaks")
153 if not (self.config.doSubtractBackground
and self.config.doWriteBackground):
154 self.outputs.remove(
"differenceBackground")
155 if not (self.config.writeGlintInfo):
156 self.outputs.remove(
"glintTrailInfo")
159class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
160 pipelineConnections=DetectAndMeasureConnections):
161 """Config for DetectAndMeasureTask
163 doMerge = pexConfig.Field(
166 doc=
"Merge positive and negative diaSources with grow radius "
167 "set by growFootprint"
169 doForcedMeasurement = pexConfig.Field(
172 doc=
"Force photometer diaSource locations on PVI?"
174 doAddMetrics = pexConfig.Field(
177 doc=
"Add columns to the source table to hold analysis metrics?"
179 doSubtractBackground = pexConfig.Field(
181 doc=
"Subtract a background model from the image before detection?",
184 doWriteBackground = pexConfig.Field(
186 doc=
"Persist the fitted background model?",
189 subtractInitialBackground = pexConfig.ConfigurableField(
191 doc=
"Task to perform intial background subtraction, before first detection pass.",
193 subtractFinalBackground = pexConfig.ConfigurableField(
195 doc=
"Task to perform final background subtraction, after first detection pass.",
197 detection = pexConfig.ConfigurableField(
198 target=SourceDetectionTask,
199 doc=
"Final source detection for diaSource measurement",
201 streakDetection = pexConfig.ConfigurableField(
202 target=SourceDetectionTask,
203 doc=
"Separate source detection used only for streak masking",
205 doDeblend = pexConfig.Field(
208 doc=
"Deblend DIASources after detection?"
210 deblend = pexConfig.ConfigurableField(
212 doc=
"Task to split blended sources into their components."
214 measurement = pexConfig.ConfigurableField(
215 target=DipoleFitTask,
216 doc=
"Task to measure sources on the difference image.",
221 doc=
"Run subtask to apply aperture corrections"
224 target=ApplyApCorrTask,
225 doc=
"Task to apply aperture corrections"
227 forcedMeasurement = pexConfig.ConfigurableField(
228 target=ForcedMeasurementTask,
229 doc=
"Task to force photometer science image at diaSource locations.",
231 growFootprint = pexConfig.Field(
234 doc=
"Grow positive and negative footprints by this many pixels before merging"
236 diaSourceMatchRadius = pexConfig.Field(
239 doc=
"Match radius (in arcseconds) for DiaSource to Source association"
241 doSkySources = pexConfig.Field(
244 doc=
"Generate sky sources?",
246 skySources = pexConfig.ConfigurableField(
247 target=SkyObjectsTask,
248 doc=
"Generate sky sources",
250 doMaskStreaks = pexConfig.Field(
253 doc=
"Turn on streak masking",
255 maskStreaks = pexConfig.ConfigurableField(
256 target=MaskStreaksTask,
257 doc=
"Subtask for masking streaks. Only used if doMaskStreaks is True. "
258 "Adds a mask plane to an exposure, with the mask plane name set by streakMaskName.",
260 streakBinFactor = pexConfig.Field(
263 doc=
"Bin scale factor to use when rerunning detection for masking streaks. "
264 "Only used if doMaskStreaks is True.",
266 writeStreakInfo = pexConfig.Field(
269 doc=
"Record the parameters of any detected streaks. For LSST, this should be turned off except for "
272 findGlints = pexConfig.ConfigurableField(
273 target=FindGlintTrailsTask,
274 doc=
"Subtask for finding glint trails, usually caused by satellites or debris."
276 writeGlintInfo = pexConfig.Field(
279 doc=
"Record the parameters of any detected glint trails."
281 setPrimaryFlags = pexConfig.ConfigurableField(
282 target=SetPrimaryFlagsTask,
283 doc=
"Task to add isPrimary and deblending-related flags to the catalog."
287 doc=
"Sources with any of these flags set are removed before writing the output catalog.",
288 default=(
"base_PixelFlags_flag_offimage",
289 "base_PixelFlags_flag_edge",
290 "base_PixelFlags_flag_interpolatedCenterAll",
291 "base_PixelFlags_flag_badCenterAll",
292 "base_PixelFlags_flag_edgeCenterAll",
293 "base_PixelFlags_flag_nodataCenterAll",
294 "base_PixelFlags_flag_saturatedCenterAll",
295 "base_PixelFlags_flag_saturated_templateCenterAll",
300 doc=
"Mask planes to clear before running detection.",
301 default=(
"DETECTED",
"DETECTED_NEGATIVE",
"NOT_DEBLENDED",
"STREAK"),
303 raiseOnBadSubtractionRatio = pexConfig.Field(
306 doc=
"Raise an error if the ratio of power in detected footprints"
307 " on the difference image to the power in footprints on the science"
308 " image exceeds ``badSubtractionRatioThreshold``",
310 badSubtractionRatioThreshold = pexConfig.Field(
313 doc=
"Maximum ratio of power in footprints on the difference image to"
314 " the same footprints on the science image."
315 "Only used if ``raiseOnBadSubtractionRatio`` is set",
317 badSubtractionVariationThreshold = pexConfig.Field(
320 doc=
"Maximum standard deviation of the ratio of power in footprints on"
321 " the difference image to the same footprints on the science image."
322 "Only used if ``raiseOnBadSubtractionRatio`` is set",
324 raiseOnNoDiaSources = pexConfig.Field(
327 doc=
"Raise an algorithm error if no diaSources are detected.",
329 run_sattle = pexConfig.Field(
332 doc=
"If true, dia source bounding boxes will be sent for verification"
333 "to the sattle service."
335 sattle_historical = pexConfig.Field(
338 doc=
"If re-running a pipeline that requires sattle, this should be set "
339 "to True. This will populate sattle's cache with the historic data "
340 "closest in time to the exposure."
342 idGenerator = DetectorVisitIdGeneratorConfig.make_field()
344 def setDefaults(self):
349 self.subtractInitialBackground.binSize = 8
350 self.subtractInitialBackground.useApprox =
False
351 self.subtractInitialBackground.statisticsProperty =
"MEDIAN"
352 self.subtractInitialBackground.doFilterSuperPixels =
True
353 self.subtractInitialBackground.ignoredPixelMask = [
"BAD",
361 self.subtractFinalBackground.binSize = 40
362 self.subtractFinalBackground.useApprox =
False
363 self.subtractFinalBackground.statisticsProperty =
"MEDIAN"
364 self.subtractFinalBackground.doFilterSuperPixels =
True
365 self.subtractFinalBackground.ignoredPixelMask = [
"BAD",
372 self.detection.thresholdPolarity =
"both"
373 self.detection.thresholdValue = 5.0
374 self.detection.reEstimateBackground =
False
375 self.detection.thresholdType =
"pixel_stdev"
376 self.detection.excludeMaskPlanes = []
379 self.streakDetection.thresholdType = self.detection.thresholdType
380 self.streakDetection.reEstimateBackground =
False
381 self.streakDetection.excludeMaskPlanes = self.detection.excludeMaskPlanes
382 self.streakDetection.thresholdValue = self.detection.thresholdValue
384 self.streakDetection.thresholdPolarity =
"positive"
386 self.streakDetection.nSigmaToGrow = 0
389 self.maskStreaks.onlyMaskDetected =
False
391 self.maskStreaks.maxStreakWidth = 100
394 self.maskStreaks.maxFitIter = 10
396 self.maskStreaks.nSigmaMask = 2
399 self.maskStreaks.absMinimumKernelHeight = 2
401 self.measurement.plugins.names |= [
"ext_trailedSources_Naive",
402 "base_LocalPhotoCalib",
404 "ext_shapeHSM_HsmSourceMoments",
405 "ext_shapeHSM_HsmPsfMoments",
406 "base_ClassificationSizeExtendedness",
408 self.measurement.slots.psfShape =
"ext_shapeHSM_HsmPsfMoments"
409 self.measurement.slots.shape =
"ext_shapeHSM_HsmSourceMoments"
410 self.measurement.plugins[
"base_SdssCentroid"].maxDistToPeak = 5.0
411 self.forcedMeasurement.plugins = [
"base_TransformedCentroid",
"base_PsfFlux"]
412 self.forcedMeasurement.copyColumns = {
413 "id":
"objectId",
"parent":
"parentObjectId",
"coord_ra":
"coord_ra",
"coord_dec":
"coord_dec"}
414 self.forcedMeasurement.slots.centroid =
"base_TransformedCentroid"
415 self.forcedMeasurement.slots.shape =
None
418 self.measurement.plugins[
"base_PixelFlags"].masksFpAnywhere = [
419 "STREAK",
"INJECTED",
"INJECTED_TEMPLATE",
"HIGH_VARIANCE",
"SATURATED_TEMPLATE"]
420 self.measurement.plugins[
"base_PixelFlags"].masksFpCenter = [
421 "STREAK",
"INJECTED",
"INJECTED_TEMPLATE",
"HIGH_VARIANCE",
"SATURATED_TEMPLATE"]
422 self.skySources.avoidMask = [
"DETECTED",
"DETECTED_NEGATIVE",
"BAD",
"NO_DATA",
"EDGE"]
428 if not os.getenv(
"SATTLE_URI_BASE"):
429 raise pexConfig.FieldValidationError(DetectAndMeasureConfig.run_sattle, self,
430 "Sattle requested but SATTLE_URI_BASE "
431 "environment variable not set.")
434class DetectAndMeasureTask(lsst.pipe.base.PipelineTask):
435 """Detect and measure sources on a difference image.
437 ConfigClass = DetectAndMeasureConfig
438 _DefaultName =
"detectAndMeasure"
440 def __init__(self, **kwargs):
441 super().__init__(**kwargs)
442 self.schema = afwTable.SourceTable.makeMinimalSchema()
445 if self.config.doSubtractBackground:
446 self.makeSubtask(
"subtractInitialBackground")
447 self.makeSubtask(
"subtractFinalBackground")
448 self.makeSubtask(
"detection", schema=self.schema)
449 if self.config.doDeblend:
450 self.makeSubtask(
"deblend", schema=self.schema)
451 self.makeSubtask(
"setPrimaryFlags", schema=self.schema, isSingleFrame=
True)
452 self.makeSubtask(
"measurement", schema=self.schema,
453 algMetadata=self.algMetadata)
454 if self.config.doApCorr:
455 self.makeSubtask(
"applyApCorr", schema=self.measurement.schema)
456 if self.config.doForcedMeasurement:
457 self.schema.addField(
458 "ip_diffim_forced_PsfFlux_instFlux",
"D",
459 "Forced PSF flux measured on the direct image.",
461 self.schema.addField(
462 "ip_diffim_forced_PsfFlux_instFluxErr",
"D",
463 "Forced PSF flux error measured on the direct image.",
465 self.schema.addField(
466 "ip_diffim_forced_PsfFlux_area",
"F",
467 "Forced PSF flux effective area of PSF.",
469 self.schema.addField(
470 "ip_diffim_forced_PsfFlux_flag",
"Flag",
471 "Forced PSF flux general failure flag.")
472 self.schema.addField(
473 "ip_diffim_forced_PsfFlux_flag_noGoodPixels",
"Flag",
474 "Forced PSF flux not enough non-rejected pixels in data to attempt the fit.")
475 self.schema.addField(
476 "ip_diffim_forced_PsfFlux_flag_edge",
"Flag",
477 "Forced PSF flux object was too close to the edge of the image to use the full PSF model.")
478 self.schema.addField(
479 "ip_diffim_forced_template_PsfFlux_instFlux",
"D",
480 "Forced PSF flux measured on the template image.",
482 self.schema.addField(
483 "ip_diffim_forced_template_PsfFlux_instFluxErr",
"D",
484 "Forced PSF flux error measured on the template image.",
486 self.schema.addField(
487 "ip_diffim_forced_template_PsfFlux_area",
"F",
488 "Forced template PSF flux effective area of PSF.",
490 self.schema.addField(
491 "ip_diffim_forced_template_PsfFlux_flag",
"Flag",
492 "Forced template PSF flux general failure flag.")
493 self.schema.addField(
494 "ip_diffim_forced_template_PsfFlux_flag_noGoodPixels",
"Flag",
495 "Forced template PSF flux not enough non-rejected pixels in data to attempt the fit.")
496 self.schema.addField(
497 "ip_diffim_forced_template_PsfFlux_flag_edge",
"Flag",
498 """Forced template PSF flux object was too close to the edge of the image """
499 """to use the full PSF model.""")
500 self.makeSubtask(
"forcedMeasurement", refSchema=self.schema)
502 self.schema.addField(
"refMatchId",
"L",
"unique id of reference catalog match")
503 self.schema.addField(
"srcMatchId",
"L",
"unique id of source match")
506 self.makeSubtask(
"skySources", schema=self.schema)
507 if self.config.doMaskStreaks:
508 self.makeSubtask(
"maskStreaks")
509 self.makeSubtask(
"streakDetection")
510 self.makeSubtask(
"findGlints")
511 self.schema.addField(
"glint_trail",
"Flag",
"DiaSource is part of a glint trail.")
512 self.schema.addField(
"reliability", type=
"F", doc=
"Reliability score of the DiaSource")
519 for flag
in self.config.badSourceFlags:
520 if flag
not in self.schema:
521 raise pipeBase.InvalidQuantumError(
"Field %s not in schema" % flag)
525 self.outputSchema.getTable().setMetadata(self.algMetadata)
527 def runQuantum(self, butlerQC: pipeBase.QuantumContext,
528 inputRefs: pipeBase.InputQuantizedConnection,
529 outputRefs: pipeBase.OutputQuantizedConnection):
530 inputs = butlerQC.get(inputRefs)
531 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId)
532 idFactory = idGenerator.make_table_id_factory()
535 measurementResults = pipeBase.Struct(
536 subtractedMeasuredExposure=
None,
539 differenceBackground=
None,
542 self.run(**inputs, idFactory=idFactory, measurementResults=measurementResults)
543 except pipeBase.AlgorithmError
as e:
544 error = pipeBase.AnnotatedPartialOutputsError.annotate(
547 measurementResults.subtractedMeasuredExposure,
548 measurementResults.diaSources,
549 measurementResults.maskedStreaks,
552 butlerQC.put(measurementResults, outputRefs)
554 butlerQC.put(measurementResults, outputRefs)
557 def run(self, science, matchedTemplate, difference, kernelSources,
558 idFactory=None, measurementResults=None):
559 """Detect and measure sources on a difference image.
561 The difference image will be convolved with a gaussian approximation of
562 the PSF to form a maximum likelihood image for detection.
563 Close positive and negative detections will optionally be merged into
565 Sky sources, or forced detections in background regions, will optionally
566 be added, and the configured measurement algorithm will be run on all
571 science : `lsst.afw.image.ExposureF`
572 Science exposure that the template was subtracted from.
573 matchedTemplate : `lsst.afw.image.ExposureF`
574 Warped and PSF-matched template that was used produce the
576 difference : `lsst.afw.image.ExposureF`
577 Result of subtracting template from the science image.
578 kernelSources : `lsst.afw.table.SourceCatalog`
579 Final selection of sources that was used for psf matching.
580 idFactory : `lsst.afw.table.IdFactory`, optional
581 Generator object used to assign ids to detected sources in the
582 difference image. Ids from this generator are not set until after
583 deblending and merging positive/negative peaks.
584 measurementResults : `lsst.pipe.base.Struct`, optional
585 Result struct that is modified to allow saving of partial outputs
586 for some failure conditions. If the task completes successfully,
587 this is also returned.
591 measurementResults : `lsst.pipe.base.Struct`
593 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF`
594 Subtracted exposure with detection mask applied.
595 ``diaSources`` : `lsst.afw.table.SourceCatalog`
596 The catalog of detected sources.
597 ``differenceBackground`` : `lsst.afw.math.BackgroundList`
598 Background that was subtracted from the difference image.
600 if measurementResults
is None:
601 measurementResults = pipeBase.Struct()
602 if idFactory
is None:
605 if self.config.doSubtractBackground:
607 detectionExposure = difference.clone()
608 background = self.subtractInitialBackground.run(detectionExposure).background
610 detectionExposure = difference
613 self._prepareInputs(detectionExposure)
618 table = afwTable.SourceTable.make(self.schema)
619 results = self.detection.run(
621 exposure=detectionExposure,
623 background=background
626 if self.config.doSubtractBackground:
630 difference.setMask(detectionExposure.mask)
631 background = self.subtractFinalBackground.run(difference).background
634 table = afwTable.SourceTable.make(self.schema)
635 results = self.detection.run(
639 background=background
641 measurementResults.differenceBackground = background
643 if self.config.doDeblend:
644 sources, positives, negatives = self._deblend(difference,
650 results.positive.makeSources(positives)
652 results.negative.makeSources(negatives)
653 sources = results.sources
655 self.processResults(science, matchedTemplate, difference,
656 sources, idFactory, kernelSources,
659 measurementResults=measurementResults)
660 return measurementResults
662 def _prepareInputs(self, difference):
663 """Ensure that we start with an empty detection and deblended mask.
667 difference : `lsst.afw.image.ExposureF`
668 The difference image that will be used for detecting diaSources.
669 The mask plane will be modified in place.
673 lsst.pipe.base.UpstreamFailureNoWorkFound
674 If the PSF is not usable for measurement.
677 sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius()
679 raise pipeBase.UpstreamFailureNoWorkFound(
"Invalid PSF detected! PSF width evaluates to NaN.")
681 mask = difference.mask
682 for mp
in self.config.clearMaskPlanes:
683 if mp
not in mask.getMaskPlaneDict():
684 mask.addMaskPlane(mp)
685 mask &= ~mask.getPlaneBitMask(self.config.clearMaskPlanes)
687 def processResults(self, science, matchedTemplate, difference, sources, idFactory, kernelSources,
688 positives=None, negatives=None, measurementResults=None):
689 """Measure and process the results of source detection.
693 science : `lsst.afw.image.ExposureF`
694 Science exposure that the template was subtracted from.
695 matchedTemplate : `lsst.afw.image.ExposureF`
696 Warped and PSF-matched template that was used produce the
698 difference : `lsst.afw.image.ExposureF`
699 Result of subtracting template from the science image.
700 sources : `lsst.afw.table.SourceCatalog`
701 Detected sources on the difference exposure.
702 idFactory : `lsst.afw.table.IdFactory`
703 Generator object used to assign ids to detected sources in the
705 kernelSources : `lsst.afw.table.SourceCatalog`
706 Final selection of sources that was used for psf matching.
707 positives : `lsst.afw.table.SourceCatalog`, optional
708 Positive polarity footprints.
709 negatives : `lsst.afw.table.SourceCatalog`, optional
710 Negative polarity footprints.
711 measurementResults : `lsst.pipe.base.Struct`, optional
712 Result struct that is modified to allow saving of partial outputs
713 for some failure conditions. If the task completes successfully,
714 this is also returned.
718 measurementResults : `lsst.pipe.base.Struct`
720 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF`
721 Subtracted exposure with detection mask applied.
722 ``diaSources`` : `lsst.afw.table.SourceCatalog`
723 The catalog of detected sources.
725 if measurementResults
is None:
726 measurementResults = pipeBase.Struct()
727 self.metadata[
"nUnmergedDiaSources"] = len(sources)
728 if self.config.doMerge:
730 if len(positives) > 0:
731 peakSchema = positives[0].getFootprint().peaks.schema
732 elif len(negatives) > 0:
733 peakSchema = negatives[0].getFootprint().peaks.schema
735 peakSchema = afwDetection.PeakTable.makeMinimalSchema()
736 mergeList = afwDetection.FootprintMergeList(self.schema,
737 [
"positive",
"negative"], peakSchema)
742 mergeList.addCatalog(initialDiaSources.table, positives,
"positive", minNewPeakDist=0)
743 mergeList.addCatalog(initialDiaSources.table, negatives,
"negative", minNewPeakDist=0)
744 mergeList.getFinalSources(initialDiaSources)
747 initialDiaSources[
"is_negative"] = initialDiaSources[
"merge_footprint_negative"] & \
748 ~initialDiaSources[
"merge_footprint_positive"]
749 self.log.info(
"Merging detections into %d sources", len(initialDiaSources))
751 initialDiaSources = sources
755 for source
in initialDiaSources:
758 initialDiaSources.getTable().setIdFactory(idFactory)
759 initialDiaSources.setMetadata(self.algMetadata)
761 self.metadata[
"nMergedDiaSources"] = len(initialDiaSources)
763 if self.config.doMaskStreaks:
764 streakInfo = self._runStreakMasking(difference)
766 if self.config.doSkySources:
767 self.addSkySources(initialDiaSources, difference.mask, difference.info.id)
769 if not initialDiaSources.isContiguous():
770 initialDiaSources = initialDiaSources.copy(deep=
True)
772 self.measureDiaSources(initialDiaSources, science, difference, matchedTemplate)
775 diaSources = self._removeBadSources(initialDiaSources)
777 if self.config.run_sattle:
778 diaSources = self.filterSatellites(diaSources, science)
781 diaSources, trail_parameters = self._find_glint_trails(diaSources)
782 if self.config.writeGlintInfo:
783 measurementResults.mergeItems(trail_parameters,
'glintTrailInfo')
785 if self.config.doForcedMeasurement:
786 self.measureForcedSources(diaSources, science, difference.getWcs())
787 self.measureForcedSources(diaSources, matchedTemplate, difference.getWcs(),
794 difference.image.array[difference.mask.array & difference.mask.getPlaneBitMask(
'NO_DATA') > 0] = 0
796 measurementResults.subtractedMeasuredExposure = difference
798 if self.config.doMaskStreaks
and self.config.writeStreakInfo:
799 measurementResults.mergeItems(streakInfo,
'maskedStreaks')
801 self.calculateMetrics(science, difference, diaSources, kernelSources)
803 if np.count_nonzero(~diaSources[
"sky_source"]) > 0:
804 measurementResults.diaSources = diaSources
805 elif self.config.raiseOnNoDiaSources:
807 elif len(diaSources) > 0:
810 measurementResults.diaSources = diaSources
811 self.log.info(
"Measured %d diaSources and %d sky sources",
812 np.count_nonzero(~diaSources[
"sky_source"]),
813 np.count_nonzero(diaSources[
"sky_source"])
815 return measurementResults
817 def _deblend(self, difference, positiveFootprints, negativeFootprints):
818 """Deblend the positive and negative footprints and return a catalog
819 containing just the children, and the deblended footprints.
823 difference : `lsst.afw.image.Exposure`
824 Result of subtracting template from the science image.
825 positiveFootprints, negativeFootprints : `lsst.afw.detection.FootprintSet`
826 Positive and negative polarity footprints measured on
827 ``difference`` to be deblended separately.
831 sources : `lsst.afw.table.SourceCatalog`
832 Positive and negative deblended children.
833 positives, negatives : `lsst.afw.table.SourceCatalog`
834 Deblended positive and negative polarity sources with footprints
835 detected on ``difference``.
838 def deblend(footprints, negative=False):
839 """Deblend a positive or negative footprint set,
840 and return the deblended children.
844 footprints : `lsst.afw.detection.FootprintSet`
846 Set True if the footprints contain negative fluxes
850 sources : `lsst.afw.table.SourceCatalog`
853 footprints.makeSources(sources)
856 difference_inverted = difference.clone()
857 difference_inverted.image *= -1
858 self.deblend.run(exposure=difference_inverted, sources=sources)
859 children = sources[sources[
"parent"] != 0]
861 for child
in children:
862 footprint = child.getFootprint()
863 array = footprint.getImageArray()
866 self.deblend.run(exposure=difference, sources=sources)
867 self.setPrimaryFlags.run(sources)
868 children = sources[
"detect_isDeblendedSource"] == 1
869 sources = sources[children].copy(deep=
True)
871 sources[
'parent'] = 0
872 return sources.copy(deep=
True)
874 positives =
deblend(positiveFootprints)
875 negatives =
deblend(negativeFootprints, negative=
True)
878 sources.reserve(len(positives) + len(negatives))
879 sources.extend(positives, deep=
True)
880 sources.extend(negatives, deep=
True)
881 if len(negatives) > 0:
882 sources[-len(negatives):][
"is_negative"] =
True
883 return sources, positives, negatives
885 def _removeBadSources(self, diaSources):
886 """Remove unphysical diaSources from the catalog.
890 diaSources : `lsst.afw.table.SourceCatalog`
891 The catalog of detected sources.
895 diaSources : `lsst.afw.table.SourceCatalog`
896 The updated catalog of detected sources, with any source that has a
897 flag in ``config.badSourceFlags`` set removed.
899 selector = np.ones(len(diaSources), dtype=bool)
900 for flag
in self.config.badSourceFlags:
901 flags = diaSources[flag]
902 nBad = np.count_nonzero(flags)
904 self.log.debug(
"Found %d unphysical sources with flag %s.", nBad, flag)
906 nBadTotal = np.count_nonzero(~selector)
907 self.metadata[
"nRemovedBadFlaggedSources"] = nBadTotal
908 self.log.info(
"Removed %d unphysical sources.", nBadTotal)
909 return diaSources[selector].copy(deep=
True)
911 def _find_glint_trails(self, diaSources):
912 """Define a new flag column for diaSources that are in a glint trail.
916 diaSources : `lsst.afw.table.SourceCatalog`
917 The catalog of detected sources.
921 diaSources : `lsst.afw.table.SourceCatalog`
922 The updated catalog of detected sources, with a new bool column
923 called 'glint_trail' added.
925 trail_parameters : `dict`
926 Parameters of all the trails that were found.
928 if self.config.doSkySources:
930 candidateDiaSources = diaSources[~diaSources[
"sky_source"]].copy(deep=
True)
932 candidateDiaSources = diaSources
933 trailed_glints = self.findGlints.run(candidateDiaSources)
934 glint_mask = [
True if id
in trailed_glints.trailed_ids
else False for id
in diaSources[
'id']]
935 if np.any(glint_mask):
936 diaSources[
'glint_trail'] = np.array(glint_mask)
938 slopes = np.array([trail.slope
for trail
in trailed_glints.parameters])
939 intercepts = np.array([trail.intercept
for trail
in trailed_glints.parameters])
940 stderrs = np.array([trail.stderr
for trail
in trailed_glints.parameters])
941 lengths = np.array([trail.length
for trail
in trailed_glints.parameters])
942 angles = np.array([trail.angle
for trail
in trailed_glints.parameters])
943 parameters = {
'slopes': slopes,
'intercepts': intercepts,
'stderrs': stderrs,
'lengths': lengths,
946 trail_parameters = pipeBase.Struct(glintTrailInfo=parameters)
948 return diaSources, trail_parameters
950 def addSkySources(self, diaSources, mask, seed,
952 """Add sources in empty regions of the difference image
953 for measuring the background.
957 diaSources : `lsst.afw.table.SourceCatalog`
958 The catalog of detected sources.
959 mask : `lsst.afw.image.Mask`
960 Mask plane for determining regions where Sky sources can be added.
962 Seed value to initialize the random number generator.
965 subtask = self.skySources
966 if subtask.config.nSources <= 0:
967 self.metadata[f
"n_{subtask.getName()}"] = 0
969 skySourceFootprints = subtask.run(mask=mask, seed=seed, catalog=diaSources)
970 self.metadata[f
"n_{subtask.getName()}"] = len(skySourceFootprints)
972 def measureDiaSources(self, diaSources, science, difference, matchedTemplate):
973 """Use (matched) template and science image to constrain dipole fitting.
977 diaSources : `lsst.afw.table.SourceCatalog`
978 The catalog of detected sources.
979 science : `lsst.afw.image.ExposureF`
980 Science exposure that the template was subtracted from.
981 difference : `lsst.afw.image.ExposureF`
982 Result of subtracting template from the science image.
983 matchedTemplate : `lsst.afw.image.ExposureF`
984 Warped and PSF-matched template that was used produce the
988 for mp
in self.config.measurement.plugins[
"base_PixelFlags"].masksFpAnywhere:
989 difference.mask.addMaskPlane(mp)
992 self.measurement.run(diaSources, difference, science, matchedTemplate)
993 if self.config.doApCorr:
994 apCorrMap = difference.getInfo().getApCorrMap()
995 if apCorrMap
is None:
996 self.log.warning(
"Difference image does not have valid aperture correction; skipping.")
998 self.applyApCorr.run(
1000 apCorrMap=apCorrMap,
1003 def measureForcedSources(self, diaSources, image, wcs, template=False):
1004 """Perform forced measurement of the diaSources on a direct image.
1008 diaSources : `lsst.afw.table.SourceCatalog`
1009 The catalog of detected sources.
1010 image: `lsst.afw.image.ExposureF`
1011 Exposure that the forced measurement is being performed on
1012 wcs : `lsst.afw.geom.SkyWcs`
1013 Coordinate system definition (wcs) for the exposure.
1015 Is the forced measurement being performed on the template?
1019 forcedSources = self.forcedMeasurement.generateMeasCat(image, diaSources, wcs)
1020 self.forcedMeasurement.run(forcedSources, image, diaSources, wcs)
1023 base_key =
'ip_diffim_forced_template_PsfFlux'
1025 base_key =
'ip_diffim_forced_PsfFlux'
1027 mapper.addMapping(forcedSources.schema.find(
"base_PsfFlux_instFlux")[0],
1028 f
"{base_key}_instFlux",
True)
1029 mapper.addMapping(forcedSources.schema.find(
"base_PsfFlux_instFluxErr")[0],
1030 f
"{base_key}_instFluxErr",
True)
1031 mapper.addMapping(forcedSources.schema.find(
"base_PsfFlux_area")[0],
1032 f
"{base_key}_area",
True)
1033 mapper.addMapping(forcedSources.schema.find(
"base_PsfFlux_flag")[0],
1034 f
"{base_key}_flag",
True)
1035 mapper.addMapping(forcedSources.schema.find(
"base_PsfFlux_flag_noGoodPixels")[0],
1036 f
"{base_key}_flag_noGoodPixels",
True)
1037 mapper.addMapping(forcedSources.schema.find(
"base_PsfFlux_flag_edge")[0],
1038 f
"{base_key}_flag_edge",
True)
1039 for diaSource, forcedSource
in zip(diaSources, forcedSources):
1040 diaSource.assign(forcedSource, mapper)
1042 def calculateMetrics(self, science, difference, diaSources, kernelSources):
1043 """Add difference image QA metrics to the Task metadata.
1045 This may be used to produce corresponding metrics (see
1046 lsst.analysis.tools.tasks.diffimTaskDetectorVisitMetricAnalysis).
1050 science : `lsst.afw.image.ExposureF`
1051 Science exposure that was subtracted.
1052 difference : `lsst.afw.image.Exposure`
1053 The target difference image to calculate metrics for.
1054 diaSources : `lsst.afw.table.SourceCatalog`
1055 The catalog of detected sources.
1056 kernelSources : `lsst.afw.table.SourceCatalog`
1057 Final selection of sources that was used for psf matching.
1059 mask = difference.mask
1060 badPix = (mask.array & mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) > 0
1061 self.metadata[
"nGoodPixels"] = np.sum(~badPix)
1062 self.metadata[
"nBadPixels"] = np.sum(badPix)
1063 detPosPix = (mask.array & mask.getPlaneBitMask(
"DETECTED")) > 0
1064 detNegPix = (mask.array & mask.getPlaneBitMask(
"DETECTED_NEGATIVE")) > 0
1065 self.metadata[
"nPixelsDetectedPositive"] = np.sum(detPosPix)
1066 self.metadata[
"nPixelsDetectedNegative"] = np.sum(detNegPix)
1069 self.metadata[
"nBadPixelsDetectedPositive"] = np.sum(detPosPix)
1070 self.metadata[
"nBadPixelsDetectedNegative"] = np.sum(detNegPix)
1072 metricsMaskPlanes = list(mask.getMaskPlaneDict().keys())
1073 for maskPlane
in metricsMaskPlanes:
1075 self.metadata[
"%s_mask_fraction"%maskPlane.lower()] = evaluateMaskFraction(mask, maskPlane)
1076 except InvalidParameterError:
1077 self.metadata[
"%s_mask_fraction"%maskPlane.lower()] = -1
1078 self.log.info(
"Unable to calculate metrics for mask plane %s: not in image"%maskPlane)
1080 if self.config.doSkySources:
1081 skySources = diaSources[diaSources[
"sky_source"]]
1084 metrics = computeDifferenceImageMetrics(science, difference, kernelSources, sky_sources=skySources)
1086 self.metadata[
"residualFootprintRatioMean"] = metrics.differenceFootprintRatioMean
1087 self.metadata[
"residualFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev
1088 self.metadata[
"differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean
1089 self.metadata[
"differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev
1090 self.log.info(
"Mean, stdev of ratio of difference to science "
1091 "pixels in star footprints: %5.4f, %5.4f",
1092 self.metadata[
"residualFootprintRatioMean"],
1093 self.metadata[
"residualFootprintRatioStdev"])
1094 if self.config.raiseOnBadSubtractionRatio:
1095 if metrics.differenceFootprintRatioMean > self.config.badSubtractionRatioThreshold:
1097 threshold=self.config.badSubtractionRatioThreshold)
1098 if metrics.differenceFootprintRatioStdev > self.config.badSubtractionVariationThreshold:
1100 threshold=self.config.badSubtractionVariationThreshold)
1102 def getSattleDiaSourceAllowlist(self, diaSources, science):
1103 """Query the sattle service and determine which diaSources are allowed.
1107 diaSources : `lsst.afw.table.SourceCatalog`
1108 The catalog of detected sources.
1109 science : `lsst.afw.image.ExposureF`
1110 Science exposure that was subtracted.
1114 allow_list : `list` of `int`
1115 diaSourceIds of diaSources that can be made public.
1120 Raised if sattle call does not return success.
1122 wcs = science.getWcs()
1123 visit_info = science.getInfo().getVisitInfo()
1124 visit_id = visit_info.getId()
1125 sattle_uri_base = os.getenv(
'SATTLE_URI_BASE')
1127 dia_sources_json = []
1128 for source
in diaSources:
1129 source_bbox = source.getFootprint().getBBox()
1130 corners = wcs.pixelToSky([
lsst.geom.Point2D(c)
for c
in source_bbox.getCorners()])
1131 bbox_radec = [[pt.getRa().asDegrees(), pt.getDec().asDegrees()]
for pt
in corners]
1132 dia_sources_json.append({
"diasource_id": source[
"id"],
"bbox": bbox_radec})
1134 payload = {
"visit_id": visit_id,
"detector_id": science.getDetector().getId(),
1135 "diasources": dia_sources_json,
"historical": self.config.sattle_historical}
1137 sattle_output = requests.put(f
'{sattle_uri_base}/diasource_allow_list',
1141 if sattle_output.status_code == 404:
1142 self.log.warning(f
'Visit {visit_id} not found in sattle cache, re-sending')
1143 populate_sattle_visit_cache(visit_info, historical=self.config.sattle_historical)
1144 sattle_output = requests.put(f
'{sattle_uri_base}/diasource_allow_list', json=payload)
1146 sattle_output.raise_for_status()
1148 return sattle_output.json()[
'allow_list']
1150 def filterSatellites(self, diaSources, science):
1151 """Remove diaSources overlapping predicted satellite positions.
1155 diaSources : `lsst.afw.table.SourceCatalog`
1156 The catalog of detected sources.
1157 science : `lsst.afw.image.ExposureF`
1158 Science exposure that was subtracted.
1162 filterdDiaSources : `lsst.afw.table.SourceCatalog`
1163 Filtered catalog of diaSources
1166 allow_list = self.getSattleDiaSourceAllowlist(diaSources, science)
1169 allow_set = set(allow_list)
1170 allowed_ids = [source[
'id']
in allow_set
for source
in diaSources]
1171 diaSources = diaSources[np.array(allowed_ids)].copy(deep=
True)
1173 self.log.warning(
'Sattle allowlist is empty, all diaSources removed')
1174 diaSources = diaSources[0:0].copy(deep=
True)
1177 def _runStreakMasking(self, difference):
1178 """Do streak masking and optionally save the resulting streak
1179 fit parameters in a catalog.
1181 Only returns non-empty streakInfo if self.config.writeStreakInfo
1182 is set. The difference image is binned by self.config.streakBinFactor
1183 (and detection is run a second time) so that regions with lower
1184 surface brightness streaks are more likely to fall above the
1185 detection threshold.
1189 difference: `lsst.afw.image.Exposure`
1190 The exposure in which to search for streaks. Must have a detection
1195 streakInfo: `lsst.pipe.base.Struct`
1196 ``rho`` : `np.ndarray`
1197 Angle of detected streak.
1198 ``theta`` : `np.ndarray`
1199 Distance from center of detected streak.
1200 ``sigma`` : `np.ndarray`
1201 Width of streak profile.
1202 ``reducedChi2`` : `np.ndarray`
1203 Reduced chi2 of the best-fit streak profile.
1204 ``modelMaximum`` : `np.ndarray`
1205 Peak value of the fit line profile.
1207 maskedImage = difference.maskedImage
1210 self.config.streakBinFactor,
1211 self.config.streakBinFactor)
1212 binnedExposure = afwImage.ExposureF(binnedMaskedImage.getBBox())
1213 binnedExposure.setMaskedImage(binnedMaskedImage)
1215 binnedExposure.mask &= ~binnedExposure.mask.getPlaneBitMask(
'DETECTED')
1217 sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius()
1218 _table = afwTable.SourceTable.make(afwTable.SourceTable.makeMinimalSchema())
1219 self.streakDetection.run(table=_table, exposure=binnedExposure, doSmooth=
True,
1220 sigma=sigma/self.config.streakBinFactor)
1221 binnedDetectedMaskPlane = binnedExposure.mask.array & binnedExposure.mask.getPlaneBitMask(
'DETECTED')
1222 rescaledDetectedMaskPlane = binnedDetectedMaskPlane.repeat(self.config.streakBinFactor,
1223 axis=0).repeat(self.config.streakBinFactor,
1226 streakMaskedImage = maskedImage.clone()
1227 ysize, xsize = rescaledDetectedMaskPlane.shape
1228 streakMaskedImage.mask.array[:ysize, :xsize] |= rescaledDetectedMaskPlane
1230 streaks = self.maskStreaks.run(streakMaskedImage)
1231 streakMaskPlane = streakMaskedImage.mask.array & streakMaskedImage.mask.getPlaneBitMask(
'STREAK')
1233 maskedImage.mask.array |= streakMaskPlane
1235 if self.config.writeStreakInfo:
1236 rhos = np.array([line.rho
for line
in streaks.lines])
1237 thetas = np.array([line.theta
for line
in streaks.lines])
1238 sigmas = np.array([line.sigma
for line
in streaks.lines])
1239 chi2s = np.array([line.reducedChi2
for line
in streaks.lines])
1240 modelMaximums = np.array([line.modelMaximum
for line
in streaks.lines])
1241 streakInfo = {
'rho': rhos,
'theta': thetas,
'sigma': sigmas,
'reducedChi2': chi2s,
1242 'modelMaximum': modelMaximums}
1244 streakInfo = {
'rho': np.array([]),
'theta': np.array([]),
'sigma': np.array([]),
1245 'reducedChi2': np.array([]),
'modelMaximum': np.array([])}
1246 return pipeBase.Struct(maskedStreaks=streakInfo)
1250 scoreExposure = pipeBase.connectionTypes.Input(
1251 doc=
"Maximum likelihood image for detection.",
1252 dimensions=(
"instrument",
"visit",
"detector"),
1253 storageClass=
"ExposureF",
1254 name=
"{fakesType}{coaddName}Diff_scoreExp",
1258class DetectAndMeasureScoreConfig(DetectAndMeasureConfig,
1259 pipelineConnections=DetectAndMeasureScoreConnections):
1263class DetectAndMeasureScoreTask(DetectAndMeasureTask):
1264 """Detect DIA sources using a score image,
1265 and measure the detections on the difference image.
1267 Source detection is run on the supplied score, or maximum likelihood,
1268 image. Note that no additional convolution will be done in this case.
1269 Close positive and negative detections will optionally be merged into
1271 Sky sources, or forced detections in background regions, will optionally
1272 be added, and the configured measurement algorithm will be run on all
1275 ConfigClass = DetectAndMeasureScoreConfig
1276 _DefaultName =
"detectAndMeasureScore"
1279 def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources,
1281 """Detect and measure sources on a score image.
1285 science : `lsst.afw.image.ExposureF`
1286 Science exposure that the template was subtracted from.
1287 matchedTemplate : `lsst.afw.image.ExposureF`
1288 Warped and PSF-matched template that was used produce the
1290 difference : `lsst.afw.image.ExposureF`
1291 Result of subtracting template from the science image.
1292 scoreExposure : `lsst.afw.image.ExposureF`
1293 Score or maximum likelihood difference image
1294 kernelSources : `lsst.afw.table.SourceCatalog`
1295 Final selection of sources that was used for psf matching.
1296 idFactory : `lsst.afw.table.IdFactory`, optional
1297 Generator object used to assign ids to detected sources in the
1298 difference image. Ids from this generator are not set until after
1299 deblending and merging positive/negative peaks.
1303 measurementResults : `lsst.pipe.base.Struct`
1305 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF`
1306 Subtracted exposure with detection mask applied.
1307 ``diaSources`` : `lsst.afw.table.SourceCatalog`
1308 The catalog of detected sources.
1310 if idFactory
is None:
1313 self._prepareInputs(scoreExposure)
1318 table = afwTable.SourceTable.make(self.schema)
1319 results = self.detection.run(
1321 exposure=scoreExposure,
1325 difference.mask.assign(scoreExposure.mask, scoreExposure.getBBox())
1327 if self.config.doDeblend:
1328 sources, positives, negatives = self._deblend(difference,
1332 return self.processResults(science, matchedTemplate, difference,
1333 sources, idFactory, kernelSources,
1334 positives=positives,
1335 negatives=negatives)
1339 results.positive.makeSources(positives)
1341 results.negative.makeSources(negatives)
1342 return self.processResults(science, matchedTemplate, difference,
1343 results.sources, idFactory, kernelSources,
1344 positives=positives,
1345 negatives=negatives)
A mapping between the keys of two Schemas, used to copy data between them.
Class for storing ordered metadata with comments.
__init__(self, *, ratio, threshold)
std::shared_ptr< ImageT > binImage(ImageT const &inImage, int const binX, int const binY, lsst::afw::math::Property const flags=lsst::afw::math::MEAN)
run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter)