27 import astropy.units
as u
31 from lsst.afw.image import abMagErrFromFluxErr, makePhotoCalibFromCalibZeroPoint
36 from lsst.utils.timer
import timeMethod
37 from .colorterms
import ColortermLibrary
39 __all__ = [
"PhotoCalTask",
"PhotoCalConfig"]
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(
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."),
53 applyColorTerms = pexConf.Field(
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."),
64 sigmaMax = pexConf.Field(
67 doc=
"maximum sigma to use when clipping",
70 nSigma = pexConf.Field(
75 useMedian = pexConf.Field(
78 doc=
"use median instead of mean to compute zeropoint",
80 nIter = pexConf.Field(
83 doc=
"number of iterations",
85 colorterms = pexConf.ConfigField(
86 dtype=ColortermLibrary,
87 doc=
"Library of photometric reference catalog name: color term dict",
89 photoCatName = pexConf.Field(
92 doc=(
"Name of photometric reference catalog; used to select a color term dict in colorterms."
93 " see also applyColorTerms"),
95 magErrFloor = pexConf.RangeField(
98 doc=
"Additional magnitude uncertainty to be added in quadrature with measurement errors.",
103 pexConf.Config.validate(self)
105 raise RuntimeError(
"applyColorTerms=True requires photoCatName is non-None")
107 raise RuntimeError(
"applyColorTerms=True requires colorterms be provided")
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",
117 self.
matchmatch.sourceSelection.doUnresolved =
True
129 @anchor PhotoCalTask_
131 @brief Calculate the zero point of an exposure given a lsst.afw.table.ReferenceMatchVector.
133 @section pipe_tasks_photocal_Contents Contents
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
142 @section pipe_tasks_photocal_Purpose Description
144 @copybrief PhotoCalTask
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.
149 The algorithm clips outliers iteratively, with parameters set in the configuration.
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
154 @section pipe_tasks_photocal_Initialize Task initialisation
156 @copydoc \_\_init\_\_
158 @section pipe_tasks_photocal_IO Inputs/Outputs to the run method
162 @section pipe_tasks_photocal_Config Configuration parameters
164 See @ref PhotoCalConfig
166 @section pipe_tasks_photocal_Debug Debug variables
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.
171 The available variables in PhotoCalTask are:
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.
179 <DD> Reserved objects
181 <DD> Objects used in the photometric calibration
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)
190 @section pipe_tasks_photocal_Example A complete example of using PhotoCalTask
192 This code is in @link examples/photoCalTask.py@endlink, and can be run as @em e.g.
194 examples/photoCalTask.py
196 @dontinclude photoCalTask.py
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
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:
207 @skip AstrometryTask.ConfigClass
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)
216 If the schema has indeed changed we need to add the new columns to the source table
217 (yes; this should be easier!)
221 We're now ready to process the data (we could loop over multiple exposures/catalogues using the same
226 We can then unpack and use the results:
231 To investigate the @ref pipe_tasks_photocal_Debug, put something like
235 di = lsstDebug.getInfo(name) # N.b. lsstDebug.Info(name) would call us recursively
236 if name.endswith(".PhotoCal"):
241 lsstDebug.Info = DebugInfo
243 into your debug.py file and run photoCalTask.py with the @c --debug flag.
245 ConfigClass = PhotoCalConfig
246 _DefaultName =
"photoCal"
248 def __init__(self, refObjLoader, schema=None, **kwds):
249 """!Create the photometric calibration task. See PhotoCalTask.init for documentation
251 pipeBase.Task.__init__(self, **kwds)
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")
260 name=
"match", parentTask=self)
261 self.makeSubtask(
"reserve", columnName=
"calib_photometry", schema=schema,
262 doc=
"set if source was reserved from photometric calibration")
265 """Return a struct containing the source catalog keys for fields used
271 schema : `lsst.afw.table.schema`
272 Schema of the catalog to get keys from.
276 result : `lsst.pipe.base.Struct`
277 Result struct with components:
279 - ``instFlux``: Instrument flux key.
280 - ``instFluxErr``: Instrument flux error key.
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)
288 """!Extract magnitude and magnitude error arrays from the given matches.
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()
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
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)
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)):
307 self.log.
warning(
"Source catalog does not have flux uncertainties; using sqrt(flux).")
308 srcInstFluxErrArr = np.sqrt(srcInstFluxArr)
311 referenceFlux = (0*u.ABmag).to_value(u.nJy)
312 srcInstFluxArr = srcInstFluxArr * referenceFlux
313 srcInstFluxErrArr = srcInstFluxErrArr * referenceFlux
316 raise RuntimeError(
"No reference stars are available")
317 refSchema = matches[0].first.schema
319 applyColorTerms = self.config.applyColorTerms
320 applyCTReason =
"config.applyColorTerms is %s" % (self.config.applyColorTerms,)
321 if self.config.applyColorTerms
is None:
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
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,
338 refCat.reserve(len(matches))
340 record = refCat.addNew()
341 record.assign(x.first)
343 refMagArr, refMagErrArr = colorterm.getCorrectedMagnitudes(refCat)
344 fluxFieldList = [
getRefFluxField(refSchema, filt)
for filt
in (colorterm.primary,
345 colorterm.secondary)]
348 self.log.
info(
"Not applying color terms because %s", applyCTReason)
353 fluxKey = refSchema.find(fluxField).key
354 refFluxArr = np.array([m.first.get(fluxKey)
for m
in matches])
357 fluxErrKey = refSchema.find(fluxField +
"Err").key
358 refFluxErrArr = np.array([m.first.get(fluxErrKey)
for m
in matches])
361 self.log.
warning(
"Reference catalog does not have flux uncertainties for %s;"
362 " using sqrt(flux).", fluxField)
363 refFluxErrArr = np.sqrt(refFluxArr)
365 refMagArr = u.Quantity(refFluxArr, u.nJy).to_value(u.ABmag)
370 srcMagArr = u.Quantity(srcInstFluxArr, u.nJy).to_value(u.ABmag)
375 if self.config.magErrFloor != 0.0:
376 magErrArr = (magErrArr**2 + self.config.magErrFloor**2)**0.5
380 good = np.isfinite(srcMagArr) & np.isfinite(refMagArr)
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,
392 def run(self, exposure, sourceCat, expId=0):
393 """!Do photometric calibration - select matches to use and (possibly iteratively) compute
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.).
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
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).
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
423 The measured sources:
424 - Must include PhotoCalConfig.fluxField; the flux measurement to be used for calibration
426 @throws RuntimeError with the following strings:
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.
438 displaySources = display
and lsstDebug.Info(__name__).displaySources
442 from matplotlib
import pyplot
446 self.
figfig = pyplot.figure()
448 filterLabel = exposure.getFilterLabel()
451 matchResults = self.
matchmatch.
run(sourceCat, filterLabel.bandLabel)
452 matches = matchResults.matches
454 reserveResults = self.reserve.
run([mm.second
for mm
in matches], expId=expId)
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:
463 mm.second.set(self.
usedKeyusedKey,
True)
466 sourceKeys = self.
getSourceKeysgetSourceKeys(matches[0].second.schema)
467 arrays = self.
extractMagArraysextractMagArrays(matches, filterLabel, sourceKeys)
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)
474 flux0 = 10**(0.4*r.zp)
475 flux0err = 0.4*math.log(10)*flux0*r.sigma
478 return pipeBase.Struct(
479 photoCalib=photoCalib,
488 """Display sources we'll use for photocal
490 Sources that will be actually used will be green.
491 Sources reserved from the fit will be red.
495 exposure : `lsst.afw.image.ExposureF`
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.
502 Frame number for display.
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)
513 """!Flux calibration code, returning (ZeroPoint, Distribution Width, Number of stars)
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
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.)
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
529 sigmaMax = self.config.sigmaMax
533 indArr = np.argsort(dmag)
536 if srcErr
is not None:
537 dmagErr = srcErr[indArr]
539 dmagErr = np.ones(len(dmag))
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]
547 IQ_TO_STDEV = 0.741301109252802
552 for i
in range(self.config.nIter):
563 hist, edges = np.histogram(dmag, nhist, new=
True)
565 hist, edges = np.histogram(dmag, nhist)
566 imode = np.arange(nhist)[np.where(hist == hist.max())]
568 if imode[-1] - imode[0] + 1 == len(imode):
572 center = 0.5*(edges[imode[0]] + edges[imode[-1] + 1])
574 peak = sum(hist[imode])/len(imode)
578 while j >= 0
and hist[j] > 0.5*peak:
581 q1 = dmag[sum(hist[range(j)])]
584 while j < nhist
and hist[j] > 0.5*peak:
586 j =
min(j, nhist - 1)
587 j =
min(sum(hist[range(j)]), npt - 1)
591 q1 = dmag[int(0.25*npt)]
592 q3 = dmag[int(0.75*npt)]
599 self.log.
debug(
"Photo calibration histogram: center = %.2f, sig = %.2f", center, sig)
603 sigmaMax = dmag[-1] - dmag[0]
605 center = np.median(dmag)
606 q1 = dmag[int(0.25*npt)]
607 q3 = dmag[int(0.75*npt)]
612 if self.config.useMedian:
613 center = np.median(gdmag)
615 gdmagErr = dmagErr[good]
616 center = np.average(gdmag, weights=gdmagErr)
618 q3 = gdmag[
min(int(0.75*npt + 0.5), npt - 1)]
619 q1 = gdmag[
min(int(0.25*npt + 0.5), npt - 1)]
621 sig = IQ_TO_STDEV*(q3 - q1)
623 good =
abs(dmag - center) < self.config.nSigma*
min(sig, sigmaMax)
630 axes = self.
figfig.add_axes((0.1, 0.1, 0.85, 0.80))
632 axes.plot(ref[good], dmag[good] - center,
"b+")
633 axes.errorbar(ref[good], dmag[good] - center, yerr=dmagErr[good],
634 linestyle=
'', color=
'b')
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')
642 axes.plot((-100, 100), (0, 0),
"g-")
644 axes.plot((-100, 100), x*0.05*np.ones(2),
"g--")
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")
655 while i == 0
or reply !=
"c":
657 reply = input(
"Next iteration? [ynhpc] ")
662 print(
"Options: c[ontinue] h[elp] n[o] p[db] y[es]", file=sys.stderr)
665 if reply
in (
"",
"c",
"n",
"p",
"y"):
668 print(
"Unrecognised response: %s" % reply, file=sys.stderr)
675 except Exception
as e:
676 print(
"Error plotting in PhotoCal.getZeroPoint: %s" % e, file=sys.stderr)
683 msg =
"PhotoCal.getZeroPoint: no good stars remain"
686 center = np.average(dmag, weights=dmagErr)
687 msg +=
" on first iteration; using average of all calibration stars"
691 return pipeBase.Struct(
695 elif ngood == old_ngood:
701 dmagErr = dmagErr[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(
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
def getSourceKeys(self, schema)
def getZeroPoint(self, src, ref, srcErr=None, zp0=None)
Flux calibration code, returning (ZeroPoint, Distribution Width, Number of stars)
def extractMagArrays(self, matches, filterLabel, sourceKeys)
Extract magnitude and magnitude error arrays from the given matches.
def __init__(self, refObjLoader, schema=None, **kwds)
Create the photometric calibration task.
def displaySources(self, exposure, matches, reserved, frame=1)
def run(self, exposure, sourceCat, expId=0)
Do photometric calibration - select matches to use and (possibly iteratively) compute the zero point.
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.
std::shared_ptr< PhotoCalib > makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err)
Construct a PhotoCalib from the deprecated Calib-style instFluxMag0/instFluxMag0Err values.
def getRefFluxField(schema, filterName=None)
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.
Angle abs(Angle const &a)