LSSTApplications  18.1.0
LSSTDataManagementBasePackage
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 # @package lsst.pipe.tasks.
23 import math
24 import sys
25 
26 import numpy as np
27 import astropy.units as u
28 
29 import lsst.pex.config as pexConf
30 import lsst.pipe.base as pipeBase
31 from lsst.afw.image import abMagErrFromFluxErr, makePhotoCalibFromCalibZeroPoint
32 import lsst.afw.table as afwTable
33 from lsst.meas.astrom import DirectMatchTask, DirectMatchConfigWithoutLoader
34 import lsst.afw.display as afwDisplay
35 from lsst.meas.algorithms import getRefFluxField, ReserveSourcesTask
36 from .colorterms import ColortermLibrary
37 
38 __all__ = ["PhotoCalTask", "PhotoCalConfig"]
39 
40 
41 class PhotoCalConfig(pexConf.Config):
42  """Config for PhotoCal"""
43  match = pexConf.ConfigField("Match to reference catalog",
44  DirectMatchConfigWithoutLoader)
45  reserve = pexConf.ConfigurableField(target=ReserveSourcesTask, doc="Reserve sources from fitting")
46  fluxField = pexConf.Field(
47  dtype=str,
48  default="slot_CalibFlux_instFlux",
49  doc=("Name of the source instFlux field to use. The associated flag field\n"
50  "('<name>_flags') will be implicitly included in badFlags."),
51  )
52  applyColorTerms = pexConf.Field(
53  dtype=bool,
54  default=None,
55  doc=("Apply photometric color terms to reference stars? One of:\n"
56  "None: apply if colorterms and photoCatName are not None;\n"
57  " fail if color term data is not available for the specified ref catalog and filter.\n"
58  "True: always apply colorterms; fail if color term data is not available for the\n"
59  " specified reference catalog and filter.\n"
60  "False: do not apply."),
61  optional=True,
62  )
63  sigmaMax = pexConf.Field(
64  dtype=float,
65  default=0.25,
66  doc="maximum sigma to use when clipping",
67  optional=True,
68  )
69  nSigma = pexConf.Field(
70  dtype=float,
71  default=3.0,
72  doc="clip at nSigma",
73  )
74  useMedian = pexConf.Field(
75  dtype=bool,
76  default=True,
77  doc="use median instead of mean to compute zeropoint",
78  )
79  nIter = pexConf.Field(
80  dtype=int,
81  default=20,
82  doc="number of iterations",
83  )
84  colorterms = pexConf.ConfigField(
85  dtype=ColortermLibrary,
86  doc="Library of photometric reference catalog name: color term dict",
87  )
88  photoCatName = pexConf.Field(
89  dtype=str,
90  optional=True,
91  doc=("Name of photometric reference catalog; used to select a color term dict in colorterms."
92  " see also applyColorTerms"),
93  )
94  magErrFloor = pexConf.RangeField(
95  dtype=float,
96  default=0.0,
97  doc="Additional magnitude uncertainty to be added in quadrature with measurement errors.",
98  min=0.0,
99  )
100 
101  def validate(self):
102  pexConf.Config.validate(self)
103  if self.applyColorTerms and self.photoCatName is None:
104  raise RuntimeError("applyColorTerms=True requires photoCatName is non-None")
105  if self.applyColorTerms and len(self.colorterms.data) == 0:
106  raise RuntimeError("applyColorTerms=True requires colorterms be provided")
107 
108  def setDefaults(self):
109  pexConf.Config.setDefaults(self)
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 
119 
125 
126 class PhotoCalTask(pipeBase.Task):
127  r"""!
128 @anchor PhotoCalTask_
129 
130 @brief Calculate the zero point of an exposure given a lsst.afw.table.ReferenceMatchVector.
131 
132 @section pipe_tasks_photocal_Contents Contents
133 
134  - @ref pipe_tasks_photocal_Purpose
135  - @ref pipe_tasks_photocal_Initialize
136  - @ref pipe_tasks_photocal_IO
137  - @ref pipe_tasks_photocal_Config
138  - @ref pipe_tasks_photocal_Debug
139  - @ref pipe_tasks_photocal_Example
140 
141 @section pipe_tasks_photocal_Purpose Description
142 
143 @copybrief PhotoCalTask
144 
145 Calculate an Exposure's zero-point given a set of flux measurements of stars matched to an input catalogue.
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 @note 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 @ref pipe_tasks_photocal_Example
152 
153 @section pipe_tasks_photocal_Initialize Task initialisation
154 
155 @copydoc \_\_init\_\_
156 
157 @section pipe_tasks_photocal_IO Inputs/Outputs to the run method
158 
159 @copydoc run
160 
161 @section pipe_tasks_photocal_Config Configuration parameters
162 
163 See @ref PhotoCalConfig
164 
165 @section pipe_tasks_photocal_Debug Debug variables
166 
167 The @link lsst.pipe.base.cmdLineTask.CmdLineTask command line task@endlink interface supports a
168 flag @c -d to import @b debug.py from your @c PYTHONPATH; see @ref baseDebug for more about @b debug.py files.
169 
170 The available variables in PhotoCalTask are:
171 <DL>
172  <DT> @c display
173  <DD> If True enable other debug outputs
174  <DT> @c displaySources
175  <DD> If True, display the exposure on the display's frame 1 and overlay the source catalogue.
176  <DL>
177  <DT> red o
178  <DD> Reserved objects
179  <DT> green o
180  <DD> Objects used in the photometric calibration
181  </DL>
182  <DT> @c scatterPlot
183  <DD> Make a scatter plot of flux v. reference magnitude as a function of reference magnitude.
184  - good objects in blue
185  - rejected objects in red
186  (if @c scatterPlot is 2 or more, prompt to continue after each iteration)
187 </DL>
188 
189 @section pipe_tasks_photocal_Example A complete example of using PhotoCalTask
190 
191 This code is in @link examples/photoCalTask.py@endlink, and can be run as @em e.g.
192 @code
193 examples/photoCalTask.py
194 @endcode
195 @dontinclude photoCalTask.py
196 
197 Import the tasks (there are some other standard imports; read the file for details)
198 @skipline from lsst.pipe.tasks.astrometry
199 @skipline measPhotocal
200 
201 We need to create both our tasks before processing any data as the task constructors
202 can add extra columns to the schema which we get from the input catalogue, @c scrCat:
203 @skipline getSchema
204 
205 Astrometry first:
206 @skip AstrometryTask.ConfigClass
207 @until aTask
208 (that @c filterMap line is because our test code doesn't use a filter that the reference catalogue recognises,
209 so we tell it to use the @c r band)
210 
211 Then photometry:
212 @skip measPhotocal
213 @until pTask
214 
215 If the schema has indeed changed we need to add the new columns to the source table
216 (yes; this should be easier!)
217 @skip srcCat
218 @until srcCat = cat
219 
220 We're now ready to process the data (we could loop over multiple exposures/catalogues using the same
221 task objects):
222 @skip matches
223 @until result
224 
225 We can then unpack and use the results:
226 @skip calib
227 @until np.log
228 
229 <HR>
230 To investigate the @ref pipe_tasks_photocal_Debug, put something like
231 @code{.py}
232  import lsstDebug
233  def DebugInfo(name):
234  di = lsstDebug.getInfo(name) # N.b. lsstDebug.Info(name) would call us recursively
235  if name.endswith(".PhotoCal"):
236  di.display = 1
237 
238  return di
239 
240  lsstDebug.Info = DebugInfo
241 @endcode
242 into your debug.py file and run photoCalTask.py with the @c --debug flag.
243  """
244  ConfigClass = PhotoCalConfig
245  _DefaultName = "photoCal"
246 
247  def __init__(self, refObjLoader, schema=None, **kwds):
248  """!Create the photometric calibration task. See PhotoCalTask.init for documentation
249  """
250  pipeBase.Task.__init__(self, **kwds)
251  self.scatterPlot = None
252  self.fig = None
253  if schema is not None:
254  self.usedKey = schema.addField("calib_photometry_used", type="Flag",
255  doc="set if source was used in photometric calibration")
256  else:
257  self.usedKey = None
258  self.match = DirectMatchTask(config=self.config.match, refObjLoader=refObjLoader,
259  name="match", parentTask=self)
260  self.makeSubtask("reserve", columnName="calib_photometry", schema=schema,
261  doc="set if source was reserved from photometric calibration")
262 
263  def getSourceKeys(self, schema):
264  """Return a struct containing the source catalog keys for fields used
265  by PhotoCalTask.
266 
267 
268  Parameters
269  ----------
270  schema : `lsst.afw.table.schema`
271  Schema of the catalog to get keys from.
272 
273  Returns
274  -------
275  result : `lsst.pipe.base.Struct`
276  Result struct with components:
277 
278  - ``instFlux``: Instrument flux key.
279  - ``instFluxErr``: Instrument flux error key.
280  """
281  instFlux = schema.find(self.config.fluxField).key
282  instFluxErr = schema.find(self.config.fluxField + "Err").key
283  return pipeBase.Struct(instFlux=instFlux, instFluxErr=instFluxErr)
284 
285  @pipeBase.timeMethod
286  def extractMagArrays(self, matches, filterName, sourceKeys):
287  """!Extract magnitude and magnitude error arrays from the given matches.
288 
289  @param[in] matches Reference/source matches, a @link lsst::afw::table::ReferenceMatchVector@endlink
290  @param[in] filterName Name of filter being calibrated
291  @param[in] sourceKeys Struct of source catalog keys, as returned by getSourceKeys()
292 
293  @return Struct containing srcMag, refMag, srcMagErr, refMagErr, and magErr numpy arrays
294  where magErr is an error in the magnitude; the error in srcMag - refMag
295  If nonzero, config.magErrFloor will be added to magErr *only* (not srcMagErr or refMagErr), as
296  magErr is what is later used to determine the zero point.
297  Struct also contains refFluxFieldList: a list of field names of the reference catalog used for fluxes
298  (1 or 2 strings)
299  @note These magnitude arrays are the @em inputs to the photometric calibration, some may have been
300  discarded by clipping while estimating the calibration (https://jira.lsstcorp.org/browse/DM-813)
301  """
302  srcInstFluxArr = np.array([m.second.get(sourceKeys.instFlux) for m in matches])
303  srcInstFluxErrArr = np.array([m.second.get(sourceKeys.instFluxErr) for m in matches])
304  if not np.all(np.isfinite(srcInstFluxErrArr)):
305  # this is an unpleasant hack; see DM-2308 requesting a better solution
306  self.log.warn("Source catalog does not have flux uncertainties; using sqrt(flux).")
307  srcInstFluxErrArr = np.sqrt(srcInstFluxArr)
308 
309  # convert source instFlux from DN to an estimate of nJy
310  referenceFlux = (0*u.ABmag).to_value(u.nJy)
311  srcInstFluxArr = srcInstFluxArr * referenceFlux
312  srcInstFluxErrArr = srcInstFluxErrArr * referenceFlux
313 
314  if not matches:
315  raise RuntimeError("No reference stars are available")
316  refSchema = matches[0].first.schema
317 
318  applyColorTerms = self.config.applyColorTerms
319  applyCTReason = "config.applyColorTerms is %s" % (self.config.applyColorTerms,)
320  if self.config.applyColorTerms is None:
321  # apply color terms if color term data is available and photoCatName specified
322  ctDataAvail = len(self.config.colorterms.data) > 0
323  photoCatSpecified = self.config.photoCatName is not None
324  applyCTReason += " and data %s available" % ("is" if ctDataAvail else "is not")
325  applyCTReason += " and photoRefCat %s provided" % ("is" if photoCatSpecified else "is not")
326  applyColorTerms = ctDataAvail and photoCatSpecified
327 
328  if applyColorTerms:
329  self.log.info("Applying color terms for filterName=%r, config.photoCatName=%s because %s",
330  filterName, self.config.photoCatName, applyCTReason)
331  colorterm = self.config.colorterms.getColorterm(
332  filterName=filterName, photoCatName=self.config.photoCatName, doRaise=True)
333  refCat = afwTable.SimpleCatalog(matches[0].first.schema)
334 
335  # extract the matched refCat as a Catalog for the colorterm code
336  refCat.reserve(len(matches))
337  for x in matches:
338  record = refCat.addNew()
339  record.assign(x.first)
340 
341  refMagArr, refMagErrArr = colorterm.getCorrectedMagnitudes(refCat, filterName)
342  fluxFieldList = [getRefFluxField(refSchema, filt) for filt in (colorterm.primary,
343  colorterm.secondary)]
344  else:
345  # no colorterms to apply
346  self.log.info("Not applying color terms because %s", applyCTReason)
347  colorterm = None
348 
349  fluxFieldList = [getRefFluxField(refSchema, filterName)]
350  fluxField = getRefFluxField(refSchema, filterName)
351  fluxKey = refSchema.find(fluxField).key
352  refFluxArr = np.array([m.first.get(fluxKey) for m in matches])
353 
354  try:
355  fluxErrKey = refSchema.find(fluxField + "Err").key
356  refFluxErrArr = np.array([m.first.get(fluxErrKey) for m in matches])
357  except KeyError:
358  # Reference catalogue may not have flux uncertainties; HACK DM-2308
359  self.log.warn("Reference catalog does not have flux uncertainties for %s; using sqrt(flux).",
360  fluxField)
361  refFluxErrArr = np.sqrt(refFluxArr)
362 
363  refMagArr = u.Quantity(refFluxArr, u.nJy).to_value(u.ABmag)
364  # HACK convert to Jy until we have a replacement for this (DM-16903)
365  refMagErrArr = abMagErrFromFluxErr(refFluxErrArr*1e-9, refFluxArr*1e-9)
366 
367  # compute the source catalog magnitudes and errors
368  srcMagArr = u.Quantity(srcInstFluxArr, u.nJy).to_value(u.ABmag)
369  # Fitting with error bars in both axes is hard
370  # for now ignore reference flux error, but ticket DM-2308 is a request for a better solution
371  # HACK convert to Jy until we have a replacement for this (DM-16903)
372  magErrArr = abMagErrFromFluxErr(srcInstFluxErrArr*1e-9, srcInstFluxArr*1e-9)
373  if self.config.magErrFloor != 0.0:
374  magErrArr = (magErrArr**2 + self.config.magErrFloor**2)**0.5
375 
376  srcMagErrArr = abMagErrFromFluxErr(srcInstFluxErrArr*1e-9, srcInstFluxArr*1e-9)
377 
378  good = np.isfinite(srcMagArr) & np.isfinite(refMagArr)
379 
380  return pipeBase.Struct(
381  srcMag=srcMagArr[good],
382  refMag=refMagArr[good],
383  magErr=magErrArr[good],
384  srcMagErr=srcMagErrArr[good],
385  refMagErr=refMagErrArr[good],
386  refFluxFieldList=fluxFieldList,
387  )
388 
389  @pipeBase.timeMethod
390  def run(self, exposure, sourceCat, expId=0):
391  """!Do photometric calibration - select matches to use and (possibly iteratively) compute
392  the zero point.
393 
394  @param[in] exposure Exposure upon which the sources in the matches were detected.
395  @param[in] sourceCat A catalog of sources to use in the calibration
396  (@em i.e. a list of lsst.afw.table.Match with
397  @c first being of type lsst.afw.table.SimpleRecord and @c second type lsst.afw.table.SourceRecord ---
398  the reference object and matched object respectively).
399  (will not be modified except to set the outputField if requested.).
400 
401  @return Struct of:
402  - photoCalib -- @link lsst::afw::image::PhotoCalib@endlink object containing the calibration
403  - arrays ------ Magnitude arrays returned be PhotoCalTask.extractMagArrays
404  - matches ----- Final ReferenceMatchVector, as returned by PhotoCalTask.selectMatches.
405  - zp ---------- Photometric zero point (mag)
406  - sigma ------- Standard deviation of fit of photometric zero point (mag)
407  - ngood ------- Number of sources used to fit photometric zero point
408 
409  The exposure is only used to provide the name of the filter being calibrated (it may also be
410  used to generate debugging plots).
411 
412  The reference objects:
413  - Must include a field @c photometric; True for objects which should be considered as
414  photometric standards
415  - Must include a field @c flux; the flux used to impose a magnitude limit and also to calibrate
416  the data to (unless a color term is specified, in which case ColorTerm.primary is used;
417  See https://jira.lsstcorp.org/browse/DM-933)
418  - May include a field @c stargal; if present, True means that the object is a star
419  - May include a field @c var; if present, True means that the object is variable
420 
421  The measured sources:
422  - Must include PhotoCalConfig.fluxField; the flux measurement to be used for calibration
423 
424  @throws RuntimeError with the following strings:
425 
426  <DL>
427  <DT> No matches to use for photocal
428  <DD> No matches are available (perhaps no sources/references were selected by the matcher).
429  <DT> No reference stars are available
430  <DD> No matches are available from which to extract magnitudes.
431  </DL>
432  """
433  import lsstDebug
434 
435  display = lsstDebug.Info(__name__).display
436  displaySources = display and lsstDebug.Info(__name__).displaySources
437  self.scatterPlot = display and lsstDebug.Info(__name__).scatterPlot
438 
439  if self.scatterPlot:
440  from matplotlib import pyplot
441  try:
442  self.fig.clf()
443  except Exception:
444  self.fig = pyplot.figure()
445 
446  filterName = exposure.getFilter().getName()
447 
448  # Match sources
449  matchResults = self.match.run(sourceCat, filterName)
450  matches = matchResults.matches
451 
452  reserveResults = self.reserve.run([mm.second for mm in matches], expId=expId)
453  if displaySources:
454  self.displaySources(exposure, matches, reserveResults.reserved)
455  if reserveResults.reserved.sum() > 0:
456  matches = [mm for mm, use in zip(matches, reserveResults.use) if use]
457  if len(matches) == 0:
458  raise RuntimeError("No matches to use for photocal")
459  if self.usedKey is not None:
460  for mm in matches:
461  mm.second.set(self.usedKey, True)
462 
463  # Prepare for fitting
464  sourceKeys = self.getSourceKeys(matches[0].second.schema)
465  arrays = self.extractMagArrays(matches=matches, filterName=filterName, sourceKeys=sourceKeys)
466 
467  # Fit for zeropoint
468  r = self.getZeroPoint(arrays.srcMag, arrays.refMag, arrays.magErr)
469  self.log.info("Magnitude zero point: %f +/- %f from %d stars", r.zp, r.sigma, r.ngood)
470 
471  # Prepare the results
472  flux0 = 10**(0.4*r.zp) # Flux of mag=0 star
473  flux0err = 0.4*math.log(10)*flux0*r.sigma # Error in flux0
474  photoCalib = makePhotoCalibFromCalibZeroPoint(flux0, flux0err)
475 
476  return pipeBase.Struct(
477  photoCalib=photoCalib,
478  arrays=arrays,
479  matches=matches,
480  zp=r.zp,
481  sigma=r.sigma,
482  ngood=r.ngood,
483  )
484 
485  def displaySources(self, exposure, matches, reserved, frame=1):
486  """Display sources we'll use for photocal
487 
488  Sources that will be actually used will be green.
489  Sources reserved from the fit will be red.
490 
491  Parameters
492  ----------
493  exposure : `lsst.afw.image.ExposureF`
494  Exposure to display.
495  matches : `list` of `lsst.afw.table.RefMatch`
496  Matches used for photocal.
497  reserved : `numpy.ndarray` of type `bool`
498  Boolean array indicating sources that are reserved.
499  frame : `int`
500  Frame number for display.
501  """
502  disp = afwDisplay.getDisplay(frame=frame)
503  disp.mtv(exposure, title="photocal")
504  with disp.Buffering():
505  for mm, rr in zip(matches, reserved):
506  x, y = mm.second.getCentroid()
507  ctype = afwDisplay.RED if rr else afwDisplay.GREEN
508  disp.dot("o", x, y, size=4, ctype=ctype)
509 
510  def getZeroPoint(self, src, ref, srcErr=None, zp0=None):
511  """!Flux calibration code, returning (ZeroPoint, Distribution Width, Number of stars)
512 
513  We perform nIter iterations of a simple sigma-clipping algorithm with a couple of twists:
514  1. We use the median/interquartile range to estimate the position to clip around, and the
515  "sigma" to use.
516  2. We never allow sigma to go _above_ a critical value sigmaMax --- if we do, a sufficiently
517  large estimate will prevent the clipping from ever taking effect.
518  3. Rather than start with the median we start with a crude mode. This means that a set of magnitude
519  residuals with a tight core and asymmetrical outliers will start in the core. We use the width of
520  this core to set our maximum sigma (see 2.)
521 
522  @return Struct of:
523  - zp ---------- Photometric zero point (mag)
524  - sigma ------- Standard deviation of fit of zero point (mag)
525  - ngood ------- Number of sources used to fit zero point
526  """
527  sigmaMax = self.config.sigmaMax
528 
529  dmag = ref - src
530 
531  indArr = np.argsort(dmag)
532  dmag = dmag[indArr]
533 
534  if srcErr is not None:
535  dmagErr = srcErr[indArr]
536  else:
537  dmagErr = np.ones(len(dmag))
538 
539  # need to remove nan elements to avoid errors in stats calculation with numpy
540  ind_noNan = np.array([i for i in range(len(dmag))
541  if (not np.isnan(dmag[i]) and not np.isnan(dmagErr[i]))])
542  dmag = dmag[ind_noNan]
543  dmagErr = dmagErr[ind_noNan]
544 
545  IQ_TO_STDEV = 0.741301109252802 # 1 sigma in units of interquartile (assume Gaussian)
546 
547  npt = len(dmag)
548  ngood = npt
549  good = None # set at end of first iteration
550  for i in range(self.config.nIter):
551  if i > 0:
552  npt = sum(good)
553 
554  center = None
555  if i == 0:
556  #
557  # Start by finding the mode
558  #
559  nhist = 20
560  try:
561  hist, edges = np.histogram(dmag, nhist, new=True)
562  except TypeError:
563  hist, edges = np.histogram(dmag, nhist) # they removed new=True around numpy 1.5
564  imode = np.arange(nhist)[np.where(hist == hist.max())]
565 
566  if imode[-1] - imode[0] + 1 == len(imode): # Multiple modes, but all contiguous
567  if zp0:
568  center = zp0
569  else:
570  center = 0.5*(edges[imode[0]] + edges[imode[-1] + 1])
571 
572  peak = sum(hist[imode])/len(imode) # peak height
573 
574  # Estimate FWHM of mode
575  j = imode[0]
576  while j >= 0 and hist[j] > 0.5*peak:
577  j -= 1
578  j = max(j, 0)
579  q1 = dmag[sum(hist[range(j)])]
580 
581  j = imode[-1]
582  while j < nhist and hist[j] > 0.5*peak:
583  j += 1
584  j = min(j, nhist - 1)
585  j = min(sum(hist[range(j)]), npt - 1)
586  q3 = dmag[j]
587 
588  if q1 == q3:
589  q1 = dmag[int(0.25*npt)]
590  q3 = dmag[int(0.75*npt)]
591 
592  sig = (q3 - q1)/2.3 # estimate of standard deviation (based on FWHM; 2.358 for Gaussian)
593 
594  if sigmaMax is None:
595  sigmaMax = 2*sig # upper bound on st. dev. for clipping. multiplier is a heuristic
596 
597  self.log.debug("Photo calibration histogram: center = %.2f, sig = %.2f", center, sig)
598 
599  else:
600  if sigmaMax is None:
601  sigmaMax = dmag[-1] - dmag[0]
602 
603  center = np.median(dmag)
604  q1 = dmag[int(0.25*npt)]
605  q3 = dmag[int(0.75*npt)]
606  sig = (q3 - q1)/2.3 # estimate of standard deviation (based on FWHM; 2.358 for Gaussian)
607 
608  if center is None: # usually equivalent to (i > 0)
609  gdmag = dmag[good]
610  if self.config.useMedian:
611  center = np.median(gdmag)
612  else:
613  gdmagErr = dmagErr[good]
614  center = np.average(gdmag, weights=gdmagErr)
615 
616  q3 = gdmag[min(int(0.75*npt + 0.5), npt - 1)]
617  q1 = gdmag[min(int(0.25*npt + 0.5), npt - 1)]
618 
619  sig = IQ_TO_STDEV*(q3 - q1) # estimate of standard deviation
620 
621  good = abs(dmag - center) < self.config.nSigma*min(sig, sigmaMax) # don't clip too softly
622 
623  # =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
624  if self.scatterPlot:
625  try:
626  self.fig.clf()
627 
628  axes = self.fig.add_axes((0.1, 0.1, 0.85, 0.80))
629 
630  axes.plot(ref[good], dmag[good] - center, "b+")
631  axes.errorbar(ref[good], dmag[good] - center, yerr=dmagErr[good],
632  linestyle='', color='b')
633 
634  bad = np.logical_not(good)
635  if len(ref[bad]) > 0:
636  axes.plot(ref[bad], dmag[bad] - center, "r+")
637  axes.errorbar(ref[bad], dmag[bad] - center, yerr=dmagErr[bad],
638  linestyle='', color='r')
639 
640  axes.plot((-100, 100), (0, 0), "g-")
641  for x in (-1, 1):
642  axes.plot((-100, 100), x*0.05*np.ones(2), "g--")
643 
644  axes.set_ylim(-1.1, 1.1)
645  axes.set_xlim(24, 13)
646  axes.set_xlabel("Reference")
647  axes.set_ylabel("Reference - Instrumental")
648 
649  self.fig.show()
650 
651  if self.scatterPlot > 1:
652  reply = None
653  while i == 0 or reply != "c":
654  try:
655  reply = input("Next iteration? [ynhpc] ")
656  except EOFError:
657  reply = "n"
658 
659  if reply == "h":
660  print("Options: c[ontinue] h[elp] n[o] p[db] y[es]", file=sys.stderr)
661  continue
662 
663  if reply in ("", "c", "n", "p", "y"):
664  break
665  else:
666  print("Unrecognised response: %s" % reply, file=sys.stderr)
667 
668  if reply == "n":
669  break
670  elif reply == "p":
671  import pdb
672  pdb.set_trace()
673  except Exception as e:
674  print("Error plotting in PhotoCal.getZeroPoint: %s" % e, file=sys.stderr)
675 
676  # =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
677 
678  old_ngood = ngood
679  ngood = sum(good)
680  if ngood == 0:
681  msg = "PhotoCal.getZeroPoint: no good stars remain"
682 
683  if i == 0: # failed the first time round -- probably all fell in one bin
684  center = np.average(dmag, weights=dmagErr)
685  msg += " on first iteration; using average of all calibration stars"
686 
687  self.log.warn(msg)
688 
689  return pipeBase.Struct(
690  zp=center,
691  sigma=sig,
692  ngood=len(dmag))
693  elif ngood == old_ngood:
694  break
695 
696  if False:
697  ref = ref[good]
698  dmag = dmag[good]
699  dmagErr = dmagErr[good]
700 
701  dmag = dmag[good]
702  dmagErr = dmagErr[good]
703  zp, weightSum = np.average(dmag, weights=1/dmagErr**2, returned=True)
704  sigma = np.sqrt(1.0/weightSum)
705  return pipeBase.Struct(
706  zp=zp,
707  sigma=sigma,
708  ngood=len(dmag),
709  )
def __init__(self, refObjLoader, schema=None, kwds)
Create the photometric calibration task.
Definition: photoCal.py:247
Angle abs(Angle const &a)
Definition: Angle.h:106
std::shared_ptr< PhotoCalib > makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err)
Construct a PhotoCalib from the deprecated Calib-style instFluxMag0/instFluxMag0Err values...
Definition: PhotoCalib.cc:644
def run(self, exposure, sourceCat, expId=0)
Do photometric calibration - select matches to use and (possibly iteratively) compute the zero point...
Definition: photoCal.py:390
double abMagErrFromFluxErr(double fluxErr, double flux)
Compute AB magnitude error from flux and flux error in Janskys.
Definition: Calib.h:52
def extractMagArrays(self, matches, filterName, sourceKeys)
Extract magnitude and magnitude error arrays from the given matches.
Definition: photoCal.py:286
Calculate the zero point of an exposure given a lsst.afw.table.ReferenceMatchVector.
Definition: photoCal.py:126
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations...
int min
def show(frame=None)
Definition: ds9.py:89
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:63
def getZeroPoint(self, src, ref, srcErr=None, zp0=None)
Flux calibration code, returning (ZeroPoint, Distribution Width, Number of stars) ...
Definition: photoCal.py:510
int max
def displaySources(self, exposure, matches, reserved, frame=1)
Definition: photoCal.py:485
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...