Coverage for tests/test_computeExposureSummaryStats.py: 19%

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 compute_magnitude_limit 

38 

39# Values from syseng_throughputs notebook assuming 30s exposure 

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

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

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

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

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

45# Output magnitude limit from syseng_throughputs notebook 

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

47 

48 

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

50 

51 def testComputeExposureSummary(self): 

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

53 """ 

54 np.random.seed(12345) 

55 

56 # Make an exposure with a noise image 

57 exposure = afwImage.ExposureF(100, 100) 

58 

59 band = "i" 

60 physical_filter = "test-i" 

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

62 

63 readNoise = 5.0 

64 gain = 1.0 

65 detector = DetectorWrapper(numAmps=1).detector 

66 metadata = PropertyList() 

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

68 for amp in detector.getAmplifiers(): 

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

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

71 exposure.setDetector(detector) 

72 exposure.setMetadata(metadata) 

73 

74 skyMean = 100.0 

75 skySigma = 10.0 

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

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

78 

79 # Set the visitInfo 

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

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

82 2650.) 

83 expTime = 10.0 

84 visitInfo = afwImage.VisitInfo(exposureTime=expTime, 

85 date=date, 

86 observatory=observatory) 

87 exposure.getInfo().setVisitInfo(visitInfo) 

88 

89 # Install a Gaussian PSF 

90 psfSize = 2.0 

91 psf = GaussianPsf(5, 5, psfSize) 

92 exposure.setPsf(psf) 

93 

94 # Install a simple WCS 

95 scale = 0.2*lsst.geom.arcseconds 

96 raCenter = 300.0*lsst.geom.degrees 

97 decCenter = 0.0*lsst.geom.degrees 

98 cdMatrix = makeCdMatrix(scale=scale) 

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

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

101 cdMatrix=cdMatrix) 

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

103 exposure.setWcs(skyWcs) 

104 

105 # Install a simple photoCalib 

106 photoCalib = afwImage.PhotoCalib(calibrationMean=0.3) 

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

108 exposure.setPhotoCalib(photoCalib) 

109 

110 # Install a simple apCorrMap 

111 apCorrMap = afwImage.ApCorrMap() 

112 apCorrMap.set( 

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

114 ) 

115 exposure.setApCorrMap(apCorrMap) 

116 

117 # Compute the background image 

118 bgGridSize = 10 

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

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

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

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

123 background = afwMath.BackgroundList() 

124 background.append(backobj) 

125 

126 # Configure and run the task 

127 expSummaryTask = ComputeExposureSummaryStatsTask() 

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

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

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

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

132 # Run the task 

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

134 

135 # Test the outputs 

136 self.assertFloatsAlmostEqual(summary.psfSigma, psfSize) 

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

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

139 self.assertFloatsAlmostEqual(summary.psfIxy, 0.0) 

140 self.assertFloatsAlmostEqual(summary.psfArea, 23.088975164455444) 

141 

142 self.assertFloatsAlmostEqual(summary.psfTraceRadiusDelta, 0.0) 

143 self.assertFloatsAlmostEqual(summary.psfApFluxDelta, 0.0) 

144 self.assertFloatsAlmostEqual(summary.psfApCorrSigmaScaledDelta, 0.0) 

145 

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

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

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

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

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

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

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

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

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

155 

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

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

158 

159 self.assertFloatsAlmostEqual(summary.pixelScale, pixelScale) 

160 

161 self.assertFloatsAlmostEqual(summary.zeroPoint, zp) 

162 self.assertFloatsAlmostEqual(summary.expTime, expTime) 

163 

164 # Need to compare background level and noise 

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

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

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

168 

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

170 

171 # Effective exposure time and depth 

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

173 self.assertFloatsAlmostEqual(summary.magLim, 26.584, rtol=1e-3) 

174 

175 def testComputeMagnitudeLimit(self): 

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

177 and syseng_throughputs.""" 

178 

179 # Assumed values from SMTN-002 

180 snr = 5 

181 gain = 1.0 

182 # For testing non-unity gain 

183 gain_test = 2.0 

184 

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

186 # Translate into DM quantities 

187 psfArea = 2.266 * (fwhm_eff_fid[band] / 0.2)**2 

188 skyBg = skycounts_fid[band] 

189 zeroPoint = zeropoint_fid[band] + 2.5*np.log10(30) 

190 readNoise = readnoise_fid[band] 

191 

192 # Calculate the M5 values 

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

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

195 

196 # Changing the gain consistently should not change M5 

197 m5_gain = compute_magnitude_limit(psfArea, skyBg/gain_test, 

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

199 readNoise/gain_test, gain_test, snr) 

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

201 

202 # Check that input NaN lead to output NaN 

203 nan = float('nan') 

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

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

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

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

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

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

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

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

212 

213 

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

215 pass 

216 

217 

218def setup_module(module): 

219 lsst.utils.tests.init() 

220 

221 

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

223 lsst.utils.tests.init() 

224 unittest.main()