Coverage for tests/test_computeExposureSummaryStats.py: 15%

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 

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

41 

42 def testComputeExposureSummary(self): 

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

44 """ 

45 np.random.seed(12345) 

46 

47 # Make an exposure with a noise image 

48 exposure = afwImage.ExposureF(100, 100) 

49 

50 band = "i" 

51 physical_filter = "test-i" 

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

53 

54 readNoise = 5.0 

55 gain = 1.0 

56 detector = DetectorWrapper(numAmps=1).detector 

57 metadata = PropertyList() 

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

59 for amp in detector.getAmplifiers(): 

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

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

62 exposure.setDetector(detector) 

63 exposure.setMetadata(metadata) 

64 

65 skyMean = 100.0 

66 skySigma = 10.0 

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

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

69 

70 # Set the visitInfo 

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

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

73 2650.) 

74 expTime = 10.0 

75 visitInfo = afwImage.VisitInfo(exposureTime=expTime, 

76 date=date, 

77 observatory=observatory) 

78 exposure.getInfo().setVisitInfo(visitInfo) 

79 

80 # Install a Gaussian PSF 

81 psfSize = 2.0 

82 psf = GaussianPsf(5, 5, psfSize) 

83 exposure.setPsf(psf) 

84 

85 # Install a simple WCS 

86 scale = 0.2*lsst.geom.arcseconds 

87 raCenter = 300.0*lsst.geom.degrees 

88 decCenter = 0.0*lsst.geom.degrees 

89 cdMatrix = makeCdMatrix(scale=scale) 

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

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

92 cdMatrix=cdMatrix) 

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

94 exposure.setWcs(skyWcs) 

95 

96 # Install a simple photoCalib 

97 photoCalib = afwImage.PhotoCalib(calibrationMean=0.3) 

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

99 exposure.setPhotoCalib(photoCalib) 

100 

101 # Install a simple apCorrMap 

102 apCorrMap = afwImage.ApCorrMap() 

103 apCorrMap.set( 

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

105 ) 

106 exposure.setApCorrMap(apCorrMap) 

107 

108 # Compute the background image 

109 bgGridSize = 10 

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

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

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

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

114 background = afwMath.BackgroundList() 

115 background.append(backobj) 

116 

117 # Configure and run the task 

118 expSummaryTask = ComputeExposureSummaryStatsTask() 

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

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

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

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

123 # Run the task 

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

125 

126 # Test the outputs 

127 self.assertFloatsAlmostEqual(summary.psfSigma, psfSize) 

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

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

130 self.assertFloatsAlmostEqual(summary.psfIxy, 0.0) 

131 self.assertFloatsAlmostEqual(summary.psfArea, 23.088975164455444) 

132 

133 self.assertFloatsAlmostEqual(summary.psfTraceRadiusDelta, 0.0) 

134 self.assertFloatsAlmostEqual(summary.psfApFluxDelta, 0.0) 

135 self.assertFloatsAlmostEqual(summary.psfApCorrSigmaScaledDelta, 0.0) 

136 

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

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

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

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

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

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

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

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

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

146 

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

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

149 

150 self.assertFloatsAlmostEqual(summary.pixelScale, pixelScale) 

151 

152 self.assertFloatsAlmostEqual(summary.zeroPoint, zp) 

153 self.assertFloatsAlmostEqual(summary.expTime, expTime) 

154 

155 # Need to compare background level and noise 

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

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

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

159 

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

161 

162 # Effective exposure time and depth 

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

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

165 

166 def testComputeMagnitudeLimit(self): 

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

168 and syseng_throughputs.""" 

169 

170 # Values from syseng_throughputs notebook assuming 30s exposure 

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

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

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

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

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

176 # Assumed values from SMTN-002 

177 snr = 5 

178 gain = 1.0 

179 # Output magnitude limit from syseng_throughputs notebook 

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

181 

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

183 # Translate into DM quantities 

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

185 skyBg = skycounts_fid[band] 

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

187 readNoise = readnoise_fid[band] 

188 

189 # Calculate the M5 values 

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

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

192 

193 # Check that input NaN lead to output NaN 

194 nan = float('nan') 

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

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

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

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

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

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

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

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

203 

204 

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

206 pass 

207 

208 

209def setup_module(module): 

210 lsst.utils.tests.init() 

211 

212 

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

214 lsst.utils.tests.init() 

215 unittest.main()