Coverage for python / lsst / meas / algorithms / testUtils.py: 21%
96 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-07 08:25 +0000
1#
2# LSST Data Management System
3#
4# Copyright 2008-2017 AURA/LSST.
5#
6# This product includes software developed by the
7# LSST Project (http://www.lsst.org/).
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 LSST License Statement and
20# the GNU General Public License along with this program. If not,
21# see <https://www.lsstcorp.org/LegalNotices/>.
22#
24__all__ = ["plantSources", "makeRandomTransmissionCurve", "makeDefectList",
25 "MockReferenceObjectLoaderFromFiles", "MockRefcatDataId",
26 "MockReferenceObjectLoaderFromMemory"]
28import numpy as np
29import esutil
31import lsst.geom
32import lsst.afw.image as afwImage
33from lsst.pipe.base import InMemoryDatasetHandle
34from lsst import sphgeom
35from . import SingleGaussianPsf
36from . import Defect
38from . import ReferenceObjectLoader
39import lsst.afw.table as afwTable
42def plantSources(bbox, kwid, sky, coordList, addPoissonNoise=True):
43 """Make an exposure with stars (modelled as Gaussians)
45 Parameters
46 ----------
47 bbox : `lsst.geom.Box2I`
48 Parent bbox of exposure
49 kwid : `int`
50 Kernal width (and height; kernal is square)
51 sky : `float`
52 Amount of sky background (counts)
53 coordList : `list [tuple]`
54 A list of [x, y, counts, sigma] where:
55 * x,y are relative to exposure origin
56 * counts is the integrated counts for the star
57 * sigma is the Gaussian sigma in pixels
58 addPoissonNoise : `bool`
59 If True: add Poisson noise to the exposure
60 """
61 # make an image with sources
62 img = afwImage.ImageD(bbox)
63 meanSigma = 0.0
64 for coord in coordList:
65 x, y, counts, sigma = coord
66 meanSigma += sigma
68 # make a single gaussian psf
69 psf = SingleGaussianPsf(kwid, kwid, sigma)
71 # make an image of it and scale to the desired number of counts
72 thisPsfImg = psf.computeImage(lsst.geom.PointD(x, y))
73 thisPsfImg *= counts
75 # bbox a window in our image and add the fake star image
76 psfBox = thisPsfImg.getBBox()
77 psfBox.clip(bbox)
78 if psfBox != thisPsfImg.getBBox():
79 thisPsfImg = thisPsfImg[psfBox, afwImage.PARENT]
80 imgSeg = img[psfBox, afwImage.PARENT]
81 imgSeg += thisPsfImg
82 meanSigma /= len(coordList)
84 img += sky
86 # add Poisson noise
87 if (addPoissonNoise):
88 # Make results predictable over different numpy versions.
89 rng = np.random.Generator(np.random.MT19937(5))
90 imgArr = img.getArray()
91 imgArr[:] = rng.poisson(imgArr)
93 # bundle into a maskedimage and an exposure
94 mask = afwImage.Mask(bbox)
95 var = img.convertFloat()
96 img -= sky
97 mimg = afwImage.MaskedImageF(img.convertFloat(), mask, var)
98 exposure = afwImage.makeExposure(mimg)
100 # insert an approximate psf
101 psf = SingleGaussianPsf(kwid, kwid, meanSigma)
102 exposure.setPsf(psf)
104 return exposure
107def makeRandomTransmissionCurve(rng, minWavelength=4000.0, maxWavelength=7000.0, nWavelengths=200,
108 maxRadius=80.0, nRadii=30, perturb=0.05):
109 """Create a random TransmissionCurve with nontrivial spatial and
110 wavelength variation.
112 Parameters
113 ----------
114 rng : numpy.random.RandomState
115 Random number generator.
116 minWavelength : float
117 Average minimum wavelength for generated TransmissionCurves (will be
118 randomly perturbed).
119 maxWavelength : float
120 Average maximum wavelength for generated TransmissionCurves (will be
121 randomly perturbed).
122 nWavelengths : int
123 Number of samples in the wavelength dimension.
124 maxRadius : float
125 Average maximum radius for spatial variation (will be perturbed).
126 nRadii : int
127 Number of samples in the radial dimension.
128 perturb: float
129 Fraction by which wavelength and radius bounds should be randomly
130 perturbed.
131 """
132 dWavelength = maxWavelength - minWavelength
134 def perturbed(x, s=perturb*dWavelength):
135 return x + 2.0*s*(rng.rand() - 0.5)
137 wavelengths = np.linspace(perturbed(minWavelength), perturbed(maxWavelength), nWavelengths)
138 radii = np.linspace(0.0, perturbed(maxRadius, perturb*maxRadius), nRadii)
139 throughput = np.zeros(wavelengths.shape + radii.shape, dtype=float)
140 # throughput will be a rectangle in wavelength, shifting to higher wavelengths and shrinking
141 # in height with radius, going to zero at all bounds.
142 peak0 = perturbed(0.9, 0.05)
143 start0 = perturbed(minWavelength + 0.25*dWavelength)
144 stop0 = perturbed(minWavelength + 0.75*dWavelength)
145 for i, r in enumerate(radii):
146 mask = np.logical_and(wavelengths >= start0 + r, wavelengths <= stop0 + r)
147 throughput[mask, i] = peak0*(1.0 - r/1000.0)
148 return afwImage.TransmissionCurve.makeRadial(throughput, wavelengths, radii)
151def makeDefectList():
152 """Create a list of defects that can be used for testing.
154 Returns
155 -------
156 defectList = `list` [`lsst.meas.algorithms.Defect`]
157 The list of defects.
158 """
159 defectList = [Defect(lsst.geom.Box2I(lsst.geom.Point2I(962, 0),
160 lsst.geom.Extent2I(2, 4611))),
161 Defect(lsst.geom.Box2I(lsst.geom.Point2I(1316, 0),
162 lsst.geom.Extent2I(2, 4611))),
163 Defect(lsst.geom.Box2I(lsst.geom.Point2I(1576, 0),
164 lsst.geom.Extent2I(4, 4611))),
165 Defect(lsst.geom.Box2I(lsst.geom.Point2I(1626, 0),
166 lsst.geom.Extent2I(2, 4611))),
167 Defect(lsst.geom.Box2I(lsst.geom.Point2I(1994, 252),
168 lsst.geom.Extent2I(2, 4359))),
169 Defect(lsst.geom.Box2I(lsst.geom.Point2I(1426, 702),
170 lsst.geom.Extent2I(2, 3909))),
171 Defect(lsst.geom.Box2I(lsst.geom.Point2I(1526, 1140),
172 lsst.geom.Extent2I(2, 3471))),
173 Defect(lsst.geom.Box2I(lsst.geom.Point2I(856, 2300),
174 lsst.geom.Extent2I(2, 2311))),
175 Defect(lsst.geom.Box2I(lsst.geom.Point2I(858, 2328),
176 lsst.geom.Extent2I(2, 65))),
177 Defect(lsst.geom.Box2I(lsst.geom.Point2I(859, 2328),
178 lsst.geom.Extent2I(1, 56))),
179 Defect(lsst.geom.Box2I(lsst.geom.Point2I(844, 2796),
180 lsst.geom.Extent2I(4, 1814))),
181 Defect(lsst.geom.Box2I(lsst.geom.Point2I(1366, 2804),
182 lsst.geom.Extent2I(2, 1806))),
183 Defect(lsst.geom.Box2I(lsst.geom.Point2I(1766, 3844),
184 lsst.geom.Extent2I(2, 766))),
185 Defect(lsst.geom.Box2I(lsst.geom.Point2I(1872, 4228),
186 lsst.geom.Extent2I(2, 382))),
187 ]
189 return defectList
192class MockRefcatDataId:
193 """Mock reference catalog dataId.
195 The reference catalog dataId is only used to retrieve a region property.
197 Parameters
198 ----------
199 region : `lsst.sphgeom.Region`
200 The region associated with this mock dataId.
201 """
202 def __init__(self, region):
203 self._region = region
205 @property
206 def region(self):
207 return self._region
210class MockReferenceObjectLoaderFromFiles(ReferenceObjectLoader):
211 """A mock of ReferenceObjectLoader using files on disk.
213 This mock ReferenceObjectLoader uses a set of files on disk to create
214 mock dataIds and data reference handles that can be accessed
215 without a butler. The files must be afw catalog files in the reference
216 catalog format, sharded with HTM pixelization.
218 Parameters
219 ----------
220 filenames : `list` [`str`]
221 Names of files to use.
222 config : `lsst.meas.astrom.LoadReferenceObjectsConfig`, optional
223 Configuration object if necessary to override defaults.
224 htmLevel : `int`, optional
225 HTM level to use for the loader.
226 """
227 def __init__(self, filenames, name='cal_ref_cat', config=None, htmLevel=4):
228 dataIds, refCats = self._createDataIdsAndRefcats(filenames, htmLevel, name)
230 super().__init__(dataIds, refCats, name=name, config=config)
232 def _createDataIdsAndRefcats(self, filenames, htmLevel, name):
233 """Create mock dataIds and refcat handles.
235 Parameters
236 ----------
237 filenames : `list` [`str`]
238 Names of files to use.
239 htmLevel : `int`
240 HTM level to use for the loader.
241 name : `str`
242 Name of reference catalog (for logging).
244 Returns
245 -------
246 dataIds : `list` [`MockRefcatDataId`]
247 List of mock dataIds.
248 refCats : `list` [`lsst.pipe.base.InMemoryDatasetHandle`]
249 List of mock deferred dataset handles.
251 Raises
252 ------
253 RuntimeError if any file contains sources that cover more than one HTM
254 pixel at level ``htmLevel``.
255 """
256 pixelization = sphgeom.HtmPixelization(htmLevel)
257 htm = esutil.htm.HTM(htmLevel)
259 dataIds = []
260 refCats = []
262 for filename in filenames:
263 cat = afwTable.BaseCatalog.readFits(filename)
265 ids = htm.lookup_id(np.rad2deg(cat['coord_ra']), np.rad2deg(cat['coord_dec']))
267 if len(np.unique(ids)) != 1:
268 raise RuntimeError(f"File {filename} contains more than one pixel at level {htmLevel}")
270 dataIds.append(MockRefcatDataId(pixelization.pixel(ids[0])))
271 refCats.append(InMemoryDatasetHandle(cat, name=name))
273 return dataIds, refCats
276class MockReferenceObjectLoaderFromMemory(ReferenceObjectLoader):
277 """A mock of ReferenceObjectLoader using catalogs in memory.
279 Parameters
280 ----------
281 catalogs : `list` [`lsst.afw.table.SimpleCatalog`]
282 In-memory catalogs to use to mock dataIds.
283 config : `lsst.meas.astrom.LoadReferenceObjectsConfig`, optional
284 Configuration object if necessary to override defaults.
285 htmLevel : `int`, optional
286 HTM level to use for the loader.
287 """
288 def __init__(self, catalogs, name='mock_ref_cat', config=None, htmLevel=4):
289 dataIds, refCats = self._createDataIdsAndRefcats(catalogs, htmLevel, name)
290 super().__init__(dataIds, refCats, name=name, config=config)
292 def _createDataIdsAndRefcats(self, catalogs, htmLevel, name):
293 pixelization = sphgeom.HtmPixelization(htmLevel)
294 htm = esutil.htm.HTM(htmLevel)
296 dataIds = []
297 refCats = []
299 for i, catalog in enumerate(catalogs):
300 ids = htm.lookup_id(np.rad2deg(catalog['coord_ra']), np.rad2deg(catalog['coord_dec']))
302 if len(np.unique(ids)) != 1:
303 raise RuntimeError(f"Catalog number {i} contains more than one pixel at level {htmLevel}")
305 dataIds.append(MockRefcatDataId(pixelization.pixel(ids[0])))
306 refCats.append(InMemoryDatasetHandle(catalog, name=name))
308 return dataIds, refCats