210 matchTolerance=None):
211 """Match sources to position reference stars
213 refCat : `lsst.afw.table.SimpleCatalog`
214 catalog of reference objects that overlap the exposure; reads
218 - the specified flux field
220 sourceCat : `lsst.afw.table.SourceCatalog`
221 Catalog of sources found on an exposure. This should already be
222 down-selected to "good"/"usable" sources in the calling Task.
223 wcs : `lsst.afw.geom.SkyWcs`
225 sourceFluxField: `str`
226 field of sourceCat to use for flux
228 field of refCat to use for flux
229 matchTolerance : `lsst.meas.astrom.MatchTolerancePessimistic`
230 is a MatchTolerance class object or `None`. This this class is used
231 to communicate state between AstrometryTask and MatcherTask.
232 AstrometryTask will also set the MatchTolerance class variable
233 maxMatchDist based on the scatter AstrometryTask has found after
238 result : `lsst.pipe.base.Struct`
239 Result struct with components:
241 - ``matches`` : source to reference matches found (`list` of
242 `lsst.afw.table.ReferenceMatch`)
243 - ``usableSourceCat`` : a catalog of sources potentially usable for
244 matching and WCS fitting (`lsst.afw.table.SourceCatalog`).
245 - ``matchTolerance`` : a MatchTolerance object containing the
246 resulting state variables from the match
247 (`lsst.meas.astrom.MatchTolerancePessimistic`).
254 if matchTolerance
is None:
259 goodSourceCat = sourceCat
261 if (numUsableSources := len(goodSourceCat)) == 0:
264 minMatchedPairs =
min(self.config.minMatchedPairs,
265 int(self.config.minFracMatchedPairs
266 *
min([len(refCat), len(goodSourceCat)])))
268 if len(goodSourceCat) <= self.config.numPointsForShape:
269 msg = (f
"Not enough catalog objects ({len(goodSourceCat)}) to make a "
270 f
"shape for the matcher (need {self.config.numPointsForShape}).")
272 if len(refCat) <= self.config.numPointsForShape:
273 msg = (f
"Not enough refcat objects ({len(refCat)}) to make a "
274 f
"shape for the matcher (need {self.config.numPointsForShape}).")
277 if len(refCat) > self.config.maxRefObjects:
279 "WARNING: Reference catalog larger than maximum allowed. "
280 "Trimming to %i", self.config.maxRefObjects)
283 trimmedRefCat = refCat
286 refCat=trimmedRefCat,
287 sourceCat=goodSourceCat,
289 refFluxField=refFluxField,
290 numUsableSources=numUsableSources,
291 minMatchedPairs=minMatchedPairs,
292 matchTolerance=matchTolerance,
293 sourceFluxField=sourceFluxField,
294 verbose=debug.verbose,
296 matches = doMatchReturn.matches
297 matchTolerance = doMatchReturn.matchTolerance
299 if (nMatches := len(matches)) == 0:
302 self.
log.info(
"Matched %d sources", nMatches)
303 if nMatches < minMatchedPairs:
304 self.
log.warning(
"Number of matches (%s) is smaller than minimum requested (%s)",
305 nMatches, minMatchedPairs)
307 return pipeBase.Struct(
309 usableSourceCat=goodSourceCat,
310 matchTolerance=matchTolerance,
347 def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources,
348 minMatchedPairs, matchTolerance, sourceFluxField, verbose):
349 """Implementation of matching sources to position reference objects
351 Unlike matchObjectsToSources, this method does not check if the sources
356 refCat : `lsst.afw.table.SimpleCatalog`
357 catalog of position reference objects that overlap an exposure
358 sourceCat : `lsst.afw.table.SourceCatalog`
359 catalog of sources found on the exposure
360 wcs : `lsst.afw.geom.SkyWcs`
361 estimated WCS of exposure
363 field of refCat to use for flux
364 numUsableSources : `int`
365 number of usable sources (sources with known centroid that are not
366 near the edge, but may be saturated)
367 minMatchedPairs : `int`
368 minimum number of matches
369 matchTolerance : `lsst.meas.astrom.MatchTolerancePessimistic`
370 a MatchTolerance object containing variables specifying matcher
371 tolerances and state from possible previous runs.
372 sourceFluxField : `str`
373 Name of the flux field in the source catalog.
375 Set true to print diagnostic information to std::cout
380 Results struct with components:
382 - ``matches`` : a list the matches found
383 (`list` of `lsst.afw.table.ReferenceMatch`).
384 - ``matchTolerance`` : MatchTolerance containing updated values from
385 this fit iteration (`lsst.meas.astrom.MatchTolerancePessimistic`)
394 src_array = np.empty((len(sourceCat), 4), dtype=np.float64)
395 for src_idx, srcObj
in enumerate(sourceCat):
396 coord = wcs.pixelToSky(srcObj.getCentroid())
397 theta = np.pi / 2 - coord.getLatitude().asRadians()
398 phi = coord.getLongitude().asRadians()
399 flux = srcObj[sourceFluxField]
400 src_array[src_idx, :] = \
403 if matchTolerance.PPMbObj
is None or \
404 matchTolerance.autoMaxMatchDist
is None:
408 ref_array = np.empty((len(refCat), 4), dtype=np.float64)
409 for ref_idx, refObj
in enumerate(refCat):
410 theta = np.pi / 2 - refObj.getDec().asRadians()
411 phi = refObj.getRa().asRadians()
412 flux = refObj[refFluxField]
413 ref_array[ref_idx, :] = \
417 ref_array[:, :3], self.
log)
418 self.
log.debug(
"Computing source statistics...")
421 self.
log.debug(
"Computing reference statistics...")
424 maxMatchDistArcSec = np.max((
425 self.config.minMatchDistPixels
426 * wcs.getPixelScale().asArcseconds(),
427 np.min((maxMatchDistArcSecSrc,
428 maxMatchDistArcSecRef))))
430 maxMatchDistArcSec, geom.arcseconds)
434 if matchTolerance.maxShift
is None:
435 maxShiftArcseconds = (self.config.maxOffsetPix
436 * wcs.getPixelScale().asArcseconds())
440 maxShiftArcseconds = np.max(
441 (matchTolerance.maxShift.asArcseconds(),
442 self.config.minMatchDistPixels
443 * wcs.getPixelScale().asArcseconds()))
449 if matchTolerance.maxMatchDist
is None:
450 matchTolerance.maxMatchDist = matchTolerance.autoMaxMatchDist
452 maxMatchDistArcSec = np.max(
453 (self.config.minMatchDistPixels
454 * wcs.getPixelScale().asArcseconds(),
455 np.min((matchTolerance.maxMatchDist.asArcseconds(),
456 matchTolerance.autoMaxMatchDist.asArcseconds()))))
462 numConsensus = self.config.numPatternConsensus
463 if len(refCat) < self.config.numRefRequireConsensus:
464 minObjectsForConsensus = \
465 self.config.numBrightStars + \
466 self.config.numPointsForShapeAttempt
467 if len(refCat) < minObjectsForConsensus
or \
468 len(sourceCat) < minObjectsForConsensus:
471 self.
log.debug(
"Current tol maxDist: %.4f arcsec",
473 self.
log.debug(
"Current shift: %.4f arcsec",
478 for soften_dist
in range(self.config.matcherIterations):
479 if soften_dist == 0
and \
480 matchTolerance.lastMatchedPattern
is not None:
489 run_n_consent = numConsensus
492 matcher_struct = matchTolerance.PPMbObj.match(
493 source_array=src_array,
494 n_check=self.config.numPointsForShapeAttempt,
495 n_match=self.config.numPointsForShape,
496 n_agree=run_n_consent,
497 max_n_patterns=self.config.numBrightStars,
498 max_shift=maxShiftArcseconds,
499 max_rotation=self.config.maxRotationDeg,
500 max_dist=maxMatchDistArcSec * 2. ** soften_dist,
501 min_matches=minMatchedPairs,
502 pattern_skip_array=np.array(
503 matchTolerance.failedPatternList)
506 if soften_dist == 0
and \
507 len(matcher_struct.match_ids) == 0
and \
508 matchTolerance.lastMatchedPattern
is not None:
515 matchTolerance.failedPatternList.append(
516 matchTolerance.lastMatchedPattern)
517 matchTolerance.lastMatchedPattern =
None
518 maxShiftArcseconds = \
519 self.config.maxOffsetPix * wcs.getPixelScale().asArcseconds()
520 elif len(matcher_struct.match_ids) > 0:
523 matchTolerance.maxShift = \
524 matcher_struct.shift * geom.arcseconds
525 matchTolerance.lastMatchedPattern = \
526 matcher_struct.pattern_idx
532 return pipeBase.Struct(
534 matchTolerance=matchTolerance,
546 distances_arcsec = np.degrees(matcher_struct.distances_rad) * 3600
547 dist_cut_arcsec = np.max(
548 (np.degrees(matcher_struct.max_dist_rad) * 3600,
549 self.config.minMatchDistPixels * wcs.getPixelScale().asArcseconds()))
554 for match_id_pair, dist_arcsec
in zip(matcher_struct.match_ids,
556 if dist_arcsec < dist_cut_arcsec:
558 match.first = refCat[int(match_id_pair[1])]
559 match.second = sourceCat[int(match_id_pair[0])]
563 match.distance = match.first.getCoord().separation(
564 match.second.getCoord()).asArcseconds()
565 matches.append(match)
567 return pipeBase.Struct(
569 matchTolerance=matchTolerance,
606 """ Compute the tolerances for the matcher automatically by comparing
607 pinwheel patterns as we would in the matcher.
609 We test how similar the patterns we can create from a given set of
610 objects by computing the spoke lengths for each pattern and sorting
611 them from smallest to largest. The match tolerance is the average
612 distance per spoke between the closest two patterns in the sorted
617 cat_array : `numpy.ndarray`, (N, 3)
618 array of 3 vectors representing the x, y, z position of catalog
619 objects on the unit sphere.
624 Suggested max match tolerance distance calculated from comparisons
625 between pinwheel patterns used in optimistic/pessimistic pattern
629 self.
log.debug(
"Starting automated tolerance calculation...")
633 pattern_array = np.empty(
634 (cat_array.shape[0] - self.config.numPointsForShape,
635 self.config.numPointsForShape - 1))
636 flux_args_array = np.argsort(cat_array[:, -1])
639 tmp_sort_array = cat_array[flux_args_array]
642 for start_idx
in range(cat_array.shape[0]
643 - self.config.numPointsForShape):
644 pattern_points = tmp_sort_array[start_idx:start_idx
645 + self.config.numPointsForShape, :-1]
646 pattern_delta = pattern_points[1:, :] - pattern_points[0, :]
647 pattern_array[start_idx, :] = np.sqrt(
648 pattern_delta[:, 0] ** 2
649 + pattern_delta[:, 1] ** 2
650 + pattern_delta[:, 2] ** 2)
655 pattern_array[start_idx, :] = pattern_array[
656 start_idx, np.argsort(pattern_array[start_idx, :])]
662 pattern_array[:, :(self.config.numPointsForShape - 1)])
663 dist_nearest_array, ids = dist_tree.query(
664 pattern_array[:, :(self.config.numPointsForShape - 1)], k=2)
665 dist_nearest_array = dist_nearest_array[:, 1]
666 dist_nearest_array.sort()
670 dist_tol = (np.degrees(dist_nearest_array[dist_idx]) * 3600.
671 / (self.config.numPointsForShape - 1.))
673 self.
log.debug(
"Automated tolerance")
674 self.
log.debug(
"\tdistance/match tol: %.4f [arcsec]", dist_tol)