LSSTApplications  20.0.0
LSSTDataManagementBasePackage
psfex.py
Go to the documentation of this file.
1 import os
2 import re
3 import sys
4 
5 import numpy as np
6 from astropy.io import fits
7 try:
8  import matplotlib.pyplot as plt
9 except ImportError:
10  plt = None
11 
12 import lsst.geom as geom
13 import lsst.afw.geom as afwGeom
14 from lsst.afw.fits import readMetadata
15 import lsst.afw.image as afwImage
16 import lsst.afw.table as afwTable
17 import lsst.afw.display as afwDisplay
18 from lsst.daf.base import PropertySet
19 from . import psfexLib
20 
21 afwDisplay.setDefaultMaskTransparency(75)
22 
23 
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 
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 
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 
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 
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
387 except NameError:
388  class _SExtractor():
389  pass
390 
391  class _LSST():
392  pass
393 
394  _dataTypes = dict(LSST=_LSST,
395  SExtractor=_SExtractor
396  )
397  _dataType = _SExtractor
398 
399 
400 def setDataType(t):
401  _dataType = _dataTypes[t] # noqa: F841
402 
403 
404 def getFlags():
405  if _dataType == _LSST:
406  return getLsstFlags(None)
407  else:
408  return getSexFlags(None)
409 
410 
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 
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 
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 
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 
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 
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 
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 
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
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
lsst::meas::extensions::psfex.psfex._SExtractor
Definition: psfex.py:388
lsst::meas::extensions::psfex.psfex._LSST
Definition: psfex.py:391
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::afw.display
Definition: __init__.py:1
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.load_samples
def load_samples(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict())
Definition: psfex.py:411
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::afw::table
Definition: table.dox:3
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
lsst::geom
Definition: geomOperators.dox:4
lsst::daf::base
Definition: Utils.h:47
items
std::vector< SchemaItem< Flag > > * items
Definition: BaseColumnView.cc:142
lsst::afw::fits
Definition: fits.h:31
list
daf::base::PropertyList * list
Definition: fits.cc:913
min
int min
Definition: BoundedField.cc:103
lsst::geom::Point< double, 2 >
lsst::daf::base::PropertySet
Class for storing generic metadata.
Definition: PropertySet.h:67
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
lsst::afw::geom
Definition: frameSetUtils.h:40
Generated on Wed Jun 24 2020 18:10:10 for LSSTApplications by   doxygen 1.8.18