Coverage for tests / test_functors.py: 10%
538 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-18 09:04 +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 )
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 def testConvertPixelToArcseconds(self):
767 """Test calculations of the pixel scale, conversions of pixel to
768 arcseconds.
769 """
770 dipoleSep = 10
771 ixx = 10
772 testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2))
773 localWcsPlugin = measBase.EvaluateLocalWcsPlugin(
774 None,
775 "base_LocalWcs",
776 afwTable.SourceTable.makeMinimalSchema(),
777 None)
778 for dec in np.linspace(-90, 90, 10):
779 for x, y in zip(np.random.uniform(2 * 1109.99981456774, size=10),
780 np.random.uniform(2 * 560.018167811613, size=10)):
781 center = geom.Point2D(x, y)
782 wcs = self._makeWcs(dec)
783 skyOrigin = wcs.pixelToSky(center)
785 linAffMatrix = localWcsPlugin.makeLocalTransformMatrix(wcs,
786 center)
787 self.dataDict["dipoleSep"] = np.full(self.nRecords, dipoleSep)
788 self.dataDict["ixx"] = np.full(self.nRecords, ixx)
789 self.dataDict["slot_Centroid_x"] = np.full(self.nRecords, x)
790 self.dataDict["slot_Centroid_y"] = np.full(self.nRecords, y)
791 self.dataDict["someCentroid_x"] = x + testPixelDeltas[:, 0]
792 self.dataDict["someCentroid_y"] = y + testPixelDeltas[:, 1]
794 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(self.nRecords,
795 linAffMatrix[0, 0])
796 self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(self.nRecords,
797 linAffMatrix[0, 1])
798 self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(self.nRecords,
799 linAffMatrix[1, 0])
800 self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(self.nRecords,
801 linAffMatrix[1, 1])
802 df = self.getMultiIndexDataFrame(self.dataDict)
803 func = LocalWcs("base_LocalWcs_CDMatrix_1_1",
804 "base_LocalWcs_CDMatrix_1_2",
805 "base_LocalWcs_CDMatrix_2_1",
806 "base_LocalWcs_CDMatrix_2_2")
808 # Exercise the full set of functions in LocalWcs.
809 sepRadians = func.getSkySeparationFromPixel(
810 df[("meas", "g", "someCentroid_x")] - df[("meas", "g", "slot_Centroid_x")],
811 df[("meas", "g", "someCentroid_y")] - df[("meas", "g", "slot_Centroid_y")],
812 0.0,
813 0.0,
814 df[("meas", "g", "base_LocalWcs_CDMatrix_1_1")],
815 df[("meas", "g", "base_LocalWcs_CDMatrix_1_2")],
816 df[("meas", "g", "base_LocalWcs_CDMatrix_2_1")],
817 df[("meas", "g", "base_LocalWcs_CDMatrix_2_2")])
819 # Test functor values against afw SkyWcs computations.
820 for centX, centY, sep in zip(testPixelDeltas[:, 0],
821 testPixelDeltas[:, 1],
822 sepRadians.values):
823 afwSepRadians = skyOrigin.separation(
824 wcs.pixelToSky(x + centX, y + centY)).asRadians()
825 self.assertAlmostEqual(1 - sep / afwSepRadians, 0, places=5)
827 # Test the pixel scale computation.
828 func = ComputePixelScale("base_LocalWcs_CDMatrix_1_1",
829 "base_LocalWcs_CDMatrix_1_2",
830 "base_LocalWcs_CDMatrix_2_1",
831 "base_LocalWcs_CDMatrix_2_2")
832 pixelScale = self._funcVal(func, df)
833 self.assertTrue(np.allclose(
834 wcs.getPixelScale(center).asArcseconds(),
835 pixelScale.values,
836 rtol=1e-6,
837 atol=0))
839 # Test pixel -> arcsec conversion.
840 func = ConvertPixelToArcseconds("dipoleSep",
841 "base_LocalWcs_CDMatrix_1_1",
842 "base_LocalWcs_CDMatrix_1_2",
843 "base_LocalWcs_CDMatrix_2_1",
844 "base_LocalWcs_CDMatrix_2_2")
845 val = self._funcVal(func, df)
846 self.assertTrue(np.allclose(pixelScale.values * dipoleSep,
847 val.values,
848 atol=1e-16,
849 rtol=1e-16))
851 # Test pixel^2 -> arcsec^2 conversion.
852 func = ConvertPixelSqToArcsecondsSq("ixx",
853 "base_LocalWcs_CDMatrix_1_1",
854 "base_LocalWcs_CDMatrix_1_2",
855 "base_LocalWcs_CDMatrix_2_1",
856 "base_LocalWcs_CDMatrix_2_2")
857 val = self._funcVal(func, df)
858 self.assertTrue(np.allclose(pixelScale.values ** 2 * dipoleSep,
859 val.values,
860 atol=1e-16,
861 rtol=1e-16))
863 def _makeWcs(self, dec=53.1595451514076, theta=0):
864 """Create a wcs from real CFHT values, rotated by an optional theta.
866 dec : `float`
867 Set reference declination of CRVAL2 [degrees]
868 theta : `float`
869 Rotate CD matrix by theta [degrees]
871 Returns
872 -------
873 wcs : `lsst.afw.geom`
874 Created wcs.
875 """
876 metadata = dafBase.PropertySet()
878 metadata.set("SIMPLE", "T")
879 metadata.set("BITPIX", -32)
880 metadata.set("NAXIS", 2)
881 metadata.set("NAXIS1", 1024)
882 metadata.set("NAXIS2", 1153)
883 metadata.set("RADECSYS", 'FK5')
884 metadata.set("EQUINOX", 2000.)
886 metadata.setDouble("CRVAL1", 215.604025685476)
887 metadata.setDouble("CRVAL2", dec)
888 metadata.setDouble("CRPIX1", 1109.99981456774)
889 metadata.setDouble("CRPIX2", 560.018167811613)
890 metadata.set("CTYPE1", 'RA---SIN')
891 metadata.set("CTYPE2", 'DEC--SIN')
893 cd_matrix = np.array(
894 [
895 [5.10808596133527E-05, 1.85579539217196E-07],
896 [-8.27440751733828E-07, -5.10281493481982E-05]
897 ]
898 )
899 # rotate CD matrix
900 theta_rad = np.deg2rad(theta)
901 rotation_matrix = np.array(
902 [
903 [np.cos(theta_rad), -np.sin(theta_rad)],
904 [np.sin(theta_rad), np.cos(theta_rad)],
905 ]
906 )
907 cd_matrix = rotation_matrix @ cd_matrix
909 metadata.setDouble("CD1_1", cd_matrix[0, 0])
910 metadata.setDouble("CD1_2", cd_matrix[0, 1])
911 metadata.setDouble("CD2_1", cd_matrix[1, 0])
912 metadata.setDouble("CD2_2", cd_matrix[1, 1])
914 return afwGeom.makeSkyWcs(metadata)
916 def testHtm(self):
917 """Test that HtmIndxes are created as expected.
918 """
919 df = self.getMultiIndexDataFrame(self.dataDict)
920 func = HtmIndex20("coord_ra", "coord_dec")
922 val = self._funcVal(func, df)
923 # Test that the HtmIds come out as the ra/dec in dataDict.
924 self.assertTrue(np.all(np.equal(
925 val.values,
926 [14924528684992, 14924528689697, 14924528501716, 14924526434259,
927 14924526433879])))
929 def testEbv(self):
930 """Test that EBV works.
931 """
932 df = self.getMultiIndexDataFrame(self.dataDict)
933 func = Ebv()
935 val = self._funcVal(func, df)
936 np.testing.assert_array_almost_equal(
937 val.values,
938 [0.029100, 0.029013, 0.028857, 0.028802, 0.028797]
939 )
941 def testSkyMoments(self):
942 # TODO: We should vectorize the afwGeom functions for the conversions instead of just using
943 # them for tests: DM-54015
945 self.columns.extend([
946 "slot_Shape_xx",
947 "slot_Shape_yy",
948 "slot_Shape_xy",
949 "base_LocalWcs_CDMatrix_1_1",
950 "base_LocalWcs_CDMatrix_1_2",
951 "base_LocalWcs_CDMatrix_2_1",
952 "base_LocalWcs_CDMatrix_1_1",
953 ])
955 # CD Matrix from a ComCam exposure.
956 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = \
957 np.full(self.nRecords, -9.695088e-07)
958 self.dataDict["base_LocalWcs_CDMatrix_1_2"] = \
959 np.full(self.nRecords, 3.950301849959342e-09)
960 self.dataDict["base_LocalWcs_CDMatrix_2_1"] = \
961 np.full(self.nRecords, 3.8766915166433014e-09)
962 self.dataDict["base_LocalWcs_CDMatrix_2_2"] = \
963 np.full(self.nRecords, 9.695092484727074e-07)
964 self.dataDict["slot_Shape_xx"] = \
965 np.array([6.52675084, 74.17426471, 6.45283335, 36.82870958, 6.45685472])
966 self.dataDict["slot_Shape_yy"] = \
967 np.array([6.12848637, 80.99510036, 6.05671667, 35.79219613, 5.97778765])
968 self.dataDict["slot_Shape_xy"] = \
969 np.array([-0.10281872, 0.88788384, -0.1261287, -1.60504171, 0.11974515])
970 self.dataDict["slot_Shape_sigma_x"] = np.sqrt(self.dataDict["slot_Shape_xx"])
971 self.dataDict["slot_Shape_sigma_y"] = np.sqrt(self.dataDict["slot_Shape_yy"])
972 self.dataDict["slot_Shape_rho"] = self.dataDict["slot_Shape_xy"]/(
973 self.dataDict["slot_Shape_sigma_x"]*self.dataDict["slot_Shape_sigma_y"]
974 )
976 args_cd = [
977 "base_LocalWcs_CDMatrix_1_1", "base_LocalWcs_CDMatrix_1_2",
978 "base_LocalWcs_CDMatrix_2_1", "base_LocalWcs_CDMatrix_2_2",
979 ]
980 args = ["slot_Shape_xx", "slot_Shape_yy", "slot_Shape_xy"] + args_cd
981 args_corr = ["slot_Shape_sigma_x", "slot_Shape_sigma_y", "slot_Shape_rho"] + args_cd
983 skyXX_functor = MomentsIuuSky(*args, filt="g")
984 skyYY_functor = MomentsIvvSky(*args, filt="g")
985 skyXY_functor = MomentsIuvSky(*args, filt="g")
987 axesA_functor = SemimajorAxisFromMoments(*args, filt="g")
988 axesB_functor = SemiminorAxisFromMoments(*args, filt="g")
989 axesTheta_functor = PositionAngleFromMoments(*args, filt="g")
991 skyXX_corr_functor = CorrelationIuuSky(*args_corr, filt="g")
992 skyYY_corr_functor = CorrelationIvvSky(*args_corr, filt="g")
993 skyXY_corr_functor = CorrelationIuvSky(*args_corr, filt="g")
995 axesA_corr_functor = SemimajorAxisFromCorrelation(*args_corr, filt="g")
996 axesB_corr_functor = SemiminorAxisFromCorrelation(*args_corr, filt="g")
997 axesTheta_corr_functor = PositionAngleFromCorrelation(*args_corr, filt="g")
999 moments_g1_functor = MomentsG1Sky(*args, filt="g")
1000 moments_g2_functor = MomentsG2Sky(*args, filt="g")
1001 moments_trace_functor = MomentsTraceSky(*args, filt="g")
1003 df = self.getMultiIndexDataFrame(self.dataDict)
1004 output_sky_xx = skyXX_functor(df)
1005 output_sky_yy = skyYY_functor(df)
1006 output_sky_xy = skyXY_functor(df)
1008 output_axes_a = axesA_functor(df)
1009 output_axes_b = axesB_functor(df)
1010 output_axes_theta = axesTheta_functor(df)
1012 output_sky_xx_corr = skyXX_corr_functor(df)
1013 output_sky_yy_corr = skyYY_corr_functor(df)
1014 output_sky_xy_corr = skyXY_corr_functor(df)
1016 output_axes_a_corr = axesA_corr_functor(df)
1017 output_axes_b_corr = axesB_corr_functor(df)
1018 output_axes_theta_corr = axesTheta_corr_functor(df)
1020 output_g1 = moments_g1_functor(df)
1021 output_g2 = moments_g2_functor(df)
1022 output_trace = moments_trace_functor(df)
1024 transformed_xx = []
1025 transformed_yy = []
1026 transformed_xy = []
1027 axes_a = []
1028 axes_b = []
1029 axes_theta = []
1031 transformed_g1 = []
1032 transformed_g2 = []
1033 transformed_trace = []
1035 for n in range(5):
1036 Ixx = df[('meas', 'g', 'slot_Shape_xx')].iloc[n]
1037 Iyy = df[('meas', 'g', 'slot_Shape_yy')].iloc[n]
1038 Ixy = df[('meas', 'g', 'slot_Shape_xy')].iloc[n]
1039 localWCS_CD_1_1 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_1_1')].iloc[n]
1040 localWCS_CD_2_1 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_2_1')].iloc[n]
1041 localWCS_CD_1_2 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_1_2')].iloc[n]
1042 localWCS_CD_2_2 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_2_2')].iloc[n]
1043 CD_matrix = np.array([[localWCS_CD_1_1, localWCS_CD_1_2],
1044 [localWCS_CD_2_1, localWCS_CD_2_2]])
1046 q = afwGeom.ellipses.Quadrupole(Ixx, Iyy, Ixy)
1047 lt = geom.LinearTransform(CD_matrix)
1048 transformed_q = q.transform(lt)
1049 transformed_q.scale((180/np.pi) * (3600))
1051 axes = afwGeom.ellipses.Axes(transformed_q)
1053 transformed_xx.append(transformed_q.getIxx())
1054 transformed_yy.append(transformed_q.getIyy())
1055 transformed_xy.append(transformed_q.getIxy())
1057 axes_a.append(axes.getA())
1058 axes_b.append(axes.getB())
1059 axes_theta.append(np.degrees(axes.getTheta()))
1061 reduced_shear = afwGeom.ellipses.SeparableReducedShearTraceRadius(transformed_q)
1063 transformed_g1.append(reduced_shear.getE1())
1064 # Sign flip for consistency with HSM E2 sign convention.
1065 transformed_g2.append(-1*reduced_shear.getE2())
1066 transformed_trace.append(reduced_shear.getTraceRadius())
1068 self.assertFloatsAlmostEqual(output_sky_xx, np.array(transformed_xx), rtol=1e-7)
1069 self.assertFloatsAlmostEqual(output_sky_yy, np.array(transformed_yy), rtol=1e-7)
1070 self.assertFloatsAlmostEqual(output_sky_xy, np.array(transformed_xy), rtol=1e-7)
1071 self.assertFloatsAlmostEqual(output_sky_xx_corr, np.array(transformed_xx), rtol=1e-7)
1072 self.assertFloatsAlmostEqual(output_sky_yy_corr, np.array(transformed_yy), rtol=1e-7)
1073 self.assertFloatsAlmostEqual(output_sky_xy_corr, np.array(transformed_xy), rtol=1e-7)
1075 self.assertFloatsAlmostEqual(output_axes_a, np.array(axes_a), rtol=1e-7)
1076 self.assertFloatsAlmostEqual(output_axes_b, np.array(axes_b), rtol=1e-7)
1077 self.assertFloatsAlmostEqual(output_axes_theta, np.array(axes_theta), rtol=1e-7)
1078 self.assertFloatsAlmostEqual(output_axes_a_corr, np.array(axes_a), rtol=1e-7)
1079 self.assertFloatsAlmostEqual(output_axes_b_corr, np.array(axes_b), rtol=1e-7)
1080 self.assertFloatsAlmostEqual(output_axes_theta_corr, np.array(axes_theta), rtol=1e-7)
1082 self.assertFloatsAlmostEqual(output_g1, np.array(transformed_g1), rtol=1e-7)
1083 self.assertFloatsAlmostEqual(output_g2, np.array(transformed_g2), rtol=1e-7)
1084 self.assertFloatsAlmostEqual(output_trace, np.array(transformed_trace), rtol=2e-7)
1086 def _dropLevels(self, df):
1087 levelsToDrop = [n for lev, n in zip(df.columns.levels, df.columns.names) if len(lev) == 1]
1089 # Prevent error when trying to drop *all* columns
1090 if len(levelsToDrop) == len(df.columns.names):
1091 levelsToDrop.remove(df.columns.names[-1])
1093 df.columns = df.columns.droplevel(levelsToDrop)
1096class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase):
1097 pass
1100def setup_module(module):
1101 lsst.utils.tests.init()
1104if __name__ == "__main__": 1104 ↛ 1105line 1104 didn't jump to line 1105 because the condition on line 1104 was never true
1105 lsst.utils.tests.init()
1106 unittest.main()