Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DM-41285: Add sky sources to CalibrateImage #868

Merged
merged 1 commit into from
Feb 16, 2024
Merged

DM-41285: Add sky sources to CalibrateImage #868

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 30 additions & 8 deletions python/lsst/pipe/tasks/calibrateImage.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.",
Expand Down Expand Up @@ -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.",
Expand Down Expand Up @@ -317,15 +324,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.
Expand Down Expand Up @@ -378,6 +390,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)
Expand All @@ -396,6 +409,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],
Expand All @@ -411,12 +425,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.

Expand All @@ -427,6 +441,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
-------
Expand Down Expand Up @@ -456,13 +472,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)
Expand Down Expand Up @@ -602,7 +621,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.

Expand All @@ -613,6 +632,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
-------
Expand All @@ -625,6 +646,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)
Expand Down
63 changes: 42 additions & 21 deletions tests/test_calibrateImage.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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

Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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)

Expand All @@ -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)
Expand Down Expand Up @@ -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()
Expand All @@ -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()
Expand All @@ -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.
Expand Down
Loading