LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+1b3060144d,g18429d2f64+f642bf4753,g199a45376c+0ba108daf9,g1fd858c14a+2dcf163641,g262e1987ae+7b8c96d2ca,g29ae962dfc+3bd6ecb08a,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53e1a9e7c5,g4595892280+fef73a337f,g47891489e3+2efcf17695,g4d44eb3520+642b70b07e,g53246c7159+8c5ae1fdc5,g67b6fd64d1+2efcf17695,g67fd3c3899+b70e05ef52,g74acd417e5+317eb4c7d4,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+2efcf17695,g8d7436a09f+3be3c13596,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+a4e7b16d9b,g97be763408+ad77d7208f,g9dd6db0277+b70e05ef52,ga681d05dcb+a3f46e7fff,gabf8522325+735880ea63,gac2eed3f23+2efcf17695,gb89ab40317+2efcf17695,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+b70e05ef52,gdab6d2f7ff+317eb4c7d4,gdc713202bf+b70e05ef52,gdfd2d52018+b10e285e0f,ge365c994fd+310e8507c4,ge410e46f29+2efcf17695,geaed405ab2+562b3308c0,gffca2db377+8c5ae1fdc5,w.2025.35
LSST Data Management Base Package
Loading...
Searching...
No Matches
photoCal.py
Go to the documentation of this file.
1# This file is part of pipe_tasks.
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__ = ["PhotoCalTask", "PhotoCalConfig"]
23
24import math
25import sys
26
27import numpy as np
28import astropy.units as u
29
30import lsst.pex.config as pexConf
31import lsst.pipe.base as pipeBase
32from lsst.afw.image import abMagErrFromFluxErr, makePhotoCalibFromCalibZeroPoint
33import lsst.afw.table as afwTable
34from lsst.meas.astrom import DirectMatchTask, DirectMatchConfigWithoutLoader, MatcherFailure
35import lsst.afw.display as afwDisplay
36from lsst.meas.algorithms import getRefFluxField, ReserveSourcesTask
37from lsst.utils.timer import timeMethod
38from .colorterms import ColortermLibrary
39
40
41class PhotoCalConfig(pexConf.Config):
42 """Config for PhotoCal."""
43
44 match = pexConf.ConfigField("Match to reference catalog",
45 DirectMatchConfigWithoutLoader)
46 reserve = pexConf.ConfigurableField(target=ReserveSourcesTask, doc="Reserve sources from fitting")
47 fluxField = pexConf.Field(
48 dtype=str,
49 default="slot_CalibFlux_instFlux",
50 doc=("Name of the source instFlux field to use.\nThe associated flag field "
51 "('<name>_flags') will be implicitly included in badFlags."),
52 )
53 applyColorTerms = pexConf.Field(
54 dtype=bool,
55 default=False,
56 doc=("Apply photometric color terms to reference stars?\n"
57 "`True`: attempt to apply color terms; fail if color term data is "
58 "not available for the specified reference catalog and filter.\n"
59 "`False`: do not apply color terms."),
60 optional=True,
61 )
62 sigmaMax = pexConf.Field(
63 dtype=float,
64 default=0.25,
65 doc="maximum sigma to use when clipping",
66 optional=True,
67 )
68 nSigma = pexConf.Field(
69 dtype=float,
70 default=3.0,
71 doc="clip at nSigma",
72 )
73 useMedian = pexConf.Field(
74 dtype=bool,
75 default=True,
76 doc="use median instead of mean to compute zeropoint",
77 )
78 nIter = pexConf.Field(
79 dtype=int,
80 default=20,
81 doc="number of iterations",
82 )
83 colorterms = pexConf.ConfigField(
84 dtype=ColortermLibrary,
85 doc="Library of photometric reference catalog name: color term dict (see also applyColorTerms).",
86 )
87 photoCatName = pexConf.Field(
88 dtype=str,
89 optional=True,
90 doc=("Name of photometric reference catalog; used to select a color term dict in colorterms.\n"
91 "See also applyColorTerms."),
92 )
93 magErrFloor = pexConf.RangeField(
94 dtype=float,
95 default=0.0,
96 doc="Additional magnitude uncertainty to be added in quadrature with measurement errors.",
97 min=0.0,
98 )
99
100 def validate(self):
101 pexConf.Config.validate(self)
102 if self.applyColorTerms and self.photoCatName is None:
103 raise RuntimeError("applyColorTerms=True requires photoCatName is non-None")
104 if self.applyColorTerms and len(self.colorterms.data) == 0:
105 raise RuntimeError("applyColorTerms=True requires colorterms be provided")
106 if self.fluxField != self.match.sourceSelection.signalToNoise.fluxField:
107 raise RuntimeError("Configured flux field %s does not match source selector field %s",
108 self.fluxField, self.match.sourceSelection.signalToNoise.fluxField)
109 if self.fluxField + "Err" != self.match.sourceSelection.signalToNoise.errField:
110 raise RuntimeError("Configured flux field %sErr does not match source selector error field %s",
111 self.fluxField, self.match.sourceSelection.signalToNoise.errField)
112
113 def setDefaults(self):
114 pexConf.Config.setDefaults(self)
115 self.match.sourceSelection.doRequirePrimary = True
116 self.match.sourceSelection.doFlags = True
117 self.match.sourceSelection.flags.bad = [
118 "base_PixelFlags_flag_edge",
119 "base_PixelFlags_flag_nodata",
120 "base_PixelFlags_flag_interpolated",
121 "base_PixelFlags_flag_saturated",
122 ]
123 self.match.sourceSelection.doUnresolved = True
124 self.match.sourceSelection.doSignalToNoise = True
125 self.match.sourceSelection.signalToNoise.minimum = 10.0
126 self.match.sourceSelection.signalToNoise.fluxField = self.fluxField
127 self.match.sourceSelection.signalToNoise.errField = self.fluxField + "Err"
128
129
130class PhotoCalTask(pipeBase.Task):
131 """Calculate an Exposure's zero-point given a set of flux measurements
132 of stars matched to an input catalogue.
133
134 Parameters
135 ----------
136 refObjLoader : `lsst.meas.algorithms.ReferenceObjectLoader`
137 A reference object loader object; gen3 pipeline tasks will pass `None`
138 and call `match.setRefObjLoader` in `runQuantum`.
139 schema : `lsst.afw.table.Schema`, optional
140 The schema of the detection catalogs used as input to this task.
141 **kwds
142 Additional keyword arguments.
143
144 Notes
145 -----
146 The type of flux to use is specified by PhotoCalConfig.fluxField.
147
148 The algorithm clips outliers iteratively, with parameters set in the configuration.
149
150 This task can adds fields to the schema, so any code calling this task must ensure that
151 these columns are indeed present in the input match list; see `pipe_tasks_photocal_Example`.
152
153 Debugging:
154
155 The available `~lsst.base.lsstDebug` variables in PhotoCalTask are:
156
157 display :
158 If True enable other debug outputs.
159 displaySources :
160 If True, display the exposure on ds9's frame 1 and overlay the source catalogue.
161
162 red o :
163 Reserved objects.
164 green o :
165 Objects used in the photometric calibration.
166
167 scatterPlot :
168 Make a scatter plot of flux v. reference magnitude as a function of reference magnitude:
169
170 - good objects in blue
171 - rejected objects in red
172
173 (if scatterPlot is 2 or more, prompt to continue after each iteration)
174 """
175
176 ConfigClass = PhotoCalConfig
177 _DefaultName = "photoCal"
178
179 def __init__(self, refObjLoader=None, schema=None, **kwds):
180 pipeBase.Task.__init__(self, **kwds)
181 self.scatterPlot = None
182 self.fig = None
183 if schema is not None:
184 self.usedKey = schema.addField("calib_photometry_used", type="Flag",
185 doc="set if source was used in photometric calibration")
186 else:
187 self.usedKey = None
188 self.match = DirectMatchTask(config=self.config.match, refObjLoader=refObjLoader,
189 name="match", parentTask=self)
190 self.makeSubtask("reserve", columnName="calib_photometry", schema=schema,
191 doc="set if source was reserved from photometric calibration")
192
193 def getSourceKeys(self, schema):
194 """Return a struct containing the source catalog keys for fields used
195 by PhotoCalTask.
196
197 Parameters
198 ----------
199 schema : `lsst.afw.table.schema`
200 Schema of the catalog to get keys from.
201
202 Returns
203 -------
204 result : `lsst.pipe.base.Struct`
205 Results as a struct with attributes:
206
207 ``instFlux``
208 Instrument flux key.
209 ``instFluxErr``
210 Instrument flux error key.
211 """
212 instFlux = schema.find(self.config.fluxField).key
213 instFluxErr = schema.find(self.config.fluxField + "Err").key
214 return pipeBase.Struct(instFlux=instFlux, instFluxErr=instFluxErr)
215
216 @timeMethod
217 def extractMagArrays(self, matches, filterLabel, sourceKeys):
218 """Extract magnitude and magnitude error arrays from the given matches.
219
220 Parameters
221 ----------
222 matches : `lsst.afw.table.ReferenceMatchVector`
223 Reference/source matches.
224 filterLabel : `str`
225 Label of filter being calibrated.
226 sourceKeys : `lsst.pipe.base.Struct`
227 Struct of source catalog keys, as returned by getSourceKeys().
228
229 Returns
230 -------
231 result : `lsst.pipe.base.Struct`
232 Results as a struct with attributes:
233
234 ``srcMag``
235 Source magnitude (`np.array`).
236 ``refMag``
237 Reference magnitude (`np.array`).
238 ``srcMagErr``
239 Source magnitude error (`np.array`).
240 ``refMagErr``
241 Reference magnitude error (`np.array`).
242 ``magErr``
243 An error in the magnitude; the error in ``srcMag`` - ``refMag``.
244 If nonzero, ``config.magErrFloor`` will be added to ``magErr`` only
245 (not ``srcMagErr`` or ``refMagErr``), as
246 ``magErr`` is what is later used to determine the zero point (`np.array`).
247 ``refFluxFieldList``
248 A list of field names of the reference catalog used for fluxes (1 or 2 strings) (`list`).
249 """
250 srcInstFluxArr = np.array([m.second.get(sourceKeys.instFlux) for m in matches])
251 srcInstFluxErrArr = np.array([m.second.get(sourceKeys.instFluxErr) for m in matches])
252 if not np.all(np.isfinite(srcInstFluxErrArr)):
253 # this is an unpleasant hack; see DM-2308 requesting a better solution
254 self.log.warning("Source catalog does not have flux uncertainties; using sqrt(flux).")
255 srcInstFluxErrArr = np.sqrt(srcInstFluxArr)
256
257 # convert source instFlux from DN to an estimate of nJy
258 referenceFlux = (0*u.ABmag).to_value(u.nJy)
259 srcInstFluxArr = srcInstFluxArr * referenceFlux
260 srcInstFluxErrArr = srcInstFluxErrArr * referenceFlux
261
262 if not matches:
263 raise MatcherFailure("No reference stars are available")
264 refSchema = matches[0].first.schema
265
266 if self.config.applyColorTerms:
267 self.log.info("Applying color terms for filter=%r, config.photoCatName=%s",
268 filterLabel.physicalLabel, self.config.photoCatName)
269 colorterm = self.config.colorterms.getColorterm(filterLabel.physicalLabel,
270 self.config.photoCatName,
271 doRaise=True)
272 refCat = afwTable.SimpleCatalog(matches[0].first.schema)
273
274 # extract the matched refCat as a Catalog for the colorterm code
275 refCat.reserve(len(matches))
276 for x in matches:
277 record = refCat.addNew()
278 record.assign(x.first)
279
280 refMagArr, refMagErrArr = colorterm.getCorrectedMagnitudes(refCat)
281 fluxFieldList = [getRefFluxField(refSchema, filt) for filt in (colorterm.primary,
282 colorterm.secondary)]
283 else:
284 self.log.info("Not applying color terms.")
285 colorterm = None
286
287 fluxFieldList = [getRefFluxField(refSchema, filterLabel.bandLabel)]
288 fluxField = getRefFluxField(refSchema, filterLabel.bandLabel)
289 fluxKey = refSchema.find(fluxField).key
290 refFluxArr = np.array([m.first.get(fluxKey) for m in matches])
291
292 try:
293 fluxErrKey = refSchema.find(fluxField + "Err").key
294 refFluxErrArr = np.array([m.first.get(fluxErrKey) for m in matches])
295 except KeyError:
296 # Reference catalogue may not have flux uncertainties; HACK DM-2308
297 self.log.warning("Reference catalog does not have flux uncertainties for %s;"
298 " using sqrt(flux).", fluxField)
299 refFluxErrArr = np.sqrt(refFluxArr)
300
301 refMagArr = u.Quantity(refFluxArr, u.nJy).to_value(u.ABmag)
302 # HACK convert to Jy until we have a replacement for this (DM-16903)
303 refMagErrArr = abMagErrFromFluxErr(refFluxErrArr*1e-9, refFluxArr*1e-9)
304
305 # compute the source catalog magnitudes and errors
306 srcMagArr = u.Quantity(srcInstFluxArr, u.nJy).to_value(u.ABmag)
307 # Fitting with error bars in both axes is hard
308 # for now ignore reference flux error, but ticket DM-2308 is a request for a better solution
309 # HACK convert to Jy until we have a replacement for this (DM-16903)
310 magErrArr = abMagErrFromFluxErr(srcInstFluxErrArr*1e-9, srcInstFluxArr*1e-9)
311 if self.config.magErrFloor != 0.0:
312 magErrArr = (magErrArr**2 + self.config.magErrFloor**2)**0.5
313
314 srcMagErrArr = abMagErrFromFluxErr(srcInstFluxErrArr*1e-9, srcInstFluxArr*1e-9)
315
316 good = np.isfinite(srcMagArr) & np.isfinite(refMagArr)
317 if not np.any(good):
318 raise MatcherFailure("No matched reference stars with valid fluxes",
319 goodSrc=int(np.sum(np.isfinite(srcMagArr))),
320 goodRef=int(np.sum(np.isfinite(refMagArr))))
321
322 return pipeBase.Struct(
323 srcMag=srcMagArr[good],
324 refMag=refMagArr[good],
325 magErr=magErrArr[good],
326 srcMagErr=srcMagErrArr[good],
327 refMagErr=refMagErrArr[good],
328 refFluxFieldList=fluxFieldList,
329 )
330
331 @timeMethod
332 def run(self, exposure, sourceCat, expId=0):
333 """Do photometric calibration - select matches to use and (possibly iteratively) compute
334 the zero point.
335
336 Parameters
337 ----------
338 exposure : `lsst.afw.image.Exposure`
339 Exposure upon which the sources in the matches were detected.
340 sourceCat : `lsst.afw.table.SourceCatalog`
341 Good stars selected for use in calibration.
342 expId : `int`, optional
343 Exposure ID.
344
345 Returns
346 -------
347 result : `lsst.pipe.base.Struct`
348 Results as a struct with attributes:
349
350 ``photoCalib``
351 Object containing the zero point (`lsst.afw.image.Calib`).
352 ``arrays``
353 Magnitude arrays returned be `PhotoCalTask.extractMagArrays`.
354 ``matches``
355 ReferenceMatchVector, as returned by the matcher
356 ``matchMeta`` : metadata needed to unpersist matches, as returned
357 by the matcher (`lsst.daf.base.PropertyList`)
358 ``zp``
359 Photometric zero point (mag, `float`).
360 ``sigma``
361 Standard deviation of fit of photometric zero point (mag, `float`).
362 ``ngood``
363 Number of sources used to fit photometric zero point (`int`).
364
365 Raises
366 ------
367 RuntimeError
368 Raised if any of the following occur:
369 - No matches to use for photocal.
370 - No matches are available (perhaps no sources/references were selected by the matcher).
371 - No reference stars are available.
372 - No matches are available from which to extract magnitudes.
373
374 Notes
375 -----
376 The exposure is only used to provide the name of the filter being calibrated (it may also be
377 used to generate debugging plots).
378
379 The reference objects:
380 - Must include a field ``photometric``; True for objects which should be considered as
381 photometric standards.
382 - Must include a field ``flux``; the flux used to impose a magnitude limit and also to calibrate
383 the data to (unless a color term is specified, in which case ColorTerm.primary is used;
384 See https://jira.lsstcorp.org/browse/DM-933).
385 - May include a field ``stargal``; if present, True means that the object is a star.
386 - May include a field ``var``; if present, True means that the object is variable.
387
388 The measured sources:
389 - Must include PhotoCalConfig.fluxField; the flux measurement to be used for calibration.
390 """
391 import lsstDebug
392
393 display = lsstDebug.Info(__name__).display
394 displaySources = display and lsstDebug.Info(__name__).displaySources
395 self.scatterPlot = display and lsstDebug.Info(__name__).scatterPlot
396
397 if self.scatterPlot:
398 from matplotlib import pyplot
399 try:
400 self.fig.clf()
401 except Exception:
402 self.fig = pyplot.figure()
403
404 filterLabel = exposure.getFilter()
405
406 # Match sources
407 matchResults = self.match.run(sourceCat, filterLabel.bandLabel)
408 matches = matchResults.matches
409
410 reserveResults = self.reserve.run([mm.second for mm in matches], expId=expId)
411 if displaySources:
412 self.displaySources(exposure, matches, reserveResults.reserved)
413 if reserveResults.reserved.sum() > 0:
414 matches = [mm for mm, use in zip(matches, reserveResults.use) if use]
415 if len(matches) == 0:
416 raise MatcherFailure("No matches to use for photocal.")
417 if self.usedKey is not None:
418 for mm in matches:
419 mm.second.set(self.usedKey, True)
420
421 # Prepare for fitting
422 sourceKeys = self.getSourceKeys(matches[0].second.schema)
423 arrays = self.extractMagArrays(matches, filterLabel, sourceKeys)
424
425 # Fit for zeropoint
426 r = self.getZeroPoint(arrays.srcMag, arrays.refMag, arrays.magErr)
427 self.log.info("Magnitude zero point: %f +/- %f from %d stars", r.zp, r.sigma, r.ngood)
428
429 # Prepare the results
430 flux0 = 10**(0.4*r.zp) # Flux of mag=0 star
431 flux0err = 0.4*math.log(10)*flux0*r.sigma # Error in flux0
432 photoCalib = makePhotoCalibFromCalibZeroPoint(flux0, flux0err)
433 self.log.info("Photometric calibration factor (nJy/ADU): %f +/- %f",
434 photoCalib.getCalibrationMean(),
435 photoCalib.getCalibrationErr())
436
437 return pipeBase.Struct(
438 photoCalib=photoCalib,
439 arrays=arrays,
440 matches=matches,
441 matchMeta=matchResults.matchMeta,
442 zp=r.zp,
443 sigma=r.sigma,
444 ngood=r.ngood,
445 )
446
447 def displaySources(self, exposure, matches, reserved, frame=1):
448 """Display sources we'll use for photocal.
449
450 Sources that will be actually used will be green.
451 Sources reserved from the fit will be red.
452
453 Parameters
454 ----------
455 exposure : `lsst.afw.image.ExposureF`
456 Exposure to display.
457 matches : `list` of `lsst.afw.table.RefMatch`
458 Matches used for photocal.
459 reserved : `numpy.ndarray` of type `bool`
460 Boolean array indicating sources that are reserved.
461 frame : `int`, optional
462 Frame number for display.
463 """
464 disp = afwDisplay.getDisplay(frame=frame)
465 disp.mtv(exposure, title="photocal")
466 with disp.Buffering():
467 for mm, rr in zip(matches, reserved):
468 x, y = mm.second.getCentroid()
469 ctype = afwDisplay.RED if rr else afwDisplay.GREEN
470 disp.dot("o", x, y, size=4, ctype=ctype)
471
472 def getZeroPoint(self, src, ref, srcErr=None, zp0=None):
473 """Flux calibration code, returning (ZeroPoint, Distribution Width, Number of stars).
474
475 Returns
476 -------
477 result : `lsst.pipe.base.Struct`
478 Results as a struct with attributes:
479
480 ``zp``
481 Photometric zero point (mag, `float`).
482 ``sigma``
483 Standard deviation of fit of photometric zero point (mag, `float`).
484 ``ngood``
485 Number of sources used to fit photometric zero point (`int`).
486
487 Notes
488 -----
489 We perform nIter iterations of a simple sigma-clipping algorithm with a couple of twists:
490 - We use the median/interquartile range to estimate the position to clip around, and the
491 "sigma" to use.
492 - We never allow sigma to go _above_ a critical value sigmaMax --- if we do, a sufficiently
493 large estimate will prevent the clipping from ever taking effect.
494 - Rather than start with the median we start with a crude mode. This means that a set of magnitude
495 residuals with a tight core and asymmetrical outliers will start in the core. We use the width of
496 this core to set our maximum sigma (see second bullet).
497 """
498 sigmaMax = self.config.sigmaMax
499
500 dmag = ref - src
501
502 indArr = np.argsort(dmag)
503 dmag = dmag[indArr]
504
505 if srcErr is not None:
506 dmagErr = srcErr[indArr]
507 else:
508 dmagErr = np.ones(len(dmag))
509
510 # need to remove nan elements to avoid errors in stats calculation with numpy
511 ind_noNan = np.array([i for i in range(len(dmag))
512 if (not np.isnan(dmag[i]) and not np.isnan(dmagErr[i]))])
513 dmag = dmag[ind_noNan]
514 dmagErr = dmagErr[ind_noNan]
515
516 IQ_TO_STDEV = 0.741301109252802 # 1 sigma in units of interquartile (assume Gaussian)
517
518 npt = len(dmag)
519 ngood = npt
520 good = None # set at end of first iteration
521 for i in range(self.config.nIter):
522 if i > 0:
523 npt = sum(good)
524
525 center = None
526 if i == 0:
527 #
528 # Start by finding the mode
529 #
530 nhist = 20
531 try:
532 hist, edges = np.histogram(dmag, nhist, new=True)
533 except TypeError:
534 hist, edges = np.histogram(dmag, nhist) # they removed new=True around numpy 1.5
535 imode = np.arange(nhist)[np.where(hist == hist.max())]
536
537 if imode[-1] - imode[0] + 1 == len(imode): # Multiple modes, but all contiguous
538 if zp0:
539 center = zp0
540 else:
541 center = 0.5*(edges[imode[0]] + edges[imode[-1] + 1])
542
543 peak = sum(hist[imode])/len(imode) # peak height
544
545 # Estimate FWHM of mode
546 j = imode[0]
547 while j >= 0 and hist[j] > 0.5*peak:
548 j -= 1
549 j = max(j, 0)
550 q1 = dmag[sum(hist[range(j)])]
551
552 j = imode[-1]
553 while j < nhist and hist[j] > 0.5*peak:
554 j += 1
555 j = min(j, nhist - 1)
556 j = min(sum(hist[range(j)]), npt - 1)
557 q3 = dmag[j]
558
559 if q1 == q3:
560 q1 = dmag[int(0.25*npt)]
561 q3 = dmag[int(0.75*npt)]
562
563 sig = (q3 - q1)/2.3 # estimate of standard deviation (based on FWHM; 2.358 for Gaussian)
564
565 if sigmaMax is None:
566 sigmaMax = 2*sig # upper bound on st. dev. for clipping. multiplier is a heuristic
567
568 self.log.debug("Photo calibration histogram: center = %.2f, sig = %.2f", center, sig)
569
570 else:
571 if sigmaMax is None:
572 sigmaMax = dmag[-1] - dmag[0]
573
574 center = np.median(dmag)
575 q1 = dmag[int(0.25*npt)]
576 q3 = dmag[int(0.75*npt)]
577 sig = (q3 - q1)/2.3 # estimate of standard deviation (based on FWHM; 2.358 for Gaussian)
578
579 if center is None: # usually equivalent to (i > 0)
580 gdmag = dmag[good]
581 if self.config.useMedian:
582 center = np.median(gdmag)
583 else:
584 gdmagErr = dmagErr[good]
585 center = np.average(gdmag, weights=gdmagErr)
586
587 q3 = gdmag[min(int(0.75*npt + 0.5), npt - 1)]
588 q1 = gdmag[min(int(0.25*npt + 0.5), npt - 1)]
589
590 sig = IQ_TO_STDEV*(q3 - q1) # estimate of standard deviation
591
592 good = abs(dmag - center) < self.config.nSigma*min(sig, sigmaMax) # don't clip too softly
593
594 # =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
595 if self.scatterPlot:
596 try:
597 self.fig.clf()
598
599 axes = self.fig.add_axes((0.1, 0.1, 0.85, 0.80))
600
601 axes.plot(ref[good], dmag[good] - center, "b+")
602 axes.errorbar(ref[good], dmag[good] - center, yerr=dmagErr[good],
603 linestyle='', color='b')
604
605 bad = np.logical_not(good)
606 if len(ref[bad]) > 0:
607 axes.plot(ref[bad], dmag[bad] - center, "r+")
608 axes.errorbar(ref[bad], dmag[bad] - center, yerr=dmagErr[bad],
609 linestyle='', color='r')
610
611 axes.plot((-100, 100), (0, 0), "g-")
612 for x in (-1, 1):
613 axes.plot((-100, 100), x*0.05*np.ones(2), "g--")
614
615 axes.set_ylim(-1.1, 1.1)
616 axes.set_xlim(24, 13)
617 axes.set_xlabel("Reference")
618 axes.set_ylabel("Reference - Instrumental")
619
620 self.fig.show()
621
622 if self.scatterPlot > 1:
623 reply = None
624 while i == 0 or reply != "c":
625 try:
626 reply = input("Next iteration? [ynhpc] ")
627 except EOFError:
628 reply = "n"
629
630 if reply == "h":
631 print("Options: c[ontinue] h[elp] n[o] p[db] y[es]", file=sys.stderr)
632 continue
633
634 if reply in ("", "c", "n", "p", "y"):
635 break
636 else:
637 print("Unrecognised response: %s" % reply, file=sys.stderr)
638
639 if reply == "n":
640 break
641 elif reply == "p":
642 import pdb
643 pdb.set_trace()
644 except Exception as e:
645 print("Error plotting in PhotoCal.getZeroPoint: %s" % e, file=sys.stderr)
646
647 # =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
648
649 old_ngood = ngood
650 ngood = sum(good)
651 if ngood == 0:
652 msg = "PhotoCal.getZeroPoint: no good stars remain"
653
654 if i == 0: # failed the first time round -- probably all fell in one bin
655 center = np.average(dmag, weights=dmagErr)
656 msg += " on first iteration; using average of all calibration stars"
657
658 self.log.warning(msg)
659
660 return pipeBase.Struct(
661 zp=center,
662 sigma=sig,
663 ngood=len(dmag))
664 elif ngood == old_ngood:
665 break
666
667 if False:
668 ref = ref[good]
669 dmag = dmag[good]
670 dmagErr = dmagErr[good]
671
672 dmag = dmag[good]
673 dmagErr = dmagErr[good]
674 zp, weightSum = np.average(dmag, weights=1/dmagErr**2, returned=True)
675 sigma = np.sqrt(1.0/weightSum)
676 return pipeBase.Struct(
677 zp=zp,
678 sigma=sigma,
679 ngood=len(dmag),
680 )
displaySources(self, exposure, matches, reserved, frame=1)
Definition photoCal.py:447
run(self, exposure, sourceCat, expId=0)
Definition photoCal.py:332
extractMagArrays(self, matches, filterLabel, sourceKeys)
Definition photoCal.py:217
getZeroPoint(self, src, ref, srcErr=None, zp0=None)
Definition photoCal.py:472
__init__(self, refObjLoader=None, schema=None, **kwds)
Definition photoCal.py:179