From e9c18e9ccb67b854d46f12101811c425eb8e3046 Mon Sep 17 00:00:00 2001 From: BrunoSanchez Date: Mon, 2 Jun 2025 14:41:46 -0700 Subject: [PATCH 1/2] Replace legacy fake matching tasks with new ones. Deprecation was in place --- .../pipe/tasks/matchDiffimSourceInjected.py | 433 +++++++++++++++++ python/lsst/pipe/tasks/matchFakes.py | 445 ------------------ tests/test_matchFakes.py | 141 ------ tests/test_matchSourceInjected.py | 188 ++++++++ 4 files changed, 621 insertions(+), 586 deletions(-) create mode 100644 python/lsst/pipe/tasks/matchDiffimSourceInjected.py delete mode 100644 python/lsst/pipe/tasks/matchFakes.py delete mode 100644 tests/test_matchFakes.py create mode 100644 tests/test_matchSourceInjected.py diff --git a/python/lsst/pipe/tasks/matchDiffimSourceInjected.py b/python/lsst/pipe/tasks/matchDiffimSourceInjected.py new file mode 100644 index 000000000..6be9f9b22 --- /dev/null +++ b/python/lsst/pipe/tasks/matchDiffimSourceInjected.py @@ -0,0 +1,433 @@ +# This file is part of ap_pipe. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +__all__ = ["MatchInjectedToDiaSourceTask", + "MatchInjectedToDiaSourceConfig", + "MatchInjectedToAssocDiaSourceTask", + "MatchInjectedToAssocDiaSourceConfig"] + +import astropy.units as u +import numpy as np +from scipy.spatial import cKDTree + +from lsst.afw import table as afwTable +from lsst import geom as lsstGeom +import lsst.pex.config as pexConfig +from lsst.pipe.base import PipelineTask, PipelineTaskConnections, Struct +import lsst.pipe.base.connectionTypes as connTypes +from lsst.meas.base import ForcedMeasurementTask, ForcedMeasurementConfig +from lsst.source.injection import VisitInjectConfig + + +class MatchInjectedToDiaSourceConnections( + PipelineTaskConnections, + defaultTemplates={"coaddName": "deep", + "fakesType": "fakes_"}, + dimensions=("instrument", + "visit", + "detector")): + injectedCat = connTypes.Input( + doc="Catalog of sources injected in the images.", + name="{fakesType}_pvi_catalog", + storageClass="ArrowAstropy", + dimensions=("instrument", "visit", "detector"), + ) + diffIm = connTypes.Input( + doc="Difference image on which the DiaSources were detected.", + name="{fakesType}{coaddName}Diff_differenceExp", + storageClass="ExposureF", + dimensions=("instrument", "visit", "detector"), + ) + diaSources = connTypes.Input( + doc="A DiaSource catalog to match against fakeCat.", + name="{fakesType}{coaddName}Diff_diaSrc", + storageClass="SourceCatalog", + dimensions=("instrument", "visit", "detector"), + ) + matchDiaSources = connTypes.Output( + doc="A catalog of those fakeCat sources that have a match in " + "diaSrc. The schema is the union of the schemas for " + "``fakeCat`` and ``diaSrc``.", + name="{fakesType}{coaddName}Diff_matchDiaSrc", + storageClass="DataFrame", + dimensions=("instrument", "visit", "detector"), + ) + + +class MatchInjectedToDiaSourceConfig( + VisitInjectConfig, + pipelineConnections=MatchInjectedToDiaSourceConnections): + """Config for MatchFakesTask. + """ + matchDistanceArcseconds = pexConfig.RangeField( + doc="Distance in arcseconds to match within.", + dtype=float, + default=0.5, + min=0, + max=10, + ) + doMatchVisit = pexConfig.Field( + dtype=bool, + default=True, + doc="Match visit to trim the fakeCat" + ) + trimBuffer = pexConfig.Field( + doc="Size of the pixel buffer surrounding the image." + "Only those fake sources with a centroid" + "falling within the image+buffer region will be considered matches.", + dtype=int, + default=50, + ) + doForcedMeasurement = pexConfig.Field( + dtype=bool, + default=True, + doc="Force measurement of the fakes at the injection locations." + ) + forcedMeasurement = pexConfig.ConfigurableField( + target=ForcedMeasurementTask, + doc="Task to force photometer difference image at injection locations.", + ) + + +class MatchInjectedToDiaSourceTask(PipelineTask): + + _DefaultName = "matchInjectedToDiaSource" + ConfigClass = MatchInjectedToDiaSourceConfig + + def run(self, injectedCat, diffIm, diaSources): + """Match injected sources to detected diaSources within a difference image bound. + + Parameters + ---------- + injectedCat : `astropy.table.table.Table` + Table of catalog of synthetic sources to match to detected diaSources. + diffIm : `lsst.afw.image.Exposure` + Difference image where ``diaSources`` were detected. + diaSources : `afw.table.SourceCatalog` + Catalog of difference image sources detected in ``diffIm``. + assocDiaSources : `pandas.DataFrame` + Catalog of associated difference image sources detected in ``diffIm``. + Returns + ------- + result : `lsst.pipe.base.Struct` + Results struct with components. + + - ``matchedDiaSources`` : Fakes matched to input diaSources. Has + length of ``injectedCalexpCat``. (`pandas.DataFrame`) + """ + + if self.config.doMatchVisit: + fakeCat = self._trimFakeCat(injectedCat, diffIm) + else: + fakeCat = injectedCat + if self.config.doForcedMeasurement: + self._estimateFakesSNR(fakeCat, diffIm) + + return self._processFakes(fakeCat, diaSources) + + def _estimateFakesSNR(self, injectedCat, diffIm): + """Estimate the signal-to-noise ratio of the fakes in the given catalog. + + Parameters + ---------- + injectedCat : `astropy.table.Table` + Catalog of synthetic sources to estimate the S/N of. **This table + will be modified in place**. + diffIm : `lsst.afw.image.Exposure` + Difference image where the sources were detected. + """ + # Create a schema for the forced measurement task + schema = afwTable.SourceTable.makeMinimalSchema() + schema.addField("x", "D", "x position in image.", units="pixel") + schema.addField("y", "D", "y position in image.", units="pixel") + schema.addField("deblend_nChild", "I", "Need for minimal forced phot schema") + + pluginList = [ + "base_PixelFlags", + "base_SdssCentroid", + "base_CircularApertureFlux", + "base_PsfFlux", + "base_LocalBackground" + ] + forcedMeasConfig = ForcedMeasurementConfig(plugins=pluginList) + forcedMeasConfig.slots.centroid = 'base_SdssCentroid' + forcedMeasConfig.slots.shape = None + + # Create the forced measurement task + forcedMeas = ForcedMeasurementTask(schema, config=forcedMeasConfig) + + # Specify the columns to copy from the input catalog to the output catalog + forcedMeas.copyColumns = {"coord_ra": "ra", "coord_dec": "dec"} + + # Create an afw table from the input catalog + outputCatalog = afwTable.SourceCatalog(schema) + outputCatalog.reserve(len(injectedCat)) + for row in injectedCat: + outputRecord = outputCatalog.addNew() + outputRecord.setId(row['injection_id']) + outputRecord.setCoord(lsstGeom.SpherePoint(row["ra"], row["dec"], lsstGeom.degrees)) + outputRecord.set("x", row["x"]) + outputRecord.set("y", row["y"]) + + # Generate the forced measurement catalog + forcedSources = forcedMeas.generateMeasCat(diffIm, outputCatalog, diffIm.getWcs()) + # Attach the PSF shape footprints to the forced measurement catalog + forcedMeas.attachPsfShapeFootprints(forcedSources, diffIm) + + # Copy the x and y positions from the forced measurement catalog back + # to the input catalog + for src, tgt in zip(forcedSources, outputCatalog): + src.set('base_SdssCentroid_x', tgt['x']) + src.set('base_SdssCentroid_y', tgt['y']) + + # Define the centroid for the forced measurement catalog + forcedSources.defineCentroid('base_SdssCentroid') + # Run the forced measurement task + forcedMeas.run(forcedSources, diffIm, outputCatalog, diffIm.getWcs()) + # Convert the forced measurement catalog to an astropy table + forcedSources_table = forcedSources.asAstropy() + + # Add the forced measurement columns to the input catalog + for column in forcedSources_table.columns: + if "Flux" in column or "flag" in column: + injectedCat["forced_"+column] = forcedSources_table[column] + + # Add the SNR columns to the input catalog + for column in injectedCat.colnames: + if column.endswith("instFlux"): + flux = injectedCat[column] + fluxErr = injectedCat[column+"Err"].copy() + fluxErr = np.where( + (fluxErr <= 0) | (np.isnan(fluxErr)), np.nanmax(fluxErr), fluxErr) + + injectedCat[column+"_SNR"] = flux / fluxErr + + def _processFakes(self, injectedCat, diaSources): + """Match fakes to detected diaSources within a difference image bound. + + Parameters + ---------- + injectedCat : `astropy.table.table.Table` + Catalog of injected sources to match to detected diaSources. + diaSources : `afw.table.SourceCatalog` + Catalog of difference image sources detected in ``diffIm``. + associatedDiaSources : `pandas.DataFrame` + Catalog of associated difference image sources detected in ``diffIm``. + + Returns + ------- + result : `lsst.pipe.base.Struct` + Results struct with components. + + - ``matchedDiaSources`` : Fakes matched to input diaSources. Has + length of ``fakeCat``. (`pandas.DataFrame`) + """ + # First match the diaSrc to the injected fakes + injectedCat = injectedCat.to_pandas() + nPossibleFakes = len(injectedCat) + + fakeVects = self._getVectors( + np.radians(injectedCat.ra), + np.radians(injectedCat.dec)) + diaSrcVects = self._getVectors( + diaSources['coord_ra'], + diaSources['coord_dec']) + + diaSrcTree = cKDTree(diaSrcVects) + dist, idxs = diaSrcTree.query( + fakeVects, + distance_upper_bound=np.radians(self.config.matchDistanceArcseconds / 3600)) + nFakesFound = np.isfinite(dist).sum() + + self.log.info("Found %d out of %d possible in diaSources.", nFakesFound, nPossibleFakes) + + # assign diaSourceId to the matched fakes + diaSrcIds = diaSources['id'][np.where(np.isfinite(dist), idxs, 0)] + matchedFakes = injectedCat.assign(diaSourceId=np.where(np.isfinite(dist), diaSrcIds, 0)) + matchedFakes['dist_diaSrc'] = np.where(np.isfinite(dist), 3600*np.rad2deg(dist), -1) + + return Struct(matchDiaSources=matchedFakes) + + def _getVectors(self, ras, decs): + """Convert ra dec to unit vectors on the sphere. + + Parameters + ---------- + ras : `numpy.ndarray`, (N,) + RA coordinates in radians. + decs : `numpy.ndarray`, (N,) + Dec coordinates in radians. + + Returns + ------- + vectors : `numpy.ndarray`, (N, 3) + Vectors on the unit sphere for the given RA/DEC values. + """ + vectors = np.empty((len(ras), 3)) + + vectors[:, 2] = np.sin(decs) + vectors[:, 0] = np.cos(decs) * np.cos(ras) + vectors[:, 1] = np.cos(decs) * np.sin(ras) + + return vectors + + def _addPixCoords(self, fakeCat, image): + """Add pixel coordinates to the catalog of fakes. + + Parameters + ---------- + fakeCat : `astropy.table.table.Table` + The catalog of fake sources to be input + image : `lsst.afw.image.exposure.exposure.ExposureF` + The image into which the fake sources should be added + Returns + ------- + fakeCat : `astropy.table.table.Table` + """ + + wcs = image.getWcs() + + # Get x/y pixel coordinates for injected sources. + xs, ys = wcs.skyToPixelArray( + fakeCat["ra"], + fakeCat["dec"], + degrees=True + ) + fakeCat["x"] = xs + fakeCat["y"] = ys + + return fakeCat + + def _trimFakeCat(self, fakeCat, image): + """Trim the fake cat to the exact size of the input image. + + Parameters + ---------- + fakeCat : `astropy.table.table.Table` + The catalog of fake sources that was input + image : `lsst.afw.image.exposure.exposure.ExposureF` + The image into which the fake sources were added + Returns + ------- + fakeCat : `astropy.table.table.Table` + The original fakeCat trimmed to the area of the image + """ + + # fakeCat must be processed with _addPixCoords before trimming + fakeCat = self._addPixCoords(fakeCat, image) + + # Prefilter in ra/dec to avoid cases where the wcs incorrectly maps + # input fakes which are really off the chip onto it. + ras = fakeCat["ra"] * u.deg + decs = fakeCat["dec"] * u.deg + + isContainedRaDec = image.containsSkyCoords(ras, decs, padding=0) + + # now use the exact pixel BBox to filter to only fakes that were inserted + xs = fakeCat["x"] + ys = fakeCat["y"] + + bbox = lsstGeom.Box2D(image.getBBox()) + isContainedXy = xs - self.config.trimBuffer >= bbox.minX + isContainedXy &= xs + self.config.trimBuffer <= bbox.maxX + isContainedXy &= ys - self.config.trimBuffer >= bbox.minY + isContainedXy &= ys + self.config.trimBuffer <= bbox.maxY + + return fakeCat[isContainedRaDec & isContainedXy] + + +class MatchInjectedToAssocDiaSourceConnections( + PipelineTaskConnections, + defaultTemplates={"coaddName": "deep", + "fakesType": "fakes_"}, + dimensions=("instrument", + "visit", + "detector")): + + assocDiaSources = connTypes.Input( + doc="An assocDiaSource catalog to match against fakeCat from the" + "diaPipe run. Assumed to be SDMified.", + name="{fakesType}{coaddName}Diff_assocDiaSrc", + storageClass="DataFrame", + dimensions=("instrument", "visit", "detector"), + ) + matchDiaSources = connTypes.Input( + doc="A catalog of those fakeCat sources that have a match in " + "diaSrc. The schema is the union of the schemas for " + "``fakeCat`` and ``diaSrc``.", + name="{fakesType}{coaddName}Diff_matchDiaSrc", + storageClass="DataFrame", + dimensions=("instrument", "visit", "detector"), + ) + matchAssocDiaSources = connTypes.Output( + doc="A catalog of those fakeCat sources that have a match in " + "associatedDiaSources. The schema is the union of the schemas for " + "``fakeCat`` and ``associatedDiaSources``.", + name="{fakesType}{coaddName}Diff_matchAssocDiaSrc", + storageClass="DataFrame", + dimensions=("instrument", "visit", "detector"), + ) + + +class MatchInjectedToAssocDiaSourceConfig( + VisitInjectConfig, + pipelineConnections=MatchInjectedToAssocDiaSourceConnections): + """Config for MatchFakesTask. + """ + + +class MatchInjectedToAssocDiaSourceTask(PipelineTask): + + _DefaultName = "matchInjectedToAssocDiaSource" + ConfigClass = MatchInjectedToAssocDiaSourceConfig + + def run(self, assocDiaSources, matchDiaSources): + """Tag matched injected sources to associated diaSources. + + Parameters + ---------- + matchDiaSources : `pandas.DataFrame` + Catalog of matched diaSrc to injected sources + assocDiaSources : `pandas.DataFrame` + Catalog of associated difference image sources detected in ``diffIm``. + Returns + ------- + result : `lsst.pipe.base.Struct` + Results struct with components. + + - ``matchAssocDiaSources`` : Fakes matched to associated diaSources. Has + length of ``matchDiaSources``. (`pandas.DataFrame`) + """ + # Match the fakes to the associated sources. For this we don't use the coordinates + # but instead check for the diaSources. Since they were present in the table already + nPossibleFakes = len(matchDiaSources) + matchDiaSources['isAssocDiaSource'] = matchDiaSources.diaSourceId.isin(assocDiaSources.diaSourceId) + assocNFakesFound = matchDiaSources.isAssocDiaSource.sum() + self.log.info("Found %d out of %d possible in assocDiaSources."%(assocNFakesFound, nPossibleFakes)) + + return Struct( + matchAssocDiaSources=matchDiaSources.merge( + assocDiaSources.reset_index(drop=True), + on="diaSourceId", + how="left", + suffixes=('_ssi', '_diaSrc') + ) + ) diff --git a/python/lsst/pipe/tasks/matchFakes.py b/python/lsst/pipe/tasks/matchFakes.py deleted file mode 100644 index c9cd68d12..000000000 --- a/python/lsst/pipe/tasks/matchFakes.py +++ /dev/null @@ -1,445 +0,0 @@ -# This file is part of pipe_tasks. -# -# Developed for the LSST Data Management System. -# This product includes software developed by the LSST Project -# (https://www.lsst.org). -# See the COPYRIGHT file at the top-level directory of this distribution -# for details of code ownership. -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -__all__ = ["MatchFakesTask", - "MatchFakesConfig", - "MatchVariableFakesConfig", - "MatchVariableFakesTask"] - -import astropy.units as u -import numpy as np -import pandas as pd -from scipy.spatial import cKDTree - -from lsst.geom import Box2D -import lsst.pex.config as pexConfig -from lsst.pipe.base import PipelineTask, PipelineTaskConnections, Struct -import lsst.pipe.base.connectionTypes as connTypes -from lsst.skymap import BaseSkyMap - -from lsst.pipe.tasks.insertFakes import InsertFakesConfig - -from deprecated.sphinx import deprecated - - -class MatchFakesConnections(PipelineTaskConnections, - defaultTemplates={"coaddName": "deep", - "fakesType": "fakes_"}, - dimensions=("instrument", - "visit", - "detector")): - skyMap = connTypes.Input( - doc="Input definition of geometry/bbox and projection/wcs for " - "template exposures. Needed to test which tract to generate ", - name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, - dimensions=("skymap",), - storageClass="SkyMap", - ) - fakeCats = connTypes.Input( - doc="Catalog of fake sources inserted into an image.", - name="{fakesType}fakeSourceCat", - storageClass="DataFrame", - dimensions=("tract", "skymap"), - deferLoad=True, - multiple=True - ) - diffIm = connTypes.Input( - doc="Difference image on which the DiaSources were detected.", - name="{fakesType}{coaddName}Diff_differenceExp", - storageClass="ExposureF", - dimensions=("instrument", "visit", "detector"), - ) - associatedDiaSources = connTypes.Input( - doc="A DiaSource catalog to match against fakeCat. Assumed " - "to be SDMified.", - name="{fakesType}{coaddName}Diff_assocDiaSrc", - storageClass="DataFrame", - dimensions=("instrument", "visit", "detector"), - ) - matchedDiaSources = connTypes.Output( - doc="A catalog of those fakeCat sources that have a match in " - "associatedDiaSources. The schema is the union of the schemas for " - "``fakeCat`` and ``associatedDiaSources``.", - name="{fakesType}{coaddName}Diff_matchDiaSrc", - storageClass="DataFrame", - dimensions=("instrument", "visit", "detector"), - ) - - -@deprecated( - reason="This task will be removed in v28.0 as it is replaced by `source_injection` tasks.", - version="v28.0", - category=FutureWarning, -) -class MatchFakesConfig( - InsertFakesConfig, - pipelineConnections=MatchFakesConnections): - """Config for MatchFakesTask. - """ - matchDistanceArcseconds = pexConfig.RangeField( - doc="Distance in arcseconds to match within.", - dtype=float, - default=0.5, - min=0, - max=10, - ) - - doMatchVisit = pexConfig.Field( - dtype=bool, - default=False, - doc="Match visit to trim the fakeCat" - ) - - trimBuffer = pexConfig.Field( - doc="Size of the pixel buffer surrounding the image. Only those fake sources with a centroid" - "falling within the image+buffer region will be considered matches.", - dtype=int, - default=100, - ) - - -@deprecated( - reason="This task will be removed in v28.0 as it is replaced by `source_injection` tasks.", - version="v28.0", - category=FutureWarning, -) -class MatchFakesTask(PipelineTask): - """Match a pre-existing catalog of fakes to a catalog of detections on - a difference image. - - This task is generally for injected sources that cannot be easily - identified by their footprints such as in the case of detector sources - post image differencing. - """ - - _DefaultName = "matchFakes" - ConfigClass = MatchFakesConfig - - def run(self, fakeCats, skyMap, diffIm, associatedDiaSources): - """Compose fakes into a single catalog and match fakes to detected - diaSources within a difference image bound. - - Parameters - ---------- - fakeCats : `pandas.DataFrame` - List of catalog of fakes to match to detected diaSources. - skyMap : `lsst.skymap.SkyMap` - SkyMap defining the tracts and patches the fakes are stored over. - diffIm : `lsst.afw.image.Exposure` - Difference image where ``associatedDiaSources`` were detected. - associatedDiaSources : `pandas.DataFrame` - Catalog of difference image sources detected in ``diffIm``. - - Returns - ------- - result : `lsst.pipe.base.Struct` - Results struct with components. - - - ``matchedDiaSources`` : Fakes matched to input diaSources. Has - length of ``fakeCat``. (`pandas.DataFrame`) - """ - fakeCat = self.composeFakeCat(fakeCats, skyMap) - - if self.config.doMatchVisit: - fakeCat = self.getVisitMatchedFakeCat(fakeCat, diffIm) - - return self._processFakes(fakeCat, diffIm, associatedDiaSources) - - def _processFakes(self, fakeCat, diffIm, associatedDiaSources): - """Match fakes to detected diaSources within a difference image bound. - - Parameters - ---------- - fakeCat : `pandas.DataFrame` - Catalog of fakes to match to detected diaSources. - diffIm : `lsst.afw.image.Exposure` - Difference image where ``associatedDiaSources`` were detected. - associatedDiaSources : `pandas.DataFrame` - Catalog of difference image sources detected in ``diffIm``. - - Returns - ------- - result : `lsst.pipe.base.Struct` - Results struct with components. - - - ``matchedDiaSources`` : Fakes matched to input diaSources. Has - length of ``fakeCat``. (`pandas.DataFrame`) - """ - trimmedFakes = self._trimFakeCat(fakeCat, diffIm) - nPossibleFakes = len(trimmedFakes) - - fakeVects = self._getVectors(trimmedFakes[self.config.ra_col], - trimmedFakes[self.config.dec_col]) - diaSrcVects = self._getVectors( - np.radians(associatedDiaSources.loc[:, "ra"]), - np.radians(associatedDiaSources.loc[:, "dec"])) - - diaSrcTree = cKDTree(diaSrcVects) - dist, idxs = diaSrcTree.query( - fakeVects, - distance_upper_bound=np.radians(self.config.matchDistanceArcseconds / 3600)) - nFakesFound = np.isfinite(dist).sum() - - self.log.info("Found %d out of %d possible.", nFakesFound, nPossibleFakes) - diaSrcIds = associatedDiaSources.iloc[np.where(np.isfinite(dist), idxs, 0)]["diaSourceId"].to_numpy() - matchedFakes = trimmedFakes.assign(diaSourceId=np.where(np.isfinite(dist), diaSrcIds, 0)) - - return Struct( - matchedDiaSources=matchedFakes.merge( - associatedDiaSources.reset_index(drop=True), on="diaSourceId", how="left") - ) - - def composeFakeCat(self, fakeCats, skyMap): - """Concatenate the fakeCats from tracts that may cover the exposure. - - Parameters - ---------- - fakeCats : `list` of `lsst.daf.butler.DeferredDatasetHandle` - Set of fake cats to concatenate. - skyMap : `lsst.skymap.SkyMap` - SkyMap defining the geometry of the tracts and patches. - - Returns - ------- - combinedFakeCat : `pandas.DataFrame` - All fakes that cover the inner polygon of the tracts in this - quantum. - """ - if len(fakeCats) == 1: - return fakeCats[0].get() - outputCat = [] - for fakeCatRef in fakeCats: - cat = fakeCatRef.get() - tractId = fakeCatRef.dataId["tract"] - # Make sure all data is within the inner part of the tract. - outputCat.append(cat[ - skyMap.findTractIdArray(cat[self.config.ra_col], - cat[self.config.dec_col], - degrees=False) - == tractId]) - - return pd.concat(outputCat) - - def getVisitMatchedFakeCat(self, fakeCat, exposure): - """Trim the fakeCat to select particular visit - - Parameters - ---------- - fakeCat : `pandas.core.frame.DataFrame` - The catalog of fake sources to add to the exposure - exposure : `lsst.afw.image.exposure.exposure.ExposureF` - The exposure to add the fake sources to - - Returns - ------- - movingFakeCat : `pandas.DataFrame` - All fakes that belong to the visit - """ - selected = exposure.getInfo().getVisitInfo().getId() == fakeCat["visit"] - - return fakeCat[selected] - - def _addPixCoords(self, fakeCat, image): - - """Add pixel coordinates to the catalog of fakes. - - Parameters - ---------- - fakeCat : `pandas.core.frame.DataFrame` - The catalog of fake sources to be input - image : `lsst.afw.image.exposure.exposure.ExposureF` - The image into which the fake sources should be added - Returns - ------- - fakeCat : `pandas.core.frame.DataFrame` - """ - wcs = image.getWcs() - ras = fakeCat[self.config.ra_col].values - decs = fakeCat[self.config.dec_col].values - xs, ys = wcs.skyToPixelArray(ras, decs) - fakeCat["x"] = xs - fakeCat["y"] = ys - - return fakeCat - - def _trimFakeCat(self, fakeCat, image): - """Trim the fake cat to the exact size of the input image. - - Parameters - ---------- - fakeCat : `pandas.core.frame.DataFrame` - The catalog of fake sources that was input - image : `lsst.afw.image.exposure.exposure.ExposureF` - The image into which the fake sources were added - Returns - ------- - fakeCat : `pandas.core.frame.DataFrame` - The original fakeCat trimmed to the area of the image - """ - - # fakeCat must be processed with _addPixCoords before trimming - if ('x' not in fakeCat.columns) or ('y' not in fakeCat.columns): - fakeCat = self._addPixCoords(fakeCat, image) - - # Prefilter in ra/dec to avoid cases where the wcs incorrectly maps - # input fakes which are really off the chip onto it. - ras = fakeCat[self.config.ra_col].values * u.rad - decs = fakeCat[self.config.dec_col].values * u.rad - - isContainedRaDec = image.containsSkyCoords(ras, decs, padding=0) - - # now use the exact pixel BBox to filter to only fakes that were inserted - xs = fakeCat["x"].values - ys = fakeCat["y"].values - - bbox = Box2D(image.getBBox()) - isContainedXy = xs >= bbox.minX - isContainedXy &= xs <= bbox.maxX - isContainedXy &= ys >= bbox.minY - isContainedXy &= ys <= bbox.maxY - - return fakeCat[isContainedRaDec & isContainedXy] - - def _getVectors(self, ras, decs): - """Convert ra dec to unit vectors on the sphere. - - Parameters - ---------- - ras : `numpy.ndarray`, (N,) - RA coordinates in radians. - decs : `numpy.ndarray`, (N,) - Dec coordinates in radians. - - Returns - ------- - vectors : `numpy.ndarray`, (N, 3) - Vectors on the unit sphere for the given RA/DEC values. - """ - vectors = np.empty((len(ras), 3)) - - vectors[:, 2] = np.sin(decs) - vectors[:, 0] = np.cos(decs) * np.cos(ras) - vectors[:, 1] = np.cos(decs) * np.sin(ras) - - return vectors - - -class MatchVariableFakesConnections(MatchFakesConnections): - ccdVisitFakeMagnitudes = connTypes.Input( - doc="Catalog of fakes with magnitudes scattered for this ccdVisit.", - name="{fakesType}ccdVisitFakeMagnitudes", - storageClass="DataFrame", - dimensions=("instrument", "visit", "detector"), - ) - - -@deprecated( - reason="This task will be removed in v28.0 as it is replaced by `source_injection` tasks.", - version="v28.0", - category=FutureWarning, -) -class MatchVariableFakesConfig(MatchFakesConfig, - pipelineConnections=MatchVariableFakesConnections): - """Config for MatchFakesTask. - """ - pass - - -@deprecated( - reason="This task will be removed in v28.0 as it is replaced by `source_injection` tasks.", - version="v28.0", - category=FutureWarning, -) -class MatchVariableFakesTask(MatchFakesTask): - """Match injected fakes to their detected sources in the catalog and - compute their expected brightness in a difference image assuming perfect - subtraction. - - This task is generally for injected sources that cannot be easily - identified by their footprints such as in the case of detector sources - post image differencing. - """ - _DefaultName = "matchVariableFakes" - ConfigClass = MatchVariableFakesConfig - - def runQuantum(self, butlerQC, inputRefs, outputRefs): - inputs = butlerQC.get(inputRefs) - inputs["band"] = butlerQC.quantum.dataId["band"] - - outputs = self.run(**inputs) - butlerQC.put(outputs, outputRefs) - - def run(self, fakeCats, ccdVisitFakeMagnitudes, skyMap, diffIm, associatedDiaSources, band): - """Match fakes to detected diaSources within a difference image bound. - - Parameters - ---------- - fakeCat : `pandas.DataFrame` - Catalog of fakes to match to detected diaSources. - diffIm : `lsst.afw.image.Exposure` - Difference image where ``associatedDiaSources`` were detected in. - associatedDiaSources : `pandas.DataFrame` - Catalog of difference image sources detected in ``diffIm``. - - Returns - ------- - result : `lsst.pipe.base.Struct` - Results struct with components. - - - ``matchedDiaSources`` : Fakes matched to input diaSources. Has - length of ``fakeCat``. (`pandas.DataFrame`) - """ - fakeCat = self.composeFakeCat(fakeCats, skyMap) - self.computeExpectedDiffMag(fakeCat, ccdVisitFakeMagnitudes, band) - return self._processFakes(fakeCat, diffIm, associatedDiaSources) - - def computeExpectedDiffMag(self, fakeCat, ccdVisitFakeMagnitudes, band): - """Compute the magnitude expected in the difference image for this - detector/visit. Modify fakeCat in place. - - Negative magnitudes indicate that the source should be detected as - a negative source. - - Parameters - ---------- - fakeCat : `pandas.DataFrame` - Catalog of fake sources. - ccdVisitFakeMagnitudes : `pandas.DataFrame` - Magnitudes for variable sources in this specific ccdVisit. - band : `str` - Band that this ccdVisit was observed in. - """ - magName = self.config.mag_col % band - magnitudes = fakeCat[magName].to_numpy() - visitMags = ccdVisitFakeMagnitudes["variableMag"].to_numpy() - diffFlux = (visitMags * u.ABmag).to_value(u.nJy) - (magnitudes * u.ABmag).to_value(u.nJy) - diffMag = np.where(diffFlux > 0, - (diffFlux * u.nJy).to_value(u.ABmag), - -(-diffFlux * u.nJy).to_value(u.ABmag)) - - noVisit = ~fakeCat["isVisitSource"] - noTemplate = ~fakeCat["isTemplateSource"] - both = np.logical_and(fakeCat["isVisitSource"], - fakeCat["isTemplateSource"]) - - fakeCat.loc[noVisit, magName] = -magnitudes[noVisit] - fakeCat.loc[noTemplate, magName] = visitMags[noTemplate] - fakeCat.loc[both, magName] = diffMag[both] diff --git a/tests/test_matchFakes.py b/tests/test_matchFakes.py deleted file mode 100644 index 4e878d1ba..000000000 --- a/tests/test_matchFakes.py +++ /dev/null @@ -1,141 +0,0 @@ -# -# This file is part of pipe_tasks. -# -# Developed for the LSST Data Management System. -# This product includes software developed by the LSST Project -# (http://www.lsst.org). -# See the COPYRIGHT file at the top-level directory of this distribution -# for details of code ownership. -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . -# - -import numpy as np -import pandas as pd -import unittest -import uuid - -import lsst.sphgeom as sphgeom -import lsst.geom as geom -import lsst.meas.base.tests as measTests -import lsst.skymap as skyMap -import lsst.utils.tests - -from lsst.pipe.tasks.matchFakes import MatchFakesTask, MatchFakesConfig - - -class TestMatchFakes(lsst.utils.tests.TestCase): - - def setUp(self): - """Create fake data to use in the tests. - """ - self.bbox = geom.Box2I(geom.Point2I(0, 0), - geom.Extent2I(1024, 1153)) - dataset = measTests.TestDataset(self.bbox) - self.exposure = dataset.exposure - - simpleMapConfig = skyMap.discreteSkyMap.DiscreteSkyMapConfig() - simpleMapConfig.raList = [dataset.exposure.getWcs().getSkyOrigin().getRa().asDegrees()] - simpleMapConfig.decList = [dataset.exposure.getWcs().getSkyOrigin().getDec().asDegrees()] - simpleMapConfig.radiusList = [0.1] - - self.simpleMap = skyMap.DiscreteSkyMap(simpleMapConfig) - self.tractId = 0 - bCircle = self.simpleMap.generateTract(self.tractId).getOuterSkyPolygon().getBoundingCircle() - bCenter = sphgeom.LonLat(bCircle.getCenter()) - bRadius = bCircle.getOpeningAngle().asRadians() - targetSources = 10000 - self.sourceDensity = (targetSources - / (bCircle.getArea() * (180 / np.pi) ** 2)) - self.rng = np.random.default_rng(1234) - - self.fakeCat = pd.DataFrame({ - "fakeId": [uuid.uuid4().int & (1 << 64) - 1 for n in range(targetSources)], - # Quick-and-dirty values for testing - "ra": bCenter.getLon().asRadians() + bRadius * (2.0 * self.rng.random(targetSources) - 1.0), - "dec": bCenter.getLat().asRadians() + bRadius * (2.0 * self.rng.random(targetSources) - 1.0), - "isVisitSource": np.concatenate([np.ones(targetSources//2, dtype="bool"), - np.zeros(targetSources - targetSources//2, dtype="bool")]), - "isTemplateSource": np.concatenate([np.zeros(targetSources//2, dtype="bool"), - np.ones(targetSources - targetSources//2, dtype="bool")]), - **{band: self.rng.uniform(20, 30, size=targetSources) - for band in {"u", "g", "r", "i", "z", "y"}}, - "DiskHalfLightRadius": np.ones(targetSources, dtype="float"), - "BulgeHalfLightRadius": np.ones(targetSources, dtype="float"), - "disk_n": np.ones(targetSources, dtype="float"), - "bulge_n": np.ones(targetSources, dtype="float"), - "a_d": np.ones(targetSources, dtype="float"), - "a_b": np.ones(targetSources, dtype="float"), - "b_d": np.ones(targetSources, dtype="float"), - "b_b": np.ones(targetSources, dtype="float"), - "pa_disk": np.ones(targetSources, dtype="float"), - "pa_bulge": np.ones(targetSources, dtype="float"), - "sourceType": targetSources * ["star"], - }) - - self.inExp = np.zeros(len(self.fakeCat), dtype=bool) - bbox = geom.Box2D(self.exposure.getBBox()) - for idx, row in self.fakeCat.iterrows(): - coord = geom.SpherePoint(row["ra"], - row["dec"], - geom.radians) - cent = self.exposure.getWcs().skyToPixel(coord) - self.inExp[idx] = bbox.contains(cent) - - tmpCat = self.fakeCat[self.inExp].iloc[:int(self.inExp.sum() / 2)] - extraColumnData = self.rng.integers(0, 100, size=len(tmpCat)) - self.sourceCat = pd.DataFrame( - data={"ra": np.degrees(tmpCat["ra"]), - "dec": np.degrees(tmpCat["dec"]), - "diaObjectId": np.arange(1, len(tmpCat) + 1, dtype=int), - "band": "g", - "diaSourceId": np.arange(1, len(tmpCat) + 1, dtype=int), - "extraColumn": extraColumnData}) - self.sourceCat.set_index(["diaObjectId", "band", "extraColumn"], - drop=False, - inplace=True) - - def testProcessFakes(self): - """Test the run method. - """ - matchFakesConfig = MatchFakesConfig() - matchFakesConfig.matchDistanceArcseconds = 0.1 - matchFakes = MatchFakesTask(config=matchFakesConfig) - result = matchFakes._processFakes(self.fakeCat, - self.exposure, - self.sourceCat) - self.assertEqual(self.inExp.sum(), len(result.matchedDiaSources)) - self.assertEqual( - len(self.sourceCat), - np.sum(np.isfinite(result.matchedDiaSources["extraColumn"]))) - - def testTrimCat(self): - """Test that the correct number of sources are in the ccd area. - """ - matchTask = MatchFakesTask() - result = matchTask._trimFakeCat(self.fakeCat, self.exposure) - self.assertEqual(len(result), self.inExp.sum()) - - -class MemoryTester(lsst.utils.tests.MemoryTestCase): - pass - - -def setup_module(module): - lsst.utils.tests.init() - - -if __name__ == "__main__": - lsst.utils.tests.init() - unittest.main() diff --git a/tests/test_matchSourceInjected.py b/tests/test_matchSourceInjected.py new file mode 100644 index 000000000..9d2ee58f6 --- /dev/null +++ b/tests/test_matchSourceInjected.py @@ -0,0 +1,188 @@ +# +# This file is part of ap_pipe. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (http://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +import numpy as np +import pandas as pd +import unittest + +from astropy.table import Table +from lsst.afw.table import SourceCatalog, SourceTable +from lsst.afw.image import ExposureF +from lsst.afw.geom import makeSkyWcs +from lsst.geom import SpherePoint, degrees, Point2D, Extent2I +import lsst.utils.tests + +from lsst.pipe.tasks.matchDiffimSourceInjected import ( + MatchInjectedToDiaSourceTask, + MatchInjectedToAssocDiaSourceTask, + MatchInjectedToDiaSourceConfig, + MatchInjectedToAssocDiaSourceConfig, +) + + +class BaseTestMatchInjected(lsst.utils.tests.TestCase): + def setUp(self): + rng = np.random.RandomState(6) + # 0.35 arcsec = 0.5/np.sqrt(2) arcsec (half of the diagonal of a 0.5 arcsec square) + offsetFactor = 1./3600 * 0.35 + + # Create a mock injected catalog + self.injectedCat = Table( + { + "injection_id": [1, 2, 3, 5, 6, 7, 12, 21, 31, 49], + "injection_flag": np.repeat(0, 10), + # random positions with 10 arcmin size + # ra centered at 0, dec centered at -30 + "ra": (1/6.) * rng.random(size=10), + "dec": -30 + (1/6.) * rng.random(size=10), + "mag": 20.25 - 0.5 * rng.random(size=10), # random magnitudes + "source_type": np.repeat("DeltaFunction", 10) + } + ) + + # Create a mock diaSources catalog + schema = SourceTable.makeMinimalSchema() + self.diaSources = SourceCatalog(schema) + for i in range(5): + record = self.diaSources.addNew() + record.setId(100 + i) + record.setCoord( + # Random posisions of diaSources + SpherePoint((1/6.) * rng.random(), -30 + (1/6.) * rng.random(), degrees) + ) + + for i, (ra, dec) in enumerate(self.injectedCat[['ra', 'dec']][:8]): + record = self.diaSources.addNew() + sign = rng.choice([-1, 1], size=2) + record.setId(i + 200) + record.setCoord( + SpherePoint( + ra+sign[0]*rng.random()*offsetFactor, + dec+sign[1]*rng.random()*offsetFactor, + degrees + ) + ) + + # Create a mock difference image + self.diffIm = ExposureF(Extent2I(4096, 4096)) + crpix = Point2D(0, 0) + crval = SpherePoint(0, -30, degrees) + cdMatrix = np.array([[5.19513851e-05, 2.81124812e-07], + [3.25186974e-07, 5.19112119e-05]]) + wcs = makeSkyWcs(crpix, crval, cdMatrix) + self.diffIm.setWcs(wcs) + + # add a fake source outside of the image + self.injectedCatForTrimming = self.injectedCat.copy() + self.injectedCatForTrimming.add_row( + { + 'injection_id': 50, + 'injection_flag': 0, + 'ra': 360. - 0.1/6., + 'dec': -30 - 0.1/6., + 'mag': 20.0, + 'source_type': 'DeltaFunction' + } + ) + + # only 4 injected sources are associated + self.assocDiaSources = pd.DataFrame( + { + "diaSourceId": [101, 102, 103, 201, 202, 205, 207], + "band": np.repeat("r", 7), + "visit": np.repeat(410001, 7), + "detector": np.repeat(0, 7), + "diaObjectId": np.arange(7), + } + ) + + +class TestMatchInjectedToDiaSourceTask(BaseTestMatchInjected): + + def test_run(self): + config = MatchInjectedToDiaSourceConfig() + config.matchDistanceArcseconds = 0.5 + config.doMatchVisit = False + config.doForcedMeasurement = False + + task = MatchInjectedToDiaSourceTask(config=config) + + result = task.run(self.injectedCat, self.diffIm, self.diaSources) + self.assertEqual(len(result.matchDiaSources), len(self.injectedCat)) + self.assertEqual(np.sum(result.matchDiaSources['diaSourceId'] > 0), 8) + self.assertEqual(np.sum(result.matchDiaSources['dist_diaSrc'] > 0), 8) + self.assertEqual( + np.sum(np.abs(result.matchDiaSources['dist_diaSrc']) < config.matchDistanceArcseconds), 8 + ) + + def test_run_trimming(self): + config = MatchInjectedToDiaSourceConfig() + config.matchDistanceArcseconds = 0.5 + config.doMatchVisit = True + config.doForcedMeasurement = False + + task = MatchInjectedToDiaSourceTask(config=config) + result = task.run(self.injectedCatForTrimming, self.diffIm, self.diaSources) + + self.assertEqual(len(result.matchDiaSources), len(self.injectedCatForTrimming) - 1) + self.assertEqual(np.sum(result.matchDiaSources['diaSourceId'] > 0), 8) + self.assertEqual(np.sum(result.matchDiaSources['dist_diaSrc'] > 0), 8) + self.assertEqual( + np.sum(np.abs(result.matchDiaSources['dist_diaSrc']) < config.matchDistanceArcseconds), 8 + ) + + def test_getVectors(self): + config = MatchInjectedToDiaSourceConfig() + config.matchDistanceArcseconds = 0.5 + config.doMatchVisit = False + config.doForcedMeasurement = False + + task = MatchInjectedToDiaSourceTask(config=config) + + ras = np.radians(self.injectedCat['ra']) + decs = np.radians(self.injectedCat['dec']) + vectors = task._getVectors(ras, decs) + self.assertEqual(vectors.shape, (10, 3)) + + +class TestMatchInjectedToAssocDiaSourceTask(BaseTestMatchInjected): + + def test_run(self): + config = MatchInjectedToDiaSourceConfig() + config.matchDistanceArcseconds = 0.5 + config.doMatchVisit = False + config.doForcedMeasurement = False + + task = MatchInjectedToDiaSourceTask(config=config) + + result = task.run(self.injectedCat, self.diffIm, self.diaSources) + + configAssoc = MatchInjectedToAssocDiaSourceConfig() + taskAssoc = MatchInjectedToAssocDiaSourceTask(config=configAssoc) + resultAssoc = taskAssoc.run(self.assocDiaSources, result.matchDiaSources) + self.assertEqual(len(resultAssoc.matchAssocDiaSources), len(self.injectedCat)) + self.assertEqual(np.sum(resultAssoc.matchAssocDiaSources['isAssocDiaSource']), 4) + + +if __name__ == "__main__": + unittest.main() From 39d16821cd3235f018c9c7735a3ec7ba06d76c59 Mon Sep 17 00:00:00 2001 From: BrunoSanchez Date: Tue, 3 Jun 2025 01:00:07 -0700 Subject: [PATCH 2/2] Remove not needed dependency --- python/lsst/pipe/tasks/matchDiffimSourceInjected.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/python/lsst/pipe/tasks/matchDiffimSourceInjected.py b/python/lsst/pipe/tasks/matchDiffimSourceInjected.py index 6be9f9b22..d7656ab40 100644 --- a/python/lsst/pipe/tasks/matchDiffimSourceInjected.py +++ b/python/lsst/pipe/tasks/matchDiffimSourceInjected.py @@ -31,10 +31,9 @@ from lsst.afw import table as afwTable from lsst import geom as lsstGeom import lsst.pex.config as pexConfig -from lsst.pipe.base import PipelineTask, PipelineTaskConnections, Struct +from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct import lsst.pipe.base.connectionTypes as connTypes from lsst.meas.base import ForcedMeasurementTask, ForcedMeasurementConfig -from lsst.source.injection import VisitInjectConfig class MatchInjectedToDiaSourceConnections( @@ -73,7 +72,7 @@ class MatchInjectedToDiaSourceConnections( class MatchInjectedToDiaSourceConfig( - VisitInjectConfig, + PipelineTaskConfig, pipelineConnections=MatchInjectedToDiaSourceConnections): """Config for MatchFakesTask. """ @@ -388,7 +387,7 @@ class MatchInjectedToAssocDiaSourceConnections( class MatchInjectedToAssocDiaSourceConfig( - VisitInjectConfig, + PipelineTaskConfig, pipelineConnections=MatchInjectedToAssocDiaSourceConnections): """Config for MatchFakesTask. """