From 2b4a1b65e0f1fa1b58cd2a060b33d62645456810 Mon Sep 17 00:00:00 2001 From: John Parejko Date: Wed, 6 Dec 2023 12:11:15 -0800 Subject: [PATCH] Add sky sources to CalibrateImage Add id generator for subtasks to use. --- python/lsst/pipe/tasks/calibrateImage.py | 38 +++++++++++--- tests/test_calibrateImage.py | 63 ++++++++++++++++-------- 2 files changed, 72 insertions(+), 29 deletions(-) diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py index ca20a19b68..ac77ebf66c 100644 --- a/python/lsst/pipe/tasks/calibrateImage.py +++ b/python/lsst/pipe/tasks/calibrateImage.py @@ -28,7 +28,7 @@ import lsst.meas.algorithms 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 +148,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 +195,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.SkyObjectsTask, + 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.", @@ -319,15 +326,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. @@ -380,6 +392,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) @@ -398,6 +411,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], @@ -413,12 +427,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. @@ -429,6 +443,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 ------- @@ -458,13 +474,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) @@ -604,7 +623,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. @@ -615,6 +634,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 ------- @@ -627,6 +648,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, id_generator.catalog_id, sources) # Mask streaks self.star_mask_streaks.run(exposure) diff --git a/tests/test_calibrateImage.py b/tests/test_calibrateImage.py index a02ca3ed9c..f9ea8afe81 100644 --- a/tests/test_calibrateImage.py +++ b/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,18 +236,19 @@ 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) - self.assertTrue(psf_stars.isContiguous()) + # Only 5 psf-like sources with S/N>10 should be in the output catalog, + # plus two sky sources. + self.assertEqual(len(stars), 7) + self.assertTrue(stars.isContiguous()) # Sort in order of brightness, to easily compare with expected positions. - psf_stars.sort(psf_stars.getPsfFluxSlot().getMeasKey()) - for record, flux, center in zip(psf_stars[::-1], self.fluxes, self.centroids[self.fluxes > 50]): + stars.sort(stars.getPsfFluxSlot().getMeasKey()) + for record, flux, center in zip(stars[::-1], self.fluxes, self.centroids[self.fluxes > 50]): self.assertFloatsAlmostEqual(record.getX(), center[0], rtol=0.01) self.assertFloatsAlmostEqual(record.getY(), center[1], rtol=0.01) self.assertFloatsAlmostEqual(record["slot_PsfFlux_instFlux"], flux, rtol=0.01) @@ -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,21 @@ 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 +317,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 +326,14 @@ 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, with + # two sky sources as the faintest. 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 +353,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 +439,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 +466,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 +491,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.