Coverage for python / lsst / meas / astrom / astrometry.py: 16%
221 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-26 09:11 +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, exposure_region=None):
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 The exposure whose WCS is to be fit.
203 The following are read only:
204 - bbox
205 - filter (may be unset)
206 - detector (if wcs is pure tangent; may be absent)
208 The following are updated:
209 - wcs (the initial value is used as an initial guess, and is
210 required)
211 sourceCat : `lsst.afw.table.SourceCatalog`
212 The catalog of sources detected on the exposure.
213 exposure_region : `lsst.sphgeom.Region`, optional
214 The exposure region to use for the for the reference catalog
215 filtering. If `None`, this region will be set as a padded bbox +
216 current WCS of the exposure.
219 Returns
220 -------
221 result : `lsst.pipe.base.Struct`
222 with these fields:
224 - ``refCat`` : reference object catalog of objects that overlap the
225 exposure (with some margin) (`lsst.afw.table.SimpleCatalog`).
226 - ``matches`` : astrometric matches
227 (`list` of `lsst.afw.table.ReferenceMatch`).
228 - ``scatterOnSky`` : median on-sky separation between reference
229 objects and sources in "matches"
230 (`lsst.afw.geom.Angle`) or `None` if config.forceKnownWcs True
231 - ``matchMeta`` : metadata needed to unpersist matches
232 (`lsst.daf.base.PropertyList`)
233 - ``sip`` : a nested struct returned by
234 `FitSipApproximationTask.run`.
235 """
236 if self.refObjLoader is None:
237 raise RuntimeError("Running matcher task with no refObjLoader set in __init__ or setRefObjLoader")
238 if self.config.forceKnownWcs:
239 res = self.loadAndMatch(exposure=exposure, sourceCat=sourceCat)
240 res.scatterOnSky = None
241 else:
242 res = self.solve(exposure=exposure, sourceCat=sourceCat, exposure_region=exposure_region)
243 if self.config.doFitSipApproximation:
244 res.sip = self.fitSipApproximation.run(wcs=exposure.getWcs(), bbox=exposure.getBBox())
245 exposure.setWcs(res.sip.wcs)
246 return res
248 @timeMethod
249 def solve(self, exposure, sourceCat, exposure_region=None):
250 """Load reference objects overlapping an exposure, match to sources and
251 fit a WCS
253 Parameters
254 ----------
255 exposure : `lsst.afw.image.Exposure`
256 The exposure whose WCS is to be fit
257 The following are read only:
258 - bbox
259 - filter (may be unset)
260 - detector (if wcs is pure tangent; may be absent)
262 The following are updated:
263 - wcs (the initial value is used as an initial guess, and is
264 required)
265 sourceCat : `lsst.afw.table.SourceCatalog`
266 The catalog of sources detected on the exposure.
267 exposure_region : `lsst.sphgeom.Region`, optional
268 The exposure region to use for the for the reference catalog
269 filtering. If `None`, this region will be set as a padded bbox +
270 current WCS of the exposure.
272 Returns
273 -------
274 result : `lsst.pipe.base.Struct`
275 Result struct with components:
277 - ``refCat`` : reference object catalog of objects that overlap the
278 exposure (with some margin) (`lsst::afw::table::SimpleCatalog`).
279 - ``matches`` : astrometric matches
280 (`list` of `lsst.afw.table.ReferenceMatch`).
281 - ``scatterOnSky`` : median on-sky separation between reference
282 objects and sources in "matches" (`lsst.geom.Angle`)
283 - ``matchMeta`` : metadata needed to unpersist matches
284 (`lsst.daf.base.PropertyList`)
286 Raises
287 ------
288 BadAstrometryFit
289 If the measured mean on-sky distance between the matched source and
290 reference objects is greater than
291 ``self.config.maxMeanDistanceArcsec``.
293 Notes
294 -----
295 ignores config.forceKnownWcs
296 """
297 if self.refObjLoader is None:
298 raise RuntimeError("Running matcher task with no refObjLoader set in __init__ or setRefObjLoader")
299 import lsstDebug
300 debug = lsstDebug.Info(__name__)
302 epoch = exposure.visitInfo.date.toAstropy()
303 band = exposure.filter.bandLabel
305 sourceSelection = self.sourceSelector.run(sourceCat)
307 self.log.info("Purged %d sources, leaving %d good sources",
308 len(sourceCat) - len(sourceSelection.sourceCat),
309 len(sourceSelection.sourceCat))
310 if len(sourceSelection.sourceCat) == 0:
311 raise exceptions.AstrometryError(
312 "No good sources selected for astrometry.",
313 lenSourceSelectionCat=len(sourceSelection.sourceCat)
314 )
316 if exposure_region is not None:
317 loadResult = self.refObjLoader.loadRegion(
318 region=exposure_region,
319 filterName=band,
320 epoch=epoch,
321 wcsForCentroids=exposure.wcs,
322 )
323 if self.refObjLoader.config.pixelMargin > 0:
324 self.log.warning("Note that the astrometry_ref_loader.pixelMargin (currently "
325 "set to %d) is ignored when loading the reference catalog "
326 "with an exposure_region (i.e. the region is used as is, with "
327 "no additional padding).", self.refObjLoader.config.pixelMargin)
328 else:
329 loadResult = self.refObjLoader.loadPixelBox(
330 bbox=exposure.getBBox(),
331 wcs=exposure.wcs,
332 filterName=band,
333 epoch=epoch,
334 )
336 refSelection = self.referenceSelector.run(loadResult.refCat, exposure=exposure)
337 refCat = refSelection.sourceCat
339 if self.config.doFiducialZeroPointCull:
340 refCat, sourceSelection.sourceCat = self._do_fiducial_zeropoint_culling(
341 band,
342 loadResult.fluxField,
343 refCat, sourceSelection.sourceCat,
344 exposure.visitInfo.getExposureTime()
345 )
347 # Some operations below require catalog contiguity, which is not
348 # guaranteed from the source selector.
349 if not refCat.isContiguous():
350 refCat = refCat.copy(deep=True)
352 if debug.display:
353 frame = int(debug.frame)
354 displayAstrometry(
355 refCat=refCat,
356 sourceCat=sourceSelection.sourceCat,
357 exposure=exposure,
358 bbox=exposure.getBBox(),
359 frame=frame,
360 title="Reference catalog",
361 )
363 result = pipeBase.Struct(matchTolerance=None)
364 maxMatchDistance = np.inf
365 i = 0
366 while (maxMatchDistance > self.config.minMatchDistanceArcSec and i < self.config.maxIter):
367 i += 1
368 try:
369 result = self._matchAndFitWcs(
370 refCat=refCat,
371 sourceCat=sourceCat,
372 goodSourceCat=sourceSelection.sourceCat,
373 refFluxField=loadResult.fluxField,
374 bbox=exposure.getBBox(),
375 wcs=exposure.wcs,
376 exposure=exposure,
377 matchTolerance=result.matchTolerance,
378 )
379 exposure.setWcs(result.wcs)
380 except exceptions.AstrometryError as e:
381 e._metadata['iterations'] = i
382 sourceCat["coord_ra"] = np.nan
383 sourceCat["coord_dec"] = np.nan
384 exposure.setWcs(None)
385 self.log.error("Failure fitting astrometry. %s: %s", type(e).__name__, e)
386 raise
388 result.stats = self._computeMatchStatsOnSky(result.matches)
389 maxMatchDistance = result.stats.maxMatchDist.asArcseconds()
390 distMean = result.stats.distMean.asArcseconds()
391 distStdDev = result.stats.distStdDev.asArcseconds()
392 self.log.info("Astrometric fit iteration %d: found %d matches with mean separation "
393 "= %0.3f +- %0.3f arcsec; max match distance = %0.3f arcsec.",
394 i, len(result.matches), distMean, distStdDev, maxMatchDistance)
396 # If fitter converged, record the scatter in the exposure metadata
397 # even if the fit was deemed a failure according to the value of
398 # the maxMeanDistanceArcsec config.
399 md = exposure.getMetadata()
400 md['SFM_ASTROM_OFFSET_MEAN'] = distMean
401 md['SFM_ASTROM_OFFSET_STD'] = distStdDev
402 md['SFM_ASTROM_OFFSET_MEDIAN'] = result.scatterOnSky.asArcseconds()
404 if self.usedKey:
405 for m in result.matches:
406 m.second.set(self.usedKey, True)
408 matchMeta = self.refObjLoader.getMetadataBox(
409 bbox=exposure.getBBox(),
410 wcs=exposure.wcs,
411 filterName=band,
412 epoch=epoch,
413 )
415 return pipeBase.Struct(
416 refCat=refCat,
417 matches=result.matches,
418 scatterOnSky=result.scatterOnSky,
419 matchMeta=matchMeta,
420 )
422 def check(self, exposure, sourceCat, nMatches):
423 """Validate the astrometric fit against the maxMeanDistance threshold.
425 If the distMean metric does not satisfy the requirement of being less
426 than the value set in config.maxMeanDistanceArcsec, the WCS on the
427 exposure will be set to None and the coordinate values in the
428 source catalog will be set to NaN to reflect a failed astrometric fit.
430 Parameters
431 ----------
432 exposure : `lsst.afw.image.Exposure`
433 The exposure whose astrometric fit is being evaluated.
434 sourceCat : `lsst.afw.table.SourceCatalog`
435 The catalog of sources associated with the exposure.
436 nMatches : `int`
437 The number of matches that were found and used during
438 the astrometric fit (for logging purposes only).
440 Raises
441 ------
442 BadAstrometryFit
443 If the measured mean on-sky distance between the matched source and
444 reference objects is greater than
445 ``self.config.maxMeanDistanceArcsec``.
446 """
447 # Poor quality fits are a failure.
448 md = exposure.getMetadata()
449 distMean = md['SFM_ASTROM_OFFSET_MEAN']
450 distMedian = md['SFM_ASTROM_OFFSET_MEDIAN']
451 if distMean > self.config.maxMeanDistanceArcsec:
452 exception = exceptions.BadAstrometryFit(nMatches=nMatches, distMean=distMean,
453 maxMeanDist=self.config.maxMeanDistanceArcsec,
454 distMedian=distMedian)
455 exposure.setWcs(None)
456 sourceCat["coord_ra"] = np.nan
457 sourceCat["coord_dec"] = np.nan
458 self.log.error(exception)
459 raise exception
460 return
462 @timeMethod
463 def _matchAndFitWcs(self, refCat, sourceCat, goodSourceCat, refFluxField, bbox, wcs, matchTolerance,
464 exposure=None):
465 """Match sources to reference objects and fit a WCS.
467 Parameters
468 ----------
469 refCat : `lsst.afw.table.SimpleCatalog`
470 catalog of reference objects
471 sourceCat : `lsst.afw.table.SourceCatalog`
472 catalog of sources detected on the exposure
473 goodSourceCat : `lsst.afw.table.SourceCatalog`
474 catalog of down-selected good sources detected on the exposure
475 refFluxField : 'str'
476 field of refCat to use for flux
477 bbox : `lsst.geom.Box2I`
478 bounding box of exposure
479 wcs : `lsst.afw.geom.SkyWcs`
480 initial guess for WCS of exposure
481 matchTolerance : `lsst.meas.astrom.MatchTolerance`
482 a MatchTolerance object (or None) specifying
483 internal tolerances to the matcher. See the MatchTolerance
484 definition in the respective matcher for the class definition.
485 exposure : `lsst.afw.image.Exposure`, optional
486 exposure whose WCS is to be fit, or None; used only for the debug
487 display.
489 Returns
490 -------
491 result : `lsst.pipe.base.Struct`
492 Result struct with components:
494 - ``matches``: astrometric matches
495 (`list` of `lsst.afw.table.ReferenceMatch`).
496 - ``wcs``: the fit WCS (lsst.afw.geom.SkyWcs).
497 - ``scatterOnSky`` : median on-sky separation between reference
498 objects and sources in "matches" (`lsst.afw.geom.Angle`).
499 """
500 import lsstDebug
501 debug = lsstDebug.Info(__name__)
503 sourceFluxField = "slot_%sFlux_instFlux" % (self.config.sourceFluxType)
505 matchRes = self.matcher.matchObjectsToSources(
506 refCat=refCat,
507 sourceCat=goodSourceCat,
508 wcs=wcs,
509 sourceFluxField=sourceFluxField,
510 refFluxField=refFluxField,
511 matchTolerance=matchTolerance,
512 bbox=bbox,
513 )
514 self.log.debug("Found %s matches", len(matchRes.matches))
515 if debug.display:
516 frame = int(debug.frame)
517 displayAstrometry(
518 refCat=refCat,
519 sourceCat=matchRes.usableSourceCat,
520 matches=matchRes.matches,
521 exposure=exposure,
522 bbox=bbox,
523 frame=frame + 1,
524 title="Initial WCS",
525 )
527 if self.config.doMagnitudeOutlierRejection:
528 matches = self._removeMagnitudeOutliers(sourceFluxField, refFluxField, matchRes.matches)
529 else:
530 matches = matchRes.matches
532 self.log.debug("Fitting WCS")
533 fitRes = self.wcsFitter.fitWcs(
534 matches=matches,
535 initWcs=wcs,
536 bbox=bbox,
537 refCat=refCat,
538 sourceCat=sourceCat,
539 exposure=exposure,
540 )
541 fitWcs = fitRes.wcs
542 scatterOnSky = fitRes.scatterOnSky
543 if debug.display:
544 frame = int(debug.frame)
545 displayAstrometry(
546 refCat=refCat,
547 sourceCat=matchRes.usableSourceCat,
548 matches=matches,
549 exposure=exposure,
550 bbox=bbox,
551 frame=frame + 2,
552 title=f"Fitter: {self.wcsFitter._DefaultName}",
553 )
555 return pipeBase.Struct(
556 matches=matches,
557 wcs=fitWcs,
558 scatterOnSky=scatterOnSky,
559 matchTolerance=matchRes.matchTolerance,
560 )
562 def _removeMagnitudeOutliers(self, sourceFluxField, refFluxField, matchesIn):
563 """Remove magnitude outliers, computing a simple zeropoint.
565 Parameters
566 ----------
567 sourceFluxField : `str`
568 Field in source catalog for instrumental fluxes.
569 refFluxField : `str`
570 Field in reference catalog for fluxes (nJy).
571 matchesIn : `list` [`lsst.afw.table.ReferenceMatch`]
572 List of source/reference matches input
574 Returns
575 -------
576 matchesOut : `list` [`lsst.afw.table.ReferenceMatch`]
577 List of source/reference matches with magnitude
578 outliers removed.
579 """
580 nMatch = len(matchesIn)
581 sourceMag = np.zeros(nMatch)
582 refMag = np.zeros(nMatch)
583 for i, match in enumerate(matchesIn):
584 sourceMag[i] = -2.5*np.log10(match[1][sourceFluxField])
585 refMag[i] = (match[0][refFluxField]*units.nJy).to_value(units.ABmag)
587 deltaMag = refMag - sourceMag
588 # Protect against negative fluxes and nans in the reference catalog.
589 goodDelta, = np.where(np.isfinite(deltaMag))
590 zp = np.median(deltaMag[goodDelta])
591 # Use median absolute deviation (MAD) for zpSigma.
592 # Also require a minimum scatter to prevent floating-point errors from
593 # rejecting objects in zero-noise tests.
594 zpSigma = np.clip(scipy.stats.median_abs_deviation(deltaMag[goodDelta], scale='normal'),
595 1e-3,
596 None)
598 self.log.info("Rough zeropoint from astrometry matches is %.4f +/- %.4f.",
599 zp, zpSigma)
601 goodStars = goodDelta[(np.abs(deltaMag[goodDelta] - zp)
602 <= self.config.magnitudeOutlierRejectionNSigma*zpSigma)]
604 nOutlier = nMatch - goodStars.size
605 self.log.info("Removed %d magnitude outliers out of %d total astrometry matches.",
606 nOutlier, nMatch)
608 matchesOut = [matchesIn[idx] for idx in goodStars]
610 return matchesOut
612 def _compute_ref_src_filter_diff(self, band, refFluxField, refCat):
613 """Compute the median ref flux - source filter color difference.
615 The median difference in color between the flux field used for
616 selection and that of the observations being calibrated is computed
617 from the values for each in the reference catalog.
619 Parameters
620 ----------
621 band : `str`
622 Bandpass of observed data.
623 refFluxField : `str`
624 Name of the flux field used in the reference catalog.
625 refCat : `lsst.afw.table.SimpleCatalog`
626 Catalog of reference objects from which to compute the color
627 offset.
629 Returns
630 -------
631 refColorDelta : `float`
632 The median color difference.
633 """
634 if band in self.config.filterMap:
635 refFilterNameForBand = self.config.filterMap.get(band) + "_flux"
636 refCatTemp = refCat[(np.isfinite(refCat[refFluxField])
637 & np.isfinite(refCat[refFilterNameForBand]))].copy(deep=True)
638 if len(refCatTemp) < 3:
639 refColorDeltaDefaults = self.config.refColorDeltaDefaults
640 if band in refColorDeltaDefaults:
641 refColorDelta = refColorDeltaDefaults[band]
642 else:
643 self.log.warning("Band %s was not found in refColorDeltaDefaults dict "
644 "(currently set as: %s. Will use a default of 0.0 (i.e. "
645 "no color shift).", band, self.config.refColorDeltaDefaults)
646 refColorDelta = 0.0
647 self.log.warning("Not enough reference sources with finite fluxes in %s and %s, "
648 "so can't compute color shift; a default vaulue of %.2f will "
649 "be applied.", refFluxField, refFilterNameForBand, refColorDelta)
650 else:
651 refMag = (refCatTemp[refFluxField]*units.nJy).to_value(units.ABmag)
652 refMagSrcBand = (refCatTemp[refFilterNameForBand]*units.nJy).to_value(units.ABmag)
653 refColorDelta = np.nanmedian(refMag - refMagSrcBand)
654 self.log.info("Adjusting refCat cutoffs for color shift: %s - %s = %.2f.",
655 refFluxField, refFilterNameForBand, refColorDelta)
656 else:
657 refColorDelta = 0.0
658 self.log.warning("Band %s not found in filterMap: %s. No adjustment for filter "
659 "differences between reference and source catalogs attempted.",
660 band, self.config.filterMap)
661 return refColorDelta
663 def _do_fiducial_zeropoint_culling(self, band, refFluxField, refCat, sourceCat, expTime):
664 """Perform a culling of the catalogs to attempt to match their
665 effective magnitude ranges.
667 This uses a fiducial zeropoint along with the exposure time for
668 the observations to compute our best-guess magnitudes. Also,
669 accommodation is made for the median difference in color between
670 the flux field used for selection and that of the observations being
671 calibrated, which is computed from the values for each in the reference
672 catalog.
674 Parameters
675 ----------
676 band : `str`
677 Bandpass of observed data.
678 refFluxField : `str`
679 Name of the flux field used in the reference catalog.
680 refCat : `lsst.afw.table.SimpleCatalog`
681 Catalog of reference objects to be passed to the matcher. Modified
682 in place.
683 sourceCat : `lsst.afw.table.SourceCatalog`
684 Catalog of observed sources to be passed to the matcher. Modified
685 in place.
686 expTime : `float`
687 Exposure time of the observation being calibrated.
689 Returns
690 -------
691 refColorDelta : `float`
692 The median color difference.
693 """
694 nRefCatPreCull = len(refCat)
695 nSelectedSourceCat = len(sourceCat)
696 # Compute rough limiting magnitudes of selected sources
697 psfFlux = sourceCat["base_PsfFlux_instFlux"]
698 fiducialZeroPoint = self.config.fiducialZeroPoint[band]
699 psfMag = -2.5*np.log10(psfFlux) + fiducialZeroPoint + 2.5*np.log10(expTime)
700 sourceMagMin = min(psfMag) - self.config.cullMagBuffer
701 sourceMagMax = max(psfMag) + self.config.cullMagBuffer
703 # Try to account for median ref flux - source band color difference.
704 refColorDelta = 0.0
705 refColorDelta = self._compute_ref_src_filter_diff(band, refFluxField, refCat)
706 if refColorDelta > self.config.cullMagBuffer:
707 # Start shifting by just half of the median color difference.
708 sourceMagMin += 0.5*refColorDelta
709 sourceMagMax += 0.5*refColorDelta
710 if refColorDelta > 3.0*self.config.cullMagBuffer:
711 # Shift even further if color difference is large compared
712 # with the cullMagBuffer.
713 sourceMagMin = np.nanmin((refCat[refFluxField]*units.nJy).to_value(units.ABmag))
714 sourceMagMax += 0.5*refColorDelta
715 if refColorDelta < -1.0*self.config.cullMagBuffer:
716 # If the color difference is negative (i.e. sources are fainter
717 # in the observed filter), allow full bright end, and shift to
718 # fainter limit.
719 sourceMagMin = np.nanmin((refCat[refFluxField]*units.nJy).to_value(units.ABmag))
720 sourceMagMax += refColorDelta
721 self.log.debug("Number of sources = %d; Number of refs = %d; refs/source = %.2f.",
722 nSelectedSourceCat, nRefCatPreCull, nRefCatPreCull/nSelectedSourceCat)
724 # Include full bright end of reference catalog if there are very
725 # few sources compared to references loaded (this often occurs when
726 # there is a large amount of exctinction and/or scattering from
727 # from bright stars).
728 if nRefCatPreCull/nSelectedSourceCat > self.config.maxRefToSourceRatio:
729 sourceMagMin = np.nanmin((refCat[refFluxField]*units.nJy).to_value(units.ABmag))
731 self.log.info("Source selection: sourceMag (min, max) = (%.3f, %.3f)", sourceMagMin, sourceMagMax)
732 refCat = refCat[np.isfinite(refCat[refFluxField])]
733 refMag = (refCat[refFluxField]*units.nJy).to_value(units.ABmag)
734 refCat = refCat[(refMag < sourceMagMax) & (refMag > sourceMagMin)]
735 refMagMin = np.nanmin(refMag)
736 refMagMax = np.nanmax(refMag)
737 # Now make sure source cat doesn't extend beyond refCat limits.
738 goodSources = ((psfMag < refMagMax - refColorDelta) & (psfMag > refMagMin - refColorDelta))
739 if len(goodSources) < np.sum(goodSources):
740 sourceCat = sourceCat[goodSources].copy(deep=True)
741 nSelectedSourceCat = len(sourceCat)
742 self.log.debug("Number of sources = %d; Number of refs = %d; refs/source = %.2f.",
743 nSelectedSourceCat, nRefCatPreCull, nRefCatPreCull/nSelectedSourceCat)
745 self.log.info("Final: Selected %d/%d reference sources based on fiducial zeropoint culling. "
746 "Mag range in %s = (%.2f, %.2f)", len(refCat), nRefCatPreCull,
747 refFluxField, refMagMin, refMagMax)
748 return refCat, sourceCat