Add sky sources to CalibrateImage

Add id generator for subtasks to use.

python/lsst/pipe/tasks/calibrateImage.py

@@ -29,6 +29,7 @@
 import lsst.meas.algorithms.installGaussianPsf
 import lsst.meas.algorithms.measureApCorr
 from lsst.meas.algorithms import sourceSelector
+import lsst.meas.base
 import lsst.meas.astrom
 import lsst.meas.deblender
 import lsst.meas.extensions.shapeHSM
@@ -148,6 +149,9 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali
         optional=True
     )
 
+    # To generate catalog ids consistently across subtasks.
+    id_generator = lsst.meas.base.DetectorVisitIdGeneratorConfig.make_field()
+
     snap_combine = pexConfig.ConfigurableField(
         target=snapCombine.SnapCombineTask,
         doc="Task to combine two snaps to make one exposure.",
@@ -192,6 +196,10 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali
         target=lsst.meas.algorithms.SourceDetectionTask,
         doc="Task to detect stars to return in the output catalog."
     )
+    star_sky_sources = pexConfig.ConfigurableField(
+        target=lsst.meas.algorithms.SkySourcesTask,
+        doc="Task to generate sky sources ('empty' regions where there are no detections).",
+    )
     star_mask_streaks = pexConfig.ConfigurableField(
         target=maskStreaks.MaskStreaksTask,
         doc="Task for masking streaks. Adds a STREAK mask plane to an exposure.",
@@ -317,15 +325,20 @@ def setDefaults(self):
         self.star_selector["science"].doSignalToNoise = True
         self.star_selector["science"].doIsolated = True
         self.star_selector["science"].signalToNoise.minimum = 10.0
+        # Keep sky sources in the output catalog, even though they aren't
+        # wanted for calibration.
+        self.star_selector["science"].doSkySources = True
 
         # Use the affine WCS fitter (assumes we have a good camera geometry).
         self.astrometry.wcsFitter.retarget(lsst.meas.astrom.FitAffineWcsTask)
         # phot_g_mean is the primary Gaia band for all input bands.
         self.astrometry_ref_loader.anyFilterMapsToThis = "phot_g_mean"
 
-        # Do not subselect during fitting; we already selected good stars.
-        self.astrometry.sourceSelector = "null"
-        self.photometry.match.sourceSelection.retarget(sourceSelector.NullSourceSelectorTask)
+        # Only reject sky sources; we already selected good stars.
+        self.astrometry.sourceSelector["science"].doFlags = True
+        self.astrometry.sourceSelector["science"].flags.bad = ["sky_source"]
+        self.photometry.match.sourceSelection.doFlags = True
+        self.photometry.match.sourceSelection.flags.bad = ["sky_source"]
 
         # All sources should be good for PSF summary statistics.
         # TODO: These should both be changed to calib_psf_used with DM-41640.
@@ -378,6 +391,7 @@ def __init__(self, initial_stars_schema=None, **kwargs):
 
         afwTable.CoordKey.addErrorFields(initial_stars_schema)
         self.makeSubtask("star_detection", schema=initial_stars_schema)
+        self.makeSubtask("star_sky_sources", schema=initial_stars_schema)
         self.makeSubtask("star_mask_streaks")
         self.makeSubtask("star_deblend", schema=initial_stars_schema)
         self.makeSubtask("star_measurement", schema=initial_stars_schema)
@@ -396,6 +410,7 @@ def __init__(self, initial_stars_schema=None, **kwargs):
 
     def runQuantum(self, butlerQC, inputRefs, outputRefs):
         inputs = butlerQC.get(inputRefs)
+        id_generator = self.config.id_generator.apply(butlerQC.quantum.dataId)
 
         astrometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(
             dataIds=[ref.datasetRef.dataId for ref in inputRefs.astrometry_ref_cat],
@@ -411,12 +426,12 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs):
             config=self.config.photometry_ref_loader, log=self.log)
         self.photometry.match.setRefObjLoader(photometry_loader)
 
-        outputs = self.run(**inputs)
+        outputs = self.run(id_generator=id_generator, **inputs)
 
         butlerQC.put(outputs, outputRefs)
 
     @timeMethod
-    def run(self, *, exposures):
+    def run(self, *, exposures, id_generator=None):
         """Find stars and perform psf measurement, then do a deeper detection
         and measurement and calibrate astrometry and photometry from that.
 
@@ -427,6 +442,8 @@ def run(self, *, exposures):
             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.
+        id_generator : `lsst.meas.base.IdGenerator`, optional
+            Object that generates source IDs and provides random seeds.
 
         Returns
         -------
@@ -456,13 +473,16 @@ def run(self, *, exposures):
                 Reference catalog stars matches used in the photometric fit.
                 (`list` [`lsst.afw.table.ReferenceMatch`] or `lsst.afw.table.BaseCatalog`)
         """
+        if id_generator is None:
+            id_generator = lsst.meas.base.IdGenerator()
+
         exposure = self._handle_snaps(exposures)
 
         psf_stars, background, candidates = self._compute_psf(exposure)
 
         self._measure_aperture_correction(exposure, psf_stars)
 
-        stars = self._find_stars(exposure, background)
+        stars = self._find_stars(exposure, background, id_generator)
 
         astrometry_matches, astrometry_meta = self._fit_astrometry(exposure, stars)
         stars, photometry_matches, photometry_meta, photo_calib = self._fit_photometry(exposure, stars)
@@ -602,7 +622,7 @@ def _measure_aperture_correction(self, exposure, bright_sources):
         result = self.measure_aperture_correction.run(exposure, bright_sources)
         exposure.setApCorrMap(result.apCorrMap)
 
-    def _find_stars(self, exposure, background):
+    def _find_stars(self, exposure, background, id_generator):
         """Detect stars on an exposure that has a PSF model, and measure their
         PSF, circular aperture, compensated gaussian fluxes.
 
@@ -613,6 +633,8 @@ def _find_stars(self, exposure, background):
         background : `lsst.afw.math.BackgroundList`
             Background that was fit to the exposure during detection;
             modified in-place during subsequent detection.
+        id_generator : `lsst.meas.base.IdGenerator`
+            Object that generates source IDs and provides random seeds.
 
         Returns
         -------
@@ -625,6 +647,7 @@ def _find_stars(self, exposure, background):
         # measurement uses the most accurate background-subtraction.
         detections = self.star_detection.run(table=table, exposure=exposure, background=background)
         sources = detections.sources
+        self.star_sky_sources.run(exposure.mask, sources, seed=id_generator.catalog_id)
 
         # Mask streaks
         self.star_mask_streaks.run(exposure)

tests/test_calibrateImage.py

@@ -100,6 +100,8 @@ def setUp(self):
         # We don't have many sources, so have to fit simpler models.
         self.config.psf_detection.background.approxOrderX = 1
         self.config.star_detection.background.approxOrderX = 1
+        # Only insert 2 sky sources, for simplicity.
+        self.config.star_sky_sources.nSources = 2
         # Use PCA psf fitter, as psfex fails if there are only 4 stars.
         self.config.psf_measure_psf.psfDeterminer = 'pca'
         # We don't have many test points, so can't match on complicated shapes.
@@ -111,6 +113,9 @@ def setUp(self):
         # will use those fluxes here, and hopefully can remove this.
         self.config.astrometry.magnitudeOutlierRejectionNSigma = 9.0
 
+        # find_stars needs an id generator.
+        self.id_generator = lsst.meas.base.IdGenerator()
+
         # Something about this test dataset prefers the older fluxRatio here.
         self.config.star_catalog_calculation.plugins['base_ClassificationExtendedness'].fluxRatio = 0.925
 
@@ -231,14 +236,15 @@ def test_find_stars(self):
         psf_stars, background, candidates = calibrate._compute_psf(self.exposure)
         calibrate._measure_aperture_correction(self.exposure, psf_stars)
 
-        stars = calibrate._find_stars(self.exposure, background)
+        stars = calibrate._find_stars(self.exposure, background, self.id_generator)
 
         # Background should have 4 elements: 3 from compute_psf and one from
         # re-estimation during source detection.
         self.assertEqual(len(background), 4)
 
-        # Only psf-like sources with S/N>10 should be in the output catalog.
-        self.assertEqual(len(stars), 5)
+        # Only psf-like sources with S/N>10 should be in the output catalog,
+        # plus two sky sources.
+        self.assertEqual(len(stars), 7)
         self.assertTrue(psf_stars.isContiguous())
         # Sort in order of brightness, to easily compare with expected positions.
         psf_stars.sort(psf_stars.getPsfFluxSlot().getMeasKey())
@@ -254,12 +260,13 @@ def test_astrometry(self):
         calibrate.astrometry.setRefObjLoader(self.ref_loader)
         psf_stars, background, candidates = calibrate._compute_psf(self.exposure)
         calibrate._measure_aperture_correction(self.exposure, psf_stars)
-        stars = calibrate._find_stars(self.exposure, background)
+        stars = calibrate._find_stars(self.exposure, background, self.id_generator)
 
         calibrate._fit_astrometry(self.exposure, stars)
 
         # Check that we got reliable matches with the truth coordinates.
-        fitted = SkyCoord(stars['coord_ra'], stars['coord_dec'], unit="radian")
+        sky = stars["sky_source"]
+        fitted = SkyCoord(stars[~sky]['coord_ra'], stars[~sky]['coord_dec'], unit="radian")
         truth = SkyCoord(self.truth_cat['coord_ra'], self.truth_cat['coord_dec'], unit="radian")
         idx, d2d, _ = fitted.match_to_catalog_sky(truth)
         np.testing.assert_array_less(d2d.to_value(u.milliarcsecond), 35.0)
@@ -273,7 +280,7 @@ def test_photometry(self):
         calibrate.photometry.match.setRefObjLoader(self.ref_loader)
         psf_stars, background, candidates = calibrate._compute_psf(self.exposure)
         calibrate._measure_aperture_correction(self.exposure, psf_stars)
-        stars = calibrate._find_stars(self.exposure, background)
+        stars = calibrate._find_stars(self.exposure, background, self.id_generator)
         calibrate._fit_astrometry(self.exposure, stars)
 
         stars, matches, meta, photoCalib = calibrate._fit_photometry(self.exposure, stars)
@@ -287,15 +294,17 @@ def test_photometry(self):
         # PhotoCalib on the exposure must be identically 1.
         self.assertEqual(self.exposure.photoCalib.getCalibrationMean(), 1.0)
 
-        # Check that we got reliable magnitudes and fluxes vs. truth.
-        fitted = SkyCoord(stars['coord_ra'], stars['coord_dec'], unit="radian")
+        # Check that we got reliable magnitudes and fluxes vs. truth, ignoring
+        # sky sources.
+        sky = stars["sky_source"]
+        fitted = SkyCoord(stars[~sky]['coord_ra'], stars[~sky]['coord_dec'], unit="radian")
         truth = SkyCoord(self.truth_cat['coord_ra'], self.truth_cat['coord_dec'], unit="radian")
         idx, _, _ = fitted.match_to_catalog_sky(truth)
         # Because the input variance image does not include contributions from
         # the sources, we can't use fluxErr as a bound on the measurement
         # quality here.
-        self.assertFloatsAlmostEqual(stars['slot_PsfFlux_flux'], self.truth_cat['truth_flux'][idx], rtol=0.1)
-        self.assertFloatsAlmostEqual(stars['slot_PsfFlux_mag'], self.truth_cat['truth_mag'][idx], rtol=0.01)
+        self.assertFloatsAlmostEqual(stars[~sky]['slot_PsfFlux_flux'], self.truth_cat['truth_flux'][idx], rtol=0.1)
+        self.assertFloatsAlmostEqual(stars[~sky]['slot_PsfFlux_mag'], self.truth_cat['truth_mag'][idx], rtol=0.01)
 
     def test_match_psf_stars(self):
         """Test that _match_psf_stars() flags the correct stars as psf stars
@@ -304,7 +313,7 @@ def test_match_psf_stars(self):
         calibrate = CalibrateImageTask(config=self.config)
         psf_stars, background, candidates = calibrate._compute_psf(self.exposure)
         calibrate._measure_aperture_correction(self.exposure, psf_stars)
-        stars = calibrate._find_stars(self.exposure, background)
+        stars = calibrate._find_stars(self.exposure, background, self.id_generator)
 
         # There should be no psf-related flags set at first.
         self.assertEqual(stars["calib_psf_candidate"].sum(), 0)
@@ -313,12 +322,15 @@ def test_match_psf_stars(self):
 
         calibrate._match_psf_stars(psf_stars, stars)
 
-        # Sort in order of brightness; the psf stars are the 3 brightest.
+        # Sort in order of brightness; the psf stars are the 3 brightest;
+        # first two are sky sources.
         stars.sort(stars.getPsfFluxSlot().getMeasKey())
         # sort() above leaves the catalog non-contiguous.
         stars = stars.copy(deep=True)
-        np.testing.assert_array_equal(stars["calib_psf_candidate"], [False, False, True, True, True])
-        np.testing.assert_array_equal(stars["calib_psf_used"], [False, False, True, True, True])
+        np.testing.assert_array_equal(stars["calib_psf_candidate"],
+                                      [False, False, False, False, True, True, True])
+        np.testing.assert_array_equal(stars["calib_psf_used"],
+                                      [False, False, False, False, True, True, True])
         # Too few sources to reserve any in these tests.
         self.assertEqual(stars["calib_psf_reserved"].sum(), 0)
 
@@ -338,8 +350,14 @@ def setUp(self):
         self.repo_path = tempfile.TemporaryDirectory(ignore_cleanup_errors=True)
         self.repo = butlerTests.makeTestRepo(self.repo_path.name)
 
+        # A complete instrument record is necessary for the id generator.
+        instrumentRecord = self.repo.dimensions["instrument"].RecordClass(
+            name=instrument, visit_max=1e6, exposure_max=1e6, detector_max=128,
+            class_name="lsst.obs.base.instrument_tests.DummyCam",
+        )
+        self.repo.registry.syncDimensionData("instrument", instrumentRecord)
+
         # dataIds for fake data
-        butlerTests.addDataIdValue(self.repo, "instrument", instrument)
         butlerTests.addDataIdValue(self.repo, "exposure", exposure0)
         butlerTests.addDataIdValue(self.repo, "exposure", exposure1)
         butlerTests.addDataIdValue(self.repo, "visit", visit)
@@ -418,7 +436,7 @@ def test_runQuantum(self):
         self.assertEqual(task.astrometry.refObjLoader.name, "gaia_dr3_20230707")
         self.assertEqual(task.photometry.match.refObjLoader.name, "ps1_pv3_3pi_20170110")
         # Check that the proper kwargs are passed to run().
-        self.assertEqual(mock_run.call_args.kwargs.keys(), {"exposures"})
+        self.assertEqual(mock_run.call_args.kwargs.keys(), {"exposures", "id_generator"})
 
     def test_runQuantum_2_snaps(self):
         task = CalibrateImageTask()
@@ -445,7 +463,7 @@ def test_runQuantum_2_snaps(self):
         self.assertEqual(task.astrometry.refObjLoader.name, "gaia_dr3_20230707")
         self.assertEqual(task.photometry.match.refObjLoader.name, "ps1_pv3_3pi_20170110")
         # Check that the proper kwargs are passed to run().
-        self.assertEqual(mock_run.call_args.kwargs.keys(), {"exposures"})
+        self.assertEqual(mock_run.call_args.kwargs.keys(), {"exposures", "id_generator"})
 
     def test_runQuantum_no_optional_outputs(self):
         config = CalibrateImageTask.ConfigClass()
@@ -470,7 +488,7 @@ def test_runQuantum_no_optional_outputs(self):
         self.assertEqual(task.astrometry.refObjLoader.name, "gaia_dr3_20230707")
         self.assertEqual(task.photometry.match.refObjLoader.name, "ps1_pv3_3pi_20170110")
         # Check that the proper kwargs are passed to run().
-        self.assertEqual(mock_run.call_args.kwargs.keys(), {"exposures"})
+        self.assertEqual(mock_run.call_args.kwargs.keys(), {"exposures", "id_generator"})
 
     def test_lintConnections(self):
         """Check that the connections are self-consistent.