Coverage for tests / test_computeExposureSummaryStats.py: 16%

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 

22"""Test ComputeExposureSummaryStatsTask. 

23""" 

24import unittest 

25 

26import numpy as np 

27 

28import lsst.utils.tests 

29from lsst.afw.detection import GaussianPsf 

30import lsst.afw.image as afwImage 

31import lsst.afw.math as afwMath 

32from lsst.daf.base import DateTime, PropertyList 

33from lsst.afw.coord import Observatory 

34from lsst.afw.geom import makeCdMatrix, makeSkyWcs 

35from lsst.afw.cameraGeom.testUtils import DetectorWrapper 

36from lsst.pipe.tasks.computeExposureSummaryStats import ComputeExposureSummaryStatsTask 

37from lsst.pipe.tasks.computeExposureSummaryStats import ComputeExposureSummaryStatsConfig 

38from lsst.pipe.tasks.computeExposureSummaryStats import compute_magnitude_limit 

39from lsst.pipe.tasks.computeExposureSummaryStats import compute_effective_time 

40 

41# Values from syseng_throughputs notebook assuming 30s exposure 

42# consisting of 2x15s snaps each with readnoise of 9e- 

43fwhm_eff_fid = {'g': 0.87, 'r': 0.83, 'i': 0.80} 

44skycounts_fid = {'g': 463.634122, 'r': 988.626863, 'i': 1588.280513} 

45zeropoint_fid = {'g': 28.508375, 'r': 28.360838, 'i': 28.171396} 

46readnoise_fid = {'g': 12.73, 'r': 12.73, 'i': 12.73} 

47# Output magnitude limit from syseng_throughputs notebook 

48m5_fid = {'g': 24.90, 'r': 24.48, 'i': 24.10} 

49 

50 

51class ComputeExposureSummaryTestCase(lsst.utils.tests.TestCase): 

52 

53 def testComputeExposureSummary(self): 

54 """Make a fake exposure and background and compute summary. 

55 """ 

56 np.random.seed(12345) 

57 

58 # Make an exposure with a noise image 

59 exposure = afwImage.ExposureF(100, 100) 

60 

61 band = "i" 

62 physical_filter = "test-i" 

63 exposure.setFilter(afwImage.FilterLabel(band=band, physical=physical_filter)) 

64 

65 readNoise = 5.0 

66 gain = 1.0 

67 detector = DetectorWrapper(numAmps=1).detector 

68 metadata = PropertyList() 

69 metadata.add("LSST ISR UNITS", "electron") 

70 for amp in detector.getAmplifiers(): 

71 metadata.add(f"LSST ISR READNOISE {amp.getName()}", readNoise) 

72 metadata.add(f"LSST ISR GAIN {amp.getName()}", gain) 

73 exposure.setDetector(detector) 

74 exposure.setMetadata(metadata) 

75 

76 skyMean = 100.0 

77 skySigma = 10.0 

78 exposure.getImage().getArray()[:, :] = np.random.normal(skyMean, skySigma, size=(100, 100)) 

79 exposure.getVariance().getArray()[:, :] = skySigma**2. 

80 

81 # Set the visitInfo 

82 date = DateTime(date=59234.7083333334, system=DateTime.DateSystem.MJD) 

83 observatory = Observatory(-70.7366*lsst.geom.degrees, -30.2407*lsst.geom.degrees, 

84 2650.) 

85 expTime = 10.0 

86 visitInfo = afwImage.VisitInfo(exposureTime=expTime, 

87 date=date, 

88 observatory=observatory) 

89 exposure.getInfo().setVisitInfo(visitInfo) 

90 

91 # Install a Gaussian PSF 

92 psfSize = 2.0 

93 psf = GaussianPsf(5, 5, psfSize) 

94 exposure.setPsf(psf) 

95 

96 # Install a simple WCS 

97 scale = 0.2*lsst.geom.arcseconds 

98 raCenter = 300.0*lsst.geom.degrees 

99 decCenter = 0.0*lsst.geom.degrees 

100 cdMatrix = makeCdMatrix(scale=scale) 

101 skyWcs = makeSkyWcs(crpix=exposure.getBBox().getCenter(), 

102 crval=lsst.geom.SpherePoint(raCenter, decCenter), 

103 cdMatrix=cdMatrix) 

104 pixelScale = skyWcs.getPixelScale().asArcseconds() 

105 exposure.setWcs(skyWcs) 

106 

107 # Install a simple photoCalib 

108 photoCalib = afwImage.PhotoCalib(calibrationMean=0.3) 

109 zp = 2.5*np.log10(photoCalib.getInstFluxAtZeroMagnitude()) 

110 exposure.setPhotoCalib(photoCalib) 

111 

112 # Install a simple apCorrMap 

113 apCorrMap = afwImage.ApCorrMap() 

114 apCorrMap.set( 

115 "base_PsfFlux_instFlux", afwMath.ChebyshevBoundedField(exposure.getBBox(), np.zeros((3, 3))) 

116 ) 

117 exposure.setApCorrMap(apCorrMap) 

118 

119 # Compute the background image 

120 bgGridSize = 10 

121 bctrl = afwMath.BackgroundControl(afwMath.Interpolate.NATURAL_SPLINE) 

122 bctrl.setNxSample(int(exposure.getMaskedImage().getWidth()/bgGridSize) + 1) 

123 bctrl.setNySample(int(exposure.getMaskedImage().getHeight()/bgGridSize) + 1) 

124 backobj = afwMath.makeBackground(exposure.getMaskedImage().getImage(), bctrl) 

125 background = afwMath.BackgroundList() 

126 background.append(backobj) 

127 

128 # Configure and run the task 

129 expSummaryTask = ComputeExposureSummaryStatsTask() 

130 config = expSummaryTask.config 

131 

132 # Configure nominal values for effective time calculation (normalized to 1s exposure) 

133 config.fiducialZeroPoint = {band: float(zp - 2.5*np.log10(expTime))} 

134 config.fiducialPsfSigma = {band: float(psfSize)} 

135 config.fiducialSkyBackground = {band: float(skyMean/expTime)} 

136 config.fiducialReadNoise = {band: float(readNoise)} 

137 config.fiducialExpTime = {band: expTime} 

138 config.fiducialMagLim = {band: 26.584} 

139 

140 # Run the task 

141 summary = expSummaryTask.run(exposure, None, background) 

142 

143 # Test the outputs 

144 self.assertFloatsAlmostEqual(summary.psfSigma, psfSize) 

145 self.assertFloatsAlmostEqual(summary.psfIxx, psfSize**2.) 

146 self.assertFloatsAlmostEqual(summary.psfIyy, psfSize**2.) 

147 self.assertFloatsAlmostEqual(summary.psfIxy, 0.0) 

148 self.assertFloatsAlmostEqual(summary.psfArea, 23.088975164455444) 

149 

150 self.assertFloatsAlmostEqual(summary.psfTraceRadiusDelta, 0.0) 

151 self.assertFloatsAlmostEqual(summary.psfApFluxDelta, 0.0) 

152 self.assertFloatsAlmostEqual(summary.psfApCorrSigmaScaledDelta, 0.0) 

153 nan = float("nan") 

154 self.assertFloatsAlmostEqual(summary.psfStarDeltaE1Median, nan, ignoreNaNs=True) 

155 self.assertFloatsAlmostEqual(summary.psfStarDeltaE2Median, nan, ignoreNaNs=True) 

156 self.assertFloatsAlmostEqual(summary.psfStarDeltaE1Scatter, nan, ignoreNaNs=True) 

157 self.assertFloatsAlmostEqual(summary.psfStarDeltaE2Scatter, nan, ignoreNaNs=True) 

158 self.assertFloatsAlmostEqual(summary.psfStarDeltaSizeMedian, nan, ignoreNaNs=True) 

159 self.assertFloatsAlmostEqual(summary.psfStarDeltaSizeScatter, nan, ignoreNaNs=True) 

160 self.assertFloatsAlmostEqual(summary.psfStarScaledDeltaSizeScatter, nan, ignoreNaNs=True) 

161 self.assertFloatsAlmostEqual(summary.maxDistToNearestPsf, nan, ignoreNaNs=True) 

162 self.assertFloatsAlmostEqual(summary.starEMedian, nan, ignoreNaNs=True) 

163 self.assertFloatsAlmostEqual(summary.starUnNormalizedEMedian, nan, ignoreNaNs=True) 

164 

165 delta = (scale*50).asDegrees() 

166 for a, b in zip(summary.raCorners, 

167 [raCenter.asDegrees() + delta, raCenter.asDegrees() - delta, 

168 raCenter.asDegrees() - delta, raCenter.asDegrees() + delta]): 

169 self.assertFloatsAlmostEqual(a, b, atol=1e-10) 

170 for a, b in zip(summary.decCorners, 

171 [decCenter.asDegrees() - delta, decCenter.asDegrees() - delta, 

172 decCenter.asDegrees() + delta, decCenter.asDegrees() + delta]): 

173 self.assertFloatsAlmostEqual(a, b, atol=1e-10) 

174 

175 self.assertFloatsAlmostEqual(summary.ra, raCenter.asDegrees(), atol=1e-10) 

176 self.assertFloatsAlmostEqual(summary.dec, decCenter.asDegrees(), atol=1e-10) 

177 

178 self.assertFloatsAlmostEqual(summary.pixelScale, pixelScale) 

179 

180 self.assertFloatsAlmostEqual(summary.zeroPoint, zp) 

181 self.assertFloatsAlmostEqual(summary.expTime, expTime) 

182 

183 # Need to compare background level and noise 

184 # These are only approximately 0+/-10 because of the small image 

185 self.assertFloatsAlmostEqual(summary.skyBg, skyMean, rtol=1e-3) 

186 self.assertFloatsAlmostEqual(summary.meanVar, skySigma**2.) 

187 

188 self.assertFloatsAlmostEqual(summary.zenithDistance, 30.57112, atol=1e-5) 

189 

190 # Effective exposure time and depth 

191 self.assertFloatsAlmostEqual(summary.magLim, config.fiducialMagLim[band], rtol=1e-3) 

192 self.assertFloatsAlmostEqual(summary.effTime, expTime, rtol=1e-3) 

193 

194 def testComputeMagnitudeLimit(self): 

195 """Test the magnitude limit calculation using fiducials from SMTN-002 

196 and syseng_throughputs.""" 

197 

198 # Assumed values from SMTN-002 

199 snr = 5 

200 gain = 1.0 

201 # For testing non-unity gain 

202 gain_test = 2.0 

203 exposure_time = 30 

204 pixel_scale = 0.2 

205 

206 for band in ['g', 'r', 'i']: 

207 # Translate into DM quantities 

208 psfArea = 2.266 * (fwhm_eff_fid[band] / pixel_scale)**2 

209 skyBg = skycounts_fid[band] 

210 zeroPoint = zeropoint_fid[band] + 2.5*np.log10(exposure_time) 

211 readNoise = readnoise_fid[band] 

212 

213 # Calculate the M5 values 

214 m5 = compute_magnitude_limit(psfArea, skyBg, zeroPoint, readNoise, gain, snr) 

215 self.assertFloatsAlmostEqual(m5, m5_fid[band], atol=1e-2) 

216 

217 # Changing the gain consistently should not change M5 

218 m5_gain = compute_magnitude_limit(psfArea, skyBg/gain_test, 

219 zeroPoint - 2.5*np.log10(gain_test), 

220 readNoise/gain_test, gain_test, snr) 

221 self.assertFloatsAlmostEqual(m5_gain, m5, atol=1e-6) 

222 

223 # Calculate effective time values 

224 eff_time = compute_effective_time(m5, m5_fid[band], exposure_time) 

225 self.assertFloatsAlmostEqual(eff_time, exposure_time, rtol=1e-2) 

226 eff_time_half = compute_effective_time(m5 - 0.375, m5_fid[band], exposure_time) 

227 self.assertFloatsAlmostEqual(eff_time_half, exposure_time/2.0, rtol=1e-2) 

228 

229 # Check that input NaN lead to output NaN 

230 nan = float('nan') 

231 m5 = compute_magnitude_limit(nan, skyBg, zeroPoint, readNoise, gain, snr) 

232 self.assertFloatsAlmostEqual(m5, nan, ignoreNaNs=True) 

233 m5 = compute_magnitude_limit(psfArea, nan, zeroPoint, readNoise, gain, snr) 

234 self.assertFloatsAlmostEqual(m5, nan, ignoreNaNs=True) 

235 m5 = compute_magnitude_limit(psfArea, skyBg, nan, readNoise, gain, snr) 

236 self.assertFloatsAlmostEqual(m5, nan, ignoreNaNs=True) 

237 m5 = compute_magnitude_limit(psfArea, skyBg, zeroPoint, nan, gain, snr) 

238 self.assertFloatsAlmostEqual(m5, nan, ignoreNaNs=True) 

239 

240 def testComputeFiducialMagnitudeLimit(self): 

241 """Test the magnitude limit calculation using fiducials from SMTN-002 

242 and syseng_throughputs.""" 

243 

244 # Setup the config 

245 config = ComputeExposureSummaryStatsConfig() 

246 # Configure fiducial values. 

247 # These come from obs_lsst/config/lsstCam 

248 config.fiducialZeroPoint = { 

249 "u": 26.52, "g": 28.51, "r": 28.36, "i": 28.17, "z": 27.78, "y": 26.82, 

250 } 

251 config.fiducialPsfSigma = { 

252 "u": 1.72, "g": 1.63, "r": 1.56, "i": 1.51, "z": 1.47, "y": 1.44, 

253 } 

254 config.fiducialSkyBackground = { 

255 "u": 1.51, "g": 15.45, "r": 32.95, "i": 52.94, "z": 79.40, "y": 90.57, 

256 } 

257 config.fiducialReadNoise = { 

258 "u": 9.0, "g": 9.0, "r": 9.0, "i": 9.0, "z": 9.0, "y": 9.0, 

259 } 

260 config.fiducialExpTime = { 

261 "u": 30.0, "g": 30.0, "r": 30.0, "i": 30.0, "z": 30.0, "y": 30.0, 

262 } 

263 

264 # Other properties 

265 pixelScale = 0.2 # arcsec/pix 

266 gain = 1.0 # [e-/ADU] 

267 

268 # Fiducial m5 from SMTN-002 

269 fiducialMagLim = { 

270 "u": 23.70, "g": 24.97, "r": 24.52, "i": 24.13, "z": 23.56, "y": 22.55 

271 } 

272 

273 # Compare fiducial m5 calculated from config to SMTN-002 values 

274 for band in fiducialMagLim.keys(): 

275 print(f"\n{band}") 

276 m5fid = config.fiducialMagnitudeLimit(band, pixelScale, gain) 

277 print(f"{fiducialMagLim[band]:.2f} {m5fid:.3f}") 

278 self.assertFloatsAlmostEqual(m5fid, fiducialMagLim[band], atol=1e-2) 

279 

280 # Check that passing an unknown band outputs NaN 

281 nan = float('nan') 

282 m5fid = config.fiducialMagnitudeLimit("block", pixelScale, gain) 

283 self.assertFloatsAlmostEqual(m5fid, nan, ignoreNaNs=True) 

284 

285 

286class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase): 

287 pass 

288 

289 

290def setup_module(module): 

291 lsst.utils.tests.init() 

292 

293 

294if __name__ == "__main__": 294 ↛ 295line 294 didn't jump to line 295 because the condition on line 294 was never true

295 lsst.utils.tests.init() 

296 unittest.main()