LSSTApplications  20.0.0
LSSTDataManagementBasePackage
Classes | Functions | Variables
lsst.meas.extensions.psfex.psfex Namespace Reference

Classes

class  _LSST
 
class  _SExtractor
 

Functions

def splitFitsCard (line)
 
def compute_fwhmrange (fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False))
 
def read_samples (prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, plot=dict(showFlags=False, showRejection=False))
 
def getSexFlags (*args)
 
def select_candidates (set, prefs, frmin, frmax, flags, flux, fluxerr, rmsSize, elong, vignet, plot=dict(), title="")
 
def setDataType (t)
 
def getFlags ()
 
def load_samples (prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict())
 
def showPsf (psf, set, ext=None, wcsData=None, trim=0, nspot=5, diagnostics=False, outDir="", frame=None, title=None)
 
def getLsstFlags (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())
 
def makeit (prefs, context, saveWcs=False, plot=dict())
 

Variables

 plt = None
 
 LSST
 
 SExtractor
 

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 56 of file psfex.py.

56 def compute_fwhmrange(fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False)):
57  """Compute the FWHM range associated to a series of FWHM measurements.
58  AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
59  VERSION 20/03/2008
60 
61  Parameters
62  ----------
63  fwhm: iterable of `float`
64  Iterable of full width half-maxima.
65  maxvar: `float`
66  Maximum allowed FWHM variation.
67  minin: `float`
68  Minimum allowed FWHM.
69  maxin: `float`
70  Maximum allowed FWHM.
71  plot: `dict`, optional
72  Dict of plotting options.
73 
74  Returns
75  -------
76  fmin: `float`
77  FWHM mode.
78  minout: `float`
79  Lower FWHM range.
80  maxout: `float`
81  Upper FWHM range.
82  """
83 
84  nfwhm = len(fwhm)
85  fwhm.sort()
86 
87  # Find the mode
88  nw = nfwhm//4
89  if nw < 4:
90  nw = 1
91  dfmin = psfexLib.BIG
92  fmin = 0.0
93  for i in range(nfwhm - nw):
94  df = fwhm[i + nw] - fwhm[i]
95  if df < dfmin:
96  dfmin = df
97  fmin = (fwhm[i + nw] + fwhm[i])/2.0
98 
99  if nfwhm < 2:
100  fmin = fwhm[0]
101 
102  dfmin = (maxvar + 1.0)**0.3333333
103  minout = fmin/dfmin if dfmin > 0.0 else 0.0
104  if minout < minin:
105  minout = minin
106 
107  maxout = fmin*dfmin**2
108  if maxout > maxin:
109  maxout = maxin
110 
111  if plt and plot.get("fwhmHistogram"):
112  plt.clf()
113  plt.hist(fwhm, nfwhm//10 + 1, normed=1, facecolor='g', alpha=0.75)
114  plt.xlabel("FWHM")
115  plt.axvline(fmin, color='red')
116  [plt.axvline(_, color='blue') for _ in (minout, maxout)]
117 
118  input("Continue? ")
119 
120  return fmin, minout, maxout
121 
122 

◆ getFlags()

def lsst.meas.extensions.psfex.psfex.getFlags ( )

Definition at line 404 of file psfex.py.

404 def getFlags():
405  if _dataType == _LSST:
406  return getLsstFlags(None)
407  else:
408  return getSexFlags(None)
409 
410 

◆ getLsstFlags()

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

Definition at line 661 of file psfex.py.

661 def getLsstFlags(tab=None):
662  flagKeys = [
663  "base_PixelFlags_flag_edge",
664  # "base_PixelFlags_flag_interpolated",
665  # "base_PixelFlags_flag_interpolatedCenter",
666  # "base_PixelFlags_flag_saturated",
667  "base_PixelFlags_flag_saturatedCenter",
668  # "base_PixelFlags_flag_cr",
669  "base_PixelFlags_flag_crCenter",
670  "base_PixelFlags_flag_bad",
671  "base_PsfFlux_flag",
672  "parent",
673  ]
674 
675  if tab is None:
676  flags = {}
677  for i, k in enumerate(flagKeys):
678  flags[1 << i] = re.sub(r"\_flag", "",
679  re.sub(r"^base\_", "", re.sub(r"^base\_PixelFlags\_flag\_", "", k)))
680  else:
681  flags = 0
682  for i, k in enumerate(flagKeys):
683  if k == "parent":
684  try:
685  isSet = tab.get("deblend_nChild") > 0
686  except KeyError:
687  isSet = 0
688  else:
689  isSet = tab.get(k)
690  flags = np.bitwise_or(flags, np.where(isSet, 1 << i, 0))
691 
692  return flags
693 
694 

◆ getSexFlags()

def lsst.meas.extensions.psfex.psfex.getSexFlags ( *  args)

Definition at line 283 of file psfex.py.

283 def getSexFlags(*args):
284  return {1: "flux blended",
285  2: "blended",
286  4: "saturated",
287  8: "edge",
288  16: "bad aperture",
289  32: "bad isophotal",
290  64: "memory error (deblend)",
291  128: "memory error (extract)",
292  }
293 
294 

◆ guessCalexp()

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

Definition at line 695 of file psfex.py.

695 def guessCalexp(fileName):
696  for guess in [
697  re.sub("/src", r"", fileName),
698  re.sub("(SRC([^.]+))", r"CORR\2-exp", fileName),
699  ]:
700  if guess != fileName and os.path.exists(guess):
701  return guess
702 
703  raise RuntimeError("Unable to find a calexp to go with %s" % fileName)
704 
705 

◆ load_samples()

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

Definition at line 411 of file psfex.py.

411 def load_samples(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict()):
412  minsn = prefs.getMinsn()
413  maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100
414  min = prefs.getFwhmrange()[0]
415  max = prefs.getFwhmrange()[1]
416 
417  filenames = prefs.getCatalogs()
418 
419  ncat = len(filenames)
420  fwhmmin = np.empty(ncat)
421  fwhmmax = np.empty(ncat)
422 
423  if not prefs.getAutoselectFlag():
424  fwhmmin = np.zeros(ncat) + prefs.getFwhmrange()[0]
425  fwhmmax = np.zeros(ncat) + prefs.getFwhmrange()[1]
426  fwhmmode = (fwhmmin + fwhmmax)/2.0
427  else:
428  fwhms = {}
429 
430  # -- Try to estimate the most appropriate Half-light Radius range
431  # -- Get the Half-light radii
432  nobj = 0
433  for i, fileName in enumerate(filenames):
434  fwhms[i] = []
435 
436  if prefs.getVerboseType() != prefs.QUIET:
437  print("Examining Catalog #%d" % (i+1))
438 
439  # ---- Read input catalog
440  backnoises = []
441  with fits.open(fileName) as cat:
442  extCtr = -1
443  for tab in cat:
444  if tab.name == "LDAC_IMHEAD":
445  extCtr += 1
446 
447  if extCtr != ext and ext != prefs.ALL_EXTENSIONS:
448  if extCtr > ext:
449  break
450  continue
451 
452  if tab.name == "PRIMARY":
453  pass
454  elif tab.name == "LDAC_IMHEAD":
455  hdr = tab.data[0][0] # the fits header from the original fits image
456  for line in hdr:
457  try:
458  k, v = splitFitsCard(line)
459  except AttributeError:
460  continue
461 
462  if k == "SEXBKDEV":
463  if v < 1/psfexLib.BIG:
464  v = 1.0
465 
466  backnoises.append(v)
467  break
468  elif tab.name == "LDAC_OBJECTS":
469  # -------- Fill the FWHM array
470  rmsSize = tab.data["FLUX_RADIUS"]
471  fmax = tab.data["FLUX_MAX"]
472  flags = tab.data["FLAGS"]
473  elong = tab.data["ELONGATION"]
474  backnoise = backnoises[-1]
475 
476  good = np.logical_and(fmax/backnoise > minsn,
477  np.logical_not(flags & prefs.getFlagMask()))
478  good = np.logical_and(good, elong < maxelong)
479  fwhm = 2.0*rmsSize
480  good = np.logical_and(good, fwhm >= min)
481  good = np.logical_and(good, fwhm < max)
482  fwhms[i] = fwhm[good]
483 
484  if prefs.getVarType() == prefs.VAR_NONE:
485  if nobj:
486  fwhms_all = np.empty(sum([len(l) for l in fwhms.values()]))
487  i = 0
488  for l in fwhms.values():
489  fwhms_all[i:len(l)] = l
490  i += len(l)
491  mode, min, max = compute_fwhmrange(fwhms_all, prefs.getMaxvar(),
492  prefs.getFwhmrange()[0], prefs.getFwhmrange()[1],
493  plot=plot)
494  else:
495  raise RuntimeError("No source with appropriate FWHM found!!")
496  mode = min = max = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
497 
498  fwhmmin = np.zeros(ncat) + min
499  fwhmmax = np.zeros(ncat) + max
500  fwhmmode = np.zeros(ncat) + mode
501  else:
502  fwhmmode = np.empty(ncat)
503  fwhmmin = np.empty(ncat)
504  fwhmmax = np.empty(ncat)
505 
506  for i in range(ncat):
507  nobj = len(fwhms[i])
508  if (nobj):
509  fwhmmode[i], fwhmmin[i], fwhmmax[i] = \
510  compute_fwhmrange(fwhms[i], prefs.getMaxvar(),
511  prefs.getFwhmrange()[0], prefs.getFwhmrange()[1], plot=plot)
512  else:
513  raise RuntimeError("No source with appropriate FWHM found!!")
514  fwhmmode[i] = fwhmmin[i] = fwhmmax[i] = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
515 
516  # Read the samples
517  mode = psfexLib.BIG # mode of FWHM distribution
518 
519  sets = []
520  for i, fileName in enumerate(filenames):
521  set = None
522  if ext == prefs.ALL_EXTENSIONS:
523  extensions = list(range(len(backnoises)))
524  else:
525  extensions = [ext]
526 
527  for e in extensions:
528  set = read_samples(prefs, set, fileName, fwhmmin[i]/2.0, fwhmmax[i]/2.0,
529  e, next, i, context, context.getPc(i) if context.getNpc() else None, plot=plot)
530 
531  if fwhmmode[i] < mode:
532  mode = fwhmmode[i]
533 
534  set.setFwhm(mode)
535 
536  if prefs.getVerboseType() != prefs.QUIET:
537  if set.getNsample():
538  print("%d samples loaded." % set.getNsample())
539  else:
540  raise RuntimeError("No appropriate source found!!")
541 
542  sets.append(set)
543 
544  return sets
545 
546 

◆ load_samplesLsst()

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

Definition at line 902 of file psfex.py.

902 def load_samplesLsst(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict()):
903  minsn = prefs.getMinsn()
904  maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100
905  min = prefs.getFwhmrange()[0]
906  max = prefs.getFwhmrange()[1]
907 
908  filenames = prefs.getCatalogs()
909 
910  ncat = len(filenames)
911  fwhmmin = np.empty(ncat)
912  fwhmmax = np.empty(ncat)
913 
914  if not prefs.getAutoselectFlag():
915  fwhmmin = np.zeros(ncat) + prefs.getFwhmrange()[0]
916  fwhmmax = np.zeros(ncat) + prefs.getFwhmrange()[1]
917  fwhmmode = (fwhmmin + fwhmmax)/2.0
918  else:
919  fwhms = {}
920 
921  # -- Try to estimate the most appropriate Half-light Radius range
922  # -- Get the Half-light radii
923  nobj = 0
924  for i, fileName in enumerate(filenames):
925  fwhms[i] = []
926 
927  if prefs.getVerboseType() != prefs.QUIET:
928  print("Examining Catalog #%d" % (i+1))
929 
930  # ---- Read input catalog
931  tab = afwTable.SourceCatalog.readFits(fileName)
932 
933  # -------- Fill the FWHM array
934  shape = tab.getShapeDefinition()
935  ixx = tab.get("%s.xx" % shape)
936  iyy = tab.get("%s.yy" % shape)
937 
938  rmsSize = np.sqrt(0.5*(ixx + iyy))
939  elong = 0.5*(ixx - iyy)/(ixx + iyy)
940 
941  flux = tab.get(prefs.getPhotfluxRkey())
942  fluxErr = tab.get(prefs.getPhotfluxerrRkey())
943 
944  flags = getLsstFlags(tab)
945 
946  good = np.logical_and(flux/fluxErr > minsn,
947  np.logical_not(flags & prefs.getFlagMask()))
948  good = np.logical_and(good, elong < maxelong)
949  fwhm = 2.0*rmsSize
950  good = np.logical_and(good, fwhm >= min)
951  good = np.logical_and(good, fwhm < max)
952  fwhms[i] = fwhm[good]
953 
954  if prefs.getVarType() == prefs.VAR_NONE:
955  if nobj:
956  fwhms_all = np.empty(sum([len(l) for l in fwhms.values()]))
957  i = 0
958  for l in fwhms.values():
959  fwhms_all[i:len(l)] = l
960  i += len(l)
961  mode, min, max = compute_fwhmrange(fwhms_all, prefs.getMaxvar(),
962  prefs.getFwhmrange()[0], prefs.getFwhmrange()[1],
963  plot=plot)
964  else:
965  raise RuntimeError("No source with appropriate FWHM found!!")
966  mode = min = max = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
967 
968  fwhmmin = np.zeros(ncat) + min
969  fwhmmax = np.zeros(ncat) + max
970  fwhmmode = np.zeros(ncat) + mode
971  else:
972  fwhmmode = np.empty(ncat)
973  fwhmmin = np.empty(ncat)
974  fwhmmax = np.empty(ncat)
975 
976  for i in range(ncat):
977  nobj = len(fwhms[i])
978  if (nobj):
979  fwhmmode[i], fwhmmin[i], fwhmmax[i] = \
980  compute_fwhmrange(fwhms[i], prefs.getMaxvar(),
981  prefs.getFwhmrange()[0], prefs.getFwhmrange()[1], plot=plot)
982  else:
983  raise RuntimeError("No source with appropriate FWHM found!!")
984  fwhmmode[i] = fwhmmin[i] = fwhmmax[i] = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
985 
986  # Read the samples
987  mode = psfexLib.BIG # mode of FWHM distribution
988 
989  sets = []
990  for i, fileName in enumerate(filenames):
991  set = None
992  for ext in range(next):
993  set = read_samplesLsst(prefs, set, fileName, fwhmmin[i]/2.0, fwhmmax[i]/2.0,
994  ext, next, i, context,
995  context.getPc(i) if context.getNpc() else None, nobj, plot=plot)
996 
997  if fwhmmode[i] < mode:
998  mode = fwhmmode[i]
999 
1000  set.setFwhm(mode)
1001 
1002  if prefs.getVerboseType() != prefs.QUIET:
1003  if set.getNsample():
1004  print("%d samples loaded." % set.getNsample())
1005  else:
1006  raise RuntimeError("No appropriate source found!!")
1007 
1008  sets.append(set)
1009 
1010  return sets
1011 
1012 

◆ makeit()

def lsst.meas.extensions.psfex.psfex.makeit (   prefs,
  context,
  saveWcs = False,
  plot = dict() 
) 
This is the python wrapper for the original psfex that reads SExtractor
outputs

Definition at line 1013 of file psfex.py.

1013 def makeit(prefs, context, saveWcs=False, plot=dict()):
1014  """This is the python wrapper for the original psfex that reads SExtractor
1015  outputs"""
1016  # Create an array of PSFs (one PSF for each extension)
1017  if prefs.getVerboseType() != prefs.QUIET:
1018  print("----- %d input catalogues:" % prefs.getNcat())
1019 
1020  if saveWcs: # only needed for making plots
1021  wcssList = []
1022 
1023  fields = psfexLib.vectorField()
1024  for cat in prefs.getCatalogs():
1025  field = psfexLib.Field(cat)
1026  wcss = []
1027  wcssList.append(wcss)
1028  with fits.open(cat) as pf:
1029  for hdu in pf:
1030  if hdu.name == "PRIMARY":
1031  pass
1032  elif hdu.name == "LDAC_IMHEAD":
1033  hdr = hdu.data[0][0] # the fits header from the original fits image
1034  md = PropertySet()
1035  for line in hdr:
1036  try:
1037  md.set(*splitFitsCard(line))
1038  except AttributeError:
1039  continue
1040 
1041  if not md.exists("CRPIX1"): # no WCS; try WCSA
1042  for k in md.names():
1043  if re.search(r"A$", k):
1044  md.set(k[:-1], md.getScalar(k))
1045  wcs = afwGeom.makeSkyWcs(md)
1046  naxis1, naxis2 = md.getScalar("NAXIS1"), md.getScalar("NAXIS2")
1047  elif hdu.name == "LDAC_OBJECTS":
1048  nobj = len(hdu.data)
1049 
1050  assert wcs, "LDAC_OBJECTS comes after LDAC_IMHEAD"
1051  field.addExt(wcs, naxis1, naxis2, nobj)
1052  if saveWcs:
1053  wcss.append((wcs, naxis1, naxis2))
1054  wcs = None
1055 
1056  field.finalize()
1057  fields.append(field)
1058 
1059  sets = psfexLib.vectorSet()
1060  for set in load_samples(prefs, context, plot=plot):
1061  sets.append(set)
1062 
1063  psfexLib.makeit(fields, sets)
1064 
1065  ret = [[f.getPsfs() for f in fields], sets]
1066  if saveWcs:
1067  ret.append(wcssList)
1068 
1069  return ret

◆ 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 706 of file psfex.py.

706 def makeitLsst(prefs, context, saveWcs=False, plot=dict()):
707  """This is the python wrapper that reads lsst tables
708  """
709  # Create an array of PSFs (one PSF for each extension)
710  if prefs.getVerboseType() != prefs.QUIET:
711  print("----- %d input catalogues:" % prefs.getNcat())
712 
713  if saveWcs: # only needed for making plots
714  wcssList = []
715 
716  fields = psfexLib.vectorField()
717  for cat in prefs.getCatalogs():
718  field = psfexLib.Field(cat)
719  wcss = []
720  wcssList.append(wcss)
721  with fits.open(cat):
722  # Hack: I want the WCS so I'll guess where the calexp is to be
723  # found
724  calexpFile = guessCalexp(cat)
725  md = readMetadata(calexpFile)
726  wcs = afwGeom.makeSkyWcs(md)
727 
728  if not wcs:
729  cdMatrix = np.array([1.0, 0.0, 0.0, 1.0])
730  cdMatrix.shape = (2, 2)
731  wcs = afwGeom.makeSkyWcs(crpix=geom.PointD(0, 0),
732  crval=geom.SpherePoint(0.0, 0.0, geom.degrees),
733  cdMatrix=cdMatrix)
734 
735  naxis1, naxis2 = md.getScalar("NAXIS1"), md.getScalar("NAXIS2")
736  # Find how many rows there are in the catalogue
737  md = readMetadata(cat)
738 
739  field.addExt(wcs, naxis1, naxis2, md.getScalar("NAXIS2"))
740  if saveWcs:
741  wcss.append((wcs, naxis1, naxis2))
742 
743  field.finalize()
744  fields.append(field)
745 
746  fields[0].getNext() # number of extensions
747 
748  prefs.getPsfStep()
749 
750  sets = psfexLib.vectorSet()
751  for set in load_samplesLsst(prefs, context, plot=plot):
752  sets.append(set)
753 
754  psfexLib.makeit(fields, sets)
755 
756  ret = [[f.getPsfs() for f in fields], sets]
757  if saveWcs:
758  ret.append(wcssList)
759 
760  return ret
761 
762 

◆ read_samples()

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

Definition at line 123 of file psfex.py.

123 def read_samples(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval,
124  plot=dict(showFlags=False, showRejection=False)):
125  # allocate a new set iff set is None
126  if not set:
127  set = psfexLib.Set(context)
128 
129  cmin, cmax = None, None
130  if set.getNcontext():
131  cmin = np.empty(set.getNcontext())
132  cmax = np.empty(set.getNcontext())
133  for i in range(set.getNcontext()):
134  if set.getNsample():
135  cmin[i] = set.getContextOffset(i) - set.getContextScale(i)/2.0
136  cmax[i] = cmin[i] + set.getContextScale(i)
137  else:
138  cmin[i] = psfexLib.BIG
139  cmax[i] = -psfexLib.BIG
140  #
141  # Read data
142  #
143  with fits.open(filename) as cat:
144  extCtr = -1
145  for tab in cat:
146  if tab.name == "LDAC_IMHEAD":
147  extCtr += 1
148 
149  if extCtr < ext:
150  continue
151  elif extCtr > ext:
152  break
153 
154  if tab.name == "PRIMARY":
155  pass
156  elif tab.name == "LDAC_IMHEAD":
157  hdr = tab.data[0][0] # the fits header from the original fits image
158  foundCards = 0 # how many of our desired cards have we found?
159 
160  for line in hdr:
161  try:
162  k, v = splitFitsCard(line)
163  except AttributeError:
164  continue
165 
166  if k == "SEXBKDEV":
167  backnoise2 = v**2
168  foundCards += 1
169  elif k == "SEXGAIN":
170  gain = v
171  foundCards += 1
172 
173  if foundCards == 2:
174  break
175  elif tab.name == "LDAC_OBJECTS":
176  xm = tab.data[prefs.getCenterKey(0)]
177  ym = tab.data[prefs.getCenterKey(1)]
178  fluxrad = tab.data["FLUX_RADIUS"]
179  flux = tab.data[prefs.getPhotfluxRkey()]
180  fluxerr = tab.data[prefs.getPhotfluxerrRkey()]
181  elong = tab.data["ELONGATION"]
182  flags = tab.data["FLAGS"]
183 
184  n = prefs.getPhotfluxNum() - 1
185  if n:
186  raise RuntimeError("Code to handle e.g. FLUX_APER(3) isn't yet converted")
187  if key.naxis == 1 and n < key.naxisn[0]: # noqa: F821
188  flux += n
189  else:
190  print("Not enough apertures for %s in catalogue %s: using first aperture" %
191  (prefs.getPhotfluxRkey(), filename), file=sys.stderr)
192 
193  n = prefs.getPhotfluxerrNum() - 1
194  if n:
195  raise RuntimeError("Code for getPhotfluxerrNum is broken")
196  if key.naxis == 1 and n < key.naxisn[0]: # noqa: F821
197  fluxerr += n
198  else:
199  print("Not enough apertures for %s in catalogue %s: using first aperture" %
200  (prefs.getPhotfluxerrRkey(), filename), file=sys.stderr)
201  #
202  # Now the VIGNET data
203  #
204  vignet = tab.data["VIGNET"]
205 
206  try:
207  vigw, vigh = vignet[0].shape
208  except ValueError:
209  raise RuntimeError("*Error*: VIGNET should be a 2D vector; saw %s" % str(vignet[0].shape))
210 
211  if set.empty():
212  set.setVigSize(vigw, vigh)
213 
214  # Try to load the set of context keys
215  pc = 0
216  contextvalp = []
217  for i, key in enumerate(context.getName()):
218  if context.getPcflag(i):
219  contextvalp.append(pcval[pc])
220  pc += 1
221  elif key[0] == ':':
222  try:
223  contextvalp.append(tab.header[key[1:]])
224  except KeyError:
225  raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
226  (key[1:], filename))
227  else:
228  try:
229  contextvalp.append(tab.data[key])
230  except KeyError:
231  raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
232  (key, filename))
233  set.setContextname(i, key)
234 
235  # Now examine each vector of the shipment
236  good = select_candidates(set, prefs, frmin, frmax,
237  flags, flux, fluxerr, fluxrad, elong, vignet,
238  plot=plot, title="%s[%d]" % (filename, ext + 1))
239  #
240  # Insert the good candidates into the set
241  #
242  if not vignet.dtype.isnative:
243  # without the swap setVig fails with
244  # "ValueError: 'unaligned arrays cannot be converted to C++'"
245  vignet = vignet.byteswap()
246 
247  for i in np.where(good)[0]:
248  sample = set.newSample()
249  sample.setCatindex(catindex)
250  sample.setExtindex(ext)
251 
252  sample.setVig(vignet[i])
253 
254  sample.setNorm(float(flux[i]))
255  sample.setBacknoise2(backnoise2)
256  sample.setGain(gain)
257  sample.setX(float(xm[i]))
258  sample.setY(float(ym[i]))
259  sample.setFluxrad(float(fluxrad[i]))
260 
261  for j in range(set.getNcontext()):
262  sample.setContext(j, float(contextvalp[j][i]))
263 
264  set.finiSample(sample, prefs.getProfAccuracy())
265 
266  # ---- Update min and max
267  for j in range(set.getNcontext()):
268  cmin[j] = contextvalp[j][good].min()
269  cmax[j] = contextvalp[j][good].max()
270 
271  # Update the scaling
272  if set.getNsample():
273  for i in range(set.getNcontext()):
274  set.setContextScale(i, cmax[i] - cmin[i])
275  set.setContextOffset(i, (cmin[i] + cmax[i])/2.0)
276 
277  # Don't waste memory!
278  set.trimMemory()
279 
280  return set
281 
282 

◆ 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 763 of file psfex.py.

763 def read_samplesLsst(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, nobj,
764  plot=dict(showFlags=False, showRejection=False)):
765  # allocate a new set iff set is None
766  if not set:
767  set = psfexLib.Set(context)
768 
769  cmin, cmax = None, None
770  if set.getNcontext():
771  cmin = np.empty(set.getNcontext())
772  cmax = np.empty(set.getNcontext())
773  for i in range(set.getNcontext()):
774  if set.getNsample():
775  cmin[i] = set.getContextOffset(i) - set.getContextScale(i)/2.0
776  cmax[i] = cmin[i] + set.getContextScale(i)
777  else:
778  cmin[i] = psfexLib.BIG
779  cmax[i] = -psfexLib.BIG
780  #
781  # Read data
782  #
783  tab = afwTable.SourceCatalog.readFits(filename)
784 
785  centroid = tab.getCentroidDefinition()
786  xm = tab.get("%s.x" % centroid)
787  ym = tab.get("%s.y" % centroid)
788 
789  shape = tab.getShapeDefinition()
790  ixx = tab.get("%s.xx" % shape)
791  iyy = tab.get("%s.yy" % shape)
792 
793  rmsSize = np.sqrt(0.5*(ixx + iyy))
794  elong = 0.5*(ixx - iyy)/(ixx + iyy)
795 
796  flux = tab.get(prefs.getPhotfluxRkey())
797  fluxErr = tab.get(prefs.getPhotfluxerrRkey())
798  flags = getLsstFlags(tab)
799 
800  #
801  # Now the VIGNET data
802  #
803  vigw, vigh = 35, 35 # [prefs.getPsfsize()[i] for i in range(2)]
804  if set.empty():
805  set.setVigSize(vigw, vigh)
806 
807  vignet = np.empty(nobj*vigw*vigh, "float32").reshape(nobj, vigw, vigh)
808 
809  # Hack: I want the WCS so I'll guess where the calexp is to be found
810  calexpFile = guessCalexp(filename)
811  mi = afwImage.MaskedImageF(calexpFile)
812  backnoise2 = np.median(mi.getVariance().getArray())
813  gain = 1.0
814 
815  edgeBit = [k for k, v in getFlags().items() if v == "edge"][0]
816 
817  for i, xc, yc in zip(range(nobj), xm, ym):
818  try:
819  x, y = int(xc), int(yc)
820  except ValueError:
821  flags[i] |= edgeBit # mark star as bad
822 
823  try:
824  pstamp = mi[x - vigw//2:x + vigw//2 + 1, y - vigh//2:y + vigh//2 + 1]
825  vignet[i] = pstamp.getImage().getArray().transpose()
826  except Exception:
827  flags[i] |= edgeBit # mark star as bad
828 
829  # Try to load the set of context keys
830  pc = 0
831  contextvalp = []
832  for i, key in enumerate(context.getName()):
833  if context.getPcflag(i):
834  contextvalp.append(pcval[pc])
835  pc += 1
836  elif key[0] == ':':
837  try:
838  contextvalp.append(tab.header[key[1:]])
839  except KeyError:
840  raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
841  (key[1:], filename))
842  else:
843  try:
844  contextvalp.append(tab.get(key))
845  except KeyError:
846  raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
847  (key, filename))
848  set.setContextname(i, key)
849 
850  # Now examine each vector of the shipment
851  good = select_candidates(set, prefs, frmin, frmax,
852  flags, flux, fluxErr, rmsSize, elong, vignet,
853  plot=plot, title="%s[%d]" % (filename, ext + 1) if next > 1 else filename)
854  #
855  # Insert the good candidates into the set
856  #
857  if not vignet.dtype.isnative:
858  # without the swap setVig fails with
859  # "ValueError: 'unaligned arrays cannot be converted to C++'"
860  vignet = vignet.byteswap()
861 
862  for i in np.where(good)[0]:
863  sample = set.newSample()
864  sample.setCatindex(catindex)
865  sample.setExtindex(ext)
866 
867  sample.setVig(vignet[i])
868  if False:
869  pstamp = afwImage.ImageF(*vignet[i].shape)
870  pstamp.getArray()[:] = sample.getVig()
871  afwDisplay.Display().mtv(pstamp)
872 
873  sample.setNorm(float(flux[i]))
874  sample.setBacknoise2(backnoise2)
875  sample.setGain(gain)
876  sample.setX(float(xm[i]))
877  sample.setY(float(ym[i]))
878  sample.setFluxrad(float(rmsSize[i]))
879 
880  for j in range(set.getNcontext()):
881  sample.setContext(j, float(contextvalp[j][i]))
882 
883  set.finiSample(sample, prefs.getProfAccuracy())
884 
885  # ---- Update min and max
886  for j in range(set.getNcontext()):
887  cmin[j] = contextvalp[j][good].min()
888  cmax[j] = contextvalp[j][good].max()
889 
890  # Update the scaling
891  if set.getNsample():
892  for i in range(set.getNcontext()):
893  set.setContextScale(i, cmax[i] - cmin[i])
894  set.setContextOffset(i, (cmin[i] + cmax[i])/2.0)
895 
896  # Don't waste memory!
897  set.trimMemory()
898 
899  return set
900 
901 

◆ 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 295 of file psfex.py.

295 def select_candidates(set, prefs, frmin, frmax,
296  flags, flux, fluxerr, rmsSize, elong, vignet,
297  plot=dict(), title=""):
298  maxbad = prefs.getBadpixNmax()
299  maxbadflag = prefs.getBadpixFlag()
300  maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100.0
301  minsn = prefs.getMinsn()
302 
303  sn = flux/np.where(fluxerr > 0, fluxerr, 1)
304  sn[fluxerr <= 0] = -psfexLib.BIG
305  # ---- Apply some selection over flags, fluxes...
306  plotFlags = plot.get("showFlags") if plt else False
307  plotRejection = plot.get("showRejection") if plt else False
308 
309  bad = flags & prefs.getFlagMask() != 0
310  set.setBadFlags(int(sum(bad != 0)))
311 
312  if plotRejection:
313  selectionVectors = []
314  selectionVectors.append((bad, "flags %d" % sum(bad != 0)))
315 
316  dbad = sn < minsn
317  set.setBadSN(int(sum(dbad)))
318  bad = np.logical_or(bad, dbad)
319  if plotRejection:
320  selectionVectors.append((dbad, "S/N %d" % sum(dbad)))
321 
322  dbad = rmsSize < frmin
323  set.setBadFrmin(int(sum(dbad)))
324  bad = np.logical_or(bad, dbad)
325  if plotRejection:
326  selectionVectors.append((dbad, "frmin %d" % sum(dbad)))
327 
328  dbad = rmsSize > frmax
329  set.setBadFrmax(int(sum(dbad)))
330  bad = np.logical_or(bad, dbad)
331  if plotRejection:
332  selectionVectors.append((dbad, "frmax %d" % sum(dbad)))
333 
334  dbad = elong > maxelong
335  set.setBadElong(int(sum(dbad)))
336  bad = np.logical_or(bad, dbad)
337  if plotRejection:
338  selectionVectors.append((dbad, "elong %d" % sum(dbad)))
339 
340  # -- ... and check the integrity of the sample
341  if maxbadflag:
342  nbad = np.array([(v <= -psfexLib.BIG).sum() for v in vignet])
343  dbad = nbad > maxbad
344  set.setBadPix(int(sum(dbad)))
345  bad = np.logical_or(bad, dbad)
346  if plotRejection:
347  selectionVectors.append((dbad, "badpix %d" % sum(dbad)))
348 
349  good = np.logical_not(bad)
350  if plotFlags or plotRejection:
351  imag = -2.5*np.log10(flux)
352  plt.clf()
353 
354  alpha = 0.5
355  if plotFlags:
356  labels = getFlags()
357 
358  isSet = np.where(flags == 0x0)[0]
359  plt.plot(imag[isSet], rmsSize[isSet], 'o', alpha=alpha, label="good")
360 
361  for i in range(16):
362  mask = 1 << i
363  if mask & prefs.getFlagMask():
364  isSet = np.where(np.bitwise_and(flags, mask))[0]
365  if isSet.any():
366  plt.plot(imag[isSet], rmsSize[isSet], 'o', alpha=alpha, label=labels[mask])
367  else:
368  for bad, label in selectionVectors:
369  plt.plot(imag[bad], rmsSize[bad], 'o', alpha=alpha, label=label)
370 
371  plt.plot(imag[good], rmsSize[good], 'o', color="black", label="selected")
372  [plt.axhline(_, color='red') for _ in [frmin, frmax]]
373  plt.xlim(np.median(imag[good]) + 5*np.array([-1, 1]))
374  plt.ylim(-0.1, 2*frmax)
375  plt.legend(loc=2)
376  plt.xlabel("Instrumental Magnitude")
377  plt.ylabel("rmsSize")
378  plt.title("%s %d selected" % (title, sum(good)))
379 
380  input("Continue? ")
381 
382  return good
383 
384 
385 try:
386  _dataType

◆ setDataType()

def lsst.meas.extensions.psfex.psfex.setDataType (   t)

Definition at line 400 of file psfex.py.

400 def setDataType(t):
401  _dataType = _dataTypes[t] # noqa: F841
402 
403 

◆ 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 547 of file psfex.py.

547 def showPsf(psf, set, ext=None, wcsData=None, trim=0, nspot=5,
548  diagnostics=False, outDir="", frame=None, title=None):
549  """Show a PSF on display (e.g., ds9)
550  """
551 
552  if ext is not None:
553  psf = psf[ext]
554 
555  if wcsData:
556  if ext is not None:
557  wcsData = wcsData[ext]
558  wcs, naxis1, naxis2 = wcsData
559  else:
560  wcs, naxis1, naxis2 = None, None, None
561 
562  naxis = [naxis1, naxis2]
563  for i in range(2):
564  if naxis[i] is None:
565  # cmin, cmax are the range of input star positions
566  cmin, cmax = [set.getContextOffset(i) + d*set.getContextScale(i) for d in (-0.5, 0.5)]
567  naxis[i] = cmax + cmin # a decent guess
568 
569  if naxis[0] > naxis[1]:
570  nx, ny = int(nspot*naxis[0]/float(naxis[1]) + 0.5), nspot
571  else:
572  nx, ny = nspot, int(nspot*naxis[1]/float(naxis[0]) + 0.5)
573 
574  mos = afwDisplay.utils.Mosaic(gutter=2, background=0.02)
575 
576  xpos, ypos = np.linspace(0, naxis[0], nx), np.linspace(0, naxis[1], ny)
577  for y in ypos:
578  for x in xpos:
579  psf.build(x, y)
580 
581  im = afwImage.ImageF(*psf.getLoc().shape)
582  im.getArray()[:] = psf.getLoc()
583  im /= float(im.getArray().max())
584  if trim:
585  if trim > im.getHeight()//2:
586  trim = im.getHeight()//2
587 
588  im = im[trim:-trim, trim:-trim]
589 
590  mos.append(im)
591 
592  mosaic = mos.makeMosaic(mode=nx)
593  if frame is not None:
594  afwDisplay.Display(frame=frame).mtv(mosaic, title=title)
595  #
596  # Figure out the WCS for the mosaic
597  #
598  pos = []
599  if False:
600  for x, y, i in zip((xpos[0], xpos[-1]), (ypos[0], ypos[-1]), (0, mos.nImage - 1)):
601  bbox = mos.getBBox(i)
602  mosx = bbox.getMinX() + 0.5*(bbox.getWidth() - 1)
603  mosy = bbox.getMinY() + 0.5*(bbox.getHeight() - 1)
604  pos.append([geom.PointD(mosx, mosy), wcs.pixelToSky(geom.PointD(x, y))])
605  else:
606  pos.append([geom.PointD(0, 0), wcs.pixelToSky(geom.PointD(0, 0))])
607  pos.append([geom.PointD(*mosaic.getDimensions()), wcs.pixelToSky(geom.PointD(naxis1, naxis2))])
608 
609  CD = []
610  for i in range(2):
611  delta = pos[1][1][i].asDegrees() - pos[0][1][i].asDegrees()
612  CD.append(delta/(pos[1][0][i] - pos[0][0][i]))
613  CD = np.array(CD)
614  CD.shape = (2, 2)
615  mosWcs = afwGeom.makeSkyWcs(crval=pos[0][0], crpix=pos[0][1], cdMatrix=CD)
616 
617  if ext is not None:
618  title = "%s-%d" % (title, ext)
619 
620  if frame is not None:
621  afwDisplay.Display(frame=frame).mtv(mosaic, title=title, wcs=mosWcs)
622 
623  if diagnostics:
624  outFile = "%s-mod.fits" % title
625  if outDir:
626  outFile = os.path.join(outDir, outFile)
627  mosaic.writeFits(outFile, mosWcs.getFitsMetadata())
628 
629  mos = afwDisplay.utils.Mosaic(gutter=4, background=0.002)
630  for i in range(set.getNsample()):
631  s = set.getSample(i)
632  if ext is not None and s.getExtindex() != ext:
633  continue
634 
635  smos = afwDisplay.utils.Mosaic(gutter=2, background=-0.003)
636  for func in [s.getVig, s.getVigResi]:
637  arr = func()
638  if func == s.getVig:
639  norm = float(arr.max()) if True else s.getNorm()
640 
641  arr /= norm
642  im = afwImage.ImageF(*arr.shape)
643  im.getArray()[:] = arr
644  smos.append(im)
645 
646  mos.append(smos.makeMosaic(mode="x"))
647 
648  mosaic = mos.makeMosaic(title=title)
649 
650  if frame is not None:
651  afwDisplay.Display(frame=frame + 1).mtv(mosaic, title=title)
652 
653  if diagnostics:
654  outFile = "%s-psfstars.fits" % title
655  if outDir:
656  outFile = os.path.join(outDir, outFile)
657 
658  mosaic.writeFits(outFile)
659 
660 

◆ splitFitsCard()

def lsst.meas.extensions.psfex.psfex.splitFitsCard (   line) 
Split a FITS header, returning (key, value).

Parameters
----------
line: `string`
    A string from a FITS header.

Returns
-------
key: `string`
    The header key.
value: `string`
    The header value.

Definition at line 24 of file psfex.py.

24 def splitFitsCard(line):
25  """Split a FITS header, returning (key, value).
26 
27  Parameters
28  ----------
29  line: `string`
30  A string from a FITS header.
31 
32  Returns
33  -------
34  key: `string`
35  The header key.
36  value: `string`
37  The header value.
38  """
39 
40  try:
41  k, v = re.search(r"(\S+)\s*=\s*'?((?:\S+|'))", line).groups()
42  except AttributeError:
43  raise
44 
45  try:
46  v = int(v)
47  except ValueError:
48  try:
49  v = float(v)
50  except ValueError:
51  pass
52 
53  return k, v
54 
55 

Variable Documentation

◆ LSST

lsst.meas.extensions.psfex.psfex.LSST

Definition at line 394 of file psfex.py.

◆ plt

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

Definition at line 10 of file psfex.py.

◆ SExtractor

lsst.meas.extensions.psfex.psfex.SExtractor

Definition at line 395 of file psfex.py.

lsst::meas::extensions::psfex.psfex.read_samplesLsst
def read_samplesLsst(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, nobj, plot=dict(showFlags=False, showRejection=False))
Definition: psfex.py:763
lsst::meas::extensions::psfex.psfex.load_samplesLsst
def load_samplesLsst(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict())
Definition: psfex.py:902
lsst::meas::extensions::psfex.psfex.makeit
def makeit(prefs, context, saveWcs=False, plot=dict())
Definition: psfex.py:1013
lsst::meas::extensions::psfex.psfex.setDataType
def setDataType(t)
Definition: psfex.py:400
lsst::afw.display.ds9.mtv
def mtv(data, frame=None, title="", wcs=None, *args, **kwargs)
Definition: ds9.py:93
lsst::meas::extensions::psfex.psfex.getSexFlags
def getSexFlags(*args)
Definition: psfex.py:283
lsst::meas::extensions::psfex.psfex.makeitLsst
def makeitLsst(prefs, context, saveWcs=False, plot=dict())
Definition: psfex.py:706
lsst::meas::extensions::psfex.psfex.splitFitsCard
def splitFitsCard(line)
Definition: psfex.py:24
lsst::meas::extensions::psfex.psfex.getFlags
def getFlags()
Definition: psfex.py:404
lsst::meas::extensions::psfex.psfex.guessCalexp
def guessCalexp(fileName)
Definition: psfex.py:695
max
int max
Definition: BoundedField.cc:104
lsst::meas::extensions::psfex.psfex.compute_fwhmrange
def compute_fwhmrange(fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False))
Definition: psfex.py:56
lsst::afw::geom::makeSkyWcs
std::shared_ptr< SkyWcs > makeSkyWcs(TransformPoint2ToPoint2 const &pixelsToFieldAngle, lsst::geom::Angle const &orientation, bool flipX, lsst::geom::SpherePoint const &boresight, std::string const &projection="TAN")
Construct a FITS SkyWcs from camera geometry.
Definition: SkyWcs.cc:536
lsst::meas::extensions::psfex.psfex.select_candidates
def select_candidates(set, prefs, frmin, frmax, flags, flux, fluxerr, rmsSize, elong, vignet, plot=dict(), title="")
Definition: psfex.py:295
lsst::afw::image.readMetadata.readMetadataContinued.readMetadata
readMetadata
Definition: readMetadataContinued.py:28
items
std::vector< SchemaItem< Flag > > * items
Definition: BaseColumnView.cc:142
list
daf::base::PropertyList * list
Definition: fits.cc:913
min
int min
Definition: BoundedField.cc:103
lsst::geom::Point< double, 2 >
load_samples
setstruct * load_samples(char **filenames, int catindex, int ncat, int ext, int next, contextstruct *context)
Definition: psfexAdaptors.cc:14
lsst::meas::extensions::psfex.psfex.getLsstFlags
def getLsstFlags(tab=None)
Definition: psfex.py:661
lsst::geom::SpherePoint
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
lsst::meas::extensions::psfex.psfex.read_samples
def read_samples(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, plot=dict(showFlags=False, showRejection=False))
Definition: psfex.py:123
lsst::meas::extensions::psfex.psfex.showPsf
def showPsf(psf, set, ext=None, wcsData=None, trim=0, nspot=5, diagnostics=False, outDir="", frame=None, title=None)
Definition: psfex.py:547
Generated on Wed Jun 24 2020 18:10:40 for LSSTApplications by   doxygen 1.8.18