Coverage for tests/test_functors.py: 11%

499 statements  

« prev     ^ index     » next       coverage.py v7.11.0, created at 2025-11-04 02:26 -0800

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# (http://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 <http://www.gnu.org/licenses/>. 

21 

22import astropy.units as u 

23import copy 

24import functools 

25import numpy as np 

26import os 

27import pandas as pd 

28import unittest 

29import logging 

30 

31import lsst.daf.base as dafBase 

32import lsst.afw.geom as afwGeom 

33import lsst.afw.table as afwTable 

34import lsst.geom as geom 

35from lsst.sphgeom import HtmPixelization 

36import lsst.meas.base as measBase 

37import lsst.utils.tests 

38from lsst.pipe.base import InMemoryDatasetHandle 

39from lsst.pipe.tasks.functors import (CompositeFunctor, CustomFunctor, Column, RAColumn, 

40 DecColumn, Mag, MagDiff, Color, 

41 DeconvolvedMoments, SdssTraceSize, PsfSdssTraceSizeDiff, 

42 E1, E2, RadiusFromQuadrupole, 

43 HsmTraceSize, PsfHsmTraceSizeDiff, HsmFwhm, 

44 LocalPhotometry, LocalNanojansky, LocalNanojanskyErr, 

45 LocalDipoleMeanFlux, LocalDipoleMeanFluxErr, 

46 LocalDipoleDiffFlux, LocalDipoleDiffFluxErr, 

47 LocalWcs, ComputePixelScale, ConvertPixelToArcseconds, 

48 ConvertPixelSqToArcsecondsSq, 

49 ConvertDetectorAngleToPositionAngle, 

50 HtmIndex20, Ebv, MomentsIuuSky, MomentsIvvSky, MomentsIuvSky, 

51 SemimajorAxisFromMoments, SemiminorAxisFromMoments, 

52 PositionAngleFromMoments) 

53 

54ROOT = os.path.abspath(os.path.dirname(__file__)) 

55 

56 

57class FunctorTestCase(lsst.utils.tests.TestCase): 

58 

59 def getMultiIndexDataFrame(self, dataDict): 

60 """Create a simple test multi-index DataFrame.""" 

61 

62 simpleDF = pd.DataFrame(dataDict) 

63 dfFilterDSCombos = [] 

64 for ds in self.datasets: 

65 for band in self.bands: 

66 df = copy.copy(simpleDF) 

67 df.reindex(sorted(df.columns), axis=1) 

68 df['dataset'] = ds 

69 df['band'] = band 

70 df.columns = pd.MultiIndex.from_tuples( 

71 [(ds, band, c) for c in df.columns], 

72 names=('dataset', 'band', 'column')) 

73 dfFilterDSCombos.append(df) 

74 

75 df = functools.reduce(lambda d1, d2: d1.join(d2), dfFilterDSCombos) 

76 

77 return df 

78 

79 def getSimpleDataFrame(self, dataDict): 

80 return pd.DataFrame(dataDict) 

81 

82 def getDatasetHandle(self, df): 

83 lo, hi = HtmPixelization(7).universe().ranges()[0] 

84 value = np.random.randint(lo, hi) 

85 return InMemoryDatasetHandle(df, storageClass="DataFrame", dataId={"htm7": value}) 

86 

87 def setUp(self): 

88 np.random.seed(12345) 

89 self.datasets = ['forced_src', 'meas', 'ref'] 

90 self.bands = ['g', 'r'] 

91 self.columns = ['coord_ra', 'coord_dec'] 

92 self.nRecords = 5 

93 self.dataDict = { 

94 "coord_ra": [3.77654137, 3.77643059, 3.77621148, 3.77611944, 3.77610396], 

95 "coord_dec": [0.01127624, 0.01127787, 0.01127543, 0.01127543, 0.01127543]} 

96 

97 def _funcVal(self, functor, df): 

98 self.assertIsInstance(functor.name, str) 

99 self.assertIsInstance(functor.shortname, str) 

100 

101 handle = self.getDatasetHandle(df) 

102 

103 val = functor(df) 

104 val2 = functor(handle) 

105 self.assertTrue((val == val2).all()) 

106 self.assertIsInstance(val, pd.Series) 

107 

108 val = functor(df, dropna=True) 

109 val2 = functor(handle, dropna=True) 

110 self.assertTrue((val == val2).all()) 

111 self.assertEqual(val.isnull().sum(), 0) 

112 

113 return val 

114 

115 def _differenceVal(self, functor, df1, df2): 

116 self.assertIsInstance(functor.name, str) 

117 self.assertIsInstance(functor.shortname, str) 

118 

119 handle1 = self.getDatasetHandle(df1) 

120 handle2 = self.getDatasetHandle(df2) 

121 

122 val = functor.difference(df1, df2) 

123 val2 = functor.difference(handle1, handle2) 

124 self.assertTrue(val.equals(val2)) 

125 self.assertIsInstance(val, pd.Series) 

126 

127 val = functor.difference(df1, df2, dropna=True) 

128 val2 = functor.difference(handle1, handle2, dropna=True) 

129 self.assertTrue(val.equals(val2)) 

130 self.assertEqual(val.isnull().sum(), 0) 

131 

132 val1 = self._funcVal(functor, df1) 

133 val2 = self._funcVal(functor, df2) 

134 

135 self.assertTrue(np.allclose(val, val1 - val2)) 

136 

137 return val 

138 

139 def testColumn(self): 

140 self.columns.append("base_FootprintArea_value") 

141 self.dataDict["base_FootprintArea_value"] = \ 

142 np.full(self.nRecords, 1) 

143 df = self.getMultiIndexDataFrame(self.dataDict) 

144 func = Column('base_FootprintArea_value', filt='g') 

145 self._funcVal(func, df) 

146 

147 df = self.getSimpleDataFrame(self.dataDict) 

148 func = Column('base_FootprintArea_value') 

149 self._funcVal(func, df) 

150 

151 def testCustom(self): 

152 self.columns.append("base_FootprintArea_value") 

153 self.dataDict["base_FootprintArea_value"] = \ 

154 np.random.rand(self.nRecords) 

155 df = self.getMultiIndexDataFrame(self.dataDict) 

156 func = CustomFunctor('2*base_FootprintArea_value', filt='g') 

157 val = self._funcVal(func, df) 

158 

159 func2 = Column('base_FootprintArea_value', filt='g') 

160 

161 np.allclose(val.values, 2*func2(df).values, atol=1e-13, rtol=0) 

162 

163 df = self.getSimpleDataFrame(self.dataDict) 

164 func = CustomFunctor('2 * base_FootprintArea_value') 

165 val = self._funcVal(func, df) 

166 func2 = Column('base_FootprintArea_value') 

167 

168 np.allclose(val.values, 2*func2(df).values, atol=1e-13, rtol=0) 

169 

170 def testCoords(self): 

171 df = self.getMultiIndexDataFrame(self.dataDict) 

172 ra = self._funcVal(RAColumn(), df) 

173 dec = self._funcVal(DecColumn(), df) 

174 

175 columnDict = {'dataset': 'ref', 'band': 'g', 

176 'column': ['coord_ra', 'coord_dec']} 

177 

178 handle = InMemoryDatasetHandle(df, storageClass="DataFrame") 

179 dfSub = handle.get(parameters={"columns": columnDict}) 

180 self._dropLevels(dfSub) 

181 

182 coords = dfSub / np.pi * 180. 

183 

184 self.assertTrue(np.allclose(ra, coords[('ref', 'g', 'coord_ra')], atol=1e-13, rtol=0)) 

185 self.assertTrue(np.allclose(dec, coords[('ref', 'g', 'coord_dec')], atol=1e-13, rtol=0)) 

186 

187 # single-level column index table 

188 df = self.getSimpleDataFrame(self.dataDict) 

189 ra = self._funcVal(RAColumn(), df) 

190 dec = self._funcVal(DecColumn(), df) 

191 

192 handle = InMemoryDatasetHandle(df, storageClass="DataFrame") 

193 dfSub = handle.get(parameters={"columns": ['coord_ra', 'coord_dec']}) 

194 coords = dfSub / np.pi * 180. 

195 

196 self.assertTrue(np.allclose(ra, coords['coord_ra'], atol=1e-13, rtol=0)) 

197 self.assertTrue(np.allclose(dec, coords['coord_dec'], atol=1e-13, rtol=0)) 

198 

199 def testMag(self): 

200 self.columns.extend(["base_PsfFlux_instFlux", "base_PsfFlux_instFluxErr"]) 

201 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000) 

202 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10) 

203 df = self.getMultiIndexDataFrame(self.dataDict) 

204 # Change one dataset filter combinations value. 

205 df[("meas", "g", "base_PsfFlux_instFlux")] -= 1 

206 

207 fluxName = 'base_PsfFlux' 

208 

209 # Check that things work when you provide dataset explicitly 

210 for dataset in ['forced_src', 'meas']: 

211 psfMag_G = self._funcVal(Mag(fluxName, dataset=dataset, 

212 filt='g'), 

213 df) 

214 psfMag_R = self._funcVal(Mag(fluxName, dataset=dataset, 

215 filt='r'), 

216 df) 

217 

218 psfColor_GR = self._funcVal(Color(fluxName, 'g', 'r', 

219 dataset=dataset), 

220 df) 

221 

222 self.assertTrue(np.allclose((psfMag_G - psfMag_R).dropna(), psfColor_GR, rtol=0, atol=1e-13)) 

223 

224 # Check that behavior as expected when dataset not provided; 

225 # that is, that the color comes from forced and default Mag is meas 

226 psfMag_G = self._funcVal(Mag(fluxName, filt='g'), df) 

227 psfMag_R = self._funcVal(Mag(fluxName, filt='r'), df) 

228 

229 psfColor_GR = self._funcVal(Color(fluxName, 'g', 'r'), df) 

230 

231 # These should *not* be equal. 

232 self.assertFalse(np.allclose((psfMag_G - psfMag_R).dropna(), psfColor_GR)) 

233 

234 def testMagDiff(self): 

235 self.columns.extend(["base_PsfFlux_instFlux", "base_PsfFlux_instFluxErr", 

236 "modelfit_CModel_instFlux", "modelfit_CModel_instFluxErr"]) 

237 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000) 

238 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10) 

239 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1000) 

240 self.dataDict["modelfit_CModel_instFluxErr"] = np.full(self.nRecords, 10) 

241 df = self.getMultiIndexDataFrame(self.dataDict) 

242 

243 for filt in self.bands: 

244 filt = 'g' 

245 val = self._funcVal(MagDiff('base_PsfFlux', 'modelfit_CModel', filt=filt), df) 

246 

247 mag1 = self._funcVal(Mag('modelfit_CModel', filt=filt), df) 

248 mag2 = self._funcVal(Mag('base_PsfFlux', filt=filt), df) 

249 self.assertTrue(np.allclose((mag2 - mag1).dropna(), val, rtol=0, atol=1e-13)) 

250 

251 def testDifference(self): 

252 """Test .difference method using MagDiff as the example. 

253 """ 

254 self.columns.extend(["base_PsfFlux_instFlux", "base_PsfFlux_instFluxErr", 

255 "modelfit_CModel_instFlux", "modelfit_CModel_instFluxErr"]) 

256 

257 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000) 

258 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1000) 

259 df1 = self.getMultiIndexDataFrame(self.dataDict) 

260 

261 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 999) 

262 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 999) 

263 df2 = self.getMultiIndexDataFrame(self.dataDict) 

264 

265 magDiff = MagDiff('base_PsfFlux', 'modelfit_CModel', filt='g') 

266 

267 # Asserts that differences computed properly 

268 self._differenceVal(magDiff, df1, df2) 

269 

270 def testPixelScale(self): 

271 # Test that the pixel scale and pix->arcsec calculations perform as 

272 # expected. 

273 pass 

274 

275 def testShape(self): 

276 data = { 

277 "x": np.array([-0.3, 0.4, 0.7, -0.9, 1.4, -5.3]), 

278 "y": np.array([1.5, -0.7, -1.9, 2.8, -1.4, 0.01]), 

279 "rho": np.array([-0.9, 0.4, -0.7, 0., 0.3, -0.99]), 

280 } 

281 data["xx"] = data["x"]**2 

282 data["yy"] = data["y"]**2 

283 data["xy"] = data["x"]*data["y"]*data["rho"] 

284 

285 args = ("xx", "xy", "yy") 

286 functor_e1, functor_e2, functor_quadrupole = E1(*args), E2(*args), RadiusFromQuadrupole(*args) 

287 

288 xx_plus_yy = data["xx"] + data["yy"] 

289 data = pd.DataFrame(data) 

290 

291 np.testing.assert_allclose( 

292 functor_e1(data), 

293 ((data["xx"] - data["yy"])/xx_plus_yy).astype(np.float32), 

294 rtol=1e-12, atol=1e-12, 

295 ) 

296 np.testing.assert_allclose( 

297 functor_e2(data), 

298 (2.0*data["xy"]/xx_plus_yy).astype(np.float32), 

299 rtol=1e-12, atol=1e-12, 

300 ) 

301 np.testing.assert_allclose( 

302 functor_quadrupole(data), 

303 ((data["xx"]*data["yy"] - data["xy"]**2)**0.25).astype(np.float32), 

304 rtol=1e-12, atol=1e-12, 

305 ) 

306 

307 def testOther(self): 

308 self.columns.extend(["ext_shapeHSM_HsmSourceMoments_xx", "ext_shapeHSM_HsmSourceMoments_yy", 

309 "base_SdssShape_xx", "base_SdssShape_yy", 

310 "ext_shapeHSM_HsmPsfMoments_xx", "ext_shapeHSM_HsmPsfMoments_yy", 

311 "base_SdssShape_psf_xx", "base_SdssShape_psf_yy"]) 

312 self.dataDict["ext_shapeHSM_HsmSourceMoments_xx"] = np.full(self.nRecords, 1 / np.sqrt(2)) 

313 self.dataDict["ext_shapeHSM_HsmSourceMoments_yy"] = np.full(self.nRecords, 1 / np.sqrt(2)) 

314 self.dataDict["base_SdssShape_xx"] = np.full(self.nRecords, 1 / np.sqrt(2)) 

315 self.dataDict["base_SdssShape_yy"] = np.full(self.nRecords, 1 / np.sqrt(2)) 

316 self.dataDict["ext_shapeHSM_HsmPsfMoments_xx"] = np.full(self.nRecords, 1 / np.sqrt(2)) 

317 self.dataDict["ext_shapeHSM_HsmPsfMoments_yy"] = np.full(self.nRecords, 1 / np.sqrt(2)) 

318 self.dataDict["base_SdssShape_psf_xx"] = np.full(self.nRecords, 1) 

319 self.dataDict["base_SdssShape_psf_yy"] = np.full(self.nRecords, 1) 

320 df = self.getMultiIndexDataFrame(self.dataDict) 

321 # Covering the code is better than nothing 

322 for filt in self.bands: 

323 for Func in [DeconvolvedMoments, 

324 SdssTraceSize, 

325 PsfSdssTraceSizeDiff, 

326 HsmTraceSize, PsfHsmTraceSizeDiff, HsmFwhm]: 

327 _ = self._funcVal(Func(filt=filt), df) 

328 

329 def _compositeFuncVal(self, functor, df): 

330 self.assertIsInstance(functor, CompositeFunctor) 

331 

332 handle = self.getDatasetHandle(df) 

333 

334 fdf1 = functor(df) 

335 fdf2 = functor(handle) 

336 self.assertTrue(fdf1.equals(fdf2)) 

337 

338 self.assertIsInstance(fdf1, pd.DataFrame) 

339 self.assertTrue(np.all([k in fdf1.columns for k in functor.funcDict.keys()])) 

340 

341 fdf1 = functor(df, dropna=True) 

342 fdf2 = functor(handle, dropna=True) 

343 self.assertTrue(fdf1.equals(fdf2)) 

344 

345 # Check that there are no nulls 

346 self.assertFalse(fdf1.isnull().any(axis=None)) 

347 

348 return fdf1 

349 

350 def _compositeDifferenceVal(self, functor, df1, df2): 

351 self.assertIsInstance(functor, CompositeFunctor) 

352 

353 handle1 = self.getDatasetHandle(df1) 

354 handle2 = self.getDatasetHandle(df2) 

355 

356 fdf1 = functor.difference(df1, df2) 

357 fdf2 = functor.difference(handle1, handle2) 

358 self.assertTrue(fdf1.equals(fdf2)) 

359 

360 self.assertIsInstance(fdf1, pd.DataFrame) 

361 self.assertTrue(np.all([k in fdf1.columns for k in functor.funcDict.keys()])) 

362 

363 fdf1 = functor.difference(df1, df2, dropna=True) 

364 fdf2 = functor.difference(handle1, handle2, dropna=True) 

365 self.assertTrue(fdf1.equals(fdf2)) 

366 

367 # Check that there are no nulls 

368 self.assertFalse(fdf1.isnull().any(axis=None)) 

369 

370 df1_functored = functor(df1) 

371 df2_functored = functor(df2) 

372 

373 self.assertTrue(np.allclose(fdf1.values, df1_functored.values - df2_functored.values)) 

374 

375 return fdf1 

376 

377 def testComposite(self): 

378 self.columns.extend(["modelfit_CModel_instFlux", "base_PsfFlux_instFlux"]) 

379 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1) 

380 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1) 

381 

382 df = self.getMultiIndexDataFrame(self.dataDict) 

383 # Modify r band value slightly. 

384 df[("meas", "r", "base_PsfFlux_instFlux")] -= 0.1 

385 

386 filt = 'g' 

387 funcDict = {'psfMag_ref': Mag('base_PsfFlux', dataset='ref'), 

388 'ra': RAColumn(), 

389 'dec': DecColumn(), 

390 'psfMag': Mag('base_PsfFlux', filt=filt), 

391 'cmodel_magDiff': MagDiff('base_PsfFlux', 

392 'modelfit_CModel', filt=filt)} 

393 func = CompositeFunctor(funcDict) 

394 fdf1 = self._compositeFuncVal(func, df) 

395 

396 # Repeat same, but define filter globally instead of individually 

397 funcDict2 = {'psfMag_ref': Mag('base_PsfFlux', dataset='ref'), 

398 'ra': RAColumn(), 

399 'dec': DecColumn(), 

400 'psfMag': Mag('base_PsfFlux'), 

401 'cmodel_magDiff': MagDiff('base_PsfFlux', 

402 'modelfit_CModel')} 

403 

404 func2 = CompositeFunctor(funcDict2, filt=filt) 

405 fdf2 = self._compositeFuncVal(func2, df) 

406 self.assertTrue(fdf1.equals(fdf2)) 

407 

408 func2.filt = 'r' 

409 fdf3 = self._compositeFuncVal(func2, df) 

410 # Because we modified the R filter this should fail. 

411 self.assertFalse(fdf2.equals(fdf3)) 

412 

413 # Make sure things work with passing list instead of dict 

414 funcs = [Mag('base_PsfFlux', dataset='ref'), 

415 RAColumn(), 

416 DecColumn(), 

417 Mag('base_PsfFlux', filt=filt), 

418 MagDiff('base_PsfFlux', 'modelfit_CModel', filt=filt)] 

419 

420 _ = self._compositeFuncVal(CompositeFunctor(funcs), df) 

421 

422 def testCompositeSimple(self): 

423 """Test single-level composite functor for functionality 

424 """ 

425 self.columns.extend(["modelfit_CModel_instFlux", "base_PsfFlux_instFlux"]) 

426 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1) 

427 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1) 

428 

429 df = self.getSimpleDataFrame(self.dataDict) 

430 funcDict = {'ra': RAColumn(), 

431 'dec': DecColumn(), 

432 'psfMag': Mag('base_PsfFlux'), 

433 'cmodel_magDiff': MagDiff('base_PsfFlux', 

434 'modelfit_CModel')} 

435 func = CompositeFunctor(funcDict) 

436 _ = self._compositeFuncVal(func, df) 

437 

438 def testCompositeColor(self): 

439 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000) 

440 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10) 

441 df = self.getMultiIndexDataFrame(self.dataDict) 

442 funcDict = {'a': Mag('base_PsfFlux', dataset='meas', filt='g'), 

443 'b': Mag('base_PsfFlux', dataset='forced_src', filt='g'), 

444 'c': Color('base_PsfFlux', 'g', 'r')} 

445 # Covering the code is better than nothing 

446 _ = self._compositeFuncVal(CompositeFunctor(funcDict), df) 

447 

448 def testCompositeDifference(self): 

449 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000) 

450 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10) 

451 df1 = self.getMultiIndexDataFrame(self.dataDict) 

452 

453 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 999) 

454 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 9) 

455 df2 = self.getMultiIndexDataFrame(self.dataDict) 

456 

457 funcDict = {'a': Mag('base_PsfFlux', dataset='meas', filt='g'), 

458 'b': Mag('base_PsfFlux', dataset='forced_src', filt='g'), 

459 'c': Color('base_PsfFlux', 'g', 'r')} 

460 # Covering the code is better than nothing 

461 _ = self._compositeDifferenceVal(CompositeFunctor(funcDict), df1, df2) 

462 

463 def testCompositeFail(self): 

464 """Test a composite functor where one of the functors should be junk. 

465 """ 

466 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000) 

467 df = self.getMultiIndexDataFrame(self.dataDict) 

468 

469 funcDict = {'good': Column("base_PsfFlux_instFlux"), 

470 'bad': Column('not_a_column')} 

471 

472 with self.assertLogs(level=logging.ERROR) as cm: 

473 _ = self._compositeFuncVal(CompositeFunctor(funcDict), df) 

474 self.assertIn("Exception in CompositeFunctor (funcs: ['good', 'bad'])", cm.output[0]) 

475 

476 def testLocalPhotometry(self): 

477 """Test the local photometry functors. 

478 """ 

479 flux = 1000 

480 fluxErr = 10 

481 calib = 10 

482 calibErr = 0.0 

483 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, flux) 

484 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 

485 fluxErr) 

486 self.dataDict["base_LocalPhotoCalib"] = np.full(self.nRecords, calib) 

487 

488 df = self.getMultiIndexDataFrame(self.dataDict) 

489 func = LocalPhotometry("base_PsfFlux_instFlux", 

490 "base_PsfFlux_instFluxErr", 

491 "base_LocalPhotoCalib") 

492 

493 nanoJansky = func.instFluxToNanojansky( 

494 df[("meas", "g", "base_PsfFlux_instFlux")], 

495 df[("meas", "g", "base_LocalPhotoCalib")]) 

496 mag = func.instFluxToMagnitude( 

497 df[("meas", "g", "base_PsfFlux_instFlux")], 

498 df[("meas", "g", "base_LocalPhotoCalib")]) 

499 nanoJanskyErr = func.instFluxErrToNanojanskyErr( 

500 df[("meas", "g", "base_PsfFlux_instFlux")], 

501 df[("meas", "g", "base_PsfFlux_instFluxErr")], 

502 df[("meas", "g", "base_LocalPhotoCalib")]) 

503 magErr = func.instFluxErrToMagnitudeErr( 

504 df[("meas", "g", "base_PsfFlux_instFlux")], 

505 df[("meas", "g", "base_PsfFlux_instFluxErr")], 

506 df[("meas", "g", "base_LocalPhotoCalib")]) 

507 

508 self.assertTrue(np.allclose(nanoJansky.values, 

509 flux * calib, 

510 atol=0, 

511 rtol=1e-7)) 

512 self.assertTrue(np.allclose(mag.values, 

513 (flux * calib * u.nJy).to_value(u.ABmag), 

514 atol=0, 

515 rtol=1e-7)) 

516 self.assertTrue(np.allclose(nanoJanskyErr.values, 

517 np.hypot(fluxErr * calib, flux * calibErr), 

518 atol=0, 

519 rtol=1e-7)) 

520 self.assertTrue(np.allclose( 

521 magErr.values, 

522 2.5 / np.log(10) * nanoJanskyErr.values / nanoJansky.values, 

523 atol=0, 

524 rtol=1e-7)) 

525 

526 # Test functors against the values computed above. 

527 self._testLocalPhotometryFunctors(LocalNanojansky, 

528 df, 

529 nanoJansky) 

530 self._testLocalPhotometryFunctors(LocalNanojanskyErr, 

531 df, 

532 nanoJanskyErr) 

533 

534 def _testLocalPhotometryFunctors(self, functor, df, testValues): 

535 func = functor("base_PsfFlux_instFlux", 

536 "base_PsfFlux_instFluxErr", 

537 "base_LocalPhotoCalib") 

538 val = self._funcVal(func, df) 

539 self.assertTrue(np.allclose(testValues.values, 

540 val.values, 

541 atol=0, 

542 rtol=1e-7)) 

543 

544 def testDipPhotometry(self): 

545 """Test calibrated flux calculations for dipoles.""" 

546 fluxNeg = -100 

547 fluxPos = 101 

548 fluxErr = 10 

549 calib = 10 

550 calibErr = 0.0 

551 

552 # compute expected values. 

553 absMean = 0.5*(fluxPos - fluxNeg)*calib 

554 absDiff = (fluxNeg + fluxPos)*calib 

555 absMeanErr = 0.5*np.sqrt(2*(fluxErr*calib)**2 

556 + ((fluxPos - fluxNeg)*calibErr)**2) 

557 absDiffErr = np.sqrt(2*(fluxErr*calib)**2 

558 + ((fluxPos + fluxNeg)*calibErr)**2) 

559 

560 self.dataDict["ip_diffim_DipoleFluxNeg_instFlux"] = np.full(self.nRecords, fluxNeg) 

561 self.dataDict["ip_diffim_DipoleFluxNeg_instFluxErr"] = np.full(self.nRecords, fluxErr) 

562 self.dataDict["ip_diffim_DipoleFluxPos_instFlux"] = np.full(self.nRecords, fluxPos) 

563 self.dataDict["ip_diffim_DipoleFluxPos_instFluxErr"] = np.full(self.nRecords, fluxErr) 

564 self.dataDict["base_LocalPhotoCalib"] = np.full(self.nRecords, calib) 

565 

566 df = self.getMultiIndexDataFrame(self.dataDict) 

567 func = LocalDipoleMeanFlux("ip_diffim_DipoleFluxPos_instFlux", 

568 "ip_diffim_DipoleFluxNeg_instFlux", 

569 "ip_diffim_DipoleFluxPos_instFluxErr", 

570 "ip_diffim_DipoleFluxNeg_instFluxErr", 

571 "base_LocalPhotoCalib") 

572 val = self._funcVal(func, df) 

573 self.assertFloatsAlmostEqual(val.values, 

574 absMean, 

575 atol=1e-13, 

576 rtol=0) 

577 

578 func = LocalDipoleMeanFluxErr("ip_diffim_DipoleFluxPos_instFlux", 

579 "ip_diffim_DipoleFluxNeg_instFlux", 

580 "ip_diffim_DipoleFluxPos_instFluxErr", 

581 "ip_diffim_DipoleFluxNeg_instFluxErr", 

582 "base_LocalPhotoCalib") 

583 val = self._funcVal(func, df) 

584 self.assertFloatsAlmostEqual(val.values, 

585 absMeanErr, 

586 atol=1e-13, 

587 rtol=0) 

588 

589 func = LocalDipoleDiffFlux("ip_diffim_DipoleFluxPos_instFlux", 

590 "ip_diffim_DipoleFluxNeg_instFlux", 

591 "ip_diffim_DipoleFluxPos_instFluxErr", 

592 "ip_diffim_DipoleFluxNeg_instFluxErr", 

593 "base_LocalPhotoCalib") 

594 val = self._funcVal(func, df) 

595 self.assertFloatsAlmostEqual(val.values, 

596 absDiff, 

597 atol=1e-13, 

598 rtol=0) 

599 

600 func = LocalDipoleDiffFluxErr("ip_diffim_DipoleFluxPos_instFlux", 

601 "ip_diffim_DipoleFluxNeg_instFlux", 

602 "ip_diffim_DipoleFluxPos_instFluxErr", 

603 "ip_diffim_DipoleFluxNeg_instFluxErr", 

604 "base_LocalPhotoCalib") 

605 val = self._funcVal(func, df) 

606 self.assertFloatsAlmostEqual(val.values, 

607 absDiffErr, 

608 atol=1e-13, 

609 rtol=0) 

610 

611 def testComputePositionAngle(self, offset=0.0001): 

612 """Test computation of position angle from (RA1, Dec1) to (RA2, Dec2) 

613 

614 offset : `float` 

615 Arc length of the offset vector to set up test points. [radian] 

616 """ 

617 

618 d = offset 

619 columns = ("ra1", "dec1", "ra2", "dec2", "expected") 

620 position_angle_test_values = [ 

621 # Get 0, 0 right 

622 (0, 0, d, 0, np.pi/2), 

623 (0, 0, 0, d, 0), 

624 (0, 0, -d, 0, -np.pi/2), 

625 (0, 0, 0, -d, np.pi), 

626 # Make sure we get wrap-around to 0, 0 right 

627 (2*np.pi, 0, d, 0, np.pi/2), 

628 (2*np.pi, 0, 0, d, 0), 

629 (2*np.pi, 0, -d, 0, -np.pi/2), 

630 (2*np.pi, 0, 0, -d, np.pi), 

631 (+0.0015, 0, 2*np.pi - 0.05, 0, -np.pi/2), 

632 # Get another somewhat arbitrary location right [these are in rad] 

633 # It's not really important to rescale d by 1/cos(dec) 

634 # because we're just looking at orientation of vector 

635 # but foreshadowing to the poles... 

636 (0.0015, 1, 0.0015 + d*np.cos(-1), 1, np.pi/2), 

637 (0.0015, 1, 0.0015, 1 + d, 0), 

638 (0.0015, 1, 0.0015 - d*np.cos(-1), 1, - np.pi/2), 

639 (0.0015, 1, 0.0015, 1 - d, np.pi), 

640 # Negative dec 

641 (0.0015, -1, 0.0015 + d*np.cos(-1), -1, np.pi/2), 

642 (0.0015, -1, 0.0015, -1 + d, 0), 

643 (0.0015, -1, 0.0015 - d*np.cos(-1), -1, - np.pi/2), 

644 (0.0015, -1, 0.0015, -1 - d, np.pi), 

645 # Make sure we get wrap-around on that right 

646 (2*np.pi + 0.0015, 1, 2*np.pi + 0.0015 + d*np.cos(1), 1, np.pi/2), 

647 (2*np.pi + 0.0015, 1, 2*np.pi + 0.0015, 1 + d, 0), 

648 (2*np.pi + 0.0015, 1, 0.0015 - d*np.cos(1), 1, - np.pi/2), 

649 (2*np.pi + 0.0015, 1, 0.0015, 1 - d, np.pi), 

650 # Get relative wrap-around right 

651 (2*np.pi + 0.0015, 1, 0.0015 + d*np.cos(1), 1, np.pi/2), 

652 (0.0015, 1, 2*np.pi + 0.0015, 1 + d, 0), 

653 (0.0015, 1, 2*np.pi + 0.0015 - d*np.cos(1), 1, - np.pi/2), 

654 (2*np.pi + 0.0015, 1, 0.0015, 1 - d, np.pi), 

655 # Get either side of RA=0 right 

656 (+ d*np.cos(1) / 2, 1, 2*np.pi - d*np.cos(1) / 2, 1, - np.pi/2), 

657 (+ d*np.cos(1) / 2, 1, - d*np.cos(1) / 2, 1, -np.pi/2), 

658 (2*np.pi + d*np.cos(1) / 2, 1, - d*np.cos(1) / 2, 1, -np.pi/2), 

659 (-d*np.cos(1) / 2, 1, +0.0015, 1, np.pi/2), 

660 (0.0015, 1, 0.0015, 1 + d, 0), 

661 (0.0015, 1, 0.0015 - d*np.cos(1), 1, - np.pi/2), 

662 (0.0015, 1, 0.0015, 1 - d, np.pi), 

663 # Try it near the poles 

664 (0, np.pi/2, 0, np.pi/2 - d, np.pi), 

665 (0, np.pi/2 - d, 0, np.pi/2, 0), 

666 (0, -np.pi/2, 0, -np.pi/2 + d, 0), 

667 (0, -np.pi/2 + d, 0, -np.pi/2, np.pi), 

668 ] 

669 

670 df = pd.DataFrame(position_angle_test_values, columns=columns) 

671 

672 cd_matrix = np.array([[1, 0], [0, -1]]) # Doesn't matter because we don't use it. 

673 local_wcs = LocalWcs(*cd_matrix.flatten()) 

674 pa = local_wcs.computePositionAngle(df["ra1"], df["dec1"], df["ra2"], df["dec2"]) 

675 expected = df["expected"] 

676 tolerance_deg = 0.05 # degrees 

677 tolerance_rad = np.deg2rad(tolerance_deg) 

678 

679 # Use SphereGeom to handle wrap-around separations correctly. 

680 diff = [ 

681 geom.Angle(o, geom.radians).separation(geom.Angle(e, geom.radians)).asRadians() 

682 for o, e in zip(pa, expected) 

683 ] 

684 

685 np.testing.assert_allclose(diff, 0, rtol=0, atol=tolerance_rad) 

686 

687 # Test position angle 

688 def testConvertDetectorAngleToPositionAngle(self): 

689 """Test conversion of position angle in detector degrees to position angle on sky 

690 

691 There is overlap with testConvertPixelToArcseconds 

692 But we also test additional rotation angles so this is separate. 

693 """ 

694 dipoleSep = 10 

695 ixx = 10 

696 testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2)) 

697 # testConvertPixelToArcSecond uses a meas_base LocalWcsPlugin 

698 # but we're using a simple WCS as our example, so this doesn't really matter 

699 # and we'll just use the WCS directly 

700 for dec in np.linspace(-90, 90, 10): 

701 for theta in (-45, 0, 90): 

702 for x, y in zip(np.random.uniform(2 * 1109.99981456774, size=10), 

703 np.random.uniform(2 * 560.018167811613, size=10)): 

704 wcs = self._makeWcs(dec=dec, theta=theta) 

705 cd_matrix = wcs.getCdMatrix() 

706 

707 self.dataDict["dipoleSep"] = np.full(self.nRecords, dipoleSep) 

708 self.dataDict["ixx"] = np.full(self.nRecords, ixx) 

709 self.dataDict["slot_Centroid_x"] = np.full(self.nRecords, x) 

710 self.dataDict["slot_Centroid_y"] = np.full(self.nRecords, y) 

711 self.dataDict["someCentroid_x"] = x + testPixelDeltas[:, 0] 

712 self.dataDict["someCentroid_y"] = y + testPixelDeltas[:, 1] 

713 self.dataDict["orientation"] = np.arctan2( 

714 self.dataDict["someCentroid_y"] - self.dataDict["slot_Centroid_y"], 

715 self.dataDict["someCentroid_x"] - self.dataDict["slot_Centroid_x"], 

716 ) 

717 

718 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(self.nRecords, 

719 cd_matrix[0, 0]) 

720 self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(self.nRecords, 

721 cd_matrix[0, 1]) 

722 self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(self.nRecords, 

723 cd_matrix[1, 0]) 

724 self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(self.nRecords, 

725 cd_matrix[1, 1]) 

726 df = self.getMultiIndexDataFrame(self.dataDict) 

727 

728 # Test detector angle to position angle conversion 

729 func = ConvertDetectorAngleToPositionAngle( 

730 "orientation", 

731 "base_LocalWcs_CDMatrix_1_1", 

732 "base_LocalWcs_CDMatrix_1_2", 

733 "base_LocalWcs_CDMatrix_2_1", 

734 "base_LocalWcs_CDMatrix_2_2" 

735 ) 

736 val = self._funcVal(func, df) 

737 

738 delta_ra, delta_dec = func.computeDeltaRaDec( 

739 self.dataDict["someCentroid_x"] - self.dataDict["slot_Centroid_x"], 

740 self.dataDict["someCentroid_y"] - self.dataDict["slot_Centroid_y"], 

741 self.dataDict["base_LocalWcs_CDMatrix_1_1"], 

742 self.dataDict["base_LocalWcs_CDMatrix_1_2"], 

743 self.dataDict["base_LocalWcs_CDMatrix_2_1"], 

744 self.dataDict["base_LocalWcs_CDMatrix_2_2"], 

745 ) 

746 

747 dx = np.cos(0) * np.tan(delta_dec) - np.sin(0) * np.cos(delta_ra) 

748 dy = np.sin(delta_ra) 

749 comparison_pa = np.arctan2(dy, dx) 

750 comparison_pa = np.rad2deg(comparison_pa) 

751 

752 coord_diff = [] 

753 for v, c in zip(val.values, comparison_pa): 

754 observed_angle = geom.Angle(v, geom.degrees) 

755 expected_angle = geom.Angle(c, geom.degrees) 

756 diff = observed_angle.separation(expected_angle).asRadians() 

757 coord_diff.append(diff) 

758 

759 np.testing.assert_allclose(coord_diff, 0, rtol=0, atol=5e-6) 

760 

761 def testConvertPixelToArcseconds(self): 

762 """Test calculations of the pixel scale, conversions of pixel to 

763 arcseconds. 

764 """ 

765 dipoleSep = 10 

766 ixx = 10 

767 testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2)) 

768 localWcsPlugin = measBase.EvaluateLocalWcsPlugin( 

769 None, 

770 "base_LocalWcs", 

771 afwTable.SourceTable.makeMinimalSchema(), 

772 None) 

773 for dec in np.linspace(-90, 90, 10): 

774 for x, y in zip(np.random.uniform(2 * 1109.99981456774, size=10), 

775 np.random.uniform(2 * 560.018167811613, size=10)): 

776 center = geom.Point2D(x, y) 

777 wcs = self._makeWcs(dec) 

778 skyOrigin = wcs.pixelToSky(center) 

779 

780 linAffMatrix = localWcsPlugin.makeLocalTransformMatrix(wcs, 

781 center) 

782 self.dataDict["dipoleSep"] = np.full(self.nRecords, dipoleSep) 

783 self.dataDict["ixx"] = np.full(self.nRecords, ixx) 

784 self.dataDict["slot_Centroid_x"] = np.full(self.nRecords, x) 

785 self.dataDict["slot_Centroid_y"] = np.full(self.nRecords, y) 

786 self.dataDict["someCentroid_x"] = x + testPixelDeltas[:, 0] 

787 self.dataDict["someCentroid_y"] = y + testPixelDeltas[:, 1] 

788 

789 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(self.nRecords, 

790 linAffMatrix[0, 0]) 

791 self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(self.nRecords, 

792 linAffMatrix[0, 1]) 

793 self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(self.nRecords, 

794 linAffMatrix[1, 0]) 

795 self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(self.nRecords, 

796 linAffMatrix[1, 1]) 

797 df = self.getMultiIndexDataFrame(self.dataDict) 

798 func = LocalWcs("base_LocalWcs_CDMatrix_1_1", 

799 "base_LocalWcs_CDMatrix_1_2", 

800 "base_LocalWcs_CDMatrix_2_1", 

801 "base_LocalWcs_CDMatrix_2_2") 

802 

803 # Exercise the full set of functions in LocalWcs. 

804 sepRadians = func.getSkySeparationFromPixel( 

805 df[("meas", "g", "someCentroid_x")] - df[("meas", "g", "slot_Centroid_x")], 

806 df[("meas", "g", "someCentroid_y")] - df[("meas", "g", "slot_Centroid_y")], 

807 0.0, 

808 0.0, 

809 df[("meas", "g", "base_LocalWcs_CDMatrix_1_1")], 

810 df[("meas", "g", "base_LocalWcs_CDMatrix_1_2")], 

811 df[("meas", "g", "base_LocalWcs_CDMatrix_2_1")], 

812 df[("meas", "g", "base_LocalWcs_CDMatrix_2_2")]) 

813 

814 # Test functor values against afw SkyWcs computations. 

815 for centX, centY, sep in zip(testPixelDeltas[:, 0], 

816 testPixelDeltas[:, 1], 

817 sepRadians.values): 

818 afwSepRadians = skyOrigin.separation( 

819 wcs.pixelToSky(x + centX, y + centY)).asRadians() 

820 self.assertAlmostEqual(1 - sep / afwSepRadians, 0, places=5) 

821 

822 # Test the pixel scale computation. 

823 func = ComputePixelScale("base_LocalWcs_CDMatrix_1_1", 

824 "base_LocalWcs_CDMatrix_1_2", 

825 "base_LocalWcs_CDMatrix_2_1", 

826 "base_LocalWcs_CDMatrix_2_2") 

827 pixelScale = self._funcVal(func, df) 

828 self.assertTrue(np.allclose( 

829 wcs.getPixelScale(center).asArcseconds(), 

830 pixelScale.values, 

831 rtol=1e-6, 

832 atol=0)) 

833 

834 # Test pixel -> arcsec conversion. 

835 func = ConvertPixelToArcseconds("dipoleSep", 

836 "base_LocalWcs_CDMatrix_1_1", 

837 "base_LocalWcs_CDMatrix_1_2", 

838 "base_LocalWcs_CDMatrix_2_1", 

839 "base_LocalWcs_CDMatrix_2_2") 

840 val = self._funcVal(func, df) 

841 self.assertTrue(np.allclose(pixelScale.values * dipoleSep, 

842 val.values, 

843 atol=1e-16, 

844 rtol=1e-16)) 

845 

846 # Test pixel^2 -> arcsec^2 conversion. 

847 func = ConvertPixelSqToArcsecondsSq("ixx", 

848 "base_LocalWcs_CDMatrix_1_1", 

849 "base_LocalWcs_CDMatrix_1_2", 

850 "base_LocalWcs_CDMatrix_2_1", 

851 "base_LocalWcs_CDMatrix_2_2") 

852 val = self._funcVal(func, df) 

853 self.assertTrue(np.allclose(pixelScale.values ** 2 * dipoleSep, 

854 val.values, 

855 atol=1e-16, 

856 rtol=1e-16)) 

857 

858 def _makeWcs(self, dec=53.1595451514076, theta=0): 

859 """Create a wcs from real CFHT values, rotated by an optional theta. 

860 

861 dec : `float` 

862 Set reference declination of CRVAL2 [degrees] 

863 theta : `float` 

864 Rotate CD matrix by theta [degrees] 

865 

866 Returns 

867 ------- 

868 wcs : `lsst.afw.geom` 

869 Created wcs. 

870 """ 

871 metadata = dafBase.PropertySet() 

872 

873 metadata.set("SIMPLE", "T") 

874 metadata.set("BITPIX", -32) 

875 metadata.set("NAXIS", 2) 

876 metadata.set("NAXIS1", 1024) 

877 metadata.set("NAXIS2", 1153) 

878 metadata.set("RADECSYS", 'FK5') 

879 metadata.set("EQUINOX", 2000.) 

880 

881 metadata.setDouble("CRVAL1", 215.604025685476) 

882 metadata.setDouble("CRVAL2", dec) 

883 metadata.setDouble("CRPIX1", 1109.99981456774) 

884 metadata.setDouble("CRPIX2", 560.018167811613) 

885 metadata.set("CTYPE1", 'RA---SIN') 

886 metadata.set("CTYPE2", 'DEC--SIN') 

887 

888 cd_matrix = np.array( 

889 [ 

890 [5.10808596133527E-05, 1.85579539217196E-07], 

891 [-8.27440751733828E-07, -5.10281493481982E-05] 

892 ] 

893 ) 

894 # rotate CD matrix 

895 theta_rad = np.deg2rad(theta) 

896 rotation_matrix = np.array( 

897 [ 

898 [np.cos(theta_rad), -np.sin(theta_rad)], 

899 [np.sin(theta_rad), np.cos(theta_rad)], 

900 ] 

901 ) 

902 cd_matrix = rotation_matrix @ cd_matrix 

903 

904 metadata.setDouble("CD1_1", cd_matrix[0, 0]) 

905 metadata.setDouble("CD1_2", cd_matrix[0, 1]) 

906 metadata.setDouble("CD2_1", cd_matrix[1, 0]) 

907 metadata.setDouble("CD2_2", cd_matrix[1, 1]) 

908 

909 return afwGeom.makeSkyWcs(metadata) 

910 

911 def testHtm(self): 

912 """Test that HtmIndxes are created as expected. 

913 """ 

914 df = self.getMultiIndexDataFrame(self.dataDict) 

915 func = HtmIndex20("coord_ra", "coord_dec") 

916 

917 val = self._funcVal(func, df) 

918 # Test that the HtmIds come out as the ra/dec in dataDict. 

919 self.assertTrue(np.all(np.equal( 

920 val.values, 

921 [14924528684992, 14924528689697, 14924528501716, 14924526434259, 

922 14924526433879]))) 

923 

924 def testEbv(self): 

925 """Test that EBV works. 

926 """ 

927 df = self.getMultiIndexDataFrame(self.dataDict) 

928 func = Ebv() 

929 

930 val = self._funcVal(func, df) 

931 np.testing.assert_array_almost_equal( 

932 val.values, 

933 [0.029100, 0.029013, 0.028857, 0.028802, 0.028797] 

934 ) 

935 

936 def testSkyMoments(self): 

937 self.columns.extend([ 

938 "slot_Shape_xx", 

939 "slot_Shape_yy", 

940 "slot_Shape_xy", 

941 "base_LocalWcs_CDMatrix_1_1", 

942 "base_LocalWcs_CDMatrix_1_2", 

943 "base_LocalWcs_CDMatrix_2_1", 

944 "base_LocalWcs_CDMatrix_1_1"]) 

945 

946 # CD Matrix from a ComCam exposure. 

947 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = \ 

948 np.full(self.nRecords, -9.695088e-07) 

949 self.dataDict["base_LocalWcs_CDMatrix_1_2"] = \ 

950 np.full(self.nRecords, 3.950301849959342e-09) 

951 self.dataDict["base_LocalWcs_CDMatrix_2_1"] = \ 

952 np.full(self.nRecords, 3.8766915166433014e-09) 

953 self.dataDict["base_LocalWcs_CDMatrix_2_2"] = \ 

954 np.full(self.nRecords, 9.695092484727074e-07) 

955 self.dataDict["slot_Shape_xx"] = \ 

956 np.array([6.52675084, 74.17426471, 6.45283335, 36.82870958, 6.45685472]) 

957 self.dataDict["slot_Shape_yy"] = \ 

958 np.array([6.12848637, 80.99510036, 6.05671667, 35.79219613, 5.97778765]) 

959 self.dataDict["slot_Shape_xy"] = \ 

960 np.array([-0.10281872, 0.88788384, -0.1261287, -1.60504171, 0.11974515]) 

961 

962 col_names = ["slot_Shape_xx", "slot_Shape_yy", "slot_Shape_xy", 

963 "base_LocalWcs_CDMatrix_1_1", 

964 "base_LocalWcs_CDMatrix_1_2", 

965 "base_LocalWcs_CDMatrix_2_1", 

966 "base_LocalWcs_CDMatrix_2_2"] 

967 skyXX_functor = MomentsIuuSky(*col_names, filt="g") 

968 skyYY_functor = MomentsIvvSky(*col_names, filt="g") 

969 skyXY_functor = MomentsIuvSky(*col_names, filt="g") 

970 

971 axesA_functor = SemimajorAxisFromMoments(*col_names, filt="g") 

972 axesB_functor = SemiminorAxisFromMoments(*col_names, filt="g") 

973 axesTheta_functor = PositionAngleFromMoments(*col_names, filt="g") 

974 

975 df = self.getMultiIndexDataFrame(self.dataDict) 

976 output_sky_xx = skyXX_functor(df) 

977 output_sky_yy = skyYY_functor(df) 

978 output_sky_xy = skyXY_functor(df) 

979 

980 output_axes_a = axesA_functor(df) 

981 output_axes_b = axesB_functor(df) 

982 output_axes_theta = axesTheta_functor(df) 

983 

984 transformed_xx = [] 

985 transformed_yy = [] 

986 transformed_xy = [] 

987 axes_a = [] 

988 axes_b = [] 

989 axes_theta = [] 

990 

991 for n in range(5): 

992 Ixx = df[('meas', 'g', 'slot_Shape_xx')].iloc[n] 

993 Iyy = df[('meas', 'g', 'slot_Shape_yy')].iloc[n] 

994 Ixy = df[('meas', 'g', 'slot_Shape_xy')].iloc[n] 

995 localWCS_CD_1_1 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_1_1')].iloc[n] 

996 localWCS_CD_2_1 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_2_1')].iloc[n] 

997 localWCS_CD_1_2 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_1_2')].iloc[n] 

998 localWCS_CD_2_2 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_2_2')].iloc[n] 

999 CD_matrix = np.array([[localWCS_CD_1_1, localWCS_CD_1_2], 

1000 [localWCS_CD_2_1, localWCS_CD_2_2]]) 

1001 

1002 q = afwGeom.ellipses.Quadrupole(Ixx, Iyy, Ixy) 

1003 lt = geom.LinearTransform(CD_matrix) 

1004 transformed_q = q.transform(lt) 

1005 transformed_q.scale((180/np.pi) * (3600)) 

1006 

1007 axes = afwGeom.ellipses.Axes(transformed_q) 

1008 

1009 transformed_xx.append(transformed_q.getIxx()) 

1010 transformed_yy.append(transformed_q.getIyy()) 

1011 transformed_xy.append(transformed_q.getIxy()) 

1012 

1013 axes_a.append(axes.getA()) 

1014 axes_b.append(axes.getB()) 

1015 axes_theta.append(np.degrees(axes.getTheta())) 

1016 

1017 self.assertFloatsAlmostEqual(output_sky_xx, np.array(transformed_xx), rtol=1e-5) 

1018 self.assertFloatsAlmostEqual(output_sky_yy, np.array(transformed_yy), rtol=1e-5) 

1019 self.assertFloatsAlmostEqual(output_sky_xy, np.array(transformed_xy), rtol=1e-5) 

1020 

1021 self.assertFloatsAlmostEqual(output_axes_a, np.array(axes_a), rtol=1e-5) 

1022 self.assertFloatsAlmostEqual(output_axes_b, np.array(axes_b), rtol=1e-5) 

1023 self.assertFloatsAlmostEqual(output_axes_theta, np.array(axes_theta), rtol=1e-5) 

1024 

1025 def _dropLevels(self, df): 

1026 levelsToDrop = [n for lev, n in zip(df.columns.levels, df.columns.names) if len(lev) == 1] 

1027 

1028 # Prevent error when trying to drop *all* columns 

1029 if len(levelsToDrop) == len(df.columns.names): 

1030 levelsToDrop.remove(df.columns.names[-1]) 

1031 

1032 df.columns = df.columns.droplevel(levelsToDrop) 

1033 

1034 

1035class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase): 

1036 pass 

1037 

1038 

1039def setup_module(module): 

1040 lsst.utils.tests.init() 

1041 

1042 

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

1044 lsst.utils.tests.init() 

1045 unittest.main()