LSST Applications g0aad566f14+98bdd3b178,g0dd44d6229+e0c2448475,g2079a07aa2+86d27d4dc4,g2305ad1205+28192cebb5,g2bbee38e9b+047b288a59,g337abbeb29+047b288a59,g33d1c0ed96+047b288a59,g3a166c0a6a+047b288a59,g3d1719c13e+98bdd3b178,g487adcacf7+08349b74a1,g4be5004598+98bdd3b178,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+742623b69a,g63cd9335cc+585e252eca,g858d7b2824+98bdd3b178,g88963caddf+0cb8e002cc,g88c60e1ba4+10d1129650,g99cad8db69+26979ce5a5,g9ddcbc5298+9a081db1e4,ga1e77700b3+a912195c07,gae0086650b+585e252eca,gb0e22166c9+60f28cb32d,gb1fda7adca+90d62c6f66,gb3a676b8dc+b4feba26a1,gb4b16eec92+950562e714,gba4ed39666+c2a2e4ac27,gbb8dafda3b+2fe9a6a885,gc120e1dc64+0553034b2f,gc28159a63d+047b288a59,gc3e9b769f7+f8a3f20ab6,gcf0d15dbbd+e0c2448475,gdaeeff99f8+f9a426f77a,gdb0af172c8+f61b8374e1,ge79ae78c31+047b288a59,w.2024.20
LSST Data Management Base Package
Loading...
Searching...
No Matches
testUtils.py
Go to the documentation of this file.
1# This file is part of jointcal.
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/>.
21
22"""Functions to help create jointcal tests by generating fake data."""
23
24__all__ = ['createFakeCatalog', 'createTwoFakeCcdImages', 'getMeasuredStarsFromCatalog']
25
26import os
27import unittest
28
29import numpy as np
30
31import lsst.afw.geom
32import lsst.afw.table
33import lsst.daf.butler
34import lsst.pipe.base
35
36import lsst.jointcal
37
38
40 """Returns True if the necessary packages and files are available.
41
42 We need ``obs_cfht`` to load the test/data/cfht_minimal dataset, which
43 includes the metadata that is used to build the fake catalogs.
44 """
45 try:
46 import lsst.obs.cfht # noqa: F401
47 return True
48 except ImportError:
49 return False
50
51
52def createTwoFakeCcdImages(num1=4, num2=4, seed=100, fakeDetectorId=12,
53 photoCalibMean1=1e-2, photoCalibMean2=1.2e-2,
54 fakeWcses=(None, None),
55 fakeVisitInfos=(None, None)):
56 """Return two fake ccdImages built on CFHT Megacam metadata.
57
58 If ``num1 == num2``, the catalogs will align on-sky so each source will
59 have a match in the other catalog.
60
61 This uses the butler dataset stored in `tests/data/cfht_minimal` to
62 bootstrap the metadata.
63
64 Parameters
65 ----------
66 num1, num2 : `int`, optional
67 Number of sources to put in the first and second catalogs. Should be
68 a square, to have sqrt(num) centroids on a grid.
69 seed : `int`, optional
70 Seed value for np.random.
71 fakeDetectorId : `int`, optional
72 Sensor identifier to use for both CcdImages. The wcs, bbox, photoCalib, etc.
73 will still be drawn from the CFHT ccd=12 files, as that is the only
74 testdata that is included in this simple test dataset.
75 photoCalibMean1, photoCalibMean2: `float`, optional
76 The mean photometric calibration to pass to each ccdImage construction.
77 Note: this value is 1/instFluxMag0, so it should be less than 1.
78 fakeWcses : `list` [`lsst.afw.geom.SkyWcs`], optional
79 The SkyWcses to use instead of the ones read from disk.
80 fakeWcses : `list` [`lsst.afw.image.VisitInfo`], optional
81 The VisitInfos to use instead of the ones read from disk.
82
83 Returns
84 -------
85 struct : `lsst.pipe.base.Struct`
86 Result struct with components:
87
88 - `camera` : Camera representing these catalogs
89 (`lsst.afw.cameraGeom.Camera`).
90 - `catalogs` : Catalogs containing fake sources
91 (`list` of `lsst.afw.table.SourceCatalog`).
92 - `ccdImageList` : CcdImages containing the metadata and fake sources
93 (`list` of `lsst.jointcal.CcdImage`).
94 - `bbox` : Bounding Box of the image (`lsst.geom.Box2I`).
95 - 'fluxFieldName' : name of the instFlux field in the catalogs ('str').
96 """
97 if not canRunTests():
98 msg = "Necessary packages not available to run tests that use the cfht_minimal dataset."
99 raise unittest.SkipTest(msg)
100
101 np.random.seed(seed)
102
103 visit1 = 849375
104 visit2 = 850587
105 fluxFieldName = "SomeFlux"
106
107 # Load or fake the necessary metadata for each CcdImage
108 dataDir = lsst.utils.getPackageDir('jointcal')
109 inputDir = os.path.join(dataDir, 'tests/data/cfht_minimal/repo')
110 # Ensure this butler is not writeable, so that we don't mess up the repo accidentally.
111 butler = lsst.daf.butler.Butler(inputDir, collections=["singleFrame"], writeable=False)
112
113 # so we can access parts of the camera later (e.g. focal plane)
114 camera = butler.get('camera', instrument="MegaPrime")
115
116 struct1 = createFakeCcdImage(butler, visit1, num1, fluxFieldName,
117 photoCalibMean=photoCalibMean1, photoCalibErr=1.0,
118 fakeDetectorId=fakeDetectorId,
119 fakeWcs=fakeWcses[0], fakeVisitInfo=fakeVisitInfos[0])
120 struct2 = createFakeCcdImage(butler, visit2, num2, fluxFieldName,
121 photoCalibMean=photoCalibMean2, photoCalibErr=5.0,
122 fakeDetectorId=fakeDetectorId,
123 fakeWcs=fakeWcses[1], fakeVisitInfo=fakeVisitInfos[1])
124
125 return lsst.pipe.base.Struct(camera=camera,
126 catalogs=[struct1.catalog, struct2.catalog],
127 ccdImageList=[struct1.ccdImage, struct2.ccdImage],
128 bbox=struct1.bbox,
129 skyWcs=[struct1.skyWcs, struct2.skyWcs],
130 fluxFieldName=fluxFieldName)
131
132
133def createFakeCcdImage(butler, visit, num, fluxFieldName,
134 photoCalibMean=1e-2, photoCalibErr=1.0, fakeDetectorId=12,
135 fakeWcs=None, fakeVisitInfo=None):
136 """Create a fake CcdImage by making a fake catalog.
137
138 Parameters
139 ----------
140 butler : `lsst.daf.butler.Butler`
141 Butler to load metadata from.
142 visit : `int`
143 Visit identifier to build a butler dataId.
144 num : `int`
145 Number of sources to put in the catalogs. Should be
146 a square, to have sqrt(num) centroids on a grid.
147 fluxFieldName : `str`
148 Name of the flux field to populate in the catalog, without `_instFlux`
149 (e.g. "slot_CalibFlux").
150 photoCalibMean : `float`, optional
151 Value to set for calibrationMean in the created PhotoCalib.
152 Note: this value is 1/instFluxMag0, so it should be less than 1.
153 photoCalibErr : `float`, optional
154 Value to set for calibrationErr in the created PhotoCalib.
155 fakeDetectorId : `int`, optional
156 Use this as the detectorId in the returned CcdImage.
157 fakeWcs : `lsst.afw.geom.SkyWcs`, optional
158 A SkyWcs to use instead of one read from disk.
159 fakeVisitInfo : `lsst.afw.image.VisitInfo`, optional
160 A VisitInfo to use instead of one read from disk.
161
162 Returns
163 -------
164 struct : `lsst.pipe.base.Struct`
165 Result struct with components:
166
167 - `catalog` : Catalogs containing fake sources
168 (`lsst.afw.table.SourceCatalog`).
169 - `ccdImage` : CcdImage containing the metadata and fake sources
170 (`lsst.jointcal.CcdImage`).
171 - `bbox` : Bounding Box of the image (`lsst.geom.Box2I`).
172 - `skyWcs` : SkyWcs of the image (`lsst.afw.geom.SkyWcs`).
173 """
174 detectorId = 12 # we only have data for detector=12
175
176 dataId = dict(visit=visit, detector=detectorId, instrument="MegaPrime")
177 skyWcs = fakeWcs if fakeWcs is not None else butler.get('calexp.wcs', dataId=dataId)
178 visitInfo = fakeVisitInfo if fakeVisitInfo is not None else butler.get('calexp.visitInfo', dataId=dataId)
179 bbox = butler.get('calexp.bbox', dataId=dataId)
180 detector = butler.get('calexp.detector', dataId=dataId)
181 filt = butler.get("calexp.filter", dataId=dataId).bandLabel
182 photoCalib = lsst.afw.image.PhotoCalib(photoCalibMean, photoCalibErr)
183
184 catalog = createFakeCatalog(num, bbox, fluxFieldName, skyWcs=skyWcs)
185 ccdImage = lsst.jointcal.CcdImage(catalog, skyWcs, visitInfo, bbox, filt, photoCalib,
186 detector, visit, fakeDetectorId, fluxFieldName)
187
188 return lsst.pipe.base.Struct(catalog=catalog, ccdImage=ccdImage, bbox=bbox, skyWcs=skyWcs)
189
190
191def createFakeCatalog(num, bbox, fluxFieldName, skyWcs=None, refCat=False):
192 """Return a fake minimally-useful catalog for jointcal.
193
194 Parameters
195 ----------
196 num : `int`
197 Number of sources to put in the catalogs. Should be
198 a square, to have sqrt(num) centroids on a grid.
199 bbox : `lsst.geom.Box2I`
200 Bounding Box of the detector to populate.
201 fluxFieldName : `str`
202 Name of the flux field to populate in the catalog, without `_instFlux`
203 (e.g. "slot_CalibFlux").
204 skyWcs : `lsst.afw.geom.SkyWcs` or None, optional
205 If supplied, use this to fill in coordinates from centroids.
206 refCat : `bool`, optional
207 Return a ``SimpleCatalog`` so that it behaves like a reference catalog?
208
209 Returns
210 -------
211 catalog : `lsst.afw.table.SourceCatalog`
212 A populated source catalog.
213 """
215 # centroid
216 centroidKey = lsst.afw.table.Point2DKey.addFields(schema, "centroid", "centroid", "pixels")
217 xErrKey = schema.addField("centroid_xErr", type="F")
218 yErrKey = schema.addField("centroid_yErr", type="F")
219 # coordinate erros
221 # shape
222 shapeKey = lsst.afw.table.QuadrupoleKey.addFields(schema, "shape", "",
223 lsst.afw.table.CoordinateType.PIXEL)
224 # Put the fake sources in the minimal catalog.
225 schema.addField(fluxFieldName+"_instFlux", type="D", doc="post-ISR instFlux")
226 schema.addField(fluxFieldName+"_instFluxErr", type="D", doc="post-ISR instFlux stddev")
227 schema.addField(fluxFieldName+"_flux", type="D", doc="source flux (nJy)")
228 schema.addField(fluxFieldName+"_fluxErr", type="D", doc="flux stddev (nJy)")
229 schema.addField(fluxFieldName+"_mag", type="D", doc="magnitude")
230 schema.addField(fluxFieldName+"_magErr", type="D", doc="magnitude stddev")
231 return fillCatalog(schema, num, bbox,
232 centroidKey, xErrKey, yErrKey, shapeKey, fluxFieldName,
233 skyWcs=skyWcs, refCat=refCat)
234
235
236def fillCatalog(schema, num, bbox,
237 centroidKey, xErrKey, yErrKey, shapeKey, fluxFieldName,
238 skyWcs=None, fluxErrFraction=0.05, refCat=False):
239 """Return a catalog populated with fake, but reasonable, sources.
240
241 Centroids are placed on a uniform grid, errors are normally distributed.
242
243 Parameters
244 ----------
245 schema : `lsst.afw.table.Schema`
246 Pre-built schema to make the catalog from.
247 num : `int`
248 Number of sources to put in the catalog.
249 bbox : `lsst.geom.Box2I`
250 Bounding box of the ccd to put sources in.
251 centroidKey : `lsst.afw.table.Key`
252 Key for the centroid field to populate.
253 xErrKey : `lsst.afw.table.Key`
254 Key for the xErr field to populate.
255 yErrKey : `lsst.afw.table.Key`
256 Key for the yErr field to populate.
257 shapeKey : `lsst.afw.table.Key`
258 Key for the shape field to populate.
259 fluxFieldName : `str`
260 Name of the flux field to populate in the catalog, without `_instFlux`
261 (e.g. "slot_CalibFlux").
262 skyWcs : `lsst.afw.geom.SkyWcs` or None, optional
263 If supplied, use this to fill in coordinates from centroids.
264 fluxErrFraction : `float`, optional
265 Fraction of instFlux to use for the instFluxErr.
266 refCat : `bool`, optional
267 Return a ``SimpleCatalog`` so that it behaves like a reference catalog?
268
269 Returns
270 -------
271 catalog : `lsst.afw.table.SourceCatalog`
272 The filled catalog.
273 """
274 table = lsst.afw.table.SourceTable.make(schema)
275 table.defineCentroid('centroid')
276 table.defineShape('shape')
277 table.defineCalibFlux(fluxFieldName)
278 if refCat:
279 catalog = lsst.afw.table.SimpleCatalog(table)
280 else:
281 catalog = lsst.afw.table.SourceCatalog(table)
282
283 instFlux = np.random.random(num)*10000
284 instFluxErr = np.abs(instFlux * np.random.normal(fluxErrFraction, scale=0.1, size=num))
285 xx = np.linspace(bbox.getMinX(), bbox.getMaxX(), int(np.sqrt(num)))
286 yy = np.linspace(bbox.getMinY(), bbox.getMaxY(), int(np.sqrt(num)))
287 xv, yv = np.meshgrid(xx, yy)
288 vx = np.random.normal(scale=0.1, size=num)
289 vy = np.random.normal(scale=0.1, size=num)
290
291 # make all the sources perfectly spherical, for simplicity.
292 mxx = 1
293 myy = 1
294 mxy = 0
295
296 for i, (x, y) in enumerate(zip(xv.ravel(), yv.ravel())):
297 record = catalog.addNew()
298 record.set('id', i)
299 record.set(centroidKey, lsst.geom.Point2D(x, y))
300 record.set(shapeKey, lsst.afw.geom.ellipses.Quadrupole(mxx, myy, mxy))
301
302 if skyWcs is not None:
303 lsst.afw.table.updateSourceCoords(skyWcs, catalog)
304
305 catalog[xErrKey] = vx
306 catalog[yErrKey] = vy
307 catalog[fluxFieldName + '_instFlux'] = instFlux
308 catalog[fluxFieldName + '_instFluxErr'] = instFluxErr
309
310 return catalog
311
312
313def getMeasuredStarsFromCatalog(catalog, pixToFocal):
314 """Return a list of measuredStars built from a catalog.
315
316 Parameters
317 ----------
318 catalog : `lsst.afw.table.SourceCatalog`
319 The table to get sources from.
320 pixToFocal : `lsst.afw.geom.TransformPoint2ToPoint2`
321 Transform that goes from pixel to focal plane coordinates, to set the
322 MeasuredStar x/y focal points.
323
324 Returns
325 -------
326 stars : `list` of `lsst.jointcal.MeasuredStar`
327 MeasuredStars built from the catalog sources.
328 """
329 stars = []
330 for record in catalog:
332 star.x = record.getX()
333 star.y = record.getY()
334 star.setInstFluxAndErr(record.getCalibInstFlux(), record.getCalibInstFluxErr())
335 # TODO: cleanup after DM-4044
336 point = lsst.geom.Point2D(star.x, star.y)
337 pointFocal = pixToFocal.applyForward(point)
338 star.setXFocal(pointFocal.getX())
339 star.setYFocal(pointFocal.getY())
340 stars.append(star)
341
342 return stars
An ellipse core with quadrupole moments as parameters.
Definition Quadrupole.h:47
The photometric calibration of an exposure.
Definition PhotoCalib.h:114
static ErrorKey addErrorFields(Schema &schema)
static QuadrupoleKey addFields(Schema &schema, std::string const &name, std::string const &doc, CoordinateType coordType=CoordinateType::PIXEL)
Add a set of quadrupole subfields to a schema and return a QuadrupoleKey that points to them.
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
static std::shared_ptr< SourceTable > make(Schema const &schema, std::shared_ptr< IdFactory > const &idFactory)
Construct a new table.
Definition Source.cc:400
static Schema makeMinimalSchema()
Return a minimal schema for Source tables and records.
Definition Source.h:258
Handler of an actual image from a single CCD.
Definition CcdImage.h:64
Sources measured on images.
void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList, bool include_covariance=true)
Update sky coordinates in a collection of source objects.
Definition wcsUtils.cc:125
createTwoFakeCcdImages(num1=4, num2=4, seed=100, fakeDetectorId=12, photoCalibMean1=1e-2, photoCalibMean2=1.2e-2, fakeWcses=(None, None), fakeVisitInfos=(None, None))
Definition testUtils.py:55
fillCatalog(schema, num, bbox, centroidKey, xErrKey, yErrKey, shapeKey, fluxFieldName, skyWcs=None, fluxErrFraction=0.05, refCat=False)
Definition testUtils.py:238
createFakeCcdImage(butler, visit, num, fluxFieldName, photoCalibMean=1e-2, photoCalibErr=1.0, fakeDetectorId=12, fakeWcs=None, fakeVisitInfo=None)
Definition testUtils.py:135
getMeasuredStarsFromCatalog(catalog, pixToFocal)
Definition testUtils.py:313
createFakeCatalog(num, bbox, fluxFieldName, skyWcs=None, refCat=False)
Definition testUtils.py:191