LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
matchPessimisticB.py
Go to the documentation of this file.
1 
2 import numpy as np
3 from scipy.spatial import cKDTree
4 
5 import lsst.pex.config as pexConfig
6 import lsst.pipe.base as pipeBase
7 import lsst.geom as geom
8 import lsst.afw.table as afwTable
9 
10 from .matchOptimisticBTask import MatchTolerance
11 
12 from .pessimistic_pattern_matcher_b_3D import PessimisticPatternMatcherB
13 
14 __all__ = ["MatchPessimisticBTask", "MatchPessimisticBConfig",
15  "MatchTolerancePessimistic"]
16 
17 
19  """Stores match tolerances for use in AstrometryTask and later
20  iterations of the matcher.
21 
22  MatchPessimisticBTask relies on several state variables to be
23  preserved over different iterations in the
24  AstrometryTask.matchAndFitWcs loop of AstrometryTask.
25 
26  Parameters
27  ----------
28  maxMatchDist : `lsst.geom.Angle`
29  Maximum distance to consider a match from the previous match/fit
30  iteration.
31  autoMaxMatchDist : `lsst.geom.Angle`
32  Automated estimation of the maxMatchDist from the sky statistics of the
33  source and reference catalogs.
34  maxShift : `lsst.geom.Angle`
35  Maximum shift found in the previous match/fit cycle.
36  lastMatchedPattern : `int`
37  Index of the last source pattern that was matched into the reference
38  data.
39  failedPatternList : `list` of `int`
40  Previous matches were found to be false positives.
41  PPMbObj : `lsst.meas.astrom.PessimisticPatternMatcherB`
42  Initialized Pessimistic pattern matcher object. Storing this prevents
43  the need for recalculation of the searchable distances in the PPMB.
44  """
45 
46  def __init__(self, maxMatchDist=None, autoMaxMatchDist=None,
47  maxShift=None, lastMatchedPattern=None,
48  failedPatternList=None, PPMbObj=None):
49  self.maxMatchDistmaxMatchDistmaxMatchDist = maxMatchDist
50  self.autoMaxMatchDistautoMaxMatchDist = autoMaxMatchDist
51  self.maxShiftmaxShift = maxShift
52  self.lastMatchedPatternlastMatchedPattern = lastMatchedPattern
53  self.PPMbObjPPMbObj = PPMbObj
54  if failedPatternList is None:
55  self.failedPatternListfailedPatternList = []
56  else:
57  self.failedPatternListfailedPatternList = failedPatternList
58 
59 
60 class MatchPessimisticBConfig(pexConfig.Config):
61  """Configuration for MatchPessimisticBTask
62  """
63  numBrightStars = pexConfig.RangeField(
64  doc="Number of bright stars to use. Sets the max number of patterns "
65  "that can be tested.",
66  dtype=int,
67  default=200,
68  min=2,
69  )
70  minMatchedPairs = pexConfig.RangeField(
71  doc="Minimum number of matched pairs; see also minFracMatchedPairs.",
72  dtype=int,
73  default=30,
74  min=2,
75  )
76  minFracMatchedPairs = pexConfig.RangeField(
77  doc="Minimum number of matched pairs as a fraction of the smaller of "
78  "the number of reference stars or the number of good sources; "
79  "the actual minimum is the smaller of this value or "
80  "minMatchedPairs.",
81  dtype=float,
82  default=0.3,
83  min=0,
84  max=1,
85  )
86  matcherIterations = pexConfig.RangeField(
87  doc="Number of softening iterations in matcher.",
88  dtype=int,
89  default=5,
90  min=1,
91  )
92  maxOffsetPix = pexConfig.RangeField(
93  doc="Maximum allowed shift of WCS, due to matching (pixel). "
94  "When changing this value, the "
95  "LoadReferenceObjectsConfig.pixelMargin should also be updated.",
96  dtype=int,
97  default=300,
98  max=4000,
99  )
100  maxRotationDeg = pexConfig.RangeField(
101  doc="Rotation angle allowed between sources and position reference "
102  "objects (degrees).",
103  dtype=float,
104  default=1.0,
105  max=6.0,
106  )
107  numPointsForShape = pexConfig.Field(
108  doc="Number of points to define a shape for matching.",
109  dtype=int,
110  default=6,
111  )
112  numPointsForShapeAttempt = pexConfig.Field(
113  doc="Number of points to try for creating a shape. This value should "
114  "be greater than or equal to numPointsForShape. Besides "
115  "loosening the signal to noise cut in the 'matcher' SourceSelector, "
116  "increasing this number will solve CCDs where no match was found.",
117  dtype=int,
118  default=6,
119  )
120  minMatchDistPixels = pexConfig.RangeField(
121  doc="Distance in units of pixels to always consider a source-"
122  "reference pair a match. This prevents the astrometric fitter "
123  "from over-fitting and removing stars that should be matched and "
124  "allows for inclusion of new matches as the wcs improves.",
125  dtype=float,
126  default=1.0,
127  min=0.0,
128  max=6.0,
129  )
130  numPatternConsensus = pexConfig.Field(
131  doc="Number of implied shift/rotations from patterns that must agree "
132  "before it a given shift/rotation is accepted. This is only used "
133  "after the first softening iteration fails and if both the "
134  "number of reference and source objects is greater than "
135  "numBrightStars.",
136  dtype=int,
137  default=3,
138  )
139  numRefRequireConsensus = pexConfig.Field(
140  doc="If the available reference objects exceeds this number, "
141  "consensus/pessimistic mode will enforced regardless of the "
142  "number of available sources. Below this optimistic mode ("
143  "exit at first match rather than requiring numPatternConsensus to "
144  "be matched) can be used. If more sources are required to match, "
145  "decrease the signal to noise cut in the sourceSelector.",
146  dtype=int,
147  default=1000,
148  )
149  maxRefObjects = pexConfig.RangeField(
150  doc="Maximum number of reference objects to use for the matcher. The "
151  "absolute maximum allowed for is 2 ** 16 for memory reasons.",
152  dtype=int,
153  default=2**16,
154  min=0,
155  max=2**16 + 1,
156  )
157 
158  def validate(self):
159  pexConfig.Config.validate(self)
160  if self.numPointsForShapeAttemptnumPointsForShapeAttempt < self.numPointsForShapenumPointsForShape:
161  raise ValueError("numPointsForShapeAttempt must be greater than "
162  "or equal to numPointsForShape.")
163  if self.numPointsForShapenumPointsForShape > self.numBrightStarsnumBrightStars:
164  raise ValueError("numBrightStars must be greater than "
165  "numPointsForShape.")
166 
167 
168 # The following block adds links to this task from the Task Documentation page.
169 # \addtogroup LSST_task_documentation
170 # \{
171 # \page measAstrom_MatchPessimisticBTask
172 # \ref MatchPessimisticBTask "MatchPessimisticBTask"
173 # Match sources to reference objects
174 # \}
175 
176 
177 class MatchPessimisticBTask(pipeBase.Task):
178  """Match sources to reference objects.
179  """
180 
181  ConfigClass = MatchPessimisticBConfig
182  _DefaultName = "matchObjectsToSources"
183 
184  def __init__(self, **kwargs):
185  pipeBase.Task.__init__(self, **kwargs)
186 
187  @pipeBase.timeMethod
188  def matchObjectsToSources(self, refCat, sourceCat, wcs, sourceFluxField, refFluxField,
189  match_tolerance=None):
190  """Match sources to position reference stars
191 
192  refCat : `lsst.afw.table.SimpleCatalog`
193  catalog of reference objects that overlap the exposure; reads
194  fields for:
195 
196  - coord
197  - the specified flux field
198 
199  sourceCat : `lsst.afw.table.SourceCatalog`
200  Catalog of sources found on an exposure. This should already be
201  down-selected to "good"/"usable" sources in the calling Task.
202  wcs : `lsst.afw.geom.SkyWcs`
203  estimated WCS
204  sourceFluxField: `str`
205  field of sourceCat to use for flux
206  refFluxField : `str`
207  field of refCat to use for flux
208  match_tolerance : `lsst.meas.astrom.MatchTolerancePessimistic`
209  is a MatchTolerance class object or `None`. This this class is used
210  to communicate state between AstrometryTask and MatcherTask.
211  AstrometryTask will also set the MatchTolerance class variable
212  maxMatchDist based on the scatter AstrometryTask has found after
213  fitting for the wcs.
214 
215  Returns
216  -------
217  result : `lsst.pipe.base.Struct`
218  Result struct with components:
219 
220  - ``matches`` : source to reference matches found (`list` of
221  `lsst.afw.table.ReferenceMatch`)
222  - ``usableSourceCat`` : a catalog of sources potentially usable for
223  matching and WCS fitting (`lsst.afw.table.SourceCatalog`).
224  - ``match_tolerance`` : a MatchTolerance object containing the
225  resulting state variables from the match
226  (`lsst.meas.astrom.MatchTolerancePessimistic`).
227  """
228  import lsstDebug
229  debug = lsstDebug.Info(__name__)
230 
231  # If we get an empty tolerance struct create the variables we need for
232  # this matcher.
233  if match_tolerance is None:
234  match_tolerance = MatchTolerancePessimistic()
235 
236  # Make a name alias here for consistency with older code, and to make
237  # it clear that this is a good/usable (cleaned) source catalog.
238  goodSourceCat = sourceCat
239 
240  numUsableSources = len(goodSourceCat)
241 
242  if len(goodSourceCat) == 0:
243  raise pipeBase.TaskError("No sources are good")
244 
245  minMatchedPairs = min(self.config.minMatchedPairs,
246  int(self.config.minFracMatchedPairs
247  * min([len(refCat), len(goodSourceCat)])))
248 
249  if len(refCat) > self.config.maxRefObjects:
250  self.log.warn(
251  "WARNING: Reference catalog larger that maximum allowed. "
252  "Trimming to %i" % self.config.maxRefObjects)
253  trimmedRefCat = self._filterRefCat_filterRefCat(refCat, refFluxField)
254  else:
255  trimmedRefCat = refCat
256 
257  doMatchReturn = self._doMatch_doMatch(
258  refCat=trimmedRefCat,
259  sourceCat=goodSourceCat,
260  wcs=wcs,
261  refFluxField=refFluxField,
262  numUsableSources=numUsableSources,
263  minMatchedPairs=minMatchedPairs,
264  match_tolerance=match_tolerance,
265  sourceFluxField=sourceFluxField,
266  verbose=debug.verbose,
267  )
268  matches = doMatchReturn.matches
269  match_tolerance = doMatchReturn.match_tolerance
270 
271  if len(matches) == 0:
272  raise RuntimeError("Unable to match sources")
273 
274  self.log.info("Matched %d sources" % len(matches))
275  if len(matches) < minMatchedPairs:
276  self.log.warn("Number of matches is smaller than request")
277 
278  return pipeBase.Struct(
279  matches=matches,
280  usableSourceCat=goodSourceCat,
281  match_tolerance=match_tolerance,
282  )
283 
284  def _filterRefCat(self, refCat, refFluxField):
285  """Sub-select a number of reference objects starting from the brightest
286  and maxing out at the number specified by maxRefObjects in the config.
287 
288  No trimming is done if len(refCat) > config.maxRefObjects.
289 
290  Parameters
291  ----------
292  refCat : `lsst.afw.table.SimpleCatalog`
293  Catalog of reference objects to trim.
294  refFluxField : `str`
295  field of refCat to use for flux
296  Returns
297  -------
298  outCat : `lsst.afw.table.SimpleCatalog`
299  Catalog trimmed to the number set in the task config from the
300  brightest flux down.
301  """
302  # Find the flux cut that gives us the desired number of objects.
303  if len(refCat) <= self.config.maxRefObjects:
304  return refCat
305  fluxArray = refCat.get(refFluxField)
306  sortedFluxArray = fluxArray[fluxArray.argsort()]
307  minFlux = sortedFluxArray[-(self.config.maxRefObjects + 1)]
308 
309  selected = (refCat.get(refFluxField) > minFlux)
310 
311  outCat = afwTable.SimpleCatalog(refCat.schema)
312  outCat.reserve(self.config.maxRefObjects)
313  outCat.extend(refCat[selected])
314 
315  return outCat
316 
317  @pipeBase.timeMethod
318  def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources,
319  minMatchedPairs, match_tolerance, sourceFluxField, verbose):
320  """Implementation of matching sources to position reference objects
321 
322  Unlike matchObjectsToSources, this method does not check if the sources
323  are suitable.
324 
325  Parameters
326  ----------
327  refCat : `lsst.afw.table.SimpleCatalog`
328  catalog of position reference objects that overlap an exposure
329  sourceCat : `lsst.afw.table.SourceCatalog`
330  catalog of sources found on the exposure
331  wcs : `lsst.afw.geom.SkyWcs`
332  estimated WCS of exposure
333  refFluxField : `str`
334  field of refCat to use for flux
335  numUsableSources : `int`
336  number of usable sources (sources with known centroid that are not
337  near the edge, but may be saturated)
338  minMatchedPairs : `int`
339  minimum number of matches
340  match_tolerance : `lsst.meas.astrom.MatchTolerancePessimistic`
341  a MatchTolerance object containing variables specifying matcher
342  tolerances and state from possible previous runs.
343  sourceFluxField : `str`
344  Name of the flux field in the source catalog.
345  verbose : `bool`
346  Set true to print diagnostic information to std::cout
347 
348  Returns
349  -------
350  result :
351  Results struct with components:
352 
353  - ``matches`` : a list the matches found
354  (`list` of `lsst.afw.table.ReferenceMatch`).
355  - ``match_tolerance`` : MatchTolerance containing updated values from
356  this fit iteration (`lsst.meas.astrom.MatchTolerancePessimistic`)
357  """
358 
359  # Load the source and reference catalog as spherical points
360  # in numpy array. We do this rather than relying on internal
361  # lsst C objects for simplicity and because we require
362  # objects contiguous in memory. We need to do these slightly
363  # differently for the reference and source cats as they are
364  # different catalog objects with different fields.
365  src_array = np.empty((len(sourceCat), 4), dtype=np.float64)
366  for src_idx, srcObj in enumerate(sourceCat):
367  coord = wcs.pixelToSky(srcObj.getCentroid())
368  theta = np.pi / 2 - coord.getLatitude().asRadians()
369  phi = coord.getLongitude().asRadians()
370  flux = srcObj[sourceFluxField]
371  src_array[src_idx, :] = \
372  self._latlong_flux_to_xyz_mag_latlong_flux_to_xyz_mag(theta, phi, flux)
373 
374  if match_tolerance.PPMbObj is None or \
375  match_tolerance.autoMaxMatchDist is None:
376  # The reference catalog is fixed per AstrometryTask so we only
377  # create the data needed if this is the first step in the match
378  # fit cycle.
379  ref_array = np.empty((len(refCat), 4), dtype=np.float64)
380  for ref_idx, refObj in enumerate(refCat):
381  theta = np.pi / 2 - refObj.getDec().asRadians()
382  phi = refObj.getRa().asRadians()
383  flux = refObj[refFluxField]
384  ref_array[ref_idx, :] = \
385  self._latlong_flux_to_xyz_mag_latlong_flux_to_xyz_mag(theta, phi, flux)
386  # Create our matcher object.
387  match_tolerance.PPMbObj = PessimisticPatternMatcherB(
388  ref_array[:, :3], self.log)
389  self.log.debug("Computing source statistics...")
390  maxMatchDistArcSecSrc = self._get_pair_pattern_statistics_get_pair_pattern_statistics(
391  src_array)
392  self.log.debug("Computing reference statistics...")
393  maxMatchDistArcSecRef = self._get_pair_pattern_statistics_get_pair_pattern_statistics(
394  ref_array)
395  maxMatchDistArcSec = np.max((
396  self.config.minMatchDistPixels
397  * wcs.getPixelScale().asArcseconds(),
398  np.min((maxMatchDistArcSecSrc,
399  maxMatchDistArcSecRef))))
400  match_tolerance.autoMaxMatchDist = geom.Angle(
401  maxMatchDistArcSec, geom.arcseconds)
402 
403  # Set configurable defaults when we encounter None type or set
404  # state based on previous run of AstrometryTask._matchAndFitWcs.
405  if match_tolerance.maxShift is None:
406  maxShiftArcseconds = (self.config.maxOffsetPix
407  * wcs.getPixelScale().asArcseconds())
408  else:
409  # We don't want to clamp down too hard on the allowed shift so
410  # we test that the smallest we ever allow is the pixel scale.
411  maxShiftArcseconds = np.max(
412  (match_tolerance.maxShift.asArcseconds(),
413  self.config.minMatchDistPixels
414  * wcs.getPixelScale().asArcseconds()))
415 
416  # If our tolerances are not set from a previous run, estimate a
417  # starting tolerance guess from the statistics of patterns we can
418  # create on both the source and reference catalog. We use the smaller
419  # of the two.
420  if match_tolerance.maxMatchDist is None:
421  match_tolerance.maxMatchDist = match_tolerance.autoMaxMatchDist
422  else:
423  maxMatchDistArcSec = np.max(
424  (self.config.minMatchDistPixels
425  * wcs.getPixelScale().asArcseconds(),
426  np.min((match_tolerance.maxMatchDist.asArcseconds(),
427  match_tolerance.autoMaxMatchDist.asArcseconds()))))
428 
429  # Make sure the data we are considering is dense enough to require
430  # the consensus mode of the matcher. If not default to Optimistic
431  # pattern matcher behavior. We enforce pessimistic mode if the
432  # reference cat is sufficiently large, avoiding false positives.
433  numConsensus = self.config.numPatternConsensus
434  if len(refCat) < self.config.numRefRequireConsensus:
435  minObjectsForConsensus = \
436  self.config.numBrightStars + \
437  self.config.numPointsForShapeAttempt
438  if len(refCat) < minObjectsForConsensus or \
439  len(sourceCat) < minObjectsForConsensus:
440  numConsensus = 1
441 
442  self.log.debug("Current tol maxDist: %.4f arcsec" %
443  maxMatchDistArcSec)
444  self.log.debug("Current shift: %.4f arcsec" %
445  maxShiftArcseconds)
446 
447  match_found = False
448  # Start the iteration over our tolerances.
449  for soften_dist in range(self.config.matcherIterations):
450  if soften_dist == 0 and \
451  match_tolerance.lastMatchedPattern is not None:
452  # If we are on the first, most stringent tolerance,
453  # and have already found a match, the matcher should behave
454  # like an optimistic pattern matcher. Exiting at the first
455  # match.
456  run_n_consent = 1
457  else:
458  # If we fail or this is the first match attempt, set the
459  # pattern consensus to the specified config value.
460  run_n_consent = numConsensus
461  # We double the match dist tolerance each round and add 1 to the
462  # to the number of candidate spokes to check.
463  matcher_struct = match_tolerance.PPMbObj.match(
464  source_array=src_array,
465  n_check=self.config.numPointsForShapeAttempt,
466  n_match=self.config.numPointsForShape,
467  n_agree=run_n_consent,
468  max_n_patterns=self.config.numBrightStars,
469  max_shift=maxShiftArcseconds,
470  max_rotation=self.config.maxRotationDeg,
471  max_dist=maxMatchDistArcSec * 2. ** soften_dist,
472  min_matches=minMatchedPairs,
473  pattern_skip_array=np.array(
474  match_tolerance.failedPatternList)
475  )
476 
477  if soften_dist == 0 and \
478  len(matcher_struct.match_ids) == 0 and \
479  match_tolerance.lastMatchedPattern is not None:
480  # If we found a pattern on a previous match-fit iteration and
481  # can't find an optimistic match on our first try with the
482  # tolerances as found in the previous match-fit,
483  # the match we found in the last iteration was likely bad. We
484  # append the bad match's index to the a list of
485  # patterns/matches to skip on subsequent iterations.
486  match_tolerance.failedPatternList.append(
487  match_tolerance.lastMatchedPattern)
488  match_tolerance.lastMatchedPattern = None
489  maxShiftArcseconds = \
490  self.config.maxOffsetPix * wcs.getPixelScale().asArcseconds()
491  elif len(matcher_struct.match_ids) > 0:
492  # Match found, save a bit a state regarding this pattern
493  # in the match tolerance class object and exit.
494  match_tolerance.maxShift = \
495  matcher_struct.shift * geom.arcseconds
496  match_tolerance.lastMatchedPattern = \
497  matcher_struct.pattern_idx
498  match_found = True
499  break
500 
501  # If we didn't find a match, exit early.
502  if not match_found:
503  return pipeBase.Struct(
504  matches=[],
505  match_tolerance=match_tolerance,
506  )
507 
508  # The matcher returns all the nearest neighbors that agree between
509  # the reference and source catalog. For the current astrometric solver
510  # we need to remove as many false positives as possible before sending
511  # the matches off to the solver. The low value of 100 and high value of
512  # 2 are the low number of sigma and high respectively. The exact values
513  # were found after testing on data of various reference/source
514  # densities and astrometric distortion quality, specifically the
515  # visits: HSC (3358), DECam (406285, 410827),
516  # CFHT (793169, 896070, 980526).
517  distances_arcsec = np.degrees(matcher_struct.distances_rad) * 3600
518  dist_cut_arcsec = np.max(
519  (np.degrees(matcher_struct.max_dist_rad) * 3600,
520  self.config.minMatchDistPixels * wcs.getPixelScale().asArcseconds()))
521 
522  # A match has been found, return our list of matches and
523  # return.
524  matches = []
525  for match_id_pair, dist_arcsec in zip(matcher_struct.match_ids,
526  distances_arcsec):
527  if dist_arcsec < dist_cut_arcsec:
528  match = afwTable.ReferenceMatch()
529  match.first = refCat[int(match_id_pair[1])]
530  match.second = sourceCat[int(match_id_pair[0])]
531  # We compute the true distance along and sphere. This isn't
532  # used in the WCS fitter however it is used in the unittest
533  # to confirm the matches computed.
534  match.distance = match.first.getCoord().separation(
535  match.second.getCoord()).asArcseconds()
536  matches.append(match)
537 
538  return pipeBase.Struct(
539  matches=matches,
540  match_tolerance=match_tolerance,
541  )
542 
543  def _latlong_flux_to_xyz_mag(self, theta, phi, flux):
544  """Convert angles theta and phi and a flux into unit sphere
545  x, y, z, and a relative magnitude.
546 
547  Takes in a afw catalog object and converts the catalog object RA, DECs
548  to points on the unit sphere. Also converts the flux into a simple,
549  non-zero-pointed magnitude for relative sorting.
550 
551  Parameters
552  ----------
553  theta : `float`
554  Angle from the north pole (z axis) of the sphere
555  phi : `float`
556  Rotation around the sphere
557 
558  Return
559  ------
560  output_array : `numpy.ndarray`, (N, 4)
561  Spherical unit vector x, y, z with flux.
562  """
563  output_array = np.empty(4, dtype=np.float64)
564  output_array[0] = np.sin(theta)*np.cos(phi)
565  output_array[1] = np.sin(theta)*np.sin(phi)
566  output_array[2] = np.cos(theta)
567  if flux > 0:
568  output_array[3] = -2.5 * np.log10(flux)
569  else:
570  # Set flux to a very faint mag if its for some reason it
571  # does not exist
572  output_array[3] = 99.
573 
574  return output_array
575 
576  def _get_pair_pattern_statistics(self, cat_array):
577  """ Compute the tolerances for the matcher automatically by comparing
578  pinwheel patterns as we would in the matcher.
579 
580  We test how similar the patterns we can create from a given set of
581  objects by computing the spoke lengths for each pattern and sorting
582  them from smallest to largest. The match tolerance is the average
583  distance per spoke between the closest two patterns in the sorted
584  spoke length space.
585 
586  Parameters
587  ----------
588  cat_array : `numpy.ndarray`, (N, 3)
589  array of 3 vectors representing the x, y, z position of catalog
590  objects on the unit sphere.
591 
592  Returns
593  -------
594  dist_tol : `float`
595  Suggested max match tolerance distance calculated from comparisons
596  between pinwheel patterns used in optimistic/pessimistic pattern
597  matcher.
598  """
599 
600  self.log.debug("Starting automated tolerance calculation...")
601 
602  # Create an empty array of all the patterns we possibly make
603  # sorting from brightest to faintest.
604  pattern_array = np.empty(
605  (cat_array.shape[0] - self.config.numPointsForShape,
606  self.config.numPointsForShape - 1))
607  flux_args_array = np.argsort(cat_array[:, -1])
608 
609  # Sort our input array.
610  tmp_sort_array = cat_array[flux_args_array]
611 
612  # Start making patterns.
613  for start_idx in range(cat_array.shape[0]
614  - self.config.numPointsForShape):
615  pattern_points = tmp_sort_array[start_idx:start_idx
616  + self.config.numPointsForShape, :-1]
617  pattern_delta = pattern_points[1:, :] - pattern_points[0, :]
618  pattern_array[start_idx, :] = np.sqrt(
619  pattern_delta[:, 0] ** 2
620  + pattern_delta[:, 1] ** 2
621  + pattern_delta[:, 2] ** 2)
622 
623  # When we store the length of each spoke in our pattern we
624  # sort from shortest to longest so we have a defined space
625  # to compare them in.
626  pattern_array[start_idx, :] = pattern_array[
627  start_idx, np.argsort(pattern_array[start_idx, :])]
628 
629  # Create a searchable tree object of the patterns and find
630  # for any given pattern the closest pattern in the sorted
631  # spoke length space.
632  dist_tree = cKDTree(
633  pattern_array[:, :(self.config.numPointsForShape - 1)])
634  dist_nearest_array, ids = dist_tree.query(
635  pattern_array[:, :(self.config.numPointsForShape - 1)], k=2)
636  dist_nearest_array = dist_nearest_array[:, 1]
637  dist_nearest_array.sort()
638 
639  # We use the two closest patterns to set our tolerance.
640  dist_idx = 0
641  dist_tol = (np.degrees(dist_nearest_array[dist_idx]) * 3600.
642  / (self.config.numPointsForShape - 1.))
643 
644  self.log.debug("Automated tolerance")
645  self.log.debug("\tdistance/match tol: %.4f [arcsec]" % dist_tol)
646 
647  return dist_tol
int min
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
Definition: SortedCatalog.h:42
A class representing an angle.
Definition: Angle.h:127
def matchObjectsToSources(self, refCat, sourceCat, wcs, sourceFluxField, refFluxField, match_tolerance=None)
def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources, minMatchedPairs, match_tolerance, sourceFluxField, verbose)
def __init__(self, maxMatchDist=None, autoMaxMatchDist=None, maxShift=None, lastMatchedPattern=None, failedPatternList=None, PPMbObj=None)
Lightweight representation of a geometric match between two records.
Definition: Match.h:67