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
measureApCorr.py
Go to the documentation of this file.
1 from __future__ import absolute_import, division
2 #
3 # LSST Data Management System
4 # Copyright 2008-2015 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 
24 import numpy
25 
26 import lsst.pex.config
27 from lsst.afw.image import ApCorrMap
28 from lsst.afw.math import ChebyshevBoundedField, ChebyshevBoundedFieldConfig
29 from lsst.pipe.base import Task, Struct
30 from .apCorrRegistry import getApCorrNameSet
31 
32 __all__ = ("MeasureApCorrConfig", "MeasureApCorrTask")
33 
34 class FluxKeys(object):
35  """A collection of keys for a given flux measurement algorithm
36  """
37  __slots__ = ("flux", "err", "flag", "used") # prevent accidentally adding fields
38 
39  def __init__(self, name, schema):
40  """Construct a FluxKeys
41 
42  @parma[in] name name of flux measurement algorithm, e.g. "base_PsfFlux"
43  @param[in,out] schema catalog schema containing the flux field
44  read: {name}_flux, {name}_fluxSigma, {name}_flag
45  added: apcorr_{name}_used
46  """
47  self.flux = schema.find(name + "_flux").key
48  self.err = schema.find(name + "_fluxSigma").key
49  self.flag = schema.find(name + "_flag").key
50  self.used = schema.addField("apcorr_" + name + "_used", type="Flag",
51  doc="set if source was used in measuring aperture correction")
52 
53 # The following block adds links to these tasks from the Task Documentation page.
54 ## \addtogroup LSST_task_documentation
55 ## \{
56 ## \page measBase_measureApCorrTask
57 ## \ref MeasureApCorrTask "MeasureApCorrTask"
58 ## Task to measure aperture correction
59 ## \}
60 
61 class MeasureApCorrConfig(lsst.pex.config.Config):
62  """!Configuration for MeasureApCorrTask
63  """
64  refFluxName = lsst.pex.config.Field(
65  doc = "Field name prefix for the flux other measurements should be aperture corrected to match",
66  dtype = str,
67  default = "base_CircularApertureFlux_17_0",
68  )
69  inputFilterFlag = lsst.pex.config.Field(
70  doc = "Name of a flag field that indicates that a source should be used to constrain the" +
71  " aperture corrections",
72  dtype = str,
73  default = "calib_psfUsed",
74  )
75  minDegreesOfFreedom = lsst.pex.config.RangeField(
76  doc = "Minimum number of degrees of freedom (# of valid data points - # of parameters);" +
77  " if this is exceeded, the order of the fit is decreased (in both dimensions), and" +
78  " if we can't decrease it enough, we'll raise ValueError.",
79  dtype = int,
80  default = 1,
81  min = 1,
82  )
83  fitConfig = lsst.pex.config.ConfigField(
84  doc = "Configuration used in fitting the aperture correction fields",
85  dtype = ChebyshevBoundedFieldConfig,
86  )
87  numIter = lsst.pex.config.Field(
88  doc = "Number of iterations for sigma clipping",
89  dtype = int,
90  default = 4,
91  )
92  numSigmaClip = lsst.pex.config.Field(
93  doc = "Number of standard devisations to clip at",
94  dtype = float,
95  default = 3.0,
96  )
97 
98 class MeasureApCorrTask(Task):
99  """!Task to measure aperture correction
100 
101  \section measBase_MeasureApCorrTask_Contents Contents
102 
103  - \ref measBase_MeasureApCorrTask_Purpose
104  - \ref measBase_MeasureApCorrTask_Config
105  - \ref measBase_MeasureApCorrTask_Debug
106 
107  \section measBase_MeasureApCorrTask_Purpose Description
108 
109  \copybrief MeasureApCorrTask
110 
111  This task measures aperture correction for the flux fields returned by
112  lsst.meas.base.getApCorrNameSet()
113 
114  The main method is \ref MeasureApCorrTask.run "run".
115 
116  \section measBase_MeasureApCorrTask_Config Configuration parameters
117 
118  See \ref MeasureApCorrConfig
119 
120  \section measBase_MeasureApCorrTask_Debug Debug variables
121 
122  This task has no debug variables.
123  """
124  ConfigClass = MeasureApCorrConfig
125  _DefaultName = "measureApCorr"
126 
127  def __init__(self, schema, **kwds):
128  """!Construct a MeasureApCorrTask
129 
130  For every name in lsst.meas.base.getApCorrNameSet():
131  - If the corresponding flux fields exist in the schema:
132  - Add a new field apcorr_{name}_used
133  - Add an entry to the self.toCorrect dict
134  - Otherwise silently skip the name
135  """
136  Task.__init__(self, **kwds)
137  self.refFluxKeys = FluxKeys(self.config.refFluxName, schema)
138  self.toCorrect = {} # dict of flux field name prefix: FluxKeys instance
139  for name in getApCorrNameSet():
140  try:
141  self.toCorrect[name] = FluxKeys(name, schema)
142  except KeyError:
143  # if a field in the registry is missing, just ignore it.
144  pass
145  self.inputFilterFlag = schema.find(self.config.inputFilterFlag).key
146 
147  def run(self, bbox, catalog):
148  """!Measure aperture correction
149 
150  @return an lsst.pipe.base.Struct containing:
151  - apCorrMap: an aperture correction map (lsst.afw.image.ApCorrMap) that contains two entries
152  for each flux field:
153  - flux field (e.g. base_PsfFlux_flux): 2d model
154  - flux sigma field (e.g. base_PsfFlux_fluxSigma): 2d model of error
155  """
156  self.log.info("Measuring aperture corrections for %d flux fields" % (len(self.toCorrect),))
157  # First, create a subset of the catalog that contains only objects with inputFilterFlag set
158  # and non-flagged reference fluxes.
159  subset1 = [record for record in catalog
160  if record.get(self.inputFilterFlag) and not record.get(self.refFluxKeys.flag)]
161 
162  apCorrMap = ApCorrMap()
163 
164  # Outer loop over the fields we want to correct
165  for name, keys in self.toCorrect.iteritems():
166  fluxName = name + "_flux"
167  fluxSigmaName = name + "_fluxSigma"
168 
169  # Create a more restricted subset with only the objects where the to-be-correct flux
170  # is not flagged.
171  subset2 = [record for record in subset1 if not record.get(keys.flag)]
172 
173  # Check that we have enough data points that we have at least the minimum of degrees of
174  # freedom specified in the config.
175  if len(subset2) - 1 < self.config.minDegreesOfFreedom:
176  self.log.warn("Only %d sources for calculation of aperture correction for '%s'; "
177  "setting to 1.0" % (len(subset2), name,))
178  apCorrMap[fluxName] = ChebyshevBoundedField(bbox, numpy.ones((1,1), dtype=float))
179  apCorrMap[fluxSigmaName] = ChebyshevBoundedField(bbox, numpy.zeros((1,1), dtype=float))
180  continue
181 
182  # If we don't have enough data points to constrain the fit, reduce the order until we do
183  ctrl = self.config.fitConfig.makeControl()
184  while len(subset2) - ctrl.computeSize() < self.config.minDegreesOfFreedom:
185  if ctrl.orderX > 0:
186  ctrl.orderX -= 1
187  if ctrl.orderY > 0:
188  ctrl.orderY -= 1
189 
190  # Fill numpy arrays with positions and the ratio of the reference flux to the to-correct flux
191  x = numpy.zeros(len(subset2), dtype=float)
192  y = numpy.zeros(len(subset2), dtype=float)
193  apCorrData = numpy.zeros(len(subset2), dtype=float)
194  indices = numpy.arange(len(subset2), dtype=int)
195  for n, record in enumerate(subset2):
196  x[n] = record.getX()
197  y[n] = record.getY()
198  apCorrData[n] = record.get(self.refFluxKeys.flux)/record.get(keys.flux)
199 
200  for _i in range(self.config.numIter):
201 
202  # Do the fit, save it in the output map
203  apCorrField = ChebyshevBoundedField.fit(bbox, x, y, apCorrData, ctrl)
204 
205  # Compute errors empirically, using the RMS difference between the true reference flux and the
206  # corrected to-be-corrected flux.
207  apCorrDiffs = apCorrField.evaluate(x, y)
208  apCorrDiffs -= apCorrData
209  apCorrErr = numpy.mean(apCorrDiffs**2)**0.5
210 
211  # Clip bad data points
212  apCorrDiffLim = self.config.numSigmaClip * apCorrErr
213  keep = numpy.fabs(apCorrDiffs) < apCorrDiffLim
214  x = x[keep]
215  y = y[keep]
216  apCorrData = apCorrData[keep]
217  indices = indices[keep]
218 
219  # Final fit after clipping
220  apCorrField = ChebyshevBoundedField.fit(bbox, x, y, apCorrData, ctrl)
221 
222  self.log.info("Aperture correction for %s: RMS %f from %d" %
223  (name, numpy.mean((apCorrField.evaluate(x, y) - apCorrData)**2)**0.5, len(indices)))
224 
225  # Save the result in the output map
226  # The error is constant spatially (we could imagine being
227  # more clever, but we're not yet sure if it's worth the effort).
228  # We save the errors as a 0th-order ChebyshevBoundedField
229  apCorrMap[fluxName] = apCorrField
230  apCorrErrCoefficients = numpy.array([[apCorrErr]], dtype=float)
231  apCorrMap[fluxSigmaName] = ChebyshevBoundedField(bbox, apCorrErrCoefficients)
232 
233  # Record which sources were used
234  for i in indices:
235  subset2[i].set(keys.used, True)
236 
237  return Struct(
238  apCorrMap = apCorrMap,
239  )
Task to measure aperture correction.
A thin wrapper around std::map to allow aperture corrections to be attached to Exposures.
Definition: ApCorrMap.h:42
def run
Measure aperture correction.
def getApCorrNameSet
Return a copy of the set of field name prefixes for fluxes that should be aperture corrected...
def __init__
Construct a MeasureApCorrTask.
A BoundedField based on 2-d Chebyshev polynomials of the first kind.
Configuration for MeasureApCorrTask.