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
psfex.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
22import os
23import re
24
25import numpy as np
26from astropy.io import fits
27try:
28 import matplotlib.pyplot as plt
29except ImportError:
30 plt = None
31
32import lsst.geom as geom
33import lsst.afw.geom as afwGeom
34from lsst.afw.fits import readMetadata
35import lsst.afw.image as afwImage
36import lsst.afw.table as afwTable
37import lsst.afw.display as afwDisplay
38from . import psfexLib
39
40afwDisplay.setDefaultMaskTransparency(75)
41
42
43def compute_fwhmrange(fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False)):
44 """Compute the FWHM range associated to a series of FWHM measurements.
45 AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
46 VERSION 20/03/2008
47
48 Parameters
49 ----------
50 fwhm: iterable of `float`
51 Iterable of full width half-maxima.
52 maxvar: `float`
53 Maximum allowed FWHM variation.
54 minin: `float`
55 Minimum allowed FWHM.
56 maxin: `float`
57 Maximum allowed FWHM.
58 plot: `dict`, optional
59 Dict of plotting options.
60
61 Returns
62 -------
63 fmin: `float`
64 FWHM mode.
65 minout: `float`
66 Lower FWHM range.
67 maxout: `float`
68 Upper FWHM range.
69 """
70
71 nfwhm = len(fwhm)
72 fwhm.sort()
73
74 # Find the mode
75 nw = nfwhm//4
76 if nw < 4:
77 nw = 1
78 dfmin = psfexLib.BIG
79 fmin = 0.0
80 for i in range(nfwhm - nw):
81 df = fwhm[i + nw] - fwhm[i]
82 if df < dfmin:
83 dfmin = df
84 fmin = (fwhm[i + nw] + fwhm[i])/2.0
85
86 if nfwhm < 2:
87 fmin = fwhm[0]
88
89 dfmin = (maxvar + 1.0)**0.3333333
90 minout = fmin/dfmin if dfmin > 0.0 else 0.0
91 if minout < minin:
92 minout = minin
93
94 maxout = fmin*dfmin**2
95 if maxout > maxin:
96 maxout = maxin
97
98 if plt and plot.get("fwhmHistogram"):
99 plt.clf()
100 plt.hist(fwhm, nfwhm//10 + 1, normed=1, facecolor='g', alpha=0.75)
101 plt.xlabel("FWHM")
102 plt.axvline(fmin, color='red')
103 [plt.axvline(_, color='blue') for _ in (minout, maxout)]
104
105 input("Continue? ")
106
107 return fmin, minout, maxout
108
109
110def select_candidates(set, prefs, frmin, frmax,
111 flags, flux, fluxerr, rmsSize, elong, vignet,
112 plot=dict(), title=""):
113 maxbad = prefs.getBadpixNmax()
114 maxbadflag = prefs.getBadpixFlag()
115 maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100.0
116 minsn = prefs.getMinsn()
117
118 sn = flux/np.where(fluxerr > 0, fluxerr, 1)
119 sn[fluxerr <= 0] = -psfexLib.BIG
120 # ---- Apply some selection over flags, fluxes...
121 plotFlags = plot.get("showFlags") if plt else False
122 plotRejection = plot.get("showRejection") if plt else False
123
124 bad = flags & prefs.getFlagMask() != 0
125 set.setBadFlags(int(sum(bad != 0)))
126
127 if plotRejection:
128 selectionVectors = []
129 selectionVectors.append((bad, "flags %d" % sum(bad != 0)))
130
131 dbad = sn < minsn
132 set.setBadSN(int(sum(dbad)))
133 bad = np.logical_or(bad, dbad)
134 if plotRejection:
135 selectionVectors.append((dbad, "S/N %d" % sum(dbad)))
136
137 dbad = rmsSize < frmin
138 set.setBadFrmin(int(sum(dbad)))
139 bad = np.logical_or(bad, dbad)
140 if plotRejection:
141 selectionVectors.append((dbad, "frmin %d" % sum(dbad)))
142
143 dbad = rmsSize > frmax
144 set.setBadFrmax(int(sum(dbad)))
145 bad = np.logical_or(bad, dbad)
146 if plotRejection:
147 selectionVectors.append((dbad, "frmax %d" % sum(dbad)))
148
149 dbad = elong > maxelong
150 set.setBadElong(int(sum(dbad)))
151 bad = np.logical_or(bad, dbad)
152 if plotRejection:
153 selectionVectors.append((dbad, "elong %d" % sum(dbad)))
154
155 # -- ... and check the integrity of the sample
156 if maxbadflag:
157 nbad = np.array([(v <= -psfexLib.BIG).sum() for v in vignet])
158 dbad = nbad > maxbad
159 set.setBadPix(int(sum(dbad)))
160 bad = np.logical_or(bad, dbad)
161 if plotRejection:
162 selectionVectors.append((dbad, "badpix %d" % sum(dbad)))
163
164 good = np.logical_not(bad)
165 if plotFlags or plotRejection:
166 imag = -2.5*np.log10(flux)
167 plt.clf()
168
169 alpha = 0.5
170 if plotFlags:
171 labels = getFlags()
172
173 isSet = np.where(flags == 0x0)[0]
174 plt.plot(imag[isSet], rmsSize[isSet], 'o', alpha=alpha, label="good")
175
176 for i in range(16):
177 mask = 1 << i
178 if mask & prefs.getFlagMask():
179 isSet = np.where(np.bitwise_and(flags, mask))[0]
180 if isSet.any():
181 plt.plot(imag[isSet], rmsSize[isSet], 'o', alpha=alpha, label=labels[mask])
182 else:
183 for bad, label in selectionVectors:
184 plt.plot(imag[bad], rmsSize[bad], 'o', alpha=alpha, label=label)
185
186 plt.plot(imag[good], rmsSize[good], 'o', color="black", label="selected")
187 [plt.axhline(_, color='red') for _ in [frmin, frmax]]
188 plt.xlim(np.median(imag[good]) + 5*np.array([-1, 1]))
189 plt.ylim(-0.1, 2*frmax)
190 plt.legend(loc=2)
191 plt.xlabel("Instrumental Magnitude")
192 plt.ylabel("rmsSize")
193 plt.title("%s %d selected" % (title, sum(good)))
194
195 input("Continue? ")
196
197 return good
198
199
200def showPsf(psf, set, ext=None, wcsData=None, trim=0, nspot=5,
201 diagnostics=False, outDir="", frame=None, title=None):
202 """Show a PSF on display (e.g., ds9)
203 """
204
205 if ext is not None:
206 psf = psf[ext]
207
208 if wcsData:
209 if ext is not None:
210 wcsData = wcsData[ext]
211 wcs, naxis1, naxis2 = wcsData
212 else:
213 wcs, naxis1, naxis2 = None, None, None
214
215 naxis = [naxis1, naxis2]
216 for i in range(2):
217 if naxis[i] is None:
218 # cmin, cmax are the range of input star positions
219 cmin, cmax = [set.getContextOffset(i) + d*set.getContextScale(i) for d in (-0.5, 0.5)]
220 naxis[i] = cmax + cmin # a decent guess
221
222 if naxis[0] > naxis[1]:
223 nx, ny = int(nspot*naxis[0]/float(naxis[1]) + 0.5), nspot
224 else:
225 nx, ny = nspot, int(nspot*naxis[1]/float(naxis[0]) + 0.5)
226
227 mos = afwDisplay.utils.Mosaic(gutter=2, background=0.02)
228
229 xpos, ypos = np.linspace(0, naxis[0], nx), np.linspace(0, naxis[1], ny)
230 for y in ypos:
231 for x in xpos:
232 psf.build(x, y)
233
234 im = afwImage.ImageF(*psf.getLoc().shape)
235 im.getArray()[:] = psf.getLoc()
236 im /= float(im.getArray().max())
237 if trim:
238 if trim > im.getHeight()//2:
239 trim = im.getHeight()//2
240
241 im = im[trim:-trim, trim:-trim]
242
243 mos.append(im)
244
245 mosaic = mos.makeMosaic(mode=nx)
246 if frame is not None:
247 afwDisplay.Display(frame=frame).mtv(mosaic, title=title)
248 #
249 # Figure out the WCS for the mosaic
250 #
251 pos = []
252 pos.append([geom.PointD(0, 0), wcs.pixelToSky(geom.PointD(0, 0))])
253 pos.append([geom.PointD(*mosaic.getDimensions()), wcs.pixelToSky(geom.PointD(naxis1, naxis2))])
254
255 CD = []
256 for i in range(2):
257 delta = pos[1][1][i].asDegrees() - pos[0][1][i].asDegrees()
258 CD.append(delta/(pos[1][0][i] - pos[0][0][i]))
259 CD = np.array(CD)
260 CD.shape = (2, 2)
261 mosWcs = afwGeom.makeSkyWcs(crval=pos[0][0], crpix=pos[0][1], cdMatrix=CD)
262
263 if ext is not None:
264 title = "%s-%d" % (title, ext)
265
266 if frame is not None:
267 afwDisplay.Display(frame=frame).mtv(mosaic, title=title, wcs=mosWcs)
268
269 if diagnostics:
270 outFile = "%s-mod.fits" % title
271 if outDir:
272 outFile = os.path.join(outDir, outFile)
273 mosaic.writeFits(outFile, mosWcs.getFitsMetadata())
274
275 mos = afwDisplay.utils.Mosaic(gutter=4, background=0.002)
276 for i in range(set.getNsample()):
277 s = set.getSample(i)
278 if ext is not None and s.getExtindex() != ext:
279 continue
280
281 smos = afwDisplay.utils.Mosaic(gutter=2, background=-0.003)
282 for func in [s.getVig, s.getVigResi]:
283 arr = func()
284 if func == s.getVig:
285 norm = float(arr.max()) if True else s.getNorm()
286
287 arr /= norm
288 im = afwImage.ImageF(*arr.shape)
289 im.getArray()[:] = arr
290 smos.append(im)
291
292 mos.append(smos.makeMosaic(mode="x"))
293
294 mosaic = mos.makeMosaic(title=title)
295
296 if frame is not None:
297 afwDisplay.Display(frame=frame + 1).mtv(mosaic, title=title)
298
299 if diagnostics:
300 outFile = "%s-psfstars.fits" % title
301 if outDir:
302 outFile = os.path.join(outDir, outFile)
303
304 mosaic.writeFits(outFile)
305
306
307def getFlags(tab=None):
308 flagKeys = [
309 "base_PixelFlags_flag_edge",
310 # "base_PixelFlags_flag_interpolated",
311 # "base_PixelFlags_flag_interpolatedCenter",
312 # "base_PixelFlags_flag_saturated",
313 "base_PixelFlags_flag_saturatedCenter",
314 # "base_PixelFlags_flag_cr",
315 "base_PixelFlags_flag_crCenter",
316 "base_PixelFlags_flag_bad",
317 "base_PsfFlux_flag",
318 "parent",
319 ]
320
321 if tab is None:
322 flags = {}
323 for i, k in enumerate(flagKeys):
324 flags[1 << i] = re.sub(r"\_flag", "",
325 re.sub(r"^base\_", "", re.sub(r"^base\_PixelFlags\_flag\_", "", k)))
326 else:
327 flags = 0
328 for i, k in enumerate(flagKeys):
329 if k == "parent":
330 try:
331 isSet = tab.get("deblend_nChild") > 0
332 except KeyError:
333 isSet = 0
334 else:
335 isSet = tab.get(k)
336 flags = np.bitwise_or(flags, np.where(isSet, 1 << i, 0))
337
338 return flags
339
340
341def guessCalexp(fileName):
342 for guess in [
343 re.sub("/src", r"", fileName),
344 re.sub("(SRC([^.]+))", r"CORR\2-exp", fileName),
345 ]:
346 if guess != fileName and os.path.exists(guess):
347 return guess
348
349 raise RuntimeError("Unable to find a calexp to go with %s" % fileName)
350
351
352def makeitLsst(prefs, context, saveWcs=False, plot=dict()):
353 """This is the python wrapper that reads lsst tables
354 """
355 # Create an array of PSFs (one PSF for each extension)
356 if prefs.getVerboseType() != prefs.QUIET:
357 print("----- %d input catalogues:" % prefs.getNcat())
358
359 if saveWcs: # only needed for making plots
360 wcssList = []
361
362 fields = psfexLib.vectorField()
363 for cat in prefs.getCatalogs():
364 field = psfexLib.Field(cat)
365 wcss = []
366 wcssList.append(wcss)
367 with fits.open(cat):
368 # Hack: I want the WCS so I'll guess where the calexp is to be
369 # found
370 calexpFile = guessCalexp(cat)
371 md = readMetadata(calexpFile)
372 wcs = afwGeom.makeSkyWcs(md)
373
374 if not wcs:
375 cdMatrix = np.array([1.0, 0.0, 0.0, 1.0])
376 cdMatrix.shape = (2, 2)
377 wcs = afwGeom.makeSkyWcs(crpix=geom.PointD(0, 0),
378 crval=geom.SpherePoint(0.0, 0.0, geom.degrees),
379 cdMatrix=cdMatrix)
380
381 naxis1, naxis2 = md.getScalar("NAXIS1"), md.getScalar("NAXIS2")
382 # Find how many rows there are in the catalogue
383 md = readMetadata(cat)
384
385 field.addExt(wcs, naxis1, naxis2, md.getScalar("NAXIS2"))
386 if saveWcs:
387 wcss.append((wcs, naxis1, naxis2))
388
389 field.finalize()
390 fields.append(field)
391
392 fields[0].getNext() # number of extensions
393
394 prefs.getPsfStep()
395
396 sets = psfexLib.vectorSet()
397 for set in load_samplesLsst(prefs, context, plot=plot):
398 sets.append(set)
399
400 psfexLib.makeit(fields, sets)
401
402 ret = [[f.getPsfs() for f in fields], sets]
403 if saveWcs:
404 ret.append(wcssList)
405
406 return ret
407
408
409def read_samplesLsst(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, nobj,
410 plot=dict(showFlags=False, showRejection=False)):
411 # allocate a new set iff set is None
412 if not set:
413 set = psfexLib.Set(context)
414
415 cmin, cmax = None, None
416 if set.getNcontext():
417 cmin = np.empty(set.getNcontext())
418 cmax = np.empty(set.getNcontext())
419 for i in range(set.getNcontext()):
420 if set.getNsample():
421 cmin[i] = set.getContextOffset(i) - set.getContextScale(i)/2.0
422 cmax[i] = cmin[i] + set.getContextScale(i)
423 else:
424 cmin[i] = psfexLib.BIG
425 cmax[i] = -psfexLib.BIG
426 #
427 # Read data
428 #
429 tab = afwTable.SourceCatalog.readFits(filename)
430
431 centroid = tab.getCentroidDefinition()
432 xm = tab.get("%s.x" % centroid)
433 ym = tab.get("%s.y" % centroid)
434
435 shape = tab.getShapeDefinition()
436 ixx = tab.get("%s.xx" % shape)
437 iyy = tab.get("%s.yy" % shape)
438
439 rmsSize = np.sqrt(0.5*(ixx + iyy))
440 elong = 0.5*(ixx - iyy)/(ixx + iyy)
441
442 flux = tab.get(prefs.getPhotfluxRkey())
443 fluxErr = tab.get(prefs.getPhotfluxerrRkey())
444 flags = getFlags(tab)
445
446 #
447 # Now the VIGNET data
448 #
449 vigw, vigh = 35, 35 # [prefs.getPsfsize()[i] for i in range(2)]
450 if set.empty():
451 set.setVigSize(vigw, vigh)
452
453 vignet = np.empty(nobj*vigw*vigh, "float32").reshape(nobj, vigw, vigh)
454
455 # Hack: I want the WCS so I'll guess where the calexp is to be found
456 calexpFile = guessCalexp(filename)
457 mi = afwImage.MaskedImageF(calexpFile)
458 backnoise2 = np.median(mi.getVariance().getArray())
459 gain = 1.0
460
461 edgeBit = [k for k, v in getFlags().items() if v == "edge"][0]
462
463 for i, xc, yc in zip(range(nobj), xm, ym):
464 try:
465 x, y = int(xc), int(yc)
466 except ValueError:
467 flags[i] |= edgeBit # mark star as bad
468
469 try:
470 pstamp = mi[x - vigw//2:x + vigw//2 + 1, y - vigh//2:y + vigh//2 + 1]
471 vignet[i] = pstamp.getImage().getArray().transpose()
472 except Exception:
473 flags[i] |= edgeBit # mark star as bad
474
475 # Try to load the set of context keys
476 pc = 0
477 contextvalp = []
478 for i, key in enumerate(context.getName()):
479 if context.getPcflag(i):
480 contextvalp.append(pcval[pc])
481 pc += 1
482 elif key[0] == ':':
483 try:
484 contextvalp.append(tab.header[key[1:]])
485 except KeyError:
486 raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
487 (key[1:], filename))
488 else:
489 try:
490 contextvalp.append(tab.get(key))
491 except KeyError:
492 raise RuntimeError("*Error*: %s parameter not found in the header of %s" %
493 (key, filename))
494 set.setContextname(i, key)
495
496 # Now examine each vector of the shipment
497 good = select_candidates(set, prefs, frmin, frmax,
498 flags, flux, fluxErr, rmsSize, elong, vignet,
499 plot=plot, title="%s[%d]" % (filename, ext + 1) if next > 1 else filename)
500 #
501 # Insert the good candidates into the set
502 #
503 if not vignet.dtype.isnative:
504 # without the swap setVig fails with
505 # "ValueError: 'unaligned arrays cannot be converted to C++'"
506 vignet = vignet.byteswap()
507
508 for i in np.where(good)[0]:
509 sample = set.newSample()
510 sample.setCatindex(catindex)
511 sample.setExtindex(ext)
512
513 sample.setVig(vignet[i])
514 sample.setNorm(float(flux[i]))
515 sample.setBacknoise2(backnoise2)
516 sample.setGain(gain)
517 sample.setX(float(xm[i]))
518 sample.setY(float(ym[i]))
519 sample.setFluxrad(float(rmsSize[i]))
520
521 for j in range(set.getNcontext()):
522 sample.setContext(j, float(contextvalp[j][i]))
523
524 set.finiSample(sample, prefs.getProfAccuracy())
525
526 # ---- Update min and max
527 for j in range(set.getNcontext()):
528 cmin[j] = contextvalp[j][good].min()
529 cmax[j] = contextvalp[j][good].max()
530
531 # Update the scaling
532 if set.getNsample():
533 for i in range(set.getNcontext()):
534 set.setContextScale(i, cmax[i] - cmin[i])
535 set.setContextOffset(i, (cmin[i] + cmax[i])/2.0)
536
537 # Don't waste memory!
538 set.trimMemory()
539
540 return set
541
542
543def load_samplesLsst(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict()):
544 minsn = prefs.getMinsn()
545 maxelong = (prefs.getMaxellip() + 1.0)/(1.0 - prefs.getMaxellip()) if prefs.getMaxellip() < 1.0 else 100
546 min = prefs.getFwhmrange()[0]
547 max = prefs.getFwhmrange()[1]
548
549 filenames = prefs.getCatalogs()
550
551 ncat = len(filenames)
552 fwhmmin = np.empty(ncat)
553 fwhmmax = np.empty(ncat)
554
555 if not prefs.getAutoselectFlag():
556 fwhmmin = np.zeros(ncat) + prefs.getFwhmrange()[0]
557 fwhmmax = np.zeros(ncat) + prefs.getFwhmrange()[1]
558 fwhmmode = (fwhmmin + fwhmmax)/2.0
559 else:
560 fwhms = {}
561
562 # -- Try to estimate the most appropriate Half-light Radius range
563 # -- Get the Half-light radii
564 nobj = 0
565 for i, fileName in enumerate(filenames):
566 fwhms[i] = []
567
568 if prefs.getVerboseType() != prefs.QUIET:
569 print("Examining Catalog #%d" % (i+1))
570
571 # ---- Read input catalog
572 tab = afwTable.SourceCatalog.readFits(fileName)
573
574 # -------- Fill the FWHM array
575 shape = tab.getShapeDefinition()
576 ixx = tab.get("%s.xx" % shape)
577 iyy = tab.get("%s.yy" % shape)
578
579 rmsSize = np.sqrt(0.5*(ixx + iyy))
580 elong = 0.5*(ixx - iyy)/(ixx + iyy)
581
582 flux = tab.get(prefs.getPhotfluxRkey())
583 fluxErr = tab.get(prefs.getPhotfluxerrRkey())
584
585 flags = getFlags(tab)
586
587 good = np.logical_and(flux/fluxErr > minsn,
588 np.logical_not(flags & prefs.getFlagMask()))
589 good = np.logical_and(good, elong < maxelong)
590 fwhm = 2.0*rmsSize
591 good = np.logical_and(good, fwhm >= min)
592 good = np.logical_and(good, fwhm < max)
593 fwhms[i] = fwhm[good]
594
595 if prefs.getVarType() == prefs.VAR_NONE:
596 if nobj:
597 fwhms_all = np.empty(sum([len(f) for f in fwhms.values()]))
598 i = 0
599 for f in fwhms.values():
600 fwhms_all[i:len(f)] = f
601 i += len(f)
602 mode, min, max = compute_fwhmrange(fwhms_all, prefs.getMaxvar(),
603 prefs.getFwhmrange()[0], prefs.getFwhmrange()[1],
604 plot=plot)
605 else:
606 raise RuntimeError("No source with appropriate FWHM found!!")
607 mode = min = max = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
608
609 fwhmmin = np.zeros(ncat) + min
610 fwhmmax = np.zeros(ncat) + max
611 fwhmmode = np.zeros(ncat) + mode
612 else:
613 fwhmmode = np.empty(ncat)
614 fwhmmin = np.empty(ncat)
615 fwhmmax = np.empty(ncat)
616
617 for i in range(ncat):
618 nobj = len(fwhms[i])
619 if (nobj):
620 fwhmmode[i], fwhmmin[i], fwhmmax[i] = \
621 compute_fwhmrange(fwhms[i], prefs.getMaxvar(),
622 prefs.getFwhmrange()[0], prefs.getFwhmrange()[1], plot=plot)
623 else:
624 raise RuntimeError("No source with appropriate FWHM found!!")
625 fwhmmode[i] = fwhmmin[i] = fwhmmax[i] = 2.35/(1.0 - 1.0/psfexLib.cvar.INTERPFAC)
626
627 # Read the samples
628 mode = psfexLib.BIG # mode of FWHM distribution
629
630 sets = []
631 for i, fileName in enumerate(filenames):
632 set = None
633 for ext in range(next):
634 set = read_samplesLsst(prefs, set, fileName, fwhmmin[i]/2.0, fwhmmax[i]/2.0,
635 ext, next, i, context,
636 context.getPc(i) if context.getNpc() else None, nobj, plot=plot)
637
638 if fwhmmode[i] < mode:
639 mode = fwhmmode[i]
640
641 set.setFwhm(mode)
642
643 if prefs.getVerboseType() != prefs.QUIET:
644 if set.getNsample():
645 print("%d samples loaded." % set.getNsample())
646 else:
647 raise RuntimeError("No appropriate source found!!")
648
649 sets.append(set)
650
651 return sets
std::vector< SchemaItem< Flag > > * items
int min
int max
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition SkyWcs.cc:521
select_candidates(set, prefs, frmin, frmax, flags, flux, fluxerr, rmsSize, elong, vignet, plot=dict(), title="")
Definition psfex.py:112
showPsf(psf, set, ext=None, wcsData=None, trim=0, nspot=5, diagnostics=False, outDir="", frame=None, title=None)
Definition psfex.py:201
read_samplesLsst(prefs, set, filename, frmin, frmax, ext, next, catindex, context, pcval, nobj, plot=dict(showFlags=False, showRejection=False))
Definition psfex.py:410
compute_fwhmrange(fwhm, maxvar, minin, maxin, plot=dict(fwhmHistogram=False))
Definition psfex.py:43
makeitLsst(prefs, context, saveWcs=False, plot=dict())
Definition psfex.py:352
load_samplesLsst(prefs, context, ext=psfexLib.Prefs.ALL_EXTENSIONS, next=1, plot=dict())
Definition psfex.py:543