Coverage for tests/test_utils.py: 20%
234 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-08 02:12 -0800
1# This file is part of summit_utils.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://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 <https://www.gnu.org/licenses/>.
22"""Test cases for utils."""
24import copy
25import datetime
26import itertools
27import unittest
28from typing import Iterable
30import astropy.time
31import astropy.units as u
32import numpy as np
33from astro_metadata_translator import makeObservationInfo
34from astropy.coordinates import HeliocentricEclipticIAU76, SkyCoord
35from astropy.time import Time
37import lsst.afw.detection as afwDetect
38import lsst.afw.geom as afwGeom
39import lsst.afw.image as afwImage
40import lsst.geom as geom
41import lsst.utils.tests
42from lsst.obs.base import createInitialSkyWcsFromBoresight
43from lsst.obs.base.makeRawVisitInfoViaObsInfo import MakeRawVisitInfoViaObsInfo
44from lsst.obs.lsst import Latiss
45from lsst.obs.lsst.translators.latiss import AUXTEL_LOCATION
46from lsst.summit.utils.utils import getFilterSeeingCorrection # deprecated
47from lsst.summit.utils.utils import (
48 calcEclipticCoords,
49 computeCcdExposureId,
50 computeExposureId,
51 fluxesFromFootprints,
52 getAirmassSeeingCorrection,
53 getBandpassSeeingCorrection,
54 getCurrentDayObs_datetime,
55 getCurrentDayObs_humanStr,
56 getCurrentDayObs_int,
57 getExpPositionOffset,
58 getFieldNameAndTileNumber,
59 getQuantiles,
60 getSunAngle,
61 quickSmooth,
62)
65class ExpSkyPositionOffsetTestCase(lsst.utils.tests.TestCase):
66 """A test case for testing sky position offsets for exposures."""
68 def setUp(self):
69 camera = Latiss.getCamera()
70 self.assertTrue(len(camera) == 1)
71 self.detector = camera[0]
73 self.viMaker = MakeRawVisitInfoViaObsInfo()
74 self.mi = afwImage.MaskedImageF(0, 0)
75 self.baseHeader = dict(
76 boresight_airmass=1.5,
77 temperature=15 * u.deg_C,
78 observation_type="science",
79 exposure_time=5 * u.ks,
80 detector_num=32,
81 location=AUXTEL_LOCATION,
82 )
84 def test_getExpPositionOffset(self):
85 epsilon = 0.0001
86 ra1s = [0, 45, 90]
87 ra2s = copy.copy(ra1s)
88 ra2s.extend([r + epsilon for r in ra1s])
89 ra1s = np.deg2rad(ra1s)
90 ra2s = np.deg2rad(ra2s)
92 epsilon = 0.0001
93 dec1s = [0, 45, 90]
94 dec2s = copy.copy(dec1s)
95 dec2s.extend([d + epsilon for d in dec1s[:-1]]) # skip last point as >90 not allowed for dec
97 rotAngle1 = geom.Angle(43.2, geom.degrees) # arbitrary non-zero
98 rotAngle2 = geom.Angle(56.7, geom.degrees)
100 t1 = astropy.time.Time("2021-09-15T12:00:00", format="isot", scale="utc")
101 t2 = astropy.time.Time("2021-09-15T12:01:00", format="isot", scale="utc")
102 expTime = astropy.time.TimeDelta(20, format="sec")
104 header1 = copy.copy(self.baseHeader)
105 header2 = copy.copy(self.baseHeader)
106 header1["datetime_begin"] = astropy.time.Time(t1, format="isot", scale="utc")
107 header2["datetime_begin"] = astropy.time.Time(t2, format="isot", scale="utc")
109 header1["datetime_end"] = astropy.time.Time(t1 + expTime, format="isot", scale="utc")
110 header2["datetime_end"] = astropy.time.Time(t2 + expTime, format="isot", scale="utc")
112 obsInfo1 = makeObservationInfo(**header1)
113 obsInfo2 = makeObservationInfo(**header2)
115 vi1 = self.viMaker.observationInfo2visitInfo(obsInfo1)
116 vi2 = self.viMaker.observationInfo2visitInfo(obsInfo2)
117 expInfo1 = afwImage.ExposureInfo()
118 expInfo1.setVisitInfo(vi1)
119 expInfo2 = afwImage.ExposureInfo()
120 expInfo2.setVisitInfo(vi2)
122 for ra1, dec1, ra2, dec2 in itertools.product(ra1s, dec1s, ra2s, dec2s):
123 pos1 = geom.SpherePoint(ra1, dec1, geom.degrees)
124 pos2 = geom.SpherePoint(ra2, dec2, geom.degrees)
126 wcs1 = createInitialSkyWcsFromBoresight(pos1, rotAngle1, self.detector, flipX=True)
127 wcs2 = createInitialSkyWcsFromBoresight(pos2, rotAngle2, self.detector, flipX=True)
129 exp1 = afwImage.ExposureF(self.mi, expInfo1)
130 exp2 = afwImage.ExposureF(self.mi, expInfo2)
132 exp1.setWcs(wcs1)
133 exp2.setWcs(wcs2)
135 result = getExpPositionOffset(exp1, exp2)
137 deltaRa = ra1 - ra2
138 deltaDec = dec1 - dec2
140 self.assertAlmostEqual(result.deltaRa.asDegrees(), deltaRa, 6)
141 self.assertAlmostEqual(result.deltaDec.asDegrees(), deltaDec, 6)
144class MiscUtilsTestCase(lsst.utils.tests.TestCase):
145 def setUp(self) -> None:
146 return super().setUp()
148 def test_getFieldNameAndTileNumber(self):
149 field, num = getFieldNameAndTileNumber("simple")
150 self.assertEqual(field, "simple")
151 self.assertIsNone(num)
153 field, num = getFieldNameAndTileNumber("_simple")
154 self.assertEqual(field, "_simple")
155 self.assertIsNone(num)
157 field, num = getFieldNameAndTileNumber("simple_321")
158 self.assertEqual(field, "simple")
159 self.assertEqual(num, 321)
161 field, num = getFieldNameAndTileNumber("_simple_321")
162 self.assertEqual(field, "_simple")
163 self.assertEqual(num, 321)
165 field, num = getFieldNameAndTileNumber("test_321a_123")
166 self.assertEqual(field, "test_321a")
167 self.assertEqual(num, 123)
169 field, num = getFieldNameAndTileNumber("test_321a_123_")
170 self.assertEqual(field, "test_321a_123_")
171 self.assertIsNone(num)
173 field, num = getFieldNameAndTileNumber("test_321a_123a")
174 self.assertEqual(field, "test_321a_123a")
175 self.assertIsNone(num)
177 field, num = getFieldNameAndTileNumber("test_321a:asd_asd-dsa_321")
178 self.assertEqual(field, "test_321a:asd_asd-dsa")
179 self.assertEqual(num, 321)
181 def test_getAirmassSeeingCorrection(self):
182 for airmass in (1.1, 2.0, 20.0):
183 correction = getAirmassSeeingCorrection(airmass)
184 self.assertGreater(correction, 0.01)
185 self.assertLess(correction, 1.0)
187 correction = getAirmassSeeingCorrection(1)
188 self.assertEqual(correction, 1.0)
190 with self.assertRaises(ValueError):
191 getAirmassSeeingCorrection(0.5)
193 def test_getFilterSeeingCorrection(self):
194 for filterName in ("SDSSg_65mm", "SDSSr_65mm", "SDSSi_65mm"):
195 with self.assertWarns(FutureWarning):
196 correction = getFilterSeeingCorrection(filterName)
197 self.assertGreater(correction, 0.5)
198 self.assertLess(correction, 1.5)
200 def test_getBandpassSeeingCorrection(self):
201 for filterName in (
202 "SDSSg_65mm",
203 "SDSSr_65mm",
204 "SDSSi_65mm",
205 "u_02",
206 "g_01",
207 "r_03",
208 "i_06",
209 "z_03",
210 "y_04",
211 ):
212 correction = getBandpassSeeingCorrection(filterName)
213 self.assertGreater(correction, 0.5)
214 self.assertLess(correction, 1.5)
216 def test_quickSmooth(self):
217 # just test that it runs and returns the right shape. It's a wrapper on
218 # scipy.ndimage.gaussian_filter we can trust that it does what it
219 # should, and we just test the interface hasn't bitrotted on either end
220 data = np.zeros((100, 100), dtype=np.float32)
221 data = quickSmooth(data, 5.0)
222 self.assertEqual(data.shape, (100, 100))
224 def test_getCurrentDayObs_datetime(self):
225 """Just a type check and a basic sanity check on the range.
227 Setting days=3 as the tolerance just because of timezones and who knows
228 what really.
229 """
230 dt = getCurrentDayObs_datetime()
231 self.assertIsInstance(dt, datetime.date)
232 self.assertLess(dt, datetime.date.today() + datetime.timedelta(days=3))
233 self.assertGreater(dt, datetime.date.today() - datetime.timedelta(days=3))
235 def test_getCurrentDayObs_int(self):
236 """Just a type check and a basic sanity check on the range."""
237 dayObs = getCurrentDayObs_int()
238 self.assertIsInstance(dayObs, int)
239 self.assertLess(dayObs, 21000101)
240 self.assertGreater(dayObs, 19700101)
242 def test_getCurrentDayObs_humanStr(self):
243 """Just a basic formatting check."""
244 dateStr = getCurrentDayObs_humanStr()
245 self.assertIsInstance(dateStr, str)
246 self.assertEqual(len(dateStr), 10)
247 self.assertRegex(dateStr, r"\d{4}-\d{2}-\d{2}")
249 def test_getSunAngle(self):
250 """Just a basic sanity check on the range."""
251 testTime = Time("2021-09-15T16:00:00", format="isot", scale="utc")
252 sunAngle = getSunAngle(testTime)
253 self.assertIsInstance(sunAngle, float)
254 self.assertLess(sunAngle, 90.01)
255 self.assertGreater(sunAngle, -90.01)
257 sunAngle = getSunAngle()
258 self.assertIsInstance(sunAngle, float)
259 self.assertLess(sunAngle, 90.01)
260 self.assertGreater(sunAngle, -90.01)
263class QuantileTestCase(lsst.utils.tests.TestCase):
264 def setUp(self) -> None:
265 return super().setUp()
267 def test_quantiles(self):
268 # We understand that our algorithm gives very large rounding error
269 # compared to the generic numpy method. But still test it.
270 np.random.seed(1234)
271 dataRanges = [(50, 1, -1), (100_000, 5_000, -2), (5_000_000, 10_000, -2), (50_000, 100_000, -3)]
272 colorRanges = [2, 256, 999] # [very few, nominal, lots and an odd number]
273 for nColors, (mean, width, decimal) in itertools.product(colorRanges, dataRanges):
274 data = np.random.normal(mean, width, (100, 100))
275 data[10, 10] = np.nan # check we're still nan-safe
276 if np.nanmax(data) - np.nanmin(data) > 300_000:
277 with self.assertLogs(level="WARNING") as cm:
278 edges1 = getQuantiles(data, nColors)
279 self.assertIn("Data range is very large", cm.output[0])
280 else:
281 with self.assertNoLogs(level="WARNING") as cm:
282 edges1 = getQuantiles(data, nColors)
283 edges2 = np.nanquantile(data, np.linspace(0, 1, nColors + 1)) # must check with nanquantile
284 np.testing.assert_almost_equal(edges1, edges2, decimal=decimal)
287class ImageBasedTestCase(lsst.utils.tests.TestCase):
288 def test_fluxFromFootprint(self):
289 image = afwImage.Image(
290 np.arange(8100, dtype=np.int32).reshape(90, 90), xy0=lsst.geom.Point2I(10, 12), dtype="I"
291 )
293 radius = 3
294 spans = afwGeom.SpanSet.fromShape(radius, afwGeom.Stencil.CIRCLE, offset=(27, 30))
295 footprint1 = afwDetect.Footprint(spans)
297 # The extracted footprint should be the same as the product of the
298 # spans and the overlapped bow with the image
299 truth1 = spans.asArray() * image.array[15:22, 14:21]
301 radius = 3
302 spans = afwGeom.SpanSet.fromShape(radius, afwGeom.Stencil.CIRCLE, offset=(44, 49))
303 footprint2 = afwDetect.Footprint(spans)
304 truth2 = spans.asArray() * image.array[34:41, 31:38]
306 allFootprints = [footprint1, footprint2]
307 footprintSet = afwDetect.FootprintSet(image.getBBox())
308 footprintSet.setFootprints(allFootprints)
310 # check it can accept a footprintSet, and single and iterables of
311 # footprints
312 with self.assertRaises(TypeError):
313 fluxesFromFootprints(10, image)
315 with self.assertRaises(TypeError):
316 fluxesFromFootprints([8, 6, 7, 5, 3, 0, 9], image)
318 # check the footPrintSet
319 fluxes = fluxesFromFootprints(footprintSet, image)
320 expectedLength = len(footprintSet.getFootprints())
321 self.assertEqual(len(fluxes), expectedLength) # always one flux per footprint
322 self.assertIsInstance(fluxes, Iterable)
323 self.assertAlmostEqual(fluxes[0], np.sum(truth1))
324 self.assertAlmostEqual(fluxes[1], np.sum(truth2))
326 # check the list of footprints
327 fluxes = fluxesFromFootprints(allFootprints, image)
328 expectedLength = 2
329 self.assertEqual(len(fluxes), expectedLength) # always one flux per footprint
330 self.assertIsInstance(fluxes, Iterable)
331 self.assertAlmostEqual(fluxes[0], np.sum(truth1))
332 self.assertAlmostEqual(fluxes[1], np.sum(truth2))
334 # ensure that subtracting the image median from fluxes leave image
335 # pixels untouched
336 oldImageArray = copy.deepcopy(image.array)
337 fluxes = fluxesFromFootprints(footprintSet, image, subtractImageMedian=True)
338 np.testing.assert_array_equal(image.array, oldImageArray)
341class IdTestCase(lsst.utils.tests.TestCase):
342 def test_exposure_id(self):
343 self.assertEqual(computeExposureId("latiss", "O", 20240402, 35), 2024040200035)
344 self.assertEqual(computeExposureId("LATISS", "C", 20240402, 35), 2024040200035)
345 self.assertEqual(computeExposureId("LSSTComCamSim", "S", 20240402, 35), 7024040200035)
346 with self.assertRaises(ValueError):
347 computeExposureId("bad_instrument", "O", 20240402, 35)
349 def test_ccdexposure_id(self):
350 self.assertEqual(computeCcdExposureId("latiss", 2024040200035, 0), 5205 * (2**23) + 35 * 256 + 0)
351 with self.assertRaises(ValueError):
352 computeCcdExposureId("latiss", 20240402000035, 1)
353 with self.assertRaises(ValueError):
354 computeCcdExposureId("LsstComCam", 20240402000035, 9)
355 with self.assertRaises(ValueError):
356 computeCcdExposureId("LsstCam", 20240402000035, 205)
357 self.assertEqual(
358 computeCcdExposureId("LSSTComCamSim", 7024040200035, 2),
359 5 * (2**37) + 5205 * (2**23) + 35 * 256 + 2,
360 )
362 def test_calcEclipticCoords(self):
363 ras = [0, 1, 45, 90, 180, 270, 359.9]
364 decs = [-90, -80, -45, -1, 0, 1, 45, 90]
366 for ra, dec in itertools.product(ras, decs):
367 p = SkyCoord(ra=ra * u.deg, dec=dec * u.deg, distance=1 * u.au, frame="hcrs")
368 eclCoords = p.transform_to(HeliocentricEclipticIAU76)
369 _lambda, beta = calcEclipticCoords(ra, dec)
370 self.assertAlmostEqual(eclCoords.lat.deg, beta, 6)
371 self.assertAlmostEqual(eclCoords.lon.deg, _lambda, 6)
374class TestMemory(lsst.utils.tests.MemoryTestCase):
375 pass
378def setup_module(module):
379 lsst.utils.tests.init()
382if __name__ == "__main__": 382 ↛ 383line 382 didn't jump to line 383 because the condition on line 382 was never true
383 lsst.utils.tests.init()
384 unittest.main()