Coverage for tests / test_photoCalib.py: 10%

1# 

2# LSST Data Management System 

3# Copyright 2008-2016 LSST Corporation. 

4# 

5# This product includes software developed by the 

6# LSST Project (http://www.lsst.org/). 

7# 

8# This program is free software: you can redistribute it and/or modify 

9# it under the terms of the GNU General Public License as published by 

10# the Free Software Foundation, either version 3 of the License, or 

11# (at your option) any later version. 

12# 

13# This program is distributed in the hope that it will be useful, 

14# but WITHOUT ANY WARRANTY; without even the implied warranty of 

15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

16# GNU General Public License for more details. 

17# 

18# You should have received a copy of the LSST License Statement and 

19# the GNU General Public License along with this program. If not, 

20# see <http://www.lsstcorp.org/LegalNotices/>. 

21# 

22 

23import os.path 

24import unittest 

25 

26import numpy as np 

27 

28import astropy.units as u 

29 

30import lsst.utils.tests 

31import lsst.geom 

32import lsst.afw.image 

33import lsst.afw.image.testUtils 

34import lsst.afw.math 

35import lsst.daf.base 

36import lsst.pex.exceptions 

37 

38 

39def computeNanojanskyErr(instFluxErr, calibration): 

40 """Return the error on the flux (nanojansky).""" 

41 return instFluxErr * calibration 

42 

43 

44def computeMagnitudeErr(instFluxErr, instFlux): 

45 """Return the error on the magnitude.""" 

46 return 2.5/np.log(10) * instFluxErr / instFlux 

47 

48 

49def makeCalibratedMaskedImage(image, mask, variance, calibration): 

50 """Return a MaskedImage that applies the given calibration to the given 

51 image, mask, and variance. 

52 """ 

53 return lsst.afw.image.makeMaskedImageFromArrays((image * calibration).astype(np.float32), 

54 mask, 

55 (variance * calibration**2).astype(np.float32)) 

56 

57 

58class PhotoCalibTestCase(lsst.utils.tests.TestCase): 

59 

60 def setUp(self): 

61 np.random.seed(100) 

62 

63 self.point0 = lsst.geom.Point2D(0, 0) 

64 self.pointXShift = lsst.geom.Point2D(-10, 0) 

65 self.pointYShift = lsst.geom.Point2D(0, -10) 

66 self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(-100, -100), lsst.geom.Point2I(100, 100)) 

67 

68 self.calibration = 1000.0 

69 # A 1% error on the calibration. 

70 self.calibrationErr = 10.0 

71 self.instFlux1 = 1.0 

72 self.instFluxErr1 = 0.1 

73 self.flux1 = 1000.0 # nJy 

74 self.mag1 = (self.flux1*u.nJy).to_value(u.ABmag) 

75 

76 # useful reference points: 575.44 nJy ~= 24.5 mag, 3630.78 * 10^9 nJy ~= 0 mag 

77 self.flux2 = 575.44 * self.flux1 

78 self.instFlux2 = self.flux2/self.calibration 

79 self.mag2 = (self.flux2*u.nJy).to_value(u.ABmag) 

80 

81 self.schema = lsst.afw.table.SourceTable.makeMinimalSchema() 

82 self.instFluxKeyName = "SomeFlux" 

83 lsst.afw.table.Point2DKey.addFields(self.schema, "centroid", "centroid", "pixels") 

84 self.schema.addField(self.instFluxKeyName+"_instFlux", type="D", units="count", 

85 doc="post-ISR instrumental Flux") 

86 self.schema.addField(self.instFluxKeyName+"_instFluxErr", type="D", units="count", 

87 doc="post-ISR instrumental flux error") 

88 self.schema.addField(self.instFluxKeyName+"_flux", type="D", units="nJy", 

89 doc="calibrated flux") 

90 self.schema.addField(self.instFluxKeyName+"_fluxErr", type="D", units="nJy", 

91 doc="calibrated flux error") 

92 self.schema.addField(self.instFluxKeyName+"_mag", type="D", 

93 doc="calibrated magnitude") 

94 self.schema.addField(self.instFluxKeyName+"_magErr", type="D", 

95 doc="calibrated magnitude error") 

96 self.otherInstFluxKeyName = "OtherFlux" 

97 self.schema.addField(self.otherInstFluxKeyName+"_instFlux", type="D", units="count", 

98 doc="another instrumental Flux") 

99 self.schema.addField(self.otherInstFluxKeyName+"_instFluxErr", type="D", units="count", 

100 doc="another instrumental flux error") 

101 self.noErrInstFluxKeyName = "NoErrFlux" 

102 self.schema.addField(self.noErrInstFluxKeyName+"_instFlux", type="D", units="count", 

103 doc="instrumental Flux with no error") 

104 self.table = lsst.afw.table.SourceTable.make(self.schema) 

105 self.table.defineCentroid('centroid') 

106 self.catalog = lsst.afw.table.SourceCatalog(self.table) 

107 record = self.catalog.addNew() 

108 record.set('id', 1) 

109 record.set('centroid_x', self.point0[0]) 

110 record.set('centroid_y', self.point0[1]) 

111 record.set(self.instFluxKeyName+'_instFlux', self.instFlux1) 

112 record.set(self.instFluxKeyName+'_instFluxErr', self.instFluxErr1) 

113 record.set(self.otherInstFluxKeyName+'_instFlux', self.instFlux1) 

114 record.set(self.otherInstFluxKeyName+'_instFluxErr', self.instFluxErr1) 

115 record.set(self.noErrInstFluxKeyName+'_instFlux', self.instFlux1) 

116 record = self.catalog.addNew() 

117 record.set('id', 2) 

118 record.set('centroid_x', self.pointYShift[0]) 

119 record.set('centroid_y', self.pointYShift[1]) 

120 record.set(self.instFluxKeyName+'_instFlux', self.instFlux2) 

121 record.set(self.instFluxKeyName+'_instFluxErr', self.instFluxErr1) 

122 record.set(self.otherInstFluxKeyName+'_instFlux', self.instFlux2) 

123 record.set(self.otherInstFluxKeyName+'_instFluxErr', self.instFluxErr1) 

124 record.set(self.noErrInstFluxKeyName+'_instFlux', self.instFlux2) 

125 

126 self.constantCalibration = lsst.afw.math.ChebyshevBoundedField(self.bbox, 

127 np.array([[self.calibration]])) 

128 self.linearXCalibration = lsst.afw.math.ChebyshevBoundedField(self.bbox, 

129 np.array([[self.calibration, 

130 self.calibration]])) 

131 

132 def tearDown(self): 

133 del self.schema 

134 del self.table 

135 del self.catalog 

136 

137 def _testPhotoCalibCenter(self, photoCalib, calibrationErr): 

138 """ 

139 Test conversions of instFlux for the mean and (0,0) value of a photoCalib. 

140 Assumes those are the same, e.g. that the non-constant terms are all 

141 odd, and that the mean of the calib is self.calibration. 

142 

143 calibrationErr is provided as an option to allow testing of photoCalibs 

144 that have no error specified, and those that do. 

145 """ 

146 # test that the constructor set the calibrationMean and err correctly 

147 self.assertEqual(self.calibration, photoCalib.getCalibrationMean()) 

148 self.assertEqual(photoCalib.instFluxToMagnitude(photoCalib.getInstFluxAtZeroMagnitude()), 0) 

149 self.assertEqual(calibrationErr, photoCalib.getCalibrationErr()) 

150 

151 # test with a "trivial" flux 

152 self.assertEqual(self.flux1, photoCalib.instFluxToNanojansky(self.instFlux1)) 

153 self.assertEqual(self.mag1, photoCalib.instFluxToMagnitude(self.instFlux1)) 

154 

155 # a less trivial flux 

156 self.assertFloatsAlmostEqual(self.flux2, photoCalib.instFluxToNanojansky(self.instFlux2)) 

157 self.assertFloatsAlmostEqual(self.mag2, photoCalib.instFluxToMagnitude(self.instFlux2)) 

158 # test that (0,0) gives the same result as above 

159 self.assertFloatsAlmostEqual(self.flux2, photoCalib.instFluxToNanojansky(self.instFlux2, self.point0)) 

160 self.assertFloatsAlmostEqual(self.mag2, photoCalib.instFluxToMagnitude(self.instFlux2, self.point0)) 

161 

162 # test that we get a correct nJy err for the base instFlux 

163 errFlux1 = computeNanojanskyErr(self.instFluxErr1, self.calibration) 

164 result = photoCalib.instFluxToNanojansky(self.instFlux1, self.instFluxErr1) 

165 self.assertEqual(self.flux1, result.value) 

166 self.assertFloatsAlmostEqual(errFlux1, result.error) 

167 result = photoCalib.instFluxToNanojansky(self.instFlux1, self.instFluxErr1, self.point0) 

168 self.assertFloatsAlmostEqual(self.flux1, result.value) 

169 self.assertFloatsAlmostEqual(errFlux1, result.error) 

170 

171 # test that we get a correct magnitude err for the base instFlux 

172 errMag1 = computeMagnitudeErr(self.instFluxErr1, self.instFlux1) 

173 result = photoCalib.instFluxToMagnitude(self.instFlux1, self.instFluxErr1) 

174 self.assertEqual(self.mag1, result.value) 

175 self.assertFloatsAlmostEqual(errMag1, result.error) 

176 # and the same given an explicit point at the center 

177 result = photoCalib.instFluxToMagnitude(self.instFlux1, self.instFluxErr1, self.point0) 

178 self.assertFloatsAlmostEqual(self.mag1, result.value) 

179 self.assertFloatsAlmostEqual(errMag1, result.error) 

180 

181 # test that we get a correct nJy err for flux2 

182 errFlux2 = computeNanojanskyErr(self.instFluxErr1, self.calibration) 

183 result = photoCalib.instFluxToNanojansky(self.instFlux2, self.instFluxErr1) 

184 self.assertFloatsAlmostEqual(self.flux2, result.value) 

185 self.assertFloatsAlmostEqual(errFlux2, result.error) 

186 result = photoCalib.instFluxToNanojansky(self.instFlux2, self.instFluxErr1, self.point0) 

187 self.assertFloatsAlmostEqual(self.flux2, result.value) 

188 self.assertFloatsAlmostEqual(errFlux2, result.error) 

189 

190 # test that we get a correct magnitude err for 575 nJy 

191 errMag2 = computeMagnitudeErr(self.instFluxErr1, self.instFlux2) 

192 result = photoCalib.instFluxToMagnitude(self.instFlux2, self.instFluxErr1) 

193 self.assertFloatsAlmostEqual(self.mag2, result.value) 

194 self.assertFloatsAlmostEqual(errMag2, result.error) 

195 result = photoCalib.instFluxToMagnitude(self.instFlux2, self.instFluxErr1, self.point0) 

196 self.assertFloatsAlmostEqual(self.mag2, result.value) 

197 self.assertFloatsAlmostEqual(errMag2, result.error) 

198 

199 # test calculations on a single sourceRecord 

200 record = self.catalog[0] 

201 result = photoCalib.instFluxToNanojansky(record, self.instFluxKeyName) 

202 self.assertEqual(self.flux1, result.value) 

203 self.assertFloatsAlmostEqual(errFlux1, result.error) 

204 result = photoCalib.instFluxToMagnitude(record, self.instFluxKeyName) 

205 self.assertEqual(self.mag1, result.value) 

206 self.assertFloatsAlmostEqual(errMag1, result.error) 

207 

208 expectNanojansky = np.array([[self.flux1, errFlux1], [self.flux2, errFlux2]]) 

209 expectMag = np.array([[self.mag1, errMag1], [self.mag2, errMag2]]) 

210 self._testSourceCatalog(photoCalib, self.catalog, expectNanojansky, expectMag) 

211 

212 # test reverse conversion: magnitude to instFlux (no position specified) 

213 self.assertFloatsAlmostEqual(self.instFlux1, photoCalib.magnitudeToInstFlux(self.mag1)) 

214 self.assertFloatsAlmostEqual(self.instFlux2, photoCalib.magnitudeToInstFlux(self.mag2), rtol=1e-15) 

215 

216 # test round-tripping instFlux->magnitude->instFlux (position specified) 

217 mag = photoCalib.instFluxToMagnitude(self.instFlux1, self.pointXShift) 

218 self.assertFloatsAlmostEqual(self.instFlux1, 

219 photoCalib.magnitudeToInstFlux(mag, self.pointXShift), 

220 rtol=1e-15) 

221 mag2 = photoCalib.instFluxToMagnitude(self.instFlux2, self.pointXShift) 

222 self.assertFloatsAlmostEqual(self.instFlux2, 

223 photoCalib.magnitudeToInstFlux(mag2, self.pointXShift), 

224 rtol=1e-15) 

225 

226 # test reverse conversion: nanojansky to instFlux (no position specified) 

227 self.assertFloatsAlmostEqual(self.instFlux1, photoCalib.nanojanskyToInstFlux(self.flux1)) 

228 self.assertFloatsAlmostEqual(self.instFlux2, photoCalib.nanojanskyToInstFlux(self.flux2), rtol=1e-15) 

229 

230 # test round-tripping instFlux->nanojansky->instFlux (position specified) 

231 flux = photoCalib.instFluxToNanojansky(self.instFlux1, self.pointXShift) 

232 self.assertFloatsAlmostEqual(self.instFlux1, 

233 photoCalib.nanojanskyToInstFlux(flux, self.pointXShift), 

234 rtol=1e-15) 

235 flux2 = photoCalib.instFluxToNanojansky(self.instFlux2, self.pointXShift) 

236 self.assertFloatsAlmostEqual(self.instFlux2, 

237 photoCalib.nanojanskyToInstFlux(flux2, self.pointXShift), 

238 rtol=1e-15) 

239 

240 # test round-tripping arrays (position specified) 

241 instFlux1Array = np.full(10, self.instFlux1) 

242 instFlux2Array = np.full(10, self.instFlux2) 

243 pointXShiftXArray = np.full(10, self.pointXShift.getX()) 

244 pointXShiftYArray = np.full(10, self.pointXShift.getY()) 

245 

246 magArray = photoCalib.instFluxToMagnitudeArray( 

247 instFlux1Array, 

248 pointXShiftXArray, 

249 pointXShiftYArray 

250 ) 

251 self.assertFloatsAlmostEqual(magArray.value, mag) 

252 self.assertFloatsAlmostEqual(photoCalib.magnitudeToInstFluxArray(magArray, 

253 pointXShiftXArray, 

254 pointXShiftYArray 

255 ), 

256 instFlux1Array, 

257 rtol=5e-15) 

258 mag2Array = photoCalib.instFluxToMagnitudeArray( 

259 np.full(10, self.instFlux2), 

260 np.full(10, self.pointXShift.getX()), 

261 np.full(10, self.pointXShift.getY()) 

262 ) 

263 self.assertFloatsAlmostEqual(mag2Array.value, mag2) 

264 self.assertFloatsAlmostEqual(photoCalib.magnitudeToInstFluxArray(mag2Array, 

265 pointXShiftXArray, 

266 pointXShiftYArray 

267 ), 

268 instFlux2Array, 

269 rtol=5e-15) 

270 

271 # test getLocalCalibration. 

272 meas = photoCalib.instFluxToNanojansky(self.instFlux1, self.instFluxErr1, self.pointXShift) 

273 localCalib = photoCalib.getLocalCalibration(self.pointXShift) 

274 flux = localCalib * self.instFlux1 

275 self.assertAlmostEqual(meas.value, flux) 

276 

277 # test getLocalCalibrationArray 

278 localCalib2 = photoCalib.getLocalCalibrationArray( 

279 pointXShiftXArray, 

280 pointXShiftYArray 

281 ) 

282 self.assertFloatsAlmostEqual(localCalib2, localCalib) 

283 

284 def _testSourceCatalog(self, photoCalib, catalog, expectNanojansky, expectMag): 

285 """Test instFluxTo*(sourceCatalog, ...), and calibrateCatalog().""" 

286 

287 # test calculations on a sourceCatalog, returning the array 

288 result = photoCalib.instFluxToNanojansky(catalog, self.instFluxKeyName) 

289 self.assertFloatsAlmostEqual(expectNanojansky, result) 

290 result = photoCalib.instFluxToMagnitude(catalog, self.instFluxKeyName) 

291 self.assertFloatsAlmostEqual(expectMag, result) 

292 

293 # Test modifying the catalog in-place with instFluxToNanojansky/instFluxToMagnitude 

294 # The original instFluxes shouldn't change: save them to test that. 

295 origInstFlux = catalog[self.instFluxKeyName+'_instFlux'].copy() 

296 origInstFluxErr = catalog[self.instFluxKeyName+'_instFluxErr'].copy() 

297 

298 def checkCatalog(catalog, expect, keyName, outField): 

299 """Test that the fields in the catalog are correct.""" 

300 self.assertFloatsAlmostEqual(catalog[keyName+'_%s' % outField], expect[:, 0]) 

301 self.assertFloatsAlmostEqual(catalog[keyName+'_%sErr' % outField], expect[:, 1]) 

302 self.assertFloatsAlmostEqual(catalog[keyName+'_instFlux'], origInstFlux) 

303 self.assertFloatsAlmostEqual(catalog[keyName+'_instFluxErr'], origInstFluxErr) 

304 

305 testCat = catalog.copy(deep=True) 

306 photoCalib.instFluxToMagnitude(testCat, self.instFluxKeyName, self.instFluxKeyName) 

307 checkCatalog(testCat, expectMag, self.instFluxKeyName, "mag") 

308 

309 testCat = catalog.copy(deep=True) 

310 photoCalib.instFluxToNanojansky(testCat, self.instFluxKeyName, self.instFluxKeyName) 

311 checkCatalog(testCat, expectNanojansky, self.instFluxKeyName, "flux") 

312 

313 testCat = catalog.copy(deep=True) 

314 photoCalib.instFluxToMagnitude(testCat, self.instFluxKeyName, self.instFluxKeyName) 

315 checkCatalog(testCat, expectMag, self.instFluxKeyName, "mag") 

316 

317 # test returning a calibrated catalog with calibrateCatalog 

318 

319 # test that trying to calibrate a non-existent flux field raises 

320 with self.assertRaises(lsst.pex.exceptions.NotFoundError): 

321 photoCalib.calibrateCatalog(testCat, ["NotARealFluxFieldName"]) 

322 

323 # test calibrating just one flux field 

324 testCat = catalog.copy(deep=True) 

325 result = photoCalib.calibrateCatalog(testCat, [self.otherInstFluxKeyName]) 

326 checkCatalog(result, expectNanojansky, self.otherInstFluxKeyName, "flux") 

327 checkCatalog(result, expectMag, self.otherInstFluxKeyName, "mag") 

328 

329 # test calibrating all of the flux fields 

330 testCat = catalog.copy(deep=True) 

331 result = photoCalib.calibrateCatalog(testCat) 

332 checkCatalog(result, expectNanojansky, self.instFluxKeyName, "flux") 

333 checkCatalog(result, expectMag, self.instFluxKeyName, "mag") 

334 checkCatalog(result, expectNanojansky, self.otherInstFluxKeyName, "flux") 

335 checkCatalog(result, expectMag, self.otherInstFluxKeyName, "mag") 

336 self.assertFloatsAlmostEqual(result[self.noErrInstFluxKeyName+'_flux'], expectNanojansky[:, 0]) 

337 self.assertFloatsAlmostEqual(result[self.noErrInstFluxKeyName+'_mag'], expectMag[:, 0]) 

338 self.assertFloatsAlmostEqual(result[self.noErrInstFluxKeyName+'_instFlux'], origInstFlux) 

339 

340 def testNonVarying(self): 

341 """Test constructing with a constant calibration factor.""" 

342 photoCalib = lsst.afw.image.PhotoCalib(self.calibration) 

343 self._testPhotoCalibCenter(photoCalib, 0) 

344 

345 # Test _isConstant 

346 self.assertTrue(photoCalib._isConstant) 

347 

348 # test on positions off the center (position should not matter) 

349 self.assertEqual(self.flux1, photoCalib.instFluxToNanojansky(self.instFlux1, self.pointXShift)) 

350 self.assertEqual(self.mag1, photoCalib.instFluxToMagnitude(self.instFlux1, self.pointXShift)) 

351 result = photoCalib.instFluxToNanojansky(self.instFlux1, self.instFluxErr1) 

352 self.assertEqual(self.flux1, result.value) 

353 

354 photoCalib = lsst.afw.image.PhotoCalib(self.calibration, self.calibrationErr) 

355 self._testPhotoCalibCenter(photoCalib, self.calibrationErr) 

356 

357 # test converting to a photoCalib 

358 photoCalib = lsst.afw.image.PhotoCalib(self.calibration, bbox=self.bbox) 

359 self._testPhotoCalibCenter(photoCalib, 0) 

360 

361 def testConstantBoundedField(self): 

362 """Test constructing with a spatially-constant bounded field.""" 

363 photoCalib = lsst.afw.image.PhotoCalib(self.constantCalibration) 

364 self._testPhotoCalibCenter(photoCalib, 0) 

365 

366 # test on positions off the center (position should not matter) 

367 self.assertEqual(self.flux1, photoCalib.instFluxToNanojansky(self.instFlux1, self.pointYShift)) 

368 self.assertEqual(self.mag1, photoCalib.instFluxToMagnitude(self.instFlux1, self.pointYShift)) 

369 self.assertFloatsAlmostEqual(self.flux2, 

370 photoCalib.instFluxToNanojansky(self.instFlux2, self.pointXShift)) 

371 self.assertFloatsAlmostEqual(self.mag2, 

372 photoCalib.instFluxToMagnitude(self.instFlux2, self.pointXShift)) 

373 

374 # test converting to a photoCalib 

375 photoCalib = lsst.afw.image.PhotoCalib(self.constantCalibration, self.calibrationErr) 

376 self._testPhotoCalibCenter(photoCalib, self.calibrationErr) 

377 

378 # test _isConstant (bounded field is not constant) 

379 self.assertFalse(photoCalib._isConstant) 

380 

381 def testLinearXBoundedField(self): 

382 photoCalib = lsst.afw.image.PhotoCalib(self.linearXCalibration) 

383 self._testPhotoCalibCenter(photoCalib, 0) 

384 

385 # test on positions off the center (Y position should not matter) 

386 self.assertEqual(self.flux1, photoCalib.instFluxToNanojansky(self.instFlux1, self.pointYShift)) 

387 self.assertEqual(self.mag1, photoCalib.instFluxToMagnitude(self.instFlux1, self.pointYShift)) 

388 

389 # test on positions off the center (X position does matter) 

390 calibration = (self.calibration + self.pointXShift.getX()*self.calibration/(self.bbox.getWidth()/2.)) 

391 expect = self.instFlux1*calibration 

392 self.assertFloatsAlmostEqual(expect, 

393 photoCalib.instFluxToNanojansky(self.instFlux1, self.pointXShift)) 

394 self.assertFloatsAlmostEqual((expect*u.nJy).to_value(u.ABmag), 

395 photoCalib.instFluxToMagnitude(self.instFlux1, self.pointXShift)) 

396 expect2 = self.instFlux2*calibration 

397 self.assertFloatsAlmostEqual(expect2, 

398 photoCalib.instFluxToNanojansky(self.instFlux2, self.pointXShift)) 

399 self.assertFloatsAlmostEqual((expect2*u.nJy).to_value(u.ABmag), 

400 photoCalib.instFluxToMagnitude(self.instFlux2, self.pointXShift)) 

401 

402 # test converting to a photoCalib 

403 photoCalib = lsst.afw.image.PhotoCalib(self.linearXCalibration, self.calibrationErr) 

404 self._testPhotoCalibCenter(photoCalib, self.calibrationErr) 

405 

406 # New catalog with a spatial component in the varying direction, 

407 # to ensure the calculations on a catalog properly handle non-constant BF. 

408 # NOTE: only the first quantity of the result (nJy or mags) should change. 

409 catalog = self.catalog.copy(deep=True) 

410 catalog[0].set('centroid_x', self.pointXShift[0]) 

411 catalog[0].set('centroid_y', self.pointXShift[1]) 

412 errFlux1 = computeNanojanskyErr(self.instFluxErr1, calibration) 

413 errMag1 = computeMagnitudeErr(self.instFluxErr1, self.instFlux1) 

414 # re-use the same instFluxErr1 for instFlux2. 

415 errFlux2 = computeNanojanskyErr(self.instFluxErr1, self.calibration) 

416 errMag2 = computeMagnitudeErr(self.instFluxErr1, self.instFlux2) 

417 expectNanojansky = np.array([[expect, errFlux1], [self.flux2, errFlux2]]) 

418 expectMag = np.array([[(expect*u.nJy).to_value(u.ABmag), errMag1], [self.mag2, errMag2]]) 

419 self._testSourceCatalog(photoCalib, catalog, expectNanojansky, expectMag) 

420 

421 def testComputeScaledCalibration(self): 

422 photoCalib = lsst.afw.image.PhotoCalib(self.calibration, bbox=self.bbox) 

423 scaledCalib = lsst.afw.image.PhotoCalib(photoCalib.computeScaledCalibration()) 

424 self.assertEqual(self.flux1, 

425 scaledCalib.instFluxToNanojansky(self.instFlux1)*photoCalib.getCalibrationMean()) 

426 self.assertEqual(photoCalib.instFluxToNanojansky(self.instFlux1), 

427 scaledCalib.instFluxToNanojansky(self.instFlux1)*photoCalib.getCalibrationMean()) 

428 

429 photoCalib = lsst.afw.image.PhotoCalib(self.constantCalibration) 

430 scaledCalib = lsst.afw.image.PhotoCalib(photoCalib.computeScaledCalibration()) 

431 

432 self.assertEqual(self.flux1, scaledCalib.instFluxToNanojansky(self.instFlux1*self.calibration)) 

433 self.assertEqual(photoCalib.instFluxToNanojansky(self.instFlux1), 

434 scaledCalib.instFluxToNanojansky(self.instFlux1)*photoCalib.getCalibrationMean()) 

435 

436 @unittest.skip("Not yet implemented: see DM-10154") 

437 def testComputeScalingTo(self): 

438 photoCalib1 = lsst.afw.image.PhotoCalib(self.calibration, self.calibrationErr, bbox=self.bbox) 

439 photoCalib2 = lsst.afw.image.PhotoCalib(self.calibration*500, self.calibrationErr, bbox=self.bbox) 

440 scaling = photoCalib1.computeScalingTo(photoCalib2)(self.pointXShift) 

441 self.assertEqual(photoCalib1.instFluxToNanojansky(self.instFlux1, self.pointXShift)*scaling, 

442 photoCalib2.instFluxToNanojansky(self.instFlux1, self.pointXShift)) 

443 

444 photoCalib3 = lsst.afw.image.PhotoCalib(self.constantCalibration, self.calibrationErr) 

445 scaling = photoCalib1.computeScalingTo(photoCalib3)(self.pointXShift) 

446 self.assertEqual(photoCalib1.instFluxToNanojansky(self.instFlux1, self.pointXShift)*scaling, 

447 photoCalib3.instFluxToNanojansky(self.instFlux1, self.pointXShift)) 

448 scaling = photoCalib3.computeScalingTo(photoCalib1)(self.pointXShift) 

449 self.assertEqual(photoCalib3.instFluxToNanojansky(self.instFlux1, self.pointXShift)*scaling, 

450 photoCalib1.instFluxToNanojansky(self.instFlux1, self.pointXShift)) 

451 

452 photoCalib4 = lsst.afw.image.PhotoCalib(self.linearXCalibration, self.calibrationErr) 

453 scaling = photoCalib1.computeScalingTo(photoCalib4)(self.pointXShift) 

454 self.assertEqual(photoCalib1.instFluxToNanojansky(self.instFlux1, self.pointXShift)*scaling, 

455 photoCalib4.instFluxToNanojansky(self.instFlux1, self.pointXShift)) 

456 scaling = photoCalib4.computeScalingTo(photoCalib1)(self.pointXShift) 

457 self.assertEqual(photoCalib4.instFluxToNanojansky(self.instFlux1, self.pointXShift)*scaling, 

458 photoCalib1.instFluxToNanojansky(self.instFlux1, self.pointXShift)) 

459 

460 # Don't allow division of BoundedFields with different bounding boxes 

461 photoCalibNoBBox = lsst.afw.image.PhotoCalib(self.calibration, self.calibrationErr) 

462 with self.assertRaises(lsst.pex.exceptions.DomainError): 

463 scaling = photoCalibNoBBox.computeScalingTo(photoCalib1) 

464 with self.assertRaises(lsst.pex.exceptions.DomainError): 

465 scaling = photoCalibNoBBox.computeScalingTo(photoCalib4) 

466 with self.assertRaises(lsst.pex.exceptions.DomainError): 

467 scaling = photoCalib1.computeScalingTo(photoCalibNoBBox) 

468 

469 def _testPersistence(self, photoCalib): 

470 with lsst.utils.tests.getTempFilePath(".fits") as filename: 

471 photoCalib.writeFits(filename) 

472 result = lsst.afw.image.PhotoCalib.readFits(filename) 

473 self.assertEqual(result, photoCalib) 

474 

475 def testPersistence(self): 

476 photoCalib = lsst.afw.image.PhotoCalib(self.calibration) 

477 self._testPersistence(photoCalib) 

478 

479 photoCalib = lsst.afw.image.PhotoCalib(self.calibration, self.calibrationErr) 

480 self._testPersistence(photoCalib) 

481 

482 photoCalib = lsst.afw.image.PhotoCalib(self.calibration, self.calibrationErr, self.bbox) 

483 self._testPersistence(photoCalib) 

484 

485 photoCalib = lsst.afw.image.PhotoCalib(self.constantCalibration, self.calibrationErr) 

486 self._testPersistence(photoCalib) 

487 

488 photoCalib = lsst.afw.image.PhotoCalib(self.linearXCalibration, self.calibrationErr) 

489 self._testPersistence(photoCalib) 

490 

491 def testPersistenceVersions(self): 

492 """Test that different versions are handled appropriately.""" 

493 # the values that were persisted in this photoCalib 

494 mean = 123 

495 err = 45 

496 dataDir = os.path.join(os.path.split(__file__)[0], "data") 

497 

498 # implicit version 0 should raise (no longer compatible) 

499 filePath = os.path.join(dataDir, "photoCalib-noversion.fits") 

500 with self.assertRaises(RuntimeError): 

501 photoCalib = lsst.afw.image.PhotoCalib.readFits(filePath) 

502 

503 # explicit version 0 should raise (no longer compatible) 

504 filePath = os.path.join(dataDir, "photoCalib-version0.fits") 

505 with self.assertRaises(RuntimeError): 

506 photoCalib = lsst.afw.image.PhotoCalib.readFits(filePath) 

507 

508 # explicit version 1 

509 filePath = os.path.join(dataDir, "photoCalib-version1.fits") 

510 photoCalib = lsst.afw.image.PhotoCalib.readFits(filePath) 

511 self.assertEqual(photoCalib.getCalibrationMean(), mean) 

512 self.assertEqual(photoCalib.getCalibrationErr(), err) 

513 

514 def testPhotoCalibEquality(self): 

515 photoCalib1 = lsst.afw.image.PhotoCalib(self.linearXCalibration, 0.5) 

516 photoCalib2 = lsst.afw.image.PhotoCalib(self.linearXCalibration, 0.5) 

517 photoCalib3 = lsst.afw.image.PhotoCalib(5, 0.5) 

518 photoCalib4 = lsst.afw.image.PhotoCalib(5, 0.5) 

519 photoCalib5 = lsst.afw.image.PhotoCalib(5) 

520 photoCalib6 = lsst.afw.image.PhotoCalib(self.linearXCalibration) 

521 photoCalib7 = lsst.afw.image.PhotoCalib(self.calibration, 0.5) 

522 photoCalib8 = lsst.afw.image.PhotoCalib(self.constantCalibration, 0.5) 

523 

524 constantCalibration = lsst.afw.math.ChebyshevBoundedField(self.bbox, np.array([[self.calibration]])) 

525 photoCalib9 = lsst.afw.image.PhotoCalib(constantCalibration, 0.5) 

526 

527 self.assertEqual(photoCalib1, photoCalib1) 

528 self.assertEqual(photoCalib1, photoCalib2) 

529 self.assertEqual(photoCalib3, photoCalib4) 

530 self.assertEqual(photoCalib8, photoCalib9) 

531 

532 self.assertNotEqual(photoCalib1, photoCalib6) 

533 self.assertNotEqual(photoCalib1, photoCalib7) 

534 self.assertNotEqual(photoCalib1, photoCalib3) 

535 self.assertNotEqual(photoCalib3, photoCalib5) 

536 self.assertNotEqual(photoCalib1, photoCalib8) 

537 

538 self.assertFalse(photoCalib1 != photoCalib2) # using assertFalse to directly test != operator 

539 

540 def setupImage(self): 

541 dim = (5, 6) 

542 npDim = (dim[1], dim[0]) # numpy and afw have a different x/y order 

543 sigma = 10.0 

544 # An image with increasing pixel values, to more easily check the 

545 # values in specific calculations. 

546 image = np.arange(dim[0]*dim[1]).astype(np.float32).reshape(npDim) + 1000 

547 mask = np.zeros(npDim, dtype=np.int32) 

548 variance = np.ones(npDim, dtype=np.float32)*sigma 

549 maskedImage = lsst.afw.image.makeMaskedImageFromArrays(image, mask, variance) 

550 maskedImage.mask[0, 0] = True # set one mask bit to check propagation of mask bits. 

551 

552 return npDim, maskedImage, image, mask, variance 

553 

554 def testCalibrateImageConstant(self): 

555 """Test a spatially-constant calibration.""" 

556 _, maskedImage, image, mask, variance = self.setupImage() 

557 photoCalib = lsst.afw.image.PhotoCalib(self.calibration, self.calibrationErr) 

558 expect = makeCalibratedMaskedImage(image, mask, variance, self.calibration) 

559 result = photoCalib.calibrateImage(maskedImage) 

560 self.assertMaskedImagesAlmostEqual(expect, result) 

561 uncalibrated = photoCalib.uncalibrateImage(result) 

562 self.assertMaskedImagesAlmostEqual(maskedImage, uncalibrated) 

563 

564 def testCalibrateImageNonConstant(self): 

565 """Test a spatially-varying calibration.""" 

566 npDim, maskedImage, image, mask, variance = self.setupImage() 

567 xIndex, yIndex = np.indices(npDim, dtype=np.float64) 

568 # y then x, as afw order and np order are flipped 

569 calibration = self.linearXCalibration.evaluate(yIndex.flatten(), xIndex.flatten()).reshape(npDim) 

570 expect = makeCalibratedMaskedImage(image, mask, variance, calibration) 

571 photoCalib = lsst.afw.image.PhotoCalib(self.linearXCalibration, self.calibrationErr) 

572 result = photoCalib.calibrateImage(maskedImage) 

573 self.assertMaskedImagesAlmostEqual(expect, result) 

574 uncalibrated = photoCalib.uncalibrateImage(result) 

575 self.assertMaskedImagesAlmostEqual(maskedImage, uncalibrated) 

576 

577 def testCalibrateImageNonConstantSubimage(self): 

578 """Test a non-constant calibration on a sub-image, to ensure we're 

579 handling xy0 correctly. 

580 """ 

581 npDim, maskedImage, image, mask, variance = self.setupImage() 

582 xIndex, yIndex = np.indices(npDim, dtype=np.float64) 

583 calibration = self.linearXCalibration.evaluate(yIndex.flatten(), xIndex.flatten()).reshape(npDim) 

584 

585 expect = makeCalibratedMaskedImage(image, mask, variance, calibration) 

586 

587 subBox = lsst.geom.Box2I(lsst.geom.Point2I(2, 4), lsst.geom.Point2I(4, 5)) 

588 subImage = maskedImage.subset(subBox) 

589 photoCalib = lsst.afw.image.PhotoCalib(self.linearXCalibration, self.calibrationErr) 

590 result = photoCalib.calibrateImage(subImage) 

591 self.assertMaskedImagesAlmostEqual(expect.subset(subBox), result) 

592 uncalibrated = photoCalib.uncalibrateImage(result) 

593 self.assertMaskedImagesAlmostEqual(subImage, uncalibrated) 

594 

595 def testNonPositiveMeans(self): 

596 # no negative calibrations 

597 with self.assertRaises(lsst.pex.exceptions.InvalidParameterError): 

598 lsst.afw.image.PhotoCalib(-1.0) 

599 # no negative errors