LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
psfexStarSelector.py
Go to the documentation of this file.
1# This file is part of meas_extensions_psfex.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21
22__all__ = ("PsfexStarSelectorConfig", "PsfexStarSelectorTask")
23
24import re
25import sys
26
27import numpy as np
28try:
29 import matplotlib.pyplot as plt
30 fig = None
31except ImportError:
32 plt = None
33
34from lsst.pipe.base import Struct
35import lsst.pex.config as pexConfig
36import lsst.afw.display as afwDisplay
37from lsst.meas.algorithms import BaseSourceSelectorTask, sourceSelectorRegistry
38from . import psfexLib
39from .psfex import compute_fwhmrange
40
41
42class PsfexStarSelectorConfig(BaseSourceSelectorTask.ConfigClass):
43 fluxName = pexConfig.Field(
44 dtype=str,
45 doc="Name of photometric flux key ",
46 default="base_PsfFlux",
47 )
48 fluxErrName = pexConfig.Field(
49 dtype=str,
50 doc="Name of phot. flux err. key",
51 default="",
52 )
53 minFwhm = pexConfig.Field(
54 dtype=float,
55 doc="Maximum allowed FWHM ",
56 default=2,
57 )
58 maxFwhm = pexConfig.Field(
59 dtype=float,
60 doc="Minimum allowed FWHM ",
61 default=10,
62 )
63 maxFwhmVariability = pexConfig.Field(
64 dtype=float,
65 doc="Allowed FWHM variability (1.0 = 100%)",
66 default=0.2,
67 )
68 maxellip = pexConfig.Field(
69 dtype=float,
70 doc="Maximum (A-B)/(A+B) ",
71 default=0.3,
72 check=lambda x: x >= 0.0,
73 )
74 minsn = pexConfig.Field(
75 dtype=float,
76 doc="Minimum S/N for candidates",
77 default=100,
78 check=lambda x: x >= 0.0,
79 )
80
81 def validate(self):
82 pexConfig.Config.validate(self)
83
84 if self.fluxErrName == "":
85 self.fluxErrName = self.fluxName + ".err"
86 elif self.fluxErrName != self.fluxName + ".err":
87 msg = f"fluxErrName ({self.fluxErrName}) doesn't correspond to fluxName ({self.fluxName})"
88 raise pexConfig.FieldValidationError(PsfexStarSelectorConfig.fluxName, self, msg)
89
90 if self.minFwhm > self.maxFwhm:
91 raise pexConfig.FieldValidationError(PsfexStarSelectorConfig.minFwhm, self,
92 f"minFwhm ({self.minFwhm}) > maxFwhm ({self.maxFwhm})")
93
94 def setDefaults(self):
95 self.badFlags = [
96 "base_PixelFlags_flag_edge",
97 "base_PixelFlags_flag_saturatedCenter",
98 "base_PixelFlags_flag_crCenter",
99 "base_PixelFlags_flag_bad",
100 "base_PixelFlags_flag_suspectCenter",
101 "base_PsfFlux_flag",
102 # "parent", # actually this is a test on deblend_nChild
103 ]
104
105
107 """A class to handle key strokes with matplotlib displays
108 """
109
110 def __init__(self, axes, xs, ys, x, y, frames=[0]):
111 self.axes = axes
112 self.xs = xs
113 self.ys = ys
114 self.x = x
115 self.y = y
116 self.frames = frames
117
118 self.cid = self.axes.figure.canvas.mpl_connect('key_press_event', self)
119
120 def __call__(self, ev):
121 if ev.inaxes != self.axes:
122 return
123
124 if ev.key and ev.key in ("p"):
125 dist = np.hypot(self.xs - ev.xdata, self.ys - ev.ydata)
126 dist[np.where(np.isnan(dist))] = 1e30
127
128 which = np.where(dist == min(dist))
129
130 x = self.x[which][0]
131 y = self.y[which][0]
132 for frame in self.frames:
133 disp = afwDisplay.Display(frame=frame)
134 disp.pan(x, y)
135 disp.flush()
136 else:
137 pass
138
139
140def plot(mag, width, centers, clusterId, marker="o", markersize=2, markeredgewidth=0, ltype='-',
141 clear=True):
142
143 global fig
144 if not fig:
145 fig = plt.figure()
146 newFig = True
147 else:
148 newFig = False
149 if clear:
150 fig.clf()
151
152 axes = fig.add_axes((0.1, 0.1, 0.85, 0.80))
153
154 xmin = sorted(mag)[int(0.05*len(mag))]
155 xmax = sorted(mag)[int(0.95*len(mag))]
156
157 axes.set_xlim(-17.5, -13)
158 axes.set_xlim(xmin - 0.1*(xmax - xmin), xmax + 0.1*(xmax - xmin))
159 axes.set_ylim(0, 10)
160
161 colors = ["r", "g", "b", "c", "m", "k", ]
162 for k, mean in enumerate(centers):
163 if k == 0:
164 axes.plot(axes.get_xlim(), (mean, mean,), "k%s" % ltype)
165
166 ll = (clusterId == k)
167 axes.plot(mag[ll], width[ll], marker, markersize=markersize, markeredgewidth=markeredgewidth,
168 color=colors[k%len(colors)])
169
170 ll = (clusterId == -1)
171 axes.plot(mag[ll], width[ll], marker, markersize=markersize, markeredgewidth=markeredgewidth,
172 color='k')
173
174 if newFig:
175 axes.set_xlabel("model")
176 axes.set_ylabel(r"$\sqrt{I_{xx} + I_{yy}}$")
177
178 return fig
179
180
186
187
188@pexConfig.registerConfigurable("psfex", sourceSelectorRegistry)
190 r"""A star selector whose algorithm is not yet documented.
191
192 @anchor PsfexStarSelectorTask_
193
194 @section meas_extensions_psfex_psfexStarSelectorStarSelector_Contents Contents
195
196 - @ref meas_extensions_psfex_psfexStarSelectorStarSelector_Purpose
197 - @ref meas_extensions_psfex_psfexStarSelectorStarSelector_Initialize
198 - @ref meas_extensions_psfex_psfexStarSelectorStarSelector_IO
199 - @ref meas_extensions_psfex_psfexStarSelectorStarSelector_Config
200 - @ref meas_extensions_psfex_psfexStarSelectorStarSelector_Debug
201
202 @section meas_extensions_psfex_psfexStarSelectorStarSelector_Purpose Description
203
204 A star selector whose algorithm is not yet documented
205
206 @section meas_extensions_psfex_psfexStarSelectorStarSelector_Initialize Task initialisation
207
208 @copydoc \_\_init\_\_
209
210 @section meas_extensions_psfex_psfexStarSelectorStarSelector_IO Invoking the Task
211
212 Like all star selectors, the main method is `run`.
213
214 @section meas_extensions_psfex_psfexStarSelectorStarSelector_Config Configuration parameters
215
216 See @ref PsfexStarSelectorConfig
217
218 @section meas_extensions_psfex_psfexStarSelectorStarSelector_Debug Debug variables
219
220 PsfexStarSelectorTask has a debug dictionary with the following keys:
221 <dl>
222 <dt>display
223 <dd>bool; if True display debug information
224 <dt>displayExposure
225 <dd>bool; if True display the exposure and spatial cells
226 <dt>plotFwhmHistogram
227 <dd>bool; if True plot histogram of FWHM
228 <dt>plotFlags
229 <dd>bool: if True plot the sources coloured by their flags
230 <dt>plotRejection
231 <dd>bool; if True plot why sources are rejected
232 </dl>
233
234 For example, put something like:
235 @code{.py}
236 import lsstDebug
237 def DebugInfo(name):
238 di = lsstDebug.getInfo(name) # N.b. lsstDebug.Info(name) would call us recursively
239 if name.endswith("objectSizeStarSelector"):
240 di.display = True
241 di.displayExposure = True
242 di.plotFwhmHistogram = True
243
244 return di
245
246 lsstDebug.Info = DebugInfo
247 @endcode
248 into your `debug.py` file and run your task with the `--debug` flag.
249 """ # noqa: W505
250 ConfigClass = PsfexStarSelectorConfig
251 usesMatches = False # selectStars does not use its matches argument
252
253 def selectSources(self, sourceCat, matches=None, exposure=None):
254 """Return a selection of psf-like objects.
255
256 Parameters
257 ----------
258 sourceCat : `lsst.afw.table.SourceCatalog`
259 Catalog of sources to select from.
260 This catalog must be contiguous in memory.
261 matches : `list` of `lsst.afw.table.ReferenceMatch` or None
262 Ignored by this source selector.
263 exposure : `lsst.afw.image.Exposure` or None
264 The exposure the catalog was built from; used for debug display.
265
266 Return
267 ------
268 struct : `lsst.pipe.base.Struct`
269 The struct contains the following data:
270
271 - selected : `numpy.ndarray` of `bool``
272 Boolean array of sources that were selected, same length as
273 sourceCat.
274 """
275 import lsstDebug
276 display = lsstDebug.Info(__name__).display
277
278 displayExposure = display and \
279 lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
280 plotFwhmHistogram = display and plt and \
281 lsstDebug.Info(__name__).plotFwhmHistogram # Plot histogram of FWHM
282 plotFlags = display and plt and \
283 lsstDebug.Info(__name__).plotFlags # Plot the sources coloured by their flags
284 plotRejection = display and plt and \
285 lsstDebug.Info(__name__).plotRejection # Plot why sources are rejected
286 afwDisplay.setDefaultMaskTransparency(75)
287
288 fluxName = self.config.fluxName
289 fluxErrName = self.config.fluxErrName
290 minFwhm = self.config.minFwhm
291 maxFwhm = self.config.maxFwhm
292 maxFwhmVariability = self.config.maxFwhmVariability
293 maxellip = self.config.maxellip
294 minsn = self.config.minsn
295
296 maxelong = (maxellip + 1.0)/(1.0 - maxellip) if maxellip < 1.0 else 100
297
298 # Unpack the catalogue
299 shape = sourceCat.getShapeDefinition()
300 ixx = sourceCat.get("%s.xx" % shape)
301 iyy = sourceCat.get("%s.yy" % shape)
302
303 fwhm = 2*np.sqrt(2*np.log(2))*np.sqrt(0.5*(ixx + iyy))
304 elong = 0.5*(ixx - iyy)/(ixx + iyy)
305
306 flux = sourceCat.get(fluxName)
307 fluxErr = sourceCat.get(fluxErrName)
308 sn = flux/np.where(fluxErr > 0, fluxErr, 1)
309 sn[fluxErr <= 0] = -psfexLib.BIG
310
311 flags = 0x0
312 for i, f in enumerate(self.config.badFlags):
313 flags = np.bitwise_or(flags, np.where(sourceCat.get(f), 1 << i, 0))
314 #
315 # Estimate the acceptable range of source widths
316 #
317 good = np.logical_and(sn > minsn, np.logical_not(flags))
318 good = np.logical_and(good, elong < maxelong)
319 good = np.logical_and(good, fwhm >= minFwhm)
320 good = np.logical_and(good, fwhm < maxFwhm)
321
322 fwhmMode, fwhmMin, fwhmMax = compute_fwhmrange(fwhm[good], maxFwhmVariability, minFwhm, maxFwhm,
323 plot=dict(fwhmHistogram=plotFwhmHistogram))
324
325 # -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
326 #
327 # Here's select_candidates
328 #
329 # ---- Apply some selection over flags, fluxes...
330
331 bad = (flags != 0)
332 # set.setBadFlags(int(sum(bad)))
333
334 if plotRejection:
335 selectionVectors = []
336 selectionVectors.append((bad, "flags %d" % sum(bad)))
337
338 dbad = sn < minsn
339 # set.setBadSN(int(sum(dbad)))
340 bad = np.logical_or(bad, dbad)
341 if plotRejection:
342 selectionVectors.append((dbad, "S/N %d" % sum(dbad)))
343
344 dbad = fwhm < fwhmMin
345 # set.setBadFrmin(int(sum(dbad)))
346 bad = np.logical_or(bad, dbad)
347 if plotRejection:
348 selectionVectors.append((dbad, "fwhmMin %d" % sum(dbad)))
349
350 dbad = fwhm > fwhmMax
351 # set.setBadFrmax(int(sum(dbad)))
352 bad = np.logical_or(bad, dbad)
353 if plotRejection:
354 selectionVectors.append((dbad, "fwhmMax %d" % sum(dbad)))
355
356 dbad = elong > maxelong
357 # set.setBadElong(int(sum(dbad)))
358 bad = np.logical_or(bad, dbad)
359 if plotRejection:
360 selectionVectors.append((dbad, "elong %d" % sum(dbad)))
361
362 good = np.logical_not(bad)
363 #
364 # We know enough to plot, if so requested
365 #
366 frame = 0
367 if displayExposure:
368 mi = exposure.getMaskedImage()
369 disp = afwDisplay.Display(frame=frame)
370 disp.mtv(mi, title="PSF candidates")
371
372 with disp.Buffering():
373 for i, source in enumerate(sourceCat):
374 if good[i]:
375 ctype = afwDisplay.GREEN # star candidate
376 else:
377 ctype = afwDisplay.RED # not star
378
379 disp.dot("+", source.getX() - mi.getX0(), source.getY() - mi.getY0(),
380 ctype=ctype)
381
382 if plotFlags or plotRejection:
383 imag = -2.5*np.log10(flux)
384 plt.clf()
385
386 alpha = 0.5
387 if plotFlags:
388 isSet = np.where(flags == 0x0)[0]
389 plt.plot(imag[isSet], fwhm[isSet], 'o', alpha=alpha, label="good")
390
391 for i, f in enumerate(self.config.badFlags):
392 mask = 1 << i
393 isSet = np.where(np.bitwise_and(flags, mask))[0]
394 if isSet.any():
395 if np.isfinite(imag[isSet] + fwhm[isSet]).any():
396 label = re.sub(r"\_flag", "",
397 re.sub(r"^base\_", "",
398 re.sub(r"^.*base\_PixelFlags\_flag\_", "", f)))
399 plt.plot(imag[isSet], fwhm[isSet], 'o', alpha=alpha, label=label)
400 else:
401 for bad, label in selectionVectors:
402 plt.plot(imag[bad], fwhm[bad], 'o', alpha=alpha, label=label)
403
404 plt.plot(imag[good], fwhm[good], 'o', color="black", label="selected")
405 [plt.axhline(_, color='red') for _ in [fwhmMin, fwhmMax]]
406 plt.xlim(np.median(imag[good]) + 5*np.array([-1, 1]))
407 plt.ylim(fwhm[np.where(np.isfinite(fwhm + imag))].min(), 2*fwhmMax)
408 plt.legend(loc=2)
409 plt.xlabel("Instrumental %s Magnitude" % fluxName.split(".")[-1].title())
410 plt.ylabel("fwhm")
411 title = "PSFEX Star Selection"
412 plt.title("%s %d selected" % (title, sum(good)))
413
414 if displayExposure:
415 global eventHandler
416 eventHandler = EventHandler(plt.axes(), imag, fwhm, sourceCat.getX(), sourceCat.getY(),
417 frames=[frame])
418
419 if plotFlags or plotRejection:
420 while True:
421 try:
422 reply = input("continue? [y[es] h(elp) p(db) q(uit)] ").strip()
423 except EOFError:
424 reply = "y"
425
426 if not reply:
427 reply = "y"
428
429 if reply[0] == "h":
430 print("""\
431At this prompt, you can continue with almost any key; 'p' enters pdb,
432 'q' returns to the shell, and
433 'h' prints this text
434""", end=' ')
435
436 if displayExposure:
437 print("""
438If you put the cursor on a point in the matplotlib scatter plot and hit 'p' you'll see it in ds9.""")
439 elif reply[0] == "p":
440 import pdb
441 pdb.set_trace()
442 elif reply[0] == 'q':
443 sys.exit(1)
444 else:
445 break
446
447 return Struct(selected=good)
int min
bool strip
Definition fits.cc:930
plot(mag, width, centers, clusterId, marker="o", markersize=2, markeredgewidth=0, ltype='-', clear=True)