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