LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Functions | Variables
lsst.ip.diffim.utils Namespace Reference

Functions

def showSourceSet
 
def showKernelSpatialCells
 
def showDiaSources
 
def showKernelCandidates
 
def showKernelBasis
 
def plotKernelSpatialModel
 
def showKernelMosaic
 
def plotPixelResiduals
 
def calcCentroid
 
def calcWidth
 
def printSkyDiffs
 
def makeRegions
 
def showSourceSetSky
 
def plotWhisker
 

Variables

 keptPlots = False
 

Function Documentation

def lsst.ip.diffim.utils.calcCentroid (   arr) 
Calculate first moment of a (kernel) image

Definition at line 619 of file utils.py.

619 
620 def calcCentroid(arr):
621  """Calculate first moment of a (kernel) image"""
622  y, x = arr.shape
623  sarr = arr*arr
624  xarr = np.asarray([[el for el in range(x)] for el2 in range(y)])
625  yarr = np.asarray([[el2 for el in range(x)] for el2 in range(y)])
626  narr = xarr*sarr
627  sarrSum = sarr.sum()
628  centx = narr.sum()/sarrSum
629  narr = yarr*sarr
630  centy = narr.sum()/sarrSum
631  return centx, centy
def lsst.ip.diffim.utils.calcWidth (   arr,
  centx,
  centy 
) 
Calculate second moment of a (kernel) image

Definition at line 632 of file utils.py.

632 
633 def calcWidth(arr, centx, centy):
634  """Calculate second moment of a (kernel) image"""
635  y, x = arr.shape
636  #Square the flux so we don't have to deal with negatives
637  sarr = arr*arr
638  xarr = np.asarray([[el for el in range(x)] for el2 in range(y)])
639  yarr = np.asarray([[el2 for el in range(x)] for el2 in range(y)])
640  narr = sarr*np.power((xarr - centx), 2.)
641  sarrSum = sarr.sum()
642  xstd = np.sqrt(narr.sum()/sarrSum)
643  narr = sarr*np.power((yarr - centy), 2.)
644  ystd = np.sqrt(narr.sum()/sarrSum)
645  return xstd, ystd
def lsst.ip.diffim.utils.makeRegions (   sources,
  outfilename,
  wcs = None 
) 
Create regions file for ds9 from input source list

Definition at line 658 of file utils.py.

658 
659 def makeRegions(sources, outfilename, wcs=None):
660  """Create regions file for ds9 from input source list"""
661  fh = open(outfilename, "w")
662  fh.write("global color=red font=\"helvetica 10 normal\" select=1 highlite=1 edit=1 move=1 delete=1 include=1 fixed=0 source\nfk5\n")
663  for s in sources:
664  if wcs:
665  (ra, dec) = wcs.pixelToSky(s.getCentroid().getX(), s.getCentroid().getY()).getPosition()
666  else:
667  (ra, dec) = s.getCoord().getPosition()
668  if np.isfinite(ra) and np.isfinite(dec):
669  fh.write("circle(%f,%f,2\")\n"%(ra,dec))
670  fh.flush()
671  fh.close()
def lsst.ip.diffim.utils.plotKernelSpatialModel (   kernel,
  kernelCellSet,
  showBadCandidates = True,
  numSample = 128,
  keepPlots = True,
  maxCoeff = 10 
) 
Plot the Kernel spatial model.

Definition at line 259 of file utils.py.

259 
260  numSample=128, keepPlots=True, maxCoeff = 10):
261  """Plot the Kernel spatial model."""
262 
263  try:
264  import matplotlib.pyplot as plt
265  import matplotlib.colors
266  except ImportError, e:
267  print "Unable to import numpy and matplotlib: %s" % e
268  return
269 
270  x0 = kernelCellSet.getBBox().getBeginX()
271  y0 = kernelCellSet.getBBox().getBeginY()
272 
273  candPos = list()
274  candFits = list()
275  badPos = list()
276  badFits = list()
277  candAmps = list()
278  badAmps = list()
279  for cell in kernelCellSet.getCellList():
280  for cand in cell.begin(False):
281  cand = diffimLib.cast_KernelCandidateF(cand)
282  if not showBadCandidates and cand.isBad():
283  continue
284  candCenter = afwGeom.PointD(cand.getXCenter(), cand.getYCenter())
285  try:
286  im = cand.getTemplateMaskedImage()
287  except Exception, e:
288  continue
289 
290  targetFits = badFits if cand.isBad() else candFits
291  targetPos = badPos if cand.isBad() else candPos
292  targetAmps = badAmps if cand.isBad() else candAmps
293 
294  # compare original and spatial kernel coefficients
295  kp0 = np.array(cand.getKernel(diffimLib.KernelCandidateF.ORIG).getKernelParameters())
296  amp = cand.getCandidateRating()
297 
298  targetFits = badFits if cand.isBad() else candFits
299  targetPos = badPos if cand.isBad() else candPos
300  targetAmps = badAmps if cand.isBad() else candAmps
301 
302  targetFits.append(kp0)
303  targetPos.append(candCenter)
304  targetAmps.append(amp)
305 
306  xGood = np.array([pos.getX() for pos in candPos]) - x0
307  yGood = np.array([pos.getY() for pos in candPos]) - y0
308  zGood = np.array(candFits)
309 
310  xBad = np.array([pos.getX() for pos in badPos]) - x0
311  yBad = np.array([pos.getY() for pos in badPos]) - y0
312  zBad = np.array(badFits)
313  numBad = len(badPos)
314 
315  xRange = np.linspace(0, kernelCellSet.getBBox().getWidth(), num=numSample)
316  yRange = np.linspace(0, kernelCellSet.getBBox().getHeight(), num=numSample)
317 
318  if maxCoeff:
319  maxCoeff = min(maxCoeff, kernel.getNKernelParameters())
320  else:
321  maxCoeff = kernel.getNKernelParameters()
322 
323  for k in range(maxCoeff):
324  func = kernel.getSpatialFunction(k)
325  dfGood = zGood[:,k] - np.array([func(pos.getX(), pos.getY()) for pos in candPos])
326  yMin = dfGood.min()
327  yMax = dfGood.max()
328  if numBad > 0:
329  dfBad = zBad[:,k] - np.array([func(pos.getX(), pos.getY()) for pos in badPos])
330  # Can really screw up the range...
331  yMin = min([yMin, dfBad.min()])
332  yMax = max([yMax, dfBad.max()])
333  yMin -= 0.05 * (yMax - yMin)
334  yMax += 0.05 * (yMax - yMin)
335 
336  fRange = np.ndarray((len(xRange), len(yRange)))
337  for j, yVal in enumerate(yRange):
338  for i, xVal in enumerate(xRange):
339  fRange[j][i] = func(xVal, yVal)
340 
341  fig = plt.figure(k)
342 
343  fig.clf()
344  try:
345  fig.canvas._tkcanvas._root().lift() # == Tk's raise, but raise is a python reserved word
346  except Exception: # protect against API changes
347  pass
348 
349  fig.suptitle('Kernel component %d' % k)
350 
351  # LL
352  ax = fig.add_axes((0.1, 0.05, 0.35, 0.35))
353  vmin = fRange.min() # - 0.05 * np.fabs(fRange.min())
354  vmax = fRange.max() # + 0.05 * np.fabs(fRange.max())
355  norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
356  im = ax.imshow(fRange, aspect='auto', norm=norm,
357  extent=[0, kernelCellSet.getBBox().getWidth()-1,
358  0, kernelCellSet.getBBox().getHeight()-1])
359  ax.set_title('Spatial polynomial')
360  plt.colorbar(im, orientation='horizontal', ticks=[vmin, vmax])
361 
362  # UL
363  ax = fig.add_axes((0.1, 0.55, 0.35, 0.35))
364  ax.plot(-2.5*np.log10(candAmps), zGood[:,k], 'b+')
365  if numBad > 0:
366  ax.plot(-2.5*np.log10(badAmps), zBad[:,k], 'r+')
367  ax.set_title("Basis Coefficients")
368  ax.set_xlabel("Instr mag")
369  ax.set_ylabel("Coeff")
370 
371  # LR
372  ax = fig.add_axes((0.55, 0.05, 0.35, 0.35))
373  ax.set_autoscale_on(False)
374  ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getHeight())
375  ax.set_ybound(lower=yMin, upper=yMax)
376  ax.plot(yGood, dfGood, 'b+')
377  if numBad > 0:
378  ax.plot(yBad, dfBad, 'r+')
379  ax.axhline(0.0)
380  ax.set_title('dCoeff (indiv-spatial) vs. y')
381 
382  # UR
383  ax = fig.add_axes((0.55, 0.55, 0.35, 0.35))
384  ax.set_autoscale_on(False)
385  ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getWidth())
386  ax.set_ybound(lower=yMin, upper=yMax)
387  ax.plot(xGood, dfGood, 'b+')
388  if numBad > 0:
389  ax.plot(xBad, dfBad, 'r+')
390  ax.axhline(0.0)
391  ax.set_title('dCoeff (indiv-spatial) vs. x')
392 
393  fig.show()
394 
395  global keptPlots
396  if keepPlots and not keptPlots:
397  # Keep plots open when done
398  def show():
399  print "%s: Please close plots when done." % __name__
400  try:
401  plt.show()
402  except Exception:
403  pass
404  print "Plots closed, exiting..."
405  import atexit
406  atexit.register(show)
407  keptPlots = True
408 
def lsst.ip.diffim.utils.plotPixelResiduals (   exposure,
  warpedTemplateExposure,
  diffExposure,
  kernelCellSet,
  kernel,
  background,
  testSources,
  config,
  origVariance = False,
  nptsFull = 1e6,
  keepPlots = True,
  titleFs = 14 
) 
Plot diffim residuals for LOCAL and SPATIAL models

Definition at line 478 of file utils.py.

478 
479  origVariance = False, nptsFull = 1e6, keepPlots = True, titleFs=14):
480  """Plot diffim residuals for LOCAL and SPATIAL models"""
481  candidateResids = []
482  spatialResids = []
483  nonfitResids = []
484 
485  for cell in kernelCellSet.getCellList():
486  for cand in cell.begin(True): # only look at good ones
487  # Be sure
488  if not (cand.getStatus() == afwMath.SpatialCellCandidate.GOOD):
489  continue
490 
491  cand = diffimLib.cast_KernelCandidateF(cand)
492  diffim = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
493  orig = cand.getScienceMaskedImage()
494 
495  ski = afwImage.ImageD(kernel.getDimensions())
496  kernel.computeImage(ski, False, int(cand.getXCenter()), int(cand.getYCenter()))
497  sk = afwMath.FixedKernel(ski)
498  sbg = background(int(cand.getXCenter()), int(cand.getYCenter()))
499  sdiffim = cand.getDifferenceImage(sk, sbg)
500 
501  # trim edgs due to convolution
502  bbox = kernel.shrinkBBox(diffim.getBBox())
503  tdiffim = diffim.Factory(diffim, bbox)
504  torig = orig.Factory(orig, bbox)
505  tsdiffim = sdiffim.Factory(sdiffim, bbox)
506 
507  if origVariance:
508  candidateResids.append(np.ravel(tdiffim.getImage().getArray() /
509  np.sqrt(torig.getVariance().getArray())))
510  spatialResids.append(np.ravel(tsdiffim.getImage().getArray() /
511  np.sqrt(torig.getVariance().getArray())))
512  else:
513  candidateResids.append(np.ravel(tdiffim.getImage().getArray() /
514  np.sqrt(tdiffim.getVariance().getArray())))
515  spatialResids.append(np.ravel(tsdiffim.getImage().getArray() /
516  np.sqrt(tsdiffim.getVariance().getArray())))
517 
518  fullIm = diffExposure.getMaskedImage().getImage().getArray()
519  fullMask = diffExposure.getMaskedImage().getMask().getArray()
520  if origVariance:
521  fullVar = exposure.getMaskedImage().getVariance().getArray()
522  else:
523  fullVar = diffExposure.getMaskedImage().getVariance().getArray()
524 
525  bitmaskBad = 0
526  bitmaskBad |= afwImage.MaskU.getPlaneBitMask('NO_DATA')
527  bitmaskBad |= afwImage.MaskU.getPlaneBitMask('SAT')
528  idx = np.where((fullMask & bitmaskBad) == 0)
529  stride = int(len(idx[0]) // nptsFull)
530  sidx = idx[0][::stride], idx[1][::stride]
531  allResids = fullIm[sidx] / np.sqrt(fullVar[sidx])
532 
533  testFootprints = diffimTools.sourceToFootprintList(testSources, warpedTemplateExposure,
534  exposure, config, pexLog.getDefaultLog())
535  for fp in testFootprints:
536  subexp = diffExposure.Factory(diffExposure, fp["footprint"].getBBox())
537  subim = subexp.getMaskedImage().getImage()
538  if origVariance:
539  subvar = afwImage.ExposureF(exposure, fp["footprint"].getBBox()).getMaskedImage().getVariance()
540  else:
541  subvar = subexp.getMaskedImage().getVariance()
542  nonfitResids.append(np.ravel(subim.getArray() / np.sqrt(subvar.getArray())))
543 
544  candidateResids = np.ravel(np.array(candidateResids))
545  spatialResids = np.ravel(np.array(spatialResids))
546  nonfitResids = np.ravel(np.array(nonfitResids))
547 
548  try:
549  import pylab
550  from matplotlib.font_manager import FontProperties
551  except ImportError, e:
552  print "Unable to import pylab: %s" % e
553  return
554 
555  fig = pylab.figure()
556  fig.clf()
557  try:
558  fig.canvas._tkcanvas._root().lift() # == Tk's raise, but raise is a python reserved word
559  except Exception: # protect against API changes
560  pass
561  if origVariance:
562  fig.suptitle("Diffim residuals: Normalized by sqrt(input variance)", fontsize=titleFs)
563  else:
564  fig.suptitle("Diffim residuals: Normalized by sqrt(diffim variance)", fontsize=titleFs)
565 
566  sp1 = pylab.subplot(221)
567  sp2 = pylab.subplot(222, sharex=sp1, sharey=sp1)
568  sp3 = pylab.subplot(223, sharex=sp1, sharey=sp1)
569  sp4 = pylab.subplot(224, sharex=sp1, sharey=sp1)
570  xs = np.arange(-5, 5.05, 0.1)
571  ys = 1. / np.sqrt(2 * np.pi) * np.exp( -0.5 * xs**2 )
572 
573  sp1.hist(candidateResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
574  % (np.mean(candidateResids), np.var(candidateResids)))
575  sp1.plot(xs, ys, "r-", lw=2, label="N(0,1)")
576  sp1.set_title("Candidates: basis fit", fontsize=titleFs-2)
577  sp1.legend(loc=1, fancybox=True, shadow=True, prop = FontProperties(size=titleFs-6))
578 
579  sp2.hist(spatialResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
580  % (np.mean(spatialResids), np.var(spatialResids)))
581  sp2.plot(xs, ys, "r-", lw=2, label="N(0,1)")
582  sp2.set_title("Candidates: spatial fit", fontsize=titleFs-2)
583  sp2.legend(loc=1, fancybox=True, shadow=True, prop = FontProperties(size=titleFs-6))
584 
585  sp3.hist(nonfitResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
586  % (np.mean(nonfitResids), np.var(nonfitResids)))
587  sp3.plot(xs, ys, "r-", lw=2, label="N(0,1)")
588  sp3.set_title("Control sample: spatial fit", fontsize=titleFs-2)
589  sp3.legend(loc=1, fancybox=True, shadow=True, prop = FontProperties(size=titleFs-6))
590 
591  sp4.hist(allResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
592  % (np.mean(allResids), np.var(allResids)))
593  sp4.plot(xs, ys, "r-", lw=2, label="N(0,1)")
594  sp4.set_title("Full image (subsampled)", fontsize=titleFs-2)
595  sp4.legend(loc=1, fancybox=True, shadow=True, prop = FontProperties(size=titleFs-6))
596 
597  pylab.setp(sp1.get_xticklabels()+sp1.get_yticklabels(), fontsize=titleFs-4)
598  pylab.setp(sp2.get_xticklabels()+sp2.get_yticklabels(), fontsize=titleFs-4)
599  pylab.setp(sp3.get_xticklabels()+sp3.get_yticklabels(), fontsize=titleFs-4)
600  pylab.setp(sp4.get_xticklabels()+sp4.get_yticklabels(), fontsize=titleFs-4)
601 
602  sp1.set_xlim(-5, 5)
603  sp1.set_ylim(0, 0.5)
604  fig.show()
605 
606  global keptPlots
607  if keepPlots and not keptPlots:
608  # Keep plots open when done
609  def show():
610  print "%s: Please close plots when done." % __name__
611  try:
612  pylab.show()
613  except Exception:
614  pass
615  print "Plots closed, exiting..."
616  import atexit
617  atexit.register(show)
618  keptPlots = True
lsst.afw.math::FixedKernel
A kernel created from an Image.
Definition: Kernel.h:551
def lsst.ip.diffim.utils.plotWhisker (   results,
  newWcs 
) 
Plot whisker diagram of astromeric offsets between results.matches

Definition at line 681 of file utils.py.

681 
682 def plotWhisker(results, newWcs):
683  """Plot whisker diagram of astromeric offsets between results.matches"""
684  refCoordKey = results.matches[0].first.getTable().getCoordKey()
685  inCentroidKey = results.matches[0].second.getTable().getCentroidKey()
686  positions = [m.first.get(refCoordKey) for m in results.matches]
687  residuals = [m.first.get(refCoordKey).getOffsetFrom(
688  newWcs.pixelToSky(m.second.get(inCentroidKey))) for
689  m in results.matches]
690  import matplotlib.pyplot as plt
691  fig = plt.figure()
692  sp = fig.add_subplot(1, 1, 0)
693  xpos = [x[0].asDegrees() for x in positions]
694  ypos = [x[1].asDegrees() for x in positions]
695  xpos.append(0.02*(max(xpos) - min(xpos)) + min(xpos))
696  ypos.append(0.98*(max(ypos) - min(ypos)) + min(ypos))
697  xidxs = np.isfinite(xpos)
698  yidxs = np.isfinite(ypos)
699  X = np.asarray(xpos)[xidxs]
700  Y = np.asarray(ypos)[yidxs]
701  distance = [x[1].asArcseconds() for x in residuals]
702  distance.append(0.2)
703  distance = np.asarray(distance)[xidxs]
704  #NOTE: This assumes that the bearing is measured positive from +RA through North.
705  #From the documentation this is not clear.
706  bearing = [x[0].asRadians() for x in residuals]
707  bearing.append(0)
708  bearing = np.asarray(bearing)[xidxs]
709  U = (distance*np.cos(bearing))
710  V = (distance*np.sin(bearing))
711  sp.quiver(X, Y, U, V)
712  sp.set_title("WCS Residual")
713  plt.show()
714 
def lsst.ip.diffim.utils.printSkyDiffs (   sources,
  wcs 
) 
Print differences in sky coordinates between source Position and its Centroid mapped through Wcs

Definition at line 646 of file utils.py.

646 
647 def printSkyDiffs(sources, wcs):
648  """Print differences in sky coordinates between source Position and its Centroid mapped through Wcs"""
649  for s in sources:
650  sCentroid = s.getCentroid()
651  sPosition = s.getCoord().getPosition()
652  dra = 3600*(sPosition.getX() - wcs.pixelToSky(sCentroid.getX(),
653  sCentroid.getY()).getPosition().getX())/0.2
654  ddec = 3600*(sPosition.getY() - wcs.pixelToSky(sCentroid.getX(),
655  sCentroid.getY()).getPosition().getY())/0.2
656  if np.isfinite(dra) and np.isfinite(ddec):
657  print dra, ddec
def lsst.ip.diffim.utils.showDiaSources (   sources,
  exposure,
  isFlagged,
  isDipole,
  frame = None 
) 
Display Dia Sources

Definition at line 92 of file utils.py.

92 
93 def showDiaSources(sources, exposure, isFlagged, isDipole, frame=None):
94  """Display Dia Sources
95  """
96  #
97  # Show us the ccandidates
98  #
99  # Too many mask planes in diffims
100  for plane in ("BAD", "CR", "EDGE", "INTERPOlATED", "INTRP", "SAT", "SATURATED"):
101  ds9.setMaskPlaneVisibility(plane, False)
102 
103  mos = displayUtils.Mosaic()
104  for i in range(len(sources)):
105  source = sources[i]
106  badFlag = isFlagged[i]
107  dipoleFlag = isDipole[i]
108  bbox = source.getFootprint().getBBox()
109  stamp = exposure.Factory(exposure, bbox, True)
110  im = displayUtils.Mosaic(gutter=1, background=0, mode="x")
111  im.append(stamp.getMaskedImage())
112  lab = "%.1f,%.1f:" % (source.getX(), source.getY())
113  if badFlag:
114  ctype = ds9.RED
115  lab += "BAD"
116  if dipoleFlag:
117  ctype = ds9.YELLOW
118  lab += "DIPOLE"
119  if not badFlag and not dipoleFlag:
120  ctype = ds9.GREEN
121  lab += "OK"
122  mos.append(im.makeMosaic(), lab, ctype)
123  #print source.getId()
124  title = "Dia Sources"
125  mosaicImage = mos.makeMosaic(frame=frame, title=title)
126  return mosaicImage
def lsst.ip.diffim.utils.showKernelBasis (   kernel,
  frame = None 
) 
Display a Kernel's basis images

Definition at line 243 of file utils.py.

243 
244 def showKernelBasis(kernel, frame=None):
245  """Display a Kernel's basis images
246  """
247  mos = displayUtils.Mosaic()
248 
249  for k in kernel.getKernelList():
250  im = afwImage.ImageD(k.getDimensions())
251  k.computeImage(im, False)
252  mos.append(im)
253  mos.makeMosaic(frame=frame, title="Kernel Basis Images")
254 
255  return mos
def lsst.ip.diffim.utils.showKernelCandidates (   kernelCellSet,
  kernel,
  background,
  frame = None,
  showBadCandidates = True,
  resids = False,
  kernels = False 
) 
Display the Kernel candidates.
If kernel is provided include spatial model and residuals;  
If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi

Definition at line 128 of file utils.py.

128 
129  resids=False, kernels=False):
130  """Display the Kernel candidates.
131  If kernel is provided include spatial model and residuals;
132  If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi
133  """
134 
135  #
136  # Show us the ccandidates
137  #
138  if kernels:
139  mos = displayUtils.Mosaic(gutter=5, background=0)
140  else:
141  mos = displayUtils.Mosaic(gutter=5, background=-1)
142  #
143  candidateCenters = []
144  candidateCentersBad = []
145  candidateIndex = 0
146  for cell in kernelCellSet.getCellList():
147  for cand in cell.begin(False): # include bad candidates
148  cand = diffimLib.cast_KernelCandidateF(cand)
149 
150  # Original difference image; if does not exist, skip candidate
151  try:
152  resid = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
153  except Exception:
154  continue
155 
156  rchi2 = cand.getChi2()
157  if rchi2 > 1e100:
158  rchi2 = np.nan
159 
160  if not showBadCandidates and cand.isBad():
161  continue
162 
163  im_resid = displayUtils.Mosaic(gutter=1, background=-0.5, mode="x")
164 
165  try:
166  im = cand.getScienceMaskedImage()
167  im = im.Factory(im, True)
168  im.setXY0(cand.getScienceMaskedImage().getXY0())
169  except Exception:
170  continue
171  if (not resids and not kernels):
172  im_resid.append(im.Factory(im, True))
173  try:
174  im = cand.getTemplateMaskedImage()
175  im = im.Factory(im, True)
176  im.setXY0(cand.getTemplateMaskedImage().getXY0())
177  except Exception:
178  continue
179  if (not resids and not kernels):
180  im_resid.append(im.Factory(im, True))
181 
182  # Difference image with original basis
183  if resids:
184  var = resid.getVariance()
185  var = var.Factory(var, True)
186  np.sqrt(var.getArray(), var.getArray()) # inplace sqrt
187  resid = resid.getImage()
188  resid /= var
189  bbox = kernel.shrinkBBox(resid.getBBox())
190  resid = resid.Factory(resid, bbox, True)
191  elif kernels:
192  kim = cand.getKernelImage(diffimLib.KernelCandidateF.ORIG).convertF()
193  resid = kim.Factory(kim, True)
194  im_resid.append(resid)
195 
196  # residuals using spatial model
197  ski = afwImage.ImageD(kernel.getDimensions())
198  kernel.computeImage(ski, False, int(cand.getXCenter()), int(cand.getYCenter()))
199  sk = afwMath.FixedKernel(ski)
200  sbg = 0.0
201  if background:
202  sbg = background(int(cand.getXCenter()), int(cand.getYCenter()))
203  sresid = cand.getDifferenceImage(sk, sbg)
204  resid = sresid
205  if resids:
206  resid = sresid.getImage()
207  resid /= var
208  bbox = kernel.shrinkBBox(resid.getBBox())
209  resid = resid.Factory(resid, bbox, True)
210  elif kernels:
211  kim = ski.convertF()
212  resid = kim.Factory(kim, True)
213  im_resid.append(resid)
214 
215  im = im_resid.makeMosaic()
216 
217  lab = "%d chi^2 %.1f" % (cand.getId(), rchi2)
218  ctype = ds9.RED if cand.isBad() else ds9.GREEN
219 
220  mos.append(im, lab, ctype)
221 
222  if False and np.isnan(rchi2):
223  ds9.mtv(cand.getScienceMaskedImage.getImage(), title="candidate", frame=1)
224  print "rating", cand.getCandidateRating()
225 
226  im = cand.getScienceMaskedImage()
227  center = (candidateIndex, cand.getXCenter() - im.getX0(), cand.getYCenter() - im.getY0())
228  candidateIndex += 1
229  if cand.isBad():
230  candidateCentersBad.append(center)
231  else:
232  candidateCenters.append(center)
233 
234  if resids:
235  title = "chi Diffim"
236  elif kernels:
237  title = "Kernels"
238  else:
239  title = "Candidates & residuals"
240  mosaicImage = mos.makeMosaic(frame=frame, title=title)
241 
242  return mosaicImage
lsst.afw.math::FixedKernel
A kernel created from an Image.
Definition: Kernel.h:551
def lsst.ip.diffim.utils.showKernelMosaic (   bbox,
  kernel,
  nx = 7,
  ny = None,
  frame = None,
  title = None,
  showCenter = True,
  showEllipticity = True 
) 
Show a mosaic of Kernel images.

Definition at line 410 of file utils.py.

410 
411  showCenter=True, showEllipticity=True):
412  """Show a mosaic of Kernel images.
413  """
414  mos = displayUtils.Mosaic()
415 
416  x0 = bbox.getBeginX()
417  y0 = bbox.getBeginY()
418  width = bbox.getWidth()
419  height = bbox.getHeight()
420 
421  if not ny:
422  ny = int(nx*float(height)/width + 0.5)
423  if not ny:
424  ny = 1
425 
426  schema = afwTable.SourceTable.makeMinimalSchema()
427  control = measAlg.GaussianCentroidControl()
428  centroider = measAlg.MeasureSourcesBuilder().addAlgorithm(control).build(schema)
429  sdssShape = measAlg.SdssShapeControl()
430  shaper = measAlg.MeasureSourcesBuilder().addAlgorithm(sdssShape).build(schema)
431  table = afwTable.SourceTable.make(schema)
432  table.defineCentroid(control.name)
433  table.defineShape(sdssShape.name)
434 
435  centers = []
436  shapes = []
437  for iy in range(ny):
438  for ix in range(nx):
439  x = int(ix*(width-1)/(nx-1)) + x0
440  y = int(iy*(height-1)/(ny-1)) + y0
441 
442  im = afwImage.ImageD(kernel.getDimensions())
443  ksum = kernel.computeImage(im, False, x, y)
444  lab = "Kernel(%d,%d)=%.2f" % (x, y, ksum) if False else ""
445  mos.append(im, lab)
446 
447  exp = afwImage.makeExposure(afwImage.makeMaskedImage(im))
448  w, h = im.getWidth(), im.getHeight()
449  cen = afwGeom.PointD(w//2, h//2)
450  src = table.makeRecord()
451  foot = afwDet.Footprint(exp.getBBox())
452  src.setFootprint(foot)
453 
454  centroider.apply(src, exp, cen)
455  centers.append((src.getX(), src.getY()))
456 
457  shaper.apply(src, exp, cen)
458  shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
459 
460  mos.makeMosaic(frame=frame, title=title if title else "Model Kernel", mode=nx)
461 
462  if centers and frame is not None:
463  i = 0
464  with ds9.Buffering():
465  for cen, shape in zip(centers, shapes):
466  bbox = mos.getBBox(i); i += 1
467  xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
468  if showCenter:
469  ds9.dot("+", xc, yc, ctype=ds9.BLUE, frame=frame)
470 
471  if showEllipticity:
472  ixx, ixy, iyy = shape
473  ds9.dot("@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, frame=frame, ctype=ds9.RED)
474 
475  return mos
lsst.afw.image::makeMaskedImage
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename Image< ImagePixelT >::Ptr image, typename Mask< MaskPixelT >::Ptr mask=typename Mask< MaskPixelT >::Ptr(), typename Image< VariancePixelT >::Ptr variance=typename Image< VariancePixelT >::Ptr())
Definition: MaskedImage.h:1067
lsst.afw.detection::Footprint
A set of pixels in an Image.
Definition: Footprint.h:62
lsst.afw.image::makeExposure
Exposure< ImagePixelT, MaskPixelT, VariancePixelT >::Ptr makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, boost::shared_ptr< Wcs const > wcs=boost::shared_ptr< Wcs const >())
Definition: Exposure.h:308
def lsst.ip.diffim.utils.showKernelSpatialCells (   maskedIm,
  kernelCellSet,
  showChi2 = False,
  symb = "o",
  ctype = None,
  ctypeUnused = None,
  ctypeBad = None,
  size = 3,
  frame = None,
  title = "Spatial Cells" 
) 
Show the SpatialCells.  If symb is something that ds9.dot
understands (e.g. "o"), the top nMaxPerCell candidates will be
indicated with that symbol, using ctype and size

Definition at line 57 of file utils.py.

57 
58  frame=None, title="Spatial Cells"):
59  """Show the SpatialCells. If symb is something that ds9.dot
60  understands (e.g. "o"), the top nMaxPerCell candidates will be
61  indicated with that symbol, using ctype and size"""
62 
63  ds9.mtv(maskedIm, frame=frame, title=title)
64  with ds9.Buffering():
65  origin = [-maskedIm.getX0(), -maskedIm.getY0()]
66  for cell in kernelCellSet.getCellList():
67  displayUtils.drawBBox(cell.getBBox(), origin=origin, frame=frame)
68 
69  goodies = ctypeBad is None
70  for cand in cell.begin(goodies):
71  cand = diffimLib.cast_KernelCandidateF(cand)
72  xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
73  if cand.getStatus() == afwMath.SpatialCellCandidate.BAD:
74  color = ctypeBad
75  elif cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
76  color = ctype
77  elif cand.getStatus() == afwMath.SpatialCellCandidate.UNKNOWN:
78  color = ctypeUnused
79  else:
80  continue
81 
82  if color:
83  ds9.dot(symb, xc, yc, frame=frame, ctype=color, size=size)
84 
85  if showChi2:
86  rchi2 = cand.getChi2()
87  if rchi2 > 1e100:
88  rchi2 = np.nan
89  ds9.dot("%d %.1f" % (cand.getId(), rchi2),
90  xc - size, yc - size - 4, frame=frame, ctype=color, size=size)
91 
def lsst.ip.diffim.utils.showSourceSet (   sSet,
  xy0 = (0, 0,
  frame = 0,
  ctype = ds9.GREEN,
  symb = "+",
  size = 2 
) 
Draw the (XAstrom, YAstrom) positions of a set of Sources.  Image has the given XY0

Definition at line 39 of file utils.py.

39 
40 def showSourceSet(sSet, xy0=(0, 0), frame=0, ctype=ds9.GREEN, symb="+", size=2):
41  """Draw the (XAstrom, YAstrom) positions of a set of Sources. Image has the given XY0"""
42 
43  with ds9.Buffering():
44  for s in sSet:
45  xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
46 
47  if symb == "id":
48  ds9.dot(str(s.getId()), xc, yc, frame=frame, ctype=ctype, size=size)
49  else:
50  ds9.dot(symb, xc, yc, frame=frame, ctype=ctype, size=size)
51 
52 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
53 #
54 # Kernel display utilities
#
def lsst.ip.diffim.utils.showSourceSetSky (   sSet,
  wcs,
  xy0,
  frame = 0,
  ctype = ds9.GREEN,
  symb = "+",
  size = 2 
) 
Draw the (RA, Dec) positions of a set of Sources. Image has the XY0.

Definition at line 672 of file utils.py.

672 
673 def showSourceSetSky(sSet, wcs, xy0, frame=0, ctype=ds9.GREEN, symb="+", size=2):
674  """Draw the (RA, Dec) positions of a set of Sources. Image has the XY0."""
675  with ds9.Buffering():
676  for s in sSet:
677  (xc, yc) = wcs.skyToPixel(s.getCoord().getRa(), s.getCoord().getDec())
678  xc -= xy0[0]
679  yc -= xy0[1]
680  ds9.dot(symb, xc, yc, frame=frame, ctype=ctype, size=size)

Variable Documentation

lsst.ip.diffim.utils.keptPlots = False

Definition at line 37 of file utils.py.

Generated on Wed Sep 16 2015 13:35:49 for LSSTApplications by   doxygen 1.8.5