LSST Applications g042eb84c57+730a74494b,g04e9c324dd+8c5ae1fdc5,g134cb467dc+1f1e3e7524,g199a45376c+0ba108daf9,g1fd858c14a+fa7d31856b,g210f2d0738+f66ac109ec,g262e1987ae+83a3acc0e5,g29ae962dfc+d856a2cb1f,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+a1e0c9f713,g47891489e3+0d594cb711,g4d44eb3520+c57ec8f3ed,g4d7b6aa1c5+f66ac109ec,g53246c7159+8c5ae1fdc5,g56a1a4eaf3+fd7ad03fde,g64539dfbff+f66ac109ec,g67b6fd64d1+0d594cb711,g67fd3c3899+f66ac109ec,g6985122a63+0d594cb711,g74acd417e5+3098891321,g786e29fd12+668abc6043,g81db2e9a8d+98e2ab9f28,g87389fa792+8856018cbb,g89139ef638+0d594cb711,g8d7436a09f+80fda9ce03,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+a8a29bda4b,g99822b682c+e3ec3c61f9,g9d5c6a246b+0d5dac0c3d,ga41d0fce20+9243b26dd2,gbf99507273+8c5ae1fdc5,gd7ef33dd92+0d594cb711,gdab6d2f7ff+3098891321,ge410e46f29+0d594cb711,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.38
LSST Data Management Base Package
Loading...
Searching...
No Matches
detectAndMeasure.py
Go to the documentation of this file.
1# This file is part of ip_diffim.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21
22import numpy as np
23import requests
24import os
25
26import lsst.afw.detection as afwDetection
27import lsst.afw.image as afwImage
28import lsst.afw.math as afwMath
29import lsst.afw.table as afwTable
30import lsst.daf.base as dafBase
31import lsst.geom
32from lsst.ip.diffim.utils import (evaluateMaskFraction, computeDifferenceImageMetrics,
33 populate_sattle_visit_cache)
34from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask
35from lsst.meas.algorithms import FindGlintTrailsTask
36from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig
39import lsst.meas.extensions.shapeHSM
40import lsst.pex.config as pexConfig
41from lsst.pex.exceptions import InvalidParameterError
42import lsst.pipe.base as pipeBase
43import lsst.utils
44from lsst.utils.timer import timeMethod
45
46from . import DipoleFitTask
47
48__all__ = ["DetectAndMeasureConfig", "DetectAndMeasureTask",
49 "DetectAndMeasureScoreConfig", "DetectAndMeasureScoreTask"]
50
51
52class BadSubtractionError(pipeBase.AlgorithmError):
53 """Raised when the residuals in footprints of stars used to compute the
54 psf-matching kernel exceeds the configured maximum.
55 """
56 def __init__(self, *, ratio, threshold):
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}")
61 super().__init__(msg)
62 self.ratio = ratio
63 self.threshold = threshold
64
65 @property
66 def metadata(self):
67 return {"ratio": self.ratio,
68 "threshold": self.threshold
69 }
70
71
72class NoDiaSourcesError(pipeBase.AlgorithmError):
73 """Raised when there are no diaSources detected on an image difference.
74 """
75 def __init__(self):
76 msg = ("No diaSources detected!")
77 super().__init__(msg)
78
79 @property
80 def metadata(self):
81 return {}
82
83
84class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections,
85 dimensions=("instrument", "visit", "detector"),
86 defaultTemplates={"coaddName": "deep",
87 "warpTypeSuffix": "",
88 "fakesType": ""}):
89 science = pipeBase.connectionTypes.Input(
90 doc="Input science exposure.",
91 dimensions=("instrument", "visit", "detector"),
92 storageClass="ExposureF",
93 name="{fakesType}calexp"
94 )
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",
100 )
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",
106 )
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"
112 )
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",
117 )
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",
123 )
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",
129 )
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",
135 )
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",
141 )
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",
147 )
148
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")
157
158
159class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
160 pipelineConnections=DetectAndMeasureConnections):
161 """Config for DetectAndMeasureTask
162 """
163 doMerge = pexConfig.Field(
164 dtype=bool,
165 default=True,
166 doc="Merge positive and negative diaSources with grow radius "
167 "set by growFootprint"
168 )
169 doForcedMeasurement = pexConfig.Field(
170 dtype=bool,
171 default=True,
172 doc="Force photometer diaSource locations on PVI?"
173 )
174 doAddMetrics = pexConfig.Field(
175 dtype=bool,
176 default=False,
177 doc="Add columns to the source table to hold analysis metrics?"
178 )
179 doSubtractBackground = pexConfig.Field(
180 dtype=bool,
181 doc="Subtract a background model from the image before detection?",
182 default=True,
183 )
184 doWriteBackground = pexConfig.Field(
185 dtype=bool,
186 doc="Persist the fitted background model?",
187 default=False,
188 )
189 subtractInitialBackground = pexConfig.ConfigurableField(
191 doc="Task to perform intial background subtraction, before first detection pass.",
192 )
193 subtractFinalBackground = pexConfig.ConfigurableField(
195 doc="Task to perform final background subtraction, after first detection pass.",
196 )
197 detection = pexConfig.ConfigurableField(
198 target=SourceDetectionTask,
199 doc="Final source detection for diaSource measurement",
200 )
201 streakDetection = pexConfig.ConfigurableField(
202 target=SourceDetectionTask,
203 doc="Separate source detection used only for streak masking",
204 )
205 doDeblend = pexConfig.Field(
206 dtype=bool,
207 default=False,
208 doc="Deblend DIASources after detection?"
209 )
210 deblend = pexConfig.ConfigurableField(
212 doc="Task to split blended sources into their components."
213 )
214 measurement = pexConfig.ConfigurableField(
215 target=DipoleFitTask,
216 doc="Task to measure sources on the difference image.",
217 )
218 doApCorr = lsst.pex.config.Field(
219 dtype=bool,
220 default=True,
221 doc="Run subtask to apply aperture corrections"
222 )
224 target=ApplyApCorrTask,
225 doc="Task to apply aperture corrections"
226 )
227 forcedMeasurement = pexConfig.ConfigurableField(
228 target=ForcedMeasurementTask,
229 doc="Task to force photometer science image at diaSource locations.",
230 )
231 growFootprint = pexConfig.Field(
232 dtype=int,
233 default=2,
234 doc="Grow positive and negative footprints by this many pixels before merging"
235 )
236 diaSourceMatchRadius = pexConfig.Field(
237 dtype=float,
238 default=0.5,
239 doc="Match radius (in arcseconds) for DiaSource to Source association"
240 )
241 doSkySources = pexConfig.Field(
242 dtype=bool,
243 default=False,
244 doc="Generate sky sources?",
245 )
246 skySources = pexConfig.ConfigurableField(
247 target=SkyObjectsTask,
248 doc="Generate sky sources",
249 )
250 doMaskStreaks = pexConfig.Field(
251 dtype=bool,
252 default=True,
253 doc="Turn on streak masking",
254 )
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.",
259 )
260 streakBinFactor = pexConfig.Field(
261 dtype=int,
262 default=4,
263 doc="Bin scale factor to use when rerunning detection for masking streaks. "
264 "Only used if doMaskStreaks is True.",
265 )
266 writeStreakInfo = pexConfig.Field(
267 dtype=bool,
268 default=False,
269 doc="Record the parameters of any detected streaks. For LSST, this should be turned off except for "
270 "development work."
271 )
272 findGlints = pexConfig.ConfigurableField(
273 target=FindGlintTrailsTask,
274 doc="Subtask for finding glint trails, usually caused by satellites or debris."
275 )
276 writeGlintInfo = pexConfig.Field(
277 dtype=bool,
278 default=True,
279 doc="Record the parameters of any detected glint trails."
280 )
281 setPrimaryFlags = pexConfig.ConfigurableField(
282 target=SetPrimaryFlagsTask,
283 doc="Task to add isPrimary and deblending-related flags to the catalog."
284 )
285 badSourceFlags = lsst.pex.config.ListField(
286 dtype=str,
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",
296 ),
297 )
298 clearMaskPlanes = lsst.pex.config.ListField(
299 dtype=str,
300 doc="Mask planes to clear before running detection.",
301 default=("DETECTED", "DETECTED_NEGATIVE", "NOT_DEBLENDED", "STREAK"),
302 )
303 raiseOnBadSubtractionRatio = pexConfig.Field(
304 dtype=bool,
305 default=True,
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``",
309 )
310 badSubtractionRatioThreshold = pexConfig.Field(
311 dtype=float,
312 default=0.2,
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",
316 )
317 badSubtractionVariationThreshold = pexConfig.Field(
318 dtype=float,
319 default=0.4,
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",
323 )
324 raiseOnNoDiaSources = pexConfig.Field(
325 dtype=bool,
326 default=True,
327 doc="Raise an algorithm error if no diaSources are detected.",
328 )
329 run_sattle = pexConfig.Field(
330 dtype=bool,
331 default=False,
332 doc="If true, dia source bounding boxes will be sent for verification"
333 "to the sattle service."
334 )
335 sattle_historical = pexConfig.Field(
336 dtype=bool,
337 default=False,
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."
341 )
342 idGenerator = DetectorVisitIdGeneratorConfig.make_field()
343
344 def setDefaults(self):
345 # Background subtraction
346 # Use a small binsize for the first pass to reduce detections on glints
347 # and extended structures. Should not affect the detectability of
348 # faint diaSources
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",
354 "EDGE",
355 "DETECTED",
356 "DETECTED_NEGATIVE",
357 "NO_DATA",
358 ]
359 # Use a larger binsize for the final background subtraction, to reduce
360 # over-subtraction of bright objects.
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",
366 "EDGE",
367 "DETECTED",
368 "DETECTED_NEGATIVE",
369 "NO_DATA",
370 ]
371 # DiaSource Detection
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 = []
377
378 # Copy configs for binned streak detection from the base detection task
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
383 # Only detect positive streaks
384 self.streakDetection.thresholdPolarity = "positive"
385 # Do not grow detected mask for streaks
386 self.streakDetection.nSigmaToGrow = 0
387 # Set the streak mask along the entire fit line, not only where the
388 # detected mask is set.
389 self.maskStreaks.onlyMaskDetected = False
390 # Restrict streak masking from growing too large
391 self.maskStreaks.maxStreakWidth = 100
392 # Restrict the number of iterations allowed for fitting streaks
393 # When the fit is good it should solve quickly, and exit a bad fit quickly
394 self.maskStreaks.maxFitIter = 10
395 # Only mask to 2 sigma in width
396 self.maskStreaks.nSigmaMask = 2
397 # Threshold for including streaks after the Hough Transform.
398 # A lower value will detect more features that are less linear.
399 self.maskStreaks.absMinimumKernelHeight = 2
400
401 self.measurement.plugins.names |= ["ext_trailedSources_Naive",
402 "base_LocalPhotoCalib",
403 "base_LocalWcs",
404 "ext_shapeHSM_HsmSourceMoments",
405 "ext_shapeHSM_HsmPsfMoments",
406 "base_ClassificationSizeExtendedness",
407 ]
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
416
417 # Keep track of which footprints contain streaks
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"]
423
424 def validate(self):
425 super().validate()
426
427 if self.run_sattle:
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.")
432
433
434class DetectAndMeasureTask(lsst.pipe.base.PipelineTask):
435 """Detect and measure sources on a difference image.
436 """
437 ConfigClass = DetectAndMeasureConfig
438 _DefaultName = "detectAndMeasure"
439
440 def __init__(self, **kwargs):
441 super().__init__(**kwargs)
442 self.schema = afwTable.SourceTable.makeMinimalSchema()
443
444 self.algMetadata = dafBase.PropertyList()
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.",
460 units="count")
461 self.schema.addField(
462 "ip_diffim_forced_PsfFlux_instFluxErr", "D",
463 "Forced PSF flux error measured on the direct image.",
464 units="count")
465 self.schema.addField(
466 "ip_diffim_forced_PsfFlux_area", "F",
467 "Forced PSF flux effective area of PSF.",
468 units="pixel")
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.",
481 units="count")
482 self.schema.addField(
483 "ip_diffim_forced_template_PsfFlux_instFluxErr", "D",
484 "Forced PSF flux error measured on the template image.",
485 units="count")
486 self.schema.addField(
487 "ip_diffim_forced_template_PsfFlux_area", "F",
488 "Forced template PSF flux effective area of PSF.",
489 units="pixel")
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)
501
502 self.schema.addField("refMatchId", "L", "unique id of reference catalog match")
503 self.schema.addField("srcMatchId", "L", "unique id of source match")
504 # Create the sky source task for use by metrics,
505 # even if sky sources are not added to the diaSource catalog
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")
513
514 # To get the "merge_*" fields in the schema; have to re-initialize
515 # this later, once we have a peak schema post-detection.
516 lsst.afw.detection.FootprintMergeList(self.schema, ["positive", "negative"])
517
518 # Check that the schema and config are consistent
519 for flag in self.config.badSourceFlags:
520 if flag not in self.schema:
521 raise pipeBase.InvalidQuantumError("Field %s not in schema" % flag)
522
523 # initialize InitOutputs
524 self.outputSchema = afwTable.SourceCatalog(self.schema)
525 self.outputSchema.getTable().setMetadata(self.algMetadata)
526
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()
533 # Specify the fields that `annotate` needs below, to ensure they
534 # exist, even as None.
535 measurementResults = pipeBase.Struct(
536 subtractedMeasuredExposure=None,
537 diaSources=None,
538 maskedStreaks=None,
539 differenceBackground=None,
540 )
541 try:
542 self.run(**inputs, idFactory=idFactory, measurementResults=measurementResults)
543 except pipeBase.AlgorithmError as e:
544 error = pipeBase.AnnotatedPartialOutputsError.annotate(
545 e,
546 self,
547 measurementResults.subtractedMeasuredExposure,
548 measurementResults.diaSources,
549 measurementResults.maskedStreaks,
550 log=self.log
551 )
552 butlerQC.put(measurementResults, outputRefs)
553 raise error from e
554 butlerQC.put(measurementResults, outputRefs)
555
556 @timeMethod
557 def run(self, science, matchedTemplate, difference, kernelSources,
558 idFactory=None, measurementResults=None):
559 """Detect and measure sources on a difference image.
560
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
564 dipole diaSources.
565 Sky sources, or forced detections in background regions, will optionally
566 be added, and the configured measurement algorithm will be run on all
567 detections.
568
569 Parameters
570 ----------
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
575 difference image.
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.
588
589 Returns
590 -------
591 measurementResults : `lsst.pipe.base.Struct`
592
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.
599 """
600 if measurementResults is None:
601 measurementResults = pipeBase.Struct()
602 if idFactory is None:
603 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
604
605 if self.config.doSubtractBackground:
606 # Run background subtraction before clearing the mask planes
607 detectionExposure = difference.clone()
608 background = self.subtractInitialBackground.run(detectionExposure).background
609 else:
610 detectionExposure = difference
611 background = afwMath.BackgroundList()
612
613 self._prepareInputs(detectionExposure)
614
615 # Don't use the idFactory until after deblend+merge, so that we aren't
616 # generating ids that just get thrown away (footprint merge doesn't
617 # know about past ids).
618 table = afwTable.SourceTable.make(self.schema)
619 results = self.detection.run(
620 table=table,
621 exposure=detectionExposure,
622 doSmooth=True,
623 background=background
624 )
625
626 if self.config.doSubtractBackground:
627 # Run background subtraction again after detecting peaks
628 # but before measurement
629 # First update the mask using the detection image
630 difference.setMask(detectionExposure.mask)
631 background = self.subtractFinalBackground.run(difference).background
632
633 # Re-run detection to get final footprints
634 table = afwTable.SourceTable.make(self.schema)
635 results = self.detection.run(
636 table=table,
637 exposure=difference,
638 doSmooth=True,
639 background=background
640 )
641 measurementResults.differenceBackground = background
642
643 if self.config.doDeblend:
644 sources, positives, negatives = self._deblend(difference,
645 results.positive,
646 results.negative)
647
648 else:
649 positives = afwTable.SourceCatalog(self.schema)
650 results.positive.makeSources(positives)
651 negatives = afwTable.SourceCatalog(self.schema)
652 results.negative.makeSources(negatives)
653 sources = results.sources
654
655 self.processResults(science, matchedTemplate, difference,
656 sources, idFactory, kernelSources,
657 positives=positives,
658 negatives=negatives,
659 measurementResults=measurementResults)
660 return measurementResults
661
662 def _prepareInputs(self, difference):
663 """Ensure that we start with an empty detection and deblended mask.
664
665 Parameters
666 ----------
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.
670
671 Raises
672 ------
673 lsst.pipe.base.UpstreamFailureNoWorkFound
674 If the PSF is not usable for measurement.
675 """
676 # Check that we have a valid PSF now before we do more work
677 sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius()
678 if np.isnan(sigma):
679 raise pipeBase.UpstreamFailureNoWorkFound("Invalid PSF detected! PSF width evaluates to NaN.")
680 # Ensure that we start with an empty detection and deblended mask.
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)
686
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.
690
691 Parameters
692 ----------
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
697 difference image.
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
704 difference image.
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.
715
716 Returns
717 -------
718 measurementResults : `lsst.pipe.base.Struct`
719
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.
724 """
725 if measurementResults is None:
726 measurementResults = pipeBase.Struct()
727 self.metadata["nUnmergedDiaSources"] = len(sources)
728 if self.config.doMerge:
729 # preserve peak schema, if there are any footprints
730 if len(positives) > 0:
731 peakSchema = positives[0].getFootprint().peaks.schema
732 elif len(negatives) > 0:
733 peakSchema = negatives[0].getFootprint().peaks.schema
734 else:
735 peakSchema = afwDetection.PeakTable.makeMinimalSchema()
736 mergeList = afwDetection.FootprintMergeList(self.schema,
737 ["positive", "negative"], peakSchema)
738 initialDiaSources = afwTable.SourceCatalog(self.schema)
739 # Start with positive, as FootprintMergeList will self-merge the
740 # subsequent added catalogs, and we want to try to preserve
741 # deblended positive sources.
742 mergeList.addCatalog(initialDiaSources.table, positives, "positive", minNewPeakDist=0)
743 mergeList.addCatalog(initialDiaSources.table, negatives, "negative", minNewPeakDist=0)
744 mergeList.getFinalSources(initialDiaSources)
745 # Flag as negative those sources that *only* came from the negative
746 # footprint set.
747 initialDiaSources["is_negative"] = initialDiaSources["merge_footprint_negative"] & \
748 ~initialDiaSources["merge_footprint_positive"]
749 self.log.info("Merging detections into %d sources", len(initialDiaSources))
750 else:
751 initialDiaSources = sources
752
753 # Assign source ids at the end: deblend/merge mean that we don't keep
754 # track of parents and children, we only care about the final ids.
755 for source in initialDiaSources:
756 source.setId(idFactory())
757 # Ensure sources added after this get correct ids.
758 initialDiaSources.getTable().setIdFactory(idFactory)
759 initialDiaSources.setMetadata(self.algMetadata)
760
761 self.metadata["nMergedDiaSources"] = len(initialDiaSources)
762
763 if self.config.doMaskStreaks:
764 streakInfo = self._runStreakMasking(difference)
765
766 if self.config.doSkySources:
767 self.addSkySources(initialDiaSources, difference.mask, difference.info.id)
768
769 if not initialDiaSources.isContiguous():
770 initialDiaSources = initialDiaSources.copy(deep=True)
771
772 self.measureDiaSources(initialDiaSources, science, difference, matchedTemplate)
773
774 # Remove unphysical diaSources per config.badSourceFlags
775 diaSources = self._removeBadSources(initialDiaSources)
776
777 if self.config.run_sattle:
778 diaSources = self.filterSatellites(diaSources, science)
779
780 # Flag diaSources in glint trails, but do not remove them
781 diaSources, trail_parameters = self._find_glint_trails(diaSources)
782 if self.config.writeGlintInfo:
783 measurementResults.mergeItems(trail_parameters, 'glintTrailInfo')
784
785 if self.config.doForcedMeasurement:
786 self.measureForcedSources(diaSources, science, difference.getWcs())
787 self.measureForcedSources(diaSources, matchedTemplate, difference.getWcs(),
788 template=True)
789
790 # Clear the image plane for regions with NO_DATA.
791 # These regions are most often caused by insufficient template coverage.
792 # Do this for the final difference image after detection and measurement
793 # since the subtasks should all be configured to handle NO_DATA properly
794 difference.image.array[difference.mask.array & difference.mask.getPlaneBitMask('NO_DATA') > 0] = 0
795
796 measurementResults.subtractedMeasuredExposure = difference
797
798 if self.config.doMaskStreaks and self.config.writeStreakInfo:
799 measurementResults.mergeItems(streakInfo, 'maskedStreaks')
800
801 self.calculateMetrics(science, difference, diaSources, kernelSources)
802
803 if np.count_nonzero(~diaSources["sky_source"]) > 0:
804 measurementResults.diaSources = diaSources
805 elif self.config.raiseOnNoDiaSources:
806 raise NoDiaSourcesError()
807 elif len(diaSources) > 0:
808 # This option allows returning sky sources,
809 # even if there are no diaSources
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"])
814 )
815 return measurementResults
816
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.
820
821 Parameters
822 ----------
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.
828
829 Returns
830 -------
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``.
836 """
837
838 def deblend(footprints, negative=False):
839 """Deblend a positive or negative footprint set,
840 and return the deblended children.
841
842 Parameters
843 ----------
844 footprints : `lsst.afw.detection.FootprintSet`
845 negative : `bool`
846 Set True if the footprints contain negative fluxes
847
848 Returns
849 -------
850 sources : `lsst.afw.table.SourceCatalog`
851 """
852 sources = afwTable.SourceCatalog(self.schema)
853 footprints.makeSources(sources)
854 if negative:
855 # Invert the image so the deblender can run on positive peaks
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]
860 # Set the heavy footprint pixel values back to reality
861 for child in children:
862 footprint = child.getFootprint()
863 array = footprint.getImageArray()
864 array *= -1
865 else:
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)
870 # Clear parents, so that measurement plugins behave correctly.
871 sources['parent'] = 0
872 return sources.copy(deep=True)
873
874 positives = deblend(positiveFootprints)
875 negatives = deblend(negativeFootprints, negative=True)
876
877 sources = afwTable.SourceCatalog(self.schema)
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
884
885 def _removeBadSources(self, diaSources):
886 """Remove unphysical diaSources from the catalog.
887
888 Parameters
889 ----------
890 diaSources : `lsst.afw.table.SourceCatalog`
891 The catalog of detected sources.
892
893 Returns
894 -------
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.
898 """
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)
903 if nBad > 0:
904 self.log.debug("Found %d unphysical sources with flag %s.", nBad, flag)
905 selector &= ~flags
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)
910
911 def _find_glint_trails(self, diaSources):
912 """Define a new flag column for diaSources that are in a glint trail.
913
914 Parameters
915 ----------
916 diaSources : `lsst.afw.table.SourceCatalog`
917 The catalog of detected sources.
918
919 Returns
920 -------
921 diaSources : `lsst.afw.table.SourceCatalog`
922 The updated catalog of detected sources, with a new bool column
923 called 'glint_trail' added.
924
925 trail_parameters : `dict`
926 Parameters of all the trails that were found.
927 """
928 if self.config.doSkySources:
929 # Do not include sky sources in glint detection
930 candidateDiaSources = diaSources[~diaSources["sky_source"]].copy(deep=True)
931 else:
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)
937
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,
944 'angles': angles}
945
946 trail_parameters = pipeBase.Struct(glintTrailInfo=parameters)
947
948 return diaSources, trail_parameters
949
950 def addSkySources(self, diaSources, mask, seed,
951 subtask=None):
952 """Add sources in empty regions of the difference image
953 for measuring the background.
954
955 Parameters
956 ----------
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.
961 seed : `int`
962 Seed value to initialize the random number generator.
963 """
964 if subtask is None:
965 subtask = self.skySources
966 if subtask.config.nSources <= 0:
967 self.metadata[f"n_{subtask.getName()}"] = 0
968 return
969 skySourceFootprints = subtask.run(mask=mask, seed=seed, catalog=diaSources)
970 self.metadata[f"n_{subtask.getName()}"] = len(skySourceFootprints)
971
972 def measureDiaSources(self, diaSources, science, difference, matchedTemplate):
973 """Use (matched) template and science image to constrain dipole fitting.
974
975 Parameters
976 ----------
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
985 difference image.
986 """
987 # Ensure that the required mask planes are present
988 for mp in self.config.measurement.plugins["base_PixelFlags"].masksFpAnywhere:
989 difference.mask.addMaskPlane(mp)
990 # Note that this may not be correct if we convolved the science image.
991 # In the future we may wish to persist the matchedScience image.
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.")
997 else:
998 self.applyApCorr.run(
999 catalog=diaSources,
1000 apCorrMap=apCorrMap,
1001 )
1002
1003 def measureForcedSources(self, diaSources, image, wcs, template=False):
1004 """Perform forced measurement of the diaSources on a direct image.
1005
1006 Parameters
1007 ----------
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.
1014 template : `bool`
1015 Is the forced measurement being performed on the template?
1016 """
1017 # Run forced psf photometry on the image at the diaSource locations.
1018 # Copy the measured flux and error into the diaSource.
1019 forcedSources = self.forcedMeasurement.generateMeasCat(image, diaSources, wcs)
1020 self.forcedMeasurement.run(forcedSources, image, diaSources, wcs)
1021
1022 if template:
1023 base_key = 'ip_diffim_forced_template_PsfFlux'
1024 else:
1025 base_key = 'ip_diffim_forced_PsfFlux'
1026 mapper = afwTable.SchemaMapper(forcedSources.schema, diaSources.schema)
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)
1041
1042 def calculateMetrics(self, science, difference, diaSources, kernelSources):
1043 """Add difference image QA metrics to the Task metadata.
1044
1045 This may be used to produce corresponding metrics (see
1046 lsst.analysis.tools.tasks.diffimTaskDetectorVisitMetricAnalysis).
1047
1048 Parameters
1049 ----------
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.
1058 """
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)
1067 detPosPix &= badPix
1068 detNegPix &= badPix
1069 self.metadata["nBadPixelsDetectedPositive"] = np.sum(detPosPix)
1070 self.metadata["nBadPixelsDetectedNegative"] = np.sum(detNegPix)
1071
1072 metricsMaskPlanes = list(mask.getMaskPlaneDict().keys())
1073 for maskPlane in metricsMaskPlanes:
1074 try:
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)
1079
1080 if self.config.doSkySources:
1081 skySources = diaSources[diaSources["sky_source"]]
1082 else:
1083 skySources = None
1084 metrics = computeDifferenceImageMetrics(science, difference, kernelSources, sky_sources=skySources)
1085
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:
1096 raise BadSubtractionError(ratio=metrics.differenceFootprintRatioMean,
1097 threshold=self.config.badSubtractionRatioThreshold)
1098 if metrics.differenceFootprintRatioStdev > self.config.badSubtractionVariationThreshold:
1099 raise BadSubtractionError(ratio=metrics.differenceFootprintRatioStdev,
1100 threshold=self.config.badSubtractionVariationThreshold)
1101
1102 def getSattleDiaSourceAllowlist(self, diaSources, science):
1103 """Query the sattle service and determine which diaSources are allowed.
1104
1105 Parameters
1106 ----------
1107 diaSources : `lsst.afw.table.SourceCatalog`
1108 The catalog of detected sources.
1109 science : `lsst.afw.image.ExposureF`
1110 Science exposure that was subtracted.
1111
1112 Returns
1113 ----------
1114 allow_list : `list` of `int`
1115 diaSourceIds of diaSources that can be made public.
1116
1117 Raises
1118 ------
1119 requests.HTTPError
1120 Raised if sattle call does not return success.
1121 """
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')
1126
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})
1133
1134 payload = {"visit_id": visit_id, "detector_id": science.getDetector().getId(),
1135 "diasources": dia_sources_json, "historical": self.config.sattle_historical}
1136
1137 sattle_output = requests.put(f'{sattle_uri_base}/diasource_allow_list',
1138 json=payload)
1139
1140 # retry once if visit cache is not populated
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)
1145
1146 sattle_output.raise_for_status()
1147
1148 return sattle_output.json()['allow_list']
1149
1150 def filterSatellites(self, diaSources, science):
1151 """Remove diaSources overlapping predicted satellite positions.
1152
1153 Parameters
1154 ----------
1155 diaSources : `lsst.afw.table.SourceCatalog`
1156 The catalog of detected sources.
1157 science : `lsst.afw.image.ExposureF`
1158 Science exposure that was subtracted.
1159
1160 Returns
1161 ----------
1162 filterdDiaSources : `lsst.afw.table.SourceCatalog`
1163 Filtered catalog of diaSources
1164 """
1165
1166 allow_list = self.getSattleDiaSourceAllowlist(diaSources, science)
1167
1168 if allow_list:
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)
1172 else:
1173 self.log.warning('Sattle allowlist is empty, all diaSources removed')
1174 diaSources = diaSources[0:0].copy(deep=True)
1175 return diaSources
1176
1177 def _runStreakMasking(self, difference):
1178 """Do streak masking and optionally save the resulting streak
1179 fit parameters in a catalog.
1180
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.
1186
1187 Parameters
1188 ----------
1189 difference: `lsst.afw.image.Exposure`
1190 The exposure in which to search for streaks. Must have a detection
1191 mask.
1192
1193 Returns
1194 -------
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.
1206 """
1207 maskedImage = difference.maskedImage
1208 # Bin the diffim to enhance low surface brightness streaks
1209 binnedMaskedImage = afwMath.binImage(maskedImage,
1210 self.config.streakBinFactor,
1211 self.config.streakBinFactor)
1212 binnedExposure = afwImage.ExposureF(binnedMaskedImage.getBBox())
1213 binnedExposure.setMaskedImage(binnedMaskedImage)
1214 # Clear the DETECTED mask plane before streak detection
1215 binnedExposure.mask &= ~binnedExposure.mask.getPlaneBitMask('DETECTED')
1216 # Rerun detection to set the DETECTED mask plane on binnedExposure
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,
1224 axis=1)
1225 # Create new version of a diffim with DETECTED based on binnedExposure
1226 streakMaskedImage = maskedImage.clone()
1227 ysize, xsize = rescaledDetectedMaskPlane.shape
1228 streakMaskedImage.mask.array[:ysize, :xsize] |= rescaledDetectedMaskPlane
1229 # Detect streaks on this new version of the diffim
1230 streaks = self.maskStreaks.run(streakMaskedImage)
1231 streakMaskPlane = streakMaskedImage.mask.array & streakMaskedImage.mask.getPlaneBitMask('STREAK')
1232 # Apply the new STREAK mask to the original diffim
1233 maskedImage.mask.array |= streakMaskPlane
1234
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}
1243 else:
1244 streakInfo = {'rho': np.array([]), 'theta': np.array([]), 'sigma': np.array([]),
1245 'reducedChi2': np.array([]), 'modelMaximum': np.array([])}
1246 return pipeBase.Struct(maskedStreaks=streakInfo)
1247
1248
1249class DetectAndMeasureScoreConnections(DetectAndMeasureConnections):
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",
1255 )
1256
1257
1258class DetectAndMeasureScoreConfig(DetectAndMeasureConfig,
1259 pipelineConnections=DetectAndMeasureScoreConnections):
1260 pass
1261
1262
1263class DetectAndMeasureScoreTask(DetectAndMeasureTask):
1264 """Detect DIA sources using a score image,
1265 and measure the detections on the difference image.
1266
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
1270 dipole diaSources.
1271 Sky sources, or forced detections in background regions, will optionally
1272 be added, and the configured measurement algorithm will be run on all
1273 detections.
1274 """
1275 ConfigClass = DetectAndMeasureScoreConfig
1276 _DefaultName = "detectAndMeasureScore"
1277
1278 @timeMethod
1279 def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources,
1280 idFactory=None):
1281 """Detect and measure sources on a score image.
1282
1283 Parameters
1284 ----------
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
1289 difference image.
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.
1300
1301 Returns
1302 -------
1303 measurementResults : `lsst.pipe.base.Struct`
1304
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.
1309 """
1310 if idFactory is None:
1311 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
1312
1313 self._prepareInputs(scoreExposure)
1314
1315 # Don't use the idFactory until after deblend+merge, so that we aren't
1316 # generating ids that just get thrown away (footprint merge doesn't
1317 # know about past ids).
1318 table = afwTable.SourceTable.make(self.schema)
1319 results = self.detection.run(
1320 table=table,
1321 exposure=scoreExposure,
1322 doSmooth=False,
1323 )
1324 # Copy the detection mask from the Score image to the difference image
1325 difference.mask.assign(scoreExposure.mask, scoreExposure.getBBox())
1326
1327 if self.config.doDeblend:
1328 sources, positives, negatives = self._deblend(difference,
1329 results.positive,
1330 results.negative)
1331
1332 return self.processResults(science, matchedTemplate, difference,
1333 sources, idFactory, kernelSources,
1334 positives=positives,
1335 negatives=negatives)
1336
1337 else:
1338 positives = afwTable.SourceCatalog(self.schema)
1339 results.positive.makeSources(positives)
1340 negatives = afwTable.SourceCatalog(self.schema)
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.
std::shared_ptr< ImageT > binImage(ImageT const &inImage, int const binX, int const binY, lsst::afw::math::Property const flags=lsst::afw::math::MEAN)
Definition binImage.cc:44
run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter)