Coverage for tests/test_functors.py: 10%
546 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-05-29 08:48 +0000
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/>.
22import astropy.units as u
23import copy
24import functools
25import numpy as np
26import os
27import pandas as pd
28import unittest
29import logging
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 MomentsG1Sky, MomentsG2Sky, MomentsTraceSky,
54 CorrelationIuuSky, CorrelationIvvSky, CorrelationIuvSky,
55 SemimajorAxisFromCorrelation, SemiminorAxisFromCorrelation,
56 PositionAngleFromCorrelation,
57 ConvertDetectorAngleErrToPositionAngleErr)
59ROOT = os.path.abspath(os.path.dirname(__file__))
62class FunctorTestCase(lsst.utils.tests.TestCase):
64 def getMultiIndexDataFrame(self, dataDict):
65 """Create a simple test multi-index DataFrame."""
67 simpleDF = pd.DataFrame(dataDict)
68 dfFilterDSCombos = []
69 for ds in self.datasets:
70 for band in self.bands:
71 df = copy.copy(simpleDF)
72 df.reindex(sorted(df.columns), axis=1)
73 df['dataset'] = ds
74 df['band'] = band
75 df.columns = pd.MultiIndex.from_tuples(
76 [(ds, band, c) for c in df.columns],
77 names=('dataset', 'band', 'column'))
78 dfFilterDSCombos.append(df)
80 df = functools.reduce(lambda d1, d2: d1.join(d2), dfFilterDSCombos)
82 return df
84 def getSimpleDataFrame(self, dataDict):
85 return pd.DataFrame(dataDict)
87 def getDatasetHandle(self, df):
88 lo, hi = HtmPixelization(7).universe().ranges()[0]
89 value = np.random.randint(lo, hi)
90 return InMemoryDatasetHandle(df, storageClass="DataFrame", dataId={"htm7": value})
92 def setUp(self):
93 np.random.seed(12345)
94 self.datasets = ['forced_src', 'meas', 'ref']
95 self.bands = ['g', 'r']
96 self.columns = ['coord_ra', 'coord_dec']
97 self.nRecords = 5
98 self.dataDict = {
99 "coord_ra": [3.77654137, 3.77643059, 3.77621148, 3.77611944, 3.77610396],
100 "coord_dec": [0.01127624, 0.01127787, 0.01127543, 0.01127543, 0.01127543]}
102 def _funcVal(self, functor, df):
103 self.assertIsInstance(functor.name, str)
104 self.assertIsInstance(functor.shortname, str)
106 handle = self.getDatasetHandle(df)
108 val = functor(df)
109 val2 = functor(handle)
110 self.assertTrue((val == val2).all())
111 self.assertIsInstance(val, pd.Series)
113 val = functor(df, dropna=True)
114 val2 = functor(handle, dropna=True)
115 self.assertTrue((val == val2).all())
116 self.assertEqual(val.isnull().sum(), 0)
118 return val
120 def _differenceVal(self, functor, df1, df2):
121 self.assertIsInstance(functor.name, str)
122 self.assertIsInstance(functor.shortname, str)
124 handle1 = self.getDatasetHandle(df1)
125 handle2 = self.getDatasetHandle(df2)
127 val = functor.difference(df1, df2)
128 val2 = functor.difference(handle1, handle2)
129 self.assertTrue(val.equals(val2))
130 self.assertIsInstance(val, pd.Series)
132 val = functor.difference(df1, df2, dropna=True)
133 val2 = functor.difference(handle1, handle2, dropna=True)
134 self.assertTrue(val.equals(val2))
135 self.assertEqual(val.isnull().sum(), 0)
137 val1 = self._funcVal(functor, df1)
138 val2 = self._funcVal(functor, df2)
140 self.assertTrue(np.allclose(val, val1 - val2))
142 return val
144 def testColumn(self):
145 self.columns.append("base_FootprintArea_value")
146 self.dataDict["base_FootprintArea_value"] = \
147 np.full(self.nRecords, 1)
148 df = self.getMultiIndexDataFrame(self.dataDict)
149 func = Column('base_FootprintArea_value', filt='g')
150 self._funcVal(func, df)
152 df = self.getSimpleDataFrame(self.dataDict)
153 func = Column('base_FootprintArea_value')
154 self._funcVal(func, df)
156 def testCustom(self):
157 self.columns.append("base_FootprintArea_value")
158 self.dataDict["base_FootprintArea_value"] = \
159 np.random.rand(self.nRecords)
160 df = self.getMultiIndexDataFrame(self.dataDict)
161 func = CustomFunctor('2*base_FootprintArea_value', filt='g')
162 val = self._funcVal(func, df)
164 func2 = Column('base_FootprintArea_value', filt='g')
166 np.allclose(val.values, 2*func2(df).values, atol=1e-13, rtol=0)
168 df = self.getSimpleDataFrame(self.dataDict)
169 func = CustomFunctor('2 * base_FootprintArea_value')
170 val = self._funcVal(func, df)
171 func2 = Column('base_FootprintArea_value')
173 np.allclose(val.values, 2*func2(df).values, atol=1e-13, rtol=0)
175 def testCoords(self):
176 df = self.getMultiIndexDataFrame(self.dataDict)
177 ra = self._funcVal(RAColumn(), df)
178 dec = self._funcVal(DecColumn(), df)
180 columnDict = {'dataset': 'ref', 'band': 'g',
181 'column': ['coord_ra', 'coord_dec']}
183 handle = InMemoryDatasetHandle(df, storageClass="DataFrame")
184 dfSub = handle.get(parameters={"columns": columnDict})
185 self._dropLevels(dfSub)
187 coords = dfSub / np.pi * 180.
189 self.assertTrue(np.allclose(ra, coords[('ref', 'g', 'coord_ra')], atol=1e-13, rtol=0))
190 self.assertTrue(np.allclose(dec, coords[('ref', 'g', 'coord_dec')], atol=1e-13, rtol=0))
192 # single-level column index table
193 df = self.getSimpleDataFrame(self.dataDict)
194 ra = self._funcVal(RAColumn(), df)
195 dec = self._funcVal(DecColumn(), df)
197 handle = InMemoryDatasetHandle(df, storageClass="DataFrame")
198 dfSub = handle.get(parameters={"columns": ['coord_ra', 'coord_dec']})
199 coords = dfSub / np.pi * 180.
201 self.assertTrue(np.allclose(ra, coords['coord_ra'], atol=1e-13, rtol=0))
202 self.assertTrue(np.allclose(dec, coords['coord_dec'], atol=1e-13, rtol=0))
204 def testMag(self):
205 self.columns.extend(["base_PsfFlux_instFlux", "base_PsfFlux_instFluxErr"])
206 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
207 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10)
208 df = self.getMultiIndexDataFrame(self.dataDict)
209 # Change one dataset filter combinations value.
210 df[("meas", "g", "base_PsfFlux_instFlux")] -= 1
212 fluxName = 'base_PsfFlux'
214 # Check that things work when you provide dataset explicitly
215 for dataset in ['forced_src', 'meas']:
216 psfMag_G = self._funcVal(Mag(fluxName, dataset=dataset,
217 filt='g'),
218 df)
219 psfMag_R = self._funcVal(Mag(fluxName, dataset=dataset,
220 filt='r'),
221 df)
223 psfColor_GR = self._funcVal(Color(fluxName, 'g', 'r',
224 dataset=dataset),
225 df)
227 self.assertTrue(np.allclose((psfMag_G - psfMag_R).dropna(), psfColor_GR, rtol=0, atol=1e-13))
229 # Check that behavior as expected when dataset not provided;
230 # that is, that the color comes from forced and default Mag is meas
231 psfMag_G = self._funcVal(Mag(fluxName, filt='g'), df)
232 psfMag_R = self._funcVal(Mag(fluxName, filt='r'), df)
234 psfColor_GR = self._funcVal(Color(fluxName, 'g', 'r'), df)
236 # These should *not* be equal.
237 self.assertFalse(np.allclose((psfMag_G - psfMag_R).dropna(), psfColor_GR))
239 def testMagDiff(self):
240 self.columns.extend(["base_PsfFlux_instFlux", "base_PsfFlux_instFluxErr",
241 "modelfit_CModel_instFlux", "modelfit_CModel_instFluxErr"])
242 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
243 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10)
244 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1000)
245 self.dataDict["modelfit_CModel_instFluxErr"] = np.full(self.nRecords, 10)
246 df = self.getMultiIndexDataFrame(self.dataDict)
248 for filt in self.bands:
249 filt = 'g'
250 val = self._funcVal(MagDiff('base_PsfFlux', 'modelfit_CModel', filt=filt), df)
252 mag1 = self._funcVal(Mag('modelfit_CModel', filt=filt), df)
253 mag2 = self._funcVal(Mag('base_PsfFlux', filt=filt), df)
254 self.assertTrue(np.allclose((mag2 - mag1).dropna(), val, rtol=0, atol=1e-13))
256 def testDifference(self):
257 """Test .difference method using MagDiff as the example.
258 """
259 self.columns.extend(["base_PsfFlux_instFlux", "base_PsfFlux_instFluxErr",
260 "modelfit_CModel_instFlux", "modelfit_CModel_instFluxErr"])
262 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
263 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1000)
264 df1 = self.getMultiIndexDataFrame(self.dataDict)
266 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 999)
267 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 999)
268 df2 = self.getMultiIndexDataFrame(self.dataDict)
270 magDiff = MagDiff('base_PsfFlux', 'modelfit_CModel', filt='g')
272 # Asserts that differences computed properly
273 self._differenceVal(magDiff, df1, df2)
275 def testPixelScale(self):
276 # Test that the pixel scale and pix->arcsec calculations perform as
277 # expected.
278 pass
280 def testShape(self):
281 data = {
282 "x": np.array([-0.3, 0.4, 0.7, -0.9, 1.4, -5.3]),
283 "y": np.array([1.5, -0.7, -1.9, 2.8, -1.4, 0.01]),
284 "rho": np.array([-0.9, 0.4, -0.7, 0., 0.3, -0.99]),
285 }
286 data["xx"] = data["x"]**2
287 data["yy"] = data["y"]**2
288 data["xy"] = data["x"]*data["y"]*data["rho"]
290 args = ("xx", "xy", "yy")
291 functor_e1, functor_e2, functor_quadrupole = E1(*args), E2(*args), RadiusFromQuadrupole(*args)
293 xx_plus_yy = data["xx"] + data["yy"]
294 data = pd.DataFrame(data)
296 np.testing.assert_allclose(
297 functor_e1(data),
298 ((data["xx"] - data["yy"])/xx_plus_yy).astype(np.float32),
299 rtol=1e-12, atol=1e-12,
300 )
301 np.testing.assert_allclose(
302 functor_e2(data),
303 (2.0*data["xy"]/xx_plus_yy).astype(np.float32),
304 rtol=1e-12, atol=1e-12,
305 )
306 np.testing.assert_allclose(
307 functor_quadrupole(data),
308 ((data["xx"]*data["yy"] - data["xy"]**2)**0.25).astype(np.float32),
309 rtol=1e-12, atol=1e-12,
310 )
312 def testOther(self):
313 self.columns.extend(["ext_shapeHSM_HsmSourceMoments_xx", "ext_shapeHSM_HsmSourceMoments_yy",
314 "base_SdssShape_xx", "base_SdssShape_yy",
315 "ext_shapeHSM_HsmPsfMoments_xx", "ext_shapeHSM_HsmPsfMoments_yy",
316 "base_SdssShape_psf_xx", "base_SdssShape_psf_yy"])
317 self.dataDict["ext_shapeHSM_HsmSourceMoments_xx"] = np.full(self.nRecords, 1 / np.sqrt(2))
318 self.dataDict["ext_shapeHSM_HsmSourceMoments_yy"] = np.full(self.nRecords, 1 / np.sqrt(2))
319 self.dataDict["base_SdssShape_xx"] = np.full(self.nRecords, 1 / np.sqrt(2))
320 self.dataDict["base_SdssShape_yy"] = np.full(self.nRecords, 1 / np.sqrt(2))
321 self.dataDict["ext_shapeHSM_HsmPsfMoments_xx"] = np.full(self.nRecords, 1 / np.sqrt(2))
322 self.dataDict["ext_shapeHSM_HsmPsfMoments_yy"] = np.full(self.nRecords, 1 / np.sqrt(2))
323 self.dataDict["base_SdssShape_psf_xx"] = np.full(self.nRecords, 1)
324 self.dataDict["base_SdssShape_psf_yy"] = np.full(self.nRecords, 1)
325 df = self.getMultiIndexDataFrame(self.dataDict)
326 # Covering the code is better than nothing
327 for filt in self.bands:
328 for Func in [DeconvolvedMoments,
329 SdssTraceSize,
330 PsfSdssTraceSizeDiff,
331 HsmTraceSize, PsfHsmTraceSizeDiff, HsmFwhm]:
332 _ = self._funcVal(Func(filt=filt), df)
334 def _compositeFuncVal(self, functor, df):
335 self.assertIsInstance(functor, CompositeFunctor)
337 handle = self.getDatasetHandle(df)
339 fdf1 = functor(df)
340 fdf2 = functor(handle)
341 self.assertTrue(fdf1.equals(fdf2))
343 self.assertIsInstance(fdf1, pd.DataFrame)
344 self.assertTrue(np.all([k in fdf1.columns for k in functor.funcDict.keys()]))
346 fdf1 = functor(df, dropna=True)
347 fdf2 = functor(handle, dropna=True)
348 self.assertTrue(fdf1.equals(fdf2))
350 # Check that there are no nulls
351 self.assertFalse(fdf1.isnull().any(axis=None))
353 return fdf1
355 def _compositeDifferenceVal(self, functor, df1, df2):
356 self.assertIsInstance(functor, CompositeFunctor)
358 handle1 = self.getDatasetHandle(df1)
359 handle2 = self.getDatasetHandle(df2)
361 fdf1 = functor.difference(df1, df2)
362 fdf2 = functor.difference(handle1, handle2)
363 self.assertTrue(fdf1.equals(fdf2))
365 self.assertIsInstance(fdf1, pd.DataFrame)
366 self.assertTrue(np.all([k in fdf1.columns for k in functor.funcDict.keys()]))
368 fdf1 = functor.difference(df1, df2, dropna=True)
369 fdf2 = functor.difference(handle1, handle2, dropna=True)
370 self.assertTrue(fdf1.equals(fdf2))
372 # Check that there are no nulls
373 self.assertFalse(fdf1.isnull().any(axis=None))
375 df1_functored = functor(df1)
376 df2_functored = functor(df2)
378 self.assertTrue(np.allclose(fdf1.values, df1_functored.values - df2_functored.values))
380 return fdf1
382 def testComposite(self):
383 self.columns.extend(["modelfit_CModel_instFlux", "base_PsfFlux_instFlux"])
384 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1)
385 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1)
387 df = self.getMultiIndexDataFrame(self.dataDict)
388 # Modify r band value slightly.
389 df[("meas", "r", "base_PsfFlux_instFlux")] -= 0.1
391 filt = 'g'
392 funcDict = {'psfMag_ref': Mag('base_PsfFlux', dataset='ref'),
393 'ra': RAColumn(),
394 'dec': DecColumn(),
395 'psfMag': Mag('base_PsfFlux', filt=filt),
396 'cmodel_magDiff': MagDiff('base_PsfFlux',
397 'modelfit_CModel', filt=filt)}
398 func = CompositeFunctor(funcDict)
399 fdf1 = self._compositeFuncVal(func, df)
401 # Repeat same, but define filter globally instead of individually
402 funcDict2 = {'psfMag_ref': Mag('base_PsfFlux', dataset='ref'),
403 'ra': RAColumn(),
404 'dec': DecColumn(),
405 'psfMag': Mag('base_PsfFlux'),
406 'cmodel_magDiff': MagDiff('base_PsfFlux',
407 'modelfit_CModel')}
409 func2 = CompositeFunctor(funcDict2, filt=filt)
410 fdf2 = self._compositeFuncVal(func2, df)
411 self.assertTrue(fdf1.equals(fdf2))
413 func2.filt = 'r'
414 fdf3 = self._compositeFuncVal(func2, df)
415 # Because we modified the R filter this should fail.
416 self.assertFalse(fdf2.equals(fdf3))
418 # Make sure things work with passing list instead of dict
419 funcs = [Mag('base_PsfFlux', dataset='ref'),
420 RAColumn(),
421 DecColumn(),
422 Mag('base_PsfFlux', filt=filt),
423 MagDiff('base_PsfFlux', 'modelfit_CModel', filt=filt)]
425 _ = self._compositeFuncVal(CompositeFunctor(funcs), df)
427 def testCompositeSimple(self):
428 """Test single-level composite functor for functionality
429 """
430 self.columns.extend(["modelfit_CModel_instFlux", "base_PsfFlux_instFlux"])
431 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1)
432 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1)
434 df = self.getSimpleDataFrame(self.dataDict)
435 funcDict = {'ra': RAColumn(),
436 'dec': DecColumn(),
437 'psfMag': Mag('base_PsfFlux'),
438 'cmodel_magDiff': MagDiff('base_PsfFlux',
439 'modelfit_CModel')}
440 func = CompositeFunctor(funcDict)
441 _ = self._compositeFuncVal(func, df)
443 def testCompositeColor(self):
444 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
445 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10)
446 df = self.getMultiIndexDataFrame(self.dataDict)
447 funcDict = {'a': Mag('base_PsfFlux', dataset='meas', filt='g'),
448 'b': Mag('base_PsfFlux', dataset='forced_src', filt='g'),
449 'c': Color('base_PsfFlux', 'g', 'r')}
450 # Covering the code is better than nothing
451 _ = self._compositeFuncVal(CompositeFunctor(funcDict), df)
453 def testCompositeDifference(self):
454 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
455 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10)
456 df1 = self.getMultiIndexDataFrame(self.dataDict)
458 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 999)
459 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 9)
460 df2 = self.getMultiIndexDataFrame(self.dataDict)
462 funcDict = {'a': Mag('base_PsfFlux', dataset='meas', filt='g'),
463 'b': Mag('base_PsfFlux', dataset='forced_src', filt='g'),
464 'c': Color('base_PsfFlux', 'g', 'r')}
465 # Covering the code is better than nothing
466 _ = self._compositeDifferenceVal(CompositeFunctor(funcDict), df1, df2)
468 def testCompositeFail(self):
469 """Test a composite functor where one of the functors should be junk.
470 """
471 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
472 df = self.getMultiIndexDataFrame(self.dataDict)
474 funcDict = {'good': Column("base_PsfFlux_instFlux"),
475 'bad': Column('not_a_column')}
477 with self.assertLogs(level=logging.ERROR) as cm:
478 _ = self._compositeFuncVal(CompositeFunctor(funcDict), df)
479 self.assertIn("Exception in CompositeFunctor (funcs: ['good', 'bad'])", cm.output[0])
481 def testLocalPhotometry(self):
482 """Test the local photometry functors.
483 """
484 flux = 1000
485 fluxErr = 10
486 calib = 10
487 calibErr = 0.0
488 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, flux)
489 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords,
490 fluxErr)
491 self.dataDict["base_LocalPhotoCalib"] = np.full(self.nRecords, calib)
493 df = self.getMultiIndexDataFrame(self.dataDict)
494 func = LocalPhotometry("base_PsfFlux_instFlux",
495 "base_PsfFlux_instFluxErr",
496 "base_LocalPhotoCalib")
498 nanoJansky = func.instFluxToNanojansky(
499 df[("meas", "g", "base_PsfFlux_instFlux")],
500 df[("meas", "g", "base_LocalPhotoCalib")])
501 mag = func.instFluxToMagnitude(
502 df[("meas", "g", "base_PsfFlux_instFlux")],
503 df[("meas", "g", "base_LocalPhotoCalib")])
504 nanoJanskyErr = func.instFluxErrToNanojanskyErr(
505 df[("meas", "g", "base_PsfFlux_instFlux")],
506 df[("meas", "g", "base_PsfFlux_instFluxErr")],
507 df[("meas", "g", "base_LocalPhotoCalib")])
508 magErr = func.instFluxErrToMagnitudeErr(
509 df[("meas", "g", "base_PsfFlux_instFlux")],
510 df[("meas", "g", "base_PsfFlux_instFluxErr")],
511 df[("meas", "g", "base_LocalPhotoCalib")])
513 self.assertTrue(np.allclose(nanoJansky.values,
514 flux * calib,
515 atol=0,
516 rtol=1e-7))
517 self.assertTrue(np.allclose(mag.values,
518 (flux * calib * u.nJy).to_value(u.ABmag),
519 atol=0,
520 rtol=1e-7))
521 self.assertTrue(np.allclose(nanoJanskyErr.values,
522 np.hypot(fluxErr * calib, flux * calibErr),
523 atol=0,
524 rtol=1e-7))
525 self.assertTrue(np.allclose(
526 magErr.values,
527 2.5 / np.log(10) * nanoJanskyErr.values / nanoJansky.values,
528 atol=0,
529 rtol=1e-7))
531 # Test functors against the values computed above.
532 self._testLocalPhotometryFunctors(LocalNanojansky,
533 df,
534 nanoJansky)
535 self._testLocalPhotometryFunctors(LocalNanojanskyErr,
536 df,
537 nanoJanskyErr)
539 def _testLocalPhotometryFunctors(self, functor, df, testValues):
540 func = functor("base_PsfFlux_instFlux",
541 "base_PsfFlux_instFluxErr",
542 "base_LocalPhotoCalib")
543 val = self._funcVal(func, df)
544 self.assertTrue(np.allclose(testValues.values,
545 val.values,
546 atol=0,
547 rtol=1e-7))
549 def testDipPhotometry(self):
550 """Test calibrated flux calculations for dipoles."""
551 fluxNeg = -100
552 fluxPos = 101
553 fluxErr = 10
554 calib = 10
555 calibErr = 0.0
557 # compute expected values.
558 absMean = 0.5*(fluxPos - fluxNeg)*calib
559 absDiff = (fluxNeg + fluxPos)*calib
560 absMeanErr = 0.5*np.sqrt(2*(fluxErr*calib)**2
561 + ((fluxPos - fluxNeg)*calibErr)**2)
562 absDiffErr = np.sqrt(2*(fluxErr*calib)**2
563 + ((fluxPos + fluxNeg)*calibErr)**2)
565 self.dataDict["ip_diffim_DipoleFluxNeg_instFlux"] = np.full(self.nRecords, fluxNeg)
566 self.dataDict["ip_diffim_DipoleFluxNeg_instFluxErr"] = np.full(self.nRecords, fluxErr)
567 self.dataDict["ip_diffim_DipoleFluxPos_instFlux"] = np.full(self.nRecords, fluxPos)
568 self.dataDict["ip_diffim_DipoleFluxPos_instFluxErr"] = np.full(self.nRecords, fluxErr)
569 self.dataDict["base_LocalPhotoCalib"] = np.full(self.nRecords, calib)
571 df = self.getMultiIndexDataFrame(self.dataDict)
572 func = LocalDipoleMeanFlux("ip_diffim_DipoleFluxPos_instFlux",
573 "ip_diffim_DipoleFluxNeg_instFlux",
574 "ip_diffim_DipoleFluxPos_instFluxErr",
575 "ip_diffim_DipoleFluxNeg_instFluxErr",
576 "base_LocalPhotoCalib")
577 val = self._funcVal(func, df)
578 self.assertFloatsAlmostEqual(val.values,
579 absMean,
580 atol=1e-13,
581 rtol=0)
583 func = LocalDipoleMeanFluxErr("ip_diffim_DipoleFluxPos_instFlux",
584 "ip_diffim_DipoleFluxNeg_instFlux",
585 "ip_diffim_DipoleFluxPos_instFluxErr",
586 "ip_diffim_DipoleFluxNeg_instFluxErr",
587 "base_LocalPhotoCalib")
588 val = self._funcVal(func, df)
589 self.assertFloatsAlmostEqual(val.values,
590 absMeanErr,
591 atol=1e-13,
592 rtol=0)
594 func = LocalDipoleDiffFlux("ip_diffim_DipoleFluxPos_instFlux",
595 "ip_diffim_DipoleFluxNeg_instFlux",
596 "ip_diffim_DipoleFluxPos_instFluxErr",
597 "ip_diffim_DipoleFluxNeg_instFluxErr",
598 "base_LocalPhotoCalib")
599 val = self._funcVal(func, df)
600 self.assertFloatsAlmostEqual(val.values,
601 absDiff,
602 atol=1e-13,
603 rtol=0)
605 func = LocalDipoleDiffFluxErr("ip_diffim_DipoleFluxPos_instFlux",
606 "ip_diffim_DipoleFluxNeg_instFlux",
607 "ip_diffim_DipoleFluxPos_instFluxErr",
608 "ip_diffim_DipoleFluxNeg_instFluxErr",
609 "base_LocalPhotoCalib")
610 val = self._funcVal(func, df)
611 self.assertFloatsAlmostEqual(val.values,
612 absDiffErr,
613 atol=1e-13,
614 rtol=0)
616 def testComputePositionAngle(self, offset=0.0001):
617 """Test computation of position angle from (RA1, Dec1) to (RA2, Dec2)
619 offset : `float`
620 Arc length of the offset vector to set up test points. [radian]
621 """
623 d = offset
624 columns = ("ra1", "dec1", "ra2", "dec2", "expected")
625 position_angle_test_values = [
626 # Get 0, 0 right
627 (0, 0, d, 0, np.pi/2),
628 (0, 0, 0, d, 0),
629 (0, 0, -d, 0, -np.pi/2),
630 (0, 0, 0, -d, np.pi),
631 # Make sure we get wrap-around to 0, 0 right
632 (2*np.pi, 0, d, 0, np.pi/2),
633 (2*np.pi, 0, 0, d, 0),
634 (2*np.pi, 0, -d, 0, -np.pi/2),
635 (2*np.pi, 0, 0, -d, np.pi),
636 (+0.0015, 0, 2*np.pi - 0.05, 0, -np.pi/2),
637 # Get another somewhat arbitrary location right [these are in rad]
638 # It's not really important to rescale d by 1/cos(dec)
639 # because we're just looking at orientation of vector
640 # but foreshadowing to the poles...
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 # Negative dec
646 (0.0015, -1, 0.0015 + d*np.cos(-1), -1, np.pi/2),
647 (0.0015, -1, 0.0015, -1 + d, 0),
648 (0.0015, -1, 0.0015 - d*np.cos(-1), -1, - np.pi/2),
649 (0.0015, -1, 0.0015, -1 - d, np.pi),
650 # Make sure we get wrap-around on that right
651 (2*np.pi + 0.0015, 1, 2*np.pi + 0.0015 + d*np.cos(1), 1, np.pi/2),
652 (2*np.pi + 0.0015, 1, 2*np.pi + 0.0015, 1 + d, 0),
653 (2*np.pi + 0.0015, 1, 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 relative wrap-around right
656 (2*np.pi + 0.0015, 1, 0.0015 + d*np.cos(1), 1, np.pi/2),
657 (0.0015, 1, 2*np.pi + 0.0015, 1 + d, 0),
658 (0.0015, 1, 2*np.pi + 0.0015 - d*np.cos(1), 1, - np.pi/2),
659 (2*np.pi + 0.0015, 1, 0.0015, 1 - d, np.pi),
660 # Get either side of RA=0 right
661 (+ d*np.cos(1) / 2, 1, 2*np.pi - d*np.cos(1) / 2, 1, - np.pi/2),
662 (+ d*np.cos(1) / 2, 1, - d*np.cos(1) / 2, 1, -np.pi/2),
663 (2*np.pi + d*np.cos(1) / 2, 1, - d*np.cos(1) / 2, 1, -np.pi/2),
664 (-d*np.cos(1) / 2, 1, +0.0015, 1, np.pi/2),
665 (0.0015, 1, 0.0015, 1 + d, 0),
666 (0.0015, 1, 0.0015 - d*np.cos(1), 1, - np.pi/2),
667 (0.0015, 1, 0.0015, 1 - d, np.pi),
668 # Try it near the poles
669 (0, np.pi/2, 0, np.pi/2 - d, np.pi),
670 (0, np.pi/2 - d, 0, np.pi/2, 0),
671 (0, -np.pi/2, 0, -np.pi/2 + d, 0),
672 (0, -np.pi/2 + d, 0, -np.pi/2, np.pi),
673 ]
675 df = pd.DataFrame(position_angle_test_values, columns=columns)
677 cd_matrix = np.array([[1, 0], [0, -1]]) # Doesn't matter because we don't use it.
678 local_wcs = LocalWcs(*cd_matrix.flatten())
679 pa = local_wcs.computePositionAngle(df["ra1"], df["dec1"], df["ra2"], df["dec2"])
680 expected = df["expected"]
681 tolerance_deg = 0.05 # degrees
682 tolerance_rad = np.deg2rad(tolerance_deg)
684 # Use SphereGeom to handle wrap-around separations correctly.
685 diff = [
686 geom.Angle(o, geom.radians).separation(geom.Angle(e, geom.radians)).asRadians()
687 for o, e in zip(pa, expected)
688 ]
690 np.testing.assert_allclose(diff, 0, rtol=0, atol=tolerance_rad)
692 # Test position angle
693 def testConvertDetectorAngleToPositionAngle(self):
694 """Test conversion of position angle in detector degrees to position angle on sky
696 There is overlap with testConvertPixelToArcseconds
697 But we also test additional rotation angles so this is separate.
698 """
699 dipoleSep = 10
700 ixx = 10
701 testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2))
702 # testConvertPixelToArcSecond uses a meas_base LocalWcsPlugin
703 # but we're using a simple WCS as our example, so this doesn't really matter
704 # and we'll just use the WCS directly
705 for dec in np.linspace(-90, 90, 10):
706 for theta in (-45, 0, 90):
707 for x, y in zip(np.random.uniform(2 * 1109.99981456774, size=10),
708 np.random.uniform(2 * 560.018167811613, size=10)):
709 wcs = self._makeWcs(dec=dec, theta=theta)
710 cd_matrix = wcs.getCdMatrix()
712 self.dataDict["dipoleSep"] = np.full(self.nRecords, dipoleSep)
713 self.dataDict["ixx"] = np.full(self.nRecords, ixx)
714 self.dataDict["slot_Centroid_x"] = np.full(self.nRecords, x)
715 self.dataDict["slot_Centroid_y"] = np.full(self.nRecords, y)
716 self.dataDict["someCentroid_x"] = x + testPixelDeltas[:, 0]
717 self.dataDict["someCentroid_y"] = y + testPixelDeltas[:, 1]
718 self.dataDict["orientation"] = np.arctan2(
719 self.dataDict["someCentroid_y"] - self.dataDict["slot_Centroid_y"],
720 self.dataDict["someCentroid_x"] - self.dataDict["slot_Centroid_x"],
721 )
723 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(self.nRecords,
724 cd_matrix[0, 0])
725 self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(self.nRecords,
726 cd_matrix[0, 1])
727 self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(self.nRecords,
728 cd_matrix[1, 0])
729 self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(self.nRecords,
730 cd_matrix[1, 1])
731 df = self.getMultiIndexDataFrame(self.dataDict)
733 # Test detector angle to position angle conversion
734 func = ConvertDetectorAngleToPositionAngle(
735 "orientation",
736 "base_LocalWcs_CDMatrix_1_1",
737 "base_LocalWcs_CDMatrix_1_2",
738 "base_LocalWcs_CDMatrix_2_1",
739 "base_LocalWcs_CDMatrix_2_2"
740 )
741 val = self._funcVal(func, df)
743 delta_ra, delta_dec = func.computeDeltaRaDec(
744 self.dataDict["someCentroid_x"] - self.dataDict["slot_Centroid_x"],
745 self.dataDict["someCentroid_y"] - self.dataDict["slot_Centroid_y"],
746 self.dataDict["base_LocalWcs_CDMatrix_1_1"],
747 self.dataDict["base_LocalWcs_CDMatrix_1_2"],
748 self.dataDict["base_LocalWcs_CDMatrix_2_1"],
749 self.dataDict["base_LocalWcs_CDMatrix_2_2"],
750 )
752 dx = np.cos(0) * np.tan(delta_dec) - np.sin(0) * np.cos(delta_ra)
753 dy = np.sin(delta_ra)
754 comparison_pa = np.arctan2(dy, dx)
755 comparison_pa = np.rad2deg(comparison_pa)
757 coord_diff = []
758 for v, c in zip(val.values, comparison_pa):
759 observed_angle = geom.Angle(v, geom.degrees)
760 expected_angle = geom.Angle(c, geom.degrees)
761 diff = observed_angle.separation(expected_angle).asRadians()
762 coord_diff.append(diff)
764 np.testing.assert_allclose(coord_diff, 0, rtol=0, atol=5e-6)
766 # Test position angle error
767 def testConvertDetectorAngleErrToPositionAngleErr(self):
768 """Test conversion of detector angle error in radians to position angle error in degrees."""
769 theta_err_values = np.random.uniform(0, np.pi / 4, size=self.nRecords)
770 self.dataDict["theta_err"] = theta_err_values
771 df = self.getSimpleDataFrame(self.dataDict)
773 func = ConvertDetectorAngleErrToPositionAngleErr("theta_err")
774 val = self._funcVal(func, df)
776 expected = np.rad2deg(theta_err_values)
777 np.testing.assert_allclose(val.values, expected, rtol=0, atol=1e-13)
779 def testConvertPixelToArcseconds(self):
780 """Test calculations of the pixel scale, conversions of pixel to
781 arcseconds.
782 """
783 dipoleSep = 10
784 ixx = 10
785 testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2))
786 localWcsPlugin = measBase.EvaluateLocalWcsPlugin(
787 None,
788 "base_LocalWcs",
789 afwTable.SourceTable.makeMinimalSchema(),
790 None)
791 for dec in np.linspace(-90, 90, 10):
792 for x, y in zip(np.random.uniform(2 * 1109.99981456774, size=10),
793 np.random.uniform(2 * 560.018167811613, size=10)):
794 center = geom.Point2D(x, y)
795 wcs = self._makeWcs(dec)
796 skyOrigin = wcs.pixelToSky(center)
798 linAffMatrix = localWcsPlugin.makeLocalTransformMatrix(wcs,
799 center)
800 self.dataDict["dipoleSep"] = np.full(self.nRecords, dipoleSep)
801 self.dataDict["ixx"] = np.full(self.nRecords, ixx)
802 self.dataDict["slot_Centroid_x"] = np.full(self.nRecords, x)
803 self.dataDict["slot_Centroid_y"] = np.full(self.nRecords, y)
804 self.dataDict["someCentroid_x"] = x + testPixelDeltas[:, 0]
805 self.dataDict["someCentroid_y"] = y + testPixelDeltas[:, 1]
807 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(self.nRecords,
808 linAffMatrix[0, 0])
809 self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(self.nRecords,
810 linAffMatrix[0, 1])
811 self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(self.nRecords,
812 linAffMatrix[1, 0])
813 self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(self.nRecords,
814 linAffMatrix[1, 1])
815 df = self.getMultiIndexDataFrame(self.dataDict)
816 func = LocalWcs("base_LocalWcs_CDMatrix_1_1",
817 "base_LocalWcs_CDMatrix_1_2",
818 "base_LocalWcs_CDMatrix_2_1",
819 "base_LocalWcs_CDMatrix_2_2")
821 # Exercise the full set of functions in LocalWcs.
822 sepRadians = func.getSkySeparationFromPixel(
823 df[("meas", "g", "someCentroid_x")] - df[("meas", "g", "slot_Centroid_x")],
824 df[("meas", "g", "someCentroid_y")] - df[("meas", "g", "slot_Centroid_y")],
825 0.0,
826 0.0,
827 df[("meas", "g", "base_LocalWcs_CDMatrix_1_1")],
828 df[("meas", "g", "base_LocalWcs_CDMatrix_1_2")],
829 df[("meas", "g", "base_LocalWcs_CDMatrix_2_1")],
830 df[("meas", "g", "base_LocalWcs_CDMatrix_2_2")])
832 # Test functor values against afw SkyWcs computations.
833 for centX, centY, sep in zip(testPixelDeltas[:, 0],
834 testPixelDeltas[:, 1],
835 sepRadians.values):
836 afwSepRadians = skyOrigin.separation(
837 wcs.pixelToSky(x + centX, y + centY)).asRadians()
838 self.assertAlmostEqual(1 - sep / afwSepRadians, 0, places=5)
840 # Test the pixel scale computation.
841 func = ComputePixelScale("base_LocalWcs_CDMatrix_1_1",
842 "base_LocalWcs_CDMatrix_1_2",
843 "base_LocalWcs_CDMatrix_2_1",
844 "base_LocalWcs_CDMatrix_2_2")
845 pixelScale = self._funcVal(func, df)
846 self.assertTrue(np.allclose(
847 wcs.getPixelScale(center).asArcseconds(),
848 pixelScale.values,
849 rtol=1e-6,
850 atol=0))
852 # Test pixel -> arcsec conversion.
853 func = ConvertPixelToArcseconds("dipoleSep",
854 "base_LocalWcs_CDMatrix_1_1",
855 "base_LocalWcs_CDMatrix_1_2",
856 "base_LocalWcs_CDMatrix_2_1",
857 "base_LocalWcs_CDMatrix_2_2")
858 val = self._funcVal(func, df)
859 self.assertTrue(np.allclose(pixelScale.values * dipoleSep,
860 val.values,
861 atol=1e-16,
862 rtol=1e-16))
864 # Test pixel^2 -> arcsec^2 conversion.
865 func = ConvertPixelSqToArcsecondsSq("ixx",
866 "base_LocalWcs_CDMatrix_1_1",
867 "base_LocalWcs_CDMatrix_1_2",
868 "base_LocalWcs_CDMatrix_2_1",
869 "base_LocalWcs_CDMatrix_2_2")
870 val = self._funcVal(func, df)
871 self.assertTrue(np.allclose(pixelScale.values ** 2 * dipoleSep,
872 val.values,
873 atol=1e-16,
874 rtol=1e-16))
876 def _makeWcs(self, dec=53.1595451514076, theta=0):
877 """Create a wcs from real CFHT values, rotated by an optional theta.
879 dec : `float`
880 Set reference declination of CRVAL2 [degrees]
881 theta : `float`
882 Rotate CD matrix by theta [degrees]
884 Returns
885 -------
886 wcs : `lsst.afw.geom`
887 Created wcs.
888 """
889 metadata = dafBase.PropertySet()
891 metadata.set("SIMPLE", "T")
892 metadata.set("BITPIX", -32)
893 metadata.set("NAXIS", 2)
894 metadata.set("NAXIS1", 1024)
895 metadata.set("NAXIS2", 1153)
896 metadata.set("RADECSYS", 'FK5')
897 metadata.set("EQUINOX", 2000.)
899 metadata.setDouble("CRVAL1", 215.604025685476)
900 metadata.setDouble("CRVAL2", dec)
901 metadata.setDouble("CRPIX1", 1109.99981456774)
902 metadata.setDouble("CRPIX2", 560.018167811613)
903 metadata.set("CTYPE1", 'RA---SIN')
904 metadata.set("CTYPE2", 'DEC--SIN')
906 cd_matrix = np.array(
907 [
908 [5.10808596133527E-05, 1.85579539217196E-07],
909 [-8.27440751733828E-07, -5.10281493481982E-05]
910 ]
911 )
912 # rotate CD matrix
913 theta_rad = np.deg2rad(theta)
914 rotation_matrix = np.array(
915 [
916 [np.cos(theta_rad), -np.sin(theta_rad)],
917 [np.sin(theta_rad), np.cos(theta_rad)],
918 ]
919 )
920 cd_matrix = rotation_matrix @ cd_matrix
922 metadata.setDouble("CD1_1", cd_matrix[0, 0])
923 metadata.setDouble("CD1_2", cd_matrix[0, 1])
924 metadata.setDouble("CD2_1", cd_matrix[1, 0])
925 metadata.setDouble("CD2_2", cd_matrix[1, 1])
927 return afwGeom.makeSkyWcs(metadata)
929 def testHtm(self):
930 """Test that HtmIndxes are created as expected.
931 """
932 df = self.getMultiIndexDataFrame(self.dataDict)
933 func = HtmIndex20("coord_ra", "coord_dec")
935 val = self._funcVal(func, df)
936 # Test that the HtmIds come out as the ra/dec in dataDict.
937 self.assertTrue(np.all(np.equal(
938 val.values,
939 [14924528684992, 14924528689697, 14924528501716, 14924526434259,
940 14924526433879])))
942 def testEbv(self):
943 """Test that EBV works.
944 """
945 df = self.getMultiIndexDataFrame(self.dataDict)
946 func = Ebv()
948 val = self._funcVal(func, df)
949 np.testing.assert_array_almost_equal(
950 val.values,
951 [0.029100, 0.029013, 0.028857, 0.028802, 0.028797]
952 )
954 def testSkyMoments(self):
955 # TODO: We should vectorize the afwGeom functions for the conversions instead of just using
956 # them for tests: DM-54015
958 self.columns.extend([
959 "slot_Shape_xx",
960 "slot_Shape_yy",
961 "slot_Shape_xy",
962 "base_LocalWcs_CDMatrix_1_1",
963 "base_LocalWcs_CDMatrix_1_2",
964 "base_LocalWcs_CDMatrix_2_1",
965 "base_LocalWcs_CDMatrix_1_1",
966 ])
968 # CD Matrix from a ComCam exposure.
969 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = \
970 np.full(self.nRecords, -9.695088e-07)
971 self.dataDict["base_LocalWcs_CDMatrix_1_2"] = \
972 np.full(self.nRecords, 3.950301849959342e-09)
973 self.dataDict["base_LocalWcs_CDMatrix_2_1"] = \
974 np.full(self.nRecords, 3.8766915166433014e-09)
975 self.dataDict["base_LocalWcs_CDMatrix_2_2"] = \
976 np.full(self.nRecords, 9.695092484727074e-07)
977 self.dataDict["slot_Shape_xx"] = \
978 np.array([6.52675084, 74.17426471, 6.45283335, 36.82870958, 6.45685472])
979 self.dataDict["slot_Shape_yy"] = \
980 np.array([6.12848637, 80.99510036, 6.05671667, 35.79219613, 5.97778765])
981 self.dataDict["slot_Shape_xy"] = \
982 np.array([-0.10281872, 0.88788384, -0.1261287, -1.60504171, 0.11974515])
983 self.dataDict["slot_Shape_sigma_x"] = np.sqrt(self.dataDict["slot_Shape_xx"])
984 self.dataDict["slot_Shape_sigma_y"] = np.sqrt(self.dataDict["slot_Shape_yy"])
985 self.dataDict["slot_Shape_rho"] = self.dataDict["slot_Shape_xy"]/(
986 self.dataDict["slot_Shape_sigma_x"]*self.dataDict["slot_Shape_sigma_y"]
987 )
989 args_cd = [
990 "base_LocalWcs_CDMatrix_1_1", "base_LocalWcs_CDMatrix_1_2",
991 "base_LocalWcs_CDMatrix_2_1", "base_LocalWcs_CDMatrix_2_2",
992 ]
993 args = ["slot_Shape_xx", "slot_Shape_yy", "slot_Shape_xy"] + args_cd
994 args_corr = ["slot_Shape_sigma_x", "slot_Shape_sigma_y", "slot_Shape_rho"] + args_cd
996 skyXX_functor = MomentsIuuSky(*args, filt="g")
997 skyYY_functor = MomentsIvvSky(*args, filt="g")
998 skyXY_functor = MomentsIuvSky(*args, filt="g")
1000 axesA_functor = SemimajorAxisFromMoments(*args, filt="g")
1001 axesB_functor = SemiminorAxisFromMoments(*args, filt="g")
1002 axesTheta_functor = PositionAngleFromMoments(*args, filt="g")
1004 skyXX_corr_functor = CorrelationIuuSky(*args_corr, filt="g")
1005 skyYY_corr_functor = CorrelationIvvSky(*args_corr, filt="g")
1006 skyXY_corr_functor = CorrelationIuvSky(*args_corr, filt="g")
1008 axesA_corr_functor = SemimajorAxisFromCorrelation(*args_corr, filt="g")
1009 axesB_corr_functor = SemiminorAxisFromCorrelation(*args_corr, filt="g")
1010 axesTheta_corr_functor = PositionAngleFromCorrelation(*args_corr, filt="g")
1012 moments_g1_functor = MomentsG1Sky(*args, filt="g")
1013 moments_g2_functor = MomentsG2Sky(*args, filt="g")
1014 moments_trace_functor = MomentsTraceSky(*args, filt="g")
1016 df = self.getMultiIndexDataFrame(self.dataDict)
1017 output_sky_xx = skyXX_functor(df)
1018 output_sky_yy = skyYY_functor(df)
1019 output_sky_xy = skyXY_functor(df)
1021 output_axes_a = axesA_functor(df)
1022 output_axes_b = axesB_functor(df)
1023 output_axes_theta = axesTheta_functor(df)
1025 output_sky_xx_corr = skyXX_corr_functor(df)
1026 output_sky_yy_corr = skyYY_corr_functor(df)
1027 output_sky_xy_corr = skyXY_corr_functor(df)
1029 output_axes_a_corr = axesA_corr_functor(df)
1030 output_axes_b_corr = axesB_corr_functor(df)
1031 output_axes_theta_corr = axesTheta_corr_functor(df)
1033 output_g1 = moments_g1_functor(df)
1034 output_g2 = moments_g2_functor(df)
1035 output_trace = moments_trace_functor(df)
1037 transformed_xx = []
1038 transformed_yy = []
1039 transformed_xy = []
1040 axes_a = []
1041 axes_b = []
1042 axes_theta = []
1044 transformed_g1 = []
1045 transformed_g2 = []
1046 transformed_trace = []
1048 for n in range(5):
1049 Ixx = df[('meas', 'g', 'slot_Shape_xx')].iloc[n]
1050 Iyy = df[('meas', 'g', 'slot_Shape_yy')].iloc[n]
1051 Ixy = df[('meas', 'g', 'slot_Shape_xy')].iloc[n]
1052 localWCS_CD_1_1 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_1_1')].iloc[n]
1053 localWCS_CD_2_1 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_2_1')].iloc[n]
1054 localWCS_CD_1_2 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_1_2')].iloc[n]
1055 localWCS_CD_2_2 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_2_2')].iloc[n]
1056 CD_matrix = np.array([[localWCS_CD_1_1, localWCS_CD_1_2],
1057 [localWCS_CD_2_1, localWCS_CD_2_2]])
1059 q = afwGeom.ellipses.Quadrupole(Ixx, Iyy, Ixy)
1060 lt = geom.LinearTransform(CD_matrix)
1061 transformed_q = q.transform(lt)
1062 transformed_q.scale((180/np.pi) * (3600))
1064 axes = afwGeom.ellipses.Axes(transformed_q)
1066 transformed_xx.append(transformed_q.getIxx())
1067 transformed_yy.append(transformed_q.getIyy())
1068 transformed_xy.append(transformed_q.getIxy())
1070 axes_a.append(axes.getA())
1071 axes_b.append(axes.getB())
1072 axes_theta.append(np.degrees(axes.getTheta()))
1074 reduced_shear = afwGeom.ellipses.SeparableReducedShearTraceRadius(transformed_q)
1076 transformed_g1.append(reduced_shear.getE1())
1077 # Sign flip for consistency with HSM E2 sign convention.
1078 transformed_g2.append(-1*reduced_shear.getE2())
1079 transformed_trace.append(reduced_shear.getTraceRadius())
1081 self.assertFloatsAlmostEqual(output_sky_xx, np.array(transformed_xx), rtol=1e-7)
1082 self.assertFloatsAlmostEqual(output_sky_yy, np.array(transformed_yy), rtol=1e-7)
1083 self.assertFloatsAlmostEqual(output_sky_xy, np.array(transformed_xy), rtol=1e-7)
1084 self.assertFloatsAlmostEqual(output_sky_xx_corr, np.array(transformed_xx), rtol=1e-7)
1085 self.assertFloatsAlmostEqual(output_sky_yy_corr, np.array(transformed_yy), rtol=1e-7)
1086 self.assertFloatsAlmostEqual(output_sky_xy_corr, np.array(transformed_xy), rtol=1e-7)
1088 self.assertFloatsAlmostEqual(output_axes_a, np.array(axes_a), rtol=1e-7)
1089 self.assertFloatsAlmostEqual(output_axes_b, np.array(axes_b), rtol=1e-7)
1090 self.assertFloatsAlmostEqual(output_axes_theta, np.array(axes_theta), rtol=1e-7)
1091 self.assertFloatsAlmostEqual(output_axes_a_corr, np.array(axes_a), rtol=1e-7)
1092 self.assertFloatsAlmostEqual(output_axes_b_corr, np.array(axes_b), rtol=1e-7)
1093 self.assertFloatsAlmostEqual(output_axes_theta_corr, np.array(axes_theta), rtol=1e-7)
1095 self.assertFloatsAlmostEqual(output_g1, np.array(transformed_g1), rtol=1e-7)
1096 self.assertFloatsAlmostEqual(output_g2, np.array(transformed_g2), rtol=1e-7)
1097 self.assertFloatsAlmostEqual(output_trace, np.array(transformed_trace), rtol=2e-7)
1099 def _dropLevels(self, df):
1100 levelsToDrop = [n for lev, n in zip(df.columns.levels, df.columns.names) if len(lev) == 1]
1102 # Prevent error when trying to drop *all* columns
1103 if len(levelsToDrop) == len(df.columns.names):
1104 levelsToDrop.remove(df.columns.names[-1])
1106 df.columns = df.columns.droplevel(levelsToDrop)
1109class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase):
1110 pass
1113def setup_module(module):
1114 lsst.utils.tests.init()
1117if __name__ == "__main__": 1117 ↛ 1118line 1117 didn't jump to line 1118 because the condition on line 1117 was never true
1118 lsst.utils.tests.init()
1119 unittest.main()