Coverage for tests/utils.py: 34%
690 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-14 19:10 +0000
1#
2# Developed for the LSST Data Management System.
3# This product includes software developed by the LSST Project
4# (https://www.lsst.org).
5# See the COPYRIGHT file at the top-level directory of this distribution
6# for details of code ownership.
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the GNU General Public License
19# along with this program. If not, see <https://www.gnu.org/licenses/>.
21"""Support utilities for Measuring sources"""
23# Export DipoleTestImage to expose fake image generating funcs
24__all__ = ["DipoleTestImage"]
27import astropy.time
28import numpy as np
29import lsst.geom as geom
30import lsst.afw.detection as afwDet
31import lsst.afw.display as afwDisplay
32import lsst.afw.detection as afwDetection
33import lsst.afw.geom as afwGeom
34import lsst.afw.image as afwImage
35from lsst.afw.coord import Observatory
36import lsst.daf.base as dafBase
37import lsst.afw.math as afwMath
38import lsst.afw.table as afwTable
39from lsst.daf.butler import DataCoordinate, DimensionUniverse, Timespan
40import lsst.meas.algorithms as measAlg
41import lsst.meas.base as measBase
42from lsst.meas.algorithms.testUtils import plantSources
43from lsst.ip.diffim.dipoleFitTask import DipoleFitAlgorithm
44from lsst.ip.diffim import diffimLib
46afwDisplay.setDefaultMaskTransparency(75)
47keptPlots = False # Have we arranged to keep spatial plots open?
50def showSourceSet(sSet, xy0=(0, 0), frame=0, ctype=afwDisplay.GREEN, symb="+", size=2):
51 """Draw the (XAstrom, YAstrom) positions of a set of Sources.
53 Image has the given XY0.
54 """
55 disp = afwDisplay.afwDisplay(frame=frame)
56 with disp.Buffering():
57 for s in sSet:
58 xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
60 if symb == "id":
61 disp.dot(str(s.getId()), xc, yc, ctype=ctype, size=size)
62 else:
63 disp.dot(symb, xc, yc, ctype=ctype, size=size)
66# Kernel display utilities
67#
70def showKernelSpatialCells(maskedIm, kernelCellSet, showChi2=False, symb="o",
71 ctype=None, ctypeUnused=None, ctypeBad=None, size=3,
72 frame=None, title="Spatial Cells"):
73 """Show the SpatialCells.
75 If symb is something that display.dot understands (e.g. "o"), the top
76 nMaxPerCell candidates will be indicated with that symbol, using ctype
77 and size.
78 """
79 disp = afwDisplay.Display(frame=frame)
80 disp.mtv(maskedIm, title=title)
81 with disp.Buffering():
82 origin = [-maskedIm.getX0(), -maskedIm.getY0()]
83 for cell in kernelCellSet.getCellList():
84 afwDisplay.utils.drawBBox(cell.getBBox(), origin=origin, display=disp)
86 goodies = ctypeBad is None
87 for cand in cell.begin(goodies):
88 xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
89 if cand.getStatus() == afwMath.SpatialCellCandidate.BAD:
90 color = ctypeBad
91 elif cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
92 color = ctype
93 elif cand.getStatus() == afwMath.SpatialCellCandidate.UNKNOWN:
94 color = ctypeUnused
95 else:
96 continue
98 if color:
99 disp.dot(symb, xc, yc, ctype=color, size=size)
101 if showChi2:
102 rchi2 = cand.getChi2()
103 if rchi2 > 1e100:
104 rchi2 = np.nan
105 disp.dot("%d %.1f" % (cand.getId(), rchi2),
106 xc - size, yc - size - 4, ctype=color, size=size)
109def showDiaSources(sources, exposure, isFlagged, isDipole, frame=None):
110 """Display Dia Sources.
111 """
112 #
113 # Show us the ccandidates
114 #
115 # Too many mask planes in diffims
116 disp = afwDisplay.Display(frame=frame)
117 for plane in ("BAD", "CR", "EDGE", "INTERPOlATED", "INTRP", "SAT", "SATURATED"):
118 disp.setMaskPlaneColor(plane, color="ignore")
120 mos = afwDisplay.utils.Mosaic()
121 for i in range(len(sources)):
122 source = sources[i]
123 badFlag = isFlagged[i]
124 dipoleFlag = isDipole[i]
125 bbox = source.getFootprint().getBBox()
126 stamp = exposure.Factory(exposure, bbox, True)
127 im = afwDisplay.utils.Mosaic(gutter=1, background=0, mode="x")
128 im.append(stamp.getMaskedImage())
129 lab = "%.1f,%.1f:" % (source.getX(), source.getY())
130 if badFlag:
131 ctype = afwDisplay.RED
132 lab += "BAD"
133 if dipoleFlag:
134 ctype = afwDisplay.YELLOW
135 lab += "DIPOLE"
136 if not badFlag and not dipoleFlag:
137 ctype = afwDisplay.GREEN
138 lab += "OK"
139 mos.append(im.makeMosaic(), lab, ctype)
140 title = "Dia Sources"
141 mosaicImage = mos.makeMosaic(display=disp, title=title)
142 return mosaicImage
145def showKernelCandidates(kernelCellSet, kernel, background, frame=None, showBadCandidates=True,
146 resids=False, kernels=False):
147 """Display the Kernel candidates.
149 If kernel is provided include spatial model and residuals;
150 If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi.
151 """
152 #
153 # Show us the ccandidates
154 #
155 if kernels:
156 mos = afwDisplay.utils.Mosaic(gutter=5, background=0)
157 else:
158 mos = afwDisplay.utils.Mosaic(gutter=5, background=-1)
159 #
160 candidateCenters = []
161 candidateCentersBad = []
162 candidateIndex = 0
163 for cell in kernelCellSet.getCellList():
164 for cand in cell.begin(False): # include bad candidates
165 # Original difference image; if does not exist, skip candidate
166 try:
167 resid = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
168 except Exception:
169 continue
171 rchi2 = cand.getChi2()
172 if rchi2 > 1e100:
173 rchi2 = np.nan
175 if not showBadCandidates and cand.isBad():
176 continue
178 im_resid = afwDisplay.utils.Mosaic(gutter=1, background=-0.5, mode="x")
180 try:
181 im = cand.getScienceMaskedImage()
182 im = im.Factory(im, True)
183 im.setXY0(cand.getScienceMaskedImage().getXY0())
184 except Exception:
185 continue
186 if (not resids and not kernels):
187 im_resid.append(im.Factory(im, True))
188 try:
189 im = cand.getTemplateMaskedImage()
190 im = im.Factory(im, True)
191 im.setXY0(cand.getTemplateMaskedImage().getXY0())
192 except Exception:
193 continue
194 if (not resids and not kernels):
195 im_resid.append(im.Factory(im, True))
197 # Difference image with original basis
198 if resids:
199 var = resid.variance
200 var = var.Factory(var, True)
201 np.sqrt(var.array, var.array) # inplace sqrt
202 resid = resid.image
203 resid /= var
204 bbox = kernel.shrinkBBox(resid.getBBox())
205 resid = resid.Factory(resid, bbox, deep=True)
206 elif kernels:
207 kim = cand.getKernelImage(diffimLib.KernelCandidateF.ORIG).convertF()
208 resid = kim.Factory(kim, True)
209 im_resid.append(resid)
211 # residuals using spatial model
212 ski = afwImage.ImageD(kernel.getDimensions())
213 kernel.computeImage(ski, False, int(cand.getXCenter()), int(cand.getYCenter()))
214 sk = afwMath.FixedKernel(ski)
215 sbg = 0.0
216 if background:
217 sbg = background(int(cand.getXCenter()), int(cand.getYCenter()))
218 sresid = cand.getDifferenceImage(sk, sbg)
219 resid = sresid
220 if resids:
221 resid = sresid.image
222 resid /= var
223 bbox = kernel.shrinkBBox(resid.getBBox())
224 resid = resid.Factory(resid, bbox, deep=True)
225 elif kernels:
226 kim = ski.convertF()
227 resid = kim.Factory(kim, True)
228 im_resid.append(resid)
230 im = im_resid.makeMosaic()
232 lab = "%d chi^2 %.1f" % (cand.getId(), rchi2)
233 ctype = afwDisplay.RED if cand.isBad() else afwDisplay.GREEN
235 mos.append(im, lab, ctype)
237 if False and np.isnan(rchi2):
238 disp = afwDisplay.Display(frame=1)
239 disp.mtv(cand.getScienceMaskedImage.image, title="candidate")
240 print("rating", cand.getCandidateRating())
242 im = cand.getScienceMaskedImage()
243 center = (candidateIndex, cand.getXCenter() - im.getX0(), cand.getYCenter() - im.getY0())
244 candidateIndex += 1
245 if cand.isBad():
246 candidateCentersBad.append(center)
247 else:
248 candidateCenters.append(center)
250 if resids:
251 title = "chi Diffim"
252 elif kernels:
253 title = "Kernels"
254 else:
255 title = "Candidates & residuals"
257 disp = afwDisplay.Display(frame=frame)
258 mosaicImage = mos.makeMosaic(display=disp, title=title)
260 return mosaicImage
263def showKernelBasis(kernel, frame=None):
264 """Display a Kernel's basis images.
265 """
266 mos = afwDisplay.utils.Mosaic()
268 for k in kernel.getKernelList():
269 im = afwImage.ImageD(k.getDimensions())
270 k.computeImage(im, False)
271 mos.append(im)
273 disp = afwDisplay.Display(frame=frame)
274 mos.makeMosaic(display=disp, title="Kernel Basis Images")
276 return mos
278###############
281def plotKernelSpatialModel(kernel, kernelCellSet, showBadCandidates=True,
282 numSample=128, keepPlots=True, maxCoeff=10):
283 """Plot the Kernel spatial model.
284 """
285 try:
286 import matplotlib.pyplot as plt
287 import matplotlib.colors
288 except ImportError as e:
289 print("Unable to import numpy and matplotlib: %s" % e)
290 return
292 x0 = kernelCellSet.getBBox().getBeginX()
293 y0 = kernelCellSet.getBBox().getBeginY()
295 candPos = list()
296 candFits = list()
297 badPos = list()
298 badFits = list()
299 candAmps = list()
300 badAmps = list()
301 for cell in kernelCellSet.getCellList():
302 for cand in cell.begin(False):
303 if not showBadCandidates and cand.isBad():
304 continue
305 candCenter = geom.PointD(cand.getXCenter(), cand.getYCenter())
306 try:
307 im = cand.getTemplateMaskedImage()
308 except Exception:
309 continue
311 targetFits = badFits if cand.isBad() else candFits
312 targetPos = badPos if cand.isBad() else candPos
313 targetAmps = badAmps if cand.isBad() else candAmps
315 # compare original and spatial kernel coefficients
316 kp0 = np.array(cand.getKernel(diffimLib.KernelCandidateF.ORIG).getKernelParameters())
317 amp = cand.getCandidateRating()
319 targetFits = badFits if cand.isBad() else candFits
320 targetPos = badPos if cand.isBad() else candPos
321 targetAmps = badAmps if cand.isBad() else candAmps
323 targetFits.append(kp0)
324 targetPos.append(candCenter)
325 targetAmps.append(amp)
327 xGood = np.array([pos.getX() for pos in candPos]) - x0
328 yGood = np.array([pos.getY() for pos in candPos]) - y0
329 zGood = np.array(candFits)
331 xBad = np.array([pos.getX() for pos in badPos]) - x0
332 yBad = np.array([pos.getY() for pos in badPos]) - y0
333 zBad = np.array(badFits)
334 numBad = len(badPos)
336 xRange = np.linspace(0, kernelCellSet.getBBox().getWidth(), num=numSample)
337 yRange = np.linspace(0, kernelCellSet.getBBox().getHeight(), num=numSample)
339 if maxCoeff:
340 maxCoeff = min(maxCoeff, kernel.getNKernelParameters())
341 else:
342 maxCoeff = kernel.getNKernelParameters()
344 for k in range(maxCoeff):
345 func = kernel.getSpatialFunction(k)
346 dfGood = zGood[:, k] - np.array([func(pos.getX(), pos.getY()) for pos in candPos])
347 yMin = dfGood.min()
348 yMax = dfGood.max()
349 if numBad > 0:
350 dfBad = zBad[:, k] - np.array([func(pos.getX(), pos.getY()) for pos in badPos])
351 # Can really screw up the range...
352 yMin = min([yMin, dfBad.min()])
353 yMax = max([yMax, dfBad.max()])
354 yMin -= 0.05*(yMax - yMin)
355 yMax += 0.05*(yMax - yMin)
357 fRange = np.ndarray((len(xRange), len(yRange)))
358 for j, yVal in enumerate(yRange):
359 for i, xVal in enumerate(xRange):
360 fRange[j][i] = func(xVal, yVal)
362 fig = plt.figure(k)
364 fig.clf()
365 try:
366 fig.canvas._tkcanvas._root().lift() # == Tk's raise, but raise is a python reserved word
367 except Exception: # protect against API changes
368 pass
370 fig.suptitle('Kernel component %d' % k)
372 # LL
373 ax = fig.add_axes((0.1, 0.05, 0.35, 0.35))
374 vmin = fRange.min() # - 0.05*np.fabs(fRange.min())
375 vmax = fRange.max() # + 0.05*np.fabs(fRange.max())
376 norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
377 im = ax.imshow(fRange, aspect='auto', norm=norm,
378 extent=[0, kernelCellSet.getBBox().getWidth() - 1,
379 0, kernelCellSet.getBBox().getHeight() - 1])
380 ax.set_title('Spatial polynomial')
381 plt.colorbar(im, orientation='horizontal', ticks=[vmin, vmax])
383 # UL
384 ax = fig.add_axes((0.1, 0.55, 0.35, 0.35))
385 ax.plot(-2.5*np.log10(candAmps), zGood[:, k], 'b+')
386 if numBad > 0:
387 ax.plot(-2.5*np.log10(badAmps), zBad[:, k], 'r+')
388 ax.set_title("Basis Coefficients")
389 ax.set_xlabel("Instr mag")
390 ax.set_ylabel("Coeff")
392 # LR
393 ax = fig.add_axes((0.55, 0.05, 0.35, 0.35))
394 ax.set_autoscale_on(False)
395 ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getHeight())
396 ax.set_ybound(lower=yMin, upper=yMax)
397 ax.plot(yGood, dfGood, 'b+')
398 if numBad > 0:
399 ax.plot(yBad, dfBad, 'r+')
400 ax.axhline(0.0)
401 ax.set_title('dCoeff (indiv-spatial) vs. y')
403 # UR
404 ax = fig.add_axes((0.55, 0.55, 0.35, 0.35))
405 ax.set_autoscale_on(False)
406 ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getWidth())
407 ax.set_ybound(lower=yMin, upper=yMax)
408 ax.plot(xGood, dfGood, 'b+')
409 if numBad > 0:
410 ax.plot(xBad, dfBad, 'r+')
411 ax.axhline(0.0)
412 ax.set_title('dCoeff (indiv-spatial) vs. x')
414 fig.show()
416 global keptPlots
417 if keepPlots and not keptPlots:
418 # Keep plots open when done
419 def show():
420 print("%s: Please close plots when done." % __name__)
421 try:
422 plt.show()
423 except Exception:
424 pass
425 print("Plots closed, exiting...")
426 import atexit
427 atexit.register(show)
428 keptPlots = True
431def plotKernelCoefficients(spatialKernel, kernelCellSet, showBadCandidates=False, keepPlots=True):
432 """Plot the individual kernel candidate and the spatial kernel solution coefficients.
434 Parameters
435 ----------
437 spatialKernel : `lsst.afw.math.LinearCombinationKernel`
438 The spatial spatialKernel solution model which is a spatially varying linear combination
439 of the spatialKernel basis functions.
440 Typically returned by `lsst.ip.diffim.SpatialKernelSolution.getSolutionPair()`.
442 kernelCellSet : `lsst.afw.math.SpatialCellSet`
443 The spatial cells that was used for solution for the spatialKernel. They contain the
444 local solutions of the AL kernel for the selected sources.
446 showBadCandidates : `bool`, optional
447 If True, plot the coefficient values for kernel candidates where the solution was marked
448 bad by the numerical algorithm. Defaults to False.
450 keepPlots: `bool`, optional
451 If True, sets ``plt.show()`` to be called before the task terminates, so that the plots
452 can be explored interactively. Defaults to True.
454 Notes
455 -----
456 This function produces 3 figures per image subtraction operation.
457 * A grid plot of the local solutions. Each grid cell corresponds to a proportional area in
458 the image. In each cell, local kernel solution coefficients are plotted of kernel candidates (color)
459 that fall into this area as a function of the kernel basis function number.
460 * A grid plot of the spatial solution. Each grid cell corresponds to a proportional area in
461 the image. In each cell, the spatial solution coefficients are evaluated for the center of the cell.
462 * Histogram of the local solution coefficients. Red line marks the spatial solution value at
463 center of the image.
465 This function is called if ``lsst.ip.diffim.psfMatch.plotKernelCoefficients==True`` in lsstDebug. This
466 function was implemented as part of DM-17825.
467 """
468 try:
469 import matplotlib.pyplot as plt
470 except ImportError as e:
471 print("Unable to import matplotlib: %s" % e)
472 return
474 # Image dimensions
475 imgBBox = kernelCellSet.getBBox()
476 x0 = imgBBox.getBeginX()
477 y0 = imgBBox.getBeginY()
478 wImage = imgBBox.getWidth()
479 hImage = imgBBox.getHeight()
480 imgCenterX = imgBBox.getCenterX()
481 imgCenterY = imgBBox.getCenterY()
483 # Plot the local solutions
484 # ----
486 # Grid size
487 nX = 8
488 nY = 8
489 wCell = wImage / nX
490 hCell = hImage / nY
492 fig = plt.figure()
493 fig.suptitle("Kernel candidate parameters on an image grid")
494 arrAx = fig.subplots(nrows=nY, ncols=nX, sharex=True, sharey=True, gridspec_kw=dict(
495 wspace=0, hspace=0))
497 # Bottom left panel is for bottom left part of the image
498 arrAx = arrAx[::-1, :]
500 allParams = []
501 for cell in kernelCellSet.getCellList():
502 cellBBox = geom.Box2D(cell.getBBox())
503 # Determine which panel this spatial cell belongs to
504 iX = int((cellBBox.getCenterX() - x0)//wCell)
505 iY = int((cellBBox.getCenterY() - y0)//hCell)
507 for cand in cell.begin(False):
508 try:
509 kernel = cand.getKernel(cand.ORIG)
510 except Exception:
511 continue
513 if not showBadCandidates and cand.isBad():
514 continue
516 nKernelParams = kernel.getNKernelParameters()
517 kernelParams = np.array(kernel.getKernelParameters())
518 allParams.append(kernelParams)
520 if cand.isBad():
521 color = 'red'
522 else:
523 color = None
524 arrAx[iY, iX].plot(np.arange(nKernelParams), kernelParams, '.-',
525 color=color, drawstyle='steps-mid', linewidth=0.1)
526 for ax in arrAx.ravel():
527 ax.grid(True, axis='y')
529 # Plot histogram of the local parameters and the global solution at the image center
530 # ----
532 spatialFuncs = spatialKernel.getSpatialFunctionList()
533 nKernelParams = spatialKernel.getNKernelParameters()
534 nX = 8
535 fig = plt.figure()
536 fig.suptitle("Hist. of parameters marked with spatial solution at img center")
537 arrAx = fig.subplots(nrows=int(nKernelParams//nX)+1, ncols=nX)
538 arrAx = arrAx[::-1, :]
539 allParams = np.array(allParams)
540 for k in range(nKernelParams):
541 ax = arrAx.ravel()[k]
542 ax.hist(allParams[:, k], bins=20, edgecolor='black')
543 ax.set_xlabel('P{}'.format(k))
544 valueParam = spatialFuncs[k](imgCenterX, imgCenterY)
545 ax.axvline(x=valueParam, color='red')
546 ax.text(0.1, 0.9, '{:.1f}'.format(valueParam),
547 transform=ax.transAxes, backgroundcolor='lightsteelblue')
549 # Plot grid of the spatial solution
550 # ----
552 nX = 8
553 nY = 8
554 wCell = wImage / nX
555 hCell = hImage / nY
556 x0 += wCell / 2
557 y0 += hCell / 2
559 fig = plt.figure()
560 fig.suptitle("Spatial solution of kernel parameters on an image grid")
561 arrAx = fig.subplots(nrows=nY, ncols=nX, sharex=True, sharey=True, gridspec_kw=dict(
562 wspace=0, hspace=0))
563 arrAx = arrAx[::-1, :]
564 kernelParams = np.zeros(nKernelParams, dtype=float)
566 for iX in range(nX):
567 for iY in range(nY):
568 x = x0 + iX * wCell
569 y = y0 + iY * hCell
570 # Evaluate the spatial solution functions for this x,y location
571 kernelParams = [f(x, y) for f in spatialFuncs]
572 arrAx[iY, iX].plot(np.arange(nKernelParams), kernelParams, '.-', drawstyle='steps-mid')
573 arrAx[iY, iX].grid(True, axis='y')
575 global keptPlots
576 if keepPlots and not keptPlots:
577 # Keep plots open when done
578 def show():
579 print("%s: Please close plots when done." % __name__)
580 try:
581 plt.show()
582 except Exception:
583 pass
584 print("Plots closed, exiting...")
585 import atexit
586 atexit.register(show)
587 keptPlots = True
590def showKernelMosaic(bbox, kernel, nx=7, ny=None, frame=None, title=None,
591 showCenter=True, showEllipticity=True):
592 """Show a mosaic of Kernel images.
593 """
594 mos = afwDisplay.utils.Mosaic()
596 x0 = bbox.getBeginX()
597 y0 = bbox.getBeginY()
598 width = bbox.getWidth()
599 height = bbox.getHeight()
601 if not ny:
602 ny = int(nx*float(height)/width + 0.5)
603 if not ny:
604 ny = 1
606 schema = afwTable.SourceTable.makeMinimalSchema()
607 centroidName = "base_SdssCentroid"
608 shapeName = "base_SdssShape"
609 control = measBase.SdssCentroidControl()
610 schema.getAliasMap().set("slot_Centroid", centroidName)
611 schema.getAliasMap().set("slot_Centroid_flag", centroidName + "_flag")
612 centroider = measBase.SdssCentroidAlgorithm(control, centroidName, schema)
613 sdssShape = measBase.SdssShapeControl()
614 shaper = measBase.SdssShapeAlgorithm(sdssShape, shapeName, schema)
615 table = afwTable.SourceTable.make(schema)
616 table.defineCentroid(centroidName)
617 table.defineShape(shapeName)
619 centers = []
620 shapes = []
621 for iy in range(ny):
622 for ix in range(nx):
623 x = int(ix*(width - 1)/(nx - 1)) + x0
624 y = int(iy*(height - 1)/(ny - 1)) + y0
626 im = afwImage.ImageD(kernel.getDimensions())
627 ksum = kernel.computeImage(im, False, x, y)
628 lab = "Kernel(%d,%d)=%.2f" % (x, y, ksum) if False else ""
629 mos.append(im, lab)
631 # SdssCentroidAlgorithm.measure requires an exposure of floats
632 exp = afwImage.makeExposure(afwImage.makeMaskedImage(im.convertF()))
634 w, h = im.getWidth(), im.getHeight()
635 centerX = im.getX0() + w//2
636 centerY = im.getY0() + h//2
637 src = table.makeRecord()
638 spans = afwGeom.SpanSet(exp.getBBox())
639 foot = afwDet.Footprint(spans)
640 foot.addPeak(centerX, centerY, 1)
641 src.setFootprint(foot)
643 try: # The centroider requires a psf, so this will fail if none is attached to exp
644 centroider.measure(src, exp)
645 centers.append((src.getX(), src.getY()))
647 shaper.measure(src, exp)
648 shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
649 except Exception:
650 pass
652 disp = afwDisplay.Display(frame=frame)
653 mos.makeMosaic(display=disp, title=title if title else "Model Kernel", mode=nx)
655 if centers and frame is not None:
656 disp = afwDisplay.Display(frame=frame)
657 i = 0
658 with disp.Buffering():
659 for cen, shape in zip(centers, shapes):
660 bbox = mos.getBBox(i)
661 i += 1
662 xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
663 if showCenter:
664 disp.dot("+", xc, yc, ctype=afwDisplay.BLUE)
666 if showEllipticity:
667 ixx, ixy, iyy = shape
668 disp.dot("@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, ctype=afwDisplay.RED)
670 return mos
673def plotWhisker(results, newWcs):
674 """Plot whisker diagram of astromeric offsets between results.matches.
675 """
676 refCoordKey = results.matches[0].first.getTable().getCoordKey()
677 inCentroidKey = results.matches[0].second.getTable().getCentroidSlot().getMeasKey()
678 positions = [m.first.get(refCoordKey) for m in results.matches]
679 residuals = [m.first.get(refCoordKey).getOffsetFrom(
680 newWcs.pixelToSky(m.second.get(inCentroidKey))) for
681 m in results.matches]
682 import matplotlib.pyplot as plt
683 fig = plt.figure()
684 sp = fig.add_subplot(1, 1, 0)
685 xpos = [x[0].asDegrees() for x in positions]
686 ypos = [x[1].asDegrees() for x in positions]
687 xpos.append(0.02*(max(xpos) - min(xpos)) + min(xpos))
688 ypos.append(0.98*(max(ypos) - min(ypos)) + min(ypos))
689 xidxs = np.isfinite(xpos)
690 yidxs = np.isfinite(ypos)
691 X = np.asarray(xpos)[xidxs]
692 Y = np.asarray(ypos)[yidxs]
693 distance = [x[1].asArcseconds() for x in residuals]
694 distance.append(0.2)
695 distance = np.asarray(distance)[xidxs]
696 # NOTE: This assumes that the bearing is measured positive from +RA through North.
697 # From the documentation this is not clear.
698 bearing = [x[0].asRadians() for x in residuals]
699 bearing.append(0)
700 bearing = np.asarray(bearing)[xidxs]
701 U = (distance*np.cos(bearing))
702 V = (distance*np.sin(bearing))
703 sp.quiver(X, Y, U, V)
704 sp.set_title("WCS Residual")
705 plt.show()
708class DipoleTestImage:
709 """Utility class for dipole measurement testing.
711 Generate an image with simulated dipoles and noise; store the original
712 "pre-subtraction" images and catalogs as well.
713 Used to generate test data for DMTN-007 (http://dmtn-007.lsst.io).
714 """
716 def __init__(self, w=101, h=101, xcenPos=[27.], ycenPos=[25.], xcenNeg=[23.], ycenNeg=[25.],
717 psfSigma=2., flux=[30000.], fluxNeg=None, noise=10., gradientParams=None, edgeWidth=8):
718 self.w = w
719 self.h = h
720 self.xcenPos = xcenPos
721 self.ycenPos = ycenPos
722 self.xcenNeg = xcenNeg
723 self.ycenNeg = ycenNeg
724 self.psfSigma = psfSigma
725 self.flux = flux
726 self.fluxNeg = fluxNeg
727 if fluxNeg is None: 727 ↛ 728line 727 didn't jump to line 728 because the condition on line 727 was never true
728 self.fluxNeg = self.flux
729 self.noise = noise
730 self.gradientParams = gradientParams
731 self.edgeWidth = edgeWidth
732 self._makeDipoleImage()
734 def _makeDipoleImage(self):
735 """Generate an exposure and catalog with the given dipole source(s).
736 """
737 # Must seed the pos/neg images with different values to ensure they get different noise realizations
738 posImage, posCatalog = self._makeStarImage(
739 xc=self.xcenPos, yc=self.ycenPos, flux=self.flux, randomSeed=111)
741 negImage, negCatalog = self._makeStarImage(
742 xc=self.xcenNeg, yc=self.ycenNeg, flux=self.fluxNeg, randomSeed=222)
744 dipole = posImage.clone()
745 di = dipole.getMaskedImage()
746 di -= negImage.getMaskedImage()
748 self.diffim, self.posImage, self.posCatalog, self.negImage, self.negCatalog \
749 = dipole, posImage, posCatalog, negImage, negCatalog
751 def _makeStarImage(self, xc=[15.3], yc=[18.6], flux=[2500], schema=None, randomSeed=None):
752 """Generate an exposure and catalog with the given stellar source(s).
753 """
754 from lsst.meas.base.tests import TestDataset
755 bbox = geom.Box2I(geom.Point2I(0, 0), geom.Point2I(self.w - 1, self.h - 1))
756 dataset = TestDataset(bbox, psfSigma=self.psfSigma, threshold=1.)
758 for i in range(len(xc)):
759 dataset.addSource(instFlux=flux[i], centroid=geom.Point2D(xc[i], yc[i]))
761 if schema is None: 761 ↛ 763line 761 didn't jump to line 763 because the condition on line 761 was always true
762 schema = TestDataset.makeMinimalSchema()
763 exposure, catalog = dataset.realize(noise=self.noise, schema=schema, randomSeed=randomSeed)
765 # set EDGE by masking the whole exposure and un-masking an inner bbox
766 edgeMask = exposure.mask.getPlaneBitMask('EDGE')
767 exposure.mask.array |= edgeMask
768 inner_bbox = exposure.getBBox()
769 inner_bbox.grow(-self.edgeWidth)
770 exposure[inner_bbox].mask.array &= ~edgeMask
772 if self.gradientParams is not None: 772 ↛ 781line 772 didn't jump to line 781 because the condition on line 772 was always true
773 y, x = np.mgrid[:self.w, :self.h]
774 gp = self.gradientParams
775 gradient = gp[0] + gp[1]*x + gp[2]*y
776 if len(self.gradientParams) > 3: # it includes a set of 2nd-order polynomial params 776 ↛ 777line 776 didn't jump to line 777 because the condition on line 776 was never true
777 gradient += gp[3]*x*y + gp[4]*x*x + gp[5]*y*y
778 imgArr = exposure.image.array
779 imgArr += gradient
781 return exposure, catalog
783 def fitDipoleSource(self, source, **kwds):
784 alg = DipoleFitAlgorithm(self.diffim, self.posImage, self.negImage)
785 fitResult = alg.fitDipole(source, **kwds)
786 return fitResult
788 def detectDipoleSources(self, doMerge=True, diffim=None, detectSigma=5.5, grow=3, minBinSize=32):
789 """Utility function for detecting dipoles.
791 Detect pos/neg sources in the diffim, then merge them. A
792 bigger "grow" parameter leads to a larger footprint which
793 helps with dipole measurement for faint dipoles.
795 Parameters
796 ----------
797 doMerge : `bool`
798 Whether to merge the positive and negagive detections into a single
799 source table.
800 diffim : `lsst.afw.image.exposure.exposure.ExposureF`
801 Difference image on which to perform detection.
802 detectSigma : `float`
803 Threshold for object detection.
804 grow : `int`
805 Number of pixels to grow the footprints before merging.
806 minBinSize : `int`
807 Minimum bin size for the background (re)estimation (only applies if
808 the default leads to min(nBinX, nBinY) < fit order so the default
809 config parameter needs to be decreased, but not to a value smaller
810 than ``minBinSize``, in which case the fitting algorithm will take
811 over and decrease the fit order appropriately.)
813 Returns
814 -------
815 sources : `lsst.afw.table.SourceCatalog`
816 If doMerge=True, the merged source catalog is returned OR
817 detectTask : `lsst.meas.algorithms.SourceDetectionTask`
818 schema : `lsst.afw.table.Schema`
819 If doMerge=False, the source detection task and its schema are
820 returned.
821 """
822 if diffim is None: 822 ↛ 826line 822 didn't jump to line 826 because the condition on line 822 was always true
823 diffim = self.diffim
825 # Start with a minimal schema - only the fields all SourceCatalogs need
826 schema = afwTable.SourceTable.makeMinimalSchema()
828 # Customize the detection task a bit (optional)
829 detectConfig = measAlg.SourceDetectionConfig()
830 detectConfig.returnOriginalFootprints = False # should be the default
832 # code from imageDifference.py:
833 detectConfig.thresholdPolarity = "both"
834 detectConfig.thresholdValue = detectSigma
835 # detectConfig.nSigmaToGrow = psfSigma
836 detectConfig.reEstimateBackground = True # if False, will fail often for faint sources on gradients?
837 detectConfig.thresholdType = "pixel_stdev"
838 detectConfig.excludeMaskPlanes = ["EDGE"]
839 # Test images are often quite small, so may need to adjust background binSize
840 while ((min(diffim.getWidth(), diffim.getHeight()))//detectConfig.background.binSize
841 < detectConfig.background.approxOrderX and detectConfig.background.binSize > minBinSize):
842 detectConfig.background.binSize = max(minBinSize, detectConfig.background.binSize//2)
844 # Create the detection task. We pass the schema so the task can declare a few flag fields
845 detectTask = measAlg.SourceDetectionTask(schema, config=detectConfig)
847 table = afwTable.SourceTable.make(schema)
848 catalog = detectTask.run(table, diffim)
850 # Now do the merge.
851 if doMerge:
852 fpSet = catalog.positive
853 fpSet.merge(catalog.negative, grow, grow, False)
854 sources = afwTable.SourceCatalog(table)
855 fpSet.makeSources(sources)
857 return sources
859 else:
860 return detectTask, schema
863def detectTestSources(exposure, addMaskPlanes=None):
864 """Minimal source detection wrapper suitable for unit tests.
866 Parameters
867 ----------
868 exposure : `lsst.afw.image.Exposure`
869 Exposure on which to run detection/measurement
870 The exposure is modified in place to set the 'DETECTED' mask plane.
871 addMaskPlanes : `list` of `str`, optional
872 Additional mask planes to add to the maskedImage of the exposure.
874 Returns
875 -------
876 selectSources
877 Source catalog containing candidates
878 """
879 if addMaskPlanes is None:
880 # And add empty INJECTED and INJECTED_TEMPLATE mask planes
881 addMaskPlanes = ["INJECTED", "INJECTED_TEMPLATE"]
883 schema = afwTable.SourceTable.makeMinimalSchema()
885 selectDetection = measAlg.SourceDetectionTask(schema=schema)
886 selectMeasurement = measBase.SingleFrameMeasurementTask(schema=schema)
887 table = afwTable.SourceTable.make(schema)
889 detRet = selectDetection.run(
890 table=table,
891 exposure=exposure,
892 sigma=None, # The appropriate sigma is calculated from the PSF
893 doSmooth=True
894 )
895 for mp in addMaskPlanes:
896 exposure.mask.addMaskPlane(mp)
898 selectSources = detRet.sources
899 selectMeasurement.run(measCat=selectSources, exposure=exposure)
901 return selectSources
904def makeFakeWcs():
905 """Make a fake, affine Wcs.
906 """
907 crpix = geom.Point2D(123.45, 678.9)
908 crval = geom.SpherePoint(0.1, 0.1, geom.degrees)
909 cdMatrix = np.array([[5.19513851e-05, -2.81124812e-07],
910 [-3.25186974e-07, -5.19112119e-05]])
911 return afwGeom.makeSkyWcs(crpix, crval, cdMatrix)
914def makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5.,
915 noiseSeed=6, fluxLevel=500., fluxRange=2.,
916 kernelSize=31, templateBorderSize=0,
917 background=None,
918 xSize=256,
919 ySize=256,
920 x0=12345,
921 y0=67890,
922 calibration=1.,
923 doApplyCalibration=False,
924 xLoc=None,
925 yLoc=None,
926 flux=None,
927 clearEdgeMask=False,
928 addMaskPlanes=None,
929 band="g",
930 physicalFilter="g NotACamera"
931 ):
932 """Make a reproduceable PSF-convolved exposure for testing.
934 Parameters
935 ----------
936 seed : `int`, optional
937 Seed value to initialize the random number generator for sources.
938 nSrc : `int`, optional
939 Number of sources to simulate.
940 psfSize : `float`, optional
941 Width of the PSF of the simulated sources, in pixels.
942 noiseLevel : `float`, optional
943 Standard deviation of the noise to add to each pixel.
944 noiseSeed : `int`, optional
945 Seed value to initialize the random number generator for noise.
946 fluxLevel : `float`, optional
947 Reference flux of the simulated sources.
948 fluxRange : `float`, optional
949 Range in flux amplitude of the simulated sources.
950 kernelSize : `int`, optional
951 Size in pixels of the kernel for simulating sources. Must be odd.
952 templateBorderSize : `int`, optional
953 Size in pixels of the image border used to pad the image.
954 background : `lsst.afw.math.Chebyshev1Function2D`, optional
955 Optional background to add to the output image.
956 xSize, ySize : `int`, optional
957 Size in pixels of the simulated image.
958 x0, y0 : `int`, optional
959 Origin of the image.
960 calibration : `float`, optional
961 Conversion factor between instFlux and nJy.
962 doApplyCalibration : `bool`, optional
963 Apply the photometric calibration and return the image in nJy?
964 xLoc, yLoc : `list` of `float`, optional
965 User-specified coordinates of the simulated sources.
966 If specified, must have length equal to ``nSrc``
967 flux : `list` of `float`, optional
968 User-specified fluxes of the simulated sources.
969 If specified, must have length equal to ``nSrc``
970 clearEdgeMask : `bool`, optional
971 Clear the "EDGE" mask plane after source detection.
972 addMaskPlanes : `list` of `str`, optional
973 Mask plane names to add to the image.
975 Returns
976 -------
977 modelExposure : `lsst.afw.image.Exposure`
978 The model image, with the mask and variance planes. The DETECTED
979 plane is filled in for the injected source footprints.
980 sourceCat : `lsst.afw.table.SourceCatalog`
981 Catalog of sources inserted in the model image.
983 Raises
984 ------
985 ValueError
986 If `xloc`, `yloc`, or `flux` are supplied with inconsistant lengths,
987 or if ``kernelSize`` is even.
988 """
989 if kernelSize % 2 == 0:
990 raise ValueError(f"kernelSize must be odd, got {kernelSize}.")
991 # Distance from the inner edge of the bounding box to avoid placing test
992 # sources in the model images.
993 bufferSize = kernelSize//2 + templateBorderSize + 1
995 bbox = geom.Box2I(geom.Point2I(x0, y0), geom.Extent2I(xSize, ySize))
996 if templateBorderSize > 0:
997 bbox.grow(templateBorderSize)
999 rng = np.random.RandomState(seed)
1000 rngNoise = np.random.RandomState(noiseSeed)
1001 x0, y0 = bbox.getBegin()
1002 xSize, ySize = bbox.getDimensions()
1003 if xLoc is None:
1004 xLoc = rng.rand(nSrc)*(xSize - 2*bufferSize) + bufferSize + x0
1005 else:
1006 if len(xLoc) != nSrc: 1006 ↛ 1007line 1006 didn't jump to line 1007 because the condition on line 1006 was never true
1007 raise ValueError("xLoc must have length equal to nSrc. %f supplied vs %f", len(xLoc), nSrc)
1008 if yLoc is None:
1009 yLoc = rng.rand(nSrc)*(ySize - 2*bufferSize) + bufferSize + y0
1010 else:
1011 if len(yLoc) != nSrc: 1011 ↛ 1012line 1011 didn't jump to line 1012 because the condition on line 1011 was never true
1012 raise ValueError("yLoc must have length equal to nSrc. %f supplied vs %f", len(yLoc), nSrc)
1014 if flux is None: 1014 ↛ 1017line 1014 didn't jump to line 1017 because the condition on line 1014 was always true
1015 flux = (rng.rand(nSrc)*(fluxRange - 1.) + 1.)*fluxLevel
1016 else:
1017 if len(flux) != nSrc:
1018 raise ValueError("flux must have length equal to nSrc. %f supplied vs %f", len(flux), nSrc)
1019 sigmas = [psfSize for src in range(nSrc)]
1020 injectList = list(zip(xLoc, yLoc, flux, sigmas))
1021 skyLevel = 0
1022 # Don't use the built in poisson noise: it modifies the global state of numpy random
1023 modelExposure = plantSources(bbox, kernelSize, skyLevel, injectList, addPoissonNoise=False)
1024 modelExposure.setWcs(makeFakeWcs())
1025 noise = rngNoise.randn(ySize, xSize)*noiseLevel
1026 noise -= np.mean(noise)
1027 modelExposure.variance.array = np.sqrt(np.abs(modelExposure.image.array)) + noiseLevel**2
1028 modelExposure.image.array += noise
1030 # Run source detection to set up the mask plane
1031 detectTestSources(modelExposure, addMaskPlanes=addMaskPlanes)
1032 if clearEdgeMask:
1033 modelExposure.mask &= ~modelExposure.mask.getPlaneBitMask("EDGE")
1034 modelExposure.setPhotoCalib(afwImage.PhotoCalib(calibration, 0., bbox))
1035 if background is not None:
1036 modelExposure.image += background
1037 modelExposure.maskedImage /= calibration
1038 modelExposure.setFilter(afwImage.FilterLabel(band, physicalFilter))
1039 modelExposure.info.setId(seed)
1040 if doApplyCalibration:
1041 modelExposure.maskedImage = modelExposure.photoCalib.calibrateImage(modelExposure.maskedImage)
1043 truth = _fillTruthCatalog(injectList)
1045 return modelExposure, truth
1048def _makeTruthSchema():
1049 """Make a schema for the truth catalog produced by `makeTestImage`.
1051 Returns
1052 -------
1053 keys : `dict` [`str`]
1054 Fields added to the catalog, to make it easier to set them.
1055 schema : `lsst.afw.table.Schema`
1056 Schema to use to make a "truth" SourceCatalog.
1057 Calib, Ap, and Psf flux slots all are set to ``truth_instFlux``.
1058 """
1059 schema = afwTable.SourceTable.makeMinimalSchema()
1060 keys = {}
1061 # Don't use a FluxResultKey so we can manage the flux and err separately.
1062 keys["instFlux"] = schema.addField("truth_instFlux", type=np.float64,
1063 doc="true instFlux", units="count")
1064 keys["instFluxErr"] = schema.addField("truth_instFluxErr", type=np.float64,
1065 doc="true instFluxErr", units="count")
1066 keys["centroid"] = afwTable.Point2DKey.addFields(schema, "truth", "true simulated centroid", "pixel")
1067 schema.addField("truth_flag", "Flag", "truth flux failure flag.")
1068 # Add the flag fields a source selector would need.
1069 schema.addField("sky_source", "Flag", "testing flag.")
1070 schema.addField("base_PixelFlags_flag_interpolated", "Flag", "testing flag.")
1071 schema.addField("base_PixelFlags_flag_saturated", "Flag", "testing flag.")
1072 schema.addField("base_PixelFlags_flag_bad", "Flag", "testing flag.")
1073 schema.addField("base_PixelFlags_flag_edge", "Flag", "testing flag.")
1074 schema.addField("base_PixelFlags_flag_nodata", "Flag", "testing flag.")
1075 schema.addField("base_PsfFlux_flag", "Flag", "testing flag.")
1076 schema.addField("base_ClassificationSizeExtendedness_value", "Flag", "testing flag.")
1077 schema.addField("deblend_nChild", "Flag", "testing flag.")
1078 schema.addField("detect_isPrimary", "Flag", "testing flag.")
1079 schema.addField("calib_psf_used", "Flag", "testing flag.")
1080 schema.getAliasMap().set("slot_Centroid", "truth")
1081 schema.getAliasMap().set("slot_CalibFlux", "truth")
1082 schema.getAliasMap().set("slot_ApFlux", "truth")
1083 schema.getAliasMap().set("slot_PsfFlux", "truth")
1084 schema.addField("glint_trail", "Flag", "testing flag.")
1085 return keys, schema
1088def _fillTruthCatalog(injectList):
1089 """Add injected sources to the truth catalog.
1091 Parameters
1092 ----------
1093 injectList : `list` [`float`]
1094 Sources that were injected; tuples of (x, y, flux, size).
1096 Returns
1097 -------
1098 catalog : `lsst.afw.table.SourceCatalog`
1099 Catalog with centroids and instFlux/instFluxErr values filled in and
1100 appropriate slots set.
1101 """
1102 keys, schema = _makeTruthSchema()
1103 catalog = afwTable.SourceCatalog(schema)
1104 catalog.reserve(len(injectList))
1105 for x, y, flux, size in injectList:
1106 record = catalog.addNew()
1107 keys["centroid"].set(record, geom.PointD(x, y))
1108 keys["instFlux"].set(record, flux)
1109 # Approximate injected errors
1110 keys["instFluxErr"].set(record, 20)
1111 # 5-sigma effective source width
1112 circle = afwGeom.Ellipse(afwGeom.ellipses.Axes(5*size, 5*size, 0), geom.Point2D(x, y))
1113 footprint = afwDetection.Footprint(afwGeom.SpanSet.fromShape(circle))
1114 footprint.addPeak(x, y, flux)
1115 record.setFootprint(footprint)
1117 # Set source records for isolated stars
1118 record["base_ClassificationSizeExtendedness_value"] = 0
1119 record["deblend_nChild"] = 0
1120 record["detect_isPrimary"] = True
1121 record["calib_psf_used"] = True
1123 return catalog
1126def makeStats(badMaskPlanes=None):
1127 """Create a statistics control for configuring calculations on images.
1129 Parameters
1130 ----------
1131 badMaskPlanes : `list` of `str`, optional
1132 List of mask planes to exclude from calculations.
1134 Returns
1135 -------
1136 statsControl : ` lsst.afw.math.StatisticsControl`
1137 Statistics control object for configuring calculations on images.
1138 """
1139 if badMaskPlanes is None:
1140 badMaskPlanes = ("INTRP", "EDGE", "DETECTED", "SAT", "CR",
1141 "BAD", "NO_DATA", "DETECTED_NEGATIVE")
1142 statsControl = afwMath.StatisticsControl()
1143 statsControl.setNumSigmaClip(3.)
1144 statsControl.setNumIter(3)
1145 statsControl.setAndMask(afwImage.Mask.getPlaneBitMask(badMaskPlanes))
1146 return statsControl
1149class CustomCoaddPsf(measAlg.CoaddPsf):
1150 """A custom CoaddPSF that overrides the getAveragePosition method.
1152 It intentionally moves the position off-image to cause a test failure.
1153 """
1154 def getAveragePosition(self):
1155 return geom.Point2D(-10000, -10000)
1158def makeTestVisitInfo(visit=12345):
1159 """Return a visit info for the test images.
1161 Converting an exposure to an `lsst.images` type reads the observation
1162 metadata it still carries from its visit info, so an exposure that will
1163 be converted needs one. The date matches the timespan of
1164 `makeTestExposureRecord`.
1166 Parameters
1167 ----------
1168 visit : `int`, optional
1169 Id of the visit to record.
1171 Returns
1172 -------
1173 visitInfo : `lsst.afw.image.VisitInfo`
1174 The visit info.
1175 """
1176 return afwImage.VisitInfo(
1177 id=visit,
1178 exposureTime=30.0,
1179 date=dafBase.DateTime(60815.000174, dafBase.DateTime.MJD, dafBase.DateTime.TAI),
1180 boresightRaDec=geom.SpherePoint(23.1, 73.2, geom.degrees),
1181 boresightAzAlt=geom.SpherePoint(120.0, 60.0, geom.degrees),
1182 boresightRotAngle=73.2*geom.degrees,
1183 rotType=afwImage.RotType.SKY,
1184 observatory=Observatory(11.1*geom.degrees, 22.2*geom.degrees, 0.333),
1185 )
1188def makeTestExposureRecord(universe, instrument="testCam", visit=12345):
1189 """Return an ``exposure`` dimension record for the test images.
1191 The record holds the observation metadata that converting an image to an
1192 `lsst.images` type needs.
1194 Parameters
1195 ----------
1196 universe : `lsst.daf.butler.DimensionUniverse`
1197 Dimension universe that defines the record's schema.
1198 instrument : `str`, optional
1199 Name of the instrument to record.
1200 visit : `int`, optional
1201 Id of the visit to record. The conversion takes the visit id of an
1202 image from the id of its exposure record.
1204 Returns
1205 -------
1206 record : `lsst.daf.butler.DimensionRecord`
1207 The ``exposure`` record.
1208 """
1209 begin = astropy.time.Time(60815.0, format="mjd", scale="tai")
1210 exposureTime = 30.0
1211 return universe["exposure"].RecordClass(
1212 instrument=instrument,
1213 id=visit,
1214 obs_id=f"test_{visit}",
1215 group=str(visit),
1216 day_obs=20250520,
1217 physical_filter="test filter",
1218 exposure_time=exposureTime,
1219 seq_num=1,
1220 seq_start=1,
1221 seq_end=1,
1222 can_see_sky=True,
1223 azimuth=120.0,
1224 zenith_angle=30.0,
1225 timespan=Timespan(begin=begin, end=begin + astropy.time.TimeDelta(exposureTime, format="sec")),
1226 )
1229def generate_data_id(*,
1230 tract: int = 9813,
1231 patch: int = 42,
1232 cell_x: int = 4,
1233 cell_y: int = 2,
1234 band: str = "notR",
1235 ) -> DataCoordinate:
1236 """Generate a DataCoordinate instance to use as data_id.
1238 Modified from ``generate_data_id`` in ``lsst.cell_coadds.test_utils``
1240 Parameters
1241 ----------
1242 tract : `int`, optional
1243 Tract ID for the data_id
1244 patch : `int`, optional
1245 Patch ID for the data_id
1246 cell_x : `int`, optional
1247 X index of the cell this patch corresponds to
1248 cell_y : `int`, optional
1249 Y index of the cell this patch corresponds to
1250 band : `str`, optional
1251 Band for the data_id
1253 Returns
1254 -------
1255 data_id : `lsst.daf.butler.DataCoordinate`
1256 An expanded data_id instance.
1257 """
1258 universe = DimensionUniverse()
1260 instrument = universe["instrument"]
1261 instrument_record = instrument.RecordClass(
1262 name="DummyCam",
1263 class_name="lsst.obs.base.instrument_tests.DummyCam",
1264 )
1266 skymap = universe["skymap"]
1267 skymap_record = skymap.RecordClass(name="test_skymap")
1269 band_element = universe["band"]
1270 band_record = band_element.RecordClass(name=band)
1272 physical_filter = universe["physical_filter"]
1273 physical_filter_record = physical_filter.RecordClass(name=band, instrument="test", band=band)
1275 patch_element = universe["patch"]
1276 patch_record = patch_element.RecordClass(
1277 skymap="test_skymap", tract=tract, patch=patch, cell_x=cell_x, cell_y=cell_y
1278 )
1280 # A dictionary with all the relevant records.
1281 record = {
1282 "instrument": instrument_record,
1283 "patch": patch_record,
1284 "tract": 9813,
1285 "band": band_record.name,
1286 "skymap": skymap_record.name,
1287 "physical_filter": physical_filter_record,
1288 }
1290 # A dictionary with all the relevant recordIds.
1291 record_id = record.copy()
1292 for key in (
1293 "instrument",
1294 "physical_filter",
1295 ):
1296 record_id[key] = record_id[key].name
1298 data_id = DataCoordinate.standardize(record_id, universe=universe)
1299 return data_id.expanded(record)
1302def checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True):
1303 """Exclude sources that are located on or adjacent to masked pixels.
1305 Parameters
1306 ----------
1307 mask : `lsst.afw.image.Mask`
1308 The image mask plane to use to reject sources
1309 based on the location of their centroid on the ccd.
1310 sources : `lsst.afw.table.SourceCatalog`
1311 The source catalog to evaluate.
1312 excludeMaskPlanes : `list` of `str`
1313 List of the names of the mask planes to exclude.
1315 Returns
1316 -------
1317 flags : `numpy.ndarray` of `bool`
1318 Array indicating whether each source in the catalog should be
1319 kept (True) or rejected (False) based on the value of the
1320 mask plane at its location.
1321 """
1322 setExcludeMaskPlanes = [
1323 maskPlane for maskPlane in excludeMaskPlanes if maskPlane in mask.getMaskPlaneDict()
1324 ]
1326 excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)
1328 xv = (np.rint(sources.getX() - mask.getX0())).astype(int)
1329 yv = (np.rint(sources.getY() - mask.getY0())).astype(int)
1331 flags = np.ones(len(sources), dtype=bool)
1332 if checkAdjacent: 1332 ↛ 1335line 1332 didn't jump to line 1335 because the condition on line 1332 was always true
1333 pixRange = (0, -1, 1)
1334 else:
1335 pixRange = (0,)
1336 for j in pixRange:
1337 for i in pixRange:
1338 mv = mask.array[yv + j, xv + i]
1339 flags *= np.bitwise_and(mv, excludePixelMask) == 0
1340 return flags