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
computeExposureSummaryStats.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 import warnings
22 import numpy as np
23 from scipy.stats import median_abs_deviation as sigmaMad
24 import astropy.units as units
25 from astropy.time import Time
26 from astropy.coordinates import AltAz, SkyCoord, EarthLocation
27 from lsst.daf.base import DateTime
28 
29 import lsst.pipe.base as pipeBase
30 import lsst.pex.config as pexConfig
31 import lsst.afw.math as afwMath
32 import lsst.afw.image as afwImage
33 import lsst.geom
34 from lsst.utils.timer import timeMethod
35 
36 
37 __all__ = ("ComputeExposureSummaryStatsTask", "ComputeExposureSummaryStatsConfig")
38 
39 
40 class ComputeExposureSummaryStatsConfig(pexConfig.Config):
41  """Config for ComputeExposureSummaryTask"""
42  sigmaClip = pexConfig.Field(
43  dtype=float,
44  doc="Sigma for outlier rejection for sky noise.",
45  default=3.0,
46  )
47  clipIter = pexConfig.Field(
48  dtype=int,
49  doc="Number of iterations of outlier rejection for sky noise.",
50  default=2,
51  )
52  badMaskPlanes = pexConfig.ListField(
53  dtype=str,
54  doc="Mask planes that, if set, the associated pixel should not be included sky noise calculation.",
55  default=("NO_DATA", "SUSPECT"),
56  )
57  starSelection = pexConfig.Field(
58  doc="Field to select sources to be used in the PSF statistics computation.",
59  dtype=str,
60  default="calib_psf_used"
61  )
62  starShape = pexConfig.Field(
63  doc="Base name of columns to use for the source shape in the PSF statistics computation.",
64  dtype=str,
65  default="base_SdssShape"
66  )
67  psfShape = pexConfig.Field(
68  doc="Base name of columns to use for the PSF shape in the PSF statistics computation.",
69  dtype=str,
70  default="base_SdssShape_psf"
71  )
72 
73 
74 class ComputeExposureSummaryStatsTask(pipeBase.Task):
75  """Task to compute exposure summary statistics.
76 
77  This task computes various quantities suitable for DPDD and other
78  downstream processing at the detector centers, including:
79  - psfSigma
80  - psfArea
81  - psfIxx
82  - psfIyy
83  - psfIxy
84  - ra
85  - decl
86  - zenithDistance
87  - zeroPoint
88  - skyBg
89  - skyNoise
90  - meanVar
91  - raCorners
92  - decCorners
93  - astromOffsetMean
94  - astromOffsetStd
95 
96  These additional quantities are computed from the stars in the detector:
97  - psfStarDeltaE1Median
98  - psfStarDeltaE2Median
99  - psfStarDeltaE1Scatter
100  - psfStarDeltaE2Scatter
101  - psfStarDeltaSizeMedian
102  - psfStarDeltaSizeScatter
103  - psfStarScaledDeltaSizeScatter
104  """
105  ConfigClass = ComputeExposureSummaryStatsConfig
106  _DefaultName = "computeExposureSummaryStats"
107 
108  @timeMethod
109  def run(self, exposure, sources, background):
110  """Measure exposure statistics from the exposure, sources, and background.
111 
112  Parameters
113  ----------
114  exposure : `lsst.afw.image.ExposureF`
115  sources : `lsst.afw.table.SourceCatalog`
116  background : `lsst.afw.math.BackgroundList`
117 
118  Returns
119  -------
120  summary : `lsst.afw.image.ExposureSummary`
121  """
122  self.log.info("Measuring exposure statistics")
123 
124  bbox = exposure.getBBox()
125 
126  psf = exposure.getPsf()
127  if psf is not None:
128  shape = psf.computeShape(bbox.getCenter())
129  psfSigma = shape.getDeterminantRadius()
130  psfIxx = shape.getIxx()
131  psfIyy = shape.getIyy()
132  psfIxy = shape.getIxy()
133  im = psf.computeKernelImage(bbox.getCenter())
134  # The calculation of effective psf area is taken from
135  # meas_base/src/PsfFlux.cc#L112. See
136  # https://github.com/lsst/meas_base/blob/
137  # 750bffe6620e565bda731add1509507f5c40c8bb/src/PsfFlux.cc#L112
138  psfArea = np.sum(im.array)/np.sum(im.array**2.)
139  else:
140  psfSigma = np.nan
141  psfIxx = np.nan
142  psfIyy = np.nan
143  psfIxy = np.nan
144  psfArea = np.nan
145 
146  wcs = exposure.getWcs()
147  if wcs is not None:
148  sph_pts = wcs.pixelToSky(lsst.geom.Box2D(bbox).getCorners())
149  raCorners = [float(sph.getRa().asDegrees()) for sph in sph_pts]
150  decCorners = [float(sph.getDec().asDegrees()) for sph in sph_pts]
151 
152  sph_pt = wcs.pixelToSky(bbox.getCenter())
153  ra = sph_pt.getRa().asDegrees()
154  decl = sph_pt.getDec().asDegrees()
155  else:
156  raCorners = [float(np.nan)]*4
157  decCorners = [float(np.nan)]*4
158  ra = np.nan
159  decl = np.nan
160 
161  photoCalib = exposure.getPhotoCalib()
162  if photoCalib is not None:
163  zeroPoint = 2.5*np.log10(photoCalib.getInstFluxAtZeroMagnitude())
164  else:
165  zeroPoint = np.nan
166 
167  visitInfo = exposure.getInfo().getVisitInfo()
168  date = visitInfo.getDate()
169 
170  if date.isValid():
171  # We compute the zenithDistance at the center of the detector rather
172  # than use the boresight value available via the visitInfo, because
173  # the zenithDistance may vary significantly over a large field of view.
174  observatory = visitInfo.getObservatory()
175  loc = EarthLocation(lat=observatory.getLatitude().asDegrees()*units.deg,
176  lon=observatory.getLongitude().asDegrees()*units.deg,
177  height=observatory.getElevation()*units.m)
178  obstime = Time(visitInfo.getDate().get(system=DateTime.MJD),
179  location=loc, format='mjd')
180  coord = SkyCoord(ra*units.degree, decl*units.degree, obstime=obstime, location=loc)
181  with warnings.catch_warnings():
182  warnings.simplefilter('ignore')
183  altaz = coord.transform_to(AltAz)
184 
185  zenithDistance = 90.0 - altaz.alt.degree
186  else:
187  zenithDistance = np.nan
188 
189  if background is not None:
190  bgStats = (bg[0].getStatsImage().getImage().array
191  for bg in background)
192  skyBg = sum(np.median(bg[np.isfinite(bg)]) for bg in bgStats)
193  else:
194  skyBg = np.nan
195 
196  statsCtrl = afwMath.StatisticsControl()
197  statsCtrl.setNumSigmaClip(self.config.sigmaClip)
198  statsCtrl.setNumIter(self.config.clipIter)
199  statsCtrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes))
200  statsCtrl.setNanSafe(True)
201 
202  statObj = afwMath.makeStatistics(exposure.getMaskedImage(), afwMath.STDEVCLIP,
203  statsCtrl)
204  skyNoise, _ = statObj.getResult(afwMath.STDEVCLIP)
205 
206  statObj = afwMath.makeStatistics(exposure.getMaskedImage().getVariance(),
207  exposure.getMaskedImage().getMask(),
208  afwMath.MEANCLIP, statsCtrl)
209  meanVar, _ = statObj.getResult(afwMath.MEANCLIP)
210 
211  md = exposure.getMetadata()
212  if 'SFM_ASTROM_OFFSET_MEAN' in md:
213  astromOffsetMean = md['SFM_ASTROM_OFFSET_MEAN']
214  astromOffsetStd = md['SFM_ASTROM_OFFSET_STD']
215  else:
216  astromOffsetMean = np.nan
217  astromOffsetStd = np.nan
218 
219  psfStats = self._computePsfStats_computePsfStats(sources)
220 
221  # Note that all numpy values have to be explicitly converted to
222  # python floats for yaml serialization.
224  psfSigma=float(psfSigma),
225  psfArea=float(psfArea),
226  psfIxx=float(psfIxx),
227  psfIyy=float(psfIyy),
228  psfIxy=float(psfIxy),
229  ra=float(ra),
230  decl=float(decl),
231  zenithDistance=float(zenithDistance),
232  zeroPoint=float(zeroPoint),
233  skyBg=float(skyBg),
234  skyNoise=float(skyNoise),
235  meanVar=float(meanVar),
236  raCorners=raCorners,
237  decCorners=decCorners,
238  astromOffsetMean=astromOffsetMean,
239  astromOffsetStd=astromOffsetStd,
240  nPsfStar=psfStats.nPsfStar,
241  psfStarDeltaE1Median=psfStats.psfStarDeltaE1Median,
242  psfStarDeltaE2Median=psfStats.psfStarDeltaE2Median,
243  psfStarDeltaE1Scatter=psfStats.psfStarDeltaE1Scatter,
244  psfStarDeltaE2Scatter=psfStats.psfStarDeltaE2Scatter,
245  psfStarDeltaSizeMedian=psfStats.psfStarDeltaSizeMedian,
246  psfStarDeltaSizeScatter=psfStats.psfStarDeltaSizeScatter,
247  psfStarScaledDeltaSizeScatter=psfStats.psfStarScaledDeltaSizeScatter
248  )
249 
250  return summary
251 
252  def _computePsfStats(self, sources):
253  """Compute psf residual statistics.
254 
255  All residuals are computed using median statistics on the difference
256  between the sources and the models.
257 
258  Parameters
259  ----------
260  sources : `lsst.afw.table.SourceCatalog`
261  Source catalog on which to measure the PSF statistics.
262 
263  Returns
264  -------
265  psfStats : `lsst.pipe.base.Struct`
266  Struct with various psf stats.
267  """
268  psfStats = pipeBase.Struct(nPsfStar=0,
269  psfStarDeltaE1Median=float(np.nan),
270  psfStarDeltaE2Median=float(np.nan),
271  psfStarDeltaE1Scatter=float(np.nan),
272  psfStarDeltaE2Scatter=float(np.nan),
273  psfStarDeltaSizeMedian=float(np.nan),
274  psfStarDeltaSizeScatter=float(np.nan),
275  psfStarScaledDeltaSizeScatter=float(np.nan))
276 
277  if sources is None:
278  # No sources are available (as in some tests)
279  return psfStats
280 
281  names = sources.schema.getNames()
282  if self.config.starSelection not in names or self.config.starShape + '_flag' not in names:
283  # The source catalog does not have the necessary fields (as in some tests)
284  return psfStats
285 
286  mask = sources[self.config.starSelection] & (~sources[self.config.starShape + '_flag'])
287  nPsfStar = mask.sum()
288 
289  if nPsfStar == 0:
290  # No stars to measure statistics, so we must return the defaults
291  # of 0 stars and NaN values.
292  return psfStats
293 
294  starXX = sources[self.config.starShape + '_xx'][mask]
295  starYY = sources[self.config.starShape + '_yy'][mask]
296  starXY = sources[self.config.starShape + '_xy'][mask]
297  psfXX = sources[self.config.psfShape + '_xx'][mask]
298  psfYY = sources[self.config.psfShape + '_yy'][mask]
299  psfXY = sources[self.config.psfShape + '_xy'][mask]
300 
301  starSize = (starXX*starYY - starXY**2.)**0.25
302  starE1 = (starXX - starYY)/(starXX + starYY)
303  starE2 = 2*starXY/(starXX + starYY)
304  starSizeMedian = np.median(starSize)
305 
306  psfSize = (psfXX*psfYY - psfXY**2)**0.25
307  psfE1 = (psfXX - psfYY)/(psfXX + psfYY)
308  psfE2 = 2*psfXY/(psfXX + psfYY)
309 
310  psfStarDeltaE1Median = np.median(starE1 - psfE1)
311  psfStarDeltaE1Scatter = sigmaMad(starE1 - psfE1, scale='normal')
312  psfStarDeltaE2Median = np.median(starE2 - psfE2)
313  psfStarDeltaE2Scatter = sigmaMad(starE2 - psfE2, scale='normal')
314 
315  psfStarDeltaSizeMedian = np.median(starSize - psfSize)
316  psfStarDeltaSizeScatter = sigmaMad(starSize - psfSize, scale='normal')
317  psfStarScaledDeltaSizeScatter = psfStarDeltaSizeScatter/starSizeMedian**2.
318 
319  psfStats.nPsfStar = int(nPsfStar)
320  psfStats.psfStarDeltaE1Median = float(psfStarDeltaE1Median)
321  psfStats.psfStarDeltaE2Median = float(psfStarDeltaE2Median)
322  psfStats.psfStarDeltaE1Scatter = float(psfStarDeltaE1Scatter)
323  psfStats.psfStarDeltaE2Scatter = float(psfStarDeltaE2Scatter)
324  psfStats.psfStarDeltaSizeMedian = float(psfStarDeltaSizeMedian)
325  psfStats.psfStarDeltaSizeScatter = float(psfStarDeltaSizeScatter)
326  psfStats.psfStarScaledDeltaSizeScatter = float(psfStarScaledDeltaSizeScatter)
327 
328  return psfStats
Pass parameters to a Statistics object.
Definition: Statistics.h:92
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition: Statistics.h:359