Coverage for tests / test_functors.py: 10%
522 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-20 01:17 -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/>.
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 CorrelationIuuSky, CorrelationIvvSky, CorrelationIuvSky,
54 SemimajorAxisFromCorrelation, SemiminorAxisFromCorrelation,
55 PositionAngleFromCorrelation
56 )
58ROOT = os.path.abspath(os.path.dirname(__file__))
61class FunctorTestCase(lsst.utils.tests.TestCase):
63 def getMultiIndexDataFrame(self, dataDict):
64 """Create a simple test multi-index DataFrame."""
66 simpleDF = pd.DataFrame(dataDict)
67 dfFilterDSCombos = []
68 for ds in self.datasets:
69 for band in self.bands:
70 df = copy.copy(simpleDF)
71 df.reindex(sorted(df.columns), axis=1)
72 df['dataset'] = ds
73 df['band'] = band
74 df.columns = pd.MultiIndex.from_tuples(
75 [(ds, band, c) for c in df.columns],
76 names=('dataset', 'band', 'column'))
77 dfFilterDSCombos.append(df)
79 df = functools.reduce(lambda d1, d2: d1.join(d2), dfFilterDSCombos)
81 return df
83 def getSimpleDataFrame(self, dataDict):
84 return pd.DataFrame(dataDict)
86 def getDatasetHandle(self, df):
87 lo, hi = HtmPixelization(7).universe().ranges()[0]
88 value = np.random.randint(lo, hi)
89 return InMemoryDatasetHandle(df, storageClass="DataFrame", dataId={"htm7": value})
91 def setUp(self):
92 np.random.seed(12345)
93 self.datasets = ['forced_src', 'meas', 'ref']
94 self.bands = ['g', 'r']
95 self.columns = ['coord_ra', 'coord_dec']
96 self.nRecords = 5
97 self.dataDict = {
98 "coord_ra": [3.77654137, 3.77643059, 3.77621148, 3.77611944, 3.77610396],
99 "coord_dec": [0.01127624, 0.01127787, 0.01127543, 0.01127543, 0.01127543]}
101 def _funcVal(self, functor, df):
102 self.assertIsInstance(functor.name, str)
103 self.assertIsInstance(functor.shortname, str)
105 handle = self.getDatasetHandle(df)
107 val = functor(df)
108 val2 = functor(handle)
109 self.assertTrue((val == val2).all())
110 self.assertIsInstance(val, pd.Series)
112 val = functor(df, dropna=True)
113 val2 = functor(handle, dropna=True)
114 self.assertTrue((val == val2).all())
115 self.assertEqual(val.isnull().sum(), 0)
117 return val
119 def _differenceVal(self, functor, df1, df2):
120 self.assertIsInstance(functor.name, str)
121 self.assertIsInstance(functor.shortname, str)
123 handle1 = self.getDatasetHandle(df1)
124 handle2 = self.getDatasetHandle(df2)
126 val = functor.difference(df1, df2)
127 val2 = functor.difference(handle1, handle2)
128 self.assertTrue(val.equals(val2))
129 self.assertIsInstance(val, pd.Series)
131 val = functor.difference(df1, df2, dropna=True)
132 val2 = functor.difference(handle1, handle2, dropna=True)
133 self.assertTrue(val.equals(val2))
134 self.assertEqual(val.isnull().sum(), 0)
136 val1 = self._funcVal(functor, df1)
137 val2 = self._funcVal(functor, df2)
139 self.assertTrue(np.allclose(val, val1 - val2))
141 return val
143 def testColumn(self):
144 self.columns.append("base_FootprintArea_value")
145 self.dataDict["base_FootprintArea_value"] = \
146 np.full(self.nRecords, 1)
147 df = self.getMultiIndexDataFrame(self.dataDict)
148 func = Column('base_FootprintArea_value', filt='g')
149 self._funcVal(func, df)
151 df = self.getSimpleDataFrame(self.dataDict)
152 func = Column('base_FootprintArea_value')
153 self._funcVal(func, df)
155 def testCustom(self):
156 self.columns.append("base_FootprintArea_value")
157 self.dataDict["base_FootprintArea_value"] = \
158 np.random.rand(self.nRecords)
159 df = self.getMultiIndexDataFrame(self.dataDict)
160 func = CustomFunctor('2*base_FootprintArea_value', filt='g')
161 val = self._funcVal(func, df)
163 func2 = Column('base_FootprintArea_value', filt='g')
165 np.allclose(val.values, 2*func2(df).values, atol=1e-13, rtol=0)
167 df = self.getSimpleDataFrame(self.dataDict)
168 func = CustomFunctor('2 * base_FootprintArea_value')
169 val = self._funcVal(func, df)
170 func2 = Column('base_FootprintArea_value')
172 np.allclose(val.values, 2*func2(df).values, atol=1e-13, rtol=0)
174 def testCoords(self):
175 df = self.getMultiIndexDataFrame(self.dataDict)
176 ra = self._funcVal(RAColumn(), df)
177 dec = self._funcVal(DecColumn(), df)
179 columnDict = {'dataset': 'ref', 'band': 'g',
180 'column': ['coord_ra', 'coord_dec']}
182 handle = InMemoryDatasetHandle(df, storageClass="DataFrame")
183 dfSub = handle.get(parameters={"columns": columnDict})
184 self._dropLevels(dfSub)
186 coords = dfSub / np.pi * 180.
188 self.assertTrue(np.allclose(ra, coords[('ref', 'g', 'coord_ra')], atol=1e-13, rtol=0))
189 self.assertTrue(np.allclose(dec, coords[('ref', 'g', 'coord_dec')], atol=1e-13, rtol=0))
191 # single-level column index table
192 df = self.getSimpleDataFrame(self.dataDict)
193 ra = self._funcVal(RAColumn(), df)
194 dec = self._funcVal(DecColumn(), df)
196 handle = InMemoryDatasetHandle(df, storageClass="DataFrame")
197 dfSub = handle.get(parameters={"columns": ['coord_ra', 'coord_dec']})
198 coords = dfSub / np.pi * 180.
200 self.assertTrue(np.allclose(ra, coords['coord_ra'], atol=1e-13, rtol=0))
201 self.assertTrue(np.allclose(dec, coords['coord_dec'], atol=1e-13, rtol=0))
203 def testMag(self):
204 self.columns.extend(["base_PsfFlux_instFlux", "base_PsfFlux_instFluxErr"])
205 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
206 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10)
207 df = self.getMultiIndexDataFrame(self.dataDict)
208 # Change one dataset filter combinations value.
209 df[("meas", "g", "base_PsfFlux_instFlux")] -= 1
211 fluxName = 'base_PsfFlux'
213 # Check that things work when you provide dataset explicitly
214 for dataset in ['forced_src', 'meas']:
215 psfMag_G = self._funcVal(Mag(fluxName, dataset=dataset,
216 filt='g'),
217 df)
218 psfMag_R = self._funcVal(Mag(fluxName, dataset=dataset,
219 filt='r'),
220 df)
222 psfColor_GR = self._funcVal(Color(fluxName, 'g', 'r',
223 dataset=dataset),
224 df)
226 self.assertTrue(np.allclose((psfMag_G - psfMag_R).dropna(), psfColor_GR, rtol=0, atol=1e-13))
228 # Check that behavior as expected when dataset not provided;
229 # that is, that the color comes from forced and default Mag is meas
230 psfMag_G = self._funcVal(Mag(fluxName, filt='g'), df)
231 psfMag_R = self._funcVal(Mag(fluxName, filt='r'), df)
233 psfColor_GR = self._funcVal(Color(fluxName, 'g', 'r'), df)
235 # These should *not* be equal.
236 self.assertFalse(np.allclose((psfMag_G - psfMag_R).dropna(), psfColor_GR))
238 def testMagDiff(self):
239 self.columns.extend(["base_PsfFlux_instFlux", "base_PsfFlux_instFluxErr",
240 "modelfit_CModel_instFlux", "modelfit_CModel_instFluxErr"])
241 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
242 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10)
243 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1000)
244 self.dataDict["modelfit_CModel_instFluxErr"] = np.full(self.nRecords, 10)
245 df = self.getMultiIndexDataFrame(self.dataDict)
247 for filt in self.bands:
248 filt = 'g'
249 val = self._funcVal(MagDiff('base_PsfFlux', 'modelfit_CModel', filt=filt), df)
251 mag1 = self._funcVal(Mag('modelfit_CModel', filt=filt), df)
252 mag2 = self._funcVal(Mag('base_PsfFlux', filt=filt), df)
253 self.assertTrue(np.allclose((mag2 - mag1).dropna(), val, rtol=0, atol=1e-13))
255 def testDifference(self):
256 """Test .difference method using MagDiff as the example.
257 """
258 self.columns.extend(["base_PsfFlux_instFlux", "base_PsfFlux_instFluxErr",
259 "modelfit_CModel_instFlux", "modelfit_CModel_instFluxErr"])
261 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
262 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1000)
263 df1 = self.getMultiIndexDataFrame(self.dataDict)
265 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 999)
266 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 999)
267 df2 = self.getMultiIndexDataFrame(self.dataDict)
269 magDiff = MagDiff('base_PsfFlux', 'modelfit_CModel', filt='g')
271 # Asserts that differences computed properly
272 self._differenceVal(magDiff, df1, df2)
274 def testPixelScale(self):
275 # Test that the pixel scale and pix->arcsec calculations perform as
276 # expected.
277 pass
279 def testShape(self):
280 data = {
281 "x": np.array([-0.3, 0.4, 0.7, -0.9, 1.4, -5.3]),
282 "y": np.array([1.5, -0.7, -1.9, 2.8, -1.4, 0.01]),
283 "rho": np.array([-0.9, 0.4, -0.7, 0., 0.3, -0.99]),
284 }
285 data["xx"] = data["x"]**2
286 data["yy"] = data["y"]**2
287 data["xy"] = data["x"]*data["y"]*data["rho"]
289 args = ("xx", "xy", "yy")
290 functor_e1, functor_e2, functor_quadrupole = E1(*args), E2(*args), RadiusFromQuadrupole(*args)
292 xx_plus_yy = data["xx"] + data["yy"]
293 data = pd.DataFrame(data)
295 np.testing.assert_allclose(
296 functor_e1(data),
297 ((data["xx"] - data["yy"])/xx_plus_yy).astype(np.float32),
298 rtol=1e-12, atol=1e-12,
299 )
300 np.testing.assert_allclose(
301 functor_e2(data),
302 (2.0*data["xy"]/xx_plus_yy).astype(np.float32),
303 rtol=1e-12, atol=1e-12,
304 )
305 np.testing.assert_allclose(
306 functor_quadrupole(data),
307 ((data["xx"]*data["yy"] - data["xy"]**2)**0.25).astype(np.float32),
308 rtol=1e-12, atol=1e-12,
309 )
311 def testOther(self):
312 self.columns.extend(["ext_shapeHSM_HsmSourceMoments_xx", "ext_shapeHSM_HsmSourceMoments_yy",
313 "base_SdssShape_xx", "base_SdssShape_yy",
314 "ext_shapeHSM_HsmPsfMoments_xx", "ext_shapeHSM_HsmPsfMoments_yy",
315 "base_SdssShape_psf_xx", "base_SdssShape_psf_yy"])
316 self.dataDict["ext_shapeHSM_HsmSourceMoments_xx"] = np.full(self.nRecords, 1 / np.sqrt(2))
317 self.dataDict["ext_shapeHSM_HsmSourceMoments_yy"] = np.full(self.nRecords, 1 / np.sqrt(2))
318 self.dataDict["base_SdssShape_xx"] = np.full(self.nRecords, 1 / np.sqrt(2))
319 self.dataDict["base_SdssShape_yy"] = np.full(self.nRecords, 1 / np.sqrt(2))
320 self.dataDict["ext_shapeHSM_HsmPsfMoments_xx"] = np.full(self.nRecords, 1 / np.sqrt(2))
321 self.dataDict["ext_shapeHSM_HsmPsfMoments_yy"] = np.full(self.nRecords, 1 / np.sqrt(2))
322 self.dataDict["base_SdssShape_psf_xx"] = np.full(self.nRecords, 1)
323 self.dataDict["base_SdssShape_psf_yy"] = np.full(self.nRecords, 1)
324 df = self.getMultiIndexDataFrame(self.dataDict)
325 # Covering the code is better than nothing
326 for filt in self.bands:
327 for Func in [DeconvolvedMoments,
328 SdssTraceSize,
329 PsfSdssTraceSizeDiff,
330 HsmTraceSize, PsfHsmTraceSizeDiff, HsmFwhm]:
331 _ = self._funcVal(Func(filt=filt), df)
333 def _compositeFuncVal(self, functor, df):
334 self.assertIsInstance(functor, CompositeFunctor)
336 handle = self.getDatasetHandle(df)
338 fdf1 = functor(df)
339 fdf2 = functor(handle)
340 self.assertTrue(fdf1.equals(fdf2))
342 self.assertIsInstance(fdf1, pd.DataFrame)
343 self.assertTrue(np.all([k in fdf1.columns for k in functor.funcDict.keys()]))
345 fdf1 = functor(df, dropna=True)
346 fdf2 = functor(handle, dropna=True)
347 self.assertTrue(fdf1.equals(fdf2))
349 # Check that there are no nulls
350 self.assertFalse(fdf1.isnull().any(axis=None))
352 return fdf1
354 def _compositeDifferenceVal(self, functor, df1, df2):
355 self.assertIsInstance(functor, CompositeFunctor)
357 handle1 = self.getDatasetHandle(df1)
358 handle2 = self.getDatasetHandle(df2)
360 fdf1 = functor.difference(df1, df2)
361 fdf2 = functor.difference(handle1, handle2)
362 self.assertTrue(fdf1.equals(fdf2))
364 self.assertIsInstance(fdf1, pd.DataFrame)
365 self.assertTrue(np.all([k in fdf1.columns for k in functor.funcDict.keys()]))
367 fdf1 = functor.difference(df1, df2, dropna=True)
368 fdf2 = functor.difference(handle1, handle2, dropna=True)
369 self.assertTrue(fdf1.equals(fdf2))
371 # Check that there are no nulls
372 self.assertFalse(fdf1.isnull().any(axis=None))
374 df1_functored = functor(df1)
375 df2_functored = functor(df2)
377 self.assertTrue(np.allclose(fdf1.values, df1_functored.values - df2_functored.values))
379 return fdf1
381 def testComposite(self):
382 self.columns.extend(["modelfit_CModel_instFlux", "base_PsfFlux_instFlux"])
383 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1)
384 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1)
386 df = self.getMultiIndexDataFrame(self.dataDict)
387 # Modify r band value slightly.
388 df[("meas", "r", "base_PsfFlux_instFlux")] -= 0.1
390 filt = 'g'
391 funcDict = {'psfMag_ref': Mag('base_PsfFlux', dataset='ref'),
392 'ra': RAColumn(),
393 'dec': DecColumn(),
394 'psfMag': Mag('base_PsfFlux', filt=filt),
395 'cmodel_magDiff': MagDiff('base_PsfFlux',
396 'modelfit_CModel', filt=filt)}
397 func = CompositeFunctor(funcDict)
398 fdf1 = self._compositeFuncVal(func, df)
400 # Repeat same, but define filter globally instead of individually
401 funcDict2 = {'psfMag_ref': Mag('base_PsfFlux', dataset='ref'),
402 'ra': RAColumn(),
403 'dec': DecColumn(),
404 'psfMag': Mag('base_PsfFlux'),
405 'cmodel_magDiff': MagDiff('base_PsfFlux',
406 'modelfit_CModel')}
408 func2 = CompositeFunctor(funcDict2, filt=filt)
409 fdf2 = self._compositeFuncVal(func2, df)
410 self.assertTrue(fdf1.equals(fdf2))
412 func2.filt = 'r'
413 fdf3 = self._compositeFuncVal(func2, df)
414 # Because we modified the R filter this should fail.
415 self.assertFalse(fdf2.equals(fdf3))
417 # Make sure things work with passing list instead of dict
418 funcs = [Mag('base_PsfFlux', dataset='ref'),
419 RAColumn(),
420 DecColumn(),
421 Mag('base_PsfFlux', filt=filt),
422 MagDiff('base_PsfFlux', 'modelfit_CModel', filt=filt)]
424 _ = self._compositeFuncVal(CompositeFunctor(funcs), df)
426 def testCompositeSimple(self):
427 """Test single-level composite functor for functionality
428 """
429 self.columns.extend(["modelfit_CModel_instFlux", "base_PsfFlux_instFlux"])
430 self.dataDict["modelfit_CModel_instFlux"] = np.full(self.nRecords, 1)
431 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1)
433 df = self.getSimpleDataFrame(self.dataDict)
434 funcDict = {'ra': RAColumn(),
435 'dec': DecColumn(),
436 'psfMag': Mag('base_PsfFlux'),
437 'cmodel_magDiff': MagDiff('base_PsfFlux',
438 'modelfit_CModel')}
439 func = CompositeFunctor(funcDict)
440 _ = self._compositeFuncVal(func, df)
442 def testCompositeColor(self):
443 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
444 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10)
445 df = self.getMultiIndexDataFrame(self.dataDict)
446 funcDict = {'a': Mag('base_PsfFlux', dataset='meas', filt='g'),
447 'b': Mag('base_PsfFlux', dataset='forced_src', filt='g'),
448 'c': Color('base_PsfFlux', 'g', 'r')}
449 # Covering the code is better than nothing
450 _ = self._compositeFuncVal(CompositeFunctor(funcDict), df)
452 def testCompositeDifference(self):
453 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
454 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 10)
455 df1 = self.getMultiIndexDataFrame(self.dataDict)
457 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 999)
458 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords, 9)
459 df2 = self.getMultiIndexDataFrame(self.dataDict)
461 funcDict = {'a': Mag('base_PsfFlux', dataset='meas', filt='g'),
462 'b': Mag('base_PsfFlux', dataset='forced_src', filt='g'),
463 'c': Color('base_PsfFlux', 'g', 'r')}
464 # Covering the code is better than nothing
465 _ = self._compositeDifferenceVal(CompositeFunctor(funcDict), df1, df2)
467 def testCompositeFail(self):
468 """Test a composite functor where one of the functors should be junk.
469 """
470 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, 1000)
471 df = self.getMultiIndexDataFrame(self.dataDict)
473 funcDict = {'good': Column("base_PsfFlux_instFlux"),
474 'bad': Column('not_a_column')}
476 with self.assertLogs(level=logging.ERROR) as cm:
477 _ = self._compositeFuncVal(CompositeFunctor(funcDict), df)
478 self.assertIn("Exception in CompositeFunctor (funcs: ['good', 'bad'])", cm.output[0])
480 def testLocalPhotometry(self):
481 """Test the local photometry functors.
482 """
483 flux = 1000
484 fluxErr = 10
485 calib = 10
486 calibErr = 0.0
487 self.dataDict["base_PsfFlux_instFlux"] = np.full(self.nRecords, flux)
488 self.dataDict["base_PsfFlux_instFluxErr"] = np.full(self.nRecords,
489 fluxErr)
490 self.dataDict["base_LocalPhotoCalib"] = np.full(self.nRecords, calib)
492 df = self.getMultiIndexDataFrame(self.dataDict)
493 func = LocalPhotometry("base_PsfFlux_instFlux",
494 "base_PsfFlux_instFluxErr",
495 "base_LocalPhotoCalib")
497 nanoJansky = func.instFluxToNanojansky(
498 df[("meas", "g", "base_PsfFlux_instFlux")],
499 df[("meas", "g", "base_LocalPhotoCalib")])
500 mag = func.instFluxToMagnitude(
501 df[("meas", "g", "base_PsfFlux_instFlux")],
502 df[("meas", "g", "base_LocalPhotoCalib")])
503 nanoJanskyErr = func.instFluxErrToNanojanskyErr(
504 df[("meas", "g", "base_PsfFlux_instFlux")],
505 df[("meas", "g", "base_PsfFlux_instFluxErr")],
506 df[("meas", "g", "base_LocalPhotoCalib")])
507 magErr = func.instFluxErrToMagnitudeErr(
508 df[("meas", "g", "base_PsfFlux_instFlux")],
509 df[("meas", "g", "base_PsfFlux_instFluxErr")],
510 df[("meas", "g", "base_LocalPhotoCalib")])
512 self.assertTrue(np.allclose(nanoJansky.values,
513 flux * calib,
514 atol=0,
515 rtol=1e-7))
516 self.assertTrue(np.allclose(mag.values,
517 (flux * calib * u.nJy).to_value(u.ABmag),
518 atol=0,
519 rtol=1e-7))
520 self.assertTrue(np.allclose(nanoJanskyErr.values,
521 np.hypot(fluxErr * calib, flux * calibErr),
522 atol=0,
523 rtol=1e-7))
524 self.assertTrue(np.allclose(
525 magErr.values,
526 2.5 / np.log(10) * nanoJanskyErr.values / nanoJansky.values,
527 atol=0,
528 rtol=1e-7))
530 # Test functors against the values computed above.
531 self._testLocalPhotometryFunctors(LocalNanojansky,
532 df,
533 nanoJansky)
534 self._testLocalPhotometryFunctors(LocalNanojanskyErr,
535 df,
536 nanoJanskyErr)
538 def _testLocalPhotometryFunctors(self, functor, df, testValues):
539 func = functor("base_PsfFlux_instFlux",
540 "base_PsfFlux_instFluxErr",
541 "base_LocalPhotoCalib")
542 val = self._funcVal(func, df)
543 self.assertTrue(np.allclose(testValues.values,
544 val.values,
545 atol=0,
546 rtol=1e-7))
548 def testDipPhotometry(self):
549 """Test calibrated flux calculations for dipoles."""
550 fluxNeg = -100
551 fluxPos = 101
552 fluxErr = 10
553 calib = 10
554 calibErr = 0.0
556 # compute expected values.
557 absMean = 0.5*(fluxPos - fluxNeg)*calib
558 absDiff = (fluxNeg + fluxPos)*calib
559 absMeanErr = 0.5*np.sqrt(2*(fluxErr*calib)**2
560 + ((fluxPos - fluxNeg)*calibErr)**2)
561 absDiffErr = np.sqrt(2*(fluxErr*calib)**2
562 + ((fluxPos + fluxNeg)*calibErr)**2)
564 self.dataDict["ip_diffim_DipoleFluxNeg_instFlux"] = np.full(self.nRecords, fluxNeg)
565 self.dataDict["ip_diffim_DipoleFluxNeg_instFluxErr"] = np.full(self.nRecords, fluxErr)
566 self.dataDict["ip_diffim_DipoleFluxPos_instFlux"] = np.full(self.nRecords, fluxPos)
567 self.dataDict["ip_diffim_DipoleFluxPos_instFluxErr"] = np.full(self.nRecords, fluxErr)
568 self.dataDict["base_LocalPhotoCalib"] = np.full(self.nRecords, calib)
570 df = self.getMultiIndexDataFrame(self.dataDict)
571 func = LocalDipoleMeanFlux("ip_diffim_DipoleFluxPos_instFlux",
572 "ip_diffim_DipoleFluxNeg_instFlux",
573 "ip_diffim_DipoleFluxPos_instFluxErr",
574 "ip_diffim_DipoleFluxNeg_instFluxErr",
575 "base_LocalPhotoCalib")
576 val = self._funcVal(func, df)
577 self.assertFloatsAlmostEqual(val.values,
578 absMean,
579 atol=1e-13,
580 rtol=0)
582 func = LocalDipoleMeanFluxErr("ip_diffim_DipoleFluxPos_instFlux",
583 "ip_diffim_DipoleFluxNeg_instFlux",
584 "ip_diffim_DipoleFluxPos_instFluxErr",
585 "ip_diffim_DipoleFluxNeg_instFluxErr",
586 "base_LocalPhotoCalib")
587 val = self._funcVal(func, df)
588 self.assertFloatsAlmostEqual(val.values,
589 absMeanErr,
590 atol=1e-13,
591 rtol=0)
593 func = LocalDipoleDiffFlux("ip_diffim_DipoleFluxPos_instFlux",
594 "ip_diffim_DipoleFluxNeg_instFlux",
595 "ip_diffim_DipoleFluxPos_instFluxErr",
596 "ip_diffim_DipoleFluxNeg_instFluxErr",
597 "base_LocalPhotoCalib")
598 val = self._funcVal(func, df)
599 self.assertFloatsAlmostEqual(val.values,
600 absDiff,
601 atol=1e-13,
602 rtol=0)
604 func = LocalDipoleDiffFluxErr("ip_diffim_DipoleFluxPos_instFlux",
605 "ip_diffim_DipoleFluxNeg_instFlux",
606 "ip_diffim_DipoleFluxPos_instFluxErr",
607 "ip_diffim_DipoleFluxNeg_instFluxErr",
608 "base_LocalPhotoCalib")
609 val = self._funcVal(func, df)
610 self.assertFloatsAlmostEqual(val.values,
611 absDiffErr,
612 atol=1e-13,
613 rtol=0)
615 def testComputePositionAngle(self, offset=0.0001):
616 """Test computation of position angle from (RA1, Dec1) to (RA2, Dec2)
618 offset : `float`
619 Arc length of the offset vector to set up test points. [radian]
620 """
622 d = offset
623 columns = ("ra1", "dec1", "ra2", "dec2", "expected")
624 position_angle_test_values = [
625 # Get 0, 0 right
626 (0, 0, d, 0, np.pi/2),
627 (0, 0, 0, d, 0),
628 (0, 0, -d, 0, -np.pi/2),
629 (0, 0, 0, -d, np.pi),
630 # Make sure we get wrap-around to 0, 0 right
631 (2*np.pi, 0, d, 0, np.pi/2),
632 (2*np.pi, 0, 0, d, 0),
633 (2*np.pi, 0, -d, 0, -np.pi/2),
634 (2*np.pi, 0, 0, -d, np.pi),
635 (+0.0015, 0, 2*np.pi - 0.05, 0, -np.pi/2),
636 # Get another somewhat arbitrary location right [these are in rad]
637 # It's not really important to rescale d by 1/cos(dec)
638 # because we're just looking at orientation of vector
639 # but foreshadowing to the poles...
640 (0.0015, 1, 0.0015 + d*np.cos(-1), 1, np.pi/2),
641 (0.0015, 1, 0.0015, 1 + d, 0),
642 (0.0015, 1, 0.0015 - d*np.cos(-1), 1, - np.pi/2),
643 (0.0015, 1, 0.0015, 1 - d, np.pi),
644 # Negative dec
645 (0.0015, -1, 0.0015 + d*np.cos(-1), -1, np.pi/2),
646 (0.0015, -1, 0.0015, -1 + d, 0),
647 (0.0015, -1, 0.0015 - d*np.cos(-1), -1, - np.pi/2),
648 (0.0015, -1, 0.0015, -1 - d, np.pi),
649 # Make sure we get wrap-around on that right
650 (2*np.pi + 0.0015, 1, 2*np.pi + 0.0015 + d*np.cos(1), 1, np.pi/2),
651 (2*np.pi + 0.0015, 1, 2*np.pi + 0.0015, 1 + d, 0),
652 (2*np.pi + 0.0015, 1, 0.0015 - d*np.cos(1), 1, - np.pi/2),
653 (2*np.pi + 0.0015, 1, 0.0015, 1 - d, np.pi),
654 # Get relative wrap-around right
655 (2*np.pi + 0.0015, 1, 0.0015 + d*np.cos(1), 1, np.pi/2),
656 (0.0015, 1, 2*np.pi + 0.0015, 1 + d, 0),
657 (0.0015, 1, 2*np.pi + 0.0015 - d*np.cos(1), 1, - np.pi/2),
658 (2*np.pi + 0.0015, 1, 0.0015, 1 - d, np.pi),
659 # Get either side of RA=0 right
660 (+ d*np.cos(1) / 2, 1, 2*np.pi - d*np.cos(1) / 2, 1, - np.pi/2),
661 (+ d*np.cos(1) / 2, 1, - d*np.cos(1) / 2, 1, -np.pi/2),
662 (2*np.pi + d*np.cos(1) / 2, 1, - d*np.cos(1) / 2, 1, -np.pi/2),
663 (-d*np.cos(1) / 2, 1, +0.0015, 1, np.pi/2),
664 (0.0015, 1, 0.0015, 1 + d, 0),
665 (0.0015, 1, 0.0015 - d*np.cos(1), 1, - np.pi/2),
666 (0.0015, 1, 0.0015, 1 - d, np.pi),
667 # Try it near the poles
668 (0, np.pi/2, 0, np.pi/2 - d, np.pi),
669 (0, np.pi/2 - d, 0, np.pi/2, 0),
670 (0, -np.pi/2, 0, -np.pi/2 + d, 0),
671 (0, -np.pi/2 + d, 0, -np.pi/2, np.pi),
672 ]
674 df = pd.DataFrame(position_angle_test_values, columns=columns)
676 cd_matrix = np.array([[1, 0], [0, -1]]) # Doesn't matter because we don't use it.
677 local_wcs = LocalWcs(*cd_matrix.flatten())
678 pa = local_wcs.computePositionAngle(df["ra1"], df["dec1"], df["ra2"], df["dec2"])
679 expected = df["expected"]
680 tolerance_deg = 0.05 # degrees
681 tolerance_rad = np.deg2rad(tolerance_deg)
683 # Use SphereGeom to handle wrap-around separations correctly.
684 diff = [
685 geom.Angle(o, geom.radians).separation(geom.Angle(e, geom.radians)).asRadians()
686 for o, e in zip(pa, expected)
687 ]
689 np.testing.assert_allclose(diff, 0, rtol=0, atol=tolerance_rad)
691 # Test position angle
692 def testConvertDetectorAngleToPositionAngle(self):
693 """Test conversion of position angle in detector degrees to position angle on sky
695 There is overlap with testConvertPixelToArcseconds
696 But we also test additional rotation angles so this is separate.
697 """
698 dipoleSep = 10
699 ixx = 10
700 testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2))
701 # testConvertPixelToArcSecond uses a meas_base LocalWcsPlugin
702 # but we're using a simple WCS as our example, so this doesn't really matter
703 # and we'll just use the WCS directly
704 for dec in np.linspace(-90, 90, 10):
705 for theta in (-45, 0, 90):
706 for x, y in zip(np.random.uniform(2 * 1109.99981456774, size=10),
707 np.random.uniform(2 * 560.018167811613, size=10)):
708 wcs = self._makeWcs(dec=dec, theta=theta)
709 cd_matrix = wcs.getCdMatrix()
711 self.dataDict["dipoleSep"] = np.full(self.nRecords, dipoleSep)
712 self.dataDict["ixx"] = np.full(self.nRecords, ixx)
713 self.dataDict["slot_Centroid_x"] = np.full(self.nRecords, x)
714 self.dataDict["slot_Centroid_y"] = np.full(self.nRecords, y)
715 self.dataDict["someCentroid_x"] = x + testPixelDeltas[:, 0]
716 self.dataDict["someCentroid_y"] = y + testPixelDeltas[:, 1]
717 self.dataDict["orientation"] = np.arctan2(
718 self.dataDict["someCentroid_y"] - self.dataDict["slot_Centroid_y"],
719 self.dataDict["someCentroid_x"] - self.dataDict["slot_Centroid_x"],
720 )
722 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(self.nRecords,
723 cd_matrix[0, 0])
724 self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(self.nRecords,
725 cd_matrix[0, 1])
726 self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(self.nRecords,
727 cd_matrix[1, 0])
728 self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(self.nRecords,
729 cd_matrix[1, 1])
730 df = self.getMultiIndexDataFrame(self.dataDict)
732 # Test detector angle to position angle conversion
733 func = ConvertDetectorAngleToPositionAngle(
734 "orientation",
735 "base_LocalWcs_CDMatrix_1_1",
736 "base_LocalWcs_CDMatrix_1_2",
737 "base_LocalWcs_CDMatrix_2_1",
738 "base_LocalWcs_CDMatrix_2_2"
739 )
740 val = self._funcVal(func, df)
742 delta_ra, delta_dec = func.computeDeltaRaDec(
743 self.dataDict["someCentroid_x"] - self.dataDict["slot_Centroid_x"],
744 self.dataDict["someCentroid_y"] - self.dataDict["slot_Centroid_y"],
745 self.dataDict["base_LocalWcs_CDMatrix_1_1"],
746 self.dataDict["base_LocalWcs_CDMatrix_1_2"],
747 self.dataDict["base_LocalWcs_CDMatrix_2_1"],
748 self.dataDict["base_LocalWcs_CDMatrix_2_2"],
749 )
751 dx = np.cos(0) * np.tan(delta_dec) - np.sin(0) * np.cos(delta_ra)
752 dy = np.sin(delta_ra)
753 comparison_pa = np.arctan2(dy, dx)
754 comparison_pa = np.rad2deg(comparison_pa)
756 coord_diff = []
757 for v, c in zip(val.values, comparison_pa):
758 observed_angle = geom.Angle(v, geom.degrees)
759 expected_angle = geom.Angle(c, geom.degrees)
760 diff = observed_angle.separation(expected_angle).asRadians()
761 coord_diff.append(diff)
763 np.testing.assert_allclose(coord_diff, 0, rtol=0, atol=5e-6)
765 def testConvertPixelToArcseconds(self):
766 """Test calculations of the pixel scale, conversions of pixel to
767 arcseconds.
768 """
769 dipoleSep = 10
770 ixx = 10
771 testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2))
772 localWcsPlugin = measBase.EvaluateLocalWcsPlugin(
773 None,
774 "base_LocalWcs",
775 afwTable.SourceTable.makeMinimalSchema(),
776 None)
777 for dec in np.linspace(-90, 90, 10):
778 for x, y in zip(np.random.uniform(2 * 1109.99981456774, size=10),
779 np.random.uniform(2 * 560.018167811613, size=10)):
780 center = geom.Point2D(x, y)
781 wcs = self._makeWcs(dec)
782 skyOrigin = wcs.pixelToSky(center)
784 linAffMatrix = localWcsPlugin.makeLocalTransformMatrix(wcs,
785 center)
786 self.dataDict["dipoleSep"] = np.full(self.nRecords, dipoleSep)
787 self.dataDict["ixx"] = np.full(self.nRecords, ixx)
788 self.dataDict["slot_Centroid_x"] = np.full(self.nRecords, x)
789 self.dataDict["slot_Centroid_y"] = np.full(self.nRecords, y)
790 self.dataDict["someCentroid_x"] = x + testPixelDeltas[:, 0]
791 self.dataDict["someCentroid_y"] = y + testPixelDeltas[:, 1]
793 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(self.nRecords,
794 linAffMatrix[0, 0])
795 self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(self.nRecords,
796 linAffMatrix[0, 1])
797 self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(self.nRecords,
798 linAffMatrix[1, 0])
799 self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(self.nRecords,
800 linAffMatrix[1, 1])
801 df = self.getMultiIndexDataFrame(self.dataDict)
802 func = LocalWcs("base_LocalWcs_CDMatrix_1_1",
803 "base_LocalWcs_CDMatrix_1_2",
804 "base_LocalWcs_CDMatrix_2_1",
805 "base_LocalWcs_CDMatrix_2_2")
807 # Exercise the full set of functions in LocalWcs.
808 sepRadians = func.getSkySeparationFromPixel(
809 df[("meas", "g", "someCentroid_x")] - df[("meas", "g", "slot_Centroid_x")],
810 df[("meas", "g", "someCentroid_y")] - df[("meas", "g", "slot_Centroid_y")],
811 0.0,
812 0.0,
813 df[("meas", "g", "base_LocalWcs_CDMatrix_1_1")],
814 df[("meas", "g", "base_LocalWcs_CDMatrix_1_2")],
815 df[("meas", "g", "base_LocalWcs_CDMatrix_2_1")],
816 df[("meas", "g", "base_LocalWcs_CDMatrix_2_2")])
818 # Test functor values against afw SkyWcs computations.
819 for centX, centY, sep in zip(testPixelDeltas[:, 0],
820 testPixelDeltas[:, 1],
821 sepRadians.values):
822 afwSepRadians = skyOrigin.separation(
823 wcs.pixelToSky(x + centX, y + centY)).asRadians()
824 self.assertAlmostEqual(1 - sep / afwSepRadians, 0, places=5)
826 # Test the pixel scale computation.
827 func = ComputePixelScale("base_LocalWcs_CDMatrix_1_1",
828 "base_LocalWcs_CDMatrix_1_2",
829 "base_LocalWcs_CDMatrix_2_1",
830 "base_LocalWcs_CDMatrix_2_2")
831 pixelScale = self._funcVal(func, df)
832 self.assertTrue(np.allclose(
833 wcs.getPixelScale(center).asArcseconds(),
834 pixelScale.values,
835 rtol=1e-6,
836 atol=0))
838 # Test pixel -> arcsec conversion.
839 func = ConvertPixelToArcseconds("dipoleSep",
840 "base_LocalWcs_CDMatrix_1_1",
841 "base_LocalWcs_CDMatrix_1_2",
842 "base_LocalWcs_CDMatrix_2_1",
843 "base_LocalWcs_CDMatrix_2_2")
844 val = self._funcVal(func, df)
845 self.assertTrue(np.allclose(pixelScale.values * dipoleSep,
846 val.values,
847 atol=1e-16,
848 rtol=1e-16))
850 # Test pixel^2 -> arcsec^2 conversion.
851 func = ConvertPixelSqToArcsecondsSq("ixx",
852 "base_LocalWcs_CDMatrix_1_1",
853 "base_LocalWcs_CDMatrix_1_2",
854 "base_LocalWcs_CDMatrix_2_1",
855 "base_LocalWcs_CDMatrix_2_2")
856 val = self._funcVal(func, df)
857 self.assertTrue(np.allclose(pixelScale.values ** 2 * dipoleSep,
858 val.values,
859 atol=1e-16,
860 rtol=1e-16))
862 def _makeWcs(self, dec=53.1595451514076, theta=0):
863 """Create a wcs from real CFHT values, rotated by an optional theta.
865 dec : `float`
866 Set reference declination of CRVAL2 [degrees]
867 theta : `float`
868 Rotate CD matrix by theta [degrees]
870 Returns
871 -------
872 wcs : `lsst.afw.geom`
873 Created wcs.
874 """
875 metadata = dafBase.PropertySet()
877 metadata.set("SIMPLE", "T")
878 metadata.set("BITPIX", -32)
879 metadata.set("NAXIS", 2)
880 metadata.set("NAXIS1", 1024)
881 metadata.set("NAXIS2", 1153)
882 metadata.set("RADECSYS", 'FK5')
883 metadata.set("EQUINOX", 2000.)
885 metadata.setDouble("CRVAL1", 215.604025685476)
886 metadata.setDouble("CRVAL2", dec)
887 metadata.setDouble("CRPIX1", 1109.99981456774)
888 metadata.setDouble("CRPIX2", 560.018167811613)
889 metadata.set("CTYPE1", 'RA---SIN')
890 metadata.set("CTYPE2", 'DEC--SIN')
892 cd_matrix = np.array(
893 [
894 [5.10808596133527E-05, 1.85579539217196E-07],
895 [-8.27440751733828E-07, -5.10281493481982E-05]
896 ]
897 )
898 # rotate CD matrix
899 theta_rad = np.deg2rad(theta)
900 rotation_matrix = np.array(
901 [
902 [np.cos(theta_rad), -np.sin(theta_rad)],
903 [np.sin(theta_rad), np.cos(theta_rad)],
904 ]
905 )
906 cd_matrix = rotation_matrix @ cd_matrix
908 metadata.setDouble("CD1_1", cd_matrix[0, 0])
909 metadata.setDouble("CD1_2", cd_matrix[0, 1])
910 metadata.setDouble("CD2_1", cd_matrix[1, 0])
911 metadata.setDouble("CD2_2", cd_matrix[1, 1])
913 return afwGeom.makeSkyWcs(metadata)
915 def testHtm(self):
916 """Test that HtmIndxes are created as expected.
917 """
918 df = self.getMultiIndexDataFrame(self.dataDict)
919 func = HtmIndex20("coord_ra", "coord_dec")
921 val = self._funcVal(func, df)
922 # Test that the HtmIds come out as the ra/dec in dataDict.
923 self.assertTrue(np.all(np.equal(
924 val.values,
925 [14924528684992, 14924528689697, 14924528501716, 14924526434259,
926 14924526433879])))
928 def testEbv(self):
929 """Test that EBV works.
930 """
931 df = self.getMultiIndexDataFrame(self.dataDict)
932 func = Ebv()
934 val = self._funcVal(func, df)
935 np.testing.assert_array_almost_equal(
936 val.values,
937 [0.029100, 0.029013, 0.028857, 0.028802, 0.028797]
938 )
940 def testSkyMoments(self):
941 self.columns.extend([
942 "slot_Shape_xx",
943 "slot_Shape_yy",
944 "slot_Shape_xy",
945 "base_LocalWcs_CDMatrix_1_1",
946 "base_LocalWcs_CDMatrix_1_2",
947 "base_LocalWcs_CDMatrix_2_1",
948 "base_LocalWcs_CDMatrix_1_1",
949 ])
951 # CD Matrix from a ComCam exposure.
952 self.dataDict["base_LocalWcs_CDMatrix_1_1"] = \
953 np.full(self.nRecords, -9.695088e-07)
954 self.dataDict["base_LocalWcs_CDMatrix_1_2"] = \
955 np.full(self.nRecords, 3.950301849959342e-09)
956 self.dataDict["base_LocalWcs_CDMatrix_2_1"] = \
957 np.full(self.nRecords, 3.8766915166433014e-09)
958 self.dataDict["base_LocalWcs_CDMatrix_2_2"] = \
959 np.full(self.nRecords, 9.695092484727074e-07)
960 self.dataDict["slot_Shape_xx"] = \
961 np.array([6.52675084, 74.17426471, 6.45283335, 36.82870958, 6.45685472])
962 self.dataDict["slot_Shape_yy"] = \
963 np.array([6.12848637, 80.99510036, 6.05671667, 35.79219613, 5.97778765])
964 self.dataDict["slot_Shape_xy"] = \
965 np.array([-0.10281872, 0.88788384, -0.1261287, -1.60504171, 0.11974515])
966 self.dataDict["slot_Shape_sigma_x"] = np.sqrt(self.dataDict["slot_Shape_xx"])
967 self.dataDict["slot_Shape_sigma_y"] = np.sqrt(self.dataDict["slot_Shape_yy"])
968 self.dataDict["slot_Shape_rho"] = self.dataDict["slot_Shape_xy"]/(
969 self.dataDict["slot_Shape_sigma_x"]*self.dataDict["slot_Shape_sigma_y"]
970 )
972 args_cd = [
973 "base_LocalWcs_CDMatrix_1_1", "base_LocalWcs_CDMatrix_1_2",
974 "base_LocalWcs_CDMatrix_2_1", "base_LocalWcs_CDMatrix_2_2",
975 ]
976 args = ["slot_Shape_xx", "slot_Shape_yy", "slot_Shape_xy"] + args_cd
977 args_corr = ["slot_Shape_sigma_x", "slot_Shape_sigma_y", "slot_Shape_rho"] + args_cd
979 skyXX_functor = MomentsIuuSky(*args, filt="g")
980 skyYY_functor = MomentsIvvSky(*args, filt="g")
981 skyXY_functor = MomentsIuvSky(*args, filt="g")
983 axesA_functor = SemimajorAxisFromMoments(*args, filt="g")
984 axesB_functor = SemiminorAxisFromMoments(*args, filt="g")
985 axesTheta_functor = PositionAngleFromMoments(*args, filt="g")
987 skyXX_corr_functor = CorrelationIuuSky(*args_corr, filt="g")
988 skyYY_corr_functor = CorrelationIvvSky(*args_corr, filt="g")
989 skyXY_corr_functor = CorrelationIuvSky(*args_corr, filt="g")
991 axesA_corr_functor = SemimajorAxisFromCorrelation(*args_corr, filt="g")
992 axesB_corr_functor = SemiminorAxisFromCorrelation(*args_corr, filt="g")
993 axesTheta_corr_functor = PositionAngleFromCorrelation(*args_corr, filt="g")
995 df = self.getMultiIndexDataFrame(self.dataDict)
996 output_sky_xx = skyXX_functor(df)
997 output_sky_yy = skyYY_functor(df)
998 output_sky_xy = skyXY_functor(df)
1000 output_axes_a = axesA_functor(df)
1001 output_axes_b = axesB_functor(df)
1002 output_axes_theta = axesTheta_functor(df)
1004 output_sky_xx_corr = skyXX_corr_functor(df)
1005 output_sky_yy_corr = skyYY_corr_functor(df)
1006 output_sky_xy_corr = skyXY_corr_functor(df)
1008 output_axes_a_corr = axesA_corr_functor(df)
1009 output_axes_b_corr = axesB_corr_functor(df)
1010 output_axes_theta_corr = axesTheta_corr_functor(df)
1012 transformed_xx = []
1013 transformed_yy = []
1014 transformed_xy = []
1015 axes_a = []
1016 axes_b = []
1017 axes_theta = []
1019 for n in range(5):
1020 Ixx = df[('meas', 'g', 'slot_Shape_xx')].iloc[n]
1021 Iyy = df[('meas', 'g', 'slot_Shape_yy')].iloc[n]
1022 Ixy = df[('meas', 'g', 'slot_Shape_xy')].iloc[n]
1023 localWCS_CD_1_1 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_1_1')].iloc[n]
1024 localWCS_CD_2_1 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_2_1')].iloc[n]
1025 localWCS_CD_1_2 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_1_2')].iloc[n]
1026 localWCS_CD_2_2 = df[('meas', 'g', 'base_LocalWcs_CDMatrix_2_2')].iloc[n]
1027 CD_matrix = np.array([[localWCS_CD_1_1, localWCS_CD_1_2],
1028 [localWCS_CD_2_1, localWCS_CD_2_2]])
1030 q = afwGeom.ellipses.Quadrupole(Ixx, Iyy, Ixy)
1031 lt = geom.LinearTransform(CD_matrix)
1032 transformed_q = q.transform(lt)
1033 transformed_q.scale((180/np.pi) * (3600))
1035 axes = afwGeom.ellipses.Axes(transformed_q)
1037 transformed_xx.append(transformed_q.getIxx())
1038 transformed_yy.append(transformed_q.getIyy())
1039 transformed_xy.append(transformed_q.getIxy())
1041 axes_a.append(axes.getA())
1042 axes_b.append(axes.getB())
1043 axes_theta.append(np.degrees(axes.getTheta()))
1045 self.assertFloatsAlmostEqual(output_sky_xx, np.array(transformed_xx), rtol=1e-5)
1046 self.assertFloatsAlmostEqual(output_sky_yy, np.array(transformed_yy), rtol=1e-5)
1047 self.assertFloatsAlmostEqual(output_sky_xy, np.array(transformed_xy), rtol=1e-5)
1048 self.assertFloatsAlmostEqual(output_sky_xx_corr, np.array(transformed_xx), rtol=1e-5)
1049 self.assertFloatsAlmostEqual(output_sky_yy_corr, np.array(transformed_yy), rtol=1e-5)
1050 self.assertFloatsAlmostEqual(output_sky_xy_corr, np.array(transformed_xy), rtol=1e-5)
1052 self.assertFloatsAlmostEqual(output_axes_a, np.array(axes_a), rtol=1e-5)
1053 self.assertFloatsAlmostEqual(output_axes_b, np.array(axes_b), rtol=1e-5)
1054 self.assertFloatsAlmostEqual(output_axes_theta, np.array(axes_theta), rtol=1e-5)
1055 self.assertFloatsAlmostEqual(output_axes_a_corr, np.array(axes_a), rtol=1e-5)
1056 self.assertFloatsAlmostEqual(output_axes_b_corr, np.array(axes_b), rtol=1e-5)
1057 self.assertFloatsAlmostEqual(output_axes_theta_corr, np.array(axes_theta), rtol=1e-5)
1059 def _dropLevels(self, df):
1060 levelsToDrop = [n for lev, n in zip(df.columns.levels, df.columns.names) if len(lev) == 1]
1062 # Prevent error when trying to drop *all* columns
1063 if len(levelsToDrop) == len(df.columns.names):
1064 levelsToDrop.remove(df.columns.names[-1])
1066 df.columns = df.columns.droplevel(levelsToDrop)
1069class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase):
1070 pass
1073def setup_module(module):
1074 lsst.utils.tests.init()
1077if __name__ == "__main__": 1077 ↛ 1078line 1077 didn't jump to line 1078 because the condition on line 1077 was never true
1078 lsst.utils.tests.init()
1079 unittest.main()