349 minMatchedPairs, match_tolerance, sourceFluxField, verbose):
350 """Implementation of matching sources to position reference objects
351
352 Unlike matchObjectsToSources, this method does not check if the sources
353 are suitable.
354
355 Parameters
356 ----------
357 refCat : `lsst.afw.table.SimpleCatalog`
358 catalog of position reference objects that overlap an exposure
359 sourceCat : `lsst.afw.table.SourceCatalog`
360 catalog of sources found on the exposure
361 wcs : `lsst.afw.geom.SkyWcs`
362 estimated WCS of exposure
363 refFluxField : `str`
364 field of refCat to use for flux
365 numUsableSources : `int`
366 number of usable sources (sources with known centroid that are not
367 near the edge, but may be saturated)
368 minMatchedPairs : `int`
369 minimum number of matches
370 match_tolerance : `lsst.meas.astrom.MatchTolerancePessimistic`
371 a MatchTolerance object containing variables specifying matcher
372 tolerances and state from possible previous runs.
373 sourceFluxField : `str`
374 Name of the flux field in the source catalog.
375 verbose : `bool`
376 Set true to print diagnostic information to std::cout
377
378 Returns
379 -------
380 result :
381 Results struct with components:
382
383 - ``matches`` : a list the matches found
384 (`list` of `lsst.afw.table.ReferenceMatch`).
385 - ``match_tolerance`` : MatchTolerance containing updated values from
386 this fit iteration (`lsst.meas.astrom.MatchTolerancePessimistic`)
387 """
388
389
390
391
392
393
394
395 src_array = np.empty((len(sourceCat), 4), dtype=np.float64)
396 for src_idx, srcObj in enumerate(sourceCat):
397 coord = wcs.pixelToSky(srcObj.getCentroid())
398 theta = np.pi / 2 - coord.getLatitude().asRadians()
399 phi = coord.getLongitude().asRadians()
400 flux = srcObj[sourceFluxField]
401 src_array[src_idx, :] = \
402 self._latlong_flux_to_xyz_mag(theta, phi, flux)
403
404 if match_tolerance.PPMbObj is None or \
405 match_tolerance.autoMaxMatchDist is None:
406
407
408
409 ref_array = np.empty((len(refCat), 4), dtype=np.float64)
410 for ref_idx, refObj in enumerate(refCat):
411 theta = np.pi / 2 - refObj.getDec().asRadians()
412 phi = refObj.getRa().asRadians()
413 flux = refObj[refFluxField]
414 ref_array[ref_idx, :] = \
415 self._latlong_flux_to_xyz_mag(theta, phi, flux)
416
417 match_tolerance.PPMbObj = PessimisticPatternMatcherB(
418 ref_array[:, :3], self.log)
419 self.log.debug("Computing source statistics...")
420 maxMatchDistArcSecSrc = self._get_pair_pattern_statistics(
421 src_array)
422 self.log.debug("Computing reference statistics...")
423 maxMatchDistArcSecRef = self._get_pair_pattern_statistics(
424 ref_array)
425 maxMatchDistArcSec = np.max((
426 self.config.minMatchDistPixels
427 * wcs.getPixelScale().asArcseconds(),
428 np.min((maxMatchDistArcSecSrc,
429 maxMatchDistArcSecRef))))
430 match_tolerance.autoMaxMatchDist =
geom.Angle(
431 maxMatchDistArcSec, geom.arcseconds)
432
433
434
435 if match_tolerance.maxShift is None:
436 maxShiftArcseconds = (self.config.maxOffsetPix
437 * wcs.getPixelScale().asArcseconds())
438 else:
439
440
441 maxShiftArcseconds = np.max(
442 (match_tolerance.maxShift.asArcseconds(),
443 self.config.minMatchDistPixels
444 * wcs.getPixelScale().asArcseconds()))
445
446
447
448
449
450 if match_tolerance.maxMatchDist is None:
451 match_tolerance.maxMatchDist = match_tolerance.autoMaxMatchDist
452 else:
453 maxMatchDistArcSec = np.max(
454 (self.config.minMatchDistPixels
455 * wcs.getPixelScale().asArcseconds(),
456 np.min((match_tolerance.maxMatchDist.asArcseconds(),
457 match_tolerance.autoMaxMatchDist.asArcseconds()))))
458
459
460
461
462
463 numConsensus = self.config.numPatternConsensus
464 if len(refCat) < self.config.numRefRequireConsensus:
465 minObjectsForConsensus = \
466 self.config.numBrightStars + \
467 self.config.numPointsForShapeAttempt
468 if len(refCat) < minObjectsForConsensus or \
469 len(sourceCat) < minObjectsForConsensus:
470 numConsensus = 1
471
472 self.log.debug("Current tol maxDist: %.4f arcsec",
473 maxMatchDistArcSec)
474 self.log.debug("Current shift: %.4f arcsec",
475 maxShiftArcseconds)
476
477 match_found = False
478
479 for soften_dist in range(self.config.matcherIterations):
480 if soften_dist == 0 and \
481 match_tolerance.lastMatchedPattern is not None:
482
483
484
485
486 run_n_consent = 1
487 else:
488
489
490 run_n_consent = numConsensus
491
492
493 matcher_struct = match_tolerance.PPMbObj.match(
494 source_array=src_array,
495 n_check=self.config.numPointsForShapeAttempt,
496 n_match=self.config.numPointsForShape,
497 n_agree=run_n_consent,
498 max_n_patterns=self.config.numBrightStars,
499 max_shift=maxShiftArcseconds,
500 max_rotation=self.config.maxRotationDeg,
501 max_dist=maxMatchDistArcSec * 2. ** soften_dist,
502 min_matches=minMatchedPairs,
503 pattern_skip_array=np.array(
504 match_tolerance.failedPatternList)
505 )
506
507 if soften_dist == 0 and \
508 len(matcher_struct.match_ids) == 0 and \
509 match_tolerance.lastMatchedPattern is not None:
510
511
512
513
514
515
516 match_tolerance.failedPatternList.append(
517 match_tolerance.lastMatchedPattern)
518 match_tolerance.lastMatchedPattern = None
519 maxShiftArcseconds = \
520 self.config.maxOffsetPix * wcs.getPixelScale().asArcseconds()
521 elif len(matcher_struct.match_ids) > 0:
522
523
524 match_tolerance.maxShift = \
525 matcher_struct.shift * geom.arcseconds
526 match_tolerance.lastMatchedPattern = \
527 matcher_struct.pattern_idx
528 match_found = True
529 break
530
531
532 if not match_found:
533 return pipeBase.Struct(
534 matches=[],
535 match_tolerance=match_tolerance,
536 )
537
538
539
540
541
542
543
544
545
546
547 distances_arcsec = np.degrees(matcher_struct.distances_rad) * 3600
548 dist_cut_arcsec = np.max(
549 (np.degrees(matcher_struct.max_dist_rad) * 3600,
550 self.config.minMatchDistPixels * wcs.getPixelScale().asArcseconds()))
551
552
553
554 matches = []
555 for match_id_pair, dist_arcsec in zip(matcher_struct.match_ids,
556 distances_arcsec):
557 if dist_arcsec < dist_cut_arcsec:
559 match.first = refCat[int(match_id_pair[1])]
560 match.second = sourceCat[int(match_id_pair[0])]
561
562
563
564 match.distance = match.first.getCoord().separation(
565 match.second.getCoord()).asArcseconds()
566 matches.append(match)
567
568 return pipeBase.Struct(
569 matches=matches,
570 match_tolerance=match_tolerance,
571 )
572
Tag types used to declare specialized field types.
A class representing an angle.