Coverage for python / lsst / meas / astrom / astrometry.py: 16%
217 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-14 23:59 +0000
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
35from .fit_sip_approximation import FitSipApproximationTask
38class AstrometryConfig(RefMatchConfig):
39 """Config for AstrometryTask.
40 """
41 wcsFitter = pexConfig.ConfigurableField(
42 target=FitTanSipWcsTask,
43 doc="WCS fitter",
44 )
45 forceKnownWcs = pexConfig.Field(
46 dtype=bool,
47 doc="If True then load reference objects and match sources but do not fit a WCS; "
48 "this simply controls whether 'run' calls 'solve' or 'loadAndMatch'",
49 default=False,
50 )
51 maxIter = pexConfig.RangeField(
52 doc="maximum number of iterations of match sources and fit WCS"
53 "ignored if not fitting a WCS",
54 dtype=int,
55 default=3,
56 min=1,
57 )
58 minMatchDistanceArcSec = pexConfig.RangeField(
59 doc="the match distance below which further iteration is pointless (arcsec); "
60 "ignored if not fitting a WCS",
61 dtype=float,
62 default=0.001,
63 min=0,
64 )
65 maxMeanDistanceArcsec = pexConfig.RangeField(
66 doc="Maximum mean on-sky distance (in arcsec) between matched source and reference "
67 "objects post-fit. A mean distance greater than this threshold raises BadAstrometryFit "
68 "and the WCS fit is considered a failure. The default is set to the maximum tolerated "
69 "by the external global calibration (e.g. jointcal) step for conceivable recovery; "
70 "the appropriate value will be dataset and workflow dependent.",
71 dtype=float,
72 default=0.5,
73 min=0,
74 )
75 doMagnitudeOutlierRejection = pexConfig.Field(
76 dtype=bool,
77 doc=("If True then a rough zeropoint will be computed from matched sources "
78 "and outliers will be rejected in the iterations."),
79 default=True,
80 )
81 magnitudeOutlierRejectionNSigma = pexConfig.Field(
82 dtype=float,
83 doc=("Number of sigma (measured from the distribution) in magnitude "
84 "for a potential reference/source match to be rejected during "
85 "iteration."),
86 default=3.0,
87 )
88 fiducialZeroPoint = pexConfig.DictField(
89 keytype=str,
90 itemtype=float,
91 doc="Fiducial zeropoint values, keyed by band.",
92 default=None,
93 optional=True,
94 )
95 doFiducialZeroPointCull = pexConfig.Field(
96 dtype=bool,
97 doc="If True, use the obs_package-defined fiducial zeropoint values to compute the "
98 "expected zeropoint for the current exposure. This is for use in culling reference "
99 "objects down to the approximate magnitude range of the detection source catalog "
100 "used for astrometric calibration.",
101 default=False,
102 )
103 cullMagBuffer = pexConfig.Field(
104 dtype=float,
105 doc="Generous buffer on the fiducial zero point culling to account for observing "
106 "condition variations and uncertainty of the fiducial values.",
107 default=0.3,
108 optional=True,
109 )
110 maxRefToSourceRatio = pexConfig.Field(
111 dtype=float,
112 doc="Maximum ratio of loaded reference objects to detected sources in play. If exceded "
113 "the source catalog will be trimmed to the minimum (i.e. brightest) mag of the "
114 "reference catalog.",
115 default=20.0,
116 )
117 filterMap = pexConfig.DictField(
118 doc="Mapping from physical filter label to reference filter name.",
119 keytype=str,
120 itemtype=str,
121 default={},
122 )
123 refColorDeltaDefaults = pexConfig.DictField(
124 doc="Fallback values for color differences between the reference band and the "
125 "physical filter of the observation (note that these values apply to LSSTCam "
126 "filters and may not be appropriate for others).",
127 keytype=str,
128 itemtype=float,
129 default={"u": -1.5, "g": -0.6, "r": 0.0, "i": 0.5, "z": 0.6, "y": 0.8},
130 )
131 doFitSipApproximation = pexConfig.Field(
132 "Whether to fit a TAN-SIP approximation to the true WCS for use in FITS headers.",
133 dtype=bool,
134 default=True,
135 )
136 fitSipApproximation = pexConfig.ConfigurableField(
137 "Configuration for fitting a TAN-SIP approximation to the true WCS for use in FITS headers.",
138 target=FitSipApproximationTask,
139 )
141 def setDefaults(self):
142 super().setDefaults()
143 # Astrometry needs sources to be isolated, too.
144 self.sourceSelector["science"].doRequirePrimary = True
145 self.sourceSelector["science"].doIsolated = True
146 self.sourceSelector["science"].doCentroidErrorLimit = True
147 self.referenceSelector.doCullFromMaskedRegion = True
149 def validate(self):
150 super().validate()
151 if self.doFiducialZeroPointCull and self.fiducialZeroPoint is None:
152 msg = "doFiducialZeroPointCull=True requires `fiducialZeroPoint`, a dict of the "
153 "fiducial zeropoints measured for the camera/filter, be set."
154 raise pexConfig.FieldValidationError(AstrometryConfig.fiducialZeroPoint, self, msg)
157class AstrometryTask(RefMatchTask):
158 """Match an input source catalog with objects from a reference catalog and
159 solve for the WCS.
161 This task is broken into two main subasks: matching and WCS fitting which
162 are very interactive. The matching here can be considered in part a first
163 pass WCS fitter due to the fitter's sensitivity to outliers.
165 Parameters
166 ----------
167 refObjLoader : `lsst.meas.algorithms.ReferenceLoader`
168 A reference object loader object; gen3 pipeline tasks will pass `None`
169 and call `setRefObjLoader` in `runQuantum`.
170 schema : `lsst.afw.table.Schema`
171 Used to set "calib_astrometry_used" flag in output source catalog.
172 **kwargs
173 Additional keyword arguments for pipe_base
174 `lsst.pipe.base.Task.__init__`.
175 """
176 ConfigClass = AstrometryConfig
177 _DefaultName = "astrometricSolver"
179 def __init__(self, refObjLoader=None, schema=None, **kwargs):
180 RefMatchTask.__init__(self, refObjLoader=refObjLoader, **kwargs)
182 if schema is not None:
183 self.usedKey = schema.addField("calib_astrometry_used", type="Flag",
184 doc="set if source was used in astrometric calibration")
185 else:
186 self.usedKey = None
188 self.makeSubtask("wcsFitter")
189 if self.config.doFitSipApproximation:
190 self.makeSubtask("fitSipApproximation")
192 @timeMethod
193 def run(self, sourceCat, exposure):
194 """Load reference objects, match sources and optionally fit a WCS.
196 This is a thin layer around solve or loadAndMatch, depending on
197 config.forceKnownWcs.
199 Parameters
200 ----------
201 exposure : `lsst.afw.image.Exposure`
202 exposure whose WCS is to be fit
203 The following are read only:
205 - bbox
206 - filter (may be unset)
207 - detector (if wcs is pure tangent; may be absent)
209 The following are updated:
211 - wcs (the initial value is used as an initial guess, and is
212 required)
214 sourceCat : `lsst.afw.table.SourceCatalog`
215 catalog of sources detected on the exposure
217 Returns
218 -------
219 result : `lsst.pipe.base.Struct`
220 with these fields:
222 - ``refCat`` : reference object catalog of objects that overlap the
223 exposure (with some margin) (`lsst.afw.table.SimpleCatalog`).
224 - ``matches`` : astrometric matches
225 (`list` of `lsst.afw.table.ReferenceMatch`).
226 - ``scatterOnSky`` : median on-sky separation between reference
227 objects and sources in "matches"
228 (`lsst.afw.geom.Angle`) or `None` if config.forceKnownWcs True
229 - ``matchMeta`` : metadata needed to unpersist matches
230 (`lsst.daf.base.PropertyList`)
231 - ``sip`` : a nested struct returned by
232 `FitSipApproximationTask.run`.
233 """
234 if self.refObjLoader is None:
235 raise RuntimeError("Running matcher task with no refObjLoader set in __init__ or setRefObjLoader")
236 if self.config.forceKnownWcs:
237 res = self.loadAndMatch(exposure=exposure, sourceCat=sourceCat)
238 res.scatterOnSky = None
239 else:
240 res = self.solve(exposure=exposure, sourceCat=sourceCat)
241 if self.config.doFitSipApproximation:
242 res.sip = self.fitSipApproximation.run(wcs=exposure.getWcs(), bbox=exposure.getBBox())
243 exposure.setWcs(res.sip.wcs)
244 return res
246 @timeMethod
247 def solve(self, exposure, sourceCat):
248 """Load reference objects overlapping an exposure, match to sources and
249 fit a WCS
251 Returns
252 -------
253 result : `lsst.pipe.base.Struct`
254 Result struct with components:
256 - ``refCat`` : reference object catalog of objects that overlap the
257 exposure (with some margin) (`lsst::afw::table::SimpleCatalog`).
258 - ``matches`` : astrometric matches
259 (`list` of `lsst.afw.table.ReferenceMatch`).
260 - ``scatterOnSky`` : median on-sky separation between reference
261 objects and sources in "matches" (`lsst.geom.Angle`)
262 - ``matchMeta`` : metadata needed to unpersist matches
263 (`lsst.daf.base.PropertyList`)
265 Raises
266 ------
267 BadAstrometryFit
268 If the measured mean on-sky distance between the matched source and
269 reference objects is greater than
270 ``self.config.maxMeanDistanceArcsec``.
272 Notes
273 -----
274 ignores config.forceKnownWcs
275 """
276 if self.refObjLoader is None:
277 raise RuntimeError("Running matcher task with no refObjLoader set in __init__ or setRefObjLoader")
278 import lsstDebug
279 debug = lsstDebug.Info(__name__)
281 epoch = exposure.visitInfo.date.toAstropy()
282 band = exposure.filter.bandLabel
284 sourceSelection = self.sourceSelector.run(sourceCat)
286 self.log.info("Purged %d sources, leaving %d good sources",
287 len(sourceCat) - len(sourceSelection.sourceCat),
288 len(sourceSelection.sourceCat))
289 if len(sourceSelection.sourceCat) == 0:
290 raise exceptions.AstrometryError(
291 "No good sources selected for astrometry.",
292 lenSourceSelectionCat=len(sourceSelection.sourceCat)
293 )
295 loadResult = self.refObjLoader.loadPixelBox(
296 bbox=exposure.getBBox(),
297 wcs=exposure.wcs,
298 filterName=band,
299 epoch=epoch,
300 )
302 refSelection = self.referenceSelector.run(loadResult.refCat, exposure=exposure)
303 refCat = refSelection.sourceCat
305 if self.config.doFiducialZeroPointCull:
306 refCat, sourceSelection.sourceCat = self._do_fiducial_zeropoint_culling(
307 band,
308 loadResult.fluxField,
309 refCat, sourceSelection.sourceCat,
310 exposure.visitInfo.getExposureTime()
311 )
313 # Some operations below require catalog contiguity, which is not
314 # guaranteed from the source selector.
315 if not refCat.isContiguous():
316 refCat = refCat.copy(deep=True)
318 if debug.display:
319 frame = int(debug.frame)
320 displayAstrometry(
321 refCat=refCat,
322 sourceCat=sourceSelection.sourceCat,
323 exposure=exposure,
324 bbox=exposure.getBBox(),
325 frame=frame,
326 title="Reference catalog",
327 )
329 result = pipeBase.Struct(matchTolerance=None)
330 maxMatchDistance = np.inf
331 i = 0
332 while (maxMatchDistance > self.config.minMatchDistanceArcSec and i < self.config.maxIter):
333 i += 1
334 try:
335 result = self._matchAndFitWcs(
336 refCat=refCat,
337 sourceCat=sourceCat,
338 goodSourceCat=sourceSelection.sourceCat,
339 refFluxField=loadResult.fluxField,
340 bbox=exposure.getBBox(),
341 wcs=exposure.wcs,
342 exposure=exposure,
343 matchTolerance=result.matchTolerance,
344 )
345 exposure.setWcs(result.wcs)
346 except exceptions.AstrometryError as e:
347 e._metadata['iterations'] = i
348 sourceCat["coord_ra"] = np.nan
349 sourceCat["coord_dec"] = np.nan
350 exposure.setWcs(None)
351 self.log.error("Failure fitting astrometry. %s: %s", type(e).__name__, e)
352 raise
354 result.stats = self._computeMatchStatsOnSky(result.matches)
355 maxMatchDistance = result.stats.maxMatchDist.asArcseconds()
356 distMean = result.stats.distMean.asArcseconds()
357 distStdDev = result.stats.distStdDev.asArcseconds()
358 self.log.info("Astrometric fit iteration %d: found %d matches with mean separation "
359 "= %0.3f +- %0.3f arcsec; max match distance = %0.3f arcsec.",
360 i, len(result.matches), distMean, distStdDev, maxMatchDistance)
362 # If fitter converged, record the scatter in the exposure metadata
363 # even if the fit was deemed a failure according to the value of
364 # the maxMeanDistanceArcsec config.
365 md = exposure.getMetadata()
366 md['SFM_ASTROM_OFFSET_MEAN'] = distMean
367 md['SFM_ASTROM_OFFSET_STD'] = distStdDev
368 md['SFM_ASTROM_OFFSET_MEDIAN'] = result.scatterOnSky.asArcseconds()
370 if self.usedKey:
371 for m in result.matches:
372 m.second.set(self.usedKey, True)
374 matchMeta = self.refObjLoader.getMetadataBox(
375 bbox=exposure.getBBox(),
376 wcs=exposure.wcs,
377 filterName=band,
378 epoch=epoch,
379 )
381 return pipeBase.Struct(
382 refCat=refCat,
383 matches=result.matches,
384 scatterOnSky=result.scatterOnSky,
385 matchMeta=matchMeta,
386 )
388 def check(self, exposure, sourceCat, nMatches):
389 """Validate the astrometric fit against the maxMeanDistance threshold.
391 If the distMean metric does not satisfy the requirement of being less
392 than the value set in config.maxMeanDistanceArcsec, the WCS on the
393 exposure will be set to None and the coordinate values in the
394 source catalog will be set to NaN to reflect a failed astrometric fit.
396 Parameters
397 ----------
398 exposure : `lsst.afw.image.Exposure`
399 The exposure whose astrometric fit is being evaluated.
400 sourceCat : `lsst.afw.table.SourceCatalog`
401 The catalog of sources associated with the exposure.
402 nMatches : `int`
403 The number of matches that were found and used during
404 the astrometric fit (for logging purposes only).
406 Raises
407 ------
408 BadAstrometryFit
409 If the measured mean on-sky distance between the matched source and
410 reference objects is greater than
411 ``self.config.maxMeanDistanceArcsec``.
412 """
413 # Poor quality fits are a failure.
414 md = exposure.getMetadata()
415 distMean = md['SFM_ASTROM_OFFSET_MEAN']
416 distMedian = md['SFM_ASTROM_OFFSET_MEDIAN']
417 if distMean > self.config.maxMeanDistanceArcsec:
418 exception = exceptions.BadAstrometryFit(nMatches=nMatches, distMean=distMean,
419 maxMeanDist=self.config.maxMeanDistanceArcsec,
420 distMedian=distMedian)
421 exposure.setWcs(None)
422 sourceCat["coord_ra"] = np.nan
423 sourceCat["coord_dec"] = np.nan
424 self.log.error(exception)
425 raise exception
426 return
428 @timeMethod
429 def _matchAndFitWcs(self, refCat, sourceCat, goodSourceCat, refFluxField, bbox, wcs, matchTolerance,
430 exposure=None):
431 """Match sources to reference objects and fit a WCS.
433 Parameters
434 ----------
435 refCat : `lsst.afw.table.SimpleCatalog`
436 catalog of reference objects
437 sourceCat : `lsst.afw.table.SourceCatalog`
438 catalog of sources detected on the exposure
439 goodSourceCat : `lsst.afw.table.SourceCatalog`
440 catalog of down-selected good sources detected on the exposure
441 refFluxField : 'str'
442 field of refCat to use for flux
443 bbox : `lsst.geom.Box2I`
444 bounding box of exposure
445 wcs : `lsst.afw.geom.SkyWcs`
446 initial guess for WCS of exposure
447 matchTolerance : `lsst.meas.astrom.MatchTolerance`
448 a MatchTolerance object (or None) specifying
449 internal tolerances to the matcher. See the MatchTolerance
450 definition in the respective matcher for the class definition.
451 exposure : `lsst.afw.image.Exposure`, optional
452 exposure whose WCS is to be fit, or None; used only for the debug
453 display.
455 Returns
456 -------
457 result : `lsst.pipe.base.Struct`
458 Result struct with components:
460 - ``matches``: astrometric matches
461 (`list` of `lsst.afw.table.ReferenceMatch`).
462 - ``wcs``: the fit WCS (lsst.afw.geom.SkyWcs).
463 - ``scatterOnSky`` : median on-sky separation between reference
464 objects and sources in "matches" (`lsst.afw.geom.Angle`).
465 """
466 import lsstDebug
467 debug = lsstDebug.Info(__name__)
469 sourceFluxField = "slot_%sFlux_instFlux" % (self.config.sourceFluxType)
471 matchRes = self.matcher.matchObjectsToSources(
472 refCat=refCat,
473 sourceCat=goodSourceCat,
474 wcs=wcs,
475 sourceFluxField=sourceFluxField,
476 refFluxField=refFluxField,
477 matchTolerance=matchTolerance,
478 bbox=bbox,
479 )
480 self.log.debug("Found %s matches", len(matchRes.matches))
481 if debug.display:
482 frame = int(debug.frame)
483 displayAstrometry(
484 refCat=refCat,
485 sourceCat=matchRes.usableSourceCat,
486 matches=matchRes.matches,
487 exposure=exposure,
488 bbox=bbox,
489 frame=frame + 1,
490 title="Initial WCS",
491 )
493 if self.config.doMagnitudeOutlierRejection:
494 matches = self._removeMagnitudeOutliers(sourceFluxField, refFluxField, matchRes.matches)
495 else:
496 matches = matchRes.matches
498 self.log.debug("Fitting WCS")
499 fitRes = self.wcsFitter.fitWcs(
500 matches=matches,
501 initWcs=wcs,
502 bbox=bbox,
503 refCat=refCat,
504 sourceCat=sourceCat,
505 exposure=exposure,
506 )
507 fitWcs = fitRes.wcs
508 scatterOnSky = fitRes.scatterOnSky
509 if debug.display:
510 frame = int(debug.frame)
511 displayAstrometry(
512 refCat=refCat,
513 sourceCat=matchRes.usableSourceCat,
514 matches=matches,
515 exposure=exposure,
516 bbox=bbox,
517 frame=frame + 2,
518 title=f"Fitter: {self.wcsFitter._DefaultName}",
519 )
521 return pipeBase.Struct(
522 matches=matches,
523 wcs=fitWcs,
524 scatterOnSky=scatterOnSky,
525 matchTolerance=matchRes.matchTolerance,
526 )
528 def _removeMagnitudeOutliers(self, sourceFluxField, refFluxField, matchesIn):
529 """Remove magnitude outliers, computing a simple zeropoint.
531 Parameters
532 ----------
533 sourceFluxField : `str`
534 Field in source catalog for instrumental fluxes.
535 refFluxField : `str`
536 Field in reference catalog for fluxes (nJy).
537 matchesIn : `list` [`lsst.afw.table.ReferenceMatch`]
538 List of source/reference matches input
540 Returns
541 -------
542 matchesOut : `list` [`lsst.afw.table.ReferenceMatch`]
543 List of source/reference matches with magnitude
544 outliers removed.
545 """
546 nMatch = len(matchesIn)
547 sourceMag = np.zeros(nMatch)
548 refMag = np.zeros(nMatch)
549 for i, match in enumerate(matchesIn):
550 sourceMag[i] = -2.5*np.log10(match[1][sourceFluxField])
551 refMag[i] = (match[0][refFluxField]*units.nJy).to_value(units.ABmag)
553 deltaMag = refMag - sourceMag
554 # Protect against negative fluxes and nans in the reference catalog.
555 goodDelta, = np.where(np.isfinite(deltaMag))
556 zp = np.median(deltaMag[goodDelta])
557 # Use median absolute deviation (MAD) for zpSigma.
558 # Also require a minimum scatter to prevent floating-point errors from
559 # rejecting objects in zero-noise tests.
560 zpSigma = np.clip(scipy.stats.median_abs_deviation(deltaMag[goodDelta], scale='normal'),
561 1e-3,
562 None)
564 self.log.info("Rough zeropoint from astrometry matches is %.4f +/- %.4f.",
565 zp, zpSigma)
567 goodStars = goodDelta[(np.abs(deltaMag[goodDelta] - zp)
568 <= self.config.magnitudeOutlierRejectionNSigma*zpSigma)]
570 nOutlier = nMatch - goodStars.size
571 self.log.info("Removed %d magnitude outliers out of %d total astrometry matches.",
572 nOutlier, nMatch)
574 matchesOut = [matchesIn[idx] for idx in goodStars]
576 return matchesOut
578 def _compute_ref_src_filter_diff(self, band, refFluxField, refCat):
579 """Compute the median ref flux - source filter color difference.
581 The median difference in color between the flux field used for
582 selection and that of the observations being calibrated is computed
583 from the values for each in the reference catalog.
585 Parameters
586 ----------
587 band : `str`
588 Bandpass of observed data.
589 refFluxField : `str`
590 Name of the flux field used in the reference catalog.
591 refCat : `lsst.afw.table.SimpleCatalog`
592 Catalog of reference objects from which to compute the color
593 offset.
595 Returns
596 -------
597 refColorDelta : `float`
598 The median color difference.
599 """
600 if band in self.config.filterMap:
601 refFilterNameForBand = self.config.filterMap.get(band) + "_flux"
602 refCatTemp = refCat[(np.isfinite(refCat[refFluxField])
603 & np.isfinite(refCat[refFilterNameForBand]))].copy(deep=True)
604 if len(refCatTemp) < 3:
605 refColorDeltaDefaults = self.config.refColorDeltaDefaults
606 if band in refColorDeltaDefaults:
607 refColorDelta = refColorDeltaDefaults[band]
608 else:
609 self.log.warning("Band %s was not found in refColorDeltaDefaults dict "
610 "(currently set as: %s. Will use a default of 0.0 (i.e. "
611 "no color shift).", band, self.config.refColorDeltaDefaults)
612 refColorDelta = 0.0
613 self.log.warning("Not enough reference sources with finite fluxes in %s and %s, "
614 "so can't compute color shift; a default vaulue of %.2f will "
615 "be applied.", refFluxField, refFilterNameForBand, refColorDelta)
616 else:
617 refMag = (refCatTemp[refFluxField]*units.nJy).to_value(units.ABmag)
618 refMagSrcBand = (refCatTemp[refFilterNameForBand]*units.nJy).to_value(units.ABmag)
619 refColorDelta = np.nanmedian(refMag - refMagSrcBand)
620 self.log.info("Adjusting refCat cutoffs for color shift: %s - %s = %.2f.",
621 refFluxField, refFilterNameForBand, refColorDelta)
622 else:
623 refColorDelta = 0.0
624 self.log.warning("Band %s not found in filterMap: %s. No adjustment for filter "
625 "differences between reference and source catalogs attempted.",
626 band, self.config.filterMap)
627 return refColorDelta
629 def _do_fiducial_zeropoint_culling(self, band, refFluxField, refCat, sourceCat, expTime):
630 """Perform a culling of the catalogs to attempt to match their
631 effective magnitude ranges.
633 This uses a fiducial zeropoint along with the exposure time for
634 the observations to compute our best-guess magnitudes. Also,
635 accommodation is made for the median difference in color between
636 the flux field used for selection and that of the observations being
637 calibrated, which is computed from the values for each in the reference
638 catalog.
640 Parameters
641 ----------
642 band : `str`
643 Bandpass of observed data.
644 refFluxField : `str`
645 Name of the flux field used in the reference catalog.
646 refCat : `lsst.afw.table.SimpleCatalog`
647 Catalog of reference objects to be passed to the matcher. Modified
648 in place.
649 sourceCat : `lsst.afw.table.SourceCatalog`
650 Catalog of observed sources to be passed to the matcher. Modified
651 in place.
652 expTime : `float`
653 Exposure time of the observation being calibrated.
655 Returns
656 -------
657 refColorDelta : `float`
658 The median color difference.
659 """
660 nRefCatPreCull = len(refCat)
661 nSelectedSourceCat = len(sourceCat)
662 # Compute rough limiting magnitudes of selected sources
663 psfFlux = sourceCat["base_PsfFlux_instFlux"]
664 fiducialZeroPoint = self.config.fiducialZeroPoint[band]
665 psfMag = -2.5*np.log10(psfFlux) + fiducialZeroPoint + 2.5*np.log10(expTime)
666 sourceMagMin = min(psfMag) - self.config.cullMagBuffer
667 sourceMagMax = max(psfMag) + self.config.cullMagBuffer
669 # Try to account for median ref flux - source band color difference.
670 refColorDelta = 0.0
671 refColorDelta = self._compute_ref_src_filter_diff(band, refFluxField, refCat)
672 if refColorDelta > self.config.cullMagBuffer:
673 # Start shifting by just half of the median color difference.
674 sourceMagMin += 0.5*refColorDelta
675 sourceMagMax += 0.5*refColorDelta
676 if refColorDelta > 3.0*self.config.cullMagBuffer:
677 # Shift even further if color difference is large compared
678 # with the cullMagBuffer.
679 sourceMagMin = np.nanmin((refCat[refFluxField]*units.nJy).to_value(units.ABmag))
680 sourceMagMax += 0.5*refColorDelta
681 if refColorDelta < -1.0*self.config.cullMagBuffer:
682 # If the color difference is negative (i.e. sources are fainter
683 # in the observed filter), allow full bright end, and shift to
684 # fainter limit.
685 sourceMagMin = np.nanmin((refCat[refFluxField]*units.nJy).to_value(units.ABmag))
686 sourceMagMax += refColorDelta
687 self.log.debug("Number of sources = %d; Number of refs = %d; refs/source = %.2f.",
688 nSelectedSourceCat, nRefCatPreCull, nRefCatPreCull/nSelectedSourceCat)
690 # Include full bright end of reference catalog if there are very
691 # few sources compared to references loaded (this often occurs when
692 # there is a large amount of exctinction and/or scattering from
693 # from bright stars).
694 if nRefCatPreCull/nSelectedSourceCat > self.config.maxRefToSourceRatio:
695 sourceMagMin = np.nanmin((refCat[refFluxField]*units.nJy).to_value(units.ABmag))
697 self.log.info("Source selection: sourceMag (min, max) = (%.3f, %.3f)", sourceMagMin, sourceMagMax)
698 refCat = refCat[np.isfinite(refCat[refFluxField])]
699 refMag = (refCat[refFluxField]*units.nJy).to_value(units.ABmag)
700 refCat = refCat[(refMag < sourceMagMax) & (refMag > sourceMagMin)]
701 refMagMin = np.nanmin(refMag)
702 refMagMax = np.nanmax(refMag)
703 # Now make sure source cat doesn't extend beyond refCat limits.
704 goodSources = ((psfMag < refMagMax - refColorDelta) & (psfMag > refMagMin - refColorDelta))
705 if len(goodSources) < np.sum(goodSources):
706 sourceCat = sourceCat[goodSources].copy(deep=True)
707 nSelectedSourceCat = len(sourceCat)
708 self.log.debug("Number of sources = %d; Number of refs = %d; refs/source = %.2f.",
709 nSelectedSourceCat, nRefCatPreCull, nRefCatPreCull/nSelectedSourceCat)
711 self.log.info("Final: Selected %d/%d reference sources based on fiducial zeropoint culling. "
712 "Mag range in %s = (%.2f, %.2f)", len(refCat), nRefCatPreCull,
713 refFluxField, refMagMin, refMagMax)
714 return refCat, sourceCat