Coverage for tests / test_calibrateImage.py: 13%

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.compute_shapelet_decomposition.min_n_stars = 4 

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

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

129 # Only insert 2 sky sources, for simplicity. 

130 self.config.star_sky_sources.nSources = 2 

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

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

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

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

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

136 self.config.run_sattle = False 

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

138 self.config.do_adaptive_threshold_detection = False 

139 self.config.psf_detection.reEstimateBackground = True 

140 self.config.star_detection.reEstimateBackground = True 

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

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

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

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

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

146 self.config.astrometry.magnitudeOutlierRejectionNSigma = 9.0 

147 

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

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

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

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

152 # tests to fail. 

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

154 instrument="I", 

155 visit=self.truth_exposure.visitInfo.id, 

156 detector=12, 

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

158 ) 

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

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

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

162 self.config.id_generator.n_releases = 8 

163 self.config.id_generator.release_id = 2 

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

165 

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

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

168 

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

170 expect_n_background: int = 4, expect_n_background_equal_or_greater_than: int = -1, 

171 do_shapelet_check: bool = False): 

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

173 

174 Parameters 

175 ---------- 

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

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

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

179 Result of calling calibrate.run(). 

180 expect_calibrated_pixels : `bool`, optional 

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

182 """ 

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

184 # re-estimation during source detection. 

185 if expect_n_background is not None: 

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

187 if expect_n_background_equal_or_greater_than is not None: 

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

189 

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

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

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

193 

194 # Check that the summary statistics are reasonable. 

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

196 # Values below are only correct for the base run test. 

197 if do_shapelet_check: 

198 self.assertFloatsAlmostEqual(summary.shapeletsOnlyIqScore, 8.104235164021247e-06, rtol=1e-4) 

199 self.assertFloatsAlmostEqual(summary.shapeletsIqScore, 6.651368381893731e-05, rtol=1e-4) 

200 shapeletsCoeffs = [ 

201 0.28358882665634155, -0.00037077977926209773, -0.002802249127080195, 

202 3.439385909587145e-05, -0.0012995795522691303, -0.0018000953572589512, 

203 4.4449956901327674e-05, -2.5397295072102522e-05, 3.0203139976223994e-05, 

204 0.000617273102631627, -4.456396106385861e-05, 8.730505214771256e-05, 

205 0.00011167120828758925, 0.00017521437985124066, -1.5413010141236832e-05, 

206 -0.00015419833136474156, -8.099463936604213e-05, -0.00017129145577976686, 

207 -9.799635877346543e-05, 2.7090113467370116e-05, 0.0004025635973896668, 

208 -6.346980774284659e-05, 0.00013586174545049867, -0.00011295261583971951, 

209 0.0003724458365468307, 1.5024524355547979e-05, 0.00012497524707246052, 

210 0.00041057434072064747 

211 ] 

212 for iShapelet in range(len(shapeletsCoeffs)): 

213 self.assertFloatsAlmostEqual( 

214 summary.shapeletsCoeffs[iShapelet], shapeletsCoeffs[iShapelet], rtol=2e-3) 

215 self.assertEqual(summary.nShapeletsStar, 4) 

216 self.assertFloatsAlmostEqual( 

217 summary.centroidDiffShapeletsVsSlotMedian, 0.10838090573, rtol=1e-3) 

218 

219 self.assertFloatsAlmostEqual(summary.refCatSourceDensity, 39470.0430204, rtol=1e-3) 

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

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

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

223 

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

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

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

227 

228 if expect_calibrated_pixels: 

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

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

231 photo_calib = result.applied_photo_calib 

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

233 else: 

234 self.assertIsNone(result.applied_photo_calib) 

235 photo_calib = result.exposure.photoCalib 

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

237 # comparison on just 2-3 stars. 

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

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

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

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

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

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

244 

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

246 # whether 1 or 2 snaps. 

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

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

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

250 

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

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

253 

254 # Check metadata. 

255 key = "LSST CALIB ILLUMCORR APPLIED" 

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

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

258 key = "REF_CAT_SOURCE_DENSITY" 

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

260 

261 # Check that the psf_stars cross match worked correctly. 

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

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

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

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

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

267 result.photometry_matches["src_psf_id"]) 

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

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

270 atol=3e-4) 

271 

272 def test_run(self): 

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

274 """ 

275 calibrate = CalibrateImageTask(config=self.config) 

276 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

279 

280 self._check_run(calibrate, result, do_shapelet_check=True) 

281 

282 def test_run_adaptive_threshold_detection(self): 

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

284 """ 

285 config = copy.copy(self.config) 

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

287 config.do_adaptive_threshold_detection = True 

288 config.psf_adaptive_threshold_detection.minFootprint = 4 

289 config.psf_adaptive_threshold_detection.minIsolated = 4 

290 config.psf_adaptive_threshold_detection.sufficientIsolated = 4 

291 config.psf_detection.reEstimateBackground = False 

292 config.star_detection.reEstimateBackground = False 

293 

294 calibrate = CalibrateImageTask(config=config) 

295 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

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

299 subString = "Using adaptive threshold detection " 

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

301 

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

303 # entries in the adaptive threshold code path. 

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

305 expect_n_background_equal_or_greater_than=2) 

306 

307 def test_run_downsample(self): 

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

309 """ 

310 config = copy.copy(self.config) 

311 config.do_downsample_footprints = True 

312 config.downsample_max_footprints = 5 

313 

314 calibrate = CalibrateImageTask(config=config) 

315 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

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

319 self.assertIn( 

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

321 cm.output, 

322 ) 

323 

324 self._check_run(calibrate, result) 

325 

326 def test_run_2_snaps(self): 

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

328 passed two exposures to combine as snaps. 

329 """ 

330 calibrate = CalibrateImageTask(config=self.config) 

331 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

334 self.exposure.image /= 2 

335 self.exposure.variance /= 2 

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

337 

338 self._check_run(calibrate, result) 

339 

340 def test_run_no_optionals(self): 

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

342 struct, as appropriate. 

343 """ 

344 self.config.optional_outputs = [] 

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 # These are the only optional outputs that require extra computation, 

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

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

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

355 

356 def test_run_no_calibrate_pixels(self): 

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

358 do_calibrate_pixels=False. 

359 """ 

360 self.config.do_calibrate_pixels = False 

361 calibrate = CalibrateImageTask(config=self.config) 

362 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

365 

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

367 

368 def test_run_no_astrom_errors(self): 

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

370 do_calibrate_pixels=False. 

371 """ 

372 self.config.do_include_astrometric_errors = False 

373 calibrate = CalibrateImageTask(config=self.config) 

374 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

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

377 

378 self._check_run(calibrate, result) 

379 

380 def test_compute_psf(self): 

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

382 that a PSF is assigned to the exposure. 

383 """ 

384 calibrate = CalibrateImageTask(config=self.config) 

385 compute_psf_struct = calibrate._compute_psf(self.attributes, self.id_generator) 

386 psf_stars = compute_psf_struct.detections_sources 

387 

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

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

390 

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

392 # re-estimation during the two detection passes. 

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

394 

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

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

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

398 psf_stars.sort(psf_stars.getPsfFluxSlot().getMeasKey()) 

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

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

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

402 # PsfFlux should match the values inserted. 

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

404 

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

406 # Check that we got a useable PSF. 

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

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

409 # TestDataset sources have PSF radius=2 pixels. 

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

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

412 

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

414 # import lsst.afw.display 

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

416 # display.mtv(self.exposure) 

417 

418 def test_compute_psf_bad_centroids(self): 

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

420 are bad in compute_psf. 

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

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

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

424 the sources and causing them all to be flagged. 

425 """ 

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

427 for point in self.centroids: 

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

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

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

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

432 # and the extended central source 

433 point = self.extended_source 

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

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

436 

437 # add a diagonal gradient to each source 

438 size = 5 

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

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

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

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

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

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

445 gradient[distance > size] = 0 

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

447 cutout.array += gradient / flux 

448 

449 calibrate = CalibrateImageTask(config=self.config) 

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

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

452 

453 def test_measure_aperture_correction(self): 

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

455 exposure. 

456 """ 

457 calibrate = CalibrateImageTask(config=self.config) 

458 compute_psf_struct = calibrate._compute_psf(self.attributes, self.id_generator) 

459 psf_stars = compute_psf_struct.detections_sources 

460 

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

462 self.assertIsNone(self.exposure.apCorrMap) 

463 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

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

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

466 # from other configured plugins. 

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

468 

469 def test_find_stars(self): 

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

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

472 """ 

473 calibrate = CalibrateImageTask(config=self.config) 

474 compute_psf_struct = calibrate._compute_psf(self.attributes, self.id_generator) 

475 psf_stars = compute_psf_struct.detections_sources 

476 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

477 

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

479 

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

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

482 

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

484 # re-estimation during source detection. 

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

486 

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

488 # plus two sky sources. 

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

490 self.assertTrue(stars.isContiguous()) 

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

492 stars.sort(stars.getPsfFluxSlot().getMeasKey()) 

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

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

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

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

497 

498 def test_astrometry(self): 

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

500 """ 

501 calibrate = CalibrateImageTask(config=self.config) 

502 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

503 compute_psf_struct = calibrate._compute_psf(self.attributes, self.id_generator) 

504 psf_stars = compute_psf_struct.detections_sources 

505 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

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

507 

508 calibrate._fit_astrometry(self.exposure, stars) 

509 

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

511 sky = stars["sky_source"] 

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

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

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

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

516 

517 def test_photometry(self): 

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

519 and that the exposure is calibrated. 

520 """ 

521 calibrate = CalibrateImageTask(config=self.config) 

522 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

524 compute_psf_struct = calibrate._compute_psf(self.attributes, self.id_generator) 

525 psf_stars = compute_psf_struct.detections_sources 

526 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

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

528 calibrate._fit_astrometry(self.exposure, stars) 

529 

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

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

532 

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

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

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

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

537 # and the background should be calibrated to match. 

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

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

540 uncalibrated /= self.photo_calib 

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

542 # PhotoCalib on the exposure must be identically 1. 

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

544 

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

546 # sky sources. 

547 sky = stars["sky_source"] 

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

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

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

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

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

553 # quality here. 

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

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

556 rtol=0.1) 

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

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

559 rtol=0.01) 

560 

561 def test_match_psf_stars(self): 

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

563 and candidates. 

564 """ 

565 calibrate = CalibrateImageTask(config=self.config) 

566 compute_psf_struct = calibrate._compute_psf(self.attributes, self.id_generator) 

567 psf_stars = compute_psf_struct.detections_sources 

568 calibrate._measure_aperture_correction(self.exposure, psf_stars) 

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

570 

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

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

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

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

575 

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

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

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

579 stars = stars.copy(deep=True) 

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

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

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

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

584 

585 calibrate._match_psf_stars(psf_stars, stars) 

586 

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

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

589 stars.sort(stars.getPsfFluxSlot().getMeasKey()) 

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

591 stars = stars.copy(deep=True) 

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

593 [False, False, False, True, True, True]) 

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

595 [False, False, False, True, True, True]) 

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

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

598 

599 def test_match_psf_stars_no_matches(self): 

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

601 """ 

602 calibrate = CalibrateImageTask(config=self.config) 

603 # Make two catalogs that cannot have matches. 

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

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

606 

607 with self.assertRaisesRegex(NoPsfStarsToStarsMatchError, 

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

609 calibrate._match_psf_stars(psf_stars, stars) 

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

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

612 

613 def test_calibrate_image_illumcorr(self): 

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

615 config = copy.copy(self.config) 

616 config.do_illumination_correction = True 

617 config.psf_subtract_background.doApplyFlatBackgroundRatio = True 

618 config.psf_detection.doApplyFlatBackgroundRatio = True 

619 config.star_background.doApplyFlatBackgroundRatio = True 

620 config.star_detection.doApplyFlatBackgroundRatio = True 

621 

622 calibrate = CalibrateImageTask(config=config) 

623 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

625 

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

627 background_flat = self.exposure.clone() 

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

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

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

631 

632 # And create an illumination correction of 1.1. 

633 illum_corr_value = 1.1 

634 illumination_correction = self.exposure.clone() 

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

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

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

638 

639 result = calibrate.run( 

640 exposures=self.exposure, 

641 id_generator=self.id_generator, 

642 background_flat=background_flat, 

643 illumination_correction=illumination_correction, 

644 ) 

645 

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

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

648 # the illum_corr_value. 

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

650 # correction. 

651 self.assertFloatsAlmostEqual( 

652 result.applied_photo_calib.getCalibrationMean(), 

653 self.photo_calib * illum_corr_value, 

654 rtol=1e-2, 

655 ) 

656 

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

658 self.assertFloatsAlmostEqual( 

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

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

661 rtol=1e-3, 

662 ) 

663 

664 # Check metadata. 

665 key = "LSST CALIB ILLUMCORR APPLIED" 

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

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

668 

669 def test_run_with_diffraction_spike_mask(self): 

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

671 """ 

672 config = self.config 

673 config.doMaskDiffractionSpikes = True 

674 config.diffractionSpikeMask.magnitudeThreshold = 19 

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

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

677 config.diffractionSpikeMask.saturatedMaskPlane = "FAKESATURATED" 

678 calibrate = CalibrateImageTask(config=config) 

679 calibrate.astrometry.setRefObjLoader(self.ref_loader) 

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

681 calibrate.diffractionSpikeMask.setRefObjLoader(self.ref_loader) 

682 

683 exposure = self.exposure.clone()