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