Classes | |
| class | BadSubtractionError |
| class | NoDiaSourcesError |
| class | TooManyCosmicRays |
| class | DetectAndMeasureConnections |
Variable Documentation
◆ _TEMPLATE_EPOCH_NAMES
|
protected |
Initial value:
= {
"TEMPLATE MJD MIN": "templateMjdMin",
"TEMPLATE MJD MAX": "templateMjdMax",
"TEMPLATE MJD MEAN": "templateMjdMean",
}
Definition at line 56 of file detectAndMeasure.py.
◆ allow_list
| lsst.ip.diffim.detectAndMeasure.allow_list : `list` of `int` |
Definition at line 1569 of file detectAndMeasure.py.
◆ background
| lsst.ip.diffim.detectAndMeasure.background : `lsst.afw.math.BackgroundList` |
Definition at line 897 of file detectAndMeasure.py.
◆ diaSources
| lsst.ip.diffim.detectAndMeasure.diaSources : `lsst.afw.table.SourceCatalog` |
if measurementResults is None:
measurementResults = pipeBase.Struct()
self.metadata["nUnmergedDiaSources"] = len(sources)
if self.config.doMerge:
# preserve peak schema, if there are any footprints
if len(positives) > 0:
peakSchema = positives[0].getFootprint().peaks.schema
elif len(negatives) > 0:
peakSchema = negatives[0].getFootprint().peaks.schema
else:
peakSchema = afwDetection.PeakTable.makeMinimalSchema()
mergeList = afwDetection.FootprintMergeList(self.schema,
["positive", "negative"], peakSchema)
initialDiaSources = afwTable.SourceCatalog(self.schema)
# Start with positive, as FootprintMergeList will self-merge the
# subsequent added catalogs, and we want to try to preserve
# deblended positive sources.
mergeList.addCatalog(initialDiaSources.table, positives, "positive", minNewPeakDist=0)
mergeList.addCatalog(initialDiaSources.table, negatives, "negative", minNewPeakDist=0)
mergeList.getFinalSources(initialDiaSources)
# Flag as negative those sources that *only* came from the negative
# footprint set.
initialDiaSources["is_negative"] = initialDiaSources["merge_footprint_negative"] & \
~initialDiaSources["merge_footprint_positive"]
self.log.info("Merging detections into %d sources", len(initialDiaSources))
# All positive peaks were added before any negative peaks, so
# re-order the peaks in each footprint by significance.
self._reorderPeaksBySignificance(initialDiaSources)
else:
initialDiaSources = sources
# Assign source ids at the end: deblend/merge mean that we don't keep
# track of parents and children, we only care about the final ids.
for source in initialDiaSources:
source.setId(idFactory())
# Ensure sources added after this get correct ids.
initialDiaSources.getTable().setIdFactory(idFactory)
initialDiaSources.setMetadata(self.algMetadata)
self.metadata["nMergedDiaSources"] = len(initialDiaSources)
if self.config.doMaskStreaks:
streakInfo = self._runStreakMasking(difference)
if self.config.doSkySources:
self.addSkySources(initialDiaSources, difference.mask, difference.info.id)
if self.config.doRejectBadMaskPlaneDetections:
# Save time by rejecting peaks on bad mask planes prior to
# measurement.
initialDiaSources = self._rejectBadMaskedDetections(initialDiaSources, difference.mask)
if not initialDiaSources.isContiguous():
initialDiaSources = initialDiaSources.copy(deep=True)
self.measureDiaSources(initialDiaSources, science, difference, matchedTemplate)
# Remove unphysical diaSources per config.badSourceFlags
diaSources = self._removeBadSources(initialDiaSources)
if self.config.run_sattle:
diaSources = self.filterSatellites(diaSources, science)
# Flag diaSources in glint trails, but do not remove them
diaSources, trail_parameters = self._find_glint_trails(diaSources)
if self.config.writeGlintInfo:
measurementResults.mergeItems(trail_parameters, 'glintTrailInfo')
self._measureTemplateEpochs(diaSources, templateInfo)
measurementResults.templateEpochs = self._summarizeTemplateEpochs(diaSources)
if self.config.doForcedMeasurement:
self.measureForcedSources(diaSources, science, difference.getWcs())
self.measureForcedSources(diaSources, matchedTemplate, difference.getWcs(),
template=True)
# Clear the image plane for regions with NO_DATA.
# These regions are most often caused by insufficient template coverage.
# Do this for the final difference image after detection and measurement
# since the subtasks should all be configured to handle NO_DATA properly
difference.image.array[difference.mask.array & difference.mask.getPlaneBitMask('NO_DATA') > 0] = 0
measurementResults.subtractedMeasuredExposure = difference
if self.config.doMaskStreaks and self.config.writeStreakInfo:
measurementResults.maskedStreaks = streakInfo.maskedStreaks
if kernelSources is not None:
self.calculateMetrics(science, difference, diaSources, kernelSources)
if np.count_nonzero(~diaSources["sky_source"]) > 0:
measurementResults.diaSources = diaSources
elif self.config.raiseOnNoDiaSources:
raise NoDiaSourcesError()
elif len(diaSources) > 0:
# This option allows returning sky sources,
# even if there are no diaSources
measurementResults.diaSources = diaSources
self.log.info("Measured %d diaSources and %d sky sources",
np.count_nonzero(~diaSources["sky_source"]),
np.count_nonzero(diaSources["sky_source"])
)
return measurementResults
def _reorderPeaksBySignificance(self, diaSources):
x = diaSources.getX()
y = diaSources.getY()
numerator = np.zeros(len(diaSources))
denominator = np.zeros(len(diaSources))
epochMin = np.full(len(diaSources), np.inf)
epochMax = np.full(len(diaSources), -np.inf)
for info in templateInfo or []:
inCoadd = info.bounds.contains(x=x, y=y)
if not inCoadd.any():
continue
weightSum = np.zeros(len(diaSources))
weightedSum = np.zeros(len(diaSources))
for record in info.inputs:
covers = inCoadd & record.bounds.contains(x=x, y=y)
weightSum[covers] += record.weight
weightedSum[covers] += record.weight*record.mid_mjd_tai
epochMin[covers] = np.minimum(epochMin[covers], record.mid_mjd_tai)
epochMax[covers] = np.maximum(epochMax[covers], record.mid_mjd_tai)
inputs = weightSum > 0
numerator[inputs] += weightedSum[inputs]/np.sqrt(weightSum[inputs])
denominator[inputs] += np.sqrt(weightSum[inputs])
known = denominator > 0
epoch = np.full(len(diaSources), np.nan)
epoch[known] = numerator[known]/denominator[known]
epochMin[~known] = np.nan
epochMax[~known] = np.nan
diaSources["templateMjdMean"] = epoch
diaSources["templateMjdMin"] = epochMin
diaSources["templateMjdMax"] = epochMax
def _summarizeTemplateEpochs(self, diaSources):
known = np.isfinite(diaSources["templateMjdMean"])
if not known.any():
for name in _TEMPLATE_EPOCH_NAMES.values():
self.metadata[name] = np.nan
return {}
values = {"templateMjdMin": float(np.min(diaSources["templateMjdMin"][known])),
"templateMjdMax": float(np.max(diaSources["templateMjdMax"][known])),
"templateMjdMean": float(np.mean(diaSources["templateMjdMean"][known]))}
for name, value in values.items():
self.metadata[name] = value
return {key: values[name] for key, name in _TEMPLATE_EPOCH_NAMES.items()}
def addSkySources(self, diaSources, mask, seed,
subtask=None):
if subtask is None:
subtask = self.skySources
if subtask.config.nSources <= 0:
self.metadata[f"n_{subtask.getName()}"] = 0
return
skySourceFootprints = subtask.run(mask=mask, seed=seed, catalog=diaSources)
self.metadata[f"n_{subtask.getName()}"] = len(skySourceFootprints)
def measureDiaSources(self, diaSources, science, difference, matchedTemplate):
# Ensure that the required mask planes are present
for mp in self.config.measurement.plugins["base_PixelFlags"].masksFpAnywhere:
difference.mask.addMaskPlane(mp)
# Note that this may not be correct if we convolved the science image.
# In the future we may wish to persist the matchedScience image.
self.measurement.run(diaSources, difference, science, matchedTemplate)
if self.config.doApCorr:
apCorrMap = difference.getInfo().getApCorrMap()
if apCorrMap is None:
self.log.warning("Difference image does not have valid aperture correction; skipping.")
else:
self.applyApCorr.run(
catalog=diaSources,
apCorrMap=apCorrMap,
)
def measureForcedSources(self, diaSources, image, wcs, template=False):
mask = difference.mask
badPix = (mask.array & mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) > 0
self.metadata["nGoodPixels"] = np.sum(~badPix)
self.metadata["nBadPixels"] = np.sum(badPix)
detPosPix = (mask.array & mask.getPlaneBitMask("DETECTED")) > 0
detNegPix = (mask.array & mask.getPlaneBitMask("DETECTED_NEGATIVE")) > 0
self.metadata["nPixelsDetectedPositive"] = np.sum(detPosPix)
self.metadata["nPixelsDetectedNegative"] = np.sum(detNegPix)
detPosPix &= badPix
detNegPix &= badPix
self.metadata["nBadPixelsDetectedPositive"] = np.sum(detPosPix)
self.metadata["nBadPixelsDetectedNegative"] = np.sum(detNegPix)
metricsMaskPlanes = list(mask.getMaskPlaneDict().keys())
for maskPlane in metricsMaskPlanes:
try:
self.metadata["%s_mask_fraction"%maskPlane.lower()] = evaluateMaskFraction(mask, maskPlane)
except InvalidParameterError:
self.metadata["%s_mask_fraction"%maskPlane.lower()] = -1
self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane)
if self.config.doSkySources:
skySources = diaSources[diaSources["sky_source"]]
else:
skySources = None
metrics = computeDifferenceImageMetrics(science, difference, kernelSources, sky_sources=skySources)
self.metadata["residualFootprintRatioMean"] = metrics.differenceFootprintRatioMean
self.metadata["residualFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev
self.metadata["differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean
self.metadata["differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev
self.log.info("Mean, stdev of ratio of difference to science "
"pixels in star footprints: %5.4f, %5.4f",
self.metadata["residualFootprintRatioMean"],
self.metadata["residualFootprintRatioStdev"])
if self.config.raiseOnBadSubtractionRatio:
if metrics.differenceFootprintRatioMean > self.config.badSubtractionRatioThreshold:
raise BadSubtractionError(ratio=metrics.differenceFootprintRatioMean,
threshold=self.config.badSubtractionRatioThreshold)
if metrics.differenceFootprintRatioStdev > self.config.badSubtractionVariationThreshold:
raise BadSubtractionError(ratio=metrics.differenceFootprintRatioStdev,
threshold=self.config.badSubtractionVariationThreshold)
def getSattleDiaSourceAllowlist(self, diaSources, science):
wcs = science.getWcs()
visit_info = science.getInfo().getVisitInfo()
visit_id = visit_info.getId()
sattle_uri_base = os.getenv('SATTLE_URI_BASE')
dia_sources_json = []
for source in diaSources:
source_bbox = source.getFootprint().getBBox()
corners = wcs.pixelToSky([lsst.geom.Point2D(c) for c in source_bbox.getCorners()])
bbox_radec = [[pt.getRa().asDegrees(), pt.getDec().asDegrees()] for pt in corners]
dia_sources_json.append({"diasource_id": source["id"], "bbox": bbox_radec})
payload = {"visit_id": visit_id, "detector_id": science.getDetector().getId(),
"diasources": dia_sources_json, "historical": self.config.sattle_historical}
sattle_output = requests.put(f'{sattle_uri_base}/diasource_allow_list',
json=payload)
# retry once if visit cache is not populated
if sattle_output.status_code == 404:
self.log.warning(f'Visit {visit_id} not found in sattle cache, re-sending')
populate_sattle_visit_cache(visit_info, historical=self.config.sattle_historical)
sattle_output = requests.put(f'{sattle_uri_base}/diasource_allow_list', json=payload)
sattle_output.raise_for_status()
return sattle_output.json()['allow_list']
@timeMethod
def filterSatellites(self, diaSources, science):
Definition at line 1351 of file detectAndMeasure.py.
◆ difference
| lsst.ip.diffim.detectAndMeasure.difference : `lsst.afw.image.ExposureF` |
if measurementResults is None:
measurementResults = pipeBase.Struct()
if idFactory is None:
idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
# Check image properties and clear detection mask planes.
self._prepareInputs(difference)
if self.config.doSubtractBackground:
background = self._fitAndSubtractBackground(differenceExposure=difference)
else:
background = afwMath.BackgroundList()
if self.config.doFindCosmicRays:
self.findAndMaskCosmicRays(difference)
# Don't use the idFactory until after deblend+merge, so that we aren't
# generating ids that just get thrown away (footprint merge doesn't
# know about past ids).
table = afwTable.SourceTable.make(self.schema)
results = self.detection.run(
table=table,
exposure=difference,
doSmooth=True,
background=background,
clearMask=True,
)
measurementResults.differenceBackground = background
if self.config.doDeblend:
sources, positives, negatives = self._deblend(difference,
results.positive,
results.negative)
else:
positives = afwTable.SourceCatalog(self.schema)
results.positive.makeSources(positives)
negatives = afwTable.SourceCatalog(self.schema)
results.negative.makeSources(negatives)
sources = results.sources
self.processResults(science, matchedTemplate, difference,
sources, idFactory, kernelSources,
positives=positives,
negatives=negatives,
measurementResults=measurementResults,
templateInfo=templateInfo)
return measurementResults
def _prepareInputs(self, difference):
# Run forced psf photometry on the image at the diaSource locations.
# Copy the measured flux and error into the diaSource.
forcedSources = self.forcedMeasurement.generateMeasCat(image, diaSources, wcs)
self.forcedMeasurement.run(forcedSources, image, diaSources, wcs)
if template:
base_key = 'ip_diffim_forced_template_PsfFlux'
else:
base_key = 'ip_diffim_forced_PsfFlux'
mapper = afwTable.SchemaMapper(forcedSources.schema, diaSources.schema)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFlux")[0],
f"{base_key}_instFlux", True)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFluxErr")[0],
f"{base_key}_instFluxErr", True)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_area")[0],
f"{base_key}_area", True)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag")[0],
f"{base_key}_flag", True)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_noGoodPixels")[0],
f"{base_key}_flag_noGoodPixels", True)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_edge")[0],
f"{base_key}_flag_edge", True)
for diaSource, forcedSource in zip(diaSources, forcedSources):
diaSource.assign(forcedSource, mapper)
def findAndMaskCosmicRays(self, difference):
allow_list = self.getSattleDiaSourceAllowlist(diaSources, science)
if allow_list:
allow_set = set(allow_list)
allowed_ids = [source['id'] in allow_set for source in diaSources]
diaSources = diaSources[np.array(allowed_ids)].copy(deep=True)
else:
self.log.warning('Sattle allowlist is empty, all diaSources removed')
diaSources = diaSources[0:0].copy(deep=True)
return diaSources
def _runStreakMasking(self, difference):
Definition at line 763 of file detectAndMeasure.py.
◆ differenceExposure
| lsst.ip.diffim.detectAndMeasure.differenceExposure : `lsst.afw.image.ExposureF` |
# Check that we have a valid PSF now before we do more work
sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius()
if np.isnan(sigma):
raise pipeBase.UpstreamFailureNoWorkFound("Invalid PSF detected! PSF width evaluates to NaN.")
# Ensure that we start with an empty detection and deblended mask.
mask = difference.mask
for mp in self.config.clearMaskPlanes:
if mp not in mask.getMaskPlaneDict():
mask.addMaskPlane(mp)
mask &= ~mask.getPlaneBitMask(self.config.clearMaskPlanes)
def _fitAndSubtractBackground(self, differenceExposure, scoreExposure=None):
Definition at line 888 of file detectAndMeasure.py.
◆ filterdDiaSources
| lsst.ip.diffim.detectAndMeasure.filterdDiaSources : `lsst.afw.table.SourceCatalog` |
Definition at line 1618 of file detectAndMeasure.py.
◆ idFactory
| lsst.ip.diffim.detectAndMeasure.idFactory : `lsst.afw.table.IdFactory`, optional |
Definition at line 767 of file detectAndMeasure.py.
◆ image
| lsst.ip.diffim.detectAndMeasure.image : `lsst.afw.image.ExposureF` |
Definition at line 1432 of file detectAndMeasure.py.
◆ inputImages
| lsst.ip.diffim.detectAndMeasure.inputImages : `dict` [`str`, `lsst.images.DifferenceImage`] |
Definition at line 704 of file detectAndMeasure.py.
◆ kernelSources
| lsst.ip.diffim.detectAndMeasure.kernelSources : `lsst.afw.table.SourceCatalog`, optional |
Definition at line 765 of file detectAndMeasure.py.
◆ mask
| lsst.ip.diffim.detectAndMeasure.mask : `lsst.afw.image.Mask` |
Definition at line 1381 of file detectAndMeasure.py.
◆ matchedTemplate
| lsst.ip.diffim.detectAndMeasure.matchedTemplate : `lsst.afw.image.ExposureF` |
Definition at line 760 of file detectAndMeasure.py.
◆ measurementResults
| lsst.ip.diffim.detectAndMeasure.measurementResults : `lsst.pipe.base.Struct`, optional |
Definition at line 771 of file detectAndMeasure.py.
◆ negatives
| lsst.ip.diffim.detectAndMeasure.negatives : `lsst.afw.table.SourceCatalog`, optional |
Definition at line 961 of file detectAndMeasure.py.
◆ positives
| lsst.ip.diffim.detectAndMeasure.positives : `lsst.afw.table.SourceCatalog`, optional |
Definition at line 959 of file detectAndMeasure.py.
◆ results
| lsst.ip.diffim.detectAndMeasure.results : `lsst.pipe.base.Struct` |
doMerge = pexConfig.Field(
dtype=bool,
default=True,
doc="Merge positive and negative diaSources with grow radius "
"set by growFootprint"
)
doForcedMeasurement = pexConfig.Field(
dtype=bool,
default=True,
doc="Force photometer diaSource locations on PVI?"
)
doAddMetrics = pexConfig.Field(
dtype=bool,
default=False,
doc="Add columns to the source table to hold analysis metrics?"
)
doSubtractBackground = pexConfig.Field(
dtype=bool,
doc="Subtract a background model from the image before detection?",
default=True,
)
doWriteBackground = pexConfig.Field(
dtype=bool,
doc="Persist the fitted background model? Ignored when"
" image_type='future', which records it on the output images.",
default=False,
deprecated="Use image_type='future', which records the background on the output images."
" Will be removed after v31.",
)
doCalculateResidualMetics = pexConfig.Field(
dtype=bool,
doc="Calculate metrics to assess image subtraction quality for the task"
"metadata?",
default=True,
)
subtractInitialBackground = pexConfig.ConfigurableField(
target=lsst.meas.algorithms.SubtractBackgroundTask,
doc="Task to perform intial background subtraction, before first detection pass.",
)
subtractFinalBackground = pexConfig.ConfigurableField(
target=lsst.meas.algorithms.SubtractBackgroundTask,
doc="Task to perform final background subtraction, after first detection pass.",
)
doFindCosmicRays = pexConfig.Field(
dtype=bool,
doc="Detect and mask cosmic rays on the difference image?"
"CRs can be interpolated over by setting cosmicray.keepCRs=False",
default=True,
)
cosmicray = pexConfig.ConfigField(
dtype=FindCosmicRaysConfig,
doc="Options for finding and masking cosmic rays",
)
detection = pexConfig.ConfigurableField(
target=SourceDetectionTask,
doc="Final source detection for diaSource measurement",
)
streakDetection = pexConfig.ConfigurableField(
target=SourceDetectionTask,
doc="Separate source detection used only for streak masking",
)
doDeblend = pexConfig.Field(
dtype=bool,
default=False,
doc="Deblend DIASources after detection?"
)
deblend = pexConfig.ConfigurableField(
target=lsst.meas.deblender.SourceDeblendTask,
doc="Task to split blended sources into their components."
)
measurement = pexConfig.ConfigurableField(
target=DipoleFitTask,
doc="Task to measure sources on the difference image.",
)
doApCorr = lsst.pex.config.Field(
dtype=bool,
default=True,
doc="Run subtask to apply aperture corrections"
)
applyApCorr = lsst.pex.config.ConfigurableField(
target=ApplyApCorrTask,
doc="Task to apply aperture corrections"
)
forcedMeasurement = pexConfig.ConfigurableField(
target=ForcedMeasurementTask,
doc="Task to force photometer science image at diaSource locations.",
)
growFootprint = pexConfig.Field(
dtype=int,
default=2,
doc="Grow positive and negative footprints by this many pixels before merging"
)
diaSourceMatchRadius = pexConfig.Field(
dtype=float,
default=0.5,
doc="Match radius (in arcseconds) for DiaSource to Source association"
)
doSkySources = pexConfig.Field(
dtype=bool,
default=False,
doc="Generate sky sources?",
)
skySources = pexConfig.ConfigurableField(
target=SkyObjectsTask,
doc="Generate sky sources",
)
doMaskStreaks = pexConfig.Field(
dtype=bool,
default=True,
doc="Turn on streak masking",
)
maskStreaks = pexConfig.ConfigurableField(
target=MaskStreaksTask,
doc="Subtask for masking streaks. Only used if doMaskStreaks is True. "
"Adds a mask plane to an exposure, with the mask plane name set by streakMaskName.",
)
streakBinFactor = pexConfig.Field(
dtype=int,
default=4,
doc="Bin scale factor to use when rerunning detection for masking streaks. "
"Only used if doMaskStreaks is True.",
)
writeStreakInfo = pexConfig.Field(
dtype=bool,
default=False,
doc="Record the parameters of any detected streaks. For LSST, this should be turned off except for "
"development work."
)
findGlints = pexConfig.ConfigurableField(
target=FindGlintTrailsTask,
doc="Subtask for finding glint trails, usually caused by satellites or debris."
)
writeGlintInfo = pexConfig.Field(
dtype=bool,
default=True,
doc="Record the parameters of any detected glint trails."
)
setPrimaryFlags = pexConfig.ConfigurableField(
target=SetPrimaryFlagsTask,
doc="Task to add isPrimary and deblending-related flags to the catalog."
)
badSourceFlags = lsst.pex.config.ListField(
dtype=str,
doc="Sources with any of these flags set are removed before writing the output catalog.",
default=("base_PixelFlags_flag_offimage",
"base_PixelFlags_flag_interpolatedCenterAll",
"base_PixelFlags_flag_badCenter",
"base_PixelFlags_flag_edgeCenter",
"base_PixelFlags_flag_nodataCenter",
"base_PixelFlags_flag_saturatedCenter",
"base_PixelFlags_flag_sat_templateCenter",
"base_PixelFlags_flag_spikeCenter",
),
)
clearMaskPlanes = lsst.pex.config.ListField(
dtype=str,
doc="Mask planes to clear before running detection.",
default=("DETECTED", "DETECTED_NEGATIVE", "NOT_DEBLENDED", "STREAK"),
)
doRejectBadMaskPlaneDetections = pexConfig.Field(
dtype=bool,
default=True,
doc="Reject any peaks detected on ``badMaskPlanes`` before measurement."
"These should all be rejected downstream by ``badSourceFlags``, "
"but filtering earlier can save time."
)
badMaskPlanes = lsst.pex.config.ListField(
dtype=str,
doc="Detections whose footprint peak lies on a pixel with any of these"
" mask planes set will be rejected before measurement."
" Any missing mask planes will be silently ignored.",
default=("NO_DATA", "BAD", "SAT", "EDGE"),
)
raiseOnBadSubtractionRatio = pexConfig.Field(
dtype=bool,
default=True,
doc="Raise an error if the ratio of power in detected footprints"
" on the difference image to the power in footprints on the science"
" image exceeds ``badSubtractionRatioThreshold``",
)
badSubtractionRatioThreshold = pexConfig.Field(
dtype=float,
default=0.2,
doc="Maximum ratio of power in footprints on the difference image to"
" the same footprints on the science image."
"Only used if ``raiseOnBadSubtractionRatio`` is set",
)
badSubtractionVariationThreshold = pexConfig.Field(
dtype=float,
default=0.4,
doc="Maximum standard deviation of the ratio of power in footprints on"
" the difference image to the same footprints on the science image."
"Only used if ``raiseOnBadSubtractionRatio`` is set",
)
raiseOnNoDiaSources = pexConfig.Field(
dtype=bool,
default=True,
doc="Raise an algorithm error if no diaSources are detected.",
)
run_sattle = pexConfig.Field(
dtype=bool,
default=False,
doc="If true, dia source bounding boxes will be sent for verification"
"to the sattle service."
)
sattle_historical = pexConfig.Field(
dtype=bool,
default=False,
doc="If re-running a pipeline that requires sattle, this should be set "
"to True. This will populate sattle's cache with the historic data "
"closest in time to the exposure."
)
idGenerator = DetectorVisitIdGeneratorConfig.make_field()
image_type = pexConfig.ChoiceField[str](
"Image type of the input and measured difference images (and, for"
" DetectAndMeasureScoreTask, the input and measured score images).",
allowed={
"legacy": "Read and write lsst.afw.image.ExposureF.",
"future": "Read and write lsst.images.DifferenceImage.",
},
optional=False,
default="legacy",
)
def setDefaults(self):
# Background subtraction
# Use a small binsize for the first pass to reduce detections on glints
# and extended structures. Should not affect the detectability of
# faint diaSources
self.subtractInitialBackground.binSize = 8
self.subtractInitialBackground.useApprox = False
self.subtractInitialBackground.statisticsProperty = "MEDIAN"
self.subtractInitialBackground.doFilterSuperPixels = True
self.subtractInitialBackground.ignoredPixelMask = ["BAD",
"EDGE",
"DETECTED",
"DETECTED_NEGATIVE",
"NO_DATA",
]
# Use a larger binsize for the final background subtraction, to reduce
# over-subtraction of bright objects.
self.subtractFinalBackground.binSize = 40
self.subtractFinalBackground.useApprox = False
self.subtractFinalBackground.statisticsProperty = "MEDIAN"
self.subtractFinalBackground.doFilterSuperPixels = True
self.subtractFinalBackground.ignoredPixelMask = ["BAD",
"EDGE",
"DETECTED",
"DETECTED_NEGATIVE",
"NO_DATA",
]
# Cosmic ray detection
self.cosmicray.keepCRs = True # do not interpolate over detected CRs
# DiaSource Detection
self.detection.thresholdPolarity = "both"
self.detection.thresholdValue = 5.0
self.detection.reEstimateBackground = False
self.detection.thresholdType = "pixel_stdev"
self.detection.excludeMaskPlanes = []
# Copy configs for binned streak detection from the base detection task
self.streakDetection.thresholdType = self.detection.thresholdType
self.streakDetection.reEstimateBackground = False
self.streakDetection.excludeMaskPlanes = self.detection.excludeMaskPlanes
self.streakDetection.thresholdValue = self.detection.thresholdValue
# Only detect positive streaks
self.streakDetection.thresholdPolarity = "positive"
# Do not grow detected mask for streaks
self.streakDetection.nSigmaToGrow = 0
# Set the streak mask along the entire fit line, not only where the
# detected mask is set.
self.maskStreaks.onlyMaskDetected = False
# Restrict streak masking from growing too large
self.maskStreaks.maxStreakWidth = 100
# Restrict the number of iterations allowed for fitting streaks
# When the fit is good it should solve quickly, and exit a bad fit quickly
self.maskStreaks.maxFitIter = 10
# Only mask to 2 sigma in width
self.maskStreaks.nSigmaMask = 2
# Threshold for including streaks after the Hough Transform.
# A lower value will detect more features that are less linear.
self.maskStreaks.absMinimumKernelHeight = 2
self.measurement.plugins.names |= ["ext_trailedSources_Naive",
"base_LocalPhotoCalib",
"base_LocalWcs",
"ext_shapeHSM_HsmSourceMoments",
"ext_shapeHSM_HsmPsfMoments",
"base_ClassificationSizeExtendedness",
]
self.measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments"
self.measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
self.measurement.plugins["base_SdssCentroid"].maxDistToPeak = 5.0
self.forcedMeasurement.plugins = ["base_TransformedCentroid", "base_PsfFlux"]
self.forcedMeasurement.copyColumns = {
"id": "objectId", "parent": "parentObjectId", "coord_ra": "coord_ra", "coord_dec": "coord_dec"}
self.forcedMeasurement.slots.centroid = "base_TransformedCentroid"
self.forcedMeasurement.slots.shape = None
# Keep track of which footprints contain streaks
self.measurement.plugins["base_PixelFlags"].masksFpAnywhere = [
"STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SAT_TEMPLATE", "SPIKE"]
self.measurement.plugins["base_PixelFlags"].masksFpCenter = [
"STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SAT_TEMPLATE", "SPIKE"]
self.skySources.avoidMask = ["DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA", "EDGE"]
def validate(self):
super().validate()
if self.run_sattle:
if not os.getenv("SATTLE_URI_BASE"):
raise pexConfig.FieldValidationError(DetectAndMeasureConfig.run_sattle, self,
"Sattle requested but SATTLE_URI_BASE "
"environment variable not set.")
class DetectAndMeasureTask(lsst.pipe.base.PipelineTask):
ConfigClass = DetectAndMeasureConfig
_DefaultName = "detectAndMeasure"
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.schema = afwTable.SourceTable.makeMinimalSchema()
self.algMetadata = dafBase.PropertyList()
if self.config.doSubtractBackground:
self.makeSubtask("subtractInitialBackground")
self.makeSubtask("subtractFinalBackground")
self.makeSubtask("detection", schema=self.schema)
if self.config.doDeblend:
self.makeSubtask("deblend", schema=self.schema)
self.makeSubtask("setPrimaryFlags", schema=self.schema, isSingleFrame=True)
self.makeSubtask("measurement", schema=self.schema,
algMetadata=self.algMetadata)
if self.config.doApCorr:
self.makeSubtask("applyApCorr", schema=self.measurement.schema)
if self.config.doForcedMeasurement:
self.schema.addField(
"ip_diffim_forced_PsfFlux_instFlux", "D",
"Forced PSF flux measured on the direct image.",
units="count")
self.schema.addField(
"ip_diffim_forced_PsfFlux_instFluxErr", "D",
"Forced PSF flux error measured on the direct image.",
units="count")
self.schema.addField(
"ip_diffim_forced_PsfFlux_area", "F",
"Forced PSF flux effective area of PSF.",
units="pixel")
self.schema.addField(
"ip_diffim_forced_PsfFlux_flag", "Flag",
"Forced PSF flux general failure flag.")
self.schema.addField(
"ip_diffim_forced_PsfFlux_flag_noGoodPixels", "Flag",
"Forced PSF flux not enough non-rejected pixels in data to attempt the fit.")
self.schema.addField(
"ip_diffim_forced_PsfFlux_flag_edge", "Flag",
"Forced PSF flux object was too close to the edge of the image to use the full PSF model.")
self.schema.addField(
"ip_diffim_forced_template_PsfFlux_instFlux", "D",
"Forced PSF flux measured on the template image.",
units="count")
self.schema.addField(
"ip_diffim_forced_template_PsfFlux_instFluxErr", "D",
"Forced PSF flux error measured on the template image.",
units="count")
self.schema.addField(
"ip_diffim_forced_template_PsfFlux_area", "F",
"Forced template PSF flux effective area of PSF.",
units="pixel")
self.schema.addField(
"ip_diffim_forced_template_PsfFlux_flag", "Flag",
"Forced template PSF flux general failure flag.")
self.schema.addField(
"ip_diffim_forced_template_PsfFlux_flag_noGoodPixels", "Flag",
"Forced template PSF flux not enough non-rejected pixels in data to attempt the fit.")
self.schema.addField(
"ip_diffim_forced_template_PsfFlux_flag_edge", "Flag",
)
self.makeSubtask("forcedMeasurement", refSchema=self.schema)
self.schema.addField("refMatchId", "L", "unique id of reference catalog match")
self.schema.addField("srcMatchId", "L", "unique id of source match")
# Create the sky source task for use by metrics,
# even if sky sources are not added to the diaSource catalog
self.makeSubtask("skySources", schema=self.schema)
if self.config.doMaskStreaks:
self.makeSubtask("maskStreaks")
self.makeSubtask("streakDetection")
self.makeSubtask("findGlints")
self.schema.addField("glint_trail", "Flag", "DiaSource is part of a glint trail.")
self.schema.addField("reliability", type="F", doc="Reliability score of the DiaSource")
self.schema.addField("reliabilityVersion", type=str, size=7, doc="Version of the reliability model")
self.schema.addField("templateMjdMean", type="D", units="d",
doc="Weighted mean time (MJD, TAI) of the template inputs that"
" cover the diaSource centroid.")
self.schema.addField("templateMjdMin", type="D", units="d",
doc="Earliest time (MJD, TAI) of the template inputs that"
" cover the diaSource centroid.")
self.schema.addField("templateMjdMax", type="D", units="d",
doc="Latest time (MJD, TAI) of the template inputs that"
" cover the diaSource centroid.")
# To get the "merge_*" fields in the schema; have to re-initialize
# this later, once we have a peak schema post-detection.
lsst.afw.detection.FootprintMergeList(self.schema, ["positive", "negative"])
# Check that the schema and config are consistent
for flag in self.config.badSourceFlags:
if flag not in self.schema:
raise pipeBase.InvalidQuantumError("Field %s not in schema" % flag)
# initialize InitOutputs
self.outputSchema = afwTable.SourceCatalog(self.schema)
self.outputSchema.getTable().setMetadata(self.algMetadata)
def runQuantum(self, butlerQC: pipeBase.QuantumContext,
inputRefs: pipeBase.InputQuantizedConnection,
outputRefs: pipeBase.OutputQuantizedConnection):
inputs = butlerQC.get(inputRefs)
idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId)
idFactory = idGenerator.make_table_id_factory()
# Input image that each output is measured from, keyed by output name.
inputImages = {}
templateInfo = None
if self.config.image_type == "future":
dataId = butlerQC.quantum.dataId
for output, connection in (("subtractedMeasuredExposure", "difference"),
("scoreMeasuredExposure", "scoreExposure")):
if connection not in inputs:
continue
image = inputs[connection]
# The output is written back into this image, so it must
# already carry the kernel and the template records.
get_difference_image_provenance(image, "kernel", connection, dataId)
imageTemplates = get_difference_image_provenance(image, "templates", connection, dataId)
if connection == "difference":
templateInfo = imageTemplates
inputImages[output] = image
inputs[connection] = image.to_legacy()
# Specify the fields that `annotate` and `convert_outputs_to_future`
# need below, to ensure they exist, even as None.
measurementResults = pipeBase.Struct(
subtractedMeasuredExposure=None,
scoreMeasuredExposure=None,
diaSources=None,
maskedStreaks=None,
differenceBackground=None,
templateEpochs=None,
)
try:
self.run(**inputs, idFactory=idFactory, measurementResults=measurementResults,
templateInfo=templateInfo)
except pipeBase.AlgorithmError as e:
# Annotate before converting
error = pipeBase.AnnotatedPartialOutputsError.annotate(
e,
self,
measurementResults.subtractedMeasuredExposure,
measurementResults.scoreMeasuredExposure,
measurementResults.diaSources,
measurementResults.maskedStreaks,
log=self.log
)
if self.config.image_type == "future":
self.convert_outputs_to_future(measurementResults, inputImages)
butlerQC.put(measurementResults, outputRefs)
raise error from e
if self.config.image_type == "future":
self.convert_outputs_to_future(measurementResults, inputImages)
butlerQC.put(measurementResults, outputRefs)
def convert_outputs_to_future(self, results, inputImages):
Definition at line 702 of file detectAndMeasure.py.
◆ science
| lsst.ip.diffim.detectAndMeasure.science : `lsst.afw.image.ExposureF` |
background = getattr(results, "differenceBackground", None)
templateEpochs = getattr(results, "templateEpochs", None) or {}
for name, image in inputImages.items():
exposure = getattr(results, name, None)
if exposure is None:
continue
# These are no-ops when the legacy output is a view of the input,
# as returned by `~lsst.images.DifferenceImage.to_legacy`.
image.image.array[...] = exposure.image.array
image.variance.array[...] = exposure.variance.array
image.mask = Mask.from_legacy(exposure.mask, get_legacy_difference_image_mask_planes(),
sky_projection=image.sky_projection)
image.metadata.update(templateEpochs)
if background is not None and len(background) > 0:
# The same background is subtracted from the difference image
# and the score image, so it applies to both.
image.backgrounds.add(
"subtracted",
field_from_legacy_background(background, bounds=image.bbox, unit=image.unit),
description="Background subtracted from this image before detection.",
is_subtracted=True,
)
setattr(results, name, image)
if background is not None and len(background) == 0:
# An empty BackgroundList cannot be converted to a field, and means
# that no background was subtracted.
self.log.debug("No background model to attach to the output images.")
@timeMethod
def run(self, science, matchedTemplate, difference, kernelSources=None,
idFactory=None, measurementResults=None, templateInfo=None):
if scoreExposure is None:
detectionExposure = differenceExposure.clone()
background = self.subtractInitialBackground.run(detectionExposure).background
doSmooth = True
else:
# Use a clone of differenceExposure because the background is
# subtracted in place.
background = self.subtractInitialBackground.run(differenceExposure.clone()).background
detectionExposure = scoreExposure.clone()
detectionExposure.image -= background.getImage()
doSmooth = False
table = afwTable.SourceTable.make(self.schema)
self.detection.run(
table=table,
exposure=detectionExposure,
doSmooth=doSmooth,
background=background,
clearMask=True,
)
# Use the temporary detection mask for the final background subtraction,
# then restore the original mask. The score image variant never
# clears the detection planes of the difference image.
mask = differenceExposure.mask
originalMask = mask.array.copy()
detectedBit = mask.getPlaneBitMask(["DETECTED"])
detectedPix = detectionExposure.mask.array & detectedBit > 0
detectedNegativeBit = mask.getPlaneBitMask(["DETECTED_NEGATIVE"])
detectedNegativePix = detectionExposure.mask.array & detectedNegativeBit > 0
mask.array[detectedPix] |= detectedBit
mask.array[detectedNegativePix] |= detectedNegativeBit
background = self.subtractFinalBackground.run(differenceExposure).background
mask.array[...] = originalMask
if scoreExposure is not None:
# The preconvolution kernel is normalized to 1, so the same
# background level applies to the difference and score images.
scoreExposure.image -= background.getImage()
return background
def processResults(self, science, matchedTemplate, difference, sources, idFactory,
kernelSources=None, positives=None, negatives=None, measurementResults=None,
templateInfo=None):
# This is run on a background-subtracted difference image, so the
# remaining background should be ~0.
medianBg = 0.
try:
crs = findCosmicRays(
difference.maskedImage,
difference.psf,
medianBg,
pexConfig.makePropertySet(self.config.cosmicray),
self.config.cosmicray.keepCRs,
)
except LengthError:
raise TooManyCosmicRays(self.config.cosmicray.nCrPixelMax) from None
num = 0
if crs is not None:
mask = difference.maskedImage.mask
crBit = mask.getPlaneBitMask("CR")
afwDetection.setMaskFromFootprintList(mask, crs, crBit)
num = len(crs)
text = "masked" if self.config.cosmicray.keepCRs else "interpolated over"
self.log.info("Identified and %s %s cosmic rays.", text, num)
self.metadata["cosmic_ray_count"] = num
def calculateMetrics(self, science, difference, diaSources, kernelSources):
maskedImage = difference.maskedImage
# Bin the diffim to enhance low surface brightness streaks
binnedMaskedImage = afwMath.binImage(maskedImage,
self.config.streakBinFactor,
self.config.streakBinFactor)
binnedExposure = afwImage.ExposureF(binnedMaskedImage.getBBox())
binnedExposure.setMaskedImage(binnedMaskedImage)
# Clear the DETECTED mask plane before streak detection
binnedExposure.mask &= ~binnedExposure.mask.getPlaneBitMask('DETECTED')
# Rerun detection to set the DETECTED mask plane on binnedExposure
sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius()
_table = afwTable.SourceTable.make(afwTable.SourceTable.makeMinimalSchema())
self.streakDetection.run(table=_table, exposure=binnedExposure, doSmooth=True,
sigma=sigma/self.config.streakBinFactor)
binnedDetectedMaskPlane = binnedExposure.mask.array & binnedExposure.mask.getPlaneBitMask('DETECTED')
rescaledDetectedMaskPlane = binnedDetectedMaskPlane.repeat(self.config.streakBinFactor,
axis=0).repeat(self.config.streakBinFactor,
axis=1)
# Create new version of a diffim with DETECTED based on binnedExposure
streakMaskedImage = maskedImage.clone()
ysize, xsize = rescaledDetectedMaskPlane.shape
streakMaskedImage.mask.array[:ysize, :xsize] |= rescaledDetectedMaskPlane
# Detect streaks on this new version of the diffim
streaks = self.maskStreaks.run(streakMaskedImage)
streakMaskPlane = streakMaskedImage.mask.array & streakMaskedImage.mask.getPlaneBitMask('STREAK')
# Apply the new STREAK mask to the original diffim
maskedImage.mask.array |= streakMaskPlane
if self.config.writeStreakInfo:
rhos = np.array([line.rho for line in streaks.lines])
thetas = np.array([line.theta for line in streaks.lines])
sigmas = np.array([line.sigma for line in streaks.lines])
chi2s = np.array([line.reducedChi2 for line in streaks.lines])
modelMaximums = np.array([line.modelMaximum for line in streaks.lines])
streakInfo = {'rho': rhos, 'theta': thetas, 'sigma': sigmas, 'reducedChi2': chi2s,
'modelMaximum': modelMaximums}
else:
streakInfo = {'rho': np.array([]), 'theta': np.array([]), 'sigma': np.array([]),
'reducedChi2': np.array([]), 'modelMaximum': np.array([])}
return pipeBase.Struct(maskedStreaks=streakInfo)
class DetectAndMeasureScoreConnections(DetectAndMeasureConnections):
scoreExposure = pipeBase.connectionTypes.Input(
doc="Maximum likelihood image for detection.",
dimensions=("instrument", "visit", "detector"),
storageClass="ExposureF",
name="{fakesType}{coaddName}Diff_scoreTempExp",
)
scoreMeasuredExposure = pipeBase.connectionTypes.Output(
doc="Score image after backgrond subtraction and detection.",
dimensions=("instrument", "visit", "detector"),
storageClass="ExposureF",
name="{fakesType}{coaddName}Diff_scoreExp",
)
def __init__(self, *, config):
super().__init__(config=config)
if self.config.image_type == "future":
self.scoreExposure = dataclasses.replace(self.scoreExposure, storageClass="DifferenceImage")
self.scoreMeasuredExposure = dataclasses.replace(self.scoreMeasuredExposure,
storageClass="DifferenceImage")
class DetectAndMeasureScoreConfig(DetectAndMeasureConfig,
pipelineConnections=DetectAndMeasureScoreConnections):
pass
class DetectAndMeasureScoreTask(DetectAndMeasureTask):
ConfigClass = DetectAndMeasureScoreConfig
_DefaultName = "detectAndMeasureScore"
@timeMethod
def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources,
idFactory=None, measurementResults=None, templateInfo=None):
Definition at line 758 of file detectAndMeasure.py.
◆ scoreExposure
| lsst.ip.diffim.detectAndMeasure.scoreExposure : `lsst.afw.image.ExposureF`, optional |
Definition at line 891 of file detectAndMeasure.py.
◆ seed
| lsst.ip.diffim.detectAndMeasure.seed : `int` |
Definition at line 1383 of file detectAndMeasure.py.
◆ sources
| lsst.ip.diffim.detectAndMeasure.sources : `lsst.afw.table.SourceCatalog` |
Definition at line 952 of file detectAndMeasure.py.
◆ streakInfo
| lsst.ip.diffim.detectAndMeasure.streakInfo : `lsst.pipe.base.Struct` |
Definition at line 1651 of file detectAndMeasure.py.
◆ template
| lsst.ip.diffim.detectAndMeasure.template : `bool` |
Definition at line 1436 of file detectAndMeasure.py.
◆ templateEpochs
| lsst.ip.diffim.detectAndMeasure.templateEpochs : `dict` [`str`, `float`] |
Definition at line 1356 of file detectAndMeasure.py.
◆ templateInfo
| lsst.ip.diffim.detectAndMeasure.templateInfo |
Initial value:
: `list` [`lsst.images.DifferenceImageTemplateInfo`], \
optional
Definition at line 775 of file detectAndMeasure.py.
◆ wcs
| lsst.ip.diffim.detectAndMeasure.wcs : `lsst.afw.geom.SkyWcs` |
Definition at line 1434 of file detectAndMeasure.py.
Generated on for lsst.ip.diffim by
1.18.0