LSSTApplications  18.1.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.afw.geom as afwGeom
13 from lsst.afw.fits import readMetadata
14 import lsst.afw.image as afwImage
15 import lsst.afw.table as afwTable
16 import lsst.afw.display as afwDisplay
17 from lsst.daf.base import PropertySet
18 from . import psfexLib
19 
20 afwDisplay.setDefaultMaskTransparency(75)
21 
22 
23 def splitFitsCard(line):
24  """Split a FITS header, returning (key, value).
25 
26  Parameters
27  ----------
28  line: `string`
29  A string from a FITS header.
30 
31  Returns
32  -------
33  key: `string`
34  The header key.
35  value: `string`
36  The header value.
37  """
38 
39  try:
40  k, v = re.search(r"(\S+)\s*=\s*'?((?:\S+|'))", line).groups()
41  except AttributeError:
42  raise
43 
44  try:
45  v = int(v)
46  except ValueError:
47  try:
48  v = float(v)
49  except ValueError:
50  pass
51 
52  return k, v
53 
54 
55 def compute_fwhmrange(fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False)):
56  """Compute the FWHM range associated to a series of FWHM measurements.
57  AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
58  VERSION 20/03/2008
59 
60  Parameters
61  ----------
62  fwhm: iterable of `float`
63  Iterable of full width half-maxima.
64  maxvar: `float`
65  Maximum allowed FWHM variation.
66  minin: `float`
67  Minimum allowed FWHM.
68  maxin: `float`
69  Maximum allowed FWHM.
70  plot: `dict`, optional
71  Dict of plotting options.
72 
73  Returns
74  -------
75  fmin: `float`
76  FWHM mode.
77  minout: `float`
78  Lower FWHM range.
79  maxout: `float`
80  Upper FWHM range.
81  """
82 
83  nfwhm = len(fwhm)
84  fwhm.sort()
85 
86  # Find the mode
87  nw = nfwhm//4
88  if nw < 4:
89  nw = 1
90  dfmin = psfexLib.BIG
91  fmin = 0.0
92  for i in range(nfwhm - nw):
93  df = fwhm[i + nw] - fwhm[i]
94  if df < dfmin:
95  dfmin = df
96  fmin = (fwhm[i + nw] + fwhm[i])/2.0
97 
98  if nfwhm < 2:
99  fmin = fwhm[0]
100 
101  dfmin = (maxvar + 1.0)**0.3333333
102  minout = fmin/dfmin if dfmin > 0.0 else 0.0
103  if minout < minin:
104  minout = minin
105 
106  maxout = fmin*dfmin**2
107  if maxout > maxin:
108  maxout = maxin
109 
110  if plt and plot.get("fwhmHistogram"):
111  plt.clf()
112  plt.hist(fwhm, nfwhm//10 + 1, normed=1, facecolor='g', alpha=0.75)
113  plt.xlabel("FWHM")
114  plt.axvline(fmin, color='red')
115  [plt.axvline(_, color='blue') for _ in (minout, maxout)]
116 
117  input("Continue? ")
118 
119  return fmin, minout, maxout
120 
121 
122 def read_samples(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval,
123  plot=dict(showFlags=False, showRejection=False)):
124  # allocate a new set iff set is None
125  if not set:
126  set = psfexLib.Set(context)
127 
128  cmin, cmax = None, None
129  if set.getNcontext():
130  cmin = np.empty(set.getNcontext())
131  cmax = np.empty(set.getNcontext())
132  for i in range(set.getNcontext()):
133  if set.getNsample():
134  cmin[i] = set.getContextOffset(i) - set.getContextScale(i)/2.0
135  cmax[i] = cmin[i] + set.getContextScale(i)
136  else:
137  cmin[i] = psfexLib.BIG
138  cmax[i] = -psfexLib.BIG
139  #
140  # Read data
141  #
142  with fits.open(filename) as cat:
143  extCtr = -1
144  for tab in cat:
145  if tab.name == "LDAC_IMHEAD":
146  extCtr += 1
147 
148  if extCtr < ext:
149  continue
150  elif extCtr > ext:
151  break
152 
153  if tab.name == "PRIMARY":
154  pass
155  elif tab.name == "LDAC_IMHEAD":
156  hdr = tab.data[0][0] # the fits header from the original fits image
157  foundCards = 0 # how many of our desired cards have we found?
158 
159  for line in hdr:
160  try:
161  k, v = splitFitsCard(line)
162  except AttributeError:
163  continue
164 
165  if k == "SEXBKDEV":
166  backnoise2 = v**2
167  foundCards += 1
168  elif k == "SEXGAIN":
169  gain = v
170  foundCards += 1
171 
172  if foundCards == 2:
173  break
174  elif tab.name == "LDAC_OBJECTS":
175  xm = tab.data[prefs.getCenterKey(0)]
176  ym = tab.data[prefs.getCenterKey(1)]
177  fluxrad = tab.data["FLUX_RADIUS"]
178  flux = tab.data[prefs.getPhotfluxRkey()]
179  fluxerr = tab.data[prefs.getPhotfluxerrRkey()]
180  elong = tab.data["ELONGATION"]
181  flags = tab.data["FLAGS"]
182 
183  n = prefs.getPhotfluxNum() - 1
184  if n:
185  assert False, "Code to handle e.g. FLUX_APER(3) isn't yet converted"
186  if key.naxis == 1 and n < key.naxisn[0]:
187  flux += n
188  else:
189  print("Not enough apertures for %s in catalogue %s: using first aperture" %
190  (prefs.getPhotfluxRkey(), filename), file=sys.stderr)
191 
192  n = prefs.getPhotfluxerrNum() - 1
193  if n:
194  if key.naxis == 1 and n < key.naxisn[0]:
195  fluxerr += n
196  else:
197  print("Not enough apertures for %s in catalogue %s: using first aperture" %
198  (prefs.getPhotfluxerrRkey(), filename), file=sys.stderr)
199  #
200  # Now the VIGNET data
201  #
202  vignet = tab.data["VIGNET"]
203 
204  try:
205  vigw, vigh = vignet[0].shape
206  except ValueError:
207  raise RuntimeError("*Error*: VIGNET should be a 2D vector; saw %s" % str(vignet[0].shape))
208 
209  if set.empty():
210  set.setVigSize(vigw, vigh)
211 
212  # Try to load the set of context keys
213  pc = 0
214  contextvalp = []
215  for i, key in enumerate(context.getName()):
216  if context.getPcflag(i):
217  contextvalp.append(pcval[pc])
218  pc += 1
219  elif key[0] == ':':
220  try:
221  contextvalp.append(tab.header[key[1:]])
222  except KeyError:
223  raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
224  (key[1:], filename))
225  else:
226  try:
227  contextvalp.append(tab.data[key])
228  except KeyError:
229  raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
230  (key, filename))
231  set.setContextname(i, key)
232 
233  # Now examine each vector of the shipment
234  good = select_candidates(set, prefs, frmin, frmax,
235  flags, flux, fluxerr, fluxrad, elong, vignet,
236  plot=plot, title="%s[%d]" % (filename, ext + 1))
237  #
238  # Insert the good candidates into the set
239  #
240  if not vignet.dtype.isnative:
241  # without the swap setVig fails with "ValueError: 'unaligned arrays cannot be converted to C++'"
242  vignet = vignet.byteswap()
243 
244  for i in np.where(good)[0]:
245  sample = set.newSample()
246  sample.setCatindex(catindex)
247  sample.setExtindex(ext)
248 
249  sample.setVig(vignet[i])
250 
251  sample.setNorm(float(flux[i]))
252  sample.setBacknoise2(backnoise2)
253  sample.setGain(gain)
254  sample.setX(float(xm[i]))
255  sample.setY(float(ym[i]))
256  sample.setFluxrad(float(fluxrad[i]))
257 
258  for j in range(set.getNcontext()):
259  sample.setContext(j, float(contextvalp[j][i]))
260 
261  set.finiSample(sample, prefs.getProfAccuracy())
262 
263  # ---- Update min and max
264  for j in range(set.getNcontext()):
265  cmin[j] = contextvalp[j][good].min()
266  cmax[j] = contextvalp[j][good].max()
267 
268  # Update the scaling
269  if set.getNsample():
270  for i in range(set.getNcontext()):
271  set.setContextScale(i, cmax[i] - cmin[i])
272  set.setContextOffset(i, (cmin[i] + cmax[i])/2.0)
273 
274  # Don't waste memory!
275  set.trimMemory()
276 
277  return set
278 
279 
280 def getSexFlags(*args):
281  return {1: "flux blended",
282  2: "blended",
283  4: "saturated",
284  8: "edge",
285  16: "bad aperture",
286  32: "bad isophotal",
287  64: "memory error (deblend)",
288  128: "memory error (extract)",
289  }
290 
291 
292 def select_candidates(set, prefs, frmin, frmax,
293  flags, flux, fluxerr, rmsSize, elong, vignet,
294  plot=dict(), title=""):
295  maxbad = prefs.getBadpixNmax()
296  maxbadflag = prefs.getBadpixFlag()
297  maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100.0
298  minsn = prefs.getMinsn()
299 
300  sn = flux/np.where(fluxerr > 0, fluxerr, 1)
301  sn[fluxerr <= 0] = -psfexLib.BIG
302  # ---- Apply some selection over flags, fluxes...
303  plotFlags = plot.get("showFlags") if plt else False
304  plotRejection = plot.get("showRejection") if plt else False
305 
306  bad = flags & prefs.getFlagMask() != 0
307  set.setBadFlags(int(sum(bad != 0)))
308 
309  if plotRejection:
310  selectionVectors = []
311  selectionVectors.append((bad, "flags %d" % sum(bad != 0)))
312 
313  dbad = sn < minsn
314  set.setBadSN(int(sum(dbad)))
315  bad = np.logical_or(bad, dbad)
316  if plotRejection:
317  selectionVectors.append((dbad, "S/N %d" % sum(dbad)))
318 
319  dbad = rmsSize < frmin
320  set.setBadFrmin(int(sum(dbad)))
321  bad = np.logical_or(bad, dbad)
322  if plotRejection:
323  selectionVectors.append((dbad, "frmin %d" % sum(dbad)))
324 
325  dbad = rmsSize > frmax
326  set.setBadFrmax(int(sum(dbad)))
327  bad = np.logical_or(bad, dbad)
328  if plotRejection:
329  selectionVectors.append((dbad, "frmax %d" % sum(dbad)))
330 
331  dbad = elong > maxelong
332  set.setBadElong(int(sum(dbad)))
333  bad = np.logical_or(bad, dbad)
334  if plotRejection:
335  selectionVectors.append((dbad, "elong %d" % sum(dbad)))
336 
337  # -- ... and check the integrity of the sample
338  if maxbadflag:
339  nbad = np.array([(v <= -psfexLib.BIG).sum() for v in vignet])
340  dbad = nbad > maxbad
341  set.setBadPix(int(sum(dbad)))
342  bad = np.logical_or(bad, dbad)
343  if plotRejection:
344  selectionVectors.append((dbad, "badpix %d" % sum(dbad)))
345 
346  good = np.logical_not(bad)
347  if plotFlags or plotRejection:
348  imag = -2.5*np.log10(flux)
349  plt.clf()
350 
351  alpha = 0.5
352  if plotFlags:
353  labels = getFlags()
354 
355  isSet = np.where(flags == 0x0)[0]
356  plt.plot(imag[isSet], rmsSize[isSet], 'o', alpha=alpha, label="good")
357 
358  for i in range(16):
359  mask = 1 << i
360  if mask & prefs.getFlagMask():
361  isSet = np.where(np.bitwise_and(flags, mask))[0]
362  if isSet.any():
363  plt.plot(imag[isSet], rmsSize[isSet], 'o', alpha=alpha, label=labels[mask])
364  else:
365  for bad, label in selectionVectors:
366  plt.plot(imag[bad], rmsSize[bad], 'o', alpha=alpha, label=label)
367 
368  plt.plot(imag[good], rmsSize[good], 'o', color="black", label="selected")
369  [plt.axhline(_, color='red') for _ in [frmin, frmax]]
370  plt.xlim(np.median(imag[good]) + 5*np.array([-1, 1]))
371  plt.ylim(-0.1, 2*frmax)
372  plt.legend(loc=2)
373  plt.xlabel("Instrumental Magnitude")
374  plt.ylabel("rmsSize")
375  plt.title("%s %d selected" % (title, sum(good)))
376 
377  input("Continue? ")
378 
379  return good
380 
381 
382 try:
383  _dataType
384 except NameError:
385  class _SExtractor():
386  pass
387 
388  class _LSST():
389  pass
390 
391  _dataTypes = dict(LSST=_LSST,
392  SExtractor=_SExtractor
393  )
394  _dataType = _SExtractor
395 
396 
397 def setDataType(t):
398  _dataType = _dataTypes[t]
399 
400 
401 def getFlags():
402  if _dataType == _LSST:
403  return getLsstFlags(None)
404  else:
405  return getSexFlags(None)
406 
407 
408 def load_samples(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict()):
409  minsn = prefs.getMinsn()
410  maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100
411  min = prefs.getFwhmrange()[0]
412  max = prefs.getFwhmrange()[1]
413 
414  filenames = prefs.getCatalogs()
415 
416  ncat = len(filenames)
417  fwhmmin = np.empty(ncat)
418  fwhmmax = np.empty(ncat)
419 
420  if not prefs.getAutoselectFlag():
421  fwhmmin = np.zeros(ncat) + prefs.getFwhmrange()[0]
422  fwhmmax = np.zeros(ncat) + prefs.getFwhmrange()[1]
423  fwhmmode = (fwhmmin + fwhmmax)/2.0
424  else:
425  fwhms = {}
426 
427  # -- Try to estimate the most appropriate Half-light Radius range
428  # -- Get the Half-light radii
429  nobj = 0
430  for i, fileName in enumerate(filenames):
431  fwhms[i] = []
432 
433  if prefs.getVerboseType() != prefs.QUIET:
434  print("Examining Catalog #%d" % (i+1))
435 
436  # ---- Read input catalog
437  backnoises = []
438  with fits.open(fileName) as cat:
439  extCtr = -1
440  for tab in cat:
441  if tab.name == "LDAC_IMHEAD":
442  extCtr += 1
443 
444  if extCtr != ext and ext != prefs.ALL_EXTENSIONS:
445  if extCtr > ext:
446  break
447  continue
448 
449  if tab.name == "PRIMARY":
450  pass
451  elif tab.name == "LDAC_IMHEAD":
452  hdr = tab.data[0][0] # the fits header from the original fits image
453  for line in hdr:
454  try:
455  k, v = splitFitsCard(line)
456  except AttributeError:
457  continue
458 
459  if k == "SEXBKDEV":
460  if v < 1/psfexLib.BIG:
461  v = 1.0
462 
463  backnoises.append(v)
464  break
465  elif tab.name == "LDAC_OBJECTS":
466  # -------- Fill the FWHM array
467  rmsSize = tab.data["FLUX_RADIUS"]
468  fmax = tab.data["FLUX_MAX"]
469  flags = tab.data["FLAGS"]
470  elong = tab.data["ELONGATION"]
471  backnoise = backnoises[-1]
472 
473  good = np.logical_and(fmax/backnoise > minsn,
474  np.logical_not(flags & prefs.getFlagMask()))
475  good = np.logical_and(good, elong < maxelong)
476  fwhm = 2.0*rmsSize
477  good = np.logical_and(good, fwhm >= min)
478  good = np.logical_and(good, fwhm < max)
479  fwhms[i] = fwhm[good]
480 
481  if prefs.getVarType() == prefs.VAR_NONE:
482  if nobj:
483  fwhms_all = np.empty(sum([len(l) for l in fwhms.values()]))
484  i = 0
485  for l in fwhms.values():
486  fwhms_all[i:len(l)] = l
487  i += len(l)
488  mode, min, max = compute_fwhmrange(fwhms_all, prefs.getMaxvar(),
489  prefs.getFwhmrange()[0], prefs.getFwhmrange()[1],
490  plot=plot)
491  else:
492  raise RuntimeError("No source with appropriate FWHM found!!")
493  mode = min = max = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
494 
495  fwhmmin = np.zeros(ncat) + min
496  fwhmmax = np.zeros(ncat) + max
497  fwhmmode = np.zeros(ncat) + mode
498  else:
499  fwhmmode = np.empty(ncat)
500  fwhmmin = np.empty(ncat)
501  fwhmmax = np.empty(ncat)
502 
503  for i in range(ncat):
504  nobj = len(fwhms[i])
505  if (nobj):
506  fwhmmode[i], fwhmmin[i], fwhmmax[i] = \
507  compute_fwhmrange(fwhms[i], prefs.getMaxvar(),
508  prefs.getFwhmrange()[0], prefs.getFwhmrange()[1], plot=plot)
509  else:
510  raise RuntimeError("No source with appropriate FWHM found!!")
511  fwhmmode[i] = fwhmmin[i] = fwhmmax[i] = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
512 
513  # Read the samples
514  mode = psfexLib.BIG # mode of FWHM distribution
515 
516  sets = []
517  for i, fileName in enumerate(filenames):
518  set = None
519  if ext == prefs.ALL_EXTENSIONS:
520  extensions = list(range(len(backnoises)))
521  else:
522  extensions = [ext]
523 
524  for e in extensions:
525  set = read_samples(prefs, set, fileName, fwhmmin[i]/2.0, fwhmmax[i]/2.0,
526  e, next, i, context, context.getPc(i) if context.getNpc() else None, plot=plot)
527 
528  if fwhmmode[i] < mode:
529  mode = fwhmmode[i]
530 
531  set.setFwhm(mode)
532 
533  if prefs.getVerboseType() != prefs.QUIET:
534  if set.getNsample():
535  print("%d samples loaded." % set.getNsample())
536  else:
537  raise RuntimeError("No appropriate source found!!")
538 
539  sets.append(set)
540 
541  return sets
542 
543 
544 def showPsf(psf, set, ext=None, wcsData=None, trim=0, nspot=5,
545  diagnostics=False, outDir="", frame=None, title=None):
546  """Show a PSF on display (e.g., ds9)
547  """
548 
549  if ext is not None:
550  psf = psf[ext]
551 
552  if wcsData:
553  if ext is not None:
554  wcsData = wcsData[ext]
555  wcs, naxis1, naxis2 = wcsData
556  else:
557  wcs, naxis1, naxis2 = None, None, None
558 
559  naxis = [naxis1, naxis2]
560  for i in range(2):
561  if naxis[i] is None:
562  # cmin, cmax are the range of input star positions
563  cmin, cmax = [set.getContextOffset(i) + d*set.getContextScale(i) for d in (-0.5, 0.5)]
564  naxis[i] = cmax + cmin # a decent guess
565 
566  if naxis[0] > naxis[1]:
567  nx, ny = int(nspot*naxis[0]/float(naxis[1]) + 0.5), nspot
568  else:
569  nx, ny = nspot, int(nspot*naxis[1]/float(naxis[0]) + 0.5)
570 
571  mos = afwDisplay.utils.Mosaic(gutter=2, background=0.02)
572 
573  xpos, ypos = np.linspace(0, naxis[0], nx), np.linspace(0, naxis[1], ny)
574  for y in ypos:
575  for x in xpos:
576  psf.build(x, y)
577 
578  im = afwImage.ImageF(*psf.getLoc().shape)
579  im.getArray()[:] = psf.getLoc()
580  im /= float(im.getArray().max())
581  if trim:
582  if trim > im.getHeight()//2:
583  trim = im.getHeight()//2
584 
585  im = im[trim:-trim, trim:-trim]
586 
587  mos.append(im)
588 
589  mosaic = mos.makeMosaic(mode=nx)
590  if frame is not None:
591  afwDisplay.Display(frame=frame).mtv(mosaic, title=title)
592  #
593  # Figure out the WCS for the mosaic
594  #
595  pos = []
596  if False:
597  for x, y, i in zip((xpos[0], xpos[-1]), (ypos[0], ypos[-1]), (0, mos.nImage - 1)):
598  bbox = mos.getBBox(i)
599  mosx = bbox.getMinX() + 0.5*(bbox.getWidth() - 1)
600  mosy = bbox.getMinY() + 0.5*(bbox.getHeight() - 1)
601  pos.append([afwGeom.PointD(mosx, mosy), wcs.pixelToSky(afwGeom.PointD(x, y))])
602  else:
603  pos.append([afwGeom.PointD(0, 0), wcs.pixelToSky(afwGeom.PointD(0, 0))])
604  pos.append([afwGeom.PointD(*mosaic.getDimensions()), wcs.pixelToSky(afwGeom.PointD(naxis1, naxis2))])
605 
606  CD = []
607  for i in range(2):
608  delta = pos[1][1][i].asDegrees() - pos[0][1][i].asDegrees()
609  CD.append(delta/(pos[1][0][i] - pos[0][0][i]))
610  CD = np.array(CD)
611  CD.shape = (2, 2)
612  mosWcs = afwGeom.makeSkyWcs(crval=pos[0][0], crpix=pos[0][1], cdMatrix=CD)
613 
614  if ext is not None:
615  title = "%s-%d" % (title, ext)
616 
617  if frame is not None:
618  afwDisplay.Display(frame=frame).mtv(mosaic, title=title, wcs=mosWcs)
619 
620  if diagnostics:
621  outFile = "%s-mod.fits" % title
622  if outDir:
623  outFile = os.path.join(outDir, outFile)
624  mosaic.writeFits(outFile, mosWcs.getFitsMetadata())
625 
626  mos = afwDisplay.utils.Mosaic(gutter=4, background=0.002)
627  for i in range(set.getNsample()):
628  s = set.getSample(i)
629  if ext is not None and s.getExtindex() != ext:
630  continue
631 
632  smos = afwDisplay.utils.Mosaic(gutter=2, background=-0.003)
633  for func in [s.getVig, s.getVigResi]:
634  arr = func()
635  if func == s.getVig:
636  norm = float(arr.max()) if True else s.getNorm()
637 
638  arr /= norm
639  im = afwImage.ImageF(*arr.shape)
640  im.getArray()[:] = arr
641  smos.append(im)
642 
643  mos.append(smos.makeMosaic(mode="x"))
644 
645  mosaic = mos.makeMosaic(title=title)
646 
647  if frame is not None:
648  afwDisplay.Display(frame=frame + 1).mtv(mosaic, title=title)
649 
650  if diagnostics:
651  outFile = "%s-psfstars.fits" % title
652  if outDir:
653  outFile = os.path.join(outDir, outFile)
654 
655  mosaic.writeFits(outFile)
656 
657 
658 def getLsstFlags(tab=None):
659  flagKeys = [
660  "base_PixelFlags_flag_edge",
661  # "base_PixelFlags_flag_interpolated",
662  # "base_PixelFlags_flag_interpolatedCenter",
663  # "base_PixelFlags_flag_saturated",
664  "base_PixelFlags_flag_saturatedCenter",
665  # "base_PixelFlags_flag_cr",
666  "base_PixelFlags_flag_crCenter",
667  "base_PixelFlags_flag_bad",
668  "base_PsfFlux_flag",
669  "parent",
670  ]
671 
672  if tab is None:
673  flags = {}
674  for i, k in enumerate(flagKeys):
675  flags[1 << i] = re.sub(r"\_flag", "",
676  re.sub(r"^base\_", "", re.sub(r"^base\_PixelFlags\_flag\_", "", k)))
677  else:
678  flags = 0
679  for i, k in enumerate(flagKeys):
680  if k == "parent":
681  try:
682  isSet = tab.get("deblend_nChild") > 0
683  except KeyError:
684  isSet = 0
685  else:
686  isSet = tab.get(k)
687  flags = np.bitwise_or(flags, np.where(isSet, 1 << i, 0))
688 
689  return flags
690 
691 
692 def guessCalexp(fileName):
693  for guess in [
694  re.sub("/src", r"", fileName),
695  re.sub("(SRC([^.]+))", r"CORR\2-exp", fileName),
696  ]:
697  if guess != fileName and os.path.exists(guess):
698  return guess
699 
700  raise RuntimeError("Unable to find a calexp to go with %s" % fileName)
701 
702 
703 def makeitLsst(prefs, context, saveWcs=False, plot=dict()):
704  """This is the python wrapper that reads lsst tables
705  """
706  # Create an array of PSFs (one PSF for each extension)
707  if prefs.getVerboseType() != prefs.QUIET:
708  print("----- %d input catalogues:" % prefs.getNcat())
709 
710  if saveWcs: # only needed for making plots
711  wcssList = []
712 
713  fields = psfexLib.vectorField()
714  for cat in prefs.getCatalogs():
715  field = psfexLib.Field(cat)
716  wcss = []
717  wcssList.append(wcss)
718  with fits.open(cat):
719  # Hack: I want the WCS so I'll guess where the calexp is to be found
720  calexpFile = guessCalexp(cat)
721  md = readMetadata(calexpFile)
722  wcs = afwGeom.makeSkyWcs(md)
723 
724  if not wcs:
725  cdMatrix = np.array([1.0, 0.0, 0.0, 1.0])
726  cdMatrix.shape = (2, 2)
727  wcs = afwGeom.makeSkyWcs(crpix=afwGeom.PointD(0, 0),
728  crval=afwGeom.SpherePoint(0.0, 0.0, afwGeom.degrees),
729  cdMatrix=cdMatrix)
730 
731  naxis1, naxis2 = md.getScalar("NAXIS1"), md.getScalar("NAXIS2")
732  # Find how many rows there are in the catalogue
733  md = readMetadata(cat)
734 
735  field.addExt(wcs, naxis1, naxis2, md.getScalar("NAXIS2"))
736  if saveWcs:
737  wcss.append((wcs, naxis1, naxis2))
738 
739  field.finalize()
740  fields.append(field)
741 
742  fields[0].getNext() # number of extensions
743 
744  prefs.getPsfStep()
745 
746  sets = psfexLib.vectorSet()
747  for set in load_samplesLsst(prefs, context, plot=plot):
748  sets.append(set)
749 
750  psfexLib.makeit(fields, sets)
751 
752  ret = [[f.getPsfs() for f in fields], sets]
753  if saveWcs:
754  ret.append(wcssList)
755 
756  return ret
757 
758 
759 def read_samplesLsst(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, nobj,
760  plot=dict(showFlags=False, showRejection=False)):
761  # allocate a new set iff set is None
762  if not set:
763  set = psfexLib.Set(context)
764 
765  cmin, cmax = None, None
766  if set.getNcontext():
767  cmin = np.empty(set.getNcontext())
768  cmax = np.empty(set.getNcontext())
769  for i in range(set.getNcontext()):
770  if set.getNsample():
771  cmin[i] = set.getContextOffset(i) - set.getContextScale(i)/2.0
772  cmax[i] = cmin[i] + set.getContextScale(i)
773  else:
774  cmin[i] = psfexLib.BIG
775  cmax[i] = -psfexLib.BIG
776  #
777  # Read data
778  #
779  tab = afwTable.SourceCatalog.readFits(filename)
780 
781  centroid = tab.getCentroidDefinition()
782  xm = tab.get("%s.x" % centroid)
783  ym = tab.get("%s.y" % centroid)
784 
785  shape = tab.getShapeDefinition()
786  ixx = tab.get("%s.xx" % shape)
787  iyy = tab.get("%s.yy" % shape)
788 
789  rmsSize = np.sqrt(0.5*(ixx + iyy))
790  elong = 0.5*(ixx - iyy)/(ixx + iyy)
791 
792  flux = tab.get(prefs.getPhotfluxRkey())
793  fluxErr = tab.get(prefs.getPhotfluxerrRkey())
794  flags = getLsstFlags(tab)
795 
796  #
797  # Now the VIGNET data
798  #
799  vigw, vigh = 35, 35 # [prefs.getPsfsize()[i] for i in range(2)]
800  if set.empty():
801  set.setVigSize(vigw, vigh)
802 
803  vignet = np.empty(nobj*vigw*vigh, "float32").reshape(nobj, vigw, vigh)
804 
805  # Hack: I want the WCS so I'll guess where the calexp is to be found
806  calexpFile = guessCalexp(filename)
807  mi = afwImage.MaskedImageF(calexpFile)
808  backnoise2 = np.median(mi.getVariance().getArray())
809  gain = 1.0
810 
811  edgeBit = [k for k, v in getFlags().items() if v == "edge"][0]
812 
813  for i, xc, yc in zip(range(nobj), xm, ym):
814  try:
815  x, y = int(xc), int(yc)
816  except ValueError:
817  flags[i] |= edgeBit # mark star as bad
818 
819  try:
820  pstamp = mi[x - vigw//2:x + vigw//2 + 1, y - vigh//2:y + vigh//2 + 1]
821  vignet[i] = pstamp.getImage().getArray().transpose()
822  except Exception:
823  flags[i] |= edgeBit # mark star as bad
824 
825  # Try to load the set of context keys
826  pc = 0
827  contextvalp = []
828  for i, key in enumerate(context.getName()):
829  if context.getPcflag(i):
830  contextvalp.append(pcval[pc])
831  pc += 1
832  elif key[0] == ':':
833  try:
834  contextvalp.append(tab.header[key[1:]])
835  except KeyError:
836  raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
837  (key[1:], filename))
838  else:
839  try:
840  contextvalp.append(tab.get(key))
841  except KeyError:
842  raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
843  (key, filename))
844  set.setContextname(i, key)
845 
846  # Now examine each vector of the shipment
847  good = select_candidates(set, prefs, frmin, frmax,
848  flags, flux, fluxErr, rmsSize, elong, vignet,
849  plot=plot, title="%s[%d]" % (filename, ext + 1) if next > 1 else filename)
850  #
851  # Insert the good candidates into the set
852  #
853  if not vignet.dtype.isnative:
854  # without the swap setVig fails with "ValueError: 'unaligned arrays cannot be converted to C++'"
855  vignet = vignet.byteswap()
856 
857  for i in np.where(good)[0]:
858  sample = set.newSample()
859  sample.setCatindex(catindex)
860  sample.setExtindex(ext)
861 
862  sample.setVig(vignet[i])
863  if False:
864  pstamp = afwImage.ImageF(*vignet[i].shape)
865  pstamp.getArray()[:] = sample.getVig()
866  afwDisplay.Display().mtv(pstamp)
867 
868  sample.setNorm(float(flux[i]))
869  sample.setBacknoise2(backnoise2)
870  sample.setGain(gain)
871  sample.setX(float(xm[i]))
872  sample.setY(float(ym[i]))
873  sample.setFluxrad(float(rmsSize[i]))
874 
875  for j in range(set.getNcontext()):
876  sample.setContext(j, float(contextvalp[j][i]))
877 
878  set.finiSample(sample, prefs.getProfAccuracy())
879 
880  # ---- Update min and max
881  for j in range(set.getNcontext()):
882  cmin[j] = contextvalp[j][good].min()
883  cmax[j] = contextvalp[j][good].max()
884 
885  # Update the scaling
886  if set.getNsample():
887  for i in range(set.getNcontext()):
888  set.setContextScale(i, cmax[i] - cmin[i])
889  set.setContextOffset(i, (cmin[i] + cmax[i])/2.0)
890 
891  # Don't waste memory!
892  set.trimMemory()
893 
894  return set
895 
896 
897 def load_samplesLsst(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict()):
898  minsn = prefs.getMinsn()
899  maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100
900  min = prefs.getFwhmrange()[0]
901  max = prefs.getFwhmrange()[1]
902 
903  filenames = prefs.getCatalogs()
904 
905  ncat = len(filenames)
906  fwhmmin = np.empty(ncat)
907  fwhmmax = np.empty(ncat)
908 
909  if not prefs.getAutoselectFlag():
910  fwhmmin = np.zeros(ncat) + prefs.getFwhmrange()[0]
911  fwhmmax = np.zeros(ncat) + prefs.getFwhmrange()[1]
912  fwhmmode = (fwhmmin + fwhmmax)/2.0
913  else:
914  fwhms = {}
915 
916  # -- Try to estimate the most appropriate Half-light Radius range
917  # -- Get the Half-light radii
918  nobj = 0
919  for i, fileName in enumerate(filenames):
920  fwhms[i] = []
921 
922  if prefs.getVerboseType() != prefs.QUIET:
923  print("Examining Catalog #%d" % (i+1))
924 
925  # ---- Read input catalog
926  tab = afwTable.SourceCatalog.readFits(fileName)
927 
928  # -------- Fill the FWHM array
929  shape = tab.getShapeDefinition()
930  ixx = tab.get("%s.xx" % shape)
931  iyy = tab.get("%s.yy" % shape)
932 
933  rmsSize = np.sqrt(0.5*(ixx + iyy))
934  elong = 0.5*(ixx - iyy)/(ixx + iyy)
935 
936  flux = tab.get(prefs.getPhotfluxRkey())
937  fluxErr = tab.get(prefs.getPhotfluxerrRkey())
938 
939  flags = getLsstFlags(tab)
940 
941  good = np.logical_and(flux/fluxErr > minsn,
942  np.logical_not(flags & prefs.getFlagMask()))
943  good = np.logical_and(good, elong < maxelong)
944  fwhm = 2.0*rmsSize
945  good = np.logical_and(good, fwhm >= min)
946  good = np.logical_and(good, fwhm < max)
947  fwhms[i] = fwhm[good]
948 
949  if prefs.getVarType() == prefs.VAR_NONE:
950  if nobj:
951  fwhms_all = np.empty(sum([len(l) for l in fwhms.values()]))
952  i = 0
953  for l in fwhms.values():
954  fwhms_all[i:len(l)] = l
955  i += len(l)
956  mode, min, max = compute_fwhmrange(fwhms_all, prefs.getMaxvar(),
957  prefs.getFwhmrange()[0], prefs.getFwhmrange()[1],
958  plot=plot)
959  else:
960  raise RuntimeError("No source with appropriate FWHM found!!")
961  mode = min = max = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
962 
963  fwhmmin = np.zeros(ncat) + min
964  fwhmmax = np.zeros(ncat) + max
965  fwhmmode = np.zeros(ncat) + mode
966  else:
967  fwhmmode = np.empty(ncat)
968  fwhmmin = np.empty(ncat)
969  fwhmmax = np.empty(ncat)
970 
971  for i in range(ncat):
972  nobj = len(fwhms[i])
973  if (nobj):
974  fwhmmode[i], fwhmmin[i], fwhmmax[i] = \
975  compute_fwhmrange(fwhms[i], prefs.getMaxvar(),
976  prefs.getFwhmrange()[0], prefs.getFwhmrange()[1], plot=plot)
977  else:
978  raise RuntimeError("No source with appropriate FWHM found!!")
979  fwhmmode[i] = fwhmmin[i] = fwhmmax[i] = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
980 
981  # Read the samples
982  mode = psfexLib.BIG # mode of FWHM distribution
983 
984  sets = []
985  for i, fileName in enumerate(filenames):
986  set = None
987  for ext in range(next):
988  set = read_samplesLsst(prefs, set, fileName, fwhmmin[i]/2.0, fwhmmax[i]/2.0,
989  ext, next, i, context,
990  context.getPc(i) if context.getNpc() else None, nobj, plot=plot)
991 
992  if fwhmmode[i] < mode:
993  mode = fwhmmode[i]
994 
995  set.setFwhm(mode)
996 
997  if prefs.getVerboseType() != prefs.QUIET:
998  if set.getNsample():
999  print("%d samples loaded." % set.getNsample())
1000  else:
1001  raise RuntimeError("No appropriate source found!!")
1002 
1003  sets.append(set)
1004 
1005  return sets
1006 
1007 
1008 def makeit(prefs, context, saveWcs=False, plot=dict()):
1009  """This is the python wrapper for the original psfex that reads SExtractor outputs
1010  """
1011  # Create an array of PSFs (one PSF for each extension)
1012  if prefs.getVerboseType() != prefs.QUIET:
1013  print("----- %d input catalogues:" % prefs.getNcat())
1014 
1015  if saveWcs: # only needed for making plots
1016  wcssList = []
1017 
1018  fields = psfexLib.vectorField()
1019  for cat in prefs.getCatalogs():
1020  field = psfexLib.Field(cat)
1021  wcss = []
1022  wcssList.append(wcss)
1023  with fits.open(cat) as pf:
1024  for hdu in pf:
1025  if hdu.name == "PRIMARY":
1026  pass
1027  elif hdu.name == "LDAC_IMHEAD":
1028  hdr = hdu.data[0][0] # the fits header from the original fits image
1029  md = PropertySet()
1030  for line in hdr:
1031  try:
1032  md.set(*splitFitsCard(line))
1033  except AttributeError:
1034  continue
1035 
1036  if not md.exists("CRPIX1"): # no WCS; try WCSA
1037  for k in md.names():
1038  if re.search(r"A$", k):
1039  md.set(k[:-1], md.getScalar(k))
1040  wcs = afwGeom.makeSkyWcs(md)
1041  naxis1, naxis2 = md.getScalar("NAXIS1"), md.getScalar("NAXIS2")
1042  elif hdu.name == "LDAC_OBJECTS":
1043  nobj = len(hdu.data)
1044 
1045  assert wcs, "LDAC_OBJECTS comes after LDAC_IMHEAD"
1046  field.addExt(wcs, naxis1, naxis2, nobj)
1047  if saveWcs:
1048  wcss.append((wcs, naxis1, naxis2))
1049  wcs = None
1050 
1051  field.finalize()
1052  fields.append(field)
1053 
1054  sets = psfexLib.vectorSet()
1055  for set in load_samples(prefs, context, plot=plot):
1056  sets.append(set)
1057 
1058  psfexLib.makeit(fields, sets)
1059 
1060  ret = [[f.getPsfs() for f in fields], sets]
1061  if saveWcs:
1062  ret.append(wcssList)
1063 
1064  return ret
def select_candidates(set, prefs, frmin, frmax, flags, flux, fluxerr, rmsSize, elong, vignet, plot=dict(), title="")
Definition: psfex.py:294
def read_samples(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, plot=dict(showFlags=False, showRejection=False))
Definition: psfex.py:123
def makeit(prefs, context, saveWcs=False, plot=dict())
Definition: psfex.py:1008
def showPsf(psf, set, ext=None, wcsData=None, trim=0, nspot=5, diagnostics=False, outDir="", frame=None, title=None)
Definition: psfex.py:545
def makeitLsst(prefs, context, saveWcs=False, plot=dict())
Definition: psfex.py:703
int min
def compute_fwhmrange(fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False))
Definition: psfex.py:55
int max
def mtv(data, frame=None, title="", wcs=None, args, kwargs)
Definition: ds9.py:93
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
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:496
def load_samples(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict())
Definition: psfex.py:408
Class for storing generic metadata.
Definition: PropertySet.h:68
def read_samplesLsst(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, nobj, plot=dict(showFlags=False, showRejection=False))
Definition: psfex.py:760
std::vector< SchemaItem< Flag > > * items
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
daf::base::PropertyList * list
Definition: fits.cc:885
def load_samplesLsst(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict())
Definition: psfex.py:897