Skip to content

Commit

Permalink
Add sky sources to CalibrateImage
Browse files Browse the repository at this point in the history 
Add id generator for subtasks to use.
  • Loading branch information
@parejkoj
parejkoj committed Dec 6, 2023
1 parent 97fd7b4 commit fa72e1a
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 22 deletions.
38 changes: 30 additions & 8 deletions python/lsst/pipe/tasks/calibrateImage.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,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 @@ -144,6 +144,9 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali
optional=True
)

# To generate catalog ids consistently across subtasks.
id_generator = lsst.meas.base.DetectorVisitIdGeneratorConfig.make_field()

# subtasks used during psf characterization
install_simple_psf = pexConfig.ConfigurableField(
target=lsst.meas.algorithms.installGaussianPsf.InstallGaussianPsfTask,
Expand Down Expand Up @@ -183,6 +186,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.",
Expand Down Expand Up @@ -308,15 +315,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 @@ -356,6 +368,7 @@ def __init__(self, initial_stars_schema=None, **kwargs):
initial_stars_schema = afwTable.SourceTable.makeMinimalSchema()
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 @@ -374,6 +387,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 @@ -389,12 +403,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, *, exposure):
def run(self, *, exposure, 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 @@ -403,6 +417,8 @@ def run(self, *, exposure):
exposure : `lsst.afw.image.Exposure`
Post-ISR exposure, with an initial WCS, VisitInfo, and Filter.
Modified in-place during processing.
id_generator : `lsst.meas.base.IdGenerator`, optional
Object that generates source IDs and provides random seeds.
Returns
-------
Expand Down Expand Up @@ -432,11 +448,14 @@ def run(self, *, exposure):
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()

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 @@ -546,7 +565,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 @@ -557,6 +576,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 @@ -569,6 +590,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)
Expand Down
47 changes: 33 additions & 14 deletions tests/test_calibrateImage.py
View file
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 All @@ -120,7 +125,7 @@ def test_run(self):
calibrate = CalibrateImageTask(config=self.config)
calibrate.astrometry.setRefObjLoader(self.ref_loader)
calibrate.photometry.match.setRefObjLoader(self.ref_loader)
result = calibrate.run(exposure=self.exposure)
result = calibrate.run(exposure=self.exposure, id_generator=self.id_generator)

# Background should have 4 elements: 3 from compute_psf and one from
# re-estimation during source detection.
Expand Down Expand Up @@ -196,14 +201,15 @@ def test_find_stars(self):
sources, background, candidates = calibrate._compute_psf(self.exposure)
calibrate._measure_aperture_correction(self.exposure, sources)

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), 4)
# Only psf-like sources with S/N>10 should be in the output catalog,
# plus two sky sources
self.assertEqual(len(stars), 6)
self.assertTrue(sources.isContiguous())
# Sort in order of brightness, to easily compare with expected positions.
sources.sort(sources.getPsfFluxSlot().getMeasKey())
Expand All @@ -219,12 +225,13 @@ def test_astrometry(self):
calibrate.astrometry.setRefObjLoader(self.ref_loader)
sources, background, candidates = calibrate._compute_psf(self.exposure)
calibrate._measure_aperture_correction(self.exposure, sources)
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), 30.0)
Expand All @@ -238,7 +245,7 @@ def test_photometry(self):
calibrate.photometry.match.setRefObjLoader(self.ref_loader)
sources, background, candidates = calibrate._compute_psf(self.exposure)
calibrate._measure_aperture_correction(self.exposure, sources)
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 @@ -252,15 +259,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)


class CalibrateImageTaskRunQuantumTests(lsst.utils.tests.TestCase):
Expand All @@ -277,8 +290,14 @@ def setUp(self):
self.repo_path = tempfile.TemporaryDirectory()
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", exposure)
butlerTests.addDataIdValue(self.repo, "visit", visit)
butlerTests.addDataIdValue(self.repo, "detector", detector)
Expand Down Expand Up @@ -354,7 +373,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(), {"exposure"})
self.assertEqual(mock_run.call_args.kwargs.keys(), {"exposure", "id_generator"})

def test_runQuantum_no_optional_outputs(self):
config = CalibrateImageTask.ConfigClass()
Expand All @@ -379,7 +398,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(), {"exposure"})
self.assertEqual(mock_run.call_args.kwargs.keys(), {"exposure", "id_generator"})

def test_lintConnections(self):
"""Check that the connections are self-consistent.
Expand Down

0 comments on commit fa72e1a

Please sign in to comment.