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