LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
photoCal.py
Go to the documentation of this file.
1 from __future__ import absolute_import, division
2 #
3 # LSST Data Management System
4 # Copyright 2008, 2009, 2010 LSST Corporation.
5 #
6 # This product includes software developed by the
7 # LSST Project (http://www.lsst.org/).
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 LSST License Statement and
20 # the GNU General Public License along with this program. If not,
21 # see <http://www.lsstcorp.org/LegalNotices/>.
22 #
23 # \package lsst.pipe.tasks.
24 from itertools import izip
25 import math
26 import sys
27 
28 import numpy as np
29 
30 import lsst.pex.config as pexConf
31 import lsst.pipe.base as pipeBase
32 import lsst.afw.table as afwTable
33 from lsst.afw.image import abMagFromFlux, abMagErrFromFluxErr, fluxFromABMag, Calib
34 import lsst.afw.display.ds9 as ds9
35 from lsst.meas.algorithms import getRefFluxField
36 from .colorterms import ColortermLibrary
37 
38 __all__ = ["PhotoCalTask", "PhotoCalConfig"]
39 
40 def checkSourceFlags(source, sourceKeys):
41  """!Return True if the given source has all good flags set and none of the bad flags set.
42 
43  \param[in] source SourceRecord object to process.
44  \param[in] sourceKeys Struct of source catalog keys, as returned by PhotCalTask.getSourceKeys()
45  """
46  for k in sourceKeys.goodFlags:
47  if not source.get(k): return False
48  if source.getPsfFluxFlag(): return False
49  for k in sourceKeys.badFlags:
50  if source.get(k): return False
51  return True
52 
53 class PhotoCalConfig(pexConf.Config):
54 
55  magLimit = pexConf.Field(dtype=float, doc="Don't use objects fainter than this magnitude", default=22.0)
56  doWriteOutput = pexConf.Field(
57  doc= "Write a field name astrom_usedByPhotoCal to the schema",
58  dtype=bool,
59  default=True,
60  )
61  fluxField = pexConf.Field(
62  doc="Name of the source flux field to use. The associated flag field\n"\
63  "('<name>.flags') will be implicitly included in badFlags.\n",
64  dtype=str,
65  default="slot_CalibFlux_flux",
66  )
67  applyColorTerms = pexConf.Field(
68  doc= "Apply photometric color terms to reference stars? One of: " + \
69  "None: apply if colorterms and photoCatName are not None; " + \
70  "fail if color term data is not available for the specified ref catalog and filter. " + \
71  "True: always apply colorterms; fail if color term data is not available for the " + \
72  "specified reference catalog and filter. " + \
73  "False: do not apply.",
74  dtype=bool,
75  default=None,
76  optional=True,
77  )
78  goodFlags = pexConf.ListField(
79  doc="List of source flag fields that must be set for a source to be used.",
80  dtype=str,
81  default=[],
82  )
83  badFlags = pexConf.ListField(
84  doc="List of source flag fields that will cause a source to be rejected when they are set.",
85  dtype=str,
86  default=["base_PixelFlags_flag_edge", "base_PixelFlags_flag_interpolated",
87  "base_PixelFlags_flag_saturated"],
88  )
89  sigmaMax = pexConf.Field(
90  doc="maximum sigma to use when clipping",
91  dtype=float,
92  default=0.25,
93  optional=True,
94  )
95  nSigma = pexConf.Field(
96  doc="clip at nSigma",
97  dtype=float,
98  default=3.0,
99  )
100  useMedian = pexConf.Field(
101  doc="use median instead of mean to compute zeropoint",
102  dtype=bool,
103  default=True,
104  )
105  nIter = pexConf.Field(
106  doc="number of iterations",
107  dtype=int,
108  default=20,
109  )
110  colorterms = pexConf.ConfigField(
111  doc="Library of photometric reference catalog name: color term dict",
112  dtype=ColortermLibrary,
113  )
114  photoCatName = pexConf.Field(
115  doc="Name of photometric reference catalog; used to select a color term dict in colorterms." + \
116  " see also applyColorTerms",
117  dtype=str,
118  optional=True,
119  )
120  magErrFloor = pexConf.RangeField(
121  doc="Additional magnitude uncertainty to be added in quadrature with measurement errors.",
122  dtype=float,
123  default=0.0,
124  min=0.0
125  )
126 
127  def validate(self):
128  pexConf.Config.validate(self)
129  if self.applyColorTerms and self.photoCatName is None:
130  raise RuntimeError("applyColorTerms=True requires photoCatName is non-None")
131  if self.applyColorTerms and len(self.colorterms.data) == 0:
132  raise RuntimeError("applyColorTerms=True requires colorterms be provided")
133 
134 
135 ## \addtogroup LSST_task_documentation
136 ## \{
137 ## \page photoCalTask
138 ## \ref PhotoCalTask_ "PhotoCalTask"
139 ## Detect positive and negative sources on an exposure and return a new SourceCatalog.
140 ## \}
141 
142 class PhotoCalTask(pipeBase.Task):
143  """!
144 \anchor PhotoCalTask_
145 
146 \brief Calculate the zero point of an exposure given a lsst.afw.table.ReferenceMatchVector.
147 
148 \section pipe_tasks_photocal_Contents Contents
149 
150  - \ref pipe_tasks_photocal_Purpose
151  - \ref pipe_tasks_photocal_Initialize
152  - \ref pipe_tasks_photocal_IO
153  - \ref pipe_tasks_photocal_Config
154  - \ref pipe_tasks_photocal_Debug
155  - \ref pipe_tasks_photocal_Example
156 
157 \section pipe_tasks_photocal_Purpose Description
158 
159 \copybrief PhotoCalTask
160 
161 Calculate an Exposure's zero-point given a set of flux measurements of stars matched to an input catalogue.
162 The type of flux to use is specified by PhotoCalConfig.fluxField.
163 
164 The algorithm clips outliers iteratively, with parameters set in the configuration.
165 
166 \note This task can adds fields to the schema, so any code calling this task must ensure that
167 these columns are indeed present in the input match list; see \ref pipe_tasks_photocal_Example
168 
169 \section pipe_tasks_photocal_Initialize Task initialisation
170 
171 \copydoc \_\_init\_\_
172 
173 \section pipe_tasks_photocal_IO Inputs/Outputs to the run method
174 
175 \copydoc run
176 
177 \section pipe_tasks_photocal_Config Configuration parameters
178 
179 See \ref PhotoCalConfig
180 
181 \section pipe_tasks_photocal_Debug Debug variables
182 
183 The \link lsst.pipe.base.cmdLineTask.CmdLineTask command line task\endlink interface supports a
184 flag \c -d to import \b debug.py from your \c PYTHONPATH; see \ref baseDebug for more about \b debug.py files.
185 
186 The available variables in PhotoCalTask are:
187 <DL>
188  <DT> \c display
189  <DD> If True enable other debug outputs
190  <DT> \c displaySources
191  <DD> If True, display the exposure on ds9's frame 1 and overlay the source catalogue:
192  <DL>
193  <DT> red x
194  <DD> Bad objects
195  <DT> blue +
196  <DD> Matched objects deemed unsuitable for photometric calibration.
197  Additional information is:
198  - a cyan o for galaxies
199  - a magenta o for variables
200  <DT> magenta *
201  <DD> Objects that failed the flux cut
202  <DT> green o
203  <DD> Objects used in the photometric calibration
204  </DL>
205  <DT> \c scatterPlot
206  <DD> Make a scatter plot of flux v. reference magnitude as a function of reference magnitude.
207  - good objects in blue
208  - rejected objects in red
209  (if \c scatterPlot is 2 or more, prompt to continue after each iteration)
210 </DL>
211 
212 \section pipe_tasks_photocal_Example A complete example of using PhotoCalTask
213 
214 This code is in \link examples/photoCalTask.py\endlink, and can be run as \em e.g.
215 \code
216 examples/photoCalTask.py
217 \endcode
218 \dontinclude photoCalTask.py
219 
220 Import the tasks (there are some other standard imports; read the file for details)
221 \skipline from lsst.pipe.tasks.astrometry
222 \skipline measPhotocal
223 
224 We need to create both our tasks before processing any data as the task constructors
225 can add extra columns to the schema which we get from the input catalogue, \c scrCat:
226 \skipline getSchema
227 
228 Astrometry first:
229 \skip AstrometryTask.ConfigClass
230 \until aTask
231 (that \c filterMap line is because our test code doesn't use a filter that the reference catalogue recognises,
232 so we tell it to use the \c r band)
233 
234 Then photometry:
235 \skip measPhotocal
236 \until pTask
237 
238 If the schema has indeed changed we need to add the new columns to the source table
239 (yes; this should be easier!)
240 \skip srcCat
241 \until srcCat = cat
242 
243 We're now ready to process the data (we could loop over multiple exposures/catalogues using the same
244 task objects):
245 \skip matches
246 \until result
247 
248 We can then unpack and use the results:
249 \skip calib
250 \until np.log
251 
252 <HR>
253 To investigate the \ref pipe_tasks_photocal_Debug, put something like
254 \code{.py}
255  import lsstDebug
256  def DebugInfo(name):
257  di = lsstDebug.getInfo(name) # N.b. lsstDebug.Info(name) would call us recursively
258  if name.endswith(".PhotoCal"):
259  di.display = 1
260 
261  return di
262 
263  lsstDebug.Info = DebugInfo
264 \endcode
265 into your debug.py file and run photoCalTask.py with the \c --debug flag.
266  """
267  ConfigClass = PhotoCalConfig
268  _DefaultName = "photoCal"
269 
270  def __init__(self, schema, **kwds):
271  """!Create the photometric calibration task. See PhotoCalTask.init for documentation
272  """
273  pipeBase.Task.__init__(self, **kwds)
274  self.scatterPlot = None
275  self.fig = None
276  if self.config.doWriteOutput:
277  self.outputField = schema.addField("photocal_photometricStandard", type="Flag",
278  doc="set if source was used in photometric calibration")
279  else:
280  self.outputField = None
281 
282  def getSourceKeys(self, schema):
283  """!Return a struct containing the source catalog keys for fields used by PhotoCalTask.
284 
285  Returned fields include:
286  - flux
287  - fluxErr
288  - goodFlags: a list of keys for field names in self.config.goodFlags
289  - badFlags: a list of keys for field names in self.config.badFlags
290  """
291  goodFlags = [schema.find(name).key for name in self.config.goodFlags]
292  flux = schema.find(self.config.fluxField).key
293  fluxErr = schema.find(self.config.fluxField + "Sigma").key
294  badFlags = [schema.find(name).key for name in self.config.badFlags]
295 
296  return pipeBase.Struct(flux=flux, fluxErr=fluxErr, goodFlags=goodFlags, badFlags=badFlags)
297 
298  @pipeBase.timeMethod
299  def selectMatches(self, matches, sourceKeys, filterName, frame=None):
300  """!Select reference/source matches according the criteria specified in the config.
301 
302  \param[in] matches ReferenceMatchVector (not modified)
303  \param[in] sourceKeys Struct of source catalog keys, as returned by getSourceKeys()
304  \param[in] filterName name of camera filter; used to obtain the reference flux field
305  \param[in] frame ds9 frame number to use for debugging display
306  if frame is non-None, display information about trimmed objects on that ds9 frame:
307  - Bad: red x
308  - Unsuitable objects: blue + (and a cyan o if a galaxy)
309  - Failed flux cut: magenta *
310 
311  \return a \link lsst.afw.table.ReferenceMatchVector\endlink that contains only the selected matches.
312  If a schema was passed during task construction, a flag field will be set on sources
313  in the selected matches.
314 
315  \throws ValueError There are no valid matches.
316  """
317 
318  self.log.logdebug("Number of input matches: %d" % (len(matches)))
319  if len(matches) == 0:
320  raise ValueError("No input matches")
321 
322  # Only use stars for which the flags indicate the photometry is good.
323  afterFlagCutInd = [i for i, m in enumerate(matches) if checkSourceFlags(m.second, sourceKeys)]
324  afterFlagCut = [matches[i] for i in afterFlagCutInd]
325  self.log.logdebug("Number of matches after source flag cuts: %d" % (len(afterFlagCut)))
326 
327  if len(afterFlagCut) != len(matches):
328  if frame is not None:
329  with ds9.Buffering():
330  for i, m in enumerate(matches):
331  if i not in afterFlagCutInd:
332  x, y = m.second.getCentroid()
333  ds9.dot("x", x, y, size=4, frame=frame, ctype=ds9.RED)
334 
335  matches = afterFlagCut
336 
337  if len(matches) == 0:
338  raise ValueError("All matches eliminated by source flags")
339 
340  refSchema = matches[0].first.schema
341  try:
342  photometricKey = refSchema.find("photometric").key
343  try:
344  resolvedKey = refSchema.find("resolved").key
345  except:
346  resolvedKey = None
347 
348  try:
349  variableKey = refSchema.find("variable").key
350  except:
351  variableKey = None
352  except:
353  self.log.warn("No 'photometric' flag key found in reference schema.")
354  photometricKey = None
355 
356  if photometricKey is not None:
357  afterRefCutInd = [i for i, m in enumerate(matches) if m.first.get(photometricKey)]
358  afterRefCut = [matches[i] for i in afterRefCutInd]
359 
360  if len(afterRefCut) != len(matches):
361  if frame is not None:
362  with ds9.Buffering():
363  for i, m in enumerate(matches):
364  if i not in afterRefCutInd:
365  x, y = m.second.getCentroid()
366  ds9.dot("+", x, y, size=4, frame=frame, ctype=ds9.BLUE)
367 
368  if resolvedKey and m.first.get(resolvedKey):
369  ds9.dot("o", x, y, size=6, frame=frame, ctype=ds9.CYAN)
370  if variableKey and m.first.get(variableKey):
371  ds9.dot("o", x, y, size=6, frame=frame, ctype=ds9.MAGENTA)
372 
373  matches = afterRefCut
374 
375  self.log.logdebug("Number of matches after reference catalog cuts: %d" % (len(matches)))
376  if len(matches) == 0:
377  raise RuntimeError("No sources remain in match list after reference catalog cuts.")
378  fluxName = getRefFluxField(refSchema, filterName)
379  fluxKey = refSchema.find(fluxName).key
380  if self.config.magLimit is not None:
381  fluxLimit = fluxFromABMag(self.config.magLimit)
382 
383  afterMagCutInd = [i for i, m in enumerate(matches) if (m.first.get(fluxKey) > fluxLimit
384  and m.second.getPsfFlux() > 0.0)]
385  else:
386  afterMagCutInd = [i for i, m in enumerate(matches) if m.second.getPsfFlux() > 0.0]
387 
388  afterMagCut = [matches[i] for i in afterMagCutInd]
389 
390  if len(afterMagCut) != len(matches):
391  if frame is not None:
392  with ds9.Buffering():
393  for i, m in enumerate(matches):
394  if i not in afterMagCutInd:
395  x, y = m.second.getCentroid()
396  ds9.dot("*", x, y, size=4, frame=frame, ctype=ds9.MAGENTA)
397 
398  matches = afterMagCut
399 
400  self.log.logdebug("Number of matches after magnitude limit cuts: %d" % (len(matches)))
401 
402  if len(matches) == 0:
403  raise RuntimeError("No sources remaining in match list after magnitude limit cuts.")
404 
405  if frame is not None:
406  with ds9.Buffering():
407  for m in matches:
408  x, y = m.second.getCentroid()
409  ds9.dot("o", x, y, size=4, frame=frame, ctype=ds9.GREEN)
410 
412  for m in matches:
413  if self.outputField is not None:
414  m.second.set(self.outputField, True)
415  result.append(m)
416  return result
417 
418  @pipeBase.timeMethod
419  def extractMagArrays(self, matches, filterName, sourceKeys):
420  """!Extract magnitude and magnitude error arrays from the given matches.
421 
422  \param[in] matches Reference/source matches, a \link lsst::afw::table::ReferenceMatchVector\endlink
423  \param[in] filterName Name of filter being calibrated
424  \param[in] sourceKeys Struct of source catalog keys, as returned by getSourceKeys()
425 
426  \return Struct containing srcMag, refMag, srcMagErr, refMagErr, and magErr numpy arrays
427  where magErr is an error in the magnitude; the error in srcMag - refMag
428  If nonzero, config.magErrFloor will be added to magErr *only* (not srcMagErr or refMagErr), as
429  magErr is what is later used to determine the zero point.
430  Struct also contains refFluxFieldList: a list of field names of the reference catalog used for fluxes
431  (1 or 2 strings)
432  \note These magnitude arrays are the \em inputs to the photometric calibration, some may have been
433  discarded by clipping while estimating the calibration (https://jira.lsstcorp.org/browse/DM-813)
434  """
435  srcFluxArr = np.array([m.second.get(sourceKeys.flux) for m in matches])
436  srcFluxErrArr = np.array([m.second.get(sourceKeys.fluxErr) for m in matches])
437  if not np.all(np.isfinite(srcFluxErrArr)):
438  # this is an unpleasant hack; see DM-2308 requesting a better solution
439  self.log.warn("Source catalog does not have flux uncertainties; using sqrt(flux).")
440  srcFluxErrArr = np.sqrt(srcFluxArr)
441 
442  # convert source flux from DN to an estimate of Jy
443  JanskysPerABFlux = 3631.0
444  srcFluxArr = srcFluxArr * JanskysPerABFlux
445  srcFluxErrArr = srcFluxErrArr * JanskysPerABFlux
446 
447  if not matches:
448  raise RuntimeError("No reference stars are available")
449  refSchema = matches[0].first.schema
450 
451  applyColorTerms = self.config.applyColorTerms
452  applyCTReason = "config.applyColorTerms is %s" % (self.config.applyColorTerms,)
453  if self.config.applyColorTerms is None:
454  # apply color terms if color term data is available and photoCatName specified
455  ctDataAvail = len(self.config.colorterms.data) > 0
456  photoCatSpecified = self.config.photoCatName is not None
457  applyCTReason += " and data %s available" % ("is" if ctDataAvail else "is not")
458  applyCTReason += " and photoRefCat %s None" % ("is not" if photoCatSpecified else "is")
459  applyColorTerms = ctDataAvail and photoCatSpecified
460 
461  if applyColorTerms:
462  self.log.info("Applying color terms for filterName=%r, config.photoCatName=%s because %s" %
463  (filterName, self.config.photoCatName, applyCTReason))
464  ct = self.config.colorterms.getColorterm(
465  filterName=filterName, photoCatName=self.config.photoCatName, doRaise=True)
466  else:
467  self.log.info("Not applying color terms because %s" % (applyCTReason,))
468  ct = None
469 
470  if ct: # we have a color term to worry about
471  fluxFieldList = [getRefFluxField(refSchema, filt) for filt in (ct.primary, ct.secondary)]
472  missingFluxFieldList = []
473  for fluxField in fluxFieldList:
474  try:
475  refSchema.find(fluxField).key
476  except KeyError:
477  missingFluxFieldList.append(fluxField)
478 
479  if missingFluxFieldList:
480  self.log.warn("Source catalog does not have fluxes for %s; ignoring color terms" %
481  " ".join(missingFluxFieldList))
482  ct = None
483 
484  if not ct:
485  fluxFieldList = [getRefFluxField(refSchema, filterName)]
486 
487  refFluxArrList = [] # list of ref arrays, one per flux field
488  refFluxErrArrList = [] # list of ref flux arrays, one per flux field
489  for fluxField in fluxFieldList:
490  fluxKey = refSchema.find(fluxField).key
491  refFluxArr = np.array([m.first.get(fluxKey) for m in matches])
492  try:
493  fluxErrKey = refSchema.find(fluxField + "Sigma").key
494  refFluxErrArr = np.array([m.first.get(fluxErrKey) for m in matches])
495  except KeyError:
496  # Reference catalogue may not have flux uncertainties; HACK
497  self.log.warn("Reference catalog does not have flux uncertainties for %s; using sqrt(flux)."
498  % fluxField)
499  refFluxErrArr = np.sqrt(refFluxArr)
500 
501  refFluxArrList.append(refFluxArr)
502  refFluxErrArrList.append(refFluxErrArr)
503 
504  if ct: # we have a color term to worry about
505  refMagArr1 = np.array([abMagFromFlux(rf1) for rf1 in refFluxArrList[0]]) # primary
506  refMagArr2 = np.array([abMagFromFlux(rf2) for rf2 in refFluxArrList[1]]) # secondary
507 
508  refMagArr = ct.transformMags(refMagArr1, refMagArr2)
509  refFluxErrArr = ct.propagateFluxErrors(refFluxErrArrList[0], refFluxErrArrList[1])
510  else:
511  refMagArr = np.array([abMagFromFlux(rf) for rf in refFluxArrList[0]])
512 
513  srcMagArr = np.array([abMagFromFlux(sf) for sf in srcFluxArr])
514 
515  # Fitting with error bars in both axes is hard
516  # for now ignore reference flux error, but ticket DM-2308 is a request for a better solution
517  magErrArr = np.array([abMagErrFromFluxErr(fe, sf) for fe, sf in izip(srcFluxErrArr, srcFluxArr)])
518  if self.config.magErrFloor != 0.0:
519  magErrArr = (magErrArr**2 + self.config.magErrFloor**2)**0.5
520 
521  srcMagErrArr = np.array([abMagErrFromFluxErr(sfe, sf) for sfe, sf in izip(srcFluxErrArr, srcFluxArr)])
522  refMagErrArr = np.array([abMagErrFromFluxErr(rfe, rf) for rfe, rf in izip(refFluxErrArr, refFluxArr)])
523 
524  return pipeBase.Struct(
525  srcMag = srcMagArr,
526  refMag = refMagArr,
527  magErr = magErrArr,
528  srcMagErr = srcMagErrArr,
529  refMagErr = refMagErrArr,
530  refFluxFieldList = fluxFieldList,
531  )
532 
533  @pipeBase.timeMethod
534  def run(self, exposure, matches):
535  """!Do photometric calibration - select matches to use and (possibly iteratively) compute
536  the zero point.
537 
538  \param[in] exposure Exposure upon which the sources in the matches were detected.
539  \param[in] matches Input lsst.afw.table.ReferenceMatchVector
540  (\em i.e. a list of lsst.afw.table.Match with
541  \c first being of type lsst.afw.table.SimpleRecord and \c second type lsst.afw.table.SourceRecord ---
542  the reference object and matched object respectively).
543  (will not be modified except to set the outputField if requested.).
544 
545  \return Struct of:
546  - calib ------- \link lsst::afw::image::Calib\endlink object containing the zero point
547  - arrays ------ Magnitude arrays returned be PhotoCalTask.extractMagArrays
548  - matches ----- Final ReferenceMatchVector, as returned by PhotoCalTask.selectMatches.
549 
550  The exposure is only used to provide the name of the filter being calibrated (it may also be
551  used to generate debugging plots).
552 
553  The reference objects:
554  - Must include a field \c photometric; True for objects which should be considered as
555  photometric standards
556  - Must include a field \c flux; the flux used to impose a magnitude limit and also to calibrate
557  the data to (unless a color term is specified, in which case ColorTerm.primary is used;
558  See https://jira.lsstcorp.org/browse/DM-933)
559  - May include a field \c stargal; if present, True means that the object is a star
560  - May include a field \c var; if present, True means that the object is variable
561 
562  The measured sources:
563  - Must include PhotoCalConfig.fluxField; the flux measurement to be used for calibration
564 
565  \throws RuntimeError with the following strings:
566 
567  <DL>
568  <DT> `sources' schema does not contain the calibration object flag "XXX"`
569  <DD> The constructor added fields to the schema that aren't in the Sources
570  <DT> No input matches
571  <DD> The input match vector is empty
572  <DT> All matches eliminated by source flags
573  <DD> The flags specified by \c badFlags in the config eliminated all candidate objects
574  <DT> No sources remain in match list after reference catalog cuts
575  <DD> The reference catalogue has a column "photometric", but no matched objects have it set
576  <DT> No sources remaining in match list after magnitude limit cuts
577  <DD> All surviving matches are either too faint in the catalogue or have negative or \c NaN flux
578  <DT> No reference stars are available
579  <DD> No matches survive all the checks
580  </DL>
581  """
582  import lsstDebug
583 
584  display = lsstDebug.Info(__name__).display
585  displaySources = display and lsstDebug.Info(__name__).displaySources
586  self.scatterPlot = display and lsstDebug.Info(__name__).scatterPlot
587 
588  if self.scatterPlot:
589  from matplotlib import pyplot
590  try:
591  self.fig.clf()
592  except:
593  self.fig = pyplot.figure()
594 
595  if displaySources:
596  frame = 1
597  ds9.mtv(exposure, frame=frame, title="photocal")
598  else:
599  frame = None
600 
601  filterName = exposure.getFilter().getName()
602  sourceKeys = self.getSourceKeys(matches[0].second.schema)
603  matches = self.selectMatches(matches=matches, sourceKeys=sourceKeys, filterName=filterName,
604  frame=frame)
605  arrays = self.extractMagArrays(matches=matches, filterName=filterName, sourceKeys=sourceKeys)
606 
607  if matches and self.outputField:
608  try:
609  # matches[].second is a measured source, wherein we wish to set outputField.
610  # Check that the field is present in the Sources schema.
611  matches[0].second.getSchema().find(self.outputField)
612  except:
613  raise RuntimeError("sources' schema does not contain the used-in-calibration flag \"%s\"" %
614  self.config.outputField)
615 
616  # Fit for zeropoint. We can run the code more than once, so as to
617  # give good stars that got clipped by a bad first guess a second
618  # chance.
619 
620  calib = Calib()
621  zp = None # initial guess
622  r = self.getZeroPoint(arrays.srcMag, arrays.refMag, arrays.magErr, zp0=zp)
623  zp = r.zp
624  self.log.info("Magnitude zero point: %f +/- %f from %d stars" % (r.zp, r.sigma, r.ngood))
625 
626  flux0 = 10**(0.4*r.zp) # Flux of mag=0 star
627  flux0err = 0.4*math.log(10)*flux0*r.sigma # Error in flux0
628 
629  calib.setFluxMag0(flux0, flux0err)
630 
631  return pipeBase.Struct(
632  calib = calib,
633  arrays = arrays,
634  matches = matches,
635  zp = r.zp,
636  sigma = r.sigma,
637  ngood = r.ngood,
638  )
639 
640  def getZeroPoint(self, src, ref, srcErr=None, zp0=None):
641  """!Flux calibration code, returning (ZeroPoint, Distribution Width, Number of stars)
642 
643  We perform nIter iterations of a simple sigma-clipping algorithm with a couple of twists:
644  1. We use the median/interquartile range to estimate the position to clip around, and the
645  "sigma" to use.
646  2. We never allow sigma to go _above_ a critical value sigmaMax --- if we do, a sufficiently
647  large estimate will prevent the clipping from ever taking effect.
648  3. Rather than start with the median we start with a crude mode. This means that a set of magnitude
649  residuals with a tight core and asymmetrical outliers will start in the core. We use the width of
650  this core to set our maximum sigma (see 2.)
651  """
652  sigmaMax = self.config.sigmaMax
653 
654  dmag = ref - src
655 
656  indArr = np.argsort(dmag)
657  dmag = dmag[indArr]
658 
659  if srcErr is not None:
660  dmagErr = srcErr[indArr]
661  else:
662  dmagErr = np.ones(len(dmag))
663 
664  # need to remove nan elements to avoid errors in stats calculation with numpy
665  ind_noNan = np.array([i for i in range(len(dmag))
666  if (not np.isnan(dmag[i]) and not np.isnan(dmagErr[i])) ])
667  dmag = dmag[ind_noNan]
668  dmagErr = dmagErr[ind_noNan]
669 
670  IQ_TO_STDEV = 0.741301109252802; # 1 sigma in units of interquartile (assume Gaussian)
671 
672  npt = len(dmag)
673  ngood = npt
674  good = None # set at end of first iteration
675  for i in range(self.config.nIter):
676  if i > 0:
677  npt = sum(good)
678 
679  center = None
680  if i == 0:
681  #
682  # Start by finding the mode
683  #
684  nhist = 20
685  try:
686  hist, edges = np.histogram(dmag, nhist, new=True)
687  except TypeError:
688  hist, edges = np.histogram(dmag, nhist) # they removed new=True around numpy 1.5
689  imode = np.arange(nhist)[np.where(hist == hist.max())]
690 
691  if imode[-1] - imode[0] + 1 == len(imode): # Multiple modes, but all contiguous
692  if zp0:
693  center = zp0
694  else:
695  center = 0.5*(edges[imode[0]] + edges[imode[-1] + 1])
696 
697  peak = sum(hist[imode])/len(imode) # peak height
698 
699  # Estimate FWHM of mode
700  j = imode[0]
701  while j >= 0 and hist[j] > 0.5*peak:
702  j -= 1
703  j = max(j, 0)
704  q1 = dmag[sum(hist[range(j)])]
705 
706  j = imode[-1]
707  while j < nhist and hist[j] > 0.5*peak:
708  j += 1
709  j = min(j, nhist - 1)
710  j = min(sum(hist[range(j)]), npt - 1)
711  q3 = dmag[j]
712 
713  if q1 == q3:
714  q1 = dmag[int(0.25*npt)]
715  q3 = dmag[int(0.75*npt)]
716 
717  sig = (q3 - q1)/2.3 # estimate of standard deviation (based on FWHM; 2.358 for Gaussian)
718 
719  if sigmaMax is None:
720  sigmaMax = 2*sig # upper bound on st. dev. for clipping. multiplier is a heuristic
721 
722  self.log.logdebug("Photo calibration histogram: center = %.2f, sig = %.2f"
723  % (center, sig))
724 
725  else:
726  if sigmaMax is None:
727  sigmaMax = dmag[-1] - dmag[0]
728 
729  center = np.median(dmag)
730  q1 = dmag[int(0.25*npt)]
731  q3 = dmag[int(0.75*npt)]
732  sig = (q3 - q1)/2.3 # estimate of standard deviation (based on FWHM; 2.358 for Gaussian)
733 
734  if center is None: # usually equivalent to (i > 0)
735  gdmag = dmag[good]
736  if self.config.useMedian:
737  center = np.median(gdmag)
738  else:
739  gdmagErr = dmagErr[good]
740  center = np.average(gdmag, weights=gdmagErr)
741 
742  q3 = gdmag[min(int(0.75*npt + 0.5), npt - 1)]
743  q1 = gdmag[min(int(0.25*npt + 0.5), npt - 1)]
744 
745  sig = IQ_TO_STDEV*(q3 - q1) # estimate of standard deviation
746 
747  good = abs(dmag - center) < self.config.nSigma*min(sig, sigmaMax) # don't clip too softly
748 
749  #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
750  if self.scatterPlot:
751  try:
752  self.fig.clf()
753 
754  axes = self.fig.add_axes((0.1, 0.1, 0.85, 0.80));
755 
756  axes.plot(ref[good], dmag[good] - center, "b+")
757  axes.errorbar(ref[good], dmag[good] - center, yerr=dmagErr[good],
758  linestyle='', color='b')
759 
760  bad = np.logical_not(good)
761  if len(ref[bad]) > 0:
762  axes.plot(ref[bad], dmag[bad] - center, "r+")
763  axes.errorbar(ref[bad], dmag[bad] - center, yerr=dmagErr[bad],
764  linestyle='', color='r')
765 
766  axes.plot((-100, 100), (0, 0), "g-")
767  for x in (-1, 1):
768  axes.plot((-100, 100), x*0.05*np.ones(2), "g--")
769 
770  axes.set_ylim(-1.1, 1.1)
771  axes.set_xlim(24, 13)
772  axes.set_xlabel("Reference")
773  axes.set_ylabel("Reference - Instrumental")
774 
775  self.fig.show()
776 
777  if self.scatterPlot > 1:
778  reply = None
779  while i == 0 or reply != "c":
780  try:
781  reply = raw_input("Next iteration? [ynhpc] ")
782  except EOFError:
783  reply = "n"
784 
785  if reply == "h":
786  print >> sys.stderr, "Options: c[ontinue] h[elp] n[o] p[db] y[es]"
787  continue
788 
789  if reply in ("", "c", "n", "p", "y"):
790  break
791  else:
792  print >> sys.stderr, "Unrecognised response: %s" % reply
793 
794  if reply == "n":
795  break
796  elif reply == "p":
797  import pdb; pdb.set_trace()
798  except Exception, e:
799  print >> sys.stderr, "Error plotting in PhotoCal.getZeroPoint: %s" % e
800 
801  #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
802 
803  old_ngood = ngood
804  ngood = sum(good)
805  if ngood == 0:
806  msg = "PhotoCal.getZeroPoint: no good stars remain"
807 
808  if i == 0: # failed the first time round -- probably all fell in one bin
809  center = np.average(dmag, weights=dmagErr)
810  msg += " on first iteration; using average of all calibration stars"
811 
812  self.log.log(self.log.WARN, msg)
813 
814  return pipeBase.Struct(
815  zp = center,
816  sigma = sig,
817  ngood = len(dmag)
818  )
819  elif ngood == old_ngood:
820  break
821 
822  if False:
823  ref = ref[good]
824  dmag = dmag[good]
825  dmagErr = dmagErr[good]
826 
827  dmag = dmag[good]
828  dmagErr = dmagErr[good]
829  zp, weightSum = np.average(dmag, weights=1/dmagErr**2, returned=True)
830  sigma = np.sqrt(1.0/weightSum)
831  return pipeBase.Struct(
832  zp = zp,
833  sigma = sigma,
834  ngood = len(dmag),
835  )
double abMagErrFromFluxErr(double fluxErr, double flux)
Compute AB magnitude error from flux and flux error in Janskys.
Definition: Calib.h:64
Calculate the zero point of an exposure given a lsst.afw.table.ReferenceMatchVector.
Definition: photoCal.py:142
std::vector< ReferenceMatch > ReferenceMatchVector
Definition: fwd.h:97
def checkSourceFlags
Return True if the given source has all good flags set and none of the bad flags set.
Definition: photoCal.py:40
double abMagFromFlux(double flux)
Compute AB magnitude from flux in Janskys.
Definition: Calib.h:59
boost::enable_if< typename ExpressionTraits< Scalar >::IsScalar, Scalar >::type sum(Scalar const &scalar)
Definition: operators.h:1250
double fluxFromABMag(double mag)
Compute flux in Janskys from AB magnitude.
Definition: Calib.h:69
def run
Do photometric calibration - select matches to use and (possibly iteratively) compute the zero point...
Definition: photoCal.py:534
def extractMagArrays
Extract magnitude and magnitude error arrays from the given matches.
Definition: photoCal.py:419
def getRefFluxField
Get name of flux field in schema.
def __init__
Create the photometric calibration task.
Definition: photoCal.py:270
def getZeroPoint
Flux calibration code, returning (ZeroPoint, Distribution Width, Number of stars) ...
Definition: photoCal.py:640
def getSourceKeys
Return a struct containing the source catalog keys for fields used by PhotoCalTask.
Definition: photoCal.py:282
def selectMatches
Select reference/source matches according the criteria specified in the config.
Definition: photoCal.py:299