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

Classes

class  DipoleTestImage
 

Functions

def showSourceSet (sSet, xy0=(0, 0), frame=0, ctype=afwDisplay.GREEN, symb="+", size=2)
 
def showKernelSpatialCells (maskedIm, kernelCellSet, showChi2=False, symb="o", ctype=None, ctypeUnused=None, ctypeBad=None, size=3, frame=None, title="Spatial Cells")
 
def showDiaSources (sources, exposure, isFlagged, isDipole, frame=None)
 
def showKernelCandidates (kernelCellSet, kernel, background, frame=None, showBadCandidates=True, resids=False, kernels=False)
 
def showKernelBasis (kernel, frame=None)
 
def plotKernelSpatialModel (kernel, kernelCellSet, showBadCandidates=True, numSample=128, keepPlots=True, maxCoeff=10)
 
def plotKernelCoefficients (spatialKernel, kernelCellSet, showBadCandidates=False, keepPlots=True)
 
def showKernelMosaic (bbox, kernel, nx=7, ny=None, frame=None, title=None, showCenter=True, showEllipticity=True)
 
def plotPixelResiduals (exposure, warpedTemplateExposure, diffExposure, kernelCellSet, kernel, background, testSources, config, origVariance=False, nptsFull=1e6, keepPlots=True, titleFs=14)
 
def calcCentroid (arr)
 
def calcWidth (arr, centx, centy)
 
def printSkyDiffs (sources, wcs)
 
def makeRegions (sources, outfilename, wcs=None)
 
def showSourceSetSky (sSet, wcs, xy0, frame=0, ctype=afwDisplay.GREEN, symb="+", size=2)
 
def plotWhisker (results, newWcs)
 

Variables

bool keptPlots = False
 

Function Documentation

◆ calcCentroid()

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

Definition at line 814 of file utils.py.

814 def calcCentroid(arr):
815  """Calculate first moment of a (kernel) image.
816  """
817  y, x = arr.shape
818  sarr = arr*arr
819  xarr = np.asarray([[el for el in range(x)] for el2 in range(y)])
820  yarr = np.asarray([[el2 for el in range(x)] for el2 in range(y)])
821  narr = xarr*sarr
822  sarrSum = sarr.sum()
823  centx = narr.sum()/sarrSum
824  narr = yarr*sarr
825  centy = narr.sum()/sarrSum
826  return centx, centy
827 
828 
def calcCentroid(arr)
Definition: utils.py:814

◆ calcWidth()

def lsst.ip.diffim.utils.calcWidth (   arr,
  centx,
  centy 
)
Calculate second moment of a (kernel) image.

Definition at line 829 of file utils.py.

829 def calcWidth(arr, centx, centy):
830  """Calculate second moment of a (kernel) image.
831  """
832  y, x = arr.shape
833  # Square the flux so we don't have to deal with negatives
834  sarr = arr*arr
835  xarr = np.asarray([[el for el in range(x)] for el2 in range(y)])
836  yarr = np.asarray([[el2 for el in range(x)] for el2 in range(y)])
837  narr = sarr*np.power((xarr - centx), 2.)
838  sarrSum = sarr.sum()
839  xstd = np.sqrt(narr.sum()/sarrSum)
840  narr = sarr*np.power((yarr - centy), 2.)
841  ystd = np.sqrt(narr.sum()/sarrSum)
842  return xstd, ystd
843 
844 
def calcWidth(arr, centx, centy)
Definition: utils.py:829

◆ makeRegions()

def lsst.ip.diffim.utils.makeRegions (   sources,
  outfilename,
  wcs = None 
)
Create regions file for display from input source list.

Definition at line 860 of file utils.py.

860 def makeRegions(sources, outfilename, wcs=None):
861  """Create regions file for display from input source list.
862  """
863  fh = open(outfilename, "w")
864  fh.write("global color=red font=\"helvetica 10 normal\" "
865  "select=1 highlite=1 edit=1 move=1 delete=1 include=1 fixed=0 source\nfk5\n")
866  for s in sources:
867  if wcs:
868  (ra, dec) = wcs.pixelToSky(s.getCentroid()).getPosition(geom.degrees)
869  else:
870  (ra, dec) = s.getCoord().getPosition(geom.degrees)
871  if np.isfinite(ra) and np.isfinite(dec):
872  fh.write("circle(%f,%f,2\")\n"%(ra, dec))
873  fh.flush()
874  fh.close()
875 
876 
def makeRegions(sources, outfilename, wcs=None)
Definition: utils.py:860

◆ plotKernelCoefficients()

def lsst.ip.diffim.utils.plotKernelCoefficients (   spatialKernel,
  kernelCellSet,
  showBadCandidates = False,
  keepPlots = True 
)
Plot the individual kernel candidate and the spatial kernel solution coefficients.

Parameters
----------

spatialKernel : `lsst.afw.math.LinearCombinationKernel`
    The spatial spatialKernel solution model which is a spatially varying linear combination
    of the spatialKernel basis functions.
    Typically returned by `lsst.ip.diffim.SpatialKernelSolution.getSolutionPair()`.

kernelCellSet : `lsst.afw.math.SpatialCellSet`
    The spatial cells that was used for solution for the spatialKernel. They contain the
    local solutions of the AL kernel for the selected sources.

showBadCandidates : `bool`, optional
    If True, plot the coefficient values for kernel candidates where the solution was marked
    bad by the numerical algorithm. Defaults to False.

keepPlots: `bool`, optional
    If True, sets ``plt.show()`` to be called before the task terminates, so that the plots
    can be explored interactively. Defaults to True.

Notes
-----
This function produces 3 figures per image subtraction operation.
* A grid plot of the local solutions. Each grid cell corresponds to a proportional area in
  the image. In each cell, local kernel solution coefficients are plotted of kernel candidates (color)
  that fall into this area as a function of the kernel basis function number.
* A grid plot of the spatial solution. Each grid cell corresponds to a proportional area in
  the image. In each cell, the spatial solution coefficients are evaluated for the center of the cell.
* Histogram of the local solution coefficients. Red line marks the spatial solution value at
  center of the image.

This function is called if ``lsst.ip.diffim.psfMatch.plotKernelCoefficients==True`` in lsstDebug. This
function was implemented as part of DM-17825.

Definition at line 428 of file utils.py.

428 def plotKernelCoefficients(spatialKernel, kernelCellSet, showBadCandidates=False, keepPlots=True):
429  """Plot the individual kernel candidate and the spatial kernel solution coefficients.
430 
431  Parameters
432  ----------
433 
434  spatialKernel : `lsst.afw.math.LinearCombinationKernel`
435  The spatial spatialKernel solution model which is a spatially varying linear combination
436  of the spatialKernel basis functions.
437  Typically returned by `lsst.ip.diffim.SpatialKernelSolution.getSolutionPair()`.
438 
439  kernelCellSet : `lsst.afw.math.SpatialCellSet`
440  The spatial cells that was used for solution for the spatialKernel. They contain the
441  local solutions of the AL kernel for the selected sources.
442 
443  showBadCandidates : `bool`, optional
444  If True, plot the coefficient values for kernel candidates where the solution was marked
445  bad by the numerical algorithm. Defaults to False.
446 
447  keepPlots: `bool`, optional
448  If True, sets ``plt.show()`` to be called before the task terminates, so that the plots
449  can be explored interactively. Defaults to True.
450 
451  Notes
452  -----
453  This function produces 3 figures per image subtraction operation.
454  * A grid plot of the local solutions. Each grid cell corresponds to a proportional area in
455  the image. In each cell, local kernel solution coefficients are plotted of kernel candidates (color)
456  that fall into this area as a function of the kernel basis function number.
457  * A grid plot of the spatial solution. Each grid cell corresponds to a proportional area in
458  the image. In each cell, the spatial solution coefficients are evaluated for the center of the cell.
459  * Histogram of the local solution coefficients. Red line marks the spatial solution value at
460  center of the image.
461 
462  This function is called if ``lsst.ip.diffim.psfMatch.plotKernelCoefficients==True`` in lsstDebug. This
463  function was implemented as part of DM-17825.
464  """
465  try:
466  import matplotlib.pyplot as plt
467  except ImportError as e:
468  print("Unable to import matplotlib: %s" % e)
469  return
470 
471  # Image dimensions
472  imgBBox = kernelCellSet.getBBox()
473  x0 = imgBBox.getBeginX()
474  y0 = imgBBox.getBeginY()
475  wImage = imgBBox.getWidth()
476  hImage = imgBBox.getHeight()
477  imgCenterX = imgBBox.getCenterX()
478  imgCenterY = imgBBox.getCenterY()
479 
480  # Plot the local solutions
481  # ----
482 
483  # Grid size
484  nX = 8
485  nY = 8
486  wCell = wImage / nX
487  hCell = hImage / nY
488 
489  fig = plt.figure()
490  fig.suptitle("Kernel candidate parameters on an image grid")
491  arrAx = fig.subplots(nrows=nY, ncols=nX, sharex=True, sharey=True, gridspec_kw=dict(
492  wspace=0, hspace=0))
493 
494  # Bottom left panel is for bottom left part of the image
495  arrAx = arrAx[::-1, :]
496 
497  allParams = []
498  for cell in kernelCellSet.getCellList():
499  cellBBox = afwGeom.Box2D(cell.getBBox())
500  # Determine which panel this spatial cell belongs to
501  iX = int((cellBBox.getCenterX() - x0)//wCell)
502  iY = int((cellBBox.getCenterY() - y0)//hCell)
503 
504  for cand in cell.begin(False):
505  try:
506  kernel = cand.getKernel(cand.ORIG)
507  except Exception:
508  continue
509 
510  if not showBadCandidates and cand.isBad():
511  continue
512 
513  nKernelParams = kernel.getNKernelParameters()
514  kernelParams = np.array(kernel.getKernelParameters())
515  allParams.append(kernelParams)
516 
517  if cand.isBad():
518  color = 'red'
519  else:
520  color = None
521  arrAx[iY, iX].plot(np.arange(nKernelParams), kernelParams, '.-',
522  color=color, drawstyle='steps-mid', linewidth=0.1)
523  for ax in arrAx.ravel():
524  ax.grid(True, axis='y')
525 
526  # Plot histogram of the local parameters and the global solution at the image center
527  # ----
528 
529  spatialFuncs = spatialKernel.getSpatialFunctionList()
530  nKernelParams = spatialKernel.getNKernelParameters()
531  nX = 8
532  fig = plt.figure()
533  fig.suptitle("Hist. of parameters marked with spatial solution at img center")
534  arrAx = fig.subplots(nrows=int(nKernelParams//nX)+1, ncols=nX)
535  arrAx = arrAx[::-1, :]
536  allParams = np.array(allParams)
537  for k in range(nKernelParams):
538  ax = arrAx.ravel()[k]
539  ax.hist(allParams[:, k], bins=20, edgecolor='black')
540  ax.set_xlabel('P{}'.format(k))
541  valueParam = spatialFuncs[k](imgCenterX, imgCenterY)
542  ax.axvline(x=valueParam, color='red')
543  ax.text(0.1, 0.9, '{:.1f}'.format(valueParam),
544  transform=ax.transAxes, backgroundcolor='lightsteelblue')
545 
546  # Plot grid of the spatial solution
547  # ----
548 
549  nX = 8
550  nY = 8
551  wCell = wImage / nX
552  hCell = hImage / nY
553  x0 += wCell / 2
554  y0 += hCell / 2
555 
556  fig = plt.figure()
557  fig.suptitle("Spatial solution of kernel parameters on an image grid")
558  arrAx = fig.subplots(nrows=nY, ncols=nX, sharex=True, sharey=True, gridspec_kw=dict(
559  wspace=0, hspace=0))
560  arrAx = arrAx[::-1, :]
561  kernelParams = np.zeros(nKernelParams, dtype=float)
562 
563  for iX in range(nX):
564  for iY in range(nY):
565  x = x0 + iX * wCell
566  y = y0 + iY * hCell
567  # Evaluate the spatial solution functions for this x,y location
568  kernelParams = [f(x, y) for f in spatialFuncs]
569  arrAx[iY, iX].plot(np.arange(nKernelParams), kernelParams, '.-', drawstyle='steps-mid')
570  arrAx[iY, iX].grid(True, axis='y')
571 
572  global keptPlots
573  if keepPlots and not keptPlots:
574  # Keep plots open when done
575  def show():
576  print("%s: Please close plots when done." % __name__)
577  try:
578  plt.show()
579  except Exception:
580  pass
581  print("Plots closed, exiting...")
582  import atexit
583  atexit.register(show)
584  keptPlots = True
585 
586 
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
def show(frame=None)
Definition: ds9.py:88
def plotKernelCoefficients(spatialKernel, kernelCellSet, showBadCandidates=False, keepPlots=True)
Definition: utils.py:428
def plot(mag, width, centers, clusterId, marker="o", markersize=2, markeredgewidth=0, ltype='-', magType="model", clear=True)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174

◆ plotKernelSpatialModel()

def lsst.ip.diffim.utils.plotKernelSpatialModel (   kernel,
  kernelCellSet,
  showBadCandidates = True,
  numSample = 128,
  keepPlots = True,
  maxCoeff = 10 
)
Plot the Kernel spatial model.

Definition at line 278 of file utils.py.

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

◆ plotPixelResiduals()

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 670 of file utils.py.

672  origVariance=False, nptsFull=1e6, keepPlots=True, titleFs=14):
673  """Plot diffim residuals for LOCAL and SPATIAL models.
674  """
675  candidateResids = []
676  spatialResids = []
677  nonfitResids = []
678 
679  for cell in kernelCellSet.getCellList():
680  for cand in cell.begin(True): # only look at good ones
681  # Be sure
682  if not (cand.getStatus() == afwMath.SpatialCellCandidate.GOOD):
683  continue
684 
685  diffim = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
686  orig = cand.getScienceMaskedImage()
687 
688  ski = afwImage.ImageD(kernel.getDimensions())
689  kernel.computeImage(ski, False, int(cand.getXCenter()), int(cand.getYCenter()))
690  sk = afwMath.FixedKernel(ski)
691  sbg = background(int(cand.getXCenter()), int(cand.getYCenter()))
692  sdiffim = cand.getDifferenceImage(sk, sbg)
693 
694  # trim edgs due to convolution
695  bbox = kernel.shrinkBBox(diffim.getBBox())
696  tdiffim = diffim.Factory(diffim, bbox)
697  torig = orig.Factory(orig, bbox)
698  tsdiffim = sdiffim.Factory(sdiffim, bbox)
699 
700  if origVariance:
701  candidateResids.append(np.ravel(tdiffim.getImage().getArray()
702  / np.sqrt(torig.getVariance().getArray())))
703  spatialResids.append(np.ravel(tsdiffim.getImage().getArray()
704  / np.sqrt(torig.getVariance().getArray())))
705  else:
706  candidateResids.append(np.ravel(tdiffim.getImage().getArray()
707  / np.sqrt(tdiffim.getVariance().getArray())))
708  spatialResids.append(np.ravel(tsdiffim.getImage().getArray()
709  / np.sqrt(tsdiffim.getVariance().getArray())))
710 
711  fullIm = diffExposure.getMaskedImage().getImage().getArray()
712  fullMask = diffExposure.getMaskedImage().getMask().getArray()
713  if origVariance:
714  fullVar = exposure.getMaskedImage().getVariance().getArray()
715  else:
716  fullVar = diffExposure.getMaskedImage().getVariance().getArray()
717 
718  bitmaskBad = 0
719  bitmaskBad |= afwImage.Mask.getPlaneBitMask('NO_DATA')
720  bitmaskBad |= afwImage.Mask.getPlaneBitMask('SAT')
721  idx = np.where((fullMask & bitmaskBad) == 0)
722  stride = int(len(idx[0])//nptsFull)
723  sidx = idx[0][::stride], idx[1][::stride]
724  allResids = fullIm[sidx]/np.sqrt(fullVar[sidx])
725 
726  testFootprints = diffimTools.sourceToFootprintList(testSources, warpedTemplateExposure,
727  exposure, config, Log.getDefaultLogger())
728  for fp in testFootprints:
729  subexp = diffExposure.Factory(diffExposure, fp["footprint"].getBBox())
730  subim = subexp.getMaskedImage().getImage()
731  if origVariance:
732  subvar = afwImage.ExposureF(exposure, fp["footprint"].getBBox()).getMaskedImage().getVariance()
733  else:
734  subvar = subexp.getMaskedImage().getVariance()
735  nonfitResids.append(np.ravel(subim.getArray()/np.sqrt(subvar.getArray())))
736 
737  candidateResids = np.ravel(np.array(candidateResids))
738  spatialResids = np.ravel(np.array(spatialResids))
739  nonfitResids = np.ravel(np.array(nonfitResids))
740 
741  try:
742  import pylab
743  from matplotlib.font_manager import FontProperties
744  except ImportError as e:
745  print("Unable to import pylab: %s" % e)
746  return
747 
748  fig = pylab.figure()
749  fig.clf()
750  try:
751  fig.canvas._tkcanvas._root().lift() # == Tk's raise, but raise is a python reserved word
752  except Exception: # protect against API changes
753  pass
754  if origVariance:
755  fig.suptitle("Diffim residuals: Normalized by sqrt(input variance)", fontsize=titleFs)
756  else:
757  fig.suptitle("Diffim residuals: Normalized by sqrt(diffim variance)", fontsize=titleFs)
758 
759  sp1 = pylab.subplot(221)
760  sp2 = pylab.subplot(222, sharex=sp1, sharey=sp1)
761  sp3 = pylab.subplot(223, sharex=sp1, sharey=sp1)
762  sp4 = pylab.subplot(224, sharex=sp1, sharey=sp1)
763  xs = np.arange(-5, 5.05, 0.1)
764  ys = 1./np.sqrt(2*np.pi)*np.exp(-0.5*xs**2)
765 
766  sp1.hist(candidateResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
767  % (np.mean(candidateResids), np.var(candidateResids)))
768  sp1.plot(xs, ys, "r-", lw=2, label="N(0,1)")
769  sp1.set_title("Candidates: basis fit", fontsize=titleFs - 2)
770  sp1.legend(loc=1, fancybox=True, shadow=True, prop=FontProperties(size=titleFs - 6))
771 
772  sp2.hist(spatialResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
773  % (np.mean(spatialResids), np.var(spatialResids)))
774  sp2.plot(xs, ys, "r-", lw=2, label="N(0,1)")
775  sp2.set_title("Candidates: spatial fit", fontsize=titleFs - 2)
776  sp2.legend(loc=1, fancybox=True, shadow=True, prop=FontProperties(size=titleFs - 6))
777 
778  sp3.hist(nonfitResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
779  % (np.mean(nonfitResids), np.var(nonfitResids)))
780  sp3.plot(xs, ys, "r-", lw=2, label="N(0,1)")
781  sp3.set_title("Control sample: spatial fit", fontsize=titleFs - 2)
782  sp3.legend(loc=1, fancybox=True, shadow=True, prop=FontProperties(size=titleFs - 6))
783 
784  sp4.hist(allResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
785  % (np.mean(allResids), np.var(allResids)))
786  sp4.plot(xs, ys, "r-", lw=2, label="N(0,1)")
787  sp4.set_title("Full image (subsampled)", fontsize=titleFs - 2)
788  sp4.legend(loc=1, fancybox=True, shadow=True, prop=FontProperties(size=titleFs - 6))
789 
790  pylab.setp(sp1.get_xticklabels() + sp1.get_yticklabels(), fontsize=titleFs - 4)
791  pylab.setp(sp2.get_xticklabels() + sp2.get_yticklabels(), fontsize=titleFs - 4)
792  pylab.setp(sp3.get_xticklabels() + sp3.get_yticklabels(), fontsize=titleFs - 4)
793  pylab.setp(sp4.get_xticklabels() + sp4.get_yticklabels(), fontsize=titleFs - 4)
794 
795  sp1.set_xlim(-5, 5)
796  sp1.set_ylim(0, 0.5)
797  fig.show()
798 
799  global keptPlots
800  if keepPlots and not keptPlots:
801  # Keep plots open when done
802  def show():
803  print("%s: Please close plots when done." % __name__)
804  try:
805  pylab.show()
806  except Exception:
807  pass
808  print("Plots closed, exiting...")
809  import atexit
810  atexit.register(show)
811  keptPlots = True
812 
813 
A kernel created from an Image.
Definition: Kernel.h:471

◆ plotWhisker()

def lsst.ip.diffim.utils.plotWhisker (   results,
  newWcs 
)
Plot whisker diagram of astromeric offsets between results.matches.

Definition at line 889 of file utils.py.

889 def plotWhisker(results, newWcs):
890  """Plot whisker diagram of astromeric offsets between results.matches.
891  """
892  refCoordKey = results.matches[0].first.getTable().getCoordKey()
893  inCentroidKey = results.matches[0].second.getTable().getCentroidSlot().getMeasKey()
894  positions = [m.first.get(refCoordKey) for m in results.matches]
895  residuals = [m.first.get(refCoordKey).getOffsetFrom(
896  newWcs.pixelToSky(m.second.get(inCentroidKey))) for
897  m in results.matches]
898  import matplotlib.pyplot as plt
899  fig = plt.figure()
900  sp = fig.add_subplot(1, 1, 0)
901  xpos = [x[0].asDegrees() for x in positions]
902  ypos = [x[1].asDegrees() for x in positions]
903  xpos.append(0.02*(max(xpos) - min(xpos)) + min(xpos))
904  ypos.append(0.98*(max(ypos) - min(ypos)) + min(ypos))
905  xidxs = np.isfinite(xpos)
906  yidxs = np.isfinite(ypos)
907  X = np.asarray(xpos)[xidxs]
908  Y = np.asarray(ypos)[yidxs]
909  distance = [x[1].asArcseconds() for x in residuals]
910  distance.append(0.2)
911  distance = np.asarray(distance)[xidxs]
912  # NOTE: This assumes that the bearing is measured positive from +RA through North.
913  # From the documentation this is not clear.
914  bearing = [x[0].asRadians() for x in residuals]
915  bearing.append(0)
916  bearing = np.asarray(bearing)[xidxs]
917  U = (distance*np.cos(bearing))
918  V = (distance*np.sin(bearing))
919  sp.quiver(X, Y, U, V)
920  sp.set_title("WCS Residual")
921  plt.show()
922 
923 
def plotWhisker(results, newWcs)
Definition: utils.py:889

◆ printSkyDiffs()

def lsst.ip.diffim.utils.printSkyDiffs (   sources,
  wcs 
)
Print differences in sky coordinates.

The difference is that between the source Position and its Centroid mapped
through Wcs.

Definition at line 845 of file utils.py.

845 def printSkyDiffs(sources, wcs):
846  """Print differences in sky coordinates.
847 
848  The difference is that between the source Position and its Centroid mapped
849  through Wcs.
850  """
851  for s in sources:
852  sCentroid = s.getCentroid()
853  sPosition = s.getCoord().getPosition(geom.degrees)
854  dra = 3600*(sPosition.getX() - wcs.pixelToSky(sCentroid).getPosition(geom.degrees).getX())/0.2
855  ddec = 3600*(sPosition.getY() - wcs.pixelToSky(sCentroid).getPosition(geom.degrees).getY())/0.2
856  if np.isfinite(dra) and np.isfinite(ddec):
857  print(dra, ddec)
858 
859 
def printSkyDiffs(sources, wcs)
Definition: utils.py:845

◆ showDiaSources()

def lsst.ip.diffim.utils.showDiaSources (   sources,
  exposure,
  isFlagged,
  isDipole,
  frame = None 
)
Display Dia Sources.

Definition at line 106 of file utils.py.

106 def showDiaSources(sources, exposure, isFlagged, isDipole, frame=None):
107  """Display Dia Sources.
108  """
109  #
110  # Show us the ccandidates
111  #
112  # Too many mask planes in diffims
113  disp = afwDisplay.Display(frame=frame)
114  for plane in ("BAD", "CR", "EDGE", "INTERPOlATED", "INTRP", "SAT", "SATURATED"):
115  disp.setMaskPlaneColor(plane, color="ignore")
116 
117  mos = afwDisplay.utils.Mosaic()
118  for i in range(len(sources)):
119  source = sources[i]
120  badFlag = isFlagged[i]
121  dipoleFlag = isDipole[i]
122  bbox = source.getFootprint().getBBox()
123  stamp = exposure.Factory(exposure, bbox, True)
124  im = afwDisplay.utils.Mosaic(gutter=1, background=0, mode="x")
125  im.append(stamp.getMaskedImage())
126  lab = "%.1f,%.1f:" % (source.getX(), source.getY())
127  if badFlag:
128  ctype = afwDisplay.RED
129  lab += "BAD"
130  if dipoleFlag:
131  ctype = afwDisplay.YELLOW
132  lab += "DIPOLE"
133  if not badFlag and not dipoleFlag:
134  ctype = afwDisplay.GREEN
135  lab += "OK"
136  mos.append(im.makeMosaic(), lab, ctype)
137  title = "Dia Sources"
138  mosaicImage = mos.makeMosaic(display=disp, title=title)
139  return mosaicImage
140 
141 
def showDiaSources(sources, exposure, isFlagged, isDipole, frame=None)
Definition: utils.py:106

◆ showKernelBasis()

def lsst.ip.diffim.utils.showKernelBasis (   kernel,
  frame = None 
)
Display a Kernel's basis images.

Definition at line 260 of file utils.py.

260 def showKernelBasis(kernel, frame=None):
261  """Display a Kernel's basis images.
262  """
263  mos = afwDisplay.utils.Mosaic()
264 
265  for k in kernel.getKernelList():
266  im = afwImage.ImageD(k.getDimensions())
267  k.computeImage(im, False)
268  mos.append(im)
269 
270  disp = afwDisplay.Display(frame=frame)
271  mos.makeMosaic(display=disp, title="Kernel Basis Images")
272 
273  return mos
274 
def showKernelBasis(kernel, frame=None)
Definition: utils.py:260

◆ showKernelCandidates()

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 142 of file utils.py.

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

◆ showKernelMosaic()

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 587 of file utils.py.

588  showCenter=True, showEllipticity=True):
589  """Show a mosaic of Kernel images.
590  """
591  mos = afwDisplay.utils.Mosaic()
592 
593  x0 = bbox.getBeginX()
594  y0 = bbox.getBeginY()
595  width = bbox.getWidth()
596  height = bbox.getHeight()
597 
598  if not ny:
599  ny = int(nx*float(height)/width + 0.5)
600  if not ny:
601  ny = 1
602 
603  schema = afwTable.SourceTable.makeMinimalSchema()
604  centroidName = "base_SdssCentroid"
605  shapeName = "base_SdssShape"
606  control = measBase.SdssCentroidControl()
607  schema.getAliasMap().set("slot_Centroid", centroidName)
608  schema.getAliasMap().set("slot_Centroid_flag", centroidName + "_flag")
609  centroider = measBase.SdssCentroidAlgorithm(control, centroidName, schema)
610  sdssShape = measBase.SdssShapeControl()
611  shaper = measBase.SdssShapeAlgorithm(sdssShape, shapeName, schema)
612  table = afwTable.SourceTable.make(schema)
613  table.defineCentroid(centroidName)
614  table.defineShape(shapeName)
615 
616  centers = []
617  shapes = []
618  for iy in range(ny):
619  for ix in range(nx):
620  x = int(ix*(width - 1)/(nx - 1)) + x0
621  y = int(iy*(height - 1)/(ny - 1)) + y0
622 
623  im = afwImage.ImageD(kernel.getDimensions())
624  ksum = kernel.computeImage(im, False, x, y)
625  lab = "Kernel(%d,%d)=%.2f" % (x, y, ksum) if False else ""
626  mos.append(im, lab)
627 
628  # SdssCentroidAlgorithm.measure requires an exposure of floats
629  exp = afwImage.makeExposure(afwImage.makeMaskedImage(im.convertF()))
630 
631  w, h = im.getWidth(), im.getHeight()
632  centerX = im.getX0() + w//2
633  centerY = im.getY0() + h//2
634  src = table.makeRecord()
635  spans = afwGeom.SpanSet(exp.getBBox())
636  foot = afwDet.Footprint(spans)
637  foot.addPeak(centerX, centerY, 1)
638  src.setFootprint(foot)
639 
640  try: # The centroider requires a psf, so this will fail if none is attached to exp
641  centroider.measure(src, exp)
642  centers.append((src.getX(), src.getY()))
643 
644  shaper.measure(src, exp)
645  shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
646  except Exception:
647  pass
648 
649  disp = afwDisplay.Display(frame=frame)
650  mos.makeMosaic(display=disp, title=title if title else "Model Kernel", mode=nx)
651 
652  if centers and frame is not None:
653  disp = afwDisplay.Display(frame=frame)
654  i = 0
655  with disp.Buffering():
656  for cen, shape in zip(centers, shapes):
657  bbox = mos.getBBox(i)
658  i += 1
659  xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
660  if showCenter:
661  disp.dot("+", xc, yc, ctype=afwDisplay.BLUE)
662 
663  if showEllipticity:
664  ixx, ixy, iyy = shape
665  disp.dot("@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, ctype=afwDisplay.RED)
666 
667  return mos
668 
669 
Class to describe the properties of a detected object from an image.
Definition: Footprint.h:63
A compact representation of a collection of pixels.
Definition: SpanSet.h:78
daf::base::PropertySet * set
Definition: fits.cc:912
std::shared_ptr< Exposure< ImagePixelT, MaskPixelT, VariancePixelT > > makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, std::shared_ptr< geom::SkyWcs const > wcs=std::shared_ptr< geom::SkyWcs const >())
A function to return an Exposure of the correct type (cf.
Definition: Exposure.h:462
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename std::shared_ptr< Image< ImagePixelT >> image, typename std::shared_ptr< Mask< MaskPixelT >> mask=Mask< MaskPixelT >(), typename std::shared_ptr< Image< VariancePixelT >> variance=Image< VariancePixelT >())
A function to return a MaskedImage of the correct type (cf.
Definition: MaskedImage.h:1240

◆ showKernelSpatialCells()

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 display.dot understands (e.g. "o"), the top
nMaxPerCell candidates will be indicated with that symbol, using ctype
and size.

Definition at line 67 of file utils.py.

69  frame=None, title="Spatial Cells"):
70  """Show the SpatialCells.
71 
72  If symb is something that display.dot understands (e.g. "o"), the top
73  nMaxPerCell candidates will be indicated with that symbol, using ctype
74  and size.
75  """
76  disp = afwDisplay.Display(frame=frame)
77  disp.mtv(maskedIm, title=title)
78  with disp.Buffering():
79  origin = [-maskedIm.getX0(), -maskedIm.getY0()]
80  for cell in kernelCellSet.getCellList():
81  afwDisplay.utils.drawBBox(cell.getBBox(), origin=origin, display=disp)
82 
83  goodies = ctypeBad is None
84  for cand in cell.begin(goodies):
85  xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
86  if cand.getStatus() == afwMath.SpatialCellCandidate.BAD:
87  color = ctypeBad
88  elif cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
89  color = ctype
90  elif cand.getStatus() == afwMath.SpatialCellCandidate.UNKNOWN:
91  color = ctypeUnused
92  else:
93  continue
94 
95  if color:
96  disp.dot(symb, xc, yc, ctype=color, size=size)
97 
98  if showChi2:
99  rchi2 = cand.getChi2()
100  if rchi2 > 1e100:
101  rchi2 = np.nan
102  disp.dot("%d %.1f" % (cand.getId(), rchi2),
103  xc - size, yc - size - 4, ctype=color, size=size)
104 
105 

◆ showSourceSet()

def lsst.ip.diffim.utils.showSourceSet (   sSet,
  xy0 = (0, 0),
  frame = 0,
  ctype = afwDisplay.GREEN,
  symb = "+",
  size = 2 
)
Draw the (XAstrom, YAstrom) positions of a set of Sources.

Image has the given XY0.

Definition at line 47 of file utils.py.

47 def showSourceSet(sSet, xy0=(0, 0), frame=0, ctype=afwDisplay.GREEN, symb="+", size=2):
48  """Draw the (XAstrom, YAstrom) positions of a set of Sources.
49 
50  Image has the given XY0.
51  """
52  disp = afwDisplay.afwDisplay(frame=frame)
53  with disp.Buffering():
54  for s in sSet:
55  xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
56 
57  if symb == "id":
58  disp.dot(str(s.getId()), xc, yc, ctype=ctype, size=size)
59  else:
60  disp.dot(symb, xc, yc, ctype=ctype, size=size)
61 
62 
63 # Kernel display utilities
64 #
65 
66 
def showSourceSet(sSet, xy0=(0, 0), display=None, ctype=afwDisplay.GREEN, symb="+", size=2)
Definition: utils.py:62

◆ showSourceSetSky()

def lsst.ip.diffim.utils.showSourceSetSky (   sSet,
  wcs,
  xy0,
  frame = 0,
  ctype = afwDisplay.GREEN,
  symb = "+",
  size = 2 
)
Draw the (RA, Dec) positions of a set of Sources. Image has the XY0.

Definition at line 877 of file utils.py.

877 def showSourceSetSky(sSet, wcs, xy0, frame=0, ctype=afwDisplay.GREEN, symb="+", size=2):
878  """Draw the (RA, Dec) positions of a set of Sources. Image has the XY0.
879  """
880  disp = afwDisplay.Display(frame=frame)
881  with disp.Buffering():
882  for s in sSet:
883  (xc, yc) = wcs.skyToPixel(s.getCoord().getRa(), s.getCoord().getDec())
884  xc -= xy0[0]
885  yc -= xy0[1]
886  disp.dot(symb, xc, yc, ctype=ctype, size=size)
887 
888 
def showSourceSetSky(sSet, wcs, xy0, frame=0, ctype=afwDisplay.GREEN, symb="+", size=2)
Definition: utils.py:877

Variable Documentation

◆ keptPlots

bool lsst.ip.diffim.utils.keptPlots = False

Definition at line 44 of file utils.py.