Coverage for tests/test

2# LSST Data Management System

6# This product includes software developed by the

7# LSST Project (http://www.lsst.org/).

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 LSST License Statement and

20# the GNU General Public License along with this program. If not,

21# see <https://www.lsstcorp.org/LegalNotices/>.

22#

23import math

24import unittest

25import sys

27import numpy as np

28import itertools

29import lsst.utils.tests

30import lsst.afw.detection as afwDetection

31import lsst.afw.display as afwDisplay

32import lsst.geom as geom

33import lsst.afw.geom as afwGeom

34import lsst.afw.geom.ellipses as afwEllipses

35import lsst.afw.image as afwImage

36import lsst.afw.math as afwMath

37import lsst.afw.table as afwTable

38import lsst.meas.algorithms as measAlg

39import lsst.meas.base as measBase

40# importing this package registers essential code

41import lsst.meas.extensions.photometryKron

42from lsst.daf.base import PropertyList

44try:

45 display

46except NameError:

47 display = False

49afwDisplay.setDefaultMaskTransparency(75)

52def makeGalaxy(width, height, flux, a, b, theta, dx=0.0, dy=0.0, xy0=None, xcen=None, ycen=None):

53 """Make a fake galaxy image.

54 """

55 gal = afwImage.ImageF(width, height)

56 if xcen is None:

57 xcen = 0.5*width + dx

58 if ycen is None:

59 ycen = 0.5*height + dy

60 I0 = flux/(2*math.pi*a*b)

62 if xy0 is not None:

63 gal.setXY0(xy0)

65 c, s = math.cos(math.radians(theta)), math.sin(math.radians(theta))

66 ii, iuu, ivv = 0.0, 0.0, 0.0

67 for y in range(height):

68 for x in range(width):

69 dx, dy = x + gal.getX0() - xcen, y + gal.getY0() - ycen

70 if math.hypot(dx, dy) < 10.5:

71 nsample = 5

72 subZ = np.linspace(-0.5*(1 - 1/nsample), 0.5*(1 - 1/nsample), nsample)

73 else:

74 nsample = 1

75 subZ = [0.0]

77 val = 0

78 for sx in subZ:

79 for sy in subZ:

80 u = c*(dx + sx) + s*(dy + sy)

81 v = -s*(dx + sx) + c*(dy + sy)

82 val += I0*math.exp(-0.5*((u/a)**2 + (v/b)**2))

84 if val < 0: 84 ↛ 85line 84 didn't jump to line 85 because the condition on line 84 was never true

85 val = 0

86 gal[geom.Point2I(x, y), afwImage.LOCAL] = val/nsample**2

88 ii += val

89 iuu += val*u**2

90 ivv += val*v**2

92 iuu /= ii

93 ivv /= ii

95 exp = afwImage.makeExposure(afwImage.makeMaskedImage(gal))

96 exp.getMaskedImage().getVariance().set(1.0)

97 scale = 1.0e-4*geom.degrees

98 cdMatrix = afwGeom.makeCdMatrix(scale=scale, flipX=True)

99 exp.setWcs(afwGeom.makeSkyWcs(crpix=geom.Point2D(0.0, 0.0),

100 crval=geom.SpherePoint(0.0, 0.0, geom.degrees),

101 cdMatrix=cdMatrix))

102 # add a dummy Psf. The new SdssCentroid needs one

103 exp.setPsf(afwDetection.GaussianPsf(11, 11, 0.01))

104 return exp

105

106

107def makeMeasurementConfig(forced=False, nsigma=6.0, nIterForRadius=1, kfac=2.5):

108 """Construct a (SingleFrame|Forced)MeasurementConfig with the requested parameters.

109 """

110 if forced:

111 msConfig = measBase.ForcedMeasurementConfig()

112 msConfig.algorithms.names = ["base_TransformedCentroid", "base_TransformedShape",

113 "ext_photometryKron_KronFlux"]

114 msConfig.slots.centroid = "base_TransformedCentroid"

115 msConfig.slots.shape = "base_TransformedShape"

116 msConfig.copyColumns = {"id": "objectId", "parent": "parentObjectId"}

117 else:

118 msConfig = measBase.SingleFrameMeasurementConfig()

119 msConfig.algorithms.names = ["base_SdssCentroid", "base_SdssShape",

120 "ext_photometryKron_KronFlux", "base_SkyCoord"]

121 msConfig.slots.centroid = "base_SdssCentroid"

122 msConfig.slots.shape = "base_SdssShape"

123 msConfig.slots.apFlux = "ext_photometryKron_KronFlux"

124 msConfig.slots.modelFlux = None

125 msConfig.slots.psfFlux = None

126 msConfig.slots.gaussianFlux = None

127 msConfig.slots.calibFlux = None

128 # msConfig.algorithms.names.remove("correctfluxes")

129 msConfig.plugins["base_SdssCentroid"].maxDistToPeak = -1

130 msConfig.plugins["ext_photometryKron_KronFlux"].nSigmaForRadius = nsigma

131 msConfig.plugins["ext_photometryKron_KronFlux"].nIterForRadius = nIterForRadius

132 msConfig.plugins["ext_photometryKron_KronFlux"].nRadiusForFlux = kfac

133 msConfig.plugins["ext_photometryKron_KronFlux"].enforceMinimumRadius = False

134 return msConfig

135

136

137def measureFree(exposure, center, msConfig):

138 """Unforced measurement.

139 """

140 schema = afwTable.SourceTable.makeMinimalSchema()

141 algMeta = PropertyList()

142 task = measBase.SingleFrameMeasurementTask(schema, config=msConfig, algMetadata=algMeta)

143 measCat = afwTable.SourceCatalog(schema)

144 source = measCat.addNew()

145 source.getTable().setMetadata(algMeta)

146 ss = afwDetection.FootprintSet(exposure.getMaskedImage(), afwDetection.Threshold(0.1))

147 fp = ss.getFootprints()[0]

148 source.setFootprint(fp)

149 task.run(measCat, exposure)

150 return source

151

152

153def measureForced(exposure, source, refWcs, msConfig):

154 """Forced measurement.

155 """

156 refCat = afwTable.SourceCatalog(source.table)

157 refCat.append(source)

158 schema = afwTable.SourceTable.makeMinimalSchema()

159 algMeta = PropertyList()

160 source.getTable().setMetadata(algMeta)

161 task = measBase.ForcedMeasurementTask(schema, config=msConfig, algMetadata=algMeta)

162 measCat = task.generateMeasCat(exposure, refCat, refWcs)

163 task.attachTransformedFootprints(measCat, refCat, exposure, refWcs)

164 task.run(measCat, exposure, refCat, refWcs)

165 return measCat[0]

166

167

168class KronPhotometryTestCase(lsst.utils.tests.TestCase):

169 """A test case for measuring Kron quantities.

170 """

171

172 def setUp(self):

173 self.flux = 1e5

174 self.width, self.height = 200, 200

175 self.objImg = None

176

177 def tearDown(self):

178 if self.objImg:

179 del self.objImg

180

181 def makeAndMeasure(self, measureKron, a, b, theta, dx=0.0, dy=0.0, nsigma=6, kfac=2, nIterForRadius=1,

182 xcen=None, ycen=None, makeImage=True):

183 """Make and measure an elliptical Gaussian.

184 """

185 if xcen is None:

186 xcen = 0.5*self.width + dx

187 if ycen is None:

188 ycen = 0.5*self.height + dy

189 #

190 # Make the object

191 #

192 if a < b: 192 ↛ 193line 192 didn't jump to line 193 because the condition on line 192 was never true

193 a, b = b, a

194 theta += 90

195

196 currentDisp = None

197 if self.objImg is None:

198 makeImage = True

199 if makeImage:

200 self.objImg = makeGalaxy(self.width, self.height, self.flux, a, b, theta, dx=dx, dy=dy,

201 xy0=geom.Point2I(10, 10), xcen=xcen, ycen=ycen)

202 if display: 202 ↛ 203line 202 didn't jump to line 203 because the condition on line 202 was never true

203 disp = afwDisplay.Display(frame=0)

204 disp.mtv(self.objImg, title=self._testMethodName + ": %g %g" % (a, b))

205 currentDisp = disp

206 if not makeImage and display: 206 ↛ 207line 206 didn't jump to line 207 because the condition on line 206 was never true

207 try:

208 currentDisp = disp

209 except NameError:

210 currentDisp = None

211 if currentDisp: 211 ↛ 212line 211 didn't jump to line 212 because the condition on line 211 was never true

212 disp.pan(xcen, ycen)

213 disp.dot("+", xcen, ycen, size=1, ctype=afwDisplay.RED)

214 c, s = math.cos(math.radians(theta)), math.sin(math.radians(theta))

215 # N.b. add 1/12 in quadrature to allow for pixellisation

216 disp.dot("@:%f,%f,%f" % (nsigma**2*((a**2 + 1/12.0)*c**2 + (b**2 + 1/12.0)*s**2),

217 nsigma**2*(a**2 - b**2)*c*s,

218 nsigma**2*((a**2 + 1/12.0)*s**2 + (b**2 + 1/12.0)*c**2)),

219 xcen, ycen, size=1, ctype=afwDisplay.RED)

220 #

221 # Do the measuring

222 #

223 FWHM = 5

224 ksize = 25 # size of desired kernel

225 self.objImg.setPsf(measAlg.DoubleGaussianPsf(ksize, ksize,

226 FWHM/(2*math.sqrt(2*math.log(2))), 1, 0.1))

227

228 return measureKron(self.objImg, xcen, ycen, nsigma, kfac, nIterForRadius, disp=currentDisp)

229

230 def measureKron(self, objImg, xcen, ycen, nsigma, kfac, nIterForRadius, disp=None):

231 """Measure Kron quantities using the C++ code.

232 """

233 #

234 # Now measure things

235 #

236 center = geom.Point2D(xcen, ycen)

237 msConfig = makeMeasurementConfig(False, nsigma, nIterForRadius, kfac)

238 source = measureFree(objImg, center, msConfig)

239 algMeta = source.getTable().getMetadata()

240 self.assertTrue(algMeta.exists('EXT_PHOTOMETRYKRON_KRONFLUX_NRADIUSFORFLUX'))

241

242 R_K = source.get("ext_photometryKron_KronFlux_radius")

243 flux_K = source.get("ext_photometryKron_KronFlux_instFlux")

244 fluxErr_K = source.get("ext_photometryKron_KronFlux_instFluxErr")

245 flags_K = source.get("ext_photometryKron_KronFlux_flag")

246 if not flags_K:

247 # Forced measurement on the same image should produce exactly the same result

248 msConfig = makeMeasurementConfig(True, nsigma, nIterForRadius, kfac)

249 forced = measureForced(objImg, source, objImg.getWcs(), msConfig)

250 algMeta = source.getTable().getMetadata()

251 self.assertTrue(algMeta.exists('EXT_PHOTOMETRYKRON_KRONFLUX_NRADIUSFORFLUX'))

252 for field in (

253 "ext_photometryKron_KronFlux_instFlux",

254 "ext_photometryKron_KronFlux_instFluxErr",

255 "ext_photometryKron_KronFlux_radius",

256 "ext_photometryKron_KronFlux_flag"

257 ):

258 try:

259 if np.isnan(source.get(field)): 259 ↛ 260line 259 didn't jump to line 260 because the condition on line 259 was never true

260 self.assertTrue(np.isnan(forced.get(field)))

261 else:

262 self.assertFloatsAlmostEqual(source.get(

263 field), forced.get(field), rtol=1.0e-6, atol=None)

264 except AssertionError:

265 print("Failed:", field, source.get(field), forced.get(field))

266 raise

267

268 if disp: 268 ↛ 269line 268 didn't jump to line 269 because the condition on line 268 was never true

269 disp.dot("x", xcen, ycen, ctype=afwDisplay.MAGENTA, size=1)

270 afwDisplay.utils.drawFootprint(source.getFootprint(), display=disp)

271 shape = source.getShape()

272 if True: # nsigma*shape, the radius used to estimate R_K

273 shape = shape.clone()

274 shape.scale(source.get("ext_photometryKron_KronFlux_radius_for_radius")

275 / shape.getDeterminantRadius())

276 disp.dot(shape, xcen, ycen, ctype=afwDisplay.MAGENTA)

277 # Show R_K

278 shape = shape.clone()

279 for r, ct in [(R_K, afwDisplay.BLUE), (R_K*kfac, afwDisplay.CYAN), ]:

280 shape.scale(r/shape.getDeterminantRadius())

281 disp.dot(shape, xcen, ycen, ctype=ct)

282

283 return (R_K, flux_K, fluxErr_K, flags_K,

284 source.get("ext_photometryKron_KronFlux_flag_bad_radius"),

285 source.get("ext_photometryKron_KronFlux_flag_small_radius"))

286

287 def measureKronInPython(self, objImg, xcen, ycen, nsigma, kfac, nIterForRadius, makeImage=None):

288 """Measure the Kron quantities of an elliptical Gaussian in python.

289

290 N.b. only works for XY0 == (0, 0)

291 """

292 #

293 # Measure moments using SDSS shape algorithm

294 #

295 # Note: this code was converted to the new meas_framework, but is not exercised.

296 msConfig = makeMeasurementConfig(False, nsigma, nIterForRadius, kfac)

297 center = geom.Point2D(xcen, ycen)

298 source = self.measureFree(objImg, center, msConfig)

299 algMeta = source.getTable().getMetadata()

300 self.assertTrue(algMeta.exists('EXT_PHOTOMETRYKRON_KRONFLUX_NRADIUSFORFLUX'))

301

302 Mxx = source.getIxx()

303 Mxy = source.getIxy()

304 Myy = source.getIyy()

305 #

306 # Calculate principal axes

307 #

308 Muu_p_Mvv = Mxx + Myy

309 Muu_m_Mvv = math.sqrt((Mxx - Myy)**2 + 4*Mxy**2)

310 Muu = 0.5*(Muu_p_Mvv + Muu_m_Mvv)

311 Mvv = 0.5*(Muu_p_Mvv - Muu_m_Mvv)

312 theta = 0.5*math.atan2(2*Mxy, Mxx - Myy)

313 a = math.sqrt(Muu)

314 b = math.sqrt(Mvv)

315 ab = a/b

316 #

317 # Get footprint

318 #

319 ellipse = afwEllipses.Ellipse(afwEllipses.Axes(nsigma*a, nsigma*b, theta),

320 geom.Point2D(xcen, ycen))

321 fpEllipse = afwDetection.Footprint(ellipse)

322

323 sumI = 0.0

324 sumR = (0.38259771140356325/ab*(1 + math.sqrt(2)*math.hypot(math.fmod(xcen, 1), math.fmod(ycen, 1)))

325 * objImg.image[int(xcen), int(ycen), afwImage.LOCAL])

326

327 gal = objImg.image

328

329 c, s = math.cos(theta), math.sin(theta)

330 for sp in fpEllipse.getSpans():

331 y, x0, x1 = sp.getY(), sp.getX0(), sp.getX1()

332

333 for x in range(x0, x1 + 1):

334 dx, dy = x - xcen, y - ycen

335 u = c*dx + s*dy

336 v = -s*dx + c*dy

337

338 r = math.hypot(u, v*ab)

339 try:

340 val = gal[x, y, afwImage.LOCAL]

341 except Exception:

342 continue

343

344 sumI += val

345 sumR += val*r

346

347 R_K = sumR/sumI

348

349 sumI = 0.0

350 for y in range(self.height):

351 for x in range(self.width):

352 dx, dy = x - xcen, y - ycen

353 u = c*dx + s*dy

354 v = -s*dx + c*dy

355 if math.hypot(u/a, v/b) < kfac:

356 sumI += gal[x, y, afwImage.LOCAL]

357

358 return R_K, sumI, 0, False, False, False

359

360 def testEllipticalGaussian(self):

361 """Test measuring the Kron quantities of an elliptical Gaussian.

362 """

363 ignoreTestFailures = False # if True, keep going after test failures but always generate a failure

364 #

365 # Choose function that does the measuring

366 #

367 if False: # testing only; requires XY0 == (0, 0)

368 measureKron = self.measureKronInPython

369 else:

370 measureKron = self.measureKron

371 #

372 # Make and the objects

373 #

374 ab_vals = (0.5, 1.0, 2.0, 3.0, 5.0, )

375 nIter = 2

376 for dx in (0.0, 0.5, ):

377 for dy in (0.0, 0.5, ):

378 if measureKron == self.measureKronInPython and dx + dy != 0.0: 378 ↛ 379line 378 didn't jump to line 379 because the condition on line 378 was never true

379 continue

380

381 for theta in (0.0, 20.0, 45.0, ):

382 for a in ab_vals:

383 for b in ab_vals:

384 if b > a:

385 continue

386

387 makeImage = True

388 for kfac in (1.5, 2.5, ): # multiple of R_Kron to use for Flux_Kron

389 R_K, flux_K, fluxErr_K, flags_K, flags_radius, flags_small_radius = (

390 self.makeAndMeasure(measureKron, a, b, theta, dx=dx, dy=dy, kfac=kfac,

391 nIterForRadius=nIter, makeImage=makeImage))

392 makeImage = False

393

394 self.assertFalse(flags_radius)

395 self.assertFalse(flags_small_radius)

396 #

397 # We'll have to correct for the pixelisation as we sum over the central

398 # few pixels when making models, mostly do deal with b ~ 0.5 models.

399 #

400 # See Section 5 of

401 # http://www.astro.princeton.edu/~rhl/photomisc/aperture.pdf

402 # for the source of 0.00286 etc.

403 #

404 R_truth0 = math.sqrt(math.pi/2)

405 R_truth = R_truth0*math.sqrt(1 + 0.8*1/(12.0*a*b))

406

407 flux_truth = self.flux*(1 - math.exp(-0.5*(kfac*R_truth)**2))

408 R_truth = R_truth0*math.sqrt(a*b + 1/12.0*(1 + 0.00286/min(a, b)**3.9))

409

410 failR = math.isnan(R_K) or flags_K or (

411 abs(R_truth - R_K) > 1e-2*self.getTolRad(a, b))

412 failFlux = math.isnan(flux_K) or flags_K or (

413 abs(flux_K/flux_truth - 1) > 1e-2*self.getTolFlux(a, b, kfac))

414

415 ID = ("a,b,theta %4.1f %4.1f %4.1f dx,dy = %.1f,%.1f kfac=%g" %

416 (a, b, theta, dx, dy, kfac))

417 if ((failR or failFlux) and verbose) or verbose > 1: 417 ↛ 426line 417 didn't jump to line 426 because the condition on line 417 was always true

418 print("%s R_K %10.3f %10.3f %6.3f pixels (tol %5.3f)%s" %

419 (ID, R_K, R_truth, (R_K - R_truth), 1e-2*self.getTolRad(a, b),

420 " *" if failR else ""))

421 print("%s flux_K %10.3f %10.3f %6.2f%% (tol %5.3f)%s" %

422 (ID, flux_K, flux_truth,

423 100*(flux_K/flux_truth - 1), self.getTolFlux(a, b, kfac),

424 " *" if failFlux else ""))

425

426 if ignoreTestFailures: 426 ↛ 427line 426 didn't jump to line 427 because the condition on line 426 was never true

427 continue

428

429 self.assertFalse(failR, (("%s R_Kron: %g v. exact value %g "

430 "(error %.3f pixels; limit %.3f)") %

431 (ID, R_K, R_truth, (R_K - R_truth),

432 1e-2*self.getTolRad(a, b))))

433

434 self.assertFalse(failFlux,

435 (("%s flux_Kron: %g v. exact value %g "

436 "(error %.2f%% limit %.2f%%)") %

437 (ID, flux_K, flux_truth, 100*(flux_K/flux_truth - 1),

438 self.getTolFlux(a, b, kfac))))

439

440 self.assertFalse(ignoreTestFailures, "You are ignoring possible test failures")

441

442 def testBadFlags(self):

443 a, b, theta = 5, 1, 20.0

444 #

445 # Check that we fail if too close to the edge of the image

446 #

447 for cen in (10, 20, 38, 40, 50, self.width - 10):

448 makeImage = True

449 for kfac in (2.5, 10,):

450 R_K, flux_K, fluxErr_K, flags_K, flags_radius, flags_small_radius = (

451 self.makeAndMeasure(self.measureKron, a, b, theta, makeImage=makeImage,

452 xcen=cen, ycen=cen, kfac=kfac))

453 makeImage = False

454

455 msg = "KronFlux_flag_bad_radius: cen = (%g, %g), kfac = %g" % (cen, cen, kfac)

456

457 if kfac == 2.5 and (cen <= 20 or cen > self.width - 20):

458 self.assertTrue(flags_K, msg)

459 elif kfac == 10:

460 self.assertTrue(flags_K, msg)

461 else:

462 self.assertFalse(flags_K, msg)

463

464 def getTolRad(self, a, b):

465 """Return R_K tolerance in hundredths of a pixel.

466 """

467 if b <= 0.5:

468 if a <= 0.5:

469 tol = 35

470 elif a <= 2:

471 tol = 350

472 else:

473 tol = 25*a # i.e. 0.25*a

474 elif b <= 1:

475 tol = 7.0

476 else:

477 tol = 1.0

478

479 return tol

480

481 def getTolFlux(self, a, b, kfac):

482 """Return Flux_K tolerance in percent.

483 """

484 if b <= 0.5:

485 if a <= 0.5:

486 if kfac > 2:

487 tol = 5.0

488 else:

489 tol = 10.0

490 elif a <= 1.0:

491 if kfac <= 1.5:

492 tol = 10.0

493 else:

494 tol = 4.0

495 else:

496 if kfac > 2:

497 tol = 3.0

498 elif kfac > 1.5: 498 ↛ 499line 498 didn't jump to line 499 because the condition on line 498 was never true

499 tol = 5.0

500 else:

501 tol = 10.0

502 elif b <= 1:

503 if a <= 1:

504 tol = 2.0

505 else:

506 if kfac > 2:

507 tol = 0.25

508 elif kfac > 1.5: 508 ↛ 509line 508 didn't jump to line 509 because the condition on line 508 was never true

509 tol = 0.5

510 else:

511 tol = 1.27

512 elif b <= 2:

513 if kfac > 1.5:

514 tol = 0.1

515 else:

516 tol = 0.5

517 else:

518 tol = 0.30

519

520 return tol

521

522 def testForced(self):

523 """Check that forced photometry works in the presence of rotations and translations.

524 """

525 kfac = 2.5

526 warper = afwMath.Warper("lanczos4")

527 a = 13

528 for axisRatio in (0.25, 1.0):

529 b = a*axisRatio

530 for theta in (0, 30, 45):

531 width, height = 256, 256

532 center = geom.Point2D(0.5*width, 0.5*height)

533 original = makeGalaxy(width, height, 1000.0, a, b, theta)

534 msConfig = makeMeasurementConfig(forced=False, kfac=kfac)

535 source = measureFree(original, center, msConfig)

536 algMeta = source.getTable().getMetadata()

537 self.assertTrue(algMeta.exists('EXT_PHOTOMETRYKRON_KRONFLUX_NRADIUSFORFLUX'))

538 if source.get("ext_photometryKron_KronFlux_flag"): 538 ↛ 539line 538 didn't jump to line 539 because the condition on line 538 was never true

539 continue

540

541 angleList = [val*geom.degrees for val in (45, 90)]

542 scaleList = [1.0, 0.5]

543 offsetList = [(1.23, 4.56), (12.3, 45.6)]

544

545 for angle, scale, offset in itertools.product(angleList, scaleList, offsetList):

546 dx, dy = offset

547 pixelScale = original.getWcs().getPixelScale()*scale

548 cdMatrix = afwGeom.makeCdMatrix(scale=pixelScale, orientation=angle, flipX=True)

549 wcs = afwGeom.makeSkyWcs(crpix=geom.Point2D(dx, dy),

550 crval=geom.SpherePoint(0.0, 0.0, geom.degrees),

551 cdMatrix=cdMatrix)

552

553 warped = warper.warpExposure(wcs, original)

554 # add a Psf if there is none. The new SdssCentroid needs a Psf.

555 if warped.getPsf() is None: 555 ↛ 557line 555 didn't jump to line 557 because the condition on line 555 was always true

556 warped.setPsf(afwDetection.GaussianPsf(11, 11, 0.01))

557 msConfig = makeMeasurementConfig(kfac=kfac, forced=True)

558 forced = measureForced(warped, source, original.getWcs(), msConfig)

559 algMeta = source.getTable().getMetadata()

560 self.assertTrue(algMeta.exists('EXT_PHOTOMETRYKRON_KRONFLUX_NRADIUSFORFLUX'))

561

562 if display: 562 ↛ 563line 562 didn't jump to line 563 because the condition on line 562 was never true

563 disp1 = afwDisplay.Display(frame=1)

564 disp1.mtv(original, title=self._testMethodName + ": original image")

565 shape = source.getShape().clone()

566 xc, yc = source.getCentroid()

567 radius = source.get("ext_photometryKron_KronFlux_radius")

568 for r, ct in [(radius, afwDisplay.BLUE), (radius*kfac, afwDisplay.CYAN), ]:

569 shape.scale(r/shape.getDeterminantRadius())

570 disp1.dot(shape, xc, yc, ctype=ct)

571 disp2 = afwDisplay.Display(frame=2)

572 disp2.mtv(warped, title=self._testMethodName + ": warped image")

573 transform = (wcs.linearizeSkyToPixel(source.getCoord(), lsst.geom.degrees)

574 * original.getWcs().linearizePixelToSky(source.getCoord(),

575 lsst.geom.degrees))

576 shape = shape.transform(transform.getLinear())

577 radius = source.get("ext_photometryKron_KronFlux_radius")

578 xc, yc = wcs.skyToPixel(source.getCoord())

579 for r, ct in [(radius, afwDisplay.BLUE), (radius*kfac, afwDisplay.CYAN), ]:

580 shape.scale(r/shape.getDeterminantRadius())

581 disp2.dot(shape, xc, yc, ctype=ct)

582 try:

583 self.assertFloatsAlmostEqual(source.get("ext_photometryKron_KronFlux_instFlux"),

584 forced.get("ext_photometryKron_KronFlux_instFlux"),

585 rtol=1.0e-3

586 )

587 self.assertFloatsAlmostEqual(source.get("ext_photometryKron_KronFlux_radius"),

588 scale*forced.get("ext_photometryKron_KronFlux_radius"),

589 rtol=1.0e-3

590 )

591 self.assertEqual(source.get("ext_photometryKron_KronFlux_flag"),

592 forced.get("ext_photometryKron_KronFlux_flag")

593 )

594 except Exception:

595 print(("Failed:", angle, scale, offset,

596 [(source.get(f), forced.get(f)) for f in

597 ("ext_photometryKron_KronFlux_instFlux",

598 "ext_photometryKron_KronFlux_radius",

599 "ext_photometryKron_KronFlux_flag"

600 )

601 ]))

602 raise

603

604

605class TestMemory(lsst.utils.tests.MemoryTestCase):

606 pass

607

608

609def setup_module(module):

610 lsst.utils.tests.init()

611 # Enable verbose mode for pytest. The verbose content only appears to the

612 # user on error or in the JUnit XML.

613 global verbose

614 verbose = 2

615

616

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

618 from argparse import ArgumentParser

619 parser = ArgumentParser()

620 parser.add_argument("--display", default=False, action="store_true", help="Activate display?")

621 parser.add_argument("--verbose", type=int, default=0, help="Verbosity level")

622 parser.add_argument('unittest_args', nargs='*')

623 args = parser.parse_args()

624 display = args.display

625 verbose = args.verbose

626 sys.argv[1:] = args.unittest_args

627 lsst.utils.tests.init()

628 unittest.main()

Coverage for tests/test_kron.py: 73%

391 statements