Coverage for tests/test_getTemplate.py: 84%

1# This file is part of ip_diffim.

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.

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/>.

22import collections

23import itertools

24import unittest

26import numpy as np

28import lsst.afw.geom

29import lsst.afw.image

30import lsst.afw.math

31from lsst.daf.butler import DataCoordinate, DimensionUniverse

32import lsst.geom

33import lsst.ip.diffim

34import lsst.meas.algorithms

35import lsst.meas.base.tests

36import lsst.pipe.base as pipeBase

37import lsst.skymap

38import lsst.utils.tests

40from utils import generate_data_id

42# Change this to True, `setup display_ds9`, and open ds9 (or use another afw

43# display backend) to show the tract/patch layouts on the image.

44debug = False

45if debug: 45 ↛ 46line 45 didn't jump to line 46 because the condition on line 45 was never true

46 import lsst.afw.display

47 display = lsst.afw.display.Display()

48 display.frame = 1

51def _showTemplate(box, template):

52 """Show the corners of the template we made in this test."""

53 for point in box.getCorners():

54 display.dot("+", point.x, point.y, ctype="orange", size=40)

55 display.frame = 2

56 display.image(template, "warped template")

57 display.frame = 3

58 display.image(template.variance, "warped variance")

61class GetTemplateTaskTestCase(lsst.utils.tests.TestCase):

62 """Test that GetTemplateTask works on both one tract and multiple tract

63 input coadd exposures.

65 Makes a synthetic exposure large enough to fit four small tracts with 2x2

66 (300x300 pixel) patches each, extracts pixels for those patches by warping,

67 and tests GetTemplateTask's output against boxes that overlap various

68 combinations of one or multiple tracts.

69 """

70 def setUp(self):

71 self.scale = 0.2 # arcsec/pixel

72 self.skymap = self._makeSkymap()

73 self.patches = collections.defaultdict(list)

74 self.dataIds = collections.defaultdict(list)

75 self.exposure = self._makeExposure()

77 if debug: 77 ↛ 78line 77 didn't jump to line 78 because the condition on line 77 was never true

78 display.image(self.exposure, "base exposure")

80 for tract_id in range(4):

81 tract = self.skymap.generateTract(tract_id)

82 self._makePatches(tract)

84 def _makeSkymap(self):

85 """Make a Skymap with 4 tracts with 4 patches each.

86 """

87 tractScale = 0.02 # degrees

88 # On-sky coordinates of the tract centers.

89 coords = [(0, 0),

90 (0, tractScale),

91 (tractScale, 0),

92 (tractScale, tractScale),

93 ]

94 config = lsst.skymap.DiscreteSkyMap.ConfigClass()

95 config.raList = [c[0] for c in coords]

96 config.decList = [c[1] for c in coords]

97 # Half the tract center step size, to keep the tract overlap small.

98 config.radiusList = [tractScale/2 for c in coords]

99 config.projection = "TAN"

100 config.pixelScale = self.scale

101 config.tractOverlap = 0.0005

102 config.tractBuilder = "legacy"

103 config.tractBuilder["legacy"].patchInnerDimensions = (300, 300)

104 config.tractBuilder["legacy"].patchBorder = 10

105 return lsst.skymap.DiscreteSkyMap(config=config)

106

107 def _makeExposure(self):

108 """Create a large image to break up into tracts and patches.

109

110 The image will have a source every 100 pixels in x and y, and a WCS

111 that results in the tracts all fitting in the image, with tract=0

112 in the lower left, tract=1 to the right, tract=2 above, and tract=3

113 to the upper right.

114 """

115 box = lsst.geom.Box2I(lsst.geom.Point2I(-200, -200), lsst.geom.Point2I(800, 800))

116 # This WCS was constructed so that tract 0 mostly fills the lower left

117 # quadrant of the image, and the other tracts fill the rest; slight

118 # extra rotation as a check on the final warp layout, scaled by 5%

119 # from the patch pixel scale.

120 cd_matrix = lsst.afw.geom.makeCdMatrix(1.05*self.scale*lsst.geom.arcseconds, 93*lsst.geom.degrees)

121 wcs = lsst.afw.geom.makeSkyWcs(lsst.geom.Point2D(120, 150),

122 lsst.geom.SpherePoint(0, 0, lsst.geom.radians),

123 cd_matrix)

124 dataset = lsst.meas.base.tests.TestDataset(box, wcs=wcs)

125 for x, y in itertools.product(np.arange(0, 500, 100), np.arange(0, 500, 100)):

126 dataset.addSource(1e5, lsst.geom.Point2D(x, y))

127 exposure, _ = dataset.realize(2, dataset.makeMinimalSchema())

128 exposure.setFilter(lsst.afw.image.FilterLabel("a", "a_test"))

129 return exposure

130

131 def _makePatches(self, tract):

132 """Populate the patches and dataId dicts, keyed on tract id, with the

133 warps of the main exposure and minimal dataIds, respectively.

134 """

135 if debug: 135 ↛ 136line 135 didn't jump to line 136 because the condition on line 135 was never true

136 color = ['red', 'green', 'cyan', 'yellow'][tract.tract_id]

137 point = self.exposure.wcs.skyToPixel(tract.ctr_coord)

138 # Show the tract center, colored by tract id.

139 display.dot("x", point.x, point.y, ctype=color, size=30)

140

141 # Use 5th order to minimize artifacts on the templates.

142 config = lsst.afw.math.Warper.ConfigClass()

143 config.warpingKernelName = "lanczos5"

144 warper = lsst.afw.math.Warper.fromConfig(config)

145 for patchId in range(tract.num_patches.x*tract.num_patches.y):

146 patch = tract.getPatchInfo(patchId)

147 box = patch.getOuterBBox()

148

149 if debug: 149 ↛ 151line 149 didn't jump to line 151 because the condition on line 149 was never true

150 # Show the patch corners as patch ids, colored by tract id.

151 points = self.exposure.wcs.skyToPixel(patch.wcs.pixelToSky([lsst.geom.Point2D(x)

152 for x in box.getCorners()]))

153 for p in points:

154 display.dot(patchId, p.x, p.y, ctype=color)

155

156 # This is mostly taken from drp_tasks makePsfMatchedWarp, but

157 # ip_diffim cannot depend on drp_tasks.

158 xyTransform = lsst.afw.geom.makeWcsPairTransform(self.exposure.wcs, patch.wcs)

159 warpedPsf = lsst.meas.algorithms.WarpedPsf(self.exposure.psf, xyTransform)

160 warped = warper.warpExposure(patch.wcs, self.exposure, destBBox=box)

161 warped.setPsf(warpedPsf)

162 dataRef = pipeBase.InMemoryDatasetHandle(

163 warped,

164 storageClass="ExposureF",

165 copy=True,

166 dataId=generate_data_id(

167 tract=tract,

168 patch=patch,

169 )

170 )

171 self.patches[tract.tract_id].append(dataRef)

172 dataCoordinate = DataCoordinate.standardize({"tract": tract.tract_id,

173 "patch": patchId,

174 "band": "a",

175 "skymap": "skymap"},

176 universe=DimensionUniverse())

177 self.dataIds[tract.tract_id].append(dataCoordinate)

178

179 def _checkMetadata(self, template, config, box, wcs, nPsfs):

180 """Check that the various metadata components were set correctly.

181 """

182 expectedBox = lsst.geom.Box2I(box)

183 expectedBox.grow(config.templateBorderSize)

184 self.assertEqual(template.getBBox(), expectedBox)

185 # WCS should match our exposure, not any of the coadd tracts.

186 for tract in self.patches:

187 self.assertNotEqual(template.wcs, self.patches[tract][0].get().wcs)

188 self.assertEqual(template.wcs, self.exposure.wcs)

189 self.assertEqual(template.photoCalib, self.exposure.photoCalib)

190 self.assertEqual(template.getXY0(), expectedBox.getMin())

191 self.assertEqual(template.filter.bandLabel, "a")

192 self.assertEqual(template.filter.physicalLabel, "a_test")

193 self.assertEqual(template.psf.getComponentCount(), nPsfs)

194 self.assertTrue(template.getInfo().hasCoaddInputs())

195 self.assertEqual(len(template.getInfo().getCoaddInputs().ccds), nPsfs)

196

197 def _checkPixels(self, template, config, box):

198 """Check that the pixel values in the template are close to the

199 original image.

200 """

201 # All pixels should have real values!

202 expectedBox = lsst.geom.Box2I(box)

203 expectedBox.grow(config.templateBorderSize)

204

205 if debug: 205 ↛ 206line 205 didn't jump to line 206 because the condition on line 205 was never true

206 _showTemplate(expectedBox, template)

207

208 # Check that we fully filled the template from the patches.

209 self.assertTrue(np.all(np.isfinite(template.image.array)))

210 # Because of the scale changes, there will be some ringing in the

211 # difference between the template and the original image; pick

212 # tolerances large enough to account for that.

213 self.assertImagesAlmostEqual(template.image, self.exposure[expectedBox].image,

214 rtol=.1, atol=4)

215 # Variance plane ==2 in the original image, but the warped images will

216 # have some structure due to the warping.

217 self.assertImagesAlmostEqual(template.variance, self.exposure[expectedBox].variance,

218 rtol=0.55, msg="variance planes differ")

219 # Not checking the mask, as warping changes the sizes of the masks.

220

221 def testRunOneTractInput(self):

222 """Test a bounding box that fully fits inside one tract, with only

223 that tract passed as input. This checks that the code handles a single

224 tract input correctly.

225 """

226 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))

227 task = lsst.ip.diffim.GetTemplateTask()

228 # Restrict to tract 0, since the box fits in just that tract.

229 # Task modifies the input bbox, so pass a copy.

230 result = task.run(coaddExposureHandles={0: self.patches[0]},

231 bbox=lsst.geom.Box2I(box),

232 wcs=self.exposure.wcs,

233 dataIds={0: self.dataIds[0]},

234 physical_filter="a_test")

235

236 # All 4 patches from tract 0 are included in this template.

237 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 4)

238 self._checkPixels(result.template, task.config, box)

239

240 def testRunOneTractMultipleInputs(self):

241 """Test a bounding box that fully fits inside one tract but where

242 multiple tracts were passed in. This checks that patches that are

243 mostly NaN after warping are merged correctly in the output.

244 """

245 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))

246 task = lsst.ip.diffim.GetTemplateTask()

247 # Task modifies the input bbox, so pass a copy.

248 result = task.run(coaddExposureHandles=self.patches,

249 bbox=lsst.geom.Box2I(box),

250 wcs=self.exposure.wcs,

251 dataIds=self.dataIds,

252 physical_filter="a_test")

253

254 # All 4 patches from two tracts are included in this template.

255 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 6)

256 self._checkPixels(result.template, task.config, box)

257

258 def testRunTwoTracts(self):

259 """Test a bounding box that crosses tract boundaries.

260 """

261 box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600))

262 task = lsst.ip.diffim.GetTemplateTask()

263 # Task modifies the input bbox, so pass a copy.

264 result = task.run(coaddExposureHandles=self.patches,

265 bbox=lsst.geom.Box2I(box),

266 wcs=self.exposure.wcs,

267 dataIds=self.dataIds,

268 physical_filter="a_test")

269

270 # All 4 patches from all 4 tracts are included in this template

271 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 9)

272 self._checkPixels(result.template, task.config, box)

273

274 def testRunNoTemplate(self):

275 """A bounding box that doesn't overlap the patches will raise.

276 """

277 box = lsst.geom.Box2I(lsst.geom.Point2I(1200, 1200), lsst.geom.Point2I(1600, 1600))

278 task = lsst.ip.diffim.GetTemplateTask()

279 with self.assertRaisesRegex(lsst.pipe.base.NoWorkFound, "No patches found"):

280 task.run(coaddExposureHandles=self.patches,

281 bbox=lsst.geom.Box2I(box),

282 wcs=self.exposure.wcs,

283 dataIds=self.dataIds,

284 physical_filter="a_test")

285

286 def testMissingPatches(self):

287 """Test that a missing patch results in an appropriate mask.

288

289 This fixes the bug reported on DM-44997 (image and variance were NaN

290 but the mask was not set to NO_DATA for those pixels).

291 """

292 # tract=0, patch=1 is the lower-left corner, as displayed in DS9.

293 self.patches[0].pop(1)

294 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))

295 task = lsst.ip.diffim.GetTemplateTask()

296 # Task modifies the input bbox, so pass a copy.

297 result = task.run(coaddExposureHandles=self.patches,

298 bbox=lsst.geom.Box2I(box),

299 wcs=self.exposure.wcs,

300 dataIds=self.dataIds,

301 physical_filter="a_test")

302 no_data = (result.template.mask.array & result.template.mask.getPlaneBitMask("NO_DATA")) != 0

303 self.assertTrue(np.isfinite(result.template.image.array).all())

304 self.assertTrue(np.isfinite(result.template.variance.array).all())

305 self.assertEqual(no_data.sum(), 20990)

306

307 @lsst.utils.tests.methodParameters(

308 box=[

309 lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)),

310 lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)),

311 ],

312 nInput=[8, 16],

313 )

314 def testNanInputs(self, box=None, nInput=None):

315 """Test that the template has finite values when some of the input

316 pixels have NaN as variance.

317 """

318 for tract, patchRefs in self.patches.items():

319 for patchRef in patchRefs:

320 patchCoadd = patchRef.get()

321 bbox = lsst.geom.Box2I()

322 bbox.include(lsst.geom.Point2I(patchCoadd.getBBox().getCenter()))

323 bbox.grow(3)

324 patchCoadd.variance[bbox].array *= np.nan

325

326 box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600))

327 task = lsst.ip.diffim.GetTemplateTask()

328 result = task.run(coaddExposureHandles=self.patches,

329 bbox=lsst.geom.Box2I(box),

330 wcs=self.exposure.wcs,

331 dataIds=self.dataIds,

332 physical_filter="a_test")

333 if debug: 333 ↛ 334line 333 didn't jump to line 334 because the condition on line 333 was never true

334 _showTemplate(box, result.template)

335 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 9)

336 # We just check that the pixel values are all finite. We cannot check that pixel values

337 # in the template are closer to the original anymore.

338 self.assertTrue(np.isfinite(result.template.image.array).all())

339

340

341def setup_module(module):

342 lsst.utils.tests.init()

343

344

345class MemoryTestCase(lsst.utils.tests.MemoryTestCase):

346 pass

347

348

349if __name__ == "__main__": 349 ↛ 350line 349 didn't jump to line 350 because the condition on line 349 was never true

350 lsst.utils.tests.init()

351 unittest.main()