LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Functions | Variables
lsst.meas.extensions.psfex.psfex Namespace Reference

Functions

def compute_fwhmrange (fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False))
 
def select_candidates (set, prefs, frmin, frmax, flags, flux, fluxerr, rmsSize, elong, vignet, plot=dict(), title="")
 
def showPsf (psf, set, ext=None, wcsData=None, trim=0, nspot=5, diagnostics=False, outDir="", frame=None, title=None)
 
def getFlags (tab=None)
 
def guessCalexp (fileName)
 
def makeitLsst (prefs, context, saveWcs=False, plot=dict())
 
def read_samplesLsst (prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, nobj, plot=dict(showFlags=False, showRejection=False))
 
def load_samplesLsst (prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict())
 

Variables

 plt = None
 

Function Documentation

◆ compute_fwhmrange()

def lsst.meas.extensions.psfex.psfex.compute_fwhmrange (   fwhm,
  maxvar,
  minin,
  maxin,
  plot = dict(fwhmHistogram=False) 
)
Compute the FWHM range associated to a series of FWHM measurements.
    AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
    VERSION 20/03/2008

Parameters
----------
fwhm: iterable of `float`
    Iterable of full width half-maxima.
maxvar: `float`
    Maximum allowed FWHM variation.
minin: `float`
    Minimum allowed FWHM.
maxin: `float`
    Maximum allowed FWHM.
plot: `dict`, optional
    Dict of plotting options.

Returns
-------
fmin: `float`
    FWHM mode.
minout: `float`
    Lower FWHM range.
maxout: `float`
    Upper FWHM range.

Definition at line 22 of file psfex.py.

22 def compute_fwhmrange(fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False)):
23  """Compute the FWHM range associated to a series of FWHM measurements.
24  AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
25  VERSION 20/03/2008
26 
27  Parameters
28  ----------
29  fwhm: iterable of `float`
30  Iterable of full width half-maxima.
31  maxvar: `float`
32  Maximum allowed FWHM variation.
33  minin: `float`
34  Minimum allowed FWHM.
35  maxin: `float`
36  Maximum allowed FWHM.
37  plot: `dict`, optional
38  Dict of plotting options.
39 
40  Returns
41  -------
42  fmin: `float`
43  FWHM mode.
44  minout: `float`
45  Lower FWHM range.
46  maxout: `float`
47  Upper FWHM range.
48  """
49 
50  nfwhm = len(fwhm)
51  fwhm.sort()
52 
53  # Find the mode
54  nw = nfwhm//4
55  if nw < 4:
56  nw = 1
57  dfmin = psfexLib.BIG
58  fmin = 0.0
59  for i in range(nfwhm - nw):
60  df = fwhm[i + nw] - fwhm[i]
61  if df < dfmin:
62  dfmin = df
63  fmin = (fwhm[i + nw] + fwhm[i])/2.0
64 
65  if nfwhm < 2:
66  fmin = fwhm[0]
67 
68  dfmin = (maxvar + 1.0)**0.3333333
69  minout = fmin/dfmin if dfmin > 0.0 else 0.0
70  if minout < minin:
71  minout = minin
72 
73  maxout = fmin*dfmin**2
74  if maxout > maxin:
75  maxout = maxin
76 
77  if plt and plot.get("fwhmHistogram"):
78  plt.clf()
79  plt.hist(fwhm, nfwhm//10 + 1, normed=1, facecolor='g', alpha=0.75)
80  plt.xlabel("FWHM")
81  plt.axvline(fmin, color='red')
82  [plt.axvline(_, color='blue') for _ in (minout, maxout)]
83 
84  input("Continue? ")
85 
86  return fmin, minout, maxout
87 
88 
def compute_fwhmrange(fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False))
Definition: psfex.py:22

◆ getFlags()

def lsst.meas.extensions.psfex.psfex.getFlags (   tab = None)

Definition at line 286 of file psfex.py.

286 def getFlags(tab=None):
287  flagKeys = [
288  "base_PixelFlags_flag_edge",
289  # "base_PixelFlags_flag_interpolated",
290  # "base_PixelFlags_flag_interpolatedCenter",
291  # "base_PixelFlags_flag_saturated",
292  "base_PixelFlags_flag_saturatedCenter",
293  # "base_PixelFlags_flag_cr",
294  "base_PixelFlags_flag_crCenter",
295  "base_PixelFlags_flag_bad",
296  "base_PsfFlux_flag",
297  "parent",
298  ]
299 
300  if tab is None:
301  flags = {}
302  for i, k in enumerate(flagKeys):
303  flags[1 << i] = re.sub(r"\_flag", "",
304  re.sub(r"^base\_", "", re.sub(r"^base\_PixelFlags\_flag\_", "", k)))
305  else:
306  flags = 0
307  for i, k in enumerate(flagKeys):
308  if k == "parent":
309  try:
310  isSet = tab.get("deblend_nChild") > 0
311  except KeyError:
312  isSet = 0
313  else:
314  isSet = tab.get(k)
315  flags = np.bitwise_or(flags, np.where(isSet, 1 << i, 0))
316 
317  return flags
318 
319 

◆ guessCalexp()

def lsst.meas.extensions.psfex.psfex.guessCalexp (   fileName)

Definition at line 320 of file psfex.py.

320 def guessCalexp(fileName):
321  for guess in [
322  re.sub("/src", r"", fileName),
323  re.sub("(SRC([^.]+))", r"CORR\2-exp", fileName),
324  ]:
325  if guess != fileName and os.path.exists(guess):
326  return guess
327 
328  raise RuntimeError("Unable to find a calexp to go with %s" % fileName)
329 
330 

◆ load_samplesLsst()

def lsst.meas.extensions.psfex.psfex.load_samplesLsst (   prefs,
  context,
  ext = psfexLib.Prefs.ALL_EXTENSIONS,
  next = 1,
  plot = dict() 
)

Definition at line 522 of file psfex.py.

522 def load_samplesLsst(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict()):
523  minsn = prefs.getMinsn()
524  maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100
525  min = prefs.getFwhmrange()[0]
526  max = prefs.getFwhmrange()[1]
527 
528  filenames = prefs.getCatalogs()
529 
530  ncat = len(filenames)
531  fwhmmin = np.empty(ncat)
532  fwhmmax = np.empty(ncat)
533 
534  if not prefs.getAutoselectFlag():
535  fwhmmin = np.zeros(ncat) + prefs.getFwhmrange()[0]
536  fwhmmax = np.zeros(ncat) + prefs.getFwhmrange()[1]
537  fwhmmode = (fwhmmin + fwhmmax)/2.0
538  else:
539  fwhms = {}
540 
541  # -- Try to estimate the most appropriate Half-light Radius range
542  # -- Get the Half-light radii
543  nobj = 0
544  for i, fileName in enumerate(filenames):
545  fwhms[i] = []
546 
547  if prefs.getVerboseType() != prefs.QUIET:
548  print("Examining Catalog #%d" % (i+1))
549 
550  # ---- Read input catalog
551  tab = afwTable.SourceCatalog.readFits(fileName)
552 
553  # -------- Fill the FWHM array
554  shape = tab.getShapeDefinition()
555  ixx = tab.get("%s.xx" % shape)
556  iyy = tab.get("%s.yy" % shape)
557 
558  rmsSize = np.sqrt(0.5*(ixx + iyy))
559  elong = 0.5*(ixx - iyy)/(ixx + iyy)
560 
561  flux = tab.get(prefs.getPhotfluxRkey())
562  fluxErr = tab.get(prefs.getPhotfluxerrRkey())
563 
564  flags = getFlags(tab)
565 
566  good = np.logical_and(flux/fluxErr > minsn,
567  np.logical_not(flags & prefs.getFlagMask()))
568  good = np.logical_and(good, elong < maxelong)
569  fwhm = 2.0*rmsSize
570  good = np.logical_and(good, fwhm >= min)
571  good = np.logical_and(good, fwhm < max)
572  fwhms[i] = fwhm[good]
573 
574  if prefs.getVarType() == prefs.VAR_NONE:
575  if nobj:
576  fwhms_all = np.empty(sum([len(f) for f in fwhms.values()]))
577  i = 0
578  for f in fwhms.values():
579  fwhms_all[i:len(f)] = f
580  i += len(f)
581  mode, min, max = compute_fwhmrange(fwhms_all, prefs.getMaxvar(),
582  prefs.getFwhmrange()[0], prefs.getFwhmrange()[1],
583  plot=plot)
584  else:
585  raise RuntimeError("No source with appropriate FWHM found!!")
586  mode = min = max = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
587 
588  fwhmmin = np.zeros(ncat) + min
589  fwhmmax = np.zeros(ncat) + max
590  fwhmmode = np.zeros(ncat) + mode
591  else:
592  fwhmmode = np.empty(ncat)
593  fwhmmin = np.empty(ncat)
594  fwhmmax = np.empty(ncat)
595 
596  for i in range(ncat):
597  nobj = len(fwhms[i])
598  if (nobj):
599  fwhmmode[i], fwhmmin[i], fwhmmax[i] = \
600  compute_fwhmrange(fwhms[i], prefs.getMaxvar(),
601  prefs.getFwhmrange()[0], prefs.getFwhmrange()[1], plot=plot)
602  else:
603  raise RuntimeError("No source with appropriate FWHM found!!")
604  fwhmmode[i] = fwhmmin[i] = fwhmmax[i] = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
605 
606  # Read the samples
607  mode = psfexLib.BIG # mode of FWHM distribution
608 
609  sets = []
610  for i, fileName in enumerate(filenames):
611  set = None
612  for ext in range(next):
613  set = read_samplesLsst(prefs, set, fileName, fwhmmin[i]/2.0, fwhmmax[i]/2.0,
614  ext, next, i, context,
615  context.getPc(i) if context.getNpc() else None, nobj, plot=plot)
616 
617  if fwhmmode[i] < mode:
618  mode = fwhmmode[i]
619 
620  set.setFwhm(mode)
621 
622  if prefs.getVerboseType() != prefs.QUIET:
623  if set.getNsample():
624  print("%d samples loaded." % set.getNsample())
625  else:
626  raise RuntimeError("No appropriate source found!!")
627 
628  sets.append(set)
629 
630  return sets
def load_samplesLsst(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict())
Definition: psfex.py:522
def read_samplesLsst(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, nobj, plot=dict(showFlags=False, showRejection=False))
Definition: psfex.py:389

◆ makeitLsst()

def lsst.meas.extensions.psfex.psfex.makeitLsst (   prefs,
  context,
  saveWcs = False,
  plot = dict() 
)
This is the python wrapper that reads lsst tables

Definition at line 331 of file psfex.py.

331 def makeitLsst(prefs, context, saveWcs=False, plot=dict()):
332  """This is the python wrapper that reads lsst tables
333  """
334  # Create an array of PSFs (one PSF for each extension)
335  if prefs.getVerboseType() != prefs.QUIET:
336  print("----- %d input catalogues:" % prefs.getNcat())
337 
338  if saveWcs: # only needed for making plots
339  wcssList = []
340 
341  fields = psfexLib.vectorField()
342  for cat in prefs.getCatalogs():
343  field = psfexLib.Field(cat)
344  wcss = []
345  wcssList.append(wcss)
346  with fits.open(cat):
347  # Hack: I want the WCS so I'll guess where the calexp is to be
348  # found
349  calexpFile = guessCalexp(cat)
350  md = readMetadata(calexpFile)
351  wcs = afwGeom.makeSkyWcs(md)
352 
353  if not wcs:
354  cdMatrix = np.array([1.0, 0.0, 0.0, 1.0])
355  cdMatrix.shape = (2, 2)
356  wcs = afwGeom.makeSkyWcs(crpix=geom.PointD(0, 0),
357  crval=geom.SpherePoint(0.0, 0.0, geom.degrees),
358  cdMatrix=cdMatrix)
359 
360  naxis1, naxis2 = md.getScalar("NAXIS1"), md.getScalar("NAXIS2")
361  # Find how many rows there are in the catalogue
362  md = readMetadata(cat)
363 
364  field.addExt(wcs, naxis1, naxis2, md.getScalar("NAXIS2"))
365  if saveWcs:
366  wcss.append((wcs, naxis1, naxis2))
367 
368  field.finalize()
369  fields.append(field)
370 
371  fields[0].getNext() # number of extensions
372 
373  prefs.getPsfStep()
374 
375  sets = psfexLib.vectorSet()
376  for set in load_samplesLsst(prefs, context, plot=plot):
377  sets.append(set)
378 
379  psfexLib.makeit(fields, sets)
380 
381  ret = [[f.getPsfs() for f in fields], sets]
382  if saveWcs:
383  ret.append(wcssList)
384 
385  return ret
386 
387 
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
std::shared_ptr< daf::base::PropertyList > readMetadata(std::string const &fileName, int hdu=DEFAULT_HDU, bool strip=false)
Read FITS header.
Definition: fits.cc:1657
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition: SkyWcs.cc:521
def makeitLsst(prefs, context, saveWcs=False, plot=dict())
Definition: psfex.py:331

◆ read_samplesLsst()

def lsst.meas.extensions.psfex.psfex.read_samplesLsst (   prefs,
  set,
  filename,
  frmin,
  frmax,
  ext,
  next,
  catindex,
  context,
  pcval,
  nobj,
  plot = dict(showFlags=False, showRejection=False) 
)

Definition at line 388 of file psfex.py.

389  plot=dict(showFlags=False, showRejection=False)):
390  # allocate a new set iff set is None
391  if not set:
392  set = psfexLib.Set(context)
393 
394  cmin, cmax = None, None
395  if set.getNcontext():
396  cmin = np.empty(set.getNcontext())
397  cmax = np.empty(set.getNcontext())
398  for i in range(set.getNcontext()):
399  if set.getNsample():
400  cmin[i] = set.getContextOffset(i) - set.getContextScale(i)/2.0
401  cmax[i] = cmin[i] + set.getContextScale(i)
402  else:
403  cmin[i] = psfexLib.BIG
404  cmax[i] = -psfexLib.BIG
405  #
406  # Read data
407  #
408  tab = afwTable.SourceCatalog.readFits(filename)
409 
410  centroid = tab.getCentroidDefinition()
411  xm = tab.get("%s.x" % centroid)
412  ym = tab.get("%s.y" % centroid)
413 
414  shape = tab.getShapeDefinition()
415  ixx = tab.get("%s.xx" % shape)
416  iyy = tab.get("%s.yy" % shape)
417 
418  rmsSize = np.sqrt(0.5*(ixx + iyy))
419  elong = 0.5*(ixx - iyy)/(ixx + iyy)
420 
421  flux = tab.get(prefs.getPhotfluxRkey())
422  fluxErr = tab.get(prefs.getPhotfluxerrRkey())
423  flags = getFlags(tab)
424 
425  #
426  # Now the VIGNET data
427  #
428  vigw, vigh = 35, 35 # [prefs.getPsfsize()[i] for i in range(2)]
429  if set.empty():
430  set.setVigSize(vigw, vigh)
431 
432  vignet = np.empty(nobj*vigw*vigh, "float32").reshape(nobj, vigw, vigh)
433 
434  # Hack: I want the WCS so I'll guess where the calexp is to be found
435  calexpFile = guessCalexp(filename)
436  mi = afwImage.MaskedImageF(calexpFile)
437  backnoise2 = np.median(mi.getVariance().getArray())
438  gain = 1.0
439 
440  edgeBit = [k for k, v in getFlags().items() if v == "edge"][0]
441 
442  for i, xc, yc in zip(range(nobj), xm, ym):
443  try:
444  x, y = int(xc), int(yc)
445  except ValueError:
446  flags[i] |= edgeBit # mark star as bad
447 
448  try:
449  pstamp = mi[x - vigw//2:x + vigw//2 + 1, y - vigh//2:y + vigh//2 + 1]
450  vignet[i] = pstamp.getImage().getArray().transpose()
451  except Exception:
452  flags[i] |= edgeBit # mark star as bad
453 
454  # Try to load the set of context keys
455  pc = 0
456  contextvalp = []
457  for i, key in enumerate(context.getName()):
458  if context.getPcflag(i):
459  contextvalp.append(pcval[pc])
460  pc += 1
461  elif key[0] == ':':
462  try:
463  contextvalp.append(tab.header[key[1:]])
464  except KeyError:
465  raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
466  (key[1:], filename))
467  else:
468  try:
469  contextvalp.append(tab.get(key))
470  except KeyError:
471  raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
472  (key, filename))
473  set.setContextname(i, key)
474 
475  # Now examine each vector of the shipment
476  good = select_candidates(set, prefs, frmin, frmax,
477  flags, flux, fluxErr, rmsSize, elong, vignet,
478  plot=plot, title="%s[%d]" % (filename, ext + 1) if next > 1 else filename)
479  #
480  # Insert the good candidates into the set
481  #
482  if not vignet.dtype.isnative:
483  # without the swap setVig fails with
484  # "ValueError: 'unaligned arrays cannot be converted to C++'"
485  vignet = vignet.byteswap()
486 
487  for i in np.where(good)[0]:
488  sample = set.newSample()
489  sample.setCatindex(catindex)
490  sample.setExtindex(ext)
491 
492  sample.setVig(vignet[i])
493  sample.setNorm(float(flux[i]))
494  sample.setBacknoise2(backnoise2)
495  sample.setGain(gain)
496  sample.setX(float(xm[i]))
497  sample.setY(float(ym[i]))
498  sample.setFluxrad(float(rmsSize[i]))
499 
500  for j in range(set.getNcontext()):
501  sample.setContext(j, float(contextvalp[j][i]))
502 
503  set.finiSample(sample, prefs.getProfAccuracy())
504 
505  # ---- Update min and max
506  for j in range(set.getNcontext()):
507  cmin[j] = contextvalp[j][good].min()
508  cmax[j] = contextvalp[j][good].max()
509 
510  # Update the scaling
511  if set.getNsample():
512  for i in range(set.getNcontext()):
513  set.setContextScale(i, cmax[i] - cmin[i])
514  set.setContextOffset(i, (cmin[i] + cmax[i])/2.0)
515 
516  # Don't waste memory!
517  set.trimMemory()
518 
519  return set
520 
521 
std::vector< SchemaItem< Flag > > * items
int min
int max
def select_candidates(set, prefs, frmin, frmax, flags, flux, fluxerr, rmsSize, elong, vignet, plot=dict(), title="")
Definition: psfex.py:91

◆ select_candidates()

def lsst.meas.extensions.psfex.psfex.select_candidates (   set,
  prefs,
  frmin,
  frmax,
  flags,
  flux,
  fluxerr,
  rmsSize,
  elong,
  vignet,
  plot = dict(),
  title = "" 
)

Definition at line 89 of file psfex.py.

91  plot=dict(), title=""):
92  maxbad = prefs.getBadpixNmax()
93  maxbadflag = prefs.getBadpixFlag()
94  maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100.0
95  minsn = prefs.getMinsn()
96 
97  sn = flux/np.where(fluxerr > 0, fluxerr, 1)
98  sn[fluxerr <= 0] = -psfexLib.BIG
99  # ---- Apply some selection over flags, fluxes...
100  plotFlags = plot.get("showFlags") if plt else False
101  plotRejection = plot.get("showRejection") if plt else False
102 
103  bad = flags & prefs.getFlagMask() != 0
104  set.setBadFlags(int(sum(bad != 0)))
105 
106  if plotRejection:
107  selectionVectors = []
108  selectionVectors.append((bad, "flags %d" % sum(bad != 0)))
109 
110  dbad = sn < minsn
111  set.setBadSN(int(sum(dbad)))
112  bad = np.logical_or(bad, dbad)
113  if plotRejection:
114  selectionVectors.append((dbad, "S/N %d" % sum(dbad)))
115 
116  dbad = rmsSize < frmin
117  set.setBadFrmin(int(sum(dbad)))
118  bad = np.logical_or(bad, dbad)
119  if plotRejection:
120  selectionVectors.append((dbad, "frmin %d" % sum(dbad)))
121 
122  dbad = rmsSize > frmax
123  set.setBadFrmax(int(sum(dbad)))
124  bad = np.logical_or(bad, dbad)
125  if plotRejection:
126  selectionVectors.append((dbad, "frmax %d" % sum(dbad)))
127 
128  dbad = elong > maxelong
129  set.setBadElong(int(sum(dbad)))
130  bad = np.logical_or(bad, dbad)
131  if plotRejection:
132  selectionVectors.append((dbad, "elong %d" % sum(dbad)))
133 
134  # -- ... and check the integrity of the sample
135  if maxbadflag:
136  nbad = np.array([(v <= -psfexLib.BIG).sum() for v in vignet])
137  dbad = nbad > maxbad
138  set.setBadPix(int(sum(dbad)))
139  bad = np.logical_or(bad, dbad)
140  if plotRejection:
141  selectionVectors.append((dbad, "badpix %d" % sum(dbad)))
142 
143  good = np.logical_not(bad)
144  if plotFlags or plotRejection:
145  imag = -2.5*np.log10(flux)
146  plt.clf()
147 
148  alpha = 0.5
149  if plotFlags:
150  labels = getFlags()
151 
152  isSet = np.where(flags == 0x0)[0]
153  plt.plot(imag[isSet], rmsSize[isSet], 'o', alpha=alpha, label="good")
154 
155  for i in range(16):
156  mask = 1 << i
157  if mask & prefs.getFlagMask():
158  isSet = np.where(np.bitwise_and(flags, mask))[0]
159  if isSet.any():
160  plt.plot(imag[isSet], rmsSize[isSet], 'o', alpha=alpha, label=labels[mask])
161  else:
162  for bad, label in selectionVectors:
163  plt.plot(imag[bad], rmsSize[bad], 'o', alpha=alpha, label=label)
164 
165  plt.plot(imag[good], rmsSize[good], 'o', color="black", label="selected")
166  [plt.axhline(_, color='red') for _ in [frmin, frmax]]
167  plt.xlim(np.median(imag[good]) + 5*np.array([-1, 1]))
168  plt.ylim(-0.1, 2*frmax)
169  plt.legend(loc=2)
170  plt.xlabel("Instrumental Magnitude")
171  plt.ylabel("rmsSize")
172  plt.title("%s %d selected" % (title, sum(good)))
173 
174  input("Continue? ")
175 
176  return good
177 
178 

◆ showPsf()

def lsst.meas.extensions.psfex.psfex.showPsf (   psf,
  set,
  ext = None,
  wcsData = None,
  trim = 0,
  nspot = 5,
  diagnostics = False,
  outDir = "",
  frame = None,
  title = None 
)
Show a PSF on display (e.g., ds9)

Definition at line 179 of file psfex.py.

180  diagnostics=False, outDir="", frame=None, title=None):
181  """Show a PSF on display (e.g., ds9)
182  """
183 
184  if ext is not None:
185  psf = psf[ext]
186 
187  if wcsData:
188  if ext is not None:
189  wcsData = wcsData[ext]
190  wcs, naxis1, naxis2 = wcsData
191  else:
192  wcs, naxis1, naxis2 = None, None, None
193 
194  naxis = [naxis1, naxis2]
195  for i in range(2):
196  if naxis[i] is None:
197  # cmin, cmax are the range of input star positions
198  cmin, cmax = [set.getContextOffset(i) + d*set.getContextScale(i) for d in (-0.5, 0.5)]
199  naxis[i] = cmax + cmin # a decent guess
200 
201  if naxis[0] > naxis[1]:
202  nx, ny = int(nspot*naxis[0]/float(naxis[1]) + 0.5), nspot
203  else:
204  nx, ny = nspot, int(nspot*naxis[1]/float(naxis[0]) + 0.5)
205 
206  mos = afwDisplay.utils.Mosaic(gutter=2, background=0.02)
207 
208  xpos, ypos = np.linspace(0, naxis[0], nx), np.linspace(0, naxis[1], ny)
209  for y in ypos:
210  for x in xpos:
211  psf.build(x, y)
212 
213  im = afwImage.ImageF(*psf.getLoc().shape)
214  im.getArray()[:] = psf.getLoc()
215  im /= float(im.getArray().max())
216  if trim:
217  if trim > im.getHeight()//2:
218  trim = im.getHeight()//2
219 
220  im = im[trim:-trim, trim:-trim]
221 
222  mos.append(im)
223 
224  mosaic = mos.makeMosaic(mode=nx)
225  if frame is not None:
226  afwDisplay.Display(frame=frame).mtv(mosaic, title=title)
227  #
228  # Figure out the WCS for the mosaic
229  #
230  pos = []
231  pos.append([geom.PointD(0, 0), wcs.pixelToSky(geom.PointD(0, 0))])
232  pos.append([geom.PointD(*mosaic.getDimensions()), wcs.pixelToSky(geom.PointD(naxis1, naxis2))])
233 
234  CD = []
235  for i in range(2):
236  delta = pos[1][1][i].asDegrees() - pos[0][1][i].asDegrees()
237  CD.append(delta/(pos[1][0][i] - pos[0][0][i]))
238  CD = np.array(CD)
239  CD.shape = (2, 2)
240  mosWcs = afwGeom.makeSkyWcs(crval=pos[0][0], crpix=pos[0][1], cdMatrix=CD)
241 
242  if ext is not None:
243  title = "%s-%d" % (title, ext)
244 
245  if frame is not None:
246  afwDisplay.Display(frame=frame).mtv(mosaic, title=title, wcs=mosWcs)
247 
248  if diagnostics:
249  outFile = "%s-mod.fits" % title
250  if outDir:
251  outFile = os.path.join(outDir, outFile)
252  mosaic.writeFits(outFile, mosWcs.getFitsMetadata())
253 
254  mos = afwDisplay.utils.Mosaic(gutter=4, background=0.002)
255  for i in range(set.getNsample()):
256  s = set.getSample(i)
257  if ext is not None and s.getExtindex() != ext:
258  continue
259 
260  smos = afwDisplay.utils.Mosaic(gutter=2, background=-0.003)
261  for func in [s.getVig, s.getVigResi]:
262  arr = func()
263  if func == s.getVig:
264  norm = float(arr.max()) if True else s.getNorm()
265 
266  arr /= norm
267  im = afwImage.ImageF(*arr.shape)
268  im.getArray()[:] = arr
269  smos.append(im)
270 
271  mos.append(smos.makeMosaic(mode="x"))
272 
273  mosaic = mos.makeMosaic(title=title)
274 
275  if frame is not None:
276  afwDisplay.Display(frame=frame + 1).mtv(mosaic, title=title)
277 
278  if diagnostics:
279  outFile = "%s-psfstars.fits" % title
280  if outDir:
281  outFile = os.path.join(outDir, outFile)
282 
283  mosaic.writeFits(outFile)
284 
285 
def mtv(data, frame=None, title="", wcs=None, *args, **kwargs)
Definition: ds9.py:92

Variable Documentation

◆ plt

lsst.meas.extensions.psfex.psfex.plt = None

Definition at line 9 of file psfex.py.