Coverage for python/lsst/meas/astrom/astrometry.py: 16%
207 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-08 01:22 -0800
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 )
109 maxRefToSourceRatio = pexConfig.Field(
110 dtype=float,
111 doc="Maximum ratio of loaded reference objects to detected sources in play. If exceded "
112 "the source catalog will be trimmed to the minimum (i.e. brightest) mag of the "
113 "reference catalog.",
114 default=20.0,
115 )
116 filterMap = pexConfig.DictField(
117 doc="Mapping from physical filter label to reference filter name.",
118 keytype=str,
119 itemtype=str,
120 default={},
121 )
122 refColorDeltaDefaults = pexConfig.DictField(
123 doc="Fallback values for color differences between the reference band and the "
124 "physical filter of the observation (note that these values apply to LSSTCam "
125 "filters and may not be appropriate for others).",
126 keytype=str,
127 itemtype=float,
128 default={"u": -1.5, "g": -0.6, "r": 0.0, "i": 0.5, "z": 0.6},
129 )
131 def setDefaults(self):
132 super().setDefaults()
133 # Astrometry needs sources to be isolated, too.
134 self.sourceSelector["science"].doRequirePrimary = True
135 self.sourceSelector["science"].doIsolated = True
136 self.sourceSelector["science"].doCentroidErrorLimit = True
137 self.referenceSelector.doCullFromMaskedRegion = True
139 def validate(self):
140 super().validate()
141 if self.doFiducialZeroPointCull and self.fiducialZeroPoint is None:
142 msg = "doFiducialZeroPointCull=True requires `fiducialZeroPoint`, a dict of the "
143 "fiducial zeropoints measured for the camera/filter, be set."
144 raise pexConfig.FieldValidationError(AstrometryConfig.fiducialZeroPoint, self, msg)
147class AstrometryTask(RefMatchTask):
148 """Match an input source catalog with objects from a reference catalog and
149 solve for the WCS.
151 This task is broken into two main subasks: matching and WCS fitting which
152 are very interactive. The matching here can be considered in part a first
153 pass WCS fitter due to the fitter's sensitivity to outliers.
155 Parameters
156 ----------
157 refObjLoader : `lsst.meas.algorithms.ReferenceLoader`
158 A reference object loader object; gen3 pipeline tasks will pass `None`
159 and call `setRefObjLoader` in `runQuantum`.
160 schema : `lsst.afw.table.Schema`
161 Used to set "calib_astrometry_used" flag in output source catalog.
162 **kwargs
163 Additional keyword arguments for pipe_base
164 `lsst.pipe.base.Task.__init__`.
165 """
166 ConfigClass = AstrometryConfig
167 _DefaultName = "astrometricSolver"
169 def __init__(self, refObjLoader=None, schema=None, **kwargs):
170 RefMatchTask.__init__(self, refObjLoader=refObjLoader, **kwargs)
172 if schema is not None:
173 self.usedKey = schema.addField("calib_astrometry_used", type="Flag",
174 doc="set if source was used in astrometric calibration")
175 else:
176 self.usedKey = None
178 self.makeSubtask("wcsFitter")
180 @timeMethod
181 def run(self, sourceCat, exposure):
182 """Load reference objects, match sources and optionally fit a WCS.
184 This is a thin layer around solve or loadAndMatch, depending on
185 config.forceKnownWcs.
187 Parameters
188 ----------
189 exposure : `lsst.afw.image.Exposure`
190 exposure whose WCS is to be fit
191 The following are read only:
193 - bbox
194 - filter (may be unset)
195 - detector (if wcs is pure tangent; may be absent)
197 The following are updated:
199 - wcs (the initial value is used as an initial guess, and is
200 required)
202 sourceCat : `lsst.afw.table.SourceCatalog`
203 catalog of sources detected on the exposure
205 Returns
206 -------
207 result : `lsst.pipe.base.Struct`
208 with these fields:
210 - ``refCat`` : reference object catalog of objects that overlap the
211 exposure (with some margin) (`lsst.afw.table.SimpleCatalog`).
212 - ``matches`` : astrometric matches
213 (`list` of `lsst.afw.table.ReferenceMatch`).
214 - ``scatterOnSky`` : median on-sky separation between reference
215 objects and sources in "matches"
216 (`lsst.afw.geom.Angle`) or `None` if config.forceKnownWcs True
217 - ``matchMeta`` : metadata needed to unpersist matches
218 (`lsst.daf.base.PropertyList`)
219 """
220 if self.refObjLoader is None:
221 raise RuntimeError("Running matcher task with no refObjLoader set in __init__ or setRefObjLoader")
222 if self.config.forceKnownWcs:
223 res = self.loadAndMatch(exposure=exposure, sourceCat=sourceCat)
224 res.scatterOnSky = None
225 else:
226 res = self.solve(exposure=exposure, sourceCat=sourceCat)
227 return res
229 @timeMethod
230 def solve(self, exposure, sourceCat):
231 """Load reference objects overlapping an exposure, match to sources and
232 fit a WCS
234 Returns
235 -------
236 result : `lsst.pipe.base.Struct`
237 Result struct with components:
239 - ``refCat`` : reference object catalog of objects that overlap the
240 exposure (with some margin) (`lsst::afw::table::SimpleCatalog`).
241 - ``matches`` : astrometric matches
242 (`list` of `lsst.afw.table.ReferenceMatch`).
243 - ``scatterOnSky`` : median on-sky separation between reference
244 objects and sources in "matches" (`lsst.geom.Angle`)
245 - ``matchMeta`` : metadata needed to unpersist matches
246 (`lsst.daf.base.PropertyList`)
248 Raises
249 ------
250 BadAstrometryFit
251 If the measured mean on-sky distance between the matched source and
252 reference objects is greater than
253 ``self.config.maxMeanDistanceArcsec``.
255 Notes
256 -----
257 ignores config.forceKnownWcs
258 """
259 if self.refObjLoader is None:
260 raise RuntimeError("Running matcher task with no refObjLoader set in __init__ or setRefObjLoader")
261 import lsstDebug
262 debug = lsstDebug.Info(__name__)
264 epoch = exposure.visitInfo.date.toAstropy()
265 band = exposure.filter.bandLabel
267 sourceSelection = self.sourceSelector.run(sourceCat)
269 self.log.info("Purged %d sources, leaving %d good sources",
270 len(sourceCat) - len(sourceSelection.sourceCat),
271 len(sourceSelection.sourceCat))
272 if len(sourceSelection.sourceCat) == 0:
273 raise exceptions.AstrometryError(
274 "No good sources selected for astrometry.",
275 lenSourceSelectionCat=len(sourceSelection.sourceCat)
276 )
278 loadResult = self.refObjLoader.loadPixelBox(
279 bbox=exposure.getBBox(),
280 wcs=exposure.wcs,
281 filterName=band,
282 epoch=epoch,
283 )
285 refSelection = self.referenceSelector.run(loadResult.refCat, exposure=exposure)
286 refCat = refSelection.sourceCat
288 if self.config.doFiducialZeroPointCull:
289 refCat, sourceSelection.sourceCat = self._do_fiducial_zeropoint_culling(
290 band,
291 loadResult.fluxField,
292 refCat, sourceSelection.sourceCat,
293 exposure.visitInfo.getExposureTime()
294 )
296 # Some operations below require catalog contiguity, which is not
297 # guaranteed from the source selector.
298 if not refCat.isContiguous():
299 refCat = refCat.copy(deep=True)
301 if debug.display:
302 frame = int(debug.frame)
303 displayAstrometry(
304 refCat=refCat,
305 sourceCat=sourceSelection.sourceCat,
306 exposure=exposure,
307 bbox=exposure.getBBox(),
308 frame=frame,
309 title="Reference catalog",
310 )
312 result = pipeBase.Struct(matchTolerance=None)
313 maxMatchDistance = np.inf
314 i = 0
315 while (maxMatchDistance > self.config.minMatchDistanceArcSec and i < self.config.maxIter):
316 i += 1
317 try:
318 result = self._matchAndFitWcs(
319 refCat=refCat,
320 sourceCat=sourceCat,
321 goodSourceCat=sourceSelection.sourceCat,
322 refFluxField=loadResult.fluxField,
323 bbox=exposure.getBBox(),
324 wcs=exposure.wcs,
325 exposure=exposure,
326 matchTolerance=result.matchTolerance,
327 )
328 exposure.setWcs(result.wcs)
329 except exceptions.AstrometryError as e:
330 e._metadata['iterations'] = i
331 sourceCat["coord_ra"] = np.nan
332 sourceCat["coord_dec"] = np.nan
333 exposure.setWcs(None)
334 self.log.error("Failure fitting astrometry. %s: %s", type(e).__name__, e)
335 raise
337 result.stats = self._computeMatchStatsOnSky(result.matches)
338 maxMatchDistance = result.stats.maxMatchDist.asArcseconds()
339 distMean = result.stats.distMean.asArcseconds()
340 distStdDev = result.stats.distStdDev.asArcseconds()
341 self.log.info("Astrometric fit iteration %d: found %d matches with mean separation "
342 "= %0.3f +- %0.3f arcsec; max match distance = %0.3f arcsec.",
343 i, len(result.matches), distMean, distStdDev, maxMatchDistance)
345 # If fitter converged, record the scatter in the exposure metadata
346 # even if the fit was deemed a failure according to the value of
347 # the maxMeanDistanceArcsec config.
348 md = exposure.getMetadata()
349 md['SFM_ASTROM_OFFSET_MEAN'] = distMean
350 md['SFM_ASTROM_OFFSET_STD'] = distStdDev
351 md['SFM_ASTROM_OFFSET_MEDIAN'] = result.scatterOnSky.asArcseconds()
353 if self.usedKey:
354 for m in result.matches:
355 m.second.set(self.usedKey, True)
357 matchMeta = self.refObjLoader.getMetadataBox(
358 bbox=exposure.getBBox(),
359 wcs=exposure.wcs,
360 filterName=band,
361 epoch=epoch,
362 )
364 return pipeBase.Struct(
365 refCat=refCat,
366 matches=result.matches,
367 scatterOnSky=result.scatterOnSky,
368 matchMeta=matchMeta,
369 )
371 def check(self, exposure, sourceCat, nMatches):
372 """Validate the astrometric fit against the maxMeanDistance threshold.
374 If the distMean metric does not satisfy the requirement of being less
375 than the value set in config.maxMeanDistanceArcsec, the WCS on the
376 exposure will be set to None and the coordinate values in the
377 source catalog will be set to NaN to reflect a failed astrometric fit.
379 Parameters
380 ----------
381 exposure : `lsst.afw.image.Exposure`
382 The exposure whose astrometric fit is being evaluated.
383 sourceCat : `lsst.afw.table.SourceCatalog`
384 The catalog of sources associated with the exposure.
385 nMatches : `int`
386 The number of matches that were found and used during
387 the astrometric fit (for logging purposes only).
389 Raises
390 ------
391 BadAstrometryFit
392 If the measured mean on-sky distance between the matched source and
393 reference objects is greater than
394 ``self.config.maxMeanDistanceArcsec``.
395 """
396 # Poor quality fits are a failure.
397 md = exposure.getMetadata()
398 distMean = md['SFM_ASTROM_OFFSET_MEAN']
399 distMedian = md['SFM_ASTROM_OFFSET_MEDIAN']
400 if distMean > self.config.maxMeanDistanceArcsec:
401 exception = exceptions.BadAstrometryFit(nMatches=nMatches, distMean=distMean,
402 maxMeanDist=self.config.maxMeanDistanceArcsec,
403 distMedian=distMedian)
404 exposure.setWcs(None)
405 sourceCat["coord_ra"] = np.nan
406 sourceCat["coord_dec"] = np.nan
407 self.log.error(exception)
408 raise exception
409 return
411 @timeMethod
412 def _matchAndFitWcs(self, refCat, sourceCat, goodSourceCat, refFluxField, bbox, wcs, matchTolerance,
413 exposure=None):
414 """Match sources to reference objects and fit a WCS.
416 Parameters
417 ----------
418 refCat : `lsst.afw.table.SimpleCatalog`
419 catalog of reference objects
420 sourceCat : `lsst.afw.table.SourceCatalog`
421 catalog of sources detected on the exposure
422 goodSourceCat : `lsst.afw.table.SourceCatalog`
423 catalog of down-selected good sources detected on the exposure
424 refFluxField : 'str'
425 field of refCat to use for flux
426 bbox : `lsst.geom.Box2I`
427 bounding box of exposure
428 wcs : `lsst.afw.geom.SkyWcs`
429 initial guess for WCS of exposure
430 matchTolerance : `lsst.meas.astrom.MatchTolerance`
431 a MatchTolerance object (or None) specifying
432 internal tolerances to the matcher. See the MatchTolerance
433 definition in the respective matcher for the class definition.
434 exposure : `lsst.afw.image.Exposure`, optional
435 exposure whose WCS is to be fit, or None; used only for the debug
436 display.
438 Returns
439 -------
440 result : `lsst.pipe.base.Struct`
441 Result struct with components:
443 - ``matches``: astrometric matches
444 (`list` of `lsst.afw.table.ReferenceMatch`).
445 - ``wcs``: the fit WCS (lsst.afw.geom.SkyWcs).
446 - ``scatterOnSky`` : median on-sky separation between reference
447 objects and sources in "matches" (`lsst.afw.geom.Angle`).
448 """
449 import lsstDebug
450 debug = lsstDebug.Info(__name__)
452 sourceFluxField = "slot_%sFlux_instFlux" % (self.config.sourceFluxType)
454 matchRes = self.matcher.matchObjectsToSources(
455 refCat=refCat,
456 sourceCat=goodSourceCat,
457 wcs=wcs,
458 sourceFluxField=sourceFluxField,
459 refFluxField=refFluxField,
460 matchTolerance=matchTolerance,
461 bbox=bbox,
462 )
463 self.log.debug("Found %s matches", len(matchRes.matches))
464 if debug.display:
465 frame = int(debug.frame)
466 displayAstrometry(
467 refCat=refCat,
468 sourceCat=matchRes.usableSourceCat,
469 matches=matchRes.matches,
470 exposure=exposure,
471 bbox=bbox,
472 frame=frame + 1,
473 title="Initial WCS",
474 )
476 if self.config.doMagnitudeOutlierRejection:
477 matches = self._removeMagnitudeOutliers(sourceFluxField, refFluxField, matchRes.matches)
478 else:
479 matches = matchRes.matches
481 self.log.debug("Fitting WCS")
482 fitRes = self.wcsFitter.fitWcs(
483 matches=matches,
484 initWcs=wcs,
485 bbox=bbox,
486 refCat=refCat,
487 sourceCat=sourceCat,
488 exposure=exposure,
489 )
490 fitWcs = fitRes.wcs
491 scatterOnSky = fitRes.scatterOnSky
492 if debug.display:
493 frame = int(debug.frame)
494 displayAstrometry(
495 refCat=refCat,
496 sourceCat=matchRes.usableSourceCat,
497 matches=matches,
498 exposure=exposure,
499 bbox=bbox,
500 frame=frame + 2,
501 title=f"Fitter: {self.wcsFitter._DefaultName}",
502 )
504 return pipeBase.Struct(
505 matches=matches,
506 wcs=fitWcs,
507 scatterOnSky=scatterOnSky,
508 matchTolerance=matchRes.matchTolerance,
509 )
511 def _removeMagnitudeOutliers(self, sourceFluxField, refFluxField, matchesIn):
512 """Remove magnitude outliers, computing a simple zeropoint.
514 Parameters
515 ----------
516 sourceFluxField : `str`
517 Field in source catalog for instrumental fluxes.
518 refFluxField : `str`
519 Field in reference catalog for fluxes (nJy).
520 matchesIn : `list` [`lsst.afw.table.ReferenceMatch`]
521 List of source/reference matches input
523 Returns
524 -------
525 matchesOut : `list` [`lsst.afw.table.ReferenceMatch`]
526 List of source/reference matches with magnitude
527 outliers removed.
528 """
529 nMatch = len(matchesIn)
530 sourceMag = np.zeros(nMatch)
531 refMag = np.zeros(nMatch)
532 for i, match in enumerate(matchesIn):
533 sourceMag[i] = -2.5*np.log10(match[1][sourceFluxField])
534 refMag[i] = (match[0][refFluxField]*units.nJy).to_value(units.ABmag)
536 deltaMag = refMag - sourceMag
537 # Protect against negative fluxes and nans in the reference catalog.
538 goodDelta, = np.where(np.isfinite(deltaMag))
539 zp = np.median(deltaMag[goodDelta])
540 # Use median absolute deviation (MAD) for zpSigma.
541 # Also require a minimum scatter to prevent floating-point errors from
542 # rejecting objects in zero-noise tests.
543 zpSigma = np.clip(scipy.stats.median_abs_deviation(deltaMag[goodDelta], scale='normal'),
544 1e-3,
545 None)
547 self.log.info("Rough zeropoint from astrometry matches is %.4f +/- %.4f.",
548 zp, zpSigma)
550 goodStars = goodDelta[(np.abs(deltaMag[goodDelta] - zp)
551 <= self.config.magnitudeOutlierRejectionNSigma*zpSigma)]
553 nOutlier = nMatch - goodStars.size
554 self.log.info("Removed %d magnitude outliers out of %d total astrometry matches.",
555 nOutlier, nMatch)
557 matchesOut = [matchesIn[idx] for idx in goodStars]
559 return matchesOut
561 def _compute_ref_src_filter_diff(self, band, refFluxField, refCat):
562 """Compute the median ref flux - source filter color difference.
564 The median difference in color between the flux field used for
565 selection and that of the observations being calibrated is computed
566 from the values for each in the reference catalog.
568 Parameters
569 ----------
570 band : `str`
571 Bandpass of observed data.
572 refFluxField : `str`
573 Name of the flux field used in the reference catalog.
574 refCat : `lsst.afw.table.SimpleCatalog`
575 Catalog of reference objects from which to compute the color
576 offset.
578 Returns
579 -------
580 refColorDelta : `float`
581 The median color difference.
582 """
583 if band in self.config.filterMap:
584 refFilterNameForBand = self.config.filterMap.get(band) + "_flux"
585 refCatTemp = refCat[(np.isfinite(refCat[refFluxField])
586 & np.isfinite(refCat[refFilterNameForBand]))].copy(deep=True)
587 if len(refCatTemp) < 3:
588 refColorDeltaDefaults = self.config.refColorDeltaDefaults
589 if band in refColorDeltaDefaults:
590 refColorDelta = refColorDeltaDefaults[band]
591 self.log.warning("Not enough reference sources with finite fluxes in %s and %s, "
592 "so can't compute color shift; a default vaulue of %.2f will "
593 "be applied.", refFluxField, refFilterNameForBand, refColorDelta)
594 else:
595 refMag = (refCatTemp[refFluxField]*units.nJy).to_value(units.ABmag)
596 refMagSrcBand = (refCatTemp[refFilterNameForBand]*units.nJy).to_value(units.ABmag)
597 refColorDelta = np.nanmedian(refMag - refMagSrcBand)
598 self.log.info("Adjusting refCat cutoffs for color shift: %s - %s = %.2f.",
599 refFluxField, refFilterNameForBand, refColorDelta)
600 else:
601 refColorDelta = 0.0
602 self.log.warning("Band %s not found in filterMap: %s. No adjustment for filter "
603 "differences between reference and source catalogs attempted.",
604 band, self.config.filterMap)
605 return refColorDelta
607 def _do_fiducial_zeropoint_culling(self, band, refFluxField, refCat, sourceCat, expTime):
608 """Perform a culling of the catalogs to attempt to match their
609 effective magnitude ranges.
611 This uses a fiducial zeropoint along with the exposure time for
612 the observations to compute our best-guess magnitudes. Also,
613 accommodation is made for the median difference in color between
614 the flux field used for selection and that of the observations being
615 calibrated, which is computed from the values for each in the reference
616 catalog.
618 Parameters
619 ----------
620 band : `str`
621 Bandpass of observed data.
622 refFluxField : `str`
623 Name of the flux field used in the reference catalog.
624 refCat : `lsst.afw.table.SimpleCatalog`
625 Catalog of reference objects to be passed to the matcher. Modified
626 in place.
627 sourceCat : `lsst.afw.table.SourceCatalog`
628 Catalog of observed sources to be passed to the matcher. Modified
629 in place.
630 expTime : `float`
631 Exposure time of the observation being calibrated.
633 Returns
634 -------
635 refColorDelta : `float`
636 The median color difference.
637 """
638 nRefCatPreCull = len(refCat)
639 nSelectedSourceCat = len(sourceCat)
640 # Compute rough limiting magnitudes of selected sources
641 psfFlux = sourceCat["base_PsfFlux_instFlux"]
642 fiducialZeroPoint = self.config.fiducialZeroPoint[band]
643 psfMag = -2.5*np.log10(psfFlux) + fiducialZeroPoint + 2.5*np.log10(expTime)
644 sourceMagMin = min(psfMag) - self.config.cullMagBuffer
645 sourceMagMax = max(psfMag) + self.config.cullMagBuffer
647 # Try to account for median ref flux - source band color difference.
648 refColorDelta = 0.0
649 refColorDelta = self._compute_ref_src_filter_diff(band, refFluxField, refCat)
650 if refColorDelta > self.config.cullMagBuffer:
651 # Start shifting by just half of the median color difference.
652 sourceMagMin += 0.5*refColorDelta
653 sourceMagMax += 0.5*refColorDelta
654 if refColorDelta > 3.0*self.config.cullMagBuffer:
655 # Shift even further if color difference is large compared
656 # with the cullMagBuffer.
657 sourceMagMin = np.nanmin((refCat[refFluxField]*units.nJy).to_value(units.ABmag))
658 sourceMagMax += 0.5*refColorDelta
659 if refColorDelta < -1.0*self.config.cullMagBuffer:
660 # If the color difference is negative (i.e. sources are fainter
661 # in the observed filter), allow full bright end, and shift to
662 # fainter limit.
663 sourceMagMin = np.nanmin((refCat[refFluxField]*units.nJy).to_value(units.ABmag))
664 sourceMagMax += refColorDelta
665 self.log.debug("Number of sources = %d; Number of refs = %d; refs/source = %.2f.",
666 nSelectedSourceCat, nRefCatPreCull, nRefCatPreCull/nSelectedSourceCat)
668 # Include full bright end of reference catalog if there are very
669 # few sources compared to references loaded (this often occurs when
670 # there is a large amount of exctinction and/or scattering from
671 # from bright stars).
672 if nRefCatPreCull/nSelectedSourceCat > self.config.maxRefToSourceRatio:
673 sourceMagMin = np.nanmin((refCat[refFluxField]*units.nJy).to_value(units.ABmag))
675 self.log.info("Source selection: sourceMag (min, max) = (%.3f, %.3f)", sourceMagMin, sourceMagMax)
676 refCat = refCat[np.isfinite(refCat[refFluxField])]
677 refMag = (refCat[refFluxField]*units.nJy).to_value(units.ABmag)
678 refCat = refCat[(refMag < sourceMagMax) & (refMag > sourceMagMin)]
679 refMagMin = np.nanmin(refMag)
680 refMagMax = np.nanmax(refMag)
681 # Now make sure source cat doesn't extend beyond refCat limits.
682 goodSources = ((psfMag < refMagMax - refColorDelta) & (psfMag > refMagMin - refColorDelta))
683 if len(goodSources) < np.sum(goodSources):
684 sourceCat = sourceCat[goodSources].copy(deep=True)
685 nSelectedSourceCat = len(sourceCat)
686 self.log.debug("Number of sources = %d; Number of refs = %d; refs/source = %.2f.",
687 nSelectedSourceCat, nRefCatPreCull, nRefCatPreCull/nSelectedSourceCat)
689 self.log.info("Final: Selected %d/%d reference sources based on fiducial zeropoint culling. "
690 "Mag range in %s = (%.2f, %.2f)", len(refCat), nRefCatPreCull,
691 refFluxField, refMagMin, refMagMax)
692 return refCat, sourceCat