Coverage for tests / test_computeExposureSummaryStats.py: 17%

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 nan = float("nan") 

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

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

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

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

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

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

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

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

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

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

156 

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

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

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

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

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

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

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

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

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

166 

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

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

169 

170 self.assertFloatsAlmostEqual(summary.pixelScale, pixelScale) 

171 

172 self.assertFloatsAlmostEqual(summary.zeroPoint, zp) 

173 self.assertFloatsAlmostEqual(summary.expTime, expTime) 

174 

175 # Need to compare background level and noise 

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

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

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

179 

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

181 

182 # Effective exposure time and depth 

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

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

185 

186 def testComputeMagnitudeLimit(self): 

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

188 and syseng_throughputs.""" 

189 

190 # Assumed values from SMTN-002 

191 snr = 5 

192 gain = 1.0 

193 # For testing non-unity gain 

194 gain_test = 2.0 

195 

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

197 # Translate into DM quantities 

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

199 skyBg = skycounts_fid[band] 

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

201 readNoise = readnoise_fid[band] 

202 

203 # Calculate the M5 values 

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

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

206 

207 # Changing the gain consistently should not change M5 

208 m5_gain = compute_magnitude_limit(psfArea, skyBg/gain_test, 

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

210 readNoise/gain_test, gain_test, snr) 

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

212 

213 # Check that input NaN lead to output NaN 

214 nan = float('nan') 

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

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

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

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

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

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

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

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

223 

224 

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

226 pass 

227 

228 

229def setup_module(module): 

230 lsst.utils.tests.init() 

231 

232 

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

234 lsst.utils.tests.init() 

235 unittest.main()