204 def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None):
205 """Determine a PSFEX PSF model for an exposure given a list of PSF
210 exposure: `lsst.afw.image.Exposure`
211 Exposure containing the PSF candidates.
212 psfCandidateList: iterable of `lsst.meas.algorithms.PsfCandidate`
213 Sequence of PSF candidates typically obtained by detecting sources
214 and then running them through a star selector.
215 metadata: metadata, optional
216 A home for interesting tidbits of information.
217 flagKey: `lsst.afw.table.Key`, optional
218 Schema key used to mark sources actually used in PSF determination.
222 psf: `lsst.meas.extensions.psfex.PsfexPsf`
225 psfCandidateList = self.downsampleCandidates(psfCandidateList)
229 displayExposure = display
and \
231 displayPsfComponents = display
and \
233 showBadCandidates = display
and \
235 displayResiduals = display
and \
237 displayPsfMosaic = display
and \
240 afwDisplay.setDefaultMaskTransparency(75)
243 mi = exposure.getMaskedImage()
245 nCand = len(psfCandidateList)
253 bbox = mi.getBBox(afwImage.PARENT)
258 sizes = np.empty(nCand)
259 for i, psfCandidate
in enumerate(psfCandidateList):
262 psfCellSet.insertCandidate(psfCandidate)
263 except Exception
as e:
264 self.log.error(
"Skipping PSF candidate %d of %d: %s", i, len(psfCandidateList), e)
267 source = psfCandidate.getSource()
268 quad = afwEll.Quadrupole(source.getIxx(), source.getIyy(), source.getIxy())
269 rmsSize = quad.getTraceRadius()
272 pixKernelSize = self.config.stampSize
273 actualKernelSize = int(2*np.floor(0.5*pixKernelSize/self.config.samplingSize) + 1)
276 rms = np.median(sizes)
277 self.log.debug(
"Median PSF RMS size=%.2f pixels (\"FWHM\"=%.2f)",
278 rms, 2*np.sqrt(2*np.log(2))*rms)
280 self.log.trace(
"Psfex Kernel size=%.2f, Image Kernel Size=%.2f", actualKernelSize, pixKernelSize)
286 defaultsFile = os.path.join(os.environ[
"MEAS_EXTENSIONS_PSFEX_DIR"],
"config",
"default-lsst.psfex")
288 args_md.set(
"BASIS_TYPE", str(self.config.psfexBasis))
289 args_md.set(
"PSFVAR_DEGREES", str(self.config.spatialOrder))
290 args_md.set(
"PSF_SIZE", str(actualKernelSize))
291 args_md.set(
"PSF_SAMPLING", str(self.config.samplingSize))
292 args_md.set(
"PHOTFLUX_KEY", str(self.config.photometricFluxField))
293 args_md.set(
"PHOTFLUXERR_KEY", str(self.config.photometricFluxField) +
"Err")
294 prefs = psfex.Prefs(defaultsFile, args_md)
295 prefs.setCommandLine([])
296 prefs.addCatalog(
"psfexPsfDeterminer")
299 principalComponentExclusionFlag = bool(bool(psfex.Context.REMOVEHIDDEN)
300 if False else psfex.Context.KEEPHIDDEN)
301 context = psfex.Context(prefs.getContextName(), prefs.getContextGroup(),
302 prefs.getGroupDeg(), principalComponentExclusionFlag)
303 psfSet = psfex.Set(context)
304 psfSet.setVigSize(pixKernelSize, pixKernelSize)
305 psfSet.setFwhm(2*np.sqrt(2*np.log(2))*np.median(sizes))
306 psfSet.setRecentroid(self.config.recentroid)
310 ccd = exposure.getDetector()
312 gain = np.mean(np.array([a.getGain()
for a
in ccd]))
315 self.log.warning(
"Setting gain to %g", gain)
318 for i, key
in enumerate(context.getName()):
321 contextvalp.append(exposure.getMetadata().getScalar(key[1:]))
322 except KeyError
as e:
323 raise RuntimeError(
"%s parameter not found in the header of %s" %
324 (key[1:], prefs.getContextName()))
from e
327 contextvalp.append(np.array([psfCandidateList[_].getSource().get(key)
328 for _
in range(nCand)]))
329 except KeyError
as e:
330 raise RuntimeError(
"%s parameter not found" % (key,))
from e
331 psfSet.setContextname(i, key)
336 disp = afwDisplay.Display(frame=frame)
337 disp.mtv(exposure, title=
"psf determination")
339 badBits = mi.getMask().getPlaneBitMask(self.config.badMaskBits)
340 fluxName = prefs.getPhotfluxRkey()
341 fluxFlagName = fluxName.rpartition(
"_")[0] +
"_flag"
344 for i, psfCandidate
in enumerate(psfCandidateList):
345 source = psfCandidate.getSource()
348 xc, yc = source.getX(), source.getY()
349 if not np.isfinite(xc)
or not np.isfinite(yc):
352 if source.get(fluxFlagName):
355 flux = source.get(fluxName)
356 if flux < 0
or not np.isfinite(flux):
360 pstamp = psfCandidate.getMaskedImage(pixKernelSize, pixKernelSize).clone()
361 except pexExcept.LengthError:
362 self.log.warning(
"Could not get stamp image for psfCandidate: %s with kernel size: %s",
363 psfCandidate, pixKernelSize)
371 sample = psfSet.newSample()
372 sample.setCatindex(catindex)
373 sample.setExtindex(ext)
374 sample.setObjindex(i)
376 imArray = pstamp.getImage().getArray()
377 imArray[np.where(np.bitwise_and(pstamp.getMask().getArray(), badBits))] = \
379 sample.setVig(imArray)
382 sample.setBacknoise2(backnoise2)
386 sample.setFluxrad(sizes[i])
388 for j
in range(psfSet.getNcontext()):
389 sample.setContext(j, float(contextvalp[j][i]))
390 except Exception
as e:
391 self.log.error(
"Exception when processing sample at (%f,%f): %s", xc, yc, e)
394 psfSet.finiSample(sample)
400 with disp.Buffering():
401 disp.dot(
"o", xc, yc, ctype=afwDisplay.CYAN, size=4)
403 if psfSet.getNsample() == 0:
407 for i
in range(psfSet.getNcontext()):
408 cmin = contextvalp[i].min()
409 cmax = contextvalp[i].max()
410 psfSet.setContextScale(i, cmax - cmin)
411 psfSet.setContextOffset(i, (cmin + cmax)/2.0)
421 field = psfex.Field(
"Unknown")
422 field.addExt(exposure.getWcs(), exposure.getWidth(), exposure.getHeight(), psfSet.getNsample())
430 psfex.makeit(fields, sets)
431 psfs = field.getPsfs()
435 for i
in range(sets[0].getNsample()):
436 index = sets[0].getSample(i).getObjindex()
438 good_indices.append(index)
440 if flagKey
is not None:
441 for i, psfCandidate
in enumerate(psfCandidateList):
442 source = psfCandidate.getSource()
443 if i
in good_indices:
444 source.set(flagKey,
True)
446 xpos = np.array(xpos)
447 ypos = np.array(ypos)
448 numGoodStars = len(good_indices)
449 avgX, avgY = np.mean(xpos), np.mean(ypos)
457 if (ndim := psf.getNdim()) == 0:
459 num_available_stars=nCand,
460 num_good_stars=numGoodStars,
461 poly_ndim_initial=psfSet.getNcontext(),
462 poly_ndim_final=ndim,
469 assert psfCellSet
is not None
472 maUtils.showPsfSpatialCells(exposure, psfCellSet, showChi2=
True,
473 symb=
"o", ctype=afwDisplay.YELLOW, ctypeBad=afwDisplay.RED,
474 size=8, display=disp)
476 disp4 = afwDisplay.Display(frame=4)
477 maUtils.showPsfCandidates(exposure, psfCellSet, psf=psf, display=disp4,
478 normalize=normalizeResiduals,
479 showBadCandidates=showBadCandidates)
480 if displayPsfComponents:
481 disp6 = afwDisplay.Display(frame=6)
482 maUtils.showPsf(psf, display=disp6)
484 disp7 = afwDisplay.Display(frame=7)
485 maUtils.showPsfMosaic(exposure, psf, display=disp7, showFwhm=
True)
486 disp.scale(
'linear', 0, 1)
492 if metadata
is not None:
493 metadata[
"spatialFitChi2"] = np.nan
494 metadata[
"numAvailStars"] = nCand
495 metadata[
"numGoodStars"] = numGoodStars
496 metadata[
"avgX"] = avgX
497 metadata[
"avgY"] = avgY
499 return psf, psfCellSet