Coverage for tests/test_calibrateImage.py: 14%

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 

22import unittest 

23from unittest import mock 

24import tempfile 

25 

26import astropy.units as u 

27from astropy.coordinates import SkyCoord 

28import copy 

29import numpy as np 

30import esutil 

31import os 

32import requests 

33 

34import lsst.afw.image as afwImage 

35import lsst.afw.math as afwMath 

36import lsst.afw.table as afwTable 

37import lsst.daf.base 

38import lsst.daf.butler 

39import lsst.daf.butler.tests as butlerTests 

40import lsst.geom 

41import lsst.meas.algorithms 

42from lsst.meas.algorithms import testUtils 

43import lsst.meas.extensions.psfex 

44import lsst.meas.base 

45import lsst.meas.base.tests 

46import lsst.pipe.base.testUtils 

47from lsst.pipe.tasks.calibrateImage import CalibrateImageTask, \ 

48 NoPsfStarsToStarsMatchError, AllCentroidsFlaggedError 

49import lsst.pex.config as pexConfig 

50import lsst.utils.tests 

51 

52from utils import makeTestVisitInfo 

53 

54 

55class CalibrateImageTaskTests(lsst.utils.tests.TestCase): 

56 

57 def setUp(self): 

58 # Different x/y dimensions so they're easy to distinguish in a plot, 

59 # and non-zero minimum, to help catch xy0 errors. 

60 bbox = lsst.geom.Box2I(lsst.geom.Point2I(5, 4), lsst.geom.Point2I(205, 184)) 

61 self.sky_center = lsst.geom.SpherePoint(245.0, -45.0, lsst.geom.degrees) 

62 self.photo_calib = 12.3 

63 dataset = lsst.meas.base.tests.TestDataset(bbox, crval=self.sky_center, calibration=self.photo_calib) 

64 # sqrt of area of a normalized 2d gaussian 

65 psf_scale = np.sqrt(4*np.pi*(dataset.psfShape.getDeterminantRadius())**2) 

66 noise = 10.0 # stddev of noise per pixel 

67 # Sources ordered from faintest to brightest. 

68 self.fluxes = np.array((6*noise*psf_scale, 

69 12*noise*psf_scale, 

70 45*noise*psf_scale, 

71 150*noise*psf_scale, 

72 400*noise*psf_scale, 

73 1000*noise*psf_scale)) 

74 self.centroids = np.array(((162, 22), 

75 (40, 70), 

76 (100, 160), 

77 (50, 120), 

78 (92, 35), 

79 (175, 154)), dtype=np.float32) 

80 for flux, centroid in zip(self.fluxes, self.centroids): 

81 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(centroid[0], centroid[1])) 

82 

83 # Bright extended source in the center of the image: should not appear 

84 # in any of the output catalogs. 

85 self.extended_source = lsst.geom.Point2D(100, 100) 

86 shape = lsst.afw.geom.Quadrupole(8, 9, 3) 

87 dataset.addSource(instFlux=500*noise*psf_scale, centroid=self.extended_source, shape=shape) 

88 

89 schema = dataset.makeMinimalSchema() 

90 self.truth_exposure, self.truth_cat = dataset.realize(noise=noise, schema=schema) 

91 # Add in a significant background, so we can test that the output 

92 # background is self-consistent with the calibrated exposure. 

93 self.background_level = 500.0 

94 self.truth_exposure.image += self.background_level 

95 # To make it look like a version=1 (nJy fluxes) refcat 

96 self.truth_cat = self.truth_exposure.photoCalib.calibrateCatalog(self.truth_cat) 

97 self.ref_loader = testUtils.MockReferenceObjectLoaderFromMemory([self.truth_cat]) 

98 metadata = lsst.daf.base.PropertyList() 

99 metadata.set("REFCAT_FORMAT_VERSION", 1) 

100 self.truth_cat.setMetadata(metadata) 

101 

102 # TODO: a cosmic ray (need to figure out how to insert a fake-CR) 

103 # self.truth_exposure.image.array[10, 10] = 100000 

104 # self.truth_exposure.variance.array[10, 10] = 100000/noise 

105 

106 # Copy the truth exposure, because CalibrateImage modifies the input. 

107 # Post-ISR ccds only contain: initial WCS, VisitInfo, filter 

108 self.exposure = afwImage.ExposureF(self.truth_exposure, deep=True) 

109 self.exposure.setWcs(self.truth_exposure.wcs) 

110 self.exposure.info.setVisitInfo(self.truth_exposure.visitInfo) 

111 self.exposure.info.id = 12345 

112 # "truth" filter, to match the "truth" refcat. 

113 self.exposure.setFilter(lsst.afw.image.FilterLabel(physical='truth', band="truth")) 

114 self.exposure.metadata["LSST ISR FLAT APPLIED"] = True 

115 

116 # Test-specific configuration: 

117 self.config = CalibrateImageTask.ConfigClass() 

118 # We don't have many sources, so have to fit simpler models. 

119 self.config.psf_detection.background.approxOrderX = 1 

120 self.config.star_detection.background.approxOrderX = 1 

121 # Only insert 2 sky sources, for simplicity. 

122 self.config.star_sky_sources.nSources = 2 

123 # Use PCA psf fitter, as psfex fails if there are only 4 stars. 

124 self.config.psf_measure_psf.psfDeterminer = 'pca' 

125 # We don't have many test points, so can't match on complicated shapes. 

126 self.config.astrometry.sourceSelector["science"].flags.good = [] 

127 self.config.astrometry.matcher.numPointsForShape = 3 

128 self.config.run_sattle = False 

129 # Maintain original, no adaptive threshold detection, configs values. 

130 self.config.do_adaptive_threshold_detection = False 

131 self.config.psf_detection.reEstimateBackground = True 

132 self.config.star_detection.reEstimateBackground = True 

133 # ApFlux has more noise than PsfFlux (the latter unrealistically small 

134 # in this test data), so we need to do magnitude rejection at higher 

135 # sigma, otherwise we can lose otherwise good sources. 

136 # TODO DM-39203: Once we are using Compensated Gaussian Fluxes, we 

137 # will use those fluxes here, and hopefully can remove this. 

138 self.config.astrometry.magnitudeOutlierRejectionNSigma = 9.0 

139 

140 # Make a realistic id generator so that output catalog ids are useful. 

141 # NOTE: The id generator is used to seed the noise replacer during 

142 # measurement, so changes to values here can have subtle effects on 

143 # the centroids and fluxes mesaured on the image, which might cause 

144 # tests to fail. 

145 data_id = lsst.daf.butler.DataCoordinate.standardize( 

146 instrument="I", 

147 visit=self.truth_exposure.visitInfo.id, 

148 detector=12, 

149 universe=lsst.daf.butler.DimensionUniverse(), 

150 ) 

151 self.config.id_generator.packer.name = "observation" 

152 self.config.id_generator.packer["observation"].n_observations = 10000 

153 self.config.id_generator.packer["observation"].n_detectors = 99 

154 self.config.id_generator.n_releases = 8 

155 self.config.id_generator.release_id = 2 

156 self.id_generator = self.config.id_generator.apply(data_id) 

157 

158 # Something about this test dataset prefers a larger threshold here. 

159 self.config.star_selector["science"].unresolved.maximum = 0.2 

160 

161 def _check_run(self, calibrate, result, expect_calibrated_pixels: bool = True, 

162 expect_n_background: int = 4): 

163 """Test the result of CalibrateImage.run(). 

164 

165 Parameters 

166 ---------- 

167 calibrate : `lsst.pipe.tasks.calibrateImage.CalibrateImageTask` 

168 Configured task that had `run` called on it. 

169 result : `lsst.pipe.base.Struct` 

170 Result of calling calibrate.run(). 

171 expect_calibrated_pixels : `bool`, optional 

172 Whether to expect image and background pixels to be calibrated. 

173 """ 

174 # Background should have 4 elements: 3 from compute_psf and one from 

175 # re-estimation during source detection. 

176 self.assertEqual(len(result.background), expect_n_background) 

177 

178 # Both afw and astropy psf_stars catalogs should be populated. 

179 self.assertEqual(result.psf_stars["calib_psf_used"].sum(), 3) 

180 self.assertEqual(result.psf_stars_footprints["calib_psf_used"].sum(), 3) 

181 

182 # Check that the summary statistics are reasonable. 

183 summary = result.exposure.info.getSummaryStats() 

184 self.assertFloatsAlmostEqual(summary.psfSigma, 2.0, rtol=1e-2) 

185 self.assertFloatsAlmostEqual(summary.ra, self.sky_center.getRa().asDegrees(), rtol=1e-7) 

186 self.assertFloatsAlmostEqual(summary.dec, self.sky_center.getDec().asDegrees(), rtol=1e-7) 

187 

188 # Should have finite sky coordinates in the afw and astropy catalogs. 

189 self.assertTrue(np.isfinite(result.stars_footprints["coord_ra"]).all()) 

190 self.assertTrue(np.isfinite(result.stars["coord_ra"]).all()) 

191 

192 if expect_calibrated_pixels: 

193 # Fit photoCalib should be the applied value if we calibrated 

194 # pixels, not the ==1 one on the exposure. 

195 photo_calib = result.applied_photo_calib 

196 self.assertEqual(result.exposure.photoCalib.getCalibrationMean(), 1.0) 

197 else: 

198 self.assertIsNone(result.applied_photo_calib) 

199 photo_calib = result.exposure.photoCalib 

200 # PhotoCalib comparison is very approximate because we are basing this 

201 # comparison on just 2-3 stars. 

202 self.assertFloatsAlmostEqual(photo_calib.getCalibrationMean(), self.photo_calib, rtol=1e-2) 

203 # Should have calibrated flux/magnitudes in the afw and astropy catalogs 

204 self.assertIn("slot_PsfFlux_flux", result.stars_footprints.schema) 

205 self.assertIn("slot_PsfFlux_mag", result.stars_footprints.schema) 

206 self.assertEqual(result.stars["slot_PsfFlux_flux"].unit, u.nJy) 

207 self.assertEqual(result.stars["slot_PsfFlux_mag"].unit, u.ABmag) 

208 

209 # Should have detected all S/N >= 10 sources plus 2 sky sources, whether 1 or 2 snaps. 

210 self.assertEqual(len(result.stars), 6) 

211 # Did the psf flags get propagated from the psf_stars catalog? 

212 self.assertEqual(result.stars["calib_psf_used"].sum(), 3) 

213 

214 # Check that all necessary fields are in the output. 

215 lsst.pipe.base.testUtils.assertValidOutput(calibrate, result) 

216 

217 # Check metadata. 

218 key = "LSST CALIB ILLUMCORR APPLIED" 

219 self.assertIn(key, result.exposure.metadata) 

220 self.assertEqual(result.exposure.metadata[key], False) 

221 

222 # Check that the psf_stars cross match worked correctly. 

223 matches = esutil.numpy_util.match(result.psf_stars["id"], result.stars["psf_id"]) 

224 self.assertFloatsAlmostEqual(result.psf_stars["slot_Centroid_x"][matches[0]], 

225 result.stars["slot_Centroid_x"][matches[1]], atol=3e-4) 

226 if "astrometry_matches" in self.config.optional_outputs: 

227 matches = esutil.numpy_util.match(result.astrometry_matches["src_id"], 

228 result.photometry_matches["src_psf_id"]) 

229 self.assertFloatsAlmostEqual(result.astrometry_matches["src_slot_Centroid_x"][matches[0]], 

230 result.photometry_matches["src_slot_Centroid_x"][matches[1]], 

231 atol=3e-4) 

232 

233 def test_run(self): 

234 """Test that run() returns reasonable values to be butler put. 

235 """ 

236 calibrate = CalibrateImageTask(config=self.config) 

237 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

238 calibrate.photometry.match.setRefObjLoader(self.ref_loader) 

239 result = calibrate.run(exposures=self.exposure) 

240 

241 self._check_run(calibrate, result) 

242 

243 def test_run_adaptive_threshold_deteection(self): 

244 """Test that run() runs with adaptive threshold detection turned on. 

245 """ 

246 config = copy.copy(self.config) 

247 # Set the adaptive threshold detection, config values... 

248 config.do_adaptive_threshold_detection = True 

249 config.psf_adaptive_threshold_detection.minFootprint = 4 

250 config.psf_adaptive_threshold_detection.minIsolated = 4 

251 config.psf_adaptive_threshold_detection.sufficientIsolated = 4 

252 config.psf_detection.reEstimateBackground = False 

253 config.star_detection.reEstimateBackground = False 

254 

255 calibrate = CalibrateImageTask(config=config) 

256 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

257 calibrate.photometry.match.setRefObjLoader(self.ref_loader) 

258 with self.assertLogs("lsst.calibrateImage", level="INFO") as cm: 

259 result = calibrate.run(exposures=self.exposure) 

260 subString = "Using adaptive threshold detection " 

261 self.assertTrue(any(subString in s for s in cm.output)) 

262 

263 self._check_run(calibrate, result, expect_n_background=2) 

264 

265 def test_run_downsample(self): 

266 """Test that run() runs with downsample. 

267 """ 

268 config = copy.copy(self.config) 

269 config.do_downsample_footprints = True 

270 config.downsample_max_footprints = 5 

271 

272 calibrate = CalibrateImageTask(config=config) 

273 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

274 calibrate.photometry.match.setRefObjLoader(self.ref_loader) 

275 with self.assertLogs("lsst.calibrateImage", level="INFO") as cm: 

276 result = calibrate.run(exposures=self.exposure) 

277 self.assertIn( 

278 "INFO:lsst.calibrateImage:Downsampling from 8 to 7 non-sky-source footprints.", 

279 cm.output, 

280 ) 

281 

282 self._check_run(calibrate, result) 

283 

284 def test_run_2_snaps(self): 

285 """Test that run() returns reasonable values to be butler put, when 

286 passed two exposures to combine as snaps. 

287 """ 

288 calibrate = CalibrateImageTask(config=self.config) 

289 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

290 calibrate.photometry.match.setRefObjLoader(self.ref_loader) 

291 # Halve the flux in each exposure to get the expected visit sum. 

292 self.exposure.image /= 2 

293 self.exposure.variance /= 2 

294 result = calibrate.run(exposures=[self.exposure, self.exposure]) 

295 

296 self._check_run(calibrate, result) 

297 

298 def test_run_no_optionals(self): 

299 """Test that disabling optional outputs removes them from the output 

300 struct, as appropriate. 

301 """ 

302 self.config.optional_outputs = [] 

303 calibrate = CalibrateImageTask(config=self.config) 

304 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

305 calibrate.photometry.match.setRefObjLoader(self.ref_loader) 

306 result = calibrate.run(exposures=self.exposure) 

307 

308 self._check_run(calibrate, result) 

309 # These are the only optional outputs that require extra computation, 

310 # the others are included in the output struct regardless. 

311 self.assertNotIn("astrometry_matches", result.getDict()) 

312 self.assertNotIn("photometry_matches", result.getDict()) 

313 

314 def test_run_no_calibrate_pixels(self): 

315 """Test that run() returns reasonable values when 

316 do_calibrate_pixels=False. 

317 """ 

318 self.config.do_calibrate_pixels = False 

319 calibrate = CalibrateImageTask(config=self.config) 

320 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

321 calibrate.photometry.match.setRefObjLoader(self.ref_loader) 

322 result = calibrate.run(exposures=self.exposure) 

323 

324 self._check_run(calibrate, result, expect_calibrated_pixels=False) 

325 

326 def test_run_no_astrom_errors(self): 

327 """Test that run() returns reasonable values when 

328 do_calibrate_pixels=False. 

329 """ 

330 self.config.do_include_astrometric_errors = False 

331 calibrate = CalibrateImageTask(config=self.config) 

332 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

333 calibrate.photometry.match.setRefObjLoader(self.ref_loader) 

334 result = calibrate.run(exposures=self.exposure) 

335 

336 self._check_run(calibrate, result) 

337 

338 def test_compute_psf(self): 

339 """Test that our brightest sources are found by _compute_psf(), 

340 that a PSF is assigned to the expopsure. 

341 """ 

342 calibrate = CalibrateImageTask(config=self.config) 

343 psf_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) 

344 

345 # Catalog ids should be very large from this id generator. 

346 self.assertTrue(all(psf_stars['id'] > 1000000000)) 

347 

348 # Background should have 3 elements: initial subtraction, and two from 

349 # re-estimation during the two detection passes. 

350 self.assertEqual(len(background), 3) 

351 

352 # Only the point-sources with S/N > 50 should be in this output. 

353 self.assertEqual(psf_stars["calib_psf_used"].sum(), 3) 

354 # Sort in order of brightness, to easily compare with expected positions. 

355 psf_stars.sort(psf_stars.getPsfFluxSlot().getMeasKey()) 

356 for record, flux, center in zip(psf_stars[::-1], self.fluxes, self.centroids[self.fluxes > 50]): 

357 self.assertFloatsAlmostEqual(record.getX(), center[0], rtol=0.01) 

358 self.assertFloatsAlmostEqual(record.getY(), center[1], rtol=0.01) 

359 # PsfFlux should match the values inserted. 

360 self.assertFloatsAlmostEqual(record["slot_PsfFlux_instFlux"], flux, rtol=0.01) 

361 

362 # TODO: While debugging DM-32701, we're using PCA instead of psfex. 

363 # Check that we got a useable PSF. 

364 # self.assertIsInstance(self.exposure.psf, lsst.meas.extensions.psfex.PsfexPsf) 

365 self.assertIsInstance(self.exposure.psf, lsst.meas.algorithms.PcaPsf) 

366 # TestDataset sources have PSF radius=2 pixels. 

367 radius = self.exposure.psf.computeShape(self.exposure.psf.getAveragePosition()).getDeterminantRadius() 

368 self.assertFloatsAlmostEqual(radius, 2.0, rtol=1e-2) 

369 

370 # To look at images for debugging (`setup display_ds9` and run ds9): 

371 # import lsst.afw.display 

372 # display = lsst.afw.display.getDisplay() 

373 # display.mtv(self.exposure) 

374 

375 def test_compute_psf_bad_centroids(self): 

376 """Test that we raise an appropriate error if all measured centroids 

377 are bad in compute_psf. 

378 The root cause of this is likely an interaction between the PSF fit 

379 and cosmic ray repair. An ugly true PSF can result in the centers of 

380 sources being masked and repaired as CRs, thus removing the center of 

381 the sources and causing them all to be flagged. 

382 """ 

383 # make the sources have a (not too deep) hole in the middle 

384 for point in self.centroids: 

385 for i in (-1, 0, 1): 

386 for j in (-1, 0, 1): 

387 self.exposure.image[point[0]+i, point[1]+j] -= \ 

388 (self.exposure.image[point[0]+i, point[1]+j] - self.background_level)/1.1 

389 # and the extended central source 

390 point = self.extended_source 

391 self.exposure.image[point[0], point[1]] -= \ 

392 (self.exposure.image[point[0], point[1]] - self.background_level) 

393 

394 # add a diagonal gradient to each source 

395 size = 5 

396 for point, flux in zip(self.centroids, self.fluxes): 

397 box = lsst.geom.Box2I.makeCenteredBox(lsst.geom.Point2D(point[0], point[1]), 

398 lsst.geom.Extent2I(size*2, size*2)) 

399 X, Y = np.ogrid[point[0] - size:point[0] + size, point[1] - size:point[1] + size] 

400 distance = np.sqrt((X - point[0])**2 + (Y - point[1])**2) 

401 gradient = (X - point[0] + size)*10 + (Y - point[1] + size)*10 

402 gradient[distance > size] = 0 

403 cutout = self.exposure.image.subset(box) 

404 cutout.array += gradient / flux 

405 

406 calibrate = CalibrateImageTask(config=self.config) 

407 with self.assertRaisesRegex(AllCentroidsFlaggedError, r"source centroids \(out of 4\) flagged"): 

408 calibrate.run(exposures=[self.exposure], id_generator=self.id_generator) 

409 

410 def test_measure_aperture_correction(self): 

411 """Test that _measure_aperture_correction() assigns an ApCorrMap to the 

412 exposure. 

413 """ 

414 calibrate = CalibrateImageTask(config=self.config) 

415 psf_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) 

416 

417 # First check that the exposure doesn't have an ApCorrMap. 

418 self.assertIsNone(self.exposure.apCorrMap) 

419 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

420 self.assertIsInstance(self.exposure.apCorrMap, afwImage.ApCorrMap) 

421 # We know that there are 2 fields from the normalization, plus more from 

422 # other configured plugins. 

423 self.assertGreater(len(self.exposure.apCorrMap), 2) 

424 

425 def test_find_stars(self): 

426 """Test that _find_stars() correctly identifies the S/N>10 stars 

427 in the image and returns them in the output catalog. 

428 """ 

429 calibrate = CalibrateImageTask(config=self.config) 

430 psf_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) 

431 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

432 

433 stars = calibrate._find_stars(self.exposure, background, self.id_generator) 

434 

435 # Catalog ids should be very large from this id generator. 

436 self.assertTrue(all(stars['id'] > 1000000000)) 

437 

438 # Background should have 4 elements: 3 from compute_psf and one from 

439 # re-estimation during source detection. 

440 self.assertEqual(len(background), 4) 

441 

442 # Only 5 psf-like sources with S/N>10 should be in the output catalog, 

443 # plus two sky sources. 

444 self.assertEqual(len(stars), 6) 

445 self.assertTrue(stars.isContiguous()) 

446 # Sort in order of brightness, to easily compare with expected positions. 

447 stars.sort(stars.getPsfFluxSlot().getMeasKey()) 

448 for record, flux, center in zip(stars[::-1], self.fluxes, self.centroids[self.fluxes > 50]): 

449 self.assertFloatsAlmostEqual(record.getX(), center[0], rtol=0.01) 

450 self.assertFloatsAlmostEqual(record.getY(), center[1], rtol=0.01) 

451 self.assertFloatsAlmostEqual(record["slot_PsfFlux_instFlux"], flux, rtol=0.01) 

452 

453 def test_astrometry(self): 

454 """Test that the fitted WCS gives good catalog coordinates. 

455 """ 

456 calibrate = CalibrateImageTask(config=self.config) 

457 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

458 psf_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) 

459 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

460 stars = calibrate._find_stars(self.exposure, background, self.id_generator) 

461 

462 calibrate._fit_astrometry(self.exposure, stars) 

463 

464 # Check that we got reliable matches with the truth coordinates. 

465 sky = stars["sky_source"] 

466 fitted = SkyCoord(stars[~sky]['coord_ra'], stars[~sky]['coord_dec'], unit="radian") 

467 truth = SkyCoord(self.truth_cat['coord_ra'], self.truth_cat['coord_dec'], unit="radian") 

468 idx, d2d, _ = fitted.match_to_catalog_sky(truth) 

469 np.testing.assert_array_less(d2d.to_value(u.milliarcsecond), 35.0) 

470 

471 def test_photometry(self): 

472 """Test that the fitted photoCalib matches the one we generated, 

473 and that the exposure is calibrated. 

474 """ 

475 calibrate = CalibrateImageTask(config=self.config) 

476 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

477 calibrate.photometry.match.setRefObjLoader(self.ref_loader) 

478 psf_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) 

479 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

480 stars = calibrate._find_stars(self.exposure, background, self.id_generator) 

481 calibrate._fit_astrometry(self.exposure, stars) 

482 

483 stars, matches, meta, photoCalib = calibrate._fit_photometry(self.exposure, stars) 

484 calibrate._apply_photometry(self.exposure, background) 

485 

486 # NOTE: With this test data, PhotoCalTask returns calibrationErr==0, 

487 # so we can't check that the photoCal error has been set. 

488 self.assertFloatsAlmostEqual(photoCalib.getCalibrationMean(), self.photo_calib, rtol=1e-2) 

489 # The exposure should be calibrated by the applied photoCalib, 

490 # and the background should be calibrated to match. 

491 uncalibrated = self.exposure.image.clone() 

492 uncalibrated += background.getImage() 

493 uncalibrated /= self.photo_calib 

494 self.assertFloatsAlmostEqual(uncalibrated.array, self.truth_exposure.image.array, rtol=1e-2) 

495 # PhotoCalib on the exposure must be identically 1. 

496 self.assertEqual(self.exposure.photoCalib.getCalibrationMean(), 1.0) 

497 

498 # Check that we got reliable magnitudes and fluxes vs. truth, ignoring 

499 # sky sources. 

500 sky = stars["sky_source"] 

501 fitted = SkyCoord(stars[~sky]['coord_ra'], stars[~sky]['coord_dec'], unit="radian") 

502 truth = SkyCoord(self.truth_cat['coord_ra'], self.truth_cat['coord_dec'], unit="radian") 

503 idx, _, _ = fitted.match_to_catalog_sky(truth) 

504 # Because the input variance image does not include contributions from 

505 # the sources, we can't use fluxErr as a bound on the measurement 

506 # quality here. 

507 self.assertFloatsAlmostEqual(stars[~sky]['slot_PsfFlux_flux'], 

508 self.truth_cat['truth_flux'][idx], 

509 rtol=0.1) 

510 self.assertFloatsAlmostEqual(stars[~sky]['slot_PsfFlux_mag'], 

511 self.truth_cat['truth_mag'][idx], 

512 rtol=0.01) 

513 

514 def test_match_psf_stars(self): 

515 """Test that _match_psf_stars() flags the correct stars as psf stars 

516 and candidates. 

517 """ 

518 calibrate = CalibrateImageTask(config=self.config) 

519 psf_stars, background, candidates, _ = calibrate._compute_psf(self.exposure, self.id_generator) 

520 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

521 stars = calibrate._find_stars(self.exposure, background, self.id_generator) 

522 

523 # There should be no psf-related flags set at first. 

524 self.assertEqual(stars["calib_psf_candidate"].sum(), 0) 

525 self.assertEqual(stars["calib_psf_used"].sum(), 0) 

526 self.assertEqual(stars["calib_psf_reserved"].sum(), 0) 

527 

528 # Reorder stars to be out of order with psf_stars (putting the sky 

529 # sources in front); this tests that I get the indexing right. 

530 stars.sort(stars.getCentroidSlot().getMeasKey().getX()) 

531 stars = stars.copy(deep=True) 

532 # Re-number the ids: the matcher requires sorted ids: this is always 

533 # true in the code itself, but we've permuted them by sorting on 

534 # flux. We don't care what the actual ids themselves are here. 

535 stars["id"] = np.arange(len(stars)) 

536 

537 calibrate._match_psf_stars(psf_stars, stars) 

538 

539 # Check that the three brightest stars have the psf flags transfered 

540 # from the psf_stars catalog by sorting in order of brightness. 

541 stars.sort(stars.getPsfFluxSlot().getMeasKey()) 

542 # sort() above leaves the catalog non-contiguous. 

543 stars = stars.copy(deep=True) 

544 np.testing.assert_array_equal(stars["calib_psf_candidate"], 

545 [False, False, False, True, True, True]) 

546 np.testing.assert_array_equal(stars["calib_psf_used"], 

547 [False, False, False, True, True, True]) 

548 # Too few sources to reserve any in these tests. 

549 self.assertEqual(stars["calib_psf_reserved"].sum(), 0) 

550 

551 def test_match_psf_stars_no_matches(self): 

552 """Check that _match_psf_stars handles the case of no cross-matches. 

553 """ 

554 calibrate = CalibrateImageTask(config=self.config) 

555 # Make two catalogs that cannot have matches. 

556 stars = self.truth_cat[2:].copy(deep=True) 

557 psf_stars = self.truth_cat[:2].copy(deep=True) 

558 

559 with self.assertRaisesRegex(NoPsfStarsToStarsMatchError, 

560 "No psf stars out of 2 matched 5 calib stars") as cm: 

561 calibrate._match_psf_stars(psf_stars, stars) 

562 self.assertEqual(cm.exception.metadata["n_psf_stars"], 2) 

563 self.assertEqual(cm.exception.metadata["n_stars"], 5) 

564 

565 def test_calibrate_image_illumcorr(self): 

566 """Test running through with an illumination correction.""" 

567 config = copy.copy(self.config) 

568 config.do_illumination_correction = True 

569 config.psf_subtract_background.doApplyFlatBackgroundRatio = True 

570 config.psf_detection.doApplyFlatBackgroundRatio = True 

571 config.star_background.doApplyFlatBackgroundRatio = True 

572 config.star_detection.doApplyFlatBackgroundRatio = True 

573 

574 calibrate = CalibrateImageTask(config=config) 

575 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

576 calibrate.photometry.match.setRefObjLoader(self.ref_loader) 

577 

578 # Assume that the exposure has been flattened by a flat-flat. 

579 background_flat = self.exposure.clone() 

580 background_flat.image.array[:, :] = 1.0 

581 background_flat.mask.array[:, :] = 0 

582 background_flat.variance.array[:, :] = 0.0 

583 

584 # And create an illumination correction of 1.1. 

585 illum_corr_value = 1.1 

586 illumination_correction = self.exposure.clone() 

587 illumination_correction.image.array[:, :] = illum_corr_value 

588 illumination_correction.mask.array[:, :] = 0 

589 illumination_correction.variance.array[:, :] = 0.0 

590 

591 result = calibrate.run( 

592 exposures=self.exposure, 

593 id_generator=self.id_generator, 

594 background_flat=background_flat, 

595 illumination_correction=illumination_correction, 

596 ) 

597 

598 # We divide the image by the illumination correction, but the reference 

599 # sources stay the same, so the applied photocalib will increase by 

600 # the illum_corr_value. 

601 # Tolerance is the same as the direct test with no illumination 

602 # correction. 

603 self.assertFloatsAlmostEqual( 

604 result.applied_photo_calib.getCalibrationMean(), 

605 self.photo_calib * illum_corr_value, 

606 rtol=1e-2, 

607 ) 

608 

609 self.assertEqual(len(result.background), 4) 

610 self.assertFloatsAlmostEqual( 

611 np.median(result.background.getImage().array), 

612 result.applied_photo_calib.getCalibrationMean() * self.background_level, 

613 rtol=1e-3, 

614 ) 

615 

616 # Check metadata. 

617 key = "LSST CALIB ILLUMCORR APPLIED" 

618 self.assertIn(key, result.exposure.metadata) 

619 self.assertEqual(result.exposure.metadata[key], True) 

620 

621 def test_run_with_diffraction_spike_mask(self): 

622 """Test that the diffraction spike mask subtask runs. 

623 """ 

624 config = self.config 

625 config.doMaskDiffractionSpikes = True 

626 config.diffractionSpikeMask.magnitudeThreshold = 19 

627 # Define a fake SATURATED mask plane for use by the diffraction spike 

628 # task, so that it does not affect the rest of calibrate 

629 config.diffractionSpikeMask.saturatedMaskPlane = "FAKESATURATED" 

630 calibrate = CalibrateImageTask(config=config) 

631 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

632 calibrate.photometry.match.setRefObjLoader(self.ref_loader) 

633 calibrate.diffractionSpikeMask.setRefObjLoader(self.ref_loader) 

634 

635 exposure = self.exposure.clone() 

636 

637 exposure.info.setVisitInfo(makeTestVisitInfo()) 

638 exposure.mask.addMaskPlane(config.diffractionSpikeMask.saturatedMaskPlane) 

639 

640 # Set the saturated mask plane in half of the image 

641 saturatedMaskBit = exposure.mask.getPlaneBitMask(config.diffractionSpikeMask.saturatedMaskPlane) 

642 bbox = exposure.getBBox() 

643 bbox.grow(-lsst.geom.Extent2I(0, bbox.height//4)) 

644 exposure[bbox].mask.array |= saturatedMaskBit 

645 

646 result = calibrate.run(exposures=exposure) 

647 

648 self._check_run(calibrate, result) 

649 

650 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": ""}) 

651 def test_fail_on_sattle_miconfiguration(self): 

652 """Test for failure if sattle is requested without appropriate configurations. 

653 """ 

654 self.config.run_sattle = True 

655 with self.assertRaises(pexConfig.FieldValidationError): 

656 CalibrateImageTask(config=self.config) 

657 

658 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"}) 

659 def test_continue_on_sattle_failure(self): 

660 """Processing should continue when sattle returns status codes other than 200. 

661 """ 

662 response = MockResponse({}, 500, "internal sattle error") 

663 

664 self.config.run_sattle = True 

665 calibrate = CalibrateImageTask(config=self.config) 

666 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

667 calibrate.photometry.match.setRefObjLoader(self.ref_loader) 

668 with mock.patch('requests.put', return_value=response) as mock_put: 

669 calibrate.run(exposures=self.exposure) 

670 mock_put.assert_called_once() 

671 

672 @mock.patch.dict(os.environ, {"SATTLE_URI_BASE": "fake_host:1234"}) 

673 def test_sattle(self): 

674 """Test for successful completion when sattle call returns successfully. 

675 """ 

676 response = MockResponse({}, 200, "success") 

677 

678 self.config.run_sattle = True 

679 calibrate = CalibrateImageTask(config=self.config) 

680 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

681 calibrate.photometry.match.setRefObjLoader(self.ref_loader)