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

def calcCentroid(arr):
    """Calculate first moment of a (kernel) image.
    """
    y, x = arr.shape
    sarr = arr*arr
    xarr = np.asarray([[el for el in range(x)] for el2 in range(y)])
    yarr = np.asarray([[el2 for el in range(x)] for el2 in range(y)])
    narr = xarr*sarr
    sarrSum = sarr.sum()
    centx = narr.sum()/sarrSum
    narr = yarr*sarr
    centy = narr.sum()/sarrSum
    return centx, centy

calcWidth()

def lsst.ip.diffim.utils.calcWidth	(		arr,
			centx,
			centy
	)

Calculate second moment of a (kernel) image.

Definition at line 664 of file utils.py.

def calcWidth(arr, centx, centy):
    """Calculate second moment of a (kernel) image.
    """
    y, x = arr.shape
    # Square the flux so we don't have to deal with negatives
    sarr = arr*arr
    xarr = np.asarray([[el for el in range(x)] for el2 in range(y)])
    yarr = np.asarray([[el2 for el in range(x)] for el2 in range(y)])
    narr = sarr*np.power((xarr - centx), 2.)
    sarrSum = sarr.sum()
    xstd = np.sqrt(narr.sum()/sarrSum)
    narr = sarr*np.power((yarr - centy), 2.)
    ystd = np.sqrt(narr.sum()/sarrSum)
    return xstd, ystd

makeRegions()

def lsst.ip.diffim.utils.makeRegions	(		sources,
			outfilename,
			wcs = None
	)

Create regions file for display from input source list.

Definition at line 695 of file utils.py.

def makeRegions(sources, outfilename, wcs=None):
    """Create regions file for display from input source list.
    """
    fh = open(outfilename, "w")
    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")
    for s in sources:
        if wcs:
            (ra, dec) = wcs.pixelToSky(s.getCentroid()).getPosition(afwGeom.degrees)
        else:
            (ra, dec) = s.getCoord().getPosition(afwGeom.degrees)
        if np.isfinite(ra) and np.isfinite(dec):
            fh.write("circle(%f,%f,2\")\n"%(ra, dec))
    fh.flush()
    fh.close()

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

                           numSample=128, keepPlots=True, maxCoeff=10):
    """Plot the Kernel spatial model.
    """
    try:
        import matplotlib.pyplot as plt
        import matplotlib.colors
    except ImportError as e:
        print("Unable to import numpy and matplotlib: %s" % e)
        return

    x0 = kernelCellSet.getBBox().getBeginX()
    y0 = kernelCellSet.getBBox().getBeginY()

    candPos = list()
    candFits = list()
    badPos = list()
    badFits = list()
    candAmps = list()
    badAmps = list()
    for cell in kernelCellSet.getCellList():
        for cand in cell.begin(False):
            if not showBadCandidates and cand.isBad():
                continue
            candCenter = afwGeom.PointD(cand.getXCenter(), cand.getYCenter())
            try:
                im = cand.getTemplateMaskedImage()
            except Exception:
                continue

            targetFits = badFits if cand.isBad() else candFits
            targetPos = badPos if cand.isBad() else candPos
            targetAmps = badAmps if cand.isBad() else candAmps

            # compare original and spatial kernel coefficients
            kp0 = np.array(cand.getKernel(diffimLib.KernelCandidateF.ORIG).getKernelParameters())
            amp = cand.getCandidateRating()

            targetFits = badFits if cand.isBad() else candFits
            targetPos = badPos if cand.isBad() else candPos
            targetAmps = badAmps if cand.isBad() else candAmps

            targetFits.append(kp0)
            targetPos.append(candCenter)
            targetAmps.append(amp)

    xGood = np.array([pos.getX() for pos in candPos]) - x0
    yGood = np.array([pos.getY() for pos in candPos]) - y0
    zGood = np.array(candFits)

    xBad = np.array([pos.getX() for pos in badPos]) - x0
    yBad = np.array([pos.getY() for pos in badPos]) - y0
    zBad = np.array(badFits)
    numBad = len(badPos)

    xRange = np.linspace(0, kernelCellSet.getBBox().getWidth(), num=numSample)
    yRange = np.linspace(0, kernelCellSet.getBBox().getHeight(), num=numSample)

    if maxCoeff:
        maxCoeff = min(maxCoeff, kernel.getNKernelParameters())
    else:
        maxCoeff = kernel.getNKernelParameters()

    for k in range(maxCoeff):
        func = kernel.getSpatialFunction(k)
        dfGood = zGood[:, k] - np.array([func(pos.getX(), pos.getY()) for pos in candPos])
        yMin = dfGood.min()
        yMax = dfGood.max()
        if numBad > 0:
            dfBad = zBad[:, k] - np.array([func(pos.getX(), pos.getY()) for pos in badPos])
            # Can really screw up the range...
            yMin = min([yMin, dfBad.min()])
            yMax = max([yMax, dfBad.max()])
        yMin -= 0.05*(yMax - yMin)
        yMax += 0.05*(yMax - yMin)

        fRange = np.ndarray((len(xRange), len(yRange)))
        for j, yVal in enumerate(yRange):
            for i, xVal in enumerate(xRange):
                fRange[j][i] = func(xVal, yVal)

        fig = plt.figure(k)

        fig.clf()
        try:
            fig.canvas._tkcanvas._root().lift()  # == Tk's raise, but raise is a python reserved word
        except Exception:                                 # protect against API changes
            pass

        fig.suptitle('Kernel component %d' % k)

        # LL
        ax = fig.add_axes((0.1, 0.05, 0.35, 0.35))
        vmin = fRange.min()  # - 0.05*np.fabs(fRange.min())
        vmax = fRange.max()  # + 0.05*np.fabs(fRange.max())
        norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
        im = ax.imshow(fRange, aspect='auto', norm=norm,
                       extent=[0, kernelCellSet.getBBox().getWidth() - 1,
                               0, kernelCellSet.getBBox().getHeight() - 1])
        ax.set_title('Spatial polynomial')
        plt.colorbar(im, orientation='horizontal', ticks=[vmin, vmax])

        # UL
        ax = fig.add_axes((0.1, 0.55, 0.35, 0.35))
        ax.plot(-2.5*np.log10(candAmps), zGood[:, k], 'b+')
        if numBad > 0:
            ax.plot(-2.5*np.log10(badAmps), zBad[:, k], 'r+')
        ax.set_title("Basis Coefficients")
        ax.set_xlabel("Instr mag")
        ax.set_ylabel("Coeff")

        # LR
        ax = fig.add_axes((0.55, 0.05, 0.35, 0.35))
        ax.set_autoscale_on(False)
        ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getHeight())
        ax.set_ybound(lower=yMin, upper=yMax)
        ax.plot(yGood, dfGood, 'b+')
        if numBad > 0:
            ax.plot(yBad, dfBad, 'r+')
        ax.axhline(0.0)
        ax.set_title('dCoeff (indiv-spatial) vs. y')

        # UR
        ax = fig.add_axes((0.55, 0.55, 0.35, 0.35))
        ax.set_autoscale_on(False)
        ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getWidth())
        ax.set_ybound(lower=yMin, upper=yMax)
        ax.plot(xGood, dfGood, 'b+')
        if numBad > 0:
            ax.plot(xBad, dfBad, 'r+')
        ax.axhline(0.0)
        ax.set_title('dCoeff (indiv-spatial) vs. x')

        fig.show()

    global keptPlots
    if keepPlots and not keptPlots:
        # Keep plots open when done
        def show():
            print("%s: Please close plots when done." % __name__)
            try:
                plt.show()
            except Exception:
                pass
            print("Plots closed, exiting...")
        import atexit
        atexit.register(show)
        keptPlots = True

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

                       origVariance=False, nptsFull=1e6, keepPlots=True, titleFs=14):
    """Plot diffim residuals for LOCAL and SPATIAL models.
    """
    candidateResids = []
    spatialResids = []
    nonfitResids = []

    for cell in kernelCellSet.getCellList():
        for cand in cell.begin(True):  # only look at good ones
            # Be sure
            if not (cand.getStatus() == afwMath.SpatialCellCandidate.GOOD):
                continue

            diffim = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
            orig = cand.getScienceMaskedImage()

            ski = afwImage.ImageD(kernel.getDimensions())
            kernel.computeImage(ski, False, int(cand.getXCenter()), int(cand.getYCenter()))
            sk = afwMath.FixedKernel(ski)
            sbg = background(int(cand.getXCenter()), int(cand.getYCenter()))
            sdiffim = cand.getDifferenceImage(sk, sbg)

            # trim edgs due to convolution
            bbox = kernel.shrinkBBox(diffim.getBBox())
            tdiffim = diffim.Factory(diffim, bbox)
            torig = orig.Factory(orig, bbox)
            tsdiffim = sdiffim.Factory(sdiffim, bbox)

            if origVariance:
                candidateResids.append(np.ravel(tdiffim.getImage().getArray() /
                                                np.sqrt(torig.getVariance().getArray())))
                spatialResids.append(np.ravel(tsdiffim.getImage().getArray() /
                                              np.sqrt(torig.getVariance().getArray())))
            else:
                candidateResids.append(np.ravel(tdiffim.getImage().getArray() /
                                                np.sqrt(tdiffim.getVariance().getArray())))
                spatialResids.append(np.ravel(tsdiffim.getImage().getArray() /
                                              np.sqrt(tsdiffim.getVariance().getArray())))

    fullIm = diffExposure.getMaskedImage().getImage().getArray()
    fullMask = diffExposure.getMaskedImage().getMask().getArray()
    if origVariance:
        fullVar = exposure.getMaskedImage().getVariance().getArray()
    else:
        fullVar = diffExposure.getMaskedImage().getVariance().getArray()

    bitmaskBad = 0
    bitmaskBad |= afwImage.Mask.getPlaneBitMask('NO_DATA')
    bitmaskBad |= afwImage.Mask.getPlaneBitMask('SAT')
    idx = np.where((fullMask & bitmaskBad) == 0)
    stride = int(len(idx[0])//nptsFull)
    sidx = idx[0][::stride], idx[1][::stride]
    allResids = fullIm[sidx]/np.sqrt(fullVar[sidx])

    testFootprints = diffimTools.sourceToFootprintList(testSources, warpedTemplateExposure,
                                                       exposure, config, Log.getDefaultLogger())
    for fp in testFootprints:
        subexp = diffExposure.Factory(diffExposure, fp["footprint"].getBBox())
        subim = subexp.getMaskedImage().getImage()
        if origVariance:
            subvar = afwImage.ExposureF(exposure, fp["footprint"].getBBox()).getMaskedImage().getVariance()
        else:
            subvar = subexp.getMaskedImage().getVariance()
        nonfitResids.append(np.ravel(subim.getArray()/np.sqrt(subvar.getArray())))

    candidateResids = np.ravel(np.array(candidateResids))
    spatialResids = np.ravel(np.array(spatialResids))
    nonfitResids = np.ravel(np.array(nonfitResids))

    try:
        import pylab
        from matplotlib.font_manager import FontProperties
    except ImportError as e:
        print("Unable to import pylab: %s" % e)
        return

    fig = pylab.figure()
    fig.clf()
    try:
        fig.canvas._tkcanvas._root().lift()  # == Tk's raise, but raise is a python reserved word
    except Exception:                                 # protect against API changes
        pass
    if origVariance:
        fig.suptitle("Diffim residuals: Normalized by sqrt(input variance)", fontsize=titleFs)
    else:
        fig.suptitle("Diffim residuals: Normalized by sqrt(diffim variance)", fontsize=titleFs)

    sp1 = pylab.subplot(221)
    sp2 = pylab.subplot(222, sharex=sp1, sharey=sp1)
    sp3 = pylab.subplot(223, sharex=sp1, sharey=sp1)
    sp4 = pylab.subplot(224, sharex=sp1, sharey=sp1)
    xs = np.arange(-5, 5.05, 0.1)
    ys = 1./np.sqrt(2*np.pi)*np.exp(-0.5*xs**2)

    sp1.hist(candidateResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
             % (np.mean(candidateResids), np.var(candidateResids)))
    sp1.plot(xs, ys, "r-", lw=2, label="N(0,1)")
    sp1.set_title("Candidates: basis fit", fontsize=titleFs - 2)
    sp1.legend(loc=1, fancybox=True, shadow=True, prop=FontProperties(size=titleFs - 6))

    sp2.hist(spatialResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
             % (np.mean(spatialResids), np.var(spatialResids)))
    sp2.plot(xs, ys, "r-", lw=2, label="N(0,1)")
    sp2.set_title("Candidates: spatial fit", fontsize=titleFs - 2)
    sp2.legend(loc=1, fancybox=True, shadow=True, prop=FontProperties(size=titleFs - 6))

    sp3.hist(nonfitResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
             % (np.mean(nonfitResids), np.var(nonfitResids)))
    sp3.plot(xs, ys, "r-", lw=2, label="N(0,1)")
    sp3.set_title("Control sample: spatial fit", fontsize=titleFs - 2)
    sp3.legend(loc=1, fancybox=True, shadow=True, prop=FontProperties(size=titleFs - 6))

    sp4.hist(allResids, bins=xs, normed=True, alpha=0.5, label="N(%.2f, %.2f)"
             % (np.mean(allResids), np.var(allResids)))
    sp4.plot(xs, ys, "r-", lw=2, label="N(0,1)")
    sp4.set_title("Full image (subsampled)", fontsize=titleFs - 2)
    sp4.legend(loc=1, fancybox=True, shadow=True, prop=FontProperties(size=titleFs - 6))

    pylab.setp(sp1.get_xticklabels() + sp1.get_yticklabels(), fontsize=titleFs - 4)
    pylab.setp(sp2.get_xticklabels() + sp2.get_yticklabels(), fontsize=titleFs - 4)
    pylab.setp(sp3.get_xticklabels() + sp3.get_yticklabels(), fontsize=titleFs - 4)
    pylab.setp(sp4.get_xticklabels() + sp4.get_yticklabels(), fontsize=titleFs - 4)

    sp1.set_xlim(-5, 5)
    sp1.set_ylim(0, 0.5)
    fig.show()

    global keptPlots
    if keepPlots and not keptPlots:
        # Keep plots open when done
        def show():
            print("%s: Please close plots when done." % __name__)
            try:
                pylab.show()
            except Exception:
                pass
            print("Plots closed, exiting...")
        import atexit
        atexit.register(show)
        keptPlots = True

plotWhisker()

def lsst.ip.diffim.utils.plotWhisker	(		results,
			newWcs
	)

Plot whisker diagram of astromeric offsets between results.matches.

Definition at line 724 of file utils.py.

def plotWhisker(results, newWcs):
    """Plot whisker diagram of astromeric offsets between results.matches.
    """
    refCoordKey = results.matches[0].first.getTable().getCoordKey()
    inCentroidKey = results.matches[0].second.getTable().getCentroidKey()
    positions = [m.first.get(refCoordKey) for m in results.matches]
    residuals = [m.first.get(refCoordKey).getOffsetFrom(
        newWcs.pixelToSky(m.second.get(inCentroidKey))) for
        m in results.matches]
    import matplotlib.pyplot as plt
    fig = plt.figure()
    sp = fig.add_subplot(1, 1, 0)
    xpos = [x[0].asDegrees() for x in positions]
    ypos = [x[1].asDegrees() for x in positions]
    xpos.append(0.02*(max(xpos) - min(xpos)) + min(xpos))
    ypos.append(0.98*(max(ypos) - min(ypos)) + min(ypos))
    xidxs = np.isfinite(xpos)
    yidxs = np.isfinite(ypos)
    X = np.asarray(xpos)[xidxs]
    Y = np.asarray(ypos)[yidxs]
    distance = [x[1].asArcseconds() for x in residuals]
    distance.append(0.2)
    distance = np.asarray(distance)[xidxs]
    # NOTE: This assumes that the bearing is measured positive from +RA through North.
    # From the documentation this is not clear.
    bearing = [x[0].asRadians() for x in residuals]
    bearing.append(0)
    bearing = np.asarray(bearing)[xidxs]
    U = (distance*np.cos(bearing))
    V = (distance*np.sin(bearing))
    sp.quiver(X, Y, U, V)
    sp.set_title("WCS Residual")
    plt.show()

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

def printSkyDiffs(sources, wcs):
    """Print differences in sky coordinates.

    The difference is that between the source Position and its Centroid mapped
    through Wcs.
    """
    for s in sources:
        sCentroid = s.getCentroid()
        sPosition = s.getCoord().getPosition(afwGeom.degrees)
        dra = 3600*(sPosition.getX() - wcs.pixelToSky(sCentroid).getPosition(afwGeom.degrees).getX())/0.2
        ddec = 3600*(sPosition.getY() - wcs.pixelToSky(sCentroid).getPosition(afwGeom.degrees).getY())/0.2
        if np.isfinite(dra) and np.isfinite(ddec):
            print(dra, ddec)

showDiaSources()

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

Display Dia Sources.

Definition at line 105 of file utils.py.

def showDiaSources(sources, exposure, isFlagged, isDipole, frame=None):
    """Display Dia Sources.
    """
    #
    # Show us the ccandidates
    #
    # Too many mask planes in diffims
    disp = afwDisplay.Display(frame=frame)
    for plane in ("BAD", "CR", "EDGE", "INTERPOlATED", "INTRP", "SAT", "SATURATED"):
        disp.setMaskPlaneColor(plane, color="ignore")

    mos = afwDisplay.utils.Mosaic()
    for i in range(len(sources)):
        source = sources[i]
        badFlag = isFlagged[i]
        dipoleFlag = isDipole[i]
        bbox = source.getFootprint().getBBox()
        stamp = exposure.Factory(exposure, bbox, True)
        im = afwDisplay.utils.Mosaic(gutter=1, background=0, mode="x")
        im.append(stamp.getMaskedImage())
        lab = "%.1f,%.1f:" % (source.getX(), source.getY())
        if badFlag:
            ctype = afwDisplay.RED
            lab += "BAD"
        if dipoleFlag:
            ctype = afwDisplay.YELLOW
            lab += "DIPOLE"
        if not badFlag and not dipoleFlag:
            ctype = afwDisplay.GREEN
            lab += "OK"
        mos.append(im.makeMosaic(), lab, ctype)
    title = "Dia Sources"
    mosaicImage = mos.makeMosaic(frame=frame, title=title)
    return mosaicImage

showKernelBasis()

def lsst.ip.diffim.utils.showKernelBasis	(		kernel,
			frame = None
	)

Display a Kernel's basis images.

Definition at line 257 of file utils.py.

def showKernelBasis(kernel, frame=None):
    """Display a Kernel's basis images.
    """
    mos = afwDisplay.utils.Mosaic()

    for k in kernel.getKernelList():
        im = afwImage.ImageD(k.getDimensions())
        k.computeImage(im, False)
        mos.append(im)
    mos.makeMosaic(frame=frame, title="Kernel Basis Images")

    return mos

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.

                         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.
    """
    #
    # Show us the ccandidates
    #
    if kernels:
        mos = afwDisplay.utils.Mosaic(gutter=5, background=0)
    else:
        mos = afwDisplay.utils.Mosaic(gutter=5, background=-1)
    #
    candidateCenters = []
    candidateCentersBad = []
    candidateIndex = 0
    for cell in kernelCellSet.getCellList():
        for cand in cell.begin(False):  # include bad candidates
            # Original difference image; if does not exist, skip candidate
            try:
                resid = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
            except Exception:
                continue

            rchi2 = cand.getChi2()
            if rchi2 > 1e100:
                rchi2 = np.nan

            if not showBadCandidates and cand.isBad():
                continue

            im_resid = afwDisplay.utils.Mosaic(gutter=1, background=-0.5, mode="x")

            try:
                im = cand.getScienceMaskedImage()
                im = im.Factory(im, True)
                im.setXY0(cand.getScienceMaskedImage().getXY0())
            except Exception:
                continue
            if (not resids and not kernels):
                im_resid.append(im.Factory(im, True))
            try:
                im = cand.getTemplateMaskedImage()
                im = im.Factory(im, True)
                im.setXY0(cand.getTemplateMaskedImage().getXY0())
            except Exception:
                continue
            if (not resids and not kernels):
                im_resid.append(im.Factory(im, True))

            # Difference image with original basis
            if resids:
                var = resid.getVariance()
                var = var.Factory(var, True)
                np.sqrt(var.getArray(), var.getArray())  # inplace sqrt
                resid = resid.getImage()
                resid /= var
                bbox = kernel.shrinkBBox(resid.getBBox())
                resid = resid.Factory(resid, bbox, deep=True)
            elif kernels:
                kim = cand.getKernelImage(diffimLib.KernelCandidateF.ORIG).convertF()
                resid = kim.Factory(kim, True)
            im_resid.append(resid)

            # residuals using spatial model
            ski = afwImage.ImageD(kernel.getDimensions())
            kernel.computeImage(ski, False, int(cand.getXCenter()), int(cand.getYCenter()))
            sk = afwMath.FixedKernel(ski)
            sbg = 0.0
            if background:
                sbg = background(int(cand.getXCenter()), int(cand.getYCenter()))
            sresid = cand.getDifferenceImage(sk, sbg)
            resid = sresid
            if resids:
                resid = sresid.getImage()
                resid /= var
                bbox = kernel.shrinkBBox(resid.getBBox())
                resid = resid.Factory(resid, bbox, deep=True)
            elif kernels:
                kim = ski.convertF()
                resid = kim.Factory(kim, True)
            im_resid.append(resid)

            im = im_resid.makeMosaic()

            lab = "%d chi^2 %.1f" % (cand.getId(), rchi2)
            ctype = afwDisplay.RED if cand.isBad() else afwDisplay.GREEN

            mos.append(im, lab, ctype)

            if False and np.isnan(rchi2):
                disp = afwDisplay.Display(frame=1)
                disp.mtv(cand.getScienceMaskedImage.getImage(), title="candidate")
                print("rating", cand.getCandidateRating())

            im = cand.getScienceMaskedImage()
            center = (candidateIndex, cand.getXCenter() - im.getX0(), cand.getYCenter() - im.getY0())
            candidateIndex += 1
            if cand.isBad():
                candidateCentersBad.append(center)
            else:
                candidateCenters.append(center)

    if resids:
        title = "chi Diffim"
    elif kernels:
        title = "Kernels"
    else:
        title = "Candidates & residuals"
    mosaicImage = mos.makeMosaic(frame=frame, title=title)

    return mosaicImage

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

                     showCenter=True, showEllipticity=True):
    """Show a mosaic of Kernel images.
    """
    mos = afwDisplay.utils.Mosaic()

    x0 = bbox.getBeginX()
    y0 = bbox.getBeginY()
    width = bbox.getWidth()
    height = bbox.getHeight()

    if not ny:
        ny = int(nx*float(height)/width + 0.5)
        if not ny:
            ny = 1

    schema = afwTable.SourceTable.makeMinimalSchema()
    centroidName = "base_SdssCentroid"
    shapeName = "base_SdssShape"
    control = measBase.SdssCentroidControl()
    schema.getAliasMap().set("slot_Centroid", centroidName)
    schema.getAliasMap().set("slot_Centroid_flag", centroidName + "_flag")
    centroider = measBase.SdssCentroidAlgorithm(control, centroidName, schema)
    sdssShape = measBase.SdssShapeControl()
    shaper = measBase.SdssShapeAlgorithm(sdssShape, shapeName, schema)
    table = afwTable.SourceTable.make(schema)
    table.defineCentroid(centroidName)
    table.defineShape(shapeName)

    centers = []
    shapes = []
    for iy in range(ny):
        for ix in range(nx):
            x = int(ix*(width - 1)/(nx - 1)) + x0
            y = int(iy*(height - 1)/(ny - 1)) + y0

            im = afwImage.ImageD(kernel.getDimensions())
            ksum = kernel.computeImage(im, False, x, y)
            lab = "Kernel(%d,%d)=%.2f" % (x, y, ksum) if False else ""
            mos.append(im, lab)

            # SdssCentroidAlgorithm.measure requires an exposure of floats
            exp = afwImage.makeExposure(afwImage.makeMaskedImage(im.convertF()))

            w, h = im.getWidth(), im.getHeight()
            centerX = im.getX0() + w//2
            centerY = im.getY0() + h//2
            src = table.makeRecord()
            spans = afwGeom.SpanSet(exp.getBBox())
            foot = afwDet.Footprint(spans)
            foot.addPeak(centerX, centerY, 1)
            src.setFootprint(foot)

            try:  # The centroider requires a psf, so this will fail if none is attached to exp
                centroider.measure(src, exp)
                centers.append((src.getX(), src.getY()))

                shaper.measure(src, exp)
                shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
            except Exception:
                pass

    mos.makeMosaic(frame=frame, title=title if title else "Model Kernel", mode=nx)

    if centers and frame is not None:
        disp = afwDisplay.Display(frame=frame)
        i = 0
        with disp.Buffering():
            for cen, shape in zip(centers, shapes):
                bbox = mos.getBBox(i)
                i += 1
                xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
                if showCenter:
                    disp.dot("+", xc, yc, ctype=afwDisplay.BLUE)

                if showEllipticity:
                    ixx, ixy, iyy = shape
                    disp.dot("@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, ctype=afwDisplay.RED)

    return mos

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

                           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.
    """
    disp = afwDisplay.Display(frame=frame)
    disp.mtv(maskedIm, title=title)
    with disp.Buffering():
        origin = [-maskedIm.getX0(), -maskedIm.getY0()]
        for cell in kernelCellSet.getCellList():
            afwDisplay.utils.drawBBox(cell.getBBox(), origin=origin, display=disp)

            goodies = ctypeBad is None
            for cand in cell.begin(goodies):
                xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
                if cand.getStatus() == afwMath.SpatialCellCandidate.BAD:
                    color = ctypeBad
                elif cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
                    color = ctype
                elif cand.getStatus() == afwMath.SpatialCellCandidate.UNKNOWN:
                    color = ctypeUnused
                else:
                    continue

                if color:
                    disp.dot(symb, xc, yc, ctype=color, size=size)

                    if showChi2:
                        rchi2 = cand.getChi2()
                        if rchi2 > 1e100:
                            rchi2 = np.nan
                        disp.dot("%d %.1f" % (cand.getId(), rchi2),
                                 xc - size, yc - size - 4, ctype=color, size=size)

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

def 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.
    """
    disp = afwDisplay.afwDisplay(frame=frame)
    with disp.Buffering():
        for s in sSet:
            xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]

            if symb == "id":
                disp.dot(str(s.getId()), xc, yc, ctype=ctype, size=size)
            else:
                disp.dot(symb, xc, yc, ctype=ctype, size=size)


# Kernel display utilities
#

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

def 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.
    """
    disp = afwDisplay.Display(frame=frame)
    with disp.Buffering():
        for s in sSet:
            (xc, yc) = wcs.skyToPixel(s.getCoord().getRa(), s.getCoord().getDec())
            xc -= xy0[0]
            yc -= xy0[1]
            disp.dot(symb, xc, yc, ctype=ctype, size=size)

Variable Documentation

keptPlots

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

Definition at line 43 of file utils.py.