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
34 from lsst.utils.timer
import timeMethod
37 __all__ = (
"ComputeExposureSummaryStatsTask",
"ComputeExposureSummaryStatsConfig")
41 """Config for ComputeExposureSummaryTask"""
42 sigmaClip = pexConfig.Field(
44 doc=
"Sigma for outlier rejection for sky noise.",
47 clipIter = pexConfig.Field(
49 doc=
"Number of iterations of outlier rejection for sky noise.",
52 badMaskPlanes = pexConfig.ListField(
54 doc=
"Mask planes that, if set, the associated pixel should not be included sky noise calculation.",
55 default=(
"NO_DATA",
"SUSPECT"),
57 starSelection = pexConfig.Field(
58 doc=
"Field to select sources to be used in the PSF statistics computation.",
60 default=
"calib_psf_used"
62 starShape = pexConfig.Field(
63 doc=
"Base name of columns to use for the source shape in the PSF statistics computation.",
65 default=
"base_SdssShape"
67 psfShape = pexConfig.Field(
68 doc=
"Base name of columns to use for the PSF shape in the PSF statistics computation.",
70 default=
"base_SdssShape_psf"
75 """Task to compute exposure summary statistics.
77 This task computes various quantities suitable for DPDD and other
78 downstream processing at the detector centers, including:
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
105 ConfigClass = ComputeExposureSummaryStatsConfig
106 _DefaultName =
"computeExposureSummaryStats"
109 def run(self, exposure, sources, background):
110 """Measure exposure statistics from the exposure, sources, and background.
114 exposure : `lsst.afw.image.ExposureF`
115 sources : `lsst.afw.table.SourceCatalog`
116 background : `lsst.afw.math.BackgroundList`
120 summary : `lsst.afw.image.ExposureSummary`
122 self.log.
info(
"Measuring exposure statistics")
124 bbox = exposure.getBBox()
126 psf = exposure.getPsf()
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())
138 psfArea = np.sum(im.array)/np.sum(im.array**2.)
146 wcs = exposure.getWcs()
149 raCorners = [float(sph.getRa().asDegrees())
for sph
in sph_pts]
150 decCorners = [float(sph.getDec().asDegrees())
for sph
in sph_pts]
152 sph_pt = wcs.pixelToSky(bbox.getCenter())
153 ra = sph_pt.getRa().asDegrees()
154 decl = sph_pt.getDec().asDegrees()
156 raCorners = [float(np.nan)]*4
157 decCorners = [float(np.nan)]*4
161 photoCalib = exposure.getPhotoCalib()
162 if photoCalib
is not None:
163 zeroPoint = 2.5*np.log10(photoCalib.getInstFluxAtZeroMagnitude())
167 visitInfo = exposure.getInfo().getVisitInfo()
168 date = visitInfo.getDate()
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)
185 zenithDistance = 90.0 - altaz.alt.degree
187 zenithDistance = np.nan
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)
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)
204 skyNoise, _ = statObj.getResult(afwMath.STDEVCLIP)
207 exposure.getMaskedImage().getMask(),
208 afwMath.MEANCLIP, statsCtrl)
209 meanVar, _ = statObj.getResult(afwMath.MEANCLIP)
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']
216 astromOffsetMean = np.nan
217 astromOffsetStd = np.nan
224 psfSigma=float(psfSigma),
225 psfArea=float(psfArea),
226 psfIxx=float(psfIxx),
227 psfIyy=float(psfIyy),
228 psfIxy=float(psfIxy),
231 zenithDistance=float(zenithDistance),
232 zeroPoint=float(zeroPoint),
234 skyNoise=float(skyNoise),
235 meanVar=float(meanVar),
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
252 def _computePsfStats(self, sources):
253 """Compute psf residual statistics.
255 All residuals are computed using median statistics on the difference
256 between the sources and the models.
260 sources : `lsst.afw.table.SourceCatalog`
261 Source catalog on which to measure the PSF statistics.
265 psfStats : `lsst.pipe.base.Struct`
266 Struct with various psf stats.
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))
281 names = sources.schema.getNames()
282 if self.config.starSelection
not in names
or self.config.starShape +
'_flag' not in names:
286 mask = sources[self.config.starSelection] & (~sources[self.config.starShape +
'_flag'])
287 nPsfStar = mask.sum()
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]
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)
306 psfSize = (psfXX*psfYY - psfXY**2)**0.25
307 psfE1 = (psfXX - psfYY)/(psfXX + psfYY)
308 psfE2 = 2*psfXY/(psfXX + psfYY)
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')
315 psfStarDeltaSizeMedian = np.median(starSize - psfSize)
316 psfStarDeltaSizeScatter =
sigmaMad(starSize - psfSize, scale=
'normal')
317 psfStarScaledDeltaSizeScatter = psfStarDeltaSizeScatter/starSizeMedian**2.
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)
Pass parameters to a Statistics object.
A floating-point coordinate rectangle geometry.
def run(self, exposure, sources, background)
def _computePsfStats(self, sources)
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)