Coverage for tests / test_coadds.py: 18%

1# This file is part of cell_coadds. 

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 collections.abc import Iterable, Mapping 

24from itertools import product 

25 

26import numpy as np 

27from frozendict import frozendict 

28 

29import lsst.cell_coadds.test_utils as test_utils 

30import lsst.geom as geom 

31import lsst.meas.base.tests 

32import lsst.utils.tests 

33from lsst.afw.geom import Polygon, Quadrupole 

34from lsst.afw.image import ExposureF, ImageF 

35from lsst.cell_coadds import ( 

36 CellCoaddFitsReader, 

37 CellIdentifiers, 

38 CoaddInputs, 

39 CoaddUnits, 

40 CommonComponents, 

41 ExplodedCoadd, 

42 MultipleCellCoadd, 

43 ObservationIdentifiers, 

44 OwnedImagePlanes, 

45 PatchIdentifiers, 

46 SingleCellCoadd, 

47 StitchedCoadd, 

48 UniformGrid, 

49) 

50from lsst.meas.algorithms import SingleGaussianPsf 

51from lsst.skymap import Index2D 

52 

53 

54class BaseMultipleCellCoaddTestCase(lsst.utils.tests.TestCase): 

55 """A base class that provides a common set of methods.""" 

56 

57 psf_size: int 

58 psf_sigmas: Mapping[Index2D, float] 

59 border_size: int 

60 inner_size: int 

61 outer_size: int 

62 test_positions: Iterable[tuple[geom.Point2D, Index2D]] 

63 exposures: Mapping[Index2D, ExposureF] 

64 multiple_cell_coadd: MultipleCellCoadd 

65 

66 @classmethod 

67 def setUpClass(cls) -> None: 

68 """Set up a multiple cell coadd with 2x2 cells.""" 

69 np.random.seed(42) 

70 

71 cls.nx, cls.ny = 3, 2 

72 cls.psf_sigmas = { 

73 Index2D(x=0, y=0): 1.2, 

74 Index2D(x=0, y=1): 0.7, 

75 Index2D(x=1, y=0): 0.9, 

76 Index2D(x=1, y=1): 1.1, 

77 Index2D(x=2, y=0): 1.3, 

78 Index2D(x=2, y=1): 0.8, 

79 } 

80 

81 cls.border_size = 5 

82 # In practice, we expect this to be squares. 

83 # To check for any possible x, y swap, we use different values. 

84 cls.psf_size_x, cls.psf_size_y = 21, 23 

85 cls.inner_size_x, cls.inner_size_y = 17, 15 

86 cls.outer_size_x = cls.inner_size_x + 2 * cls.border_size 

87 cls.outer_size_y = cls.inner_size_y + 2 * cls.border_size 

88 # The origin should not be at (0, 0) for robust testing. 

89 cls.x0, cls.y0 = 5, 2 

90 cls.n_noise_realizations = 2 

91 

92 patch_outer_bbox = geom.Box2I( 

93 geom.Point2I(cls.x0, cls.y0), geom.Extent2I(cls.nx * cls.inner_size_x, cls.ny * cls.inner_size_y) 

94 ) 

95 patch_outer_bbox.grow(cls.border_size) 

96 

97 # Add one star and one galaxy per quadrant. 

98 # The mapping of positions to cell indices assume inner_size = (17, 15) 

99 # and border_size = 5. If that is changed, these values need an update. 

100 sources = ( 

101 # flux, centroid, shape 

102 (1000.0, geom.Point2D(cls.x0 + 6.3, cls.y0 + 7.2), None), 

103 (2500.0, geom.Point2D(cls.x0 + 16.8, cls.y0 + 18.3), None), 

104 (1500.0, geom.Point2D(cls.x0 + 21.2, cls.y0 + 5.1), None), 

105 (3200.0, geom.Point2D(cls.x0 + 16.1, cls.y0 + 23.9), None), 

106 (1800.0, geom.Point2D(cls.x0 + 44.7, cls.y0 + 8.9), None), 

107 (2100.0, geom.Point2D(cls.x0 + 34.1, cls.y0 + 19.2), None), 

108 (900.0, geom.Point2D(cls.x0 + 9.1, cls.y0 + 13.9), Quadrupole(2.5, 1.5, 0.8)), 

109 (1250.0, geom.Point2D(cls.x0 + 19.3, cls.y0 + 11.2), Quadrupole(1.5, 2.5, 0.75)), 

110 (2100.0, geom.Point2D(cls.x0 + 5.1, cls.y0 + 21.2), Quadrupole(1.7, 1.9, 0.05)), 

111 (2800.0, geom.Point2D(cls.x0 + 24.1, cls.y0 + 19.2), Quadrupole(1.9, 1.7, 0.1)), 

112 (2350.0, geom.Point2D(cls.x0 + 40.3, cls.y0 + 13.9), Quadrupole(1.8, 1.8, -0.4)), 

113 (4999.0, geom.Point2D(cls.x0 + 45.8, cls.y0 + 22.0), Quadrupole(1.6, 1.2, 0.2)), 

114 ) 

115 

116 # The test points are chosen to cover various corner cases assuming 

117 # inner_size = (17, 15) and border_size = 5. If that is changed, the 

118 # test points should be updated to not fall outside the coadd and still 

119 # cover the description in the inline comments. This is checked below. 

120 test_points = ( 

121 geom.Point2D(cls.x0 + 5, cls.y0 + 4), # inner point in lower left 

122 geom.Point2D(cls.x0 + 6, cls.y0 + 24), # inner point in upper left 

123 geom.Point2D(cls.x0 + 25.2, cls.y0 + 7.8), # inner point in lower middle 

124 geom.Point2D(cls.x0 + 23, cls.y0 + 22), # inner point in upper middle 

125 geom.Point2D(cls.x0 + 39, cls.y0 + 9.4), # inner point in lower right 

126 geom.Point2D(cls.x0 + 44, cls.y0 + 24), # inner point in upper right 

127 # Some points that lie on the border 

128 geom.Point2D(cls.x0 + 33, cls.y0 + 24), # inner point in upper right 

129 geom.Point2D(cls.x0 + 46, cls.y0 + 0), # inner point in lower right 

130 geom.Point2D(cls.x0 + 19, cls.y0 + 16), # inner point in upper middle 

131 geom.Point2D(cls.x0 + 17, cls.y0 + 8), # inner point in lower middle 

132 geom.Point2D(cls.x0 + 0, cls.y0 + 29), # inner point in upper left 

133 geom.Point2D(cls.x0 + 0, cls.y0 + 0), # inner point in lower left 

134 ) 

135 # A list of (point, cell_index) pairs. 

136 cls.test_positions = [ 

137 ( 

138 point, 

139 Index2D( 

140 x=int((point.getX() - cls.x0) // cls.inner_size_x), 

141 y=int((point.getY() - cls.y0) // cls.inner_size_y), 

142 ), 

143 ) 

144 for point in test_points 

145 ] 

146 

147 schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema() 

148 

149 single_cell_coadds = [] 

150 cls.exposures = dict.fromkeys(cls.psf_sigmas.keys()) 

151 

152 data_id = test_utils.generate_data_id() 

153 common = CommonComponents( 

154 units=CoaddUnits.nJy, # units here are arbitrary. 

155 wcs=test_utils.generate_wcs(), 

156 band=data_id["band"], 

157 identifiers=PatchIdentifiers.from_data_id(data_id), 

158 visit_polygons={ # This is arbitrary, just to check it if it goes through FITS persistence. 

159 ObservationIdentifiers( 

160 instrument="dummy", 

161 physical_filter="dummy-I", 

162 visit=12345, 

163 detector=67, 

164 day_obs=20000101, 

165 ): Polygon(geom.Box2D(patch_outer_bbox)) 

166 }, 

167 ) 

168 

169 for x in range(cls.nx): 

170 for y in range(cls.ny): 

171 identifiers = CellIdentifiers( 

172 cell=Index2D(x=x, y=y), 

173 skymap=common.identifiers.skymap, 

174 tract=common.identifiers.tract, 

175 patch=common.identifiers.patch, 

176 band=common.identifiers.band, 

177 ) 

178 

179 outer_bbox = geom.Box2I( 

180 geom.Point2I(cls.x0 + x * cls.inner_size_x, cls.y0 + y * cls.inner_size_y), 

181 geom.Extent2I(cls.inner_size_x, cls.inner_size_y), 

182 ) 

183 outer_bbox.grow(cls.border_size) 

184 

185 dataset = lsst.meas.base.tests.TestDataset( 

186 patch_outer_bbox, psfSigma=cls.psf_sigmas[identifiers.cell] 

187 ) 

188 

189 for inst_flux, position, shape in sources: 

190 dataset.addSource(inst_flux, position, shape) 

191 

192 # Create a spatially varying variance plane. 

193 variance = ImageF( 

194 # np.random.uniform returns an array with x-y flipped. 

195 np.random.uniform( 

196 0.8, 

197 1.2, 

198 ( 

199 cls.ny * cls.inner_size_y + 2 * cls.border_size, 

200 cls.nx * cls.inner_size_x + 2 * cls.border_size, 

201 ), 

202 ).astype(np.float32), 

203 xy0=outer_bbox.getMin(), 

204 ) 

205 noise_realizations = [ 

206 ImageF( 

207 np.random.normal( 

208 0.0, 

209 1.0, 

210 ( 

211 cls.ny * cls.inner_size_y + 2 * cls.border_size, 

212 cls.nx * cls.inner_size_x + 2 * cls.border_size, 

213 ), 

214 ).astype(np.float32) 

215 * variance.array, 

216 xy0=outer_bbox.getMin(), 

217 ) 

218 for _ in range(cls.n_noise_realizations) 

219 ] 

220 mask_fraction = ImageF( 

221 np.random.uniform( 

222 0.0, 

223 1.0, 

224 ( 

225 cls.ny * cls.inner_size_y + 2 * cls.border_size, 

226 cls.nx * cls.inner_size_x + 2 * cls.border_size, 

227 ), 

228 ).astype(np.float32), 

229 xy0=outer_bbox.getMin(), 

230 ) 

231 exposure, _ = dataset.realize(variance.getArray() ** 0.5, schema, randomSeed=123456789) 

232 cls.exposures[identifiers.cell] = exposure 

233 exposure = exposure[outer_bbox] 

234 image_plane = OwnedImagePlanes( 

235 image=exposure.image, 

236 variance=exposure.variance, 

237 mask=exposure.mask, 

238 noise_realizations=noise_realizations, 

239 mask_fractions=mask_fraction, 

240 ) 

241 aperture_correction_map = frozendict( 

242 base_GaussianFlux_instFlux=0.9 * (x + 1), 

243 base_GaussianFlux_instFluxErr=0.01 * (y + 1), 

244 base_PsfFlux_instFlux=0.8 * (y + 2), 

245 base_PsfFlux_instFluxErr=0.02 * (x + 2), 

246 ) 

247 

248 single_cell_coadds.append( 

249 SingleCellCoadd( 

250 outer=image_plane, 

251 psf=SingleGaussianPsf( 

252 cls.psf_size_x, cls.psf_size_y, cls.psf_sigmas[Index2D(x=x, y=y)] 

253 ).computeKernelImage(outer_bbox.getCenter()), 

254 inner_bbox=geom.Box2I( 

255 geom.Point2I(cls.x0 + x * cls.inner_size_x, cls.y0 + y * cls.inner_size_y), 

256 geom.Extent2I(cls.inner_size_x, cls.inner_size_y), 

257 ), 

258 inputs={ 

259 ObservationIdentifiers( 

260 instrument="dummy", 

261 physical_filter="dummy-I", 

262 visit=12345, 

263 detector=67, 

264 day_obs=20000101, 

265 ): CoaddInputs( 

266 True, 

267 1.0, 

268 1.0, 

269 Quadrupole( 

270 cls.psf_sigmas[Index2D(x=x, y=y)] ** 2, 

271 cls.psf_sigmas[Index2D(x=x, y=y)] ** 2, 

272 ), 

273 False, 

274 ) 

275 }, 

276 common=common, 

277 identifiers=identifiers, 

278 aperture_correction_map=aperture_correction_map, 

279 ) 

280 ) 

281 

282 grid_bbox = geom.Box2I( 

283 geom.Point2I(cls.x0, cls.y0), geom.Extent2I(cls.nx * cls.inner_size_x, cls.ny * cls.inner_size_y) 

284 ) 

285 grid = UniformGrid.from_bbox_shape(grid_bbox, Index2D(x=cls.nx, y=cls.ny), padding=cls.border_size) 

286 

287 for test_point, test_index in cls.test_positions: 

288 assert test_index == grid.index( 

289 geom.Point2I(test_point) 

290 ), f"Test point {test_point} is not in cell {test_index}." 

291 

292 cls.multiple_cell_coadd = MultipleCellCoadd( 

293 single_cell_coadds, 

294 grid=grid, 

295 outer_cell_size=geom.Extent2I(cls.outer_size_x, cls.outer_size_y), 

296 inner_bbox=None, 

297 common=common, 

298 psf_image_size=geom.Extent2I(cls.psf_size_x, cls.psf_size_y), 

299 ) 

300 

301 @classmethod 

302 def tearDownClass(cls) -> None: # noqa: D102 

303 # Docstring inherited 

304 del cls.multiple_cell_coadd 

305 del cls.exposures 

306 super().tearDownClass() 

307 

308 def assertMultipleCellCoaddsEqual(self, mcc1: MultipleCellCoadd, mcc2: MultipleCellCoadd) -> None: 

309 """Check the equality of two instances of `MultipleCellCoadd`. 

310 

311 Parameters 

312 ---------- 

313 mcc1 : `MultipleCellCoadd` 

314 The MultipleCellCoadd created by reading a FITS file. 

315 mcc2 : `MultipleCellCoadd` 

316 The reference MultipleCellCoadd for comparison. 

317 """ 

318 self.assertEqual(mcc1.band, mcc2.band) 

319 self.assertEqual(mcc1.identifiers, mcc2.identifiers) 

320 self.assertEqual(mcc1.inner_bbox, mcc2.inner_bbox) 

321 self.assertEqual(mcc1.outer_bbox, mcc2.outer_bbox) 

322 self.assertEqual(mcc1.outer_cell_size, mcc2.outer_cell_size) 

323 self.assertEqual(mcc1.grid, mcc2.grid) 

324 self.assertEqual(mcc1.mask_fraction_names, mcc2.mask_fraction_names) 

325 self.assertEqual(mcc1.n_noise_realizations, mcc2.n_noise_realizations) 

326 self.assertEqual(mcc1.psf_image_size, mcc2.psf_image_size) 

327 self.assertEqual(mcc1.units, mcc2.units) 

328 self.assertEqual(mcc1.wcs.getFitsMetadata().toString(), mcc2.wcs.getFitsMetadata().toString()) 

329 

330 # Check that the individual cells are identical. 

331 self.assertEqual(mcc1.cells.keys(), mcc2.cells.keys()) 

332 for idx in mcc1.cells.keys(): # noqa: SIM118 

333 self.assertImagesEqual(mcc1.cells[idx].outer.image, mcc2.cells[idx].outer.image) 

334 self.assertMasksEqual(mcc1.cells[idx].outer.mask, mcc2.cells[idx].outer.mask) 

335 self.assertImagesEqual(mcc1.cells[idx].outer.variance, mcc2.cells[idx].outer.variance) 

336 self.assertImagesEqual(mcc1.cells[idx].outer.mask_fractions, mcc2.cells[idx].outer.mask_fractions) 

337 self.assertEqual( 

338 len(mcc1.cells[idx].outer.noise_realizations), len(mcc2.cells[idx].outer.noise_realizations) 

339 ) 

340 for noise_id in range(len(mcc1.cells[idx].outer.noise_realizations)): 

341 with self.subTest(noise_id): 

342 self.assertImagesEqual( 

343 mcc1.cells[idx].outer.noise_realizations[noise_id], 

344 mcc2.cells[idx].outer.noise_realizations[noise_id], 

345 ) 

346 self.assertImagesEqual(mcc1.cells[idx].psf_image, mcc2.cells[idx].psf_image) 

347 self.assertEqual(mcc1.cells[idx].inputs, mcc2.cells[idx].inputs) 

348 

349 self.assertEqual(mcc1.cells[idx].band, mcc1.band) 

350 self.assertEqual(mcc1.cells[idx].common, mcc1.common) 

351 self.assertEqual(mcc1.cells[idx].units, mcc2.units) 

352 self.assertEqual(mcc1.cells[idx].wcs, mcc1.wcs) 

353 self.assertEqual(mcc1.cells[idx].aperture_correction_map, mcc2.cells[idx].aperture_correction_map) 

354 

355 # Identifiers differ because of the ``cell`` component. 

356 # Check the other attributes within the identifiers. 

357 for attr in ("skymap", "tract", "patch", "band"): 

358 self.assertEqual(getattr(mcc1.cells[idx].identifiers, attr), getattr(mcc1.identifiers, attr)) 

359 

360 

361class MultipleCellCoaddTestCase(BaseMultipleCellCoaddTestCase): 

362 """Test the construction and interfaces of MultipleCellCoadd.""" 

363 

364 def test_fits(self): 

365 """Test that we can write a coadd to a FITS file and read it.""" 

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

367 self.multiple_cell_coadd.write_fits(filename) 

368 mcc1 = MultipleCellCoadd.read_fits(filename) # Test the readFits method. 

369 

370 # Test the reader class. 

371 reader = CellCoaddFitsReader(filename) 

372 mcc2 = reader.readAsMultipleCellCoadd() 

373 

374 wcs = reader.readWcs() 

375 

376 self.assertMultipleCellCoaddsEqual(mcc1, self.multiple_cell_coadd) 

377 self.assertMultipleCellCoaddsEqual(mcc2, self.multiple_cell_coadd) 

378 # By transititve property of equality, mcc1 == mcc2. 

379 

380 self.assertEqual(self.multiple_cell_coadd.band, self.multiple_cell_coadd.common.band) 

381 for poly in self.multiple_cell_coadd.common.visit_polygons.values(): 

382 self.assertEqual(len(poly), 4) 

383 self.assertEqual( 

384 wcs.getFitsMetadata().toString(), self.multiple_cell_coadd.wcs.getFitsMetadata().toString() 

385 ) 

386 

387 def test_visit_count(self): 

388 """Test the visit_count method.""" 

389 # Since we don't simulate coaddition from multiple warps, the cells are 

390 # all going to have just a single visit. 

391 for cellId, singleCellCoadd in self.multiple_cell_coadd.cells.items(): 

392 with self.subTest(x=cellId.x, y=cellId.y): 

393 self.assertEqual(singleCellCoadd.visit_count, 1) 

394 

395 def test_aperture_correction(self): 

396 """Test the aperture correction values are what we expect.""" 

397 for cellId, single_cell_coadd in self.multiple_cell_coadd.cells.items(): 

398 with self.subTest(x=cellId.x, y=cellId.y): 

399 ap_corr_map = single_cell_coadd.aperture_correction_map 

400 self.assertIsNotNone(ap_corr_map) 

401 self.assertEqual(ap_corr_map["base_GaussianFlux_instFlux"], 0.9 * (cellId.x + 1)) 

402 self.assertEqual(ap_corr_map["base_GaussianFlux_instFluxErr"], 0.01 * (cellId.y + 1)) 

403 self.assertEqual(ap_corr_map["base_PsfFlux_instFlux"], 0.8 * (cellId.y + 2)) 

404 self.assertEqual(ap_corr_map["base_PsfFlux_instFluxErr"], 0.02 * (cellId.x + 2)) 

405 

406 

407class ExplodedCoaddTestCase(BaseMultipleCellCoaddTestCase): 

408 """Test the construction and methods of an ExplodedCoadd instance.""" 

409 

410 exploded_coadd: ExplodedCoadd 

411 

412 @classmethod 

413 def setUpClass(cls) -> None: # noqa: D102 

414 # Docstring inherited 

415 super().setUpClass() 

416 cls.exploded_coadd = cls.multiple_cell_coadd.explode() 

417 

418 @classmethod 

419 def tearDownClass(cls) -> None: # noqa: D102 

420 # Docstring inherited 

421 del cls.exploded_coadd 

422 super().tearDownClass() 

423 

424 def test_exploded_psf_image(self): 

425 """Show that psf_image sizes are absurd.""" 

426 self.assertEqual( 

427 self.exploded_coadd.psf_image.getBBox().getDimensions(), 

428 geom.Extent2I(self.nx * self.psf_size_x, self.ny * self.psf_size_y), 

429 ) 

430 for pad_psfs_with in (-999, -4, 0, 4, 8, 21, 40, 100): 

431 exploded_coadd = self.multiple_cell_coadd.explode(pad_psfs_with=pad_psfs_with) 

432 self.assertEqual( 

433 exploded_coadd.psf_image.getBBox().getDimensions(), 

434 geom.Extent2I(self.nx * self.outer_size_x, self.ny * self.outer_size_y), 

435 ) 

436 

437 def test_asMaskedImage(self): 

438 """Test the asMaskedImage method for an ExplodedCoadd object.""" 

439 masked_image = self.exploded_coadd.asMaskedImage() 

440 masked_image.setXY0(self.multiple_cell_coadd.outer_bbox.getMin()) 

441 base_bbox = self.multiple_cell_coadd.cells.first.outer.bbox 

442 for cell_x, cell_y in product(range(self.nx), range(self.ny)): 

443 bbox = base_bbox.shiftedBy(geom.Extent2I(cell_x * self.outer_size_x, cell_y * self.outer_size_y)) 

444 with self.subTest(cell_x=cell_x, cell_y=cell_y): 

445 self.assertMaskedImagesEqual( 

446 masked_image[bbox], 

447 self.multiple_cell_coadd.cells[Index2D(cell_x, cell_y)].outer.asMaskedImage(), 

448 ) 

449 

450 

451class StitchedCoaddTestCase(BaseMultipleCellCoaddTestCase): 

452 """Test the construction and methods of a StitchedCoadd instance.""" 

453 

454 stitched_coadd: StitchedCoadd 

455 

456 @classmethod 

457 def setUpClass(cls) -> None: # noqa: D102 

458 # Docstring inherited 

459 super().setUpClass() 

460 cls.stitched_coadd = cls.multiple_cell_coadd.stitch() 

461 

462 @classmethod 

463 def tearDownClass(cls) -> None: # noqa: D102 

464 # Docstring inherited 

465 del cls.stitched_coadd 

466 super().tearDownClass() 

467 

468 def test_computeBBox(self): 

469 """Test the computeBBox method for a StitchedPsf object.""" 

470 stitched_psf = self.stitched_coadd.psf 

471 

472 psf_bbox = geom.Box2I( 

473 geom.Point2I(-(self.psf_size_x // 2), -(self.psf_size_y // 2)), 

474 geom.Extent2I(self.psf_size_x, self.psf_size_y), 

475 ) 

476 

477 for position, _ in self.test_positions: 

478 bbox = stitched_psf.computeBBox(position) 

479 self.assertEqual(bbox, psf_bbox) 

480 

481 def test_computeShape(self): 

482 """Test the computeShape method for a StitchedPsf object.""" 

483 stitched_psf = self.stitched_coadd.psf 

484 for position, cell_index in self.test_positions: 

485 psf_shape = stitched_psf.computeShape(position) # check we can compute shape 

486 self.assertIsNot(psf_shape.getIxx(), np.nan) 

487 self.assertIsNot(psf_shape.getIyy(), np.nan) 

488 self.assertIsNot(psf_shape.getIxy(), np.nan) 

489 

490 # Moments measured from pixellated images are significantly 

491 # underestimated for small PSFs. 

492 if self.psf_sigmas[cell_index] >= 1.0: 

493 self.assertAlmostEqual(psf_shape.getIxx(), self.psf_sigmas[cell_index] ** 2, delta=1e-3) 

494 self.assertAlmostEqual(psf_shape.getIyy(), self.psf_sigmas[cell_index] ** 2, delta=1e-3) 

495 self.assertAlmostEqual(psf_shape.getIxy(), 0.0) 

496 

497 def test_computeKernelImage(self): 

498 """Test the computeKernelImage method for a StitchedPsf object.""" 

499 stitched_psf = self.stitched_coadd.psf 

500 psf_bbox = geom.Box2I( 

501 geom.Point2I(-(self.psf_size_x // 2), -(self.psf_size_y // 2)), 

502 geom.Extent2I(self.psf_size_x, self.psf_size_y), 

503 ) 

504 

505 for position, cell_index in self.test_positions: 

506 image1 = stitched_psf.computeKernelImage(position) 

507 image2 = SingleGaussianPsf( 

508 self.psf_size_x, self.psf_size_y, self.psf_sigmas[cell_index] 

509 ).computeKernelImage(position) 

510 self.assertImagesEqual(image1, image2) 

511 self.assertEqual(image1.getBBox(), psf_bbox) 

512 

513 def test_computeImage(self): 

514 """Test the computeImage method for a StitchedPsf object.""" 

515 stitched_psf = self.stitched_coadd.psf 

516 psf_extent = geom.Extent2I(self.psf_size_x, self.psf_size_y) 

517 

518 for position, cell_index in self.test_positions: 

519 image1 = stitched_psf.computeImage(position) 

520 image2 = SingleGaussianPsf( 

521 self.psf_size_x, self.psf_size_y, self.psf_sigmas[cell_index] 

522 ).computeImage(position) 

523 self.assertImagesEqual(image1, image2) 

524 self.assertEqual(image1.getBBox().getDimensions(), psf_extent) 

525 

526 def test_computeImage_computeKernelImage(self): 

527 """Test that computeImage called at integer points gives the same 

528 result as calling computeKernelImage. 

529 """ 

530 stitched_psf = self.stitched_coadd.psf 

531 for position, _cell_index in self.test_positions: 

532 pos = geom.Point2D(geom.Point2I(position)) # round to integer 

533 image1 = stitched_psf.computeKernelImage(pos) 

534 image2 = stitched_psf.computeImage(pos) 

535 self.assertImagesEqual(image1, image2) 

536 

537 def test_computeApetureFlux(self): 

538 """Test the computeApertureFlux method for a StitchedPsf object.""" 

539 stitched_psf = self.stitched_coadd.psf 

540 for position, cell_index in self.test_positions: 

541 flux1sigma = stitched_psf.computeApertureFlux(self.psf_sigmas[cell_index], position=position) 

542 self.assertAlmostEqual(flux1sigma, 0.39, delta=5e-2) 

543 

544 flux3sigma = stitched_psf.computeApertureFlux( 

545 3.0 * self.psf_sigmas[cell_index], position=position 

546 ) 

547 self.assertAlmostEqual(flux3sigma, 0.97, delta=2e-2) 

548 

549 def test_asExposure(self): 

550 """Test the asExposure method for a StitchedCoadd object.""" 

551 exposure = self.stitched_coadd.asExposure() 

552 

553 # Check that the bounding box is correct. 

554 bbox = exposure.getBBox() 

555 self.assertEqual(bbox.getWidth(), self.inner_size_x * self.nx + 2 * self.border_size) 

556 self.assertEqual(bbox.getHeight(), self.inner_size_y * self.ny + 2 * self.border_size) 

557 

558 for y in range(self.ny): 

559 for x in range(self.nx): 

560 bbox = geom.Box2I( 

561 geom.Point2I(self.x0 + x * self.inner_size_x, self.y0 + y * self.inner_size_y), 

562 geom.Extent2I(self.inner_size_x, self.inner_size_y), 

563 ) 

564 index = Index2D(x=x, y=y) 

565 self.assertMaskedImagesEqual(exposure[bbox], self.exposures[index][bbox]) 

566 

567 self.assertMaskedImagesEqual(self.stitched_coadd.asExposure(noise_index=None), exposure) 

568 for noise_index in range(self.n_noise_realizations): 

569 noise_exposure = self.stitched_coadd.asExposure(noise_index=noise_index) 

570 

571 self.assertImagesEqual(noise_exposure.variance, exposure.variance) 

572 self.assertImagesEqual(noise_exposure.mask, exposure.mask) 

573 

574 for y in range(self.ny): 

575 for x in range(self.nx): 

576 bbox = geom.Box2I( 

577 geom.Point2I(self.x0 + x * self.inner_size_x, self.y0 + y * self.inner_size_y), 

578 geom.Extent2I(self.inner_size_x, self.inner_size_y), 

579 ) 

580 index = Index2D(x=x, y=y) 

581 self.assertImagesEqual( 

582 noise_exposure.image[bbox], 

583 self.multiple_cell_coadd.cells[index].outer.noise_realizations[noise_index][bbox], 

584 ) 

585 

586 with self.assertRaises(ValueError): 

587 self.stitched_coadd.asExposure(noise_index=self.n_noise_realizations) 

588 

589 def test_aperture_correction(self): 

590 """Test the aperture correction values are what we expect.""" 

591 ap_corr_map = self.stitched_coadd.ap_corr_map 

592 algorithm_names = ap_corr_map.keys() 

593 

594 for (position, cellId), algorithm_name in product(self.test_positions, algorithm_names): 

595 with self.subTest(x=cellId.x, y=cellId.y, algorithm_name=algorithm_name): 

596 field_name = algorithm_name 

597 self.assertEqual( 

598 ap_corr_map[field_name].evaluate(position), 

599 self.multiple_cell_coadd.cells[cellId].aperture_correction_map[field_name], 

600 ) 

601 

602 def test_inputs(self): 

603 """Test that the inputs are populated correctly on stitching.""" 

604 inputs = self.stitched_coadd.ccds 

605 for position, cellId in self.test_positions: 

606 with self.subTest(x=cellId.x, y=cellId.y): 

607 self.assertEqual( 

608 inputs.subsetContaining(position), 

609 tuple(self.multiple_cell_coadd.cells[cellId].inputs.keys()), 

610 ) 

611 self.assertEqual(len(self.stitched_coadd.visits), 1) 

612 

613 def test_coaddInputs(self): 

614 """Test that the inputs are populated when converted to Exposure.""" 

615 inputs = self.stitched_coadd.asExposure().getInfo().getCoaddInputs() 

616 for position, _ in self.test_positions: 

617 ccds = inputs.subset_containing_ccds(position, None) 

618 visits = inputs.subset_containing_visits(position, None) 

619 with self.subTest(x=position.x, y=position.y): 

620 self.assertEqual(len(ccds), 1) 

621 self.assertEqual(len(visits), 1) 

622 

623 def test_borders(self): 

624 """Test that the borders are populated correctly on stitching.""" 

625 mi = self.stitched_coadd.asMaskedImage() 

626 

627 # Check that the borders are not all zero. 

628 with self.assertRaises(AssertionError): 

629 np.testing.assert_array_equal(mi.image.array[: self.border_size, :], 0) 

630 with self.assertRaises(AssertionError): 

631 np.testing.assert_array_equal(mi.image.array[-self.border_size :, :], 0) 

632 with self.assertRaises(AssertionError): 

633 np.testing.assert_array_equal(mi.image.array[:, : self.border_size], 0) 

634 with self.assertRaises(AssertionError): 

635 np.testing.assert_array_equal(mi.image.array[:, -self.border_size :], 0) 

636 

637 def test_fits(self): 

638 """Test that we can write an Exposure with StitchedPsf to a FITS file 

639 and read it. 

640 """ 

641 write_exposure = self.stitched_coadd.asExposure() 

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

643 write_exposure.writeFits(filename) 

644 read_exposure = ExposureF.readFits(filename) # Test the readFits method. 

645 

646 # Test that the image planes are identical. 

647 self.assertMaskedImagesEqual(read_exposure, write_exposure) 

648 

649 # Test the PSF images in the StitchedPsf. 

650 for index in write_exposure.psf.images.indices(): 

651 self.assertImagesEqual(read_exposure.psf.images[index], write_exposure.psf.images[index])