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 as pipeBase 

47import lsst.pipe.base.testUtils 

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

49 NoPsfStarsToStarsMatchError, AllCentroidsFlaggedError 

50import lsst.pex.config as pexConfig 

51import lsst.utils.tests 

52 

53from utils import makeTestVisitInfo 

54 

55 

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

57 

58 def setUp(self): 

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

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

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

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

63 self.photo_calib = 12.3 

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

65 # sqrt of area of a normalized 2d gaussian 

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

67 noise = 10.0 # stddev of noise per pixel 

68 # Sources ordered from faintest to brightest. 

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

70 12*noise*psf_scale, 

71 45*noise*psf_scale, 

72 150*noise*psf_scale, 

73 400*noise*psf_scale, 

74 1000*noise*psf_scale)) 

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

76 (40, 70), 

77 (100, 160), 

78 (50, 120), 

79 (92, 35), 

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

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

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

83 

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

85 # in any of the output catalogs. 

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

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

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

89 

90 schema = dataset.makeMinimalSchema() 

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

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

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

94 self.background_level = 500.0 

95 self.truth_exposure.image += self.background_level 

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

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

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

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

100 metadata.set("REFCAT_FORMAT_VERSION", 1) 

101 self.truth_cat.setMetadata(metadata) 

102 

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

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

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

106 

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

108 # Post-ISR images only contain: initial WCS, VisitInfo, filter 

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

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

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

112 self.exposure.info.id = 12345 

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

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

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

116 

117 # Set up a basic results struct to hold exposure attribute data 

118 self.attributes = pipeBase.Struct() 

119 self.attributes.exposure = self.exposure 

120 self.attributes.background = None 

121 self.attributes.background_to_photometric_ratio = None 

122 

123 # Test-specific configuration: 

124 self.config = CalibrateImageTask.ConfigClass() 

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

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

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

128 # Only insert 2 sky sources, for simplicity. 

129 self.config.star_sky_sources.nSources = 2 

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

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

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

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

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

135 self.config.run_sattle = False 

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

137 self.config.do_adaptive_threshold_detection = False 

138 self.config.psf_detection.reEstimateBackground = True 

139 self.config.star_detection.reEstimateBackground = True 

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

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

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

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

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

145 self.config.astrometry.magnitudeOutlierRejectionNSigma = 9.0 

146 

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

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

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

150 # the centroids and fluxes measured on the image, which might cause 

151 # tests to fail. 

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

153 instrument="I", 

154 visit=self.truth_exposure.visitInfo.id, 

155 detector=12, 

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

157 ) 

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

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

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

161 self.config.id_generator.n_releases = 8 

162 self.config.id_generator.release_id = 2 

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

164 

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

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

167 

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

169 expect_n_background: int = 4, expect_n_background_equal_or_greater_than: int = -1): 

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

171 

172 Parameters 

173 ---------- 

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

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

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

177 Result of calling calibrate.run(). 

178 expect_calibrated_pixels : `bool`, optional 

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

180 """ 

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

182 # re-estimation during source detection. 

183 if expect_n_background is not None: 

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

185 if expect_n_background_equal_or_greater_than is not None: 

186 self.assertTrue(len(result.background) >= expect_n_background_equal_or_greater_than) 

187 

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

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

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

191 

192 # Check that the summary statistics are reasonable. 

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

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

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

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

197 

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

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

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

201 

202 if expect_calibrated_pixels: 

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

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

205 photo_calib = result.applied_photo_calib 

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

207 else: 

208 self.assertIsNone(result.applied_photo_calib) 

209 photo_calib = result.exposure.photoCalib 

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

211 # comparison on just 2-3 stars. 

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

213 # Should have calibrated flux/mags in the afw and astropy catalogs 

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

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

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

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

218 

219 # Should have detected all S/N >= 10 sources plus 2 sky sources, 

220 # whether 1 or 2 snaps. 

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

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

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

224 

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

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

227 

228 # Check metadata. 

229 key = "LSST CALIB ILLUMCORR APPLIED" 

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

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

232 

233 # Check that the psf_stars cross match worked correctly. 

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

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

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

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

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

239 result.photometry_matches["src_psf_id"]) 

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

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

242 atol=3e-4) 

243 

244 def test_run(self): 

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

246 """ 

247 calibrate = CalibrateImageTask(config=self.config) 

248 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

251 

252 self._check_run(calibrate, result) 

253 

254 def test_run_adaptive_threshold_detection(self): 

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

256 """ 

257 config = copy.copy(self.config) 

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

259 config.do_adaptive_threshold_detection = True 

260 config.psf_adaptive_threshold_detection.minFootprint = 4 

261 config.psf_adaptive_threshold_detection.minIsolated = 4 

262 config.psf_adaptive_threshold_detection.sufficientIsolated = 4 

263 config.psf_detection.reEstimateBackground = False 

264 config.star_detection.reEstimateBackground = False 

265 

266 calibrate = CalibrateImageTask(config=config) 

267 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

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

271 subString = "Using adaptive threshold detection " 

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

273 

274 # Number of backgrounds in list is only guaranteed to have at least 2 

275 # entries in the adaptive threshold code path. 

276 self._check_run(calibrate, result, expect_n_background=None, 

277 expect_n_background_equal_or_greater_than=2) 

278 

279 def test_run_downsample(self): 

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

281 """ 

282 config = copy.copy(self.config) 

283 config.do_downsample_footprints = True 

284 config.downsample_max_footprints = 5 

285 

286 calibrate = CalibrateImageTask(config=config) 

287 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

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

291 self.assertIn( 

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

293 cm.output, 

294 ) 

295 

296 self._check_run(calibrate, result) 

297 

298 def test_run_2_snaps(self): 

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

300 passed two exposures to combine as snaps. 

301 """ 

302 calibrate = CalibrateImageTask(config=self.config) 

303 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

306 self.exposure.image /= 2 

307 self.exposure.variance /= 2 

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

309 

310 self._check_run(calibrate, result) 

311 

312 def test_run_no_optionals(self): 

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

314 struct, as appropriate. 

315 """ 

316 self.config.optional_outputs = [] 

317 calibrate = CalibrateImageTask(config=self.config) 

318 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

321 

322 self._check_run(calibrate, result) 

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

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

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

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

327 

328 def test_run_no_calibrate_pixels(self): 

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

330 do_calibrate_pixels=False. 

331 """ 

332 self.config.do_calibrate_pixels = False 

333 calibrate = CalibrateImageTask(config=self.config) 

334 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

337 

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

339 

340 def test_run_no_astrom_errors(self): 

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

342 do_calibrate_pixels=False. 

343 """ 

344 self.config.do_include_astrometric_errors = False 

345 calibrate = CalibrateImageTask(config=self.config) 

346 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

349 

350 self._check_run(calibrate, result) 

351 

352 def test_compute_psf(self): 

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

354 that a PSF is assigned to the exposure. 

355 """ 

356 calibrate = CalibrateImageTask(config=self.config) 

357 psf_stars, candidates, _ = calibrate._compute_psf(self.attributes, self.id_generator) 

358 

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

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

361 

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

363 # re-estimation during the two detection passes. 

364 self.assertEqual(len(self.attributes.background), 3) 

365 

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

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

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

369 psf_stars.sort(psf_stars.getPsfFluxSlot().getMeasKey()) 

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

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

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

373 # PsfFlux should match the values inserted. 

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

375 

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

377 # Check that we got a useable PSF. 

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

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

380 # TestDataset sources have PSF radius=2 pixels. 

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

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

383 

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

385 # import lsst.afw.display 

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

387 # display.mtv(self.exposure) 

388 

389 def test_compute_psf_bad_centroids(self): 

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

391 are bad in compute_psf. 

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

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

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

395 the sources and causing them all to be flagged. 

396 """ 

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

398 for point in self.centroids: 

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

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

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

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

403 # and the extended central source 

404 point = self.extended_source 

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

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

407 

408 # add a diagonal gradient to each source 

409 size = 5 

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

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

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

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

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

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

416 gradient[distance > size] = 0 

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

418 cutout.array += gradient / flux 

419 

420 calibrate = CalibrateImageTask(config=self.config) 

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

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

423 

424 def test_measure_aperture_correction(self): 

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

426 exposure. 

427 """ 

428 calibrate = CalibrateImageTask(config=self.config) 

429 psf_stars, candidates, _ = calibrate._compute_psf(self.attributes, self.id_generator) 

430 

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

432 self.assertIsNone(self.exposure.apCorrMap) 

433 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

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

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

436 # from other configured plugins. 

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

438 

439 def test_find_stars(self): 

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

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

442 """ 

443 calibrate = CalibrateImageTask(config=self.config) 

444 psf_stars, candidates, _ = calibrate._compute_psf(self.attributes, self.id_generator) 

445 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

446 

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

448 

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

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

451 

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

453 # re-estimation during source detection. 

454 self.assertEqual(len(self.attributes.background), 4) 

455 

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

457 # plus two sky sources. 

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

459 self.assertTrue(stars.isContiguous()) 

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

461 stars.sort(stars.getPsfFluxSlot().getMeasKey()) 

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

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

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

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

466 

467 def test_astrometry(self): 

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

469 """ 

470 calibrate = CalibrateImageTask(config=self.config) 

471 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

472 psf_stars, candidates, _ = calibrate._compute_psf(self.attributes, self.id_generator) 

473 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

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

475 

476 calibrate._fit_astrometry(self.exposure, stars) 

477 

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

479 sky = stars["sky_source"] 

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

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

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

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

484 

485 def test_photometry(self): 

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

487 and that the exposure is calibrated. 

488 """ 

489 calibrate = CalibrateImageTask(config=self.config) 

490 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

492 psf_stars, candidates, _ = calibrate._compute_psf(self.attributes, self.id_generator) 

493 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

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

495 calibrate._fit_astrometry(self.exposure, stars) 

496 

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

498 calibrate._apply_photometry(self.exposure, self.attributes.background) 

499 

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

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

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

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

504 # and the background should be calibrated to match. 

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

506 uncalibrated += self.attributes.background.getImage() 

507 uncalibrated /= self.photo_calib 

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

509 # PhotoCalib on the exposure must be identically 1. 

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

511 

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

513 # sky sources. 

514 sky = stars["sky_source"] 

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

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

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

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

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

520 # quality here. 

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

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

523 rtol=0.1) 

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

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

526 rtol=0.01) 

527 

528 def test_match_psf_stars(self): 

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

530 and candidates. 

531 """ 

532 calibrate = CalibrateImageTask(config=self.config) 

533 psf_stars, candidates, _ = calibrate._compute_psf(self.attributes, self.id_generator) 

534 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

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

536 

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

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

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

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

541 

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

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

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

545 stars = stars.copy(deep=True) 

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

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

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

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

550 

551 calibrate._match_psf_stars(psf_stars, stars) 

552 

553 # Check that the three brightest stars have the psf flags transferred 

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

555 stars.sort(stars.getPsfFluxSlot().getMeasKey()) 

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

557 stars = stars.copy(deep=True) 

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

559 [False, False, False, True, True, True]) 

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

561 [False, False, False, True, True, True]) 

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

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

564 

565 def test_match_psf_stars_no_matches(self): 

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

567 """ 

568 calibrate = CalibrateImageTask(config=self.config) 

569 # Make two catalogs that cannot have matches. 

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

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

572 

573 with self.assertRaisesRegex(NoPsfStarsToStarsMatchError, 

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

575 calibrate._match_psf_stars(psf_stars, stars) 

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

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

578 

579 def test_calibrate_image_illumcorr(self): 

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

581 config = copy.copy(self.config) 

582 config.do_illumination_correction = True 

583 config.psf_subtract_background.doApplyFlatBackgroundRatio = True 

584 config.psf_detection.doApplyFlatBackgroundRatio = True 

585 config.star_background.doApplyFlatBackgroundRatio = True 

586 config.star_detection.doApplyFlatBackgroundRatio = True 

587 

588 calibrate = CalibrateImageTask(config=config) 

589 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

591 

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

593 background_flat = self.exposure.clone() 

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

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

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

597 

598 # And create an illumination correction of 1.1. 

599 illum_corr_value = 1.1 

600 illumination_correction = self.exposure.clone() 

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

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

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

604 

605 result = calibrate.run( 

606 exposures=self.exposure, 

607 id_generator=self.id_generator, 

608 background_flat=background_flat, 

609 illumination_correction=illumination_correction, 

610 ) 

611 

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

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

614 # the illum_corr_value. 

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

616 # correction. 

617 self.assertFloatsAlmostEqual( 

618 result.applied_photo_calib.getCalibrationMean(), 

619 self.photo_calib * illum_corr_value, 

620 rtol=1e-2, 

621 ) 

622 

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

624 self.assertFloatsAlmostEqual( 

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

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

627 rtol=1e-3, 

628 ) 

629 

630 # Check metadata. 

631 key = "LSST CALIB ILLUMCORR APPLIED" 

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

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

634 

635 def test_run_with_diffraction_spike_mask(self): 

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

637 """ 

638 config = self.config 

639 config.doMaskDiffractionSpikes = True 

640 config.diffractionSpikeMask.magnitudeThreshold = 19 

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

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

643 config.diffractionSpikeMask.saturatedMaskPlane = "FAKESATURATED" 

644 calibrate = CalibrateImageTask(config=config) 

645 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

647 calibrate.diffractionSpikeMask.setRefObjLoader(self.ref_loader) 

648 

649 exposure = self.exposure.clone() 

650 

651 exposure.info.setVisitInfo(makeTestVisitInfo()) 

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

653 

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

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

656 bbox = exposure.getBBox() 

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

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

659 

660 result = calibrate.run(exposures=exposure) 

661 

662 self._check_run(calibrate, result) 

663 

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

665 def test_fail_on_sattle_misconfiguration(self): 

666 """Test for failure if sattle is requested without appropriate 

667 configurations. 

668 """ 

669 self.config.run_sattle = True 

670 with self.assertRaises(pexConfig.FieldValidationError): 

671 CalibrateImageTask(config=self.config) 

672 

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

674 def test_continue_on_sattle_failure(self): 

675 """Processing should continue when sattle returns status codes other 

676 than 200. 

677 """ 

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

679 

680 self.config.run_sattle = True 

681 calibrate = CalibrateImageTask(config=self.config) 

682 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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