Coverage for python/lsst/meas/astrom/astrometry.py: 18%
162 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-14 02:05 -0700
1# This file is part of meas_astrom.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22__all__ = ["AstrometryConfig", "AstrometryTask"]
24import numpy as np
25from astropy import units
26import scipy.stats
28import lsst.pex.config as pexConfig
29import lsst.pipe.base as pipeBase
30from lsst.utils.timer import timeMethod
31from . import exceptions
32from .ref_match import RefMatchTask, RefMatchConfig
33from .fitTanSipWcs import FitTanSipWcsTask
34from .display import displayAstrometry
37class AstrometryConfig(RefMatchConfig):
38 """Config for AstrometryTask.
39 """
40 wcsFitter = pexConfig.ConfigurableField(
41 target=FitTanSipWcsTask,
42 doc="WCS fitter",
43 )
44 forceKnownWcs = pexConfig.Field(
45 dtype=bool,
46 doc="If True then load reference objects and match sources but do not fit a WCS; "
47 "this simply controls whether 'run' calls 'solve' or 'loadAndMatch'",
48 default=False,
49 )
50 maxIter = pexConfig.RangeField(
51 doc="maximum number of iterations of match sources and fit WCS"
52 "ignored if not fitting a WCS",
53 dtype=int,
54 default=3,
55 min=1,
56 )
57 minMatchDistanceArcSec = pexConfig.RangeField(
58 doc="the match distance below which further iteration is pointless (arcsec); "
59 "ignored if not fitting a WCS",
60 dtype=float,
61 default=0.001,
62 min=0,
63 )
64 maxMeanDistanceArcsec = pexConfig.RangeField(
65 doc="Maximum mean on-sky distance (in arcsec) between matched source and reference "
66 "objects post-fit. A mean distance greater than this threshold raises BadAstrometryFit "
67 "and the WCS fit is considered a failure. The default is set to the maximum tolerated "
68 "by the external global calibration (e.g. jointcal) step for conceivable recovery; "
69 "the appropriate value will be dataset and workflow dependent.",
70 dtype=float,
71 default=0.5,
72 min=0,
73 )
74 doMagnitudeOutlierRejection = pexConfig.Field(
75 dtype=bool,
76 doc=("If True then a rough zeropoint will be computed from matched sources "
77 "and outliers will be rejected in the iterations."),
78 default=True,
79 )
80 magnitudeOutlierRejectionNSigma = pexConfig.Field(
81 dtype=float,
82 doc=("Number of sigma (measured from the distribution) in magnitude "
83 "for a potential reference/source match to be rejected during "
84 "iteration."),
85 default=3.0,
86 )
87 fiducialZeroPoint = pexConfig.DictField(
88 keytype=str,
89 itemtype=float,
90 doc="Fiducial zeropoint values, keyed by band.",
91 default=None,
92 optional=True,
93 )
94 doFiducialZeroPointCull = pexConfig.Field(
95 dtype=bool,
96 doc="If True, use the obs_package-defined fiducial zeropoint values to compute the "
97 "expected zeropoint for the current exposure. This is for use in culling reference "
98 "objects down to the approximate magnitude range of the detection source catalog "
99 "used for astrometric calibration.",
100 default=False,
101 )
102 cullMagBuffer = pexConfig.Field(
103 dtype=float,
104 doc="Generous buffer on the fiducial zero point culling to account for observing "
105 "condition variations and uncertainty of the fiducial values.",
106 default=0.3,
107 optional=True,
108 )
110 def setDefaults(self):
111 super().setDefaults()
112 # Astrometry needs sources to be isolated, too.
113 self.sourceSelector["science"].doRequirePrimary = True
114 self.sourceSelector["science"].doIsolated = True
115 self.referenceSelector.doCullFromMaskedRegion = True
117 def validate(self):
118 super().validate()
119 if self.doFiducialZeroPointCull and self.fiducialZeroPoint is None:
120 msg = "doFiducialZeroPointCull=True requires `fiducialZeroPoint`, a dict of the "
121 "fiducial zeropoints measured for the camera/filter, be set."
122 raise pexConfig.FieldValidationError(AstrometryConfig.fiducialZeroPoint, self, msg)
125class AstrometryTask(RefMatchTask):
126 """Match an input source catalog with objects from a reference catalog and
127 solve for the WCS.
129 This task is broken into two main subasks: matching and WCS fitting which
130 are very interactive. The matching here can be considered in part a first
131 pass WCS fitter due to the fitter's sensitivity to outliers.
133 Parameters
134 ----------
135 refObjLoader : `lsst.meas.algorithms.ReferenceLoader`
136 A reference object loader object; gen3 pipeline tasks will pass `None`
137 and call `setRefObjLoader` in `runQuantum`.
138 schema : `lsst.afw.table.Schema`
139 Used to set "calib_astrometry_used" flag in output source catalog.
140 **kwargs
141 Additional keyword arguments for pipe_base
142 `lsst.pipe.base.Task.__init__`.
143 """
144 ConfigClass = AstrometryConfig
145 _DefaultName = "astrometricSolver"
147 def __init__(self, refObjLoader=None, schema=None, **kwargs):
148 RefMatchTask.__init__(self, refObjLoader=refObjLoader, **kwargs)
150 if schema is not None:
151 self.usedKey = schema.addField("calib_astrometry_used", type="Flag",
152 doc="set if source was used in astrometric calibration")
153 else:
154 self.usedKey = None
156 self.makeSubtask("wcsFitter")
158 @timeMethod
159 def run(self, sourceCat, exposure):
160 """Load reference objects, match sources and optionally fit a WCS.
162 This is a thin layer around solve or loadAndMatch, depending on
163 config.forceKnownWcs.
165 Parameters
166 ----------
167 exposure : `lsst.afw.image.Exposure`
168 exposure whose WCS is to be fit
169 The following are read only:
171 - bbox
172 - filter (may be unset)
173 - detector (if wcs is pure tangent; may be absent)
175 The following are updated:
177 - wcs (the initial value is used as an initial guess, and is
178 required)
180 sourceCat : `lsst.afw.table.SourceCatalog`
181 catalog of sources detected on the exposure
183 Returns
184 -------
185 result : `lsst.pipe.base.Struct`
186 with these fields:
188 - ``refCat`` : reference object catalog of objects that overlap the
189 exposure (with some margin) (`lsst.afw.table.SimpleCatalog`).
190 - ``matches`` : astrometric matches
191 (`list` of `lsst.afw.table.ReferenceMatch`).
192 - ``scatterOnSky`` : median on-sky separation between reference
193 objects and sources in "matches"
194 (`lsst.afw.geom.Angle`) or `None` if config.forceKnownWcs True
195 - ``matchMeta`` : metadata needed to unpersist matches
196 (`lsst.daf.base.PropertyList`)
197 """
198 if self.refObjLoader is None:
199 raise RuntimeError("Running matcher task with no refObjLoader set in __init__ or setRefObjLoader")
200 if self.config.forceKnownWcs:
201 res = self.loadAndMatch(exposure=exposure, sourceCat=sourceCat)
202 res.scatterOnSky = None
203 else:
204 res = self.solve(exposure=exposure, sourceCat=sourceCat)
205 return res
207 @timeMethod
208 def solve(self, exposure, sourceCat):
209 """Load reference objects overlapping an exposure, match to sources and
210 fit a WCS
212 Returns
213 -------
214 result : `lsst.pipe.base.Struct`
215 Result struct with components:
217 - ``refCat`` : reference object catalog of objects that overlap the
218 exposure (with some margin) (`lsst::afw::table::SimpleCatalog`).
219 - ``matches`` : astrometric matches
220 (`list` of `lsst.afw.table.ReferenceMatch`).
221 - ``scatterOnSky`` : median on-sky separation between reference
222 objects and sources in "matches" (`lsst.geom.Angle`)
223 - ``matchMeta`` : metadata needed to unpersist matches
224 (`lsst.daf.base.PropertyList`)
226 Raises
227 ------
228 BadAstrometryFit
229 If the measured mean on-sky distance between the matched source and
230 reference objects is greater than
231 ``self.config.maxMeanDistanceArcsec``.
233 Notes
234 -----
235 ignores config.forceKnownWcs
236 """
237 if self.refObjLoader is None:
238 raise RuntimeError("Running matcher task with no refObjLoader set in __init__ or setRefObjLoader")
239 import lsstDebug
240 debug = lsstDebug.Info(__name__)
242 epoch = exposure.visitInfo.date.toAstropy()
244 sourceSelection = self.sourceSelector.run(sourceCat)
246 self.log.info("Purged %d sources, leaving %d good sources",
247 len(sourceCat) - len(sourceSelection.sourceCat),
248 len(sourceSelection.sourceCat))
249 if len(sourceSelection.sourceCat) == 0:
250 raise exceptions.AstrometryError(
251 "No good sources selected for astrometry.",
252 lenSourceSelectionCat=len(sourceSelection.sourceCat)
253 )
255 loadResult = self.refObjLoader.loadPixelBox(
256 bbox=exposure.getBBox(),
257 wcs=exposure.wcs,
258 filterName=exposure.filter.bandLabel,
259 epoch=epoch,
260 )
261 refSelection = self.referenceSelector.run(loadResult.refCat, exposure=exposure)
262 refCat = refSelection.sourceCat
264 if self.config.doFiducialZeroPointCull:
265 nRefCatPreCull = len(refCat)
266 # Compute rough limiting magnitudes of selected sources
267 psfFlux = sourceSelection.sourceCat["base_PsfFlux_instFlux"]
268 expTime = exposure.visitInfo.getExposureTime()
269 fiducialZeroPoint = self.config.fiducialZeroPoint[exposure.filter.bandLabel]
270 psfMag = -2.5*np.log10(psfFlux) + fiducialZeroPoint + 2.5*np.log10(expTime)
271 sourceMagMin = min(psfMag) - self.config.cullMagBuffer
272 sourceMagMax = max(psfMag) + self.config.cullMagBuffer
273 self.log.info("sourceMagMin = %0.3f sourceMagMax = %0.3f", sourceMagMin, sourceMagMax)
275 refCat = refCat[np.isfinite(refCat[loadResult.fluxField])]
276 refMag = (refCat[loadResult.fluxField]*units.nJy).to_value(units.ABmag)
277 refCat = refCat[(refMag > sourceMagMin) & (refMag < sourceMagMax)]
278 self.log.info("Selected %d/%d reference sources based on fiducial zeropoint culling.",
279 len(refCat), nRefCatPreCull)
281 # Some operations below require catalog contiguity, which is not
282 # guaranteed from the source selector.
283 if not refCat.isContiguous():
284 refCat = refCat.copy(deep=True)
286 if debug.display:
287 frame = int(debug.frame)
288 displayAstrometry(
289 refCat=refCat,
290 sourceCat=sourceSelection.sourceCat,
291 exposure=exposure,
292 bbox=exposure.getBBox(),
293 frame=frame,
294 title="Reference catalog",
295 )
297 result = pipeBase.Struct(matchTolerance=None)
298 maxMatchDistance = np.inf
299 i = 0
300 while (maxMatchDistance > self.config.minMatchDistanceArcSec and i < self.config.maxIter):
301 i += 1
302 try:
303 result = self._matchAndFitWcs(
304 refCat=refCat,
305 sourceCat=sourceCat,
306 goodSourceCat=sourceSelection.sourceCat,
307 refFluxField=loadResult.fluxField,
308 bbox=exposure.getBBox(),
309 wcs=exposure.wcs,
310 exposure=exposure,
311 matchTolerance=result.matchTolerance,
312 )
313 exposure.setWcs(result.wcs)
314 except exceptions.AstrometryError as e:
315 e._metadata['iterations'] = i
316 sourceCat["coord_ra"] = np.nan
317 sourceCat["coord_dec"] = np.nan
318 exposure.setWcs(None)
319 self.log.error("Failure fitting astrometry. %s: %s", type(e).__name__, e)
320 raise
322 result.stats = self._computeMatchStatsOnSky(result.matches)
323 maxMatchDistance = result.stats.maxMatchDist.asArcseconds()
324 distMean = result.stats.distMean.asArcseconds()
325 distStdDev = result.stats.distStdDev.asArcseconds()
326 self.log.info("Astrometric fit iteration %d: found %d matches with mean separation "
327 "= %0.3f +- %0.3f arcsec; max match distance = %0.3f arcsec.",
328 i, len(result.matches), distMean, distStdDev, maxMatchDistance)
330 # If fitter converged, record the scatter in the exposure metadata
331 # even if the fit was deemed a failure according to the value of
332 # the maxMeanDistanceArcsec config.
333 md = exposure.getMetadata()
334 md['SFM_ASTROM_OFFSET_MEAN'] = distMean
335 md['SFM_ASTROM_OFFSET_STD'] = distStdDev
336 md['SFM_ASTROM_OFFSET_MEDIAN'] = result.scatterOnSky.asArcseconds()
338 if self.usedKey:
339 for m in result.matches:
340 m.second.set(self.usedKey, True)
342 matchMeta = self.refObjLoader.getMetadataBox(
343 bbox=exposure.getBBox(),
344 wcs=exposure.wcs,
345 filterName=exposure.filter.bandLabel,
346 epoch=epoch,
347 )
349 return pipeBase.Struct(
350 refCat=refCat,
351 matches=result.matches,
352 scatterOnSky=result.scatterOnSky,
353 matchMeta=matchMeta,
354 )
356 def check(self, exposure, sourceCat, nMatches):
357 """Validate the astrometric fit against the maxMeanDistance threshold.
359 If the distMean metric does not satisfy the requirement of being less
360 than the value set in config.maxMeanDistanceArcsec, the WCS on the
361 exposure will be set to None and the coordinate values in the
362 source catalog will be set to NaN to reflect a failed astrometric fit.
364 Parameters
365 ----------
366 exposure : `lsst.afw.image.Exposure`
367 The exposure whose astrometric fit is being evaluated.
368 sourceCat : `lsst.afw.table.SourceCatalog`
369 The catalog of sources associated with the exposure.
370 nMatches : `int`
371 The number of matches that were found and used during
372 the astrometric fit (for logging purposes only).
374 Raises
375 ------
376 BadAstrometryFit
377 If the measured mean on-sky distance between the matched source and
378 reference objects is greater than
379 ``self.config.maxMeanDistanceArcsec``.
380 """
381 # Poor quality fits are a failure.
382 md = exposure.getMetadata()
383 distMean = md['SFM_ASTROM_OFFSET_MEAN']
384 distMedian = md['SFM_ASTROM_OFFSET_MEDIAN']
385 if distMean > self.config.maxMeanDistanceArcsec:
386 exception = exceptions.BadAstrometryFit(nMatches=nMatches, distMean=distMean,
387 maxMeanDist=self.config.maxMeanDistanceArcsec,
388 distMedian=distMedian)
389 exposure.setWcs(None)
390 sourceCat["coord_ra"] = np.nan
391 sourceCat["coord_dec"] = np.nan
392 self.log.error(exception)
393 raise exception
394 return
396 @timeMethod
397 def _matchAndFitWcs(self, refCat, sourceCat, goodSourceCat, refFluxField, bbox, wcs, matchTolerance,
398 exposure=None):
399 """Match sources to reference objects and fit a WCS.
401 Parameters
402 ----------
403 refCat : `lsst.afw.table.SimpleCatalog`
404 catalog of reference objects
405 sourceCat : `lsst.afw.table.SourceCatalog`
406 catalog of sources detected on the exposure
407 goodSourceCat : `lsst.afw.table.SourceCatalog`
408 catalog of down-selected good sources detected on the exposure
409 refFluxField : 'str'
410 field of refCat to use for flux
411 bbox : `lsst.geom.Box2I`
412 bounding box of exposure
413 wcs : `lsst.afw.geom.SkyWcs`
414 initial guess for WCS of exposure
415 matchTolerance : `lsst.meas.astrom.MatchTolerance`
416 a MatchTolerance object (or None) specifying
417 internal tolerances to the matcher. See the MatchTolerance
418 definition in the respective matcher for the class definition.
419 exposure : `lsst.afw.image.Exposure`, optional
420 exposure whose WCS is to be fit, or None; used only for the debug
421 display.
423 Returns
424 -------
425 result : `lsst.pipe.base.Struct`
426 Result struct with components:
428 - ``matches``: astrometric matches
429 (`list` of `lsst.afw.table.ReferenceMatch`).
430 - ``wcs``: the fit WCS (lsst.afw.geom.SkyWcs).
431 - ``scatterOnSky`` : median on-sky separation between reference
432 objects and sources in "matches" (`lsst.afw.geom.Angle`).
433 """
434 import lsstDebug
435 debug = lsstDebug.Info(__name__)
437 sourceFluxField = "slot_%sFlux_instFlux" % (self.config.sourceFluxType)
439 matchRes = self.matcher.matchObjectsToSources(
440 refCat=refCat,
441 sourceCat=goodSourceCat,
442 wcs=wcs,
443 sourceFluxField=sourceFluxField,
444 refFluxField=refFluxField,
445 matchTolerance=matchTolerance,
446 bbox=bbox,
447 )
448 self.log.debug("Found %s matches", len(matchRes.matches))
449 if debug.display:
450 frame = int(debug.frame)
451 displayAstrometry(
452 refCat=refCat,
453 sourceCat=matchRes.usableSourceCat,
454 matches=matchRes.matches,
455 exposure=exposure,
456 bbox=bbox,
457 frame=frame + 1,
458 title="Initial WCS",
459 )
461 if self.config.doMagnitudeOutlierRejection:
462 matches = self._removeMagnitudeOutliers(sourceFluxField, refFluxField, matchRes.matches)
463 else:
464 matches = matchRes.matches
466 self.log.debug("Fitting WCS")
467 fitRes = self.wcsFitter.fitWcs(
468 matches=matches,
469 initWcs=wcs,
470 bbox=bbox,
471 refCat=refCat,
472 sourceCat=sourceCat,
473 exposure=exposure,
474 )
475 fitWcs = fitRes.wcs
476 scatterOnSky = fitRes.scatterOnSky
477 if debug.display:
478 frame = int(debug.frame)
479 displayAstrometry(
480 refCat=refCat,
481 sourceCat=matchRes.usableSourceCat,
482 matches=matches,
483 exposure=exposure,
484 bbox=bbox,
485 frame=frame + 2,
486 title=f"Fitter: {self.wcsFitter._DefaultName}",
487 )
489 return pipeBase.Struct(
490 matches=matches,
491 wcs=fitWcs,
492 scatterOnSky=scatterOnSky,
493 matchTolerance=matchRes.matchTolerance,
494 )
496 def _removeMagnitudeOutliers(self, sourceFluxField, refFluxField, matchesIn):
497 """Remove magnitude outliers, computing a simple zeropoint.
499 Parameters
500 ----------
501 sourceFluxField : `str`
502 Field in source catalog for instrumental fluxes.
503 refFluxField : `str`
504 Field in reference catalog for fluxes (nJy).
505 matchesIn : `list` [`lsst.afw.table.ReferenceMatch`]
506 List of source/reference matches input
508 Returns
509 -------
510 matchesOut : `list` [`lsst.afw.table.ReferenceMatch`]
511 List of source/reference matches with magnitude
512 outliers removed.
513 """
514 nMatch = len(matchesIn)
515 sourceMag = np.zeros(nMatch)
516 refMag = np.zeros(nMatch)
517 for i, match in enumerate(matchesIn):
518 sourceMag[i] = -2.5*np.log10(match[1][sourceFluxField])
519 refMag[i] = (match[0][refFluxField]*units.nJy).to_value(units.ABmag)
521 deltaMag = refMag - sourceMag
522 # Protect against negative fluxes and nans in the reference catalog.
523 goodDelta, = np.where(np.isfinite(deltaMag))
524 zp = np.median(deltaMag[goodDelta])
525 # Use median absolute deviation (MAD) for zpSigma.
526 # Also require a minimum scatter to prevent floating-point errors from
527 # rejecting objects in zero-noise tests.
528 zpSigma = np.clip(scipy.stats.median_abs_deviation(deltaMag[goodDelta], scale='normal'),
529 1e-3,
530 None)
532 self.log.info("Rough zeropoint from astrometry matches is %.4f +/- %.4f.",
533 zp, zpSigma)
535 goodStars = goodDelta[(np.abs(deltaMag[goodDelta] - zp)
536 <= self.config.magnitudeOutlierRejectionNSigma*zpSigma)]
538 nOutlier = nMatch - goodStars.size
539 self.log.info("Removed %d magnitude outliers out of %d total astrometry matches.",
540 nOutlier, nMatch)
542 matchesOut = [matchesIn[idx] for idx in goodStars]
544 return matchesOut