diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py index 8e5d0c711..d2286ecb5 100644 --- a/python/lsst/pipe/tasks/calibrateImage.py +++ b/python/lsst/pipe/tasks/calibrateImage.py @@ -458,6 +458,13 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali doc="If True, include astrometric errors in the output catalog.", ) + do_full_star_detection = pexConfig.Field( + dtype=bool, + default=True, + doc="If True, run the full star detection/deblend/measurement. Otherwise, " + "use PSF stars.", + ) + background_stats_ignored_pixel_masks = pexConfig.ListField( dtype=str, default=["SAT", "SUSPECT", "SPIKE"], @@ -533,14 +540,18 @@ def setDefaults(self): # shapeHSM/moments for star/galaxy separation. # TODO DM-39203: we can remove aperture correction from this task # once we are using the shape-based star/galaxy code. - self.psf_source_measurement.plugins = ["base_PixelFlags", - "base_SdssCentroid", - "ext_shapeHSM_HsmSourceMoments", - "base_CircularApertureFlux", - "base_GaussianFlux", - "base_PsfFlux", - "base_CompensatedTophatFlux", - ] + self.psf_source_measurement.plugins = [ + "base_PixelFlags", + "base_SdssCentroid", + "ext_shapeHSM_HsmSourceMoments", + "ext_shapeHSM_HsmPsfMoments", + "base_GaussianFlux", + "base_PsfFlux", + "base_CircularApertureFlux", + "base_ClassificationSizeExtendedness", + "base_CompensatedTophatFlux", + ] + self.psf_source_measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments" self.psf_source_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments" # Only measure apertures we need for PSF measurement. self.psf_source_measurement.plugins["base_CircularApertureFlux"].radii = [12.0] @@ -560,16 +571,17 @@ def setDefaults(self): self.star_detection.reEstimateBackground = False self.star_detection.thresholdValue = 5.0 self.star_detection.includeThresholdMultiplier = 2.0 - self.star_measurement.plugins = ["base_PixelFlags", - "base_SdssCentroid", - "ext_shapeHSM_HsmSourceMoments", - "ext_shapeHSM_HsmPsfMoments", - "base_GaussianFlux", - "base_PsfFlux", - "base_CircularApertureFlux", - "base_ClassificationSizeExtendedness", - "base_CompensatedTophatFlux", - ] + self.star_measurement.plugins = [ + "base_PixelFlags", + "base_SdssCentroid", + "ext_shapeHSM_HsmSourceMoments", + "ext_shapeHSM_HsmPsfMoments", + "base_GaussianFlux", + "base_PsfFlux", + "base_CircularApertureFlux", + "base_ClassificationSizeExtendedness", + "base_CompensatedTophatFlux", + ] self.star_measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments" self.star_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments" # Only measure the apertures we need for star selection. @@ -1076,18 +1088,29 @@ def run( ) # Run the stars_detection subtask for the photometric calibration. - result.stars_footprints = self._find_stars( - result.exposure, - result.background, - id_generator, - background_to_photometric_ratio=result.background_to_photometric_ratio, - adaptive_det_res_struct=adaptive_det_res_struct, - num_astrometry_matches=num_astrometry_matches, - ) + if self.config.do_full_star_detection: + result.stars_footprints = self._find_and_measure_stars( + result.exposure, + result.background, + id_generator, + background_to_photometric_ratio=result.background_to_photometric_ratio, + adaptive_det_res_struct=adaptive_det_res_struct, + num_astrometry_matches=num_astrometry_matches, + ) + else: + result.stars_footprints = self._measure_psf_stars( + result.exposure, + result.psf_stars_footprints, + id_generator, + num_astrometry_matches=num_astrometry_matches, + ) + psf = result.exposure.getPsf() psfSigma = psf.computeShape(result.exposure.getBBox().getCenter()).getDeterminantRadius() - self._match_psf_stars(result.psf_stars_footprints, result.stars_footprints, - psfSigma=psfSigma) + + if self.config.do_full_star_detection: + self._match_psf_stars(result.psf_stars_footprints, result.stars_footprints, + psfSigma=psfSigma) # Update the "stars" source coordinates with the current wcs. afwTable.updateSourceCoords( @@ -1161,12 +1184,16 @@ def run( stdev_bg, _ = statObj.getResult(afwMath.STDEV) self.metadata["bg_subtracted_skyPixel_instFlux_stdev"] = stdev_bg - self.metadata["bg_subtracted_skySource_flux_median"] = ( - np.median(result.stars[result.stars['sky_source']][self.config.background_stats_flux_column]) - ) - self.metadata["bg_subtracted_skySource_flux_stdev"] = ( - np.std(result.stars[result.stars['sky_source']][self.config.background_stats_flux_column]) - ) + if self.config.do_full_star_detection: + self.metadata["bg_subtracted_skySource_flux_median"] = ( + np.median(result.stars[result.stars['sky_source']][self.config.background_stats_flux_column]) + ) + self.metadata["bg_subtracted_skySource_flux_stdev"] = ( + np.std(result.stars[result.stars['sky_source']][self.config.background_stats_flux_column]) + ) + else: + self.metadata["bg_subtracted_skySource_flux_median"] = np.nan + self.metadata["bg_subtracted_skySource_flux_stdev"] = np.nan if self.config.do_calibrate_pixels: self._apply_photometry( @@ -1415,8 +1442,15 @@ def _measure_aperture_correction(self, exposure, bright_sources): exposure.info.setApCorrMap(ap_corr_map) - def _find_stars(self, exposure, background, id_generator, background_to_photometric_ratio=None, - adaptive_det_res_struct=None, num_astrometry_matches=None): + def _find_and_measure_stars( + self, + exposure, + background, + id_generator, + background_to_photometric_ratio=None, + adaptive_det_res_struct=None, + num_astrometry_matches=None, + ): """Detect stars on an exposure that has a PSF model, and measure their PSF, circular aperture, compensated gaussian fluxes. @@ -1482,7 +1516,7 @@ def _find_stars(self, exposure, background, id_generator, background_to_photomet nPeak += nPeakSrc minIsolated = min(400, max(3, 0.005*nPeak, 0.6*num_astrometry_matches)) if nFootprint > 0: - self.log.info("Adaptive threshold detection _find_stars nIter: %d; " + self.log.info("Adaptive threshold detection _find_and_measure_stars nIter: %d; " "nPeak/nFootprint = %.2f (max is 800); nIsolated = %d (min is %.1f).", adaptiveDetIter, nPeak/nFootprint, nIsolated, minIsolated) if nPeak/nFootprint > 800 or nIsolated < minIsolated: @@ -1566,6 +1600,58 @@ def _find_stars(self, exposure, background, id_generator, background_to_photomet else: return result.sourceCat + def _measure_psf_stars(self, exposure, psf_stars, id_generator, num_astrometry_matches=None): + """Do measurements on PSF stars only. + + Parameters + ---------- + exposure : `lsst.afw.image.Exposure` + Exposure to measure stars on. + psf_stars : `lsst.afw.table.SourceCatalog` + Catalog of PSF stars. + id_generator : `lsst.meas.base.IdGenerator` + Object that generates source IDs and provides random seeds. + num_astrometry_matches : `int` + Number of astrometry matches (for metadata). + + Returns + ------- + stars : `SourceCatalog` + PSF stars, with a limited set of measurements performed on them. + """ + # Create the sources table with the correct schema. + sources = afwTable.SourceCatalog(self.initial_stars_schema.schema) + sources.resize(len(psf_stars)) + + # Copy columns. + sourcesNames = sources.schema.getNames() + for column in psf_stars.schema.getNames(): + if column in sourcesNames: + sources[column] = psf_stars[column] + + # Copy footprints. + for i, row in enumerate(psf_stars): + sources[i].setFootprint(row.getFootprint()) + + self.metadata["post_deblend_source_count"] = len(sources) + self.metadata["failed_deblend_source_count"] = 0 + self.metadata["saturated_source_count"] = np.sum(sources["base_PixelFlags_flag_saturated"]) + self.metadata["bad_source_count"] = np.sum(sources["base_PixelFlags_flag_bad"]) + + # We do not need to run star_normalized_calibration_flux because it + # was already applied. + + self.star_apply_aperture_correction.run(sources, exposure.apCorrMap) + self.star_catalog_calculation.run(sources) + self.star_set_primary_flags.run(sources) + + result = self.star_selector.run(sources) + # The star selector may not produce a contiguous catalog. + if not result.sourceCat.isContiguous(): + return result.sourceCat.copy(deep=True) + else: + return result.sourceCat + def _match_psf_stars(self, psf_stars, stars, psfSigma=None): """Match calibration stars to psf stars, to identify which were psf candidates, and which were used or reserved during psf measurement diff --git a/tests/test_calibrateImage.py b/tests/test_calibrateImage.py index 5797c105a..c08ffa04d 100644 --- a/tests/test_calibrateImage.py +++ b/tests/test_calibrateImage.py @@ -166,7 +166,8 @@ def setUp(self): self.config.star_selector["science"].unresolved.maximum = 0.2 def _check_run(self, calibrate, result, expect_calibrated_pixels: bool = True, - expect_n_background: int = 4, expect_n_background_equal_or_greater_than: int = -1): + expect_n_background: int = 4, expect_n_background_equal_or_greater_than: int = -1, + expect_n_stars: int = 6): """Test the result of CalibrateImage.run(). Parameters @@ -177,6 +178,13 @@ def _check_run(self, calibrate, result, expect_calibrated_pixels: bool = True, Result of calling calibrate.run(). expect_calibrated_pixels : `bool`, optional Whether to expect image and background pixels to be calibrated. + expect_n_background : `int`, optional + Expected number of background components. + expect_n_background_equal_or_greater_than : `int`, optional + Expect the number of background components to be greater or equal + than this. + expect_n_stars : `int`, optional + Expected number of stars. """ # Background should have 4 elements: 3 from compute_psf and one from # re-estimation during source detection. @@ -218,7 +226,7 @@ def _check_run(self, calibrate, result, expect_calibrated_pixels: bool = True, # Should have detected all S/N >= 10 sources plus 2 sky sources, # whether 1 or 2 snaps. - self.assertEqual(len(result.stars), 6) + self.assertEqual(len(result.stars), expect_n_stars) # Did the psf flags get propagated from the psf_stars catalog? self.assertEqual(result.stars["calib_psf_used"].sum(), 3) @@ -251,6 +259,35 @@ def test_run(self): self._check_run(calibrate, result) + def test_run_no_full(self): + """Test running without full star detection. + """ + config = CalibrateImageTask.ConfigClass() + config.psf_detection.background.approxOrderX = 1 + config.psf_measure_psf.psfDeterminer = 'pca' + config.astrometry.sourceSelector["science"].flags.good = [] + config.astrometry.matcher.numPointsForShape = 3 + config.run_sattle = False + config.do_adaptive_threshold_detection = False + config.psf_detection.reEstimateBackground = True + config.star_detection.reEstimateBackground = True + config.astrometry.magnitudeOutlierRejectionNSigma = 9.0 + config.id_generator.packer.name = "observation" + config.id_generator.packer["observation"].n_observations = 10000 + config.id_generator.packer["observation"].n_detectors = 99 + config.id_generator.n_releases = 8 + config.id_generator.release_id = 2 + config.star_selector["science"].unresolved.maximum = 0.2 + # Extra config. + config.do_full_star_detection = False + + calibrate = CalibrateImageTask(config=config) + calibrate.astrometry.setRefObjLoader(self.ref_loader) + calibrate.photometry.match.setRefObjLoader(self.ref_loader) + result = calibrate.run(exposures=self.exposure) + + self._check_run(calibrate, result, expect_n_background=3, expect_n_stars=3) + def test_run_adaptive_threshold_detection(self): """Test that run() runs with adaptive threshold detection turned on. """ @@ -436,15 +473,19 @@ def test_measure_aperture_correction(self): # from other configured plugins. self.assertGreater(len(self.exposure.apCorrMap), 2) - def test_find_stars(self): - """Test that _find_stars() correctly identifies the S/N>10 stars + def test_find_and_measure_stars(self): + """Test that _find_and_measure_stars() correctly identifies the S/N>10 stars in the image and returns them in the output catalog. """ calibrate = CalibrateImageTask(config=self.config) psf_stars, candidates, _ = calibrate._compute_psf(self.attributes, self.id_generator) calibrate._measure_aperture_correction(self.exposure, psf_stars) - stars = calibrate._find_stars(self.exposure, self.attributes.background, self.id_generator) + stars = calibrate._find_and_measure_stars( + self.exposure, + self.attributes.background, + self.id_generator, + ) # Catalog ids should be very large from this id generator. self.assertTrue(all(stars['id'] > 1000000000)) @@ -464,6 +505,37 @@ def test_find_stars(self): self.assertFloatsAlmostEqual(record.getY(), center[1], rtol=0.01) self.assertFloatsAlmostEqual(record["slot_PsfFlux_instFlux"], flux, rtol=0.01) + def test_measure_psf_stars(self): + """Test that _measure_psf_stars() correctly measures stars + and returns them in the output catalog. + """ + calibrate = CalibrateImageTask(config=self.config) + psf_stars, candidates, _ = calibrate._compute_psf(self.attributes, self.id_generator) + calibrate._measure_aperture_correction(self.exposure, psf_stars) + + stars = calibrate._measure_psf_stars( + self.exposure, + psf_stars, + self.id_generator, + ) + + # Catalog ids should be very large from this id generator. + self.assertTrue(all(stars['id'] > 1000000000)) + + # Background should have 3 elements from compute_psf. + self.assertEqual(len(self.attributes.background), 3) + + # Only 3 psf-like sources with S/N>50 should be in the output catalog, + # plus two sky sources. + self.assertEqual(len(stars), 3) + self.assertTrue(stars.isContiguous()) + # Sort in brightness order, to easily compare with expected positions. + 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) + def test_astrometry(self): """Test that the fitted WCS gives good catalog coordinates. """ @@ -471,7 +543,11 @@ def test_astrometry(self): calibrate.astrometry.setRefObjLoader(self.ref_loader) psf_stars, candidates, _ = calibrate._compute_psf(self.attributes, self.id_generator) calibrate._measure_aperture_correction(self.exposure, psf_stars) - stars = calibrate._find_stars(self.exposure, self.attributes.background, self.id_generator) + stars = calibrate._find_and_measure_stars( + self.exposure, + self.attributes.background, + self.id_generator, + ) calibrate._fit_astrometry(self.exposure, stars) @@ -486,44 +562,68 @@ def test_photometry(self): """Test that the fitted photoCalib matches the one we generated, and that the exposure is calibrated. """ - calibrate = CalibrateImageTask(config=self.config) - calibrate.astrometry.setRefObjLoader(self.ref_loader) - calibrate.photometry.match.setRefObjLoader(self.ref_loader) - psf_stars, candidates, _ = calibrate._compute_psf(self.attributes, self.id_generator) - calibrate._measure_aperture_correction(self.exposure, psf_stars) - stars = calibrate._find_stars(self.exposure, self.attributes.background, self.id_generator) - calibrate._fit_astrometry(self.exposure, stars) - - stars, matches, meta, photoCalib = calibrate._fit_photometry(self.exposure, stars) - calibrate._apply_photometry(self.exposure, self.attributes.background) - - # NOTE: With this test data, PhotoCalTask returns calibrationErr==0, - # so we can't check that the photoCal error has been set. - self.assertFloatsAlmostEqual(photoCalib.getCalibrationMean(), self.photo_calib, rtol=1e-2) - # The exposure should be calibrated by the applied photoCalib, - # and the background should be calibrated to match. - uncalibrated = self.exposure.image.clone() - uncalibrated += self.attributes.background.getImage() - uncalibrated /= self.photo_calib - self.assertFloatsAlmostEqual(uncalibrated.array, self.truth_exposure.image.array, rtol=1e-2) - # 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, 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[~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) + for do_full_star_detection in [True, False]: + attributes = pipeBase.Struct() + attributes.exposure = self.exposure.clone() + attributes.background = None + attributes.background_to_photometric_ratio = None + + calibrate = CalibrateImageTask(config=self.config) + calibrate.astrometry.setRefObjLoader(self.ref_loader) + calibrate.photometry.match.setRefObjLoader(self.ref_loader) + psf_stars_footprints, candidates, _ = calibrate._compute_psf(attributes, self.id_generator) + calibrate._measure_aperture_correction(attributes.exposure, psf_stars_footprints) + + if do_full_star_detection: + stars = calibrate._find_and_measure_stars( + attributes.exposure, + attributes.background, + self.id_generator, + ) + # Some hidden state in this test written before the refactor + # makes it impossible to use the psf stars here. + calibrate._fit_astrometry(attributes.exposure, stars) + else: + # The psf stars work perfectly fine in this mode. + calibrate._fit_astrometry(attributes.exposure, psf_stars_footprints) + + stars = calibrate._measure_psf_stars( + attributes.exposure, + psf_stars_footprints, + self.id_generator, + ) + + stars, matches, meta, photoCalib = calibrate._fit_photometry(attributes.exposure, stars) + + calibrate._apply_photometry(attributes.exposure, attributes.background) + + # NOTE: With this test data, PhotoCalTask returns calibrationErr==0, + # so we can't check that the photoCal error has been set. + self.assertFloatsAlmostEqual(photoCalib.getCalibrationMean(), self.photo_calib, rtol=1e-2) + # The exposure should be calibrated by the applied photoCalib, + # and the background should be calibrated to match. + uncalibrated = attributes.exposure.image.clone() + uncalibrated += attributes.background.getImage() + uncalibrated /= self.photo_calib + self.assertFloatsAlmostEqual(uncalibrated.array, self.truth_exposure.image.array, rtol=1e-2) + # PhotoCalib on the exposure must be identically 1. + self.assertEqual(attributes.exposure.photoCalib.getCalibrationMean(), 1.0) + + # 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[~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 @@ -532,7 +632,11 @@ def test_match_psf_stars(self): calibrate = CalibrateImageTask(config=self.config) psf_stars, candidates, _ = calibrate._compute_psf(self.attributes, self.id_generator) calibrate._measure_aperture_correction(self.exposure, psf_stars) - stars = calibrate._find_stars(self.exposure, self.attributes.background, self.id_generator) + stars = calibrate._find_and_measure_stars( + self.exposure, + self.attributes.background, + self.id_generator, + ) # There should be no psf-related flags set at first. self.assertEqual(stars["calib_psf_candidate"].sum(), 0)