diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index ed6b9eca..734070c1 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -53,6 +53,16 @@ __all__ = ["DetectAndMeasureConfig", "DetectAndMeasureTask", "DetectAndMeasureScoreConfig", "DetectAndMeasureScoreTask"] +_TEMPLATE_EPOCH_NAMES = { + "TEMPLATE MJD MIN": "templateMjdMin", + "TEMPLATE MJD MAX": "templateMjdMax", + "TEMPLATE MJD MEAN": "templateMjdMean", +} +"""Image metadata key of each summary of the template epochs of the +diaSources, mapped to the name of the diaSource field and task metadata +entry it summarizes. +""" + class BadSubtractionError(pipeBase.AlgorithmError): """Raised when the residuals in footprints of stars used to compute the @@ -596,6 +606,15 @@ def __init__(self, **kwargs): 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. @@ -619,6 +638,7 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, # 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"), @@ -629,7 +649,9 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, # 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) - get_difference_image_provenance(image, "templates", connection, dataId) + imageTemplates = get_difference_image_provenance(image, "templates", connection, dataId) + if connection == "difference": + templateInfo = imageTemplates inputImages[output] = image inputs[connection] = image.to_legacy() @@ -641,9 +663,11 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext, diaSources=None, maskedStreaks=None, differenceBackground=None, + templateEpochs=None, ) try: - self.run(**inputs, idFactory=idFactory, measurementResults=measurementResults) + self.run(**inputs, idFactory=idFactory, measurementResults=measurementResults, + templateInfo=templateInfo) except pipeBase.AlgorithmError as e: # Annotate before converting error = pipeBase.AnnotatedPartialOutputsError.annotate( @@ -689,6 +713,7 @@ def convert_outputs_to_future(self, results, inputImages): for any pixels in the image. """ 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: @@ -699,6 +724,7 @@ def convert_outputs_to_future(self, results, inputImages): 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. @@ -716,7 +742,7 @@ def convert_outputs_to_future(self, results, inputImages): @timeMethod def run(self, science, matchedTemplate, difference, kernelSources=None, - idFactory=None, measurementResults=None): + idFactory=None, measurementResults=None, templateInfo=None): """Detect and measure sources on a difference image. The difference image will be convolved with a gaussian approximation of @@ -746,6 +772,11 @@ def run(self, science, matchedTemplate, difference, kernelSources=None, Result struct that is modified to allow saving of partial outputs for some failure conditions. If the task completes successfully, this is also returned. + templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \ + optional + Information about the coadds and their input images that the + template was built from. Without it the template epochs of every + diaSource are NaN. Returns ------- @@ -757,6 +788,9 @@ def run(self, science, matchedTemplate, difference, kernelSources=None, The catalog of detected sources. ``differenceBackground`` : `lsst.afw.math.BackgroundList` Background that was subtracted from the difference image. + ``templateEpochs`` : `dict` [`str`, `float`] + Summary of the template epochs of the diaSources, keyed by + image metadata key; empty if they are not known. """ if measurementResults is None: measurementResults = pipeBase.Struct() @@ -802,7 +836,8 @@ def run(self, science, matchedTemplate, difference, kernelSources=None, sources, idFactory, kernelSources, positives=positives, negatives=negatives, - measurementResults=measurementResults) + measurementResults=measurementResults, + templateInfo=templateInfo) return measurementResults def _prepareInputs(self, difference): @@ -901,7 +936,8 @@ def _fitAndSubtractBackground(self, differenceExposure, scoreExposure=None): return background def processResults(self, science, matchedTemplate, difference, sources, idFactory, - kernelSources=None, positives=None, negatives=None, measurementResults=None): + kernelSources=None, positives=None, negatives=None, measurementResults=None, + templateInfo=None): """Measure and process the results of source detection. Parameters @@ -928,6 +964,10 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor Result struct that is modified to allow saving of partial outputs for some failure conditions. If the task completes successfully, this is also returned. + templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \ + optional + Information about the coadds and their input images that the + template was built from. Returns ------- @@ -937,6 +977,9 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor Subtracted exposure with detection mask applied. ``diaSources`` : `lsst.afw.table.SourceCatalog` The catalog of detected sources. + ``templateEpochs`` : `dict` [`str`, `float`] + Summary of the template epochs of the diaSources, keyed by + image metadata key; empty if they are not known. """ if measurementResults is None: measurementResults = pipeBase.Struct() @@ -1006,6 +1049,9 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor 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(), @@ -1240,6 +1286,89 @@ def _find_glint_trails(self, diaSources): return diaSources, trail_parameters + def _measureTemplateEpochs(self, diaSources, templateInfo): + """Determine the observed epochs of the template for each diaSource. + + Parameters + ---------- + diaSources : `lsst.afw.table.SourceCatalog` + Contiguous catalog to fill in place. + templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`] or \ + `None` + Information about the coadds and their input images that the + template was built from. + + Notes + ----- + An input covers a diaSource if both its coadd's bounds and its own + bounds contain the centroid. Each coadd's inputs give a mean weighted + by their coadd weights, and where coadds overlap their means are + averaged with weights ``sqrt(W)``, ``W`` being the total weight of a + coadd's inputs there. `GetTemplateTask` merges the patches of a tract + with weights ``variance**-0.5``, and the variance of a coadd is about + ``1/W``. + """ + 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): + """Summarize the template epochs of a catalog in the task metadata. + + The earliest, latest and mean values over the rows with a template + epoch, sky sources included, are set in the task metadata under the + names of the diaSource fields they summarize. They are NaN if no row + has a template epoch. + + Parameters + ---------- + diaSources : `lsst.afw.table.SourceCatalog` + Contiguous catalog with the template epoch fields set. + + Returns + ------- + templateEpochs : `dict` [`str`, `float`] + The same values keyed by their image metadata key in + ``_TEMPLATE_EPOCH_NAMES``, or empty if no row has a template epoch. + """ + 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): """Add sources in empty regions of the difference image @@ -1617,7 +1746,7 @@ class DetectAndMeasureScoreTask(DetectAndMeasureTask): @timeMethod def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources, - idFactory=None, measurementResults=None): + idFactory=None, measurementResults=None, templateInfo=None): """Detect and measure sources on a score image. Parameters @@ -1641,6 +1770,11 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources Result struct that is modified to allow saving of partial outputs for some failure conditions. If the task completes successfully, this is also returned. + templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \ + optional + Information about the coadds and their input images that the + template was built from. Without it the template epochs of every + diaSource are NaN. Returns ------- @@ -1654,6 +1788,9 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources The catalog of detected sources. ``differenceBackground`` : `lsst.afw.math.BackgroundList` Background that was subtracted from the difference image. + ``templateEpochs`` : `dict` [`str`, `float`] + Summary of the template epochs of the diaSources, keyed by + image metadata key; empty if they are not known. """ if measurementResults is None: measurementResults = pipeBase.Struct() @@ -1713,7 +1850,8 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources sources, idFactory, kernelSources, positives=positives, negatives=negatives, - measurementResults=measurementResults) + measurementResults=measurementResults, + templateInfo=templateInfo) else: positives = afwTable.SourceCatalog(self.schema) @@ -1724,5 +1862,6 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources detectResults.sources, idFactory, kernelSources, positives=positives, negatives=negatives, - measurementResults=measurementResults) + measurementResults=measurementResults, + templateInfo=templateInfo) return measurementResults diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 41a47b00..71ab85a7 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -34,13 +34,20 @@ import lsst.afw.table as afwTable from lsst.afw.math._warper import computeWarpedBBox import lsst.afw.math as afwMath +from lsst.daf.base import DateTime import lsst.pex.config as pexConfig import lsst.pex.exceptions import lsst.pipe.base as pipeBase from lsst.images import ( + Box, DifferenceImage, DifferenceImageTemplateInfo, + DifferenceImageTemplateInput, + GeneralFrame, + NoOverlapError, + Polygon, + Transform, get_legacy_optional_mask_planes, get_legacy_template_mask_planes, ) @@ -65,6 +72,106 @@ unknown instrument. """ +_FOOTPRINT_MARGIN = 32.0 +"""Distance, in coadd pixels, by which the science image's footprint in a +coadd is grown before a coadd input is clipped to it. It must exceed twice +the distance by which a distorted science image edge departs from a straight +line between its corners, and limits how far outside the science image its +WCS is evaluated. +""" + + +def _pixelTransform(fromWcs, toWcs): + """Return a chainable transformation between two pixel frames. + + Parameters + ---------- + fromWcs, toWcs : `lsst.afw.geom.SkyWcs` + WCSs of the input and output pixels. + + Returns + ------- + transform : `lsst.images.Transform` + Transform from ``fromWcs`` pixels to ``toWcs`` pixels. + """ + return Transform.from_legacy(afwGeom.makeWcsPairTransform(fromWcs, toWcs), + GeneralFrame(unit=astropy.units.pix), + GeneralFrame(unit=astropy.units.pix) + ) + + +def _getImageOverlap(record, coaddWcs, scienceFootprint, coaddToScienceTransform, scienceBox, border): + """Return the region enclosing the overlap of the science image and one + input image of the coadd. + + Note that there are three different WCS's involved here: the WCS of the + original image that was warped to contribute to the coadd, the coadd WCS, + and the current science image WCS that we are building a template for. + + Parameters + ---------- + record : `lsst.afw.table.ExposureRecord` + Input record of a coadd, with the input's WCS, bounding box and + optional valid polygon. + coaddWcs : `lsst.afw.geom.SkyWcs` + WCS of the coadd. + scienceFootprint : `lsst.images.Region` + ``scienceBox`` in coadd pixels, grown by ``_FOOTPRINT_MARGIN``. + coaddToScienceTransform : `lsst.images.Transform` + Transform from coadd pixels to science image pixels. + scienceBox : `lsst.images.Box` + Bounding box of the science image. + border : `float` + Width, in input pixels, of the border of the input that warping left + out of the coadd. + + Returns + ------- + footprint : `lsst.images.Polygon` or `None` + The part of ``scienceBox`` the input's warp covers, in science image + pixels, or `None` if it covers none of it. + """ + inputBox = geom.Box2D(record.getBBox()) + inputBox.grow(-border) + inputBounds = Polygon.from_legacy(afwGeom.Polygon(inputBox)) + polygon = record.getValidPolygon() + inputToCoadd = _pixelTransform(record.getWcs(), coaddWcs) + try: + if polygon is not None: + inputBounds = inputBounds.intersection(Polygon.from_legacy(polygon)) + # Only vertices are transformed, so the edges of the input stay + # straight. + footprint = ( + inputBounds.transform(inputToCoadd) + .intersection(scienceFootprint) + .transform(coaddToScienceTransform) + ) + # The convex hull keeps the footprint in one piece when it is clipped. + footprint = footprint.convex_hull().intersection(scienceBox) + except NoOverlapError: + return None + return Polygon.from_box(footprint) if isinstance(footprint, Box) else footprint + + +def _warpBorder(kernelName): + """Return the width of the border of an image that warping leaves out. + + A warped pixel is computed only where the whole kernel lies on the image. + + Parameters + ---------- + kernelName : `str` + Name of the warping kernel, as `lsst.afw.math.Warper` takes it. + + Returns + ------- + border : `float` + Width of the border on each side, in pixels from the edge of the + image's bounding box; 2.5 for ``lanczos3``. + """ + kernel = afwMath.WarpingControl(kernelName).getWarpingKernel() + return kernel.getCtr().getX() + 0.5 + @functools.cache def _warpKernelSumSq(kernelName, nPhase=100): @@ -213,8 +320,10 @@ class GetTemplateConfig( "lanczos5": "Lanczos kernel of order 5", }, doc="Warping kernel that was used to build the input coadds. This" - " cannot be recovered from the coadds themselves, and is used only if" - " doCorrectVarianceWarpCorrelation is True.", + " cannot be recovered from the coadds themselves. It is used to" + " correct the variance if doCorrectVarianceWarpCorrelation is True," + " and to trim the coadd inputs bboxes to account for the size of the" + " warping kernel.", ) coaddPsf = pexConfig.ConfigField( doc="Configuration for CoaddPsf", @@ -369,16 +478,19 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): outputs.template = self.convert_outputs_to_future(outputs.template, inputRefs.coaddExposures, detector=detector, - exposureRecord=exposureRecord) + exposureRecord=exposureRecord, + templateInputs=outputs.templateInputs) butlerQC.put(outputs, outputRefs) - def convert_outputs_to_future(self, template, coadd_refs, detector, exposureRecord): + def convert_outputs_to_future(self, template, coadd_refs, detector, exposureRecord, + templateInputs=None): """Convert a template image to DifferenceImage format. This replaces ``template`` with an `lsst.images.DifferenceImage` instance and records the coadds contributed to it, along with the - instrument they were made from. The Gaussian PSF set by `run` is - converted to an `lsst.images.psfs.GaussianPointSpreadFunction`. + instrument they were made from and their input image identifiers. + The Gaussian PSF set by `run` is converted to an + `lsst.images.psfs.GaussianPointSpreadFunction`. Parameters ---------- @@ -395,6 +507,10 @@ def convert_outputs_to_future(self, template, coadd_refs, detector, exposureReco The ``exposure`` record of the observation the template was built for, which supplies the observation metadata recorded on the output. + templateInputs : `dict` [`tuple` [`int`, `int`], \ + `list` [`lsst.images.DifferenceImageTemplateInput`]], optional + Input image records of each coadd, keyed by tract and patch, as + returned by `run`. Coadds not listed record no inputs. Returns ------- @@ -440,6 +556,9 @@ def convert_outputs_to_future(self, template, coadd_refs, detector, exposureReco log=self.log, instrument=instrument, ) + if templateInputs: + for info in convertedTemplate.templates: + info.inputs = templateInputs.get((info.tract, info.patch), []) return convertedTemplate def _clearUnmappedMaskPlanes(self, mask, planeMap): @@ -603,6 +722,12 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi A template coadd exposure assembled out of patches, with pixels in nJy. Its PSF is a Gaussian approximation to the CoaddPsf of its inputs (`lsst.afw.image.ExposureF`). + ``templateInputs`` + Input image records of each coadd that cover part of the + science image, keyed by tract and patch (`dict` + [`tuple` [`int`, `int`], + `list` [`lsst.images.DifferenceImageTemplateInput`]]). Empty + unless ``config.output_image_type`` is ``"future"``. Raises ------ @@ -615,11 +740,13 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi """ band, photoCalib, templateInstrument = self._checkInputs(dataIds, coaddExposureHandles, instrument) + scienceBBox = bbox bbox = geom.Box2I(bbox) bbox.grow(self.config.templateBorderSize) warped = {} catalogs = [] + inputs = [] plateScaleFactor = 1.0 if instrument == templateInstrument else None for tract in coaddExposureHandles: @@ -634,7 +761,7 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi # coadd was from a different instrument. A single correction factor # is calculated from the first usable coadd and applied to all. if plateScaleFactor is None: - plateScaleFactor = self._plateScaleFactor(coaddInputs, wcs, bbox) + plateScaleFactor = self._plateScaleFactor(coaddInputs[0], wcs, bbox) warpedBox = computeWarpedBBox(catalog[0].wcs, bbox, wcs) warpedBox.grow(5) # to ensure we catch all relevant input pixels # Combine images from individual patches together. @@ -670,6 +797,7 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi tempCatalog.reserve(len(included)) for i in included: tempCatalog.append(catalog[i]) + inputs.append(coaddInputs[i]) catalogs.append(tempCatalog) warped[tract] = potentialInput.maskedImage @@ -696,6 +824,11 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi catalog.reserve(sum([len(c) for c in catalogs])) for c in catalogs: catalog.extend(c) + # Only the future output type records the inputs. + if self.config.output_image_type == "future": + templateInputs = self._collectTemplateInputs(catalog, inputs, wcs, scienceBBox) + else: + templateInputs = {} # Set a mask plane for any regions with exceptionally high variance. self.checkHighVariance(template) @@ -714,7 +847,7 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi coaddInputs = afwImage.CoaddInputs(afwTable.ExposureTable.makeMinimalSchema(), self.schema) coaddInputs.ccds.extend(catalog, deep=True) template.getInfo().setCoaddInputs(coaddInputs) - return pipeBase.Struct(template=template) + return pipeBase.Struct(template=template, templateInputs=templateInputs) def _correctVariance(self, template, plateScaleFactor): """Correct the template variance plane for the effects of warping. @@ -817,6 +950,59 @@ def _plateScaleFactor(self, coaddInputs, wcs, bbox): return factor return None + def _collectTemplateInputs(self, catalog, coaddInputs, wcs, bbox): + """Find the input images and epochs of each coadd in the template. + + Parameters + ---------- + catalog : `lsst.afw.table.ExposureCatalog` + Coadds the template was built from, with their tract, patch and + WCS. + coaddInputs : `list` [`lsst.afw.image.CoaddInputs` or `None`] + Input records of each coadd in ``catalog``, in the same order. + wcs : `lsst.afw.geom.SkyWcs` + WCS of the science image. + bbox : `lsst.geom.Box2I` + Bounding box of the science image. + + Returns + ------- + templateInputs : `dict` [`tuple` [`int`, `int`], \ + `list` [`lsst.images.DifferenceImageTemplateInput`]] + Inputs of each coadd that cover part of the science image, keyed + by tract and patch, with the part their warps cover in science + image pixels. + """ + scienceBox = Box.from_legacy(bbox) + border = _warpBorder(self.config.coaddWarpKernel) + regions = {} + templateInputs = {} + for coadd, coaddInput in zip(catalog, coaddInputs, strict=True): + coaddList = templateInputs.setdefault((coadd["tract"], coadd["patch"]), []) + if coaddInput is None: + continue + scienceToCoadd = _pixelTransform(wcs, coadd.getWcs()) + scienceFootprint = scienceBox.transform(scienceToCoadd).dilated_by(_FOOTPRINT_MARGIN) + coaddToScience = scienceToCoadd.inverted() + for record in coaddInput.ccds: + visitInfo = record.getVisitInfo() + weight = record["weight"] + if visitInfo is None or not visitInfo.getDate().isValid() or not weight > 0: + continue + key = (record["visit"], record["ccd"]) + # Record a single region for a ccd in overlapping patches. + if key not in regions: + regions[key] = _getImageOverlap(record, coadd.getWcs(), scienceFootprint, coaddToScience, + scienceBox, border) + region = regions[key] + if region is None: + continue + mjd = visitInfo.getDate().get(DateTime.MJD, DateTime.TAI) + coaddList.append(DifferenceImageTemplateInput(visit=key[0], detector=key[1], + mid_mjd_tai=mjd, weight=weight, + bounds=region)) + return templateInputs + def checkHighVariance(self, template): """Set a mask plane for regions with unusually high variance. @@ -952,16 +1138,15 @@ def _makeExposureCatalog(self, exposureRefs, dataIds, instrument=_UNKNOWN_INSTRU Catalog of metadata for each exposure totalBox : `lsst.geom.Box2I` The union of the bounding boxes of all the input exposures. - coaddInputs : `lsst.afw.image.CoaddInputs` or `None` - Input records of the first exposure kept, or `None` if none is - kept. + coaddInputs : `list` [`lsst.afw.image.CoaddInputs` or `None`] + Input records of each exposure in ``catalog``, in the same order. """ catalog = afwTable.ExposureCatalog(self.schema) catalog.reserve(len(exposureRefs)) exposures = (exposureRef.get() for exposureRef in exposureRefs) images = {} totalBox = geom.Box2I() - firstCoaddInputs = None + keptCoaddInputs = [] for coadd, dataId in zip(exposures, dataIds): coaddInputs = coadd.getInfo().getCoaddInputs() @@ -974,8 +1159,7 @@ def _makeExposureCatalog(self, exposureRefs, dataIds, instrument=_UNKNOWN_INSTRU instrument, ) continue - if not images: - firstCoaddInputs = coaddInputs + keptCoaddInputs.append(coaddInputs) images[dataId] = coadd.maskedImage bbox = coadd.getBBox() totalBox = totalBox.expandedTo(bbox) @@ -992,7 +1176,7 @@ def _makeExposureCatalog(self, exposureRefs, dataIds, instrument=_UNKNOWN_INSTRU # should be very similar, so this value mostly shouldn't matter. record.set("weight", 1) - return images, catalog, totalBox, firstCoaddInputs + return images, catalog, totalBox, keptCoaddInputs def _merge(self, maskedImages, bbox, wcs): """Merge the images that came from one tract into one larger image, diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 84b5d5aa..910f2ac2 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -34,8 +34,8 @@ import lsst.afw.table as afwTable import lsst.afw.math as afwMath import lsst.geom -from lsst.images import (Box, DifferenceImage, DifferenceImageTemplateInfo, Polygon, - get_legacy_difference_image_mask_planes, get_legacy_optional_mask_planes) +from lsst.images import (Box, DifferenceImage, DifferenceImageTemplateInfo, DifferenceImageTemplateInput, + Polygon, get_legacy_difference_image_mask_planes, get_legacy_optional_mask_planes) from lsst.images.convolution_kernels import ImageBasisConvolutionKernel from lsst.images.fields import ChebyshevField import lsst.images.psfs @@ -2032,6 +2032,205 @@ def test_run_quantum_score_task_provenance(self): self.assertIs(out.photometric_scaling, photometricScaling) self.assertIs(scoreOut.photometric_scaling, photometricScaling) + @staticmethod + def _makeEpochTemplates(bbox): + """Return two coadd records covering the left and right parts of an + image, overlapping in a strip down the middle, with three inputs + between them. + + Neither coadd covers the lowest quarter of the image. Input 1 is + listed by both coadds, and input 2 only by the left one; both cover + the whole image. Input 3 is listed only by the right coadd and covers + the upper half. + + Parameters + ---------- + bbox : `lsst.geom.Box2I` + Bounding box of the image. + + Returns + ------- + templates : `list` [`lsst.images.DifferenceImageTemplateInfo`] + The left and right coadd records. + xLeft, xRight : `float` + Right edge of the left coadd and left edge of the right one. + yMid : `float` + Lower edge of input 3. + yLow : `float` + Lower edge of both coadds. + """ + x0, x1 = bbox.getMinX() - 0.5, bbox.getMaxX() + 0.5 + y0, y1 = bbox.getMinY() - 0.5, bbox.getMaxY() + 0.5 + xLeft = x0 + 0.6*(x1 - x0) + xRight = x0 + 0.4*(x1 - x0) + yMid = 0.5*(y0 + y1) + yLow = y0 + 0.25*(y1 - y0) + + def box(xa, xb, ya, yb): + return Polygon(x_vertices=[xa, xb, xb, xa], y_vertices=[ya, ya, yb, yb]) + + def makeInput(visit, mjd, weight, bounds): + return DifferenceImageTemplateInput(visit=visit, detector=5, mid_mjd_tai=mjd, weight=weight, + bounds=bounds) + + whole = box(x0, x1, y0, y1) + templates = [] + for patch, coaddBounds, inputs in ( + (1, box(x0, xLeft, yLow, y1), [makeInput(1, 60000.0, 1.0, whole), + makeInput(2, 60010.0, 3.0, whole)]), + (2, box(xRight, x1, yLow, y1), [makeInput(1, 60000.0, 1.0, whole), + makeInput(3, 60020.0, 1.0, box(x0, x1, yMid, y1))]), + ): + info = DetectAndMeasureOutputFormatTest._makeTemplateInfo(patch=patch) + info.bounds = coaddBounds + info.inputs = inputs + templates.append(info) + return templates, xLeft, xRight, yMid, yLow + + def _checkTemplateEpochs(self, diaSources, xLeft, xRight, yMid, yLow): + """Check the template epochs of diaSources measured with the records + of `_makeEpochTemplates`. + """ + x = diaSources.getX() + y = diaSources.getY() + covered = y >= yLow + left = covered & (x < xLeft) + right = covered & (x >= xRight) + upper = y >= yMid + for region in (left & ~right, right & ~left & upper, right & ~left & ~upper, left & right & upper, + ~covered): + self.assertTrue(region.any(), msg="The test image needs a diaSource in every region.") + # The left coadd has total weight 4 everywhere, the right one 2 where + # input 3 covers it and 1 elsewhere; overlapping coadds are averaged + # with weights sqrt(total weight). + leftMean = (60000.0 + 3.0*60010.0)/4.0 + rightWeight = np.where(upper, 2.0, 1.0) + rightMean = np.where(upper, (60000.0 + 60020.0)/2.0, 60000.0) + expected = ((left*2.0*leftMean + right*np.sqrt(rightWeight)*rightMean) + / np.where(covered, left*2.0 + right*np.sqrt(rightWeight), 1.0)) + expectedMax = np.where(right & upper, 60020.0, np.where(left, 60010.0, 60000.0)) + self.assertFloatsAlmostEqual(diaSources["templateMjdMean"][covered], expected[covered], rtol=1e-12) + self.assertTrue(np.all(diaSources["templateMjdMin"][covered] == 60000.0)) + self.assertFloatsEqual(diaSources["templateMjdMax"][covered], expectedMax[covered]) + for name in ("templateMjdMean", "templateMjdMin", "templateMjdMax"): + self.assertTrue(np.all(np.isnan(diaSources[name][~covered])), msg=name) + + _EPOCH_NAMES = {"TEMPLATE MJD MIN": "templateMjdMin", + "TEMPLATE MJD MAX": "templateMjdMax", + "TEMPLATE MJD MEAN": "templateMjdMean"} + """Image metadata key of each summary of the diaSource template epochs, + mapped to the field and task metadata name it summarizes. + """ + + def _checkTemplateEpochSummary(self, task, templateEpochs, diaSources): + """Check the summary of the template epochs of a catalog in the task + metadata and in the result struct. + """ + known = np.isfinite(diaSources["templateMjdMean"]) + self.assertTrue(known.any()) + expected = {"templateMjdMin": np.min(diaSources["templateMjdMin"][known]), + "templateMjdMax": np.max(diaSources["templateMjdMax"][known]), + "templateMjdMean": np.mean(diaSources["templateMjdMean"][known])} + self.assertEqual(set(templateEpochs), set(self._EPOCH_NAMES)) + for key, name in self._EPOCH_NAMES.items(): + with self.subTest(name=name): + self.assertFloatsAlmostEqual(task.metadata[name], expected[name], rtol=1e-12) + self.assertEqual(templateEpochs[key], task.metadata[name]) + + def _checkNoTemplateEpochSummary(self, task, templateEpochs): + """Check that no template epoch summary is recorded beyond NaN task + metadata. + """ + self.assertEqual(templateEpochs, {}) + for name in self._EPOCH_NAMES.values(): + self.assertTrue(np.isnan(task.metadata[name]), msg=name) + + def test_template_epochs(self): + """Each diaSource records the observation times of the template + inputs that cover its centroid, counting an input only inside the + bounds of a coadd that lists it, and the task records their summary. + """ + science, sources, matchedTemplate, difference = self._make_images() + templates, *geometry = self._makeEpochTemplates(difference.getBBox()) + task = self._setup_task() + output = task.run(science, matchedTemplate, difference, sources, templateInfo=templates) + self._checkTemplateEpochs(output.diaSources, *geometry) + self._checkTemplateEpochSummary(task, output.templateEpochs, output.diaSources) + + # Without template records the epochs are not known. + science, sources, matchedTemplate, difference = self._make_images() + task = self._setup_task() + output = task.run(science, matchedTemplate, difference, sources) + for name in ("templateMjdMean", "templateMjdMin", "templateMjdMax"): + self.assertTrue(np.all(np.isnan(output.diaSources[name])), msg=name) + self._checkNoTemplateEpochSummary(task, output.templateEpochs) + + def test_template_epoch_summary_includes_sky_sources(self): + """Sky sources are measured and summarized like diaSources. + """ + science, sources, matchedTemplate, difference = self._make_images() + templates, *_ = self._makeEpochTemplates(difference.getBBox()) + task = self._setup_detection(doSkySources=True, nSkySources=20, badSubtractionRatioThreshold=1.) + output = task.run(science, matchedTemplate, difference, sources, templateInfo=templates) + diaSources = output.diaSources + sky = diaSources["sky_source"] + self.assertTrue(np.isfinite(diaSources["templateMjdMean"][sky]).any()) + self._checkTemplateEpochSummary(task, output.templateEpochs, diaSources) + + def test_template_epochs_score(self): + """``DetectAndMeasureScoreTask`` measures the template epochs too, and + both converted images record their summary. + """ + self.detectionTask = detectAndMeasure.DetectAndMeasureScoreTask + science, sources, matchedTemplate, difference, score = self._make_images(withScore=True) + templates, *geometry = self._makeEpochTemplates(difference.getBBox()) + task = self._setup_task(image_type="future") + inputImages = {"subtractedMeasuredExposure": self._makeFutureInput(difference), + "scoreMeasuredExposure": self._makeFutureInput(score)} + output = task.run(science, matchedTemplate, + inputImages["subtractedMeasuredExposure"].to_legacy(), + inputImages["scoreMeasuredExposure"].to_legacy(), sources, templateInfo=templates) + self._checkTemplateEpochs(output.diaSources, *geometry) + task.convert_outputs_to_future(output, inputImages) + for name in inputImages: + for key, value in output.templateEpochs.items(): + with self.subTest(image=name, key=key): + self.assertEqual(getattr(output, name).metadata.native[key], value) + + def test_run_quantum_template_epochs(self): + """In future mode ``runQuantum`` measures the template epochs from the + template records of the difference image and records their summary + on the output image; in legacy mode they are NaN. + """ + for imageType, withInputs in (("legacy", True), ("future", True), ("future", False)): + with self.subTest(imageType=imageType, withInputs=withInputs): + science, sources, matchedTemplate, difference = self._make_images() + templates, *geometry = self._makeEpochTemplates(difference.getBBox()) + if not withInputs: + for info in templates: + info.inputs = [] + if imageType == "future": + difference = self._makeFutureInput(difference, templates, self._makeKernel(), + self._makePhotometricScaling(difference)) + task = self._setup_task(image_type=imageType) + butlerQC = _RecordingQuantumContext(self.dataId) + task.runQuantum(butlerQC, + _FakeRefs(science=science, matchedTemplate=matchedTemplate, + difference=difference, kernelSources=sources), + _FakeRefs()) + diaSources = butlerQC.put_values.diaSources + out = butlerQC.put_values.subtractedMeasuredExposure + if imageType == "future" and withInputs: + self._checkTemplateEpochs(diaSources, *geometry) + self._checkTemplateEpochSummary(task, butlerQC.put_values.templateEpochs, diaSources) + for key, name in self._EPOCH_NAMES.items(): + self.assertEqual(out.metadata.native[key], task.metadata[name]) + else: + self.assertTrue(np.all(np.isnan(diaSources["templateMjdMean"]))) + self._checkNoTemplateEpochSummary(task, butlerQC.put_values.templateEpochs) + for key in self._EPOCH_NAMES: + self.assertNotIn(key, out.metadata) + def test_connections_future_reads_difference_image(self): """The inputs whose provenance is copied ask for the type that carries it, and only in future mode. diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index 38da9c53..64c6de03 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -26,6 +26,7 @@ import astropy.units as u from astropy.stats import gaussian_sigma_to_fwhm +import astshim import numpy as np import lsst.afw.cameraGeom.testUtils @@ -34,6 +35,7 @@ import lsst.afw.image import lsst.afw.math import lsst.afw.table +from lsst.daf.base import DateTime from lsst.daf.butler import DataCoordinate, DatasetRef, DatasetType, DimensionUniverse import lsst.geom import lsst.images @@ -502,7 +504,7 @@ def _makeScaledWcs(self, factor): cdMatrix) @staticmethod - def _makeCoaddInputs(records, instrument="TestCam"): + def _makeCoaddInputs(records, instrument="TestCam", epochs=None): """Make a CoaddInputs holding the given input records. Parameters @@ -516,30 +518,39 @@ def _makeCoaddInputs(records, instrument="TestCam"): Instrument label to record for every input. `None` leaves the inputs without a visit info, and an empty string gives them a visit info with no label; either way the instrument is unknown. + epochs : `list` [`tuple` [`int`, `float`, `float`]], optional + Visit, mid-exposure MJD (TAI) and coadd weight of each input. By + default every input is visit 0 with weight 1 and no observation + time. The time is not recorded if ``instrument`` is `None`. """ ccdSchema = lsst.afw.table.ExposureTable.makeMinimalSchema() weightKey = ccdSchema.addField("weight", type=float, doc="Coadd weight") + visitKey = ccdSchema.addField("visit", type=np.int64, doc="Visit of the input") + ccdSchema.addField("ccd", type=np.int32, doc="Detector of the input") coaddInputs = lsst.afw.image.CoaddInputs( lsst.afw.table.ExposureTable.makeMinimalSchema(), ccdSchema) - for wcs, bbox in records: + if epochs is None: + epochs = [(0, None, 1.0)]*len(records) + for (wcs, bbox), (visit, mjd, weight) in zip(records, epochs, strict=True): record = coaddInputs.ccds.addNew() record.setWcs(wcs) record.setBBox(bbox) + record.set(visitKey, visit) + record.set(weightKey, weight) if instrument is not None: - record.setVisitInfo(lsst.afw.image.VisitInfo(instrumentLabel=instrument)) - # Included because real coadds have it, though a single-record - # correction does not use it. - record.set(weightKey, 1.0) + date = DateTime(mjd, DateTime.MJD, DateTime.TAI) if mjd is not None else DateTime() + record.setVisitInfo(lsst.afw.image.VisitInfo(instrumentLabel=instrument, date=date)) return coaddInputs - def _patchHandles(self, tract, records, instrument="TestCam"): + def _patchHandles(self, tract, records, instrument="TestCam", epochs=None): """Return handles for a tract's patches, with their CoaddInputs - replaced by ``records`` from ``instrument``. + replaced by ``records`` from ``instrument``, observed at ``epochs``. """ handles = [] for ref in self.patches[tract]: coadd = ref.get() - coadd.getInfo().setCoaddInputs(self._makeCoaddInputs(records, instrument=instrument)) + coadd.getInfo().setCoaddInputs(self._makeCoaddInputs(records, instrument=instrument, + epochs=epochs)) handles.append(pipeBase.InMemoryDatasetHandle( coadd, storageClass="ExposureF", copy=True, dataId=ref.dataId)) return handles @@ -697,14 +708,112 @@ def testCorrectVariancePlateScaleNeedsCoaddInputs(self): with self.assertRaisesRegex(RuntimeError, "cannot be corrected for the change in pixel area"): self._runCorrection(handles=handles) - def _runTract0(self, handles): - """Build a template from the given tract 0 patch handles. + def testRunReturnsTemplateInputs(self): + """`run` returns the inputs of each coadd that overlap the science + image, with the part of the science image each one covers. Inputs + with no weight or no observation time are left out, as are inputs + that cover only the template border. + """ + wcs = self.exposure.wcs + whole = (wcs, self.exposure.getBBox()) + right = (wcs, lsst.geom.Box2I(lsst.geom.Point2I(90, -200), lsst.geom.Point2I(800, 800))) + outside = (wcs, lsst.geom.Box2I(lsst.geom.Point2I(600, 600), lsst.geom.Point2I(800, 800))) + border = (wcs, lsst.geom.Box2I(lsst.geom.Point2I(185, 0), lsst.geom.Point2I(195, 180))) + handles = (self._patchHandles(0, [whole, right, outside, whole, whole, border], + epochs=[(1, 60000.25, 1.0), (2, 60010.5, 3.0), (3, 60100.0, 1.0), + (4, 60200.0, 0.0), (5, None, 1.0), (6, 60300.0, 1.0)])[:1] + + self._patchHandles(0, [whole], epochs=[(1, 60000.25, 2.0)])[1:]) + result = self._runTract0(handles, outputImageType="future") + + templateInputs = result.templateInputs + self.assertEqual(set(templateInputs), {(0, patch) for patch in range(4)}) + first = templateInputs[0, 0] + self.assertEqual([(i.visit, i.detector, i.weight) for i in first], [(1, 0, 1.0), (2, 0, 3.0)]) + self.assertFloatsAlmostEqual(np.array([i.mid_mjd_tai for i in first]), np.array([60000.25, 60010.5]), + rtol=0, atol=1e-6) + # The box that _runTract0 requests. + scienceBox = lsst.geom.Box2D(lsst.geom.Point2D(-0.5, -0.5), lsst.geom.Point2D(180.5, 180.5)) + self.assertFloatsAlmostEqual(first[0].bounds.area, scienceBox.getArea(), rtol=1e-6) + # Warping with lanczos5 left the 4.5 pixels along the input's edges + # out of the coadd. + rightArea = (scienceBox.getMaxX() - 94.0)*scienceBox.getHeight() + self.assertFloatsAlmostEqual(first[1].bounds.area, rightArea, rtol=1e-6) + self.assertFloatsAlmostEqual(first[1].bounds.x_vertices.min(), 94.0, rtol=0, atol=1e-6) + # Each coadd lists its own record of an input shared with another. + for patch in range(1, 4): + self.assertEqual([(i.visit, i.weight) for i in templateInputs[0, patch]], [(1, 2.0)]) + + def testRunTemplateInputsOnlyForFuture(self): + """`run` records no template inputs for the legacy output type. + """ + records = [(self.exposure.wcs, self.exposure.getBBox())] + result = self._runTract0(self._patchHandles(0, records, epochs=[(1, 60000.25, 1.0)])) + self.assertEqual(result.templateInputs, {}) + + def testRunTemplateInputsSkipDroppedCoadd(self): + """A coadd dropped for its instrument has no template inputs. + """ + records = [(self.exposure.wcs, self.exposure.getBBox())] + handles = (self._patchHandles(0, records, epochs=[(1, 60000.25, 1.0)])[:3] + + self._patchHandles(0, records, instrument="OtherCam", + epochs=[(2, 60100.0, 1.0)])[3:]) + result = self._runTract0(handles, outputImageType="future") + self.assertEqual(set(result.templateInputs), {(0, 0), (0, 1), (0, 2)}) + + def testRunTemplateInputsLeaveOutWarpBorder(self): + """An input's footprint leaves out the border along its edges that + the coadd's warping kernel left out. + """ + right = (self.exposure.wcs, lsst.geom.Box2I(lsst.geom.Point2I(90, -200), lsst.geom.Point2I(800, 800))) + for kernel, border in (("bilinear", 0.5), ("lanczos3", 2.5)): + with self.subTest(coaddWarpKernel=kernel): + result = self._runTract0(self._patchHandles(0, [right], epochs=[(1, 60000.25, 1.0)]), + outputImageType="future", coaddWarpKernel=kernel) + bounds = result.templateInputs[0, 0][0].bounds + self.assertFloatsAlmostEqual(bounds.x_vertices.min(), 89.5 + border, rtol=0, atol=1e-6) + + def testRunTemplateInputsCoverScienceEdges(self): + """An input's footprint reaches the edges of a distorted science image + wherever the input covers them. + """ + # Science pixel (x, y) lands on input pixel (x, y + a*(x - xc)**2), + # which distorts the top and bottom edges of the science image by about + # 2 pixels between its corners. + xc = 90.0 + left = (self.exposure.wcs, lsst.geom.Box2I(lsst.geom.Point2I(-200, -200), lsst.geom.Point2I(89, 800))) + # The centers of the edge pixels that the input's warp covers; warping + # with lanczos5 left the 4.5 pixels along the input's edges out. + x = np.arange(0.0, 85.0) + y = np.arange(0.0, 181.0) + edgeX = np.concatenate([x, x, np.zeros_like(y)]) + edgeY = np.concatenate([np.zeros_like(x), np.full_like(x, 180.0), y]) + for a in (-2.5e-4, 2.5e-4): + with self.subTest(a=a): + forward = np.array([[1.0, 1, 1, 0], [1.0, 2, 0, 1], [a, 2, 2, 0], [-2*a*xc, 2, 1, 0], + [a*xc**2, 2, 0, 0]]) + inverse = np.array([[1.0, 1, 1, 0], [1.0, 2, 0, 1], [-a, 2, 2, 0], [2*a*xc, 2, 1, 0], + [-a*xc**2, 2, 0, 0]]) + distortion = lsst.afw.geom.TransformPoint2ToPoint2(astshim.PolyMap(forward, inverse)) + wcs = lsst.afw.geom.makeModifiedWcs(distortion, self.exposure.wcs, False) + result = self._runTract0(self._patchHandles(0, [left], epochs=[(1, 60000.25, 1.0)]), + outputImageType="future", wcs=wcs) + bounds = result.templateInputs[0, 0][0].bounds + self.assertTrue(np.all(bounds.contains(x=edgeX, y=edgeY))) + + def _runTract0(self, handles, outputImageType="legacy", wcs=None, coaddWarpKernel=None): + """Build a template from the given tract 0 patch handles, with a task + configured for ``outputImageType`` and ``coaddWarpKernel``, for a + science image with ``wcs``, by default that of the test exposure. """ box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) - task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) + config = self._makeConfig() + config.output_image_type = outputImageType + if coaddWarpKernel is not None: + config.coaddWarpKernel = coaddWarpKernel + task = lsst.ip.diffim.GetTemplateTask(config=config) return task.run(coaddExposureHandles={0: handles}, bbox=box, - wcs=self.exposure.wcs, + wcs=self.exposure.wcs if wcs is None else wcs, dataIds={0: self.dataIds[0]}, physical_filter="a_test") @@ -982,7 +1091,7 @@ def _runFuture(self, handles=None): result, box = self._runLegacyForFuture(handles=handles) result.template = self.futureTask.convert_outputs_to_future( result.template, self._coaddRefs(0), detector=detector, - exposureRecord=self._exposureRecord()) + exposureRecord=self._exposureRecord(), templateInputs=result.templateInputs) return result, box, detector def _exposureRecord(self): @@ -1159,6 +1268,19 @@ def testConvertOutputsToFutureInstrument(self): result, _, _ = self._runFuture(handles=self._patchHandles(0, records, instrument=None)) self.assertEqual({template.instrument for template in result.template.templates}, {"UNKNOWN"}) + def testConvertOutputsToFutureTemplateInputs(self): + """Each coadd recorded on the converted template lists its inputs. + """ + records = [(self.exposure.wcs, self.exposure.getBBox())] + handles = self._patchHandles(0, records, epochs=[(1, 60000.25, 1.0)]) + result, _, _ = self._runFuture(handles=handles) + templates = result.template.templates + self.assertGreater(len(templates), 0) + for template in templates: + with self.subTest(tract=template.tract, patch=template.patch): + self.assertEqual([(i.visit, i.detector, i.weight) for i in template.inputs], [(1, 0, 1.0)]) + self.assertFloatsAlmostEqual(template.inputs[0].mid_mjd_tai, 60000.25, rtol=0, atol=1e-6) + def testConvertOutputsToFuturePsf(self): """The Gaussian PSF set by `run` is kept, defined over the detector. """