diff --git a/python/lsst/drp/tasks/reprocess_visit_image.py b/python/lsst/drp/tasks/reprocess_visit_image.py index d2df3a2b..7089f690 100644 --- a/python/lsst/drp/tasks/reprocess_visit_image.py +++ b/python/lsst/drp/tasks/reprocess_visit_image.py @@ -21,6 +21,8 @@ __all__ = ["ReprocessVisitImageTask", "ReprocessVisitImageConfig", "combine_backgrounds"] +import smatch + import lsst.afw.image as afwImage import lsst.afw.table as afwTable import lsst.meas.algorithms @@ -49,6 +51,13 @@ class ReprocessVisitImageConnections( storageClass="ExposureCatalog", dimensions=["instrument", "visit"], ) + initial_photo_calib = connectionTypes.Input( + doc="Photometric calibration that was applied to exposure during the measurement of background_1." + " Used to uncalibrate the background before subtracting it from the input exposure.", + name="initial_photoCalib_detector", + storageClass="PhotoCalib", + dimensions=("instrument", "visit", "detector"), + ) background_1 = connectionTypes.Input( doc="Background models estimated during calibration.", name="initial_pvi_background", @@ -61,9 +70,13 @@ class ReprocessVisitImageConnections( dimensions=("instrument", "visit", "detector"), storageClass="Background", ) - - # TODO: we'll want the finalized catalog to get calib_psf_* from. - # TODO: eventually we'll want to pull in the STREAK mask from the diffim. + calib_sources = connectionTypes.Input( + doc="Per-visit catalog of measurements to get 'calib_*' flags from.", + name="finalized_src_table", + storageClass="ArrowAstropy", + dimensions=["instrument", "visit"], + ) + # TODO DM-45980: pull in the STREAK mask from the diffim. # outputs source_schema = connectionTypes.InitOutput( @@ -151,11 +164,15 @@ class ReprocessVisitImageConfig( target=lsst.meas.base.CatalogCalculationTask, doc="Task to compute catalog values using only the catalog entries.", ) - # TODO: hopefully we can get rid of this, once Jim sorts out SDM details. post_calculations = pexConfig.ConfigurableField( target=lsst.meas.base.SingleFrameMeasurementTask, doc="Task to compute catalog values after all other calculations have been done.", ) + calib_match_radius = pexConfig.Field( + dtype=float, + default=0.2, + doc="Radius in arcseconds to match calib_sources to the output catalog." + ) def setDefaults(self): super().setDefaults() @@ -252,6 +269,7 @@ def __init__(self, source_schema=None, **kwargs): type="Flag", doc="Set if source was used in the PSF determination by FinalizeCharacterizationTask.", ) + self.psf_fields = ("calib_psf_candidate", "calib_psf_used", "calib_psf_reserved") # These fields are only here to satisfy the SDM schema, and will # be removed from there as they are misleading (because we don't @@ -286,6 +304,8 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): detector = outputRefs.exposure.dataId["detector"] exposures = inputs.pop("exposures") visit_summary = inputs.pop("visit_summary") + calib_sources = inputs.pop("calib_sources") + initial_photo_calib = inputs.pop("initial_photo_calib") background_1 = inputs.pop("background_1") if self.config.do_use_sky_corr: background_2 = inputs.pop("background_2") @@ -297,6 +317,10 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): assert not inputs, "runQuantum got more inputs than expected" detector_summary = visit_summary.find(detector) + if detector_summary is None: + msg = f"Detector {detector} not found in visit summary table; not reprocessing this exposure." + raise pipeBase.NoWorkFound(msg) + # Specify the fields that `annotate` needs below, to ensure they # exist, even as None. result = pipeBase.Struct( @@ -306,12 +330,14 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): try: self.run( exposures=exposures, + initial_photo_calib=initial_photo_calib, psf=detector_summary.psf, background=background, ap_corr=detector_summary.apCorrMap, photo_calib=detector_summary.photoCalib, wcs=detector_summary.wcs, summary_stats=afwImage.ExposureSummaryStats.from_record(detector_summary), + calib_sources=calib_sources, result=result, id_generator=id_generator, ) @@ -328,12 +354,14 @@ def run( self, *, exposures, + initial_photo_calib, psf, background, ap_corr, photo_calib, wcs, summary_stats, + calib_sources, id_generator=None, result=None ): @@ -349,6 +377,9 @@ def run( Modified in-place during processing if only one is passed. If two exposures are passed, treat them as snaps and combine before doing further processing. + initial_photo_calib : `lsst.afw.image.PhotoCalib` + Photometric calibration that was applied to exposure during the + measurement of the background. psf : `lsst.afw.detection.Psf` PSF model for this exposure. background : `lsst.afw.math.BackgroundList` @@ -362,6 +393,8 @@ def run( World Coordinate System model for this exposure. summary_stats : `lsst.afw.image.ExposureSummaryStats` Summary statistics measured on this exposure. + calib_sources : `astropy.table.Table` + Per-visit catalog of measurements to get 'calib_*' flags from. id_generator : `lsst.meas.base.IdGenerator`, optional Object that generates source IDs and provides random seeds. result : `lsst.pipe.base.Struct`, optional @@ -393,11 +426,11 @@ def run( result.exposure = self.snap_combine.run(exposures).exposure # Calibrate the image, so it's on the same units as the background. - result.exposure.maskedImage = photo_calib.calibrateImage(result.exposure.maskedImage) + result.exposure.maskedImage = initial_photo_calib.calibrateImage(result.exposure.maskedImage) result.exposure.maskedImage -= background.getImage() # Uncalibrate so that we do the measurements in instFlux, because # we don't have a way to identify measurements as being in nJy (TODO!). - result.exposure.maskedImage /= photo_calib.getCalibrationMean() + result.exposure.maskedImage /= initial_photo_calib.getCalibrationMean() result.exposure.setPsf(psf) result.exposure.setApCorrMap(ap_corr) @@ -407,9 +440,10 @@ def run( # to calibrate the image and catalog, and thus the image will have a # photoCalib of exactly 1. - result.source_footprints = self._find_sources(result.exposure, background, id_generator) - result.source_footprints["visit"] = result.exposure.visitInfo.id - result.source_footprints["detector"] = result.exposure.info.getDetector().getId() + result.source_footprints = self._find_sources(result.exposure, + background, + calib_sources, + id_generator) result.background = background self._apply_photo_calib(result.exposure, result.source_footprints, photo_calib) @@ -419,16 +453,18 @@ def run( return result - def _find_sources(self, exposure, background, id_generator): + def _find_sources(self, exposure, background, calib_sources, id_generator): """Detect and measure sources on the exposure. Parameters ---------- - exposures : `lsst.afw.image.Exposure` + exposure : `lsst.afw.image.Exposure` Exposure to detect and measure sources on; must have a valid PSF. background : `lsst.afw.math.BackgroundList` Background that was fit to the exposure during detection; modified in-place during subsequent detection. + calib_sources : `astropy.table.Table` + Per-visit catalog of measurements to get 'calib_*' flags from. id_generator : `lsst.meas.base.IdGenerator` Object that generates source IDs and provides random seeds. @@ -456,8 +492,39 @@ def _find_sources(self, exposure, background, id_generator): self.catalog_calculation.run(source) self.set_primary_flags.run(source) + source["visit"] = exposure.visitInfo.id + source["detector"] = exposure.info.getDetector().getId() + + self._match_calib_sources(source, calib_sources, exposure.info.getDetector().getId()) + return source + def _match_calib_sources(self, source, calib_sources, detector): + """Match with calib_sources to set `calib_*` flags in the output + catalog. + + Parameters + ---------- + source : `lsst.afw.table.SourceCatalog` + Catalog that was detected and measured on the exposure. Modified + in place to set the psf_fields. + calib_sources : `astropy.table.Table` + Per-visit catalog of measurements to get 'calib_*' flags from. + detector : `int` + Id of detector for this exposure, to get the correct sources from + calib_sources for cross-matching. + """ + use = calib_sources["detector"] == detector + sky = source["sky_source"] + with smatch.Matcher(source[~sky]["coord_ra"], source[~sky]["coord_dec"]) as matcher: + _, i1, i2, _ = matcher.query_knn(calib_sources[use]["coord_ra"], + calib_sources[use]["coord_dec"], + 1, + self.config.calib_match_radius/3600.0, + return_indices=True) + for field in self.psf_fields: + source[~sky][field] = calib_sources[use][i1][field] + def _apply_photo_calib(self, exposure, source_footprints, photo_calib): """Photometrically calibrate the exposure and catalog with the supplied PhotoCalib, and set the exposure's PhotoCalib to 1. diff --git a/tests/test_reprocess_visit_image.py b/tests/test_reprocess_visit_image.py index 36d2562b..63f7cae3 100644 --- a/tests/test_reprocess_visit_image.py +++ b/tests/test_reprocess_visit_image.py @@ -22,6 +22,9 @@ import tempfile import unittest +import astropy.table +import numpy as np + import lsst.afw.image import lsst.afw.math import lsst.afw.table @@ -31,7 +34,6 @@ import lsst.meas.base.tests import lsst.pipe.base.testUtils import lsst.utils.tests -import numpy as np from lsst.drp.tasks.reprocess_visit_image import ReprocessVisitImageTask @@ -72,10 +74,12 @@ class CalibrateImageTaskTests(lsst.utils.tests.TestCase): def setUp(self): # Different x/y dimensions so they're easy to distinguish in a plot, # and non-zero minimum, to help catch xy0 errors. + detector = 42 bbox = lsst.geom.Box2I(lsst.geom.Point2I(5, 4), lsst.geom.Point2I(205, 184)) self.sky_center = lsst.geom.SpherePoint(245.0, -45.0, lsst.geom.degrees) self.photo_calib = 12.3 - dataset = lsst.meas.base.tests.TestDataset(bbox, crval=self.sky_center, calibration=self.photo_calib) + dataset = lsst.meas.base.tests.TestDataset(bbox, crval=self.sky_center, calibration=self.photo_calib, + detector=detector) # sqrt of area of a normalized 2d gaussian psf_scale = np.sqrt(4 * np.pi * (dataset.psfShape.getDeterminantRadius()) ** 2) noise = 10.0 # stddev of noise per pixel @@ -91,7 +95,7 @@ def setUp(self): ) ) self.centroids = np.array( - ((162, 22), (40, 70), (100, 160), (50, 120), (92, 35), (175, 154)), dtype=np.float32 + ((162, 25), (40, 70), (100, 160), (50, 120), (92, 35), (175, 154)), dtype=np.float32 ) for flux, centroid in zip(self.fluxes, self.centroids): dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(centroid[0], centroid[1])) @@ -104,13 +108,15 @@ def setUp(self): schema = dataset.makeMinimalSchema() self.truth_exposure, self.truth_cat = dataset.realize(noise=noise, schema=schema) - # TODO: add a background gradient? # Make an exposure that looks like a PostISRCCD, to serve as the input. self.exposure = lsst.afw.image.ExposureF() self.exposure.maskedImage = self.truth_exposure.maskedImage.clone() self.exposure.mask.clearMaskPlane(self.exposure.mask.getMaskPlane("DETECTED")) self.exposure.mask.clearMaskPlane(self.exposure.mask.getMaskPlane("DETECTED_NEGATIVE")) + # PostISRCCD will have a VisitInfo and Detector attached. + self.exposure.info.setVisitInfo(self.truth_exposure.visitInfo) + self.exposure.info.setDetector(self.truth_exposure.getDetector()) # Subtract a background from the truth exposure, so that the input # exposure still has the background in it. @@ -123,6 +129,23 @@ def setUp(self): wcs=self.truth_exposure.wcs, photo_calib=self.truth_exposure.photoCalib ) + # A catalog that looks like the output of finalizeCharacterization, + # with a value set that we can test on the output. + self.visit_catalog = self.truth_cat.asAstropy() + self.visit_catalog.add_column(astropy.table.Column(data=np.zeros(len(self.visit_catalog)), + name="calib_psf_used", + dtype=bool)) + self.visit_catalog.add_column(astropy.table.Column(data=np.zeros(len(self.visit_catalog)), + name="calib_psf_reserved", + dtype=bool)) + self.visit_catalog.add_column(astropy.table.Column(data=np.zeros(len(self.visit_catalog)), + name="calib_psf_candidate", + dtype=bool)) + self.visit_catalog.add_column(astropy.table.Column(data=[detector]*len(self.visit_catalog), + name="detector")) + # Marking faintest source, so it's easy to identify later. + self.visit_catalog["calib_psf_used"][0] = True + # Test-specific configuration: self.config = ReprocessVisitImageTask.ConfigClass() # Don't really have a background, so have to fit simpler models. @@ -149,15 +172,17 @@ def setUp(self): self.id_generator = self.config.id_generator.apply(data_id) def test_run(self): - task = ReprocessVisitImageTask() + task = ReprocessVisitImageTask(config=self.config) result = task.run( exposures=self.exposure, + initial_photo_calib=self.truth_exposure.photoCalib, psf=self.truth_exposure.psf, background=self.background, ap_corr=lsst.afw.image.ApCorrMap(), photo_calib=self.truth_exposure.photoCalib, wcs=self.truth_exposure.wcs, summary_stats=lsst.afw.image.ExposureSummaryStats(), + calib_sources=self.visit_catalog, id_generator=self.id_generator, ) @@ -170,6 +195,12 @@ def test_run(self): self.assertNotEqual(result.exposure.photoCalib, self.truth_exposure.photoCalib) self.assertFloatsAlmostEqual(result.exposure.photoCalib.getCalibrationMean(), 1) + # All sources (plus sky sources) should have been detected. + self.assertEqual(len(result.source), len(self.truth_cat) + self.config.sky_sources.nSources) + # Faintest non-sky source should be marked as used. + flux_sorted = result.source[result.source.argsort("slot_CalibFlux_instFlux")] + self.assertTrue(flux_sorted[~flux_sorted["sky_source"]]["calib_psf_used"][0]) + class ReprocessVisitImageTaskRunQuantumTests(lsst.utils.tests.TestCase): """Tests of ``ReprocessVisitImageTask.runQuantum``, which need a test @@ -199,6 +230,7 @@ def setUp(self): # dataIds for fake data butlerTests.addDataIdValue(self.repo, "detector", detector) + butlerTests.addDataIdValue(self.repo, "detector", detector+1) butlerTests.addDataIdValue(self.repo, "exposure", exposure0) butlerTests.addDataIdValue(self.repo, "exposure", exposure1) butlerTests.addDataIdValue(self.repo, "visit", visit) @@ -211,7 +243,11 @@ def setUp(self): butlerTests.addDatasetType( self.repo, "initial_pvi_background", {"instrument", "visit", "detector"}, "Background" ) + butlerTests.addDatasetType( + self.repo, "initial_photoCalib_detector", {"instrument", "visit", "detector"}, "PhotoCalib" + ) butlerTests.addDatasetType(self.repo, "skyCorr", {"instrument", "visit", "detector"}, "Background") + butlerTests.addDatasetType(self.repo, "finalized_src_table", {"instrument", "visit"}, "DataFrame") # outputs butlerTests.addDatasetType( @@ -238,6 +274,10 @@ def setUp(self): self.visit_id = self.repo.registry.expandDataId( {"instrument": instrument, "visit": visit, "detector": detector} ) + # Second id for testing on a detector that is not in visitSummary. + self.visit1_id = self.repo.registry.expandDataId( + {"instrument": instrument, "visit": visit, "detector": detector+1} + ) self.visit_only_id = self.repo.registry.expandDataId({"instrument": instrument, "visit": visit}) # put empty data @@ -246,8 +286,13 @@ def setUp(self): self.butler.put(lsst.afw.image.ExposureF(), "postISRCCD", self.exposure1_id) control = lsst.afw.math.BackgroundControl(10, 10) background = lsst.afw.math.makeBackground(lsst.afw.image.ImageF(100, 100), control) + self.butler.put(lsst.afw.image.PhotoCalib(10), "initial_photoCalib_detector", self.visit_id) + self.butler.put(lsst.afw.image.PhotoCalib(10), "initial_photoCalib_detector", self.visit1_id) self.butler.put(lsst.afw.math.BackgroundList(background), "initial_pvi_background", self.visit_id) + self.butler.put(lsst.afw.math.BackgroundList(background), "initial_pvi_background", self.visit1_id) self.butler.put(lsst.afw.math.BackgroundList(background), "skyCorr", self.visit_id) + self.butler.put(lsst.afw.math.BackgroundList(background), "skyCorr", self.visit1_id) + self.butler.put(lsst.afw.table.SourceCatalog().asAstropy(), "finalized_src_table", self.visit_only_id) self.butler.put(make_visit_summary(detector=detector), "finalVisitSummary", self.visit_only_id) def tearDown(self): @@ -269,8 +314,10 @@ def test_runQuantum(self): { "exposures": [self.exposure0_id], "visit_summary": self.visit_only_id, + "initial_photo_calib": self.visit_id, "background_1": self.visit_id, "background_2": self.visit_id, + "calib_sources": self.visit_only_id, # outputs "exposure": self.visit_id, "source": self.visit_id, @@ -284,16 +331,46 @@ def test_runQuantum(self): mock_run.call_args.kwargs.keys(), { "exposures", - "background", + "initial_photo_calib", "psf", + "background", "ap_corr", "photo_calib", "wcs", "summary_stats", - "result", + "calib_sources", "id_generator", + "result", + }, + ) + + def test_runQuantum_no_detector_in_visit_summary(self): + """Test how the task handles the detector not being in the input visit + summary. + """ + task = ReprocessVisitImageTask() + lsst.pipe.base.testUtils.assertValidInitOutput(task) + + quantum = lsst.pipe.base.testUtils.makeQuantum( + task, + self.butler, + self.visit1_id, + { + "exposures": [self.exposure0_id], + "visit_summary": self.visit_only_id, + "initial_photo_calib": self.visit1_id, + "background_1": self.visit1_id, + "background_2": self.visit1_id, + "calib_sources": self.visit_only_id, + # outputs + "exposure": self.visit1_id, + "source": self.visit1_id, + "source_footprints": self.visit1_id, + "background": self.visit1_id, }, ) + with self.assertRaisesRegex(lsst.pipe.base.NoWorkFound, "Detector 43 not found"): + lsst.pipe.base.testUtils.runTestQuantum(task, self.butler, quantum) def test_runQuantum_no_sky_corr(self): """Test that the task will run if using the sky_corr input is @@ -311,7 +388,9 @@ def test_runQuantum_no_sky_corr(self): { "exposures": [self.exposure0_id], "visit_summary": self.visit_only_id, + "initial_photo_calib": self.visit_id, "background_1": self.visit_id, + "calib_sources": self.visit_only_id, # outputs "exposure": self.visit_id, "source": self.visit_id, @@ -325,14 +404,16 @@ def test_runQuantum_no_sky_corr(self): mock_run.call_args.kwargs.keys(), { "exposures", - "background", + "initial_photo_calib", "psf", + "background", "ap_corr", "photo_calib", "wcs", "summary_stats", - "result", + "calib_sources", "id_generator", + "result", }, ) diff --git a/ups/drp_tasks.table b/ups/drp_tasks.table index c6e08565..6e6d438c 100644 --- a/ups/drp_tasks.table +++ b/ups/drp_tasks.table @@ -14,12 +14,7 @@ setupRequired(geom) setupRequired(pipe_base) setupRequired(utils) setupRequired(meas_algorithms) -#setupRequired(meas_extensions_convolved) -#setupRequired(meas_extensions_gaap) setupRequired(meas_extensions_photometryKron) -#setupRequired(meas_extensions_piff) -#setupRequired(meas_extensions_psfex) -#setupRequired(meas_extensions_scarlet) setupRequired(meas_extensions_shapeHSM) setupRequired(pipe_tasks)