266 def selectSources(self, sourceCat, matches=None, exposure=None):
267 """Return a selection of psf-like objects.
271 sourceCat : `lsst.afw.table.SourceCatalog`
272 Catalog of sources to select from.
273 This catalog must be contiguous in memory.
274 matches : `list` of `lsst.afw.table.ReferenceMatch` or None
275 Ignored by this source selector.
276 exposure : `lsst.afw.image.Exposure` or None
277 The exposure the catalog was built from; used for debug display.
281 struct : `lsst.pipe.base.Struct`
282 The struct contains the following data:
284 - selected : `numpy.ndarray` of `bool``
285 Boolean array of sources that were selected, same length as
291 displayExposure = display
and \
293 plotFwhmHistogram = display
and plt
and \
295 plotFlags = display
and plt
and \
297 plotRejection = display
and plt
and \
299 afwDisplay.setDefaultMaskTransparency(75)
301 fluxName = self.config.fluxName
302 fluxErrName = self.config.fluxErrName
303 minFwhm = self.config.minFwhm
304 maxFwhm = self.config.maxFwhm
305 maxFwhmVariability = self.config.maxFwhmVariability
306 maxellip = self.config.maxellip
307 minsn = self.config.minsn
309 maxelong = (maxellip + 1.0)/(1.0 - maxellip)
if maxellip < 1.0
else 100
312 shape = sourceCat.getShapeDefinition()
313 ixx = sourceCat.get(
"%s.xx" % shape)
314 iyy = sourceCat.get(
"%s.yy" % shape)
316 fwhm = 2*np.sqrt(2*np.log(2))*np.sqrt(0.5*(ixx + iyy))
317 elong = 0.5*(ixx - iyy)/(ixx + iyy)
319 flux = sourceCat.get(fluxName)
320 fluxErr = sourceCat.get(fluxErrName)
321 sn = flux/np.where(fluxErr > 0, fluxErr, 1)
322 sn[fluxErr <= 0] = -psfexLib.BIG
325 for i, f
in enumerate(self.config.badFlags):
326 flags = np.bitwise_or(flags, np.where(sourceCat.get(f), 1 << i, 0))
330 good = np.logical_and(sn > minsn, np.logical_not(flags))
331 good = np.logical_and(good, elong < maxelong)
332 good = np.logical_and(good, fwhm >= minFwhm)
333 good = np.logical_and(good, fwhm < maxFwhm)
335 fwhmMode, fwhmMin, fwhmMax =
compute_fwhmrange(fwhm[good], maxFwhmVariability, minFwhm, maxFwhm,
336 plot=dict(fwhmHistogram=plotFwhmHistogram))
348 selectionVectors = []
349 selectionVectors.append((bad,
"flags %d" % sum(bad)))
353 bad = np.logical_or(bad, dbad)
355 selectionVectors.append((dbad,
"S/N %d" % sum(dbad)))
357 dbad = fwhm < fwhmMin
359 bad = np.logical_or(bad, dbad)
361 selectionVectors.append((dbad,
"fwhmMin %d" % sum(dbad)))
363 dbad = fwhm > fwhmMax
365 bad = np.logical_or(bad, dbad)
367 selectionVectors.append((dbad,
"fwhmMax %d" % sum(dbad)))
369 dbad = elong > maxelong
371 bad = np.logical_or(bad, dbad)
373 selectionVectors.append((dbad,
"elong %d" % sum(dbad)))
375 good = np.logical_not(bad)
381 mi = exposure.getMaskedImage()
382 disp = afwDisplay.Display(frame=frame)
383 disp.mtv(mi, title=
"PSF candidates")
385 with disp.Buffering():
386 for i, source
in enumerate(sourceCat):
388 ctype = afwDisplay.GREEN
390 ctype = afwDisplay.RED
392 disp.dot(
"+", source.getX() - mi.getX0(), source.getY() - mi.getY0(),
395 if plotFlags
or plotRejection:
396 imag = -2.5*np.log10(flux)
401 isSet = np.where(flags == 0x0)[0]
402 plt.plot(imag[isSet], fwhm[isSet],
'o', alpha=alpha, label=
"good")
404 for i, f
in enumerate(self.config.badFlags):
406 isSet = np.where(np.bitwise_and(flags, mask))[0]
408 if np.isfinite(imag[isSet] + fwhm[isSet]).
any():
409 label = re.sub(
r"\_flag",
"",
410 re.sub(
r"^base\_",
"",
411 re.sub(
r"^.*base\_PixelFlags\_flag\_",
"", f)))
412 plt.plot(imag[isSet], fwhm[isSet],
'o', alpha=alpha, label=label)
414 for bad, label
in selectionVectors:
415 plt.plot(imag[bad], fwhm[bad],
'o', alpha=alpha, label=label)
417 plt.plot(imag[good], fwhm[good],
'o', color=
"black", label=
"selected")
418 [plt.axhline(_, color=
'red')
for _
in [fwhmMin, fwhmMax]]
419 plt.xlim(np.median(imag[good]) + 5*np.array([-1, 1]))
420 plt.ylim(fwhm[np.where(np.isfinite(fwhm + imag))].
min(), 2*fwhmMax)
422 plt.xlabel(
"Instrumental %s Magnitude" % fluxName.split(
".")[-1].title())
424 title =
"PSFEX Star Selection"
425 plt.title(
"%s %d selected" % (title, sum(good)))
429 eventHandler = EventHandler(plt.axes(), imag, fwhm, sourceCat.getX(), sourceCat.getY(),
432 if plotFlags
or plotRejection:
435 reply = input(
"continue? [y[es] h(elp) p(db) q(uit)] ").
strip()
444 At this prompt, you can continue with almost any key; 'p' enters pdb,
445 'q' returns to the shell, and
451 If you put the cursor on a point in the matplotlib scatter plot and hit 'p' you'll see it in ds9.""")
452 elif reply[0] ==
"p":
455 elif reply[0] ==
'q':
460 return Struct(selected=good)
bool any(CoordinateExpr< N > const &expr) noexcept
Return true if any elements are true.
def compute_fwhmrange(fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False))