LSST Applications g013ef56533+7c9321ec0f,g042eb84c57+c6cfa41bc3,g199a45376c+0ba108daf9,g1fd858c14a+fcad0d0313,g210f2d0738+c0f94c6586,g262e1987ae+a7e710680e,g29ae962dfc+fb55f2edb0,g2ac17093b6+61d6563b1e,g2b1d02342f+df6f932764,g2cef7863aa+aef1011c0b,g2f7ad74990+c0f94c6586,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53cf87ae69,g47891489e3+4316d04fff,g511e8cfd20+baa56acf6c,g53246c7159+8c5ae1fdc5,g54cd7ddccb+fd7ad03fde,g64539dfbff+c0f94c6586,g67b6fd64d1+4316d04fff,g67fd3c3899+c0f94c6586,g6985122a63+4316d04fff,g74acd417e5+ca833bee28,g786e29fd12+668abc6043,g81db2e9a8d+b2ec8e584f,g87389fa792+8856018cbb,g89139ef638+4316d04fff,g8d7436a09f+0a24083b20,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+11eb8fd5b8,gbf99507273+8c5ae1fdc5,gcdda8b9158+e4c84c9d5c,gce8aa8abaa+8c5ae1fdc5,gd7ef33dd92+4316d04fff,gdab6d2f7ff+ca833bee28,ge410e46f29+4316d04fff,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.40
LSST Data Management Base Package
Loading...
Searching...
No Matches
fgcmBuildFromIsolatedStars.py
Go to the documentation of this file.
1# See COPYRIGHT file at the top of the source tree.
2#
3# This file is part of fgcmcal.
4#
5# Developed for the LSST Data Management System.
6# This product includes software developed by the LSST Project
7# (https://www.lsst.org).
8# See the COPYRIGHT file at the top-level directory of this distribution
9# for details of code ownership.
10#
11# This program is free software: you can redistribute it and/or modify
12# it under the terms of the GNU General Public License as published by
13# the Free Software Foundation, either version 3 of the License, or
14# (at your option) any later version.
15#
16# This program is distributed in the hope that it will be useful,
17# but WITHOUT ANY WARRANTY; without even the implied warranty of
18# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19# GNU General Public License for more details.
20#
21# You should have received a copy of the GNU General Public License
22# along with this program. If not, see <https://www.gnu.org/licenses/>.
23"""Build star observations for input to FGCM using sourceTable_visit.
24
25This task finds all the visits and sourceTable_visits in a repository (or a
26subset based on command line parameters) and extracts all the potential
27calibration stars for input into fgcm. This task additionally uses fgcm to
28match star observations into unique stars, and performs as much cleaning of the
29input catalog as possible.
30"""
31import warnings
32import numpy as np
33import esutil
34import hpgeom as hpg
35from smatch.matcher import Matcher
36from astropy.table import Table, vstack
37
38from fgcm.fgcmUtilities import objFlagDict, histogram_rev_sorted
39
40import lsst.afw.cameraGeom as afwCameraGeom
41import lsst.pex.config as pexConfig
42import lsst.pipe.base as pipeBase
43from lsst.pipe.base import connectionTypes
44from lsst.meas.algorithms import ReferenceObjectLoader, LoadReferenceObjectsConfig
45from lsst.pipe.tasks.reserveIsolatedStars import ReserveIsolatedStarsTask
46
47from .fgcmBuildStarsBase import FgcmBuildStarsConfigBase, FgcmBuildStarsBaseTask
48from .utilities import computeApproxPixelAreaFields, computeApertureRadiusFromName
49
50__all__ = ["FgcmBuildFromIsolatedStarsConfig", "FgcmBuildFromIsolatedStarsTask"]
51
52
53class FgcmBuildFromIsolatedStarsConnections(pipeBase.PipelineTaskConnections,
54 dimensions=("instrument",),
55 defaultTemplates={}):
56 camera = connectionTypes.PrerequisiteInput(
57 doc="Camera instrument",
58 name="camera",
59 storageClass="Camera",
60 dimensions=("instrument",),
61 isCalibration=True,
62 )
63 fgcm_lookup_table = connectionTypes.PrerequisiteInput(
64 doc=("Atmosphere + instrument look-up-table for FGCM throughput and "
65 "chromatic corrections."),
66 name="fgcmLookUpTable",
67 storageClass="Catalog",
68 dimensions=("instrument",),
69 deferLoad=True,
70 )
71 ref_cat = connectionTypes.PrerequisiteInput(
72 doc="Reference catalog to use for photometric calibration.",
73 name="cal_ref_cat",
74 storageClass="SimpleCatalog",
75 dimensions=("skypix",),
76 deferLoad=True,
77 multiple=True,
78 )
79 isolated_star_cats = pipeBase.connectionTypes.Input(
80 doc=("Catalog of isolated stars with average positions, number of associated "
81 "sources, and indexes to the isolated_star_sources catalogs."),
82 name="isolated_star_presource_associations",
83 storageClass="ArrowAstropy",
84 dimensions=("instrument", "tract", "skymap"),
85 deferLoad=True,
86 multiple=True,
87 )
88 isolated_star_sources = pipeBase.connectionTypes.Input(
89 doc=("Catalog of isolated star sources with sourceIds, and indexes to the "
90 "isolated_star_cats catalogs."),
91 name="isolated_star_presources",
92 storageClass="ArrowAstropy",
93 dimensions=("instrument", "tract", "skymap"),
94 deferLoad=True,
95 multiple=True,
96 )
97 visit_summaries = connectionTypes.Input(
98 doc=("Per-visit consolidated exposure metadata. These catalogs use "
99 "detector id for the id and must be sorted for fast lookups of a "
100 "detector."),
101 name="visitSummary",
102 storageClass="ExposureCatalog",
103 dimensions=("instrument", "visit"),
104 deferLoad=True,
105 multiple=True,
106 )
107 fgcm_visit_catalog = connectionTypes.Output(
108 doc="Catalog of visit information for fgcm.",
109 name="fgcmVisitCatalog",
110 storageClass="Catalog",
111 dimensions=("instrument",),
112 )
113 fgcm_star_observations = connectionTypes.Output(
114 doc="Catalog of star observations for fgcm.",
115 name="fgcm_star_observations",
116 storageClass="ArrowAstropy",
117 dimensions=("instrument",),
118 )
119 fgcm_star_ids = connectionTypes.Output(
120 doc="Catalog of fgcm calibration star IDs.",
121 name="fgcm_star_ids",
122 storageClass="ArrowAstropy",
123 dimensions=("instrument",),
124 )
125 fgcm_reference_stars = connectionTypes.Output(
126 doc="Catalog of fgcm-matched reference stars.",
127 name="fgcm_reference_stars",
128 storageClass="ArrowAstropy",
129 dimensions=("instrument",),
130 )
131
132 def __init__(self, *, config=None):
133 super().__init__(config=config)
134
135 if not config.doReferenceMatches:
136 self.prerequisiteInputs.remove("ref_cat")
137 self.outputs.remove("fgcm_reference_stars")
138
140 return ("isolated_star_cats", "visit_summaries")
141
142
144 pipelineConnections=FgcmBuildFromIsolatedStarsConnections):
145 """Config for FgcmBuildFromIsolatedStarsTask."""
146 referenceCCD = pexConfig.Field(
147 doc="Reference detector for checking PSF and background.",
148 dtype=int,
149 default=40,
150 )
151 reserve_selection = pexConfig.ConfigurableField(
152 target=ReserveIsolatedStarsTask,
153 doc="Task to select reserved stars.",
154 )
155
156 def setDefaults(self):
157 super().setDefaults()
158
159 self.reserve_selection.reserve_name = "fgcmcal"
160 self.reserve_selection.reserve_fraction = 0.1
161
162 # The names here correspond to the isolated_star_sources.
163 self.instFluxField = 'normCompTophatFlux_instFlux'
164 self.localBackgroundFluxField = 'localBackground_instFlux'
165 self.apertureInnerInstFluxField = 'apFlux_12_0_instFlux'
166 self.apertureOuterInstFluxField = 'apFlux_17_0_instFlux'
167
168 source_selector = self.sourceSelector["science"]
169 source_selector.setDefaults()
170
171 source_selector.doFlags = False
172 source_selector.doSignalToNoise = True
173 source_selector.doUnresolved = False
174 source_selector.doIsolated = False
175 source_selector.doRequireFiniteRaDec = False
176
177 source_selector.flags.bad = []
178
179 source_selector.signalToNoise.minimum = 9.0
180 source_selector.signalToNoise.maximum = 1000.0
181 source_selector.signalToNoise.fluxField = self.instFluxField
182 source_selector.signalToNoise.errField = self.instFluxField + 'Err'
183
184
186 """Build star catalog for FGCM global calibration, using the isolated star catalogs.
187 """
188 ConfigClass = FgcmBuildFromIsolatedStarsConfig
189 _DefaultName = "fgcmBuildFromIsolatedStars"
190
191 canMultiprocess = False
192
193 def __init__(self, initInputs=None, **kwargs):
194 super().__init__(**kwargs)
195 self.makeSubtask('reserve_selection')
196
197 def runQuantum(self, butlerQC, inputRefs, outputRefs):
198 input_ref_dict = butlerQC.get(inputRefs)
199
200 isolated_star_cat_handles = input_ref_dict["isolated_star_cats"]
201 isolated_star_source_handles = input_ref_dict["isolated_star_sources"]
202
203 isolated_star_cat_handle_dict = {
204 handle.dataId["tract"]: handle for handle in isolated_star_cat_handles
205 }
206 isolated_star_source_handle_dict = {
207 handle.dataId["tract"]: handle for handle in isolated_star_source_handles
208 }
209
210 if len(isolated_star_cat_handle_dict) != len(isolated_star_source_handle_dict):
211 raise RuntimeError("isolated_star_cats and isolate_star_sources must have same length.")
212
213 for tract in isolated_star_cat_handle_dict:
214 if tract not in isolated_star_source_handle_dict:
215 raise RuntimeError(f"tract {tract} in isolated_star_cats but not isolated_star_sources")
216
217 if self.config.doReferenceMatches:
218
219 # Prepare the reference catalog loader
220 ref_config = LoadReferenceObjectsConfig()
221 ref_config.filterMap = self.config.fgcmLoadReferenceCatalog.filterMap
222 ref_obj_loader = ReferenceObjectLoader(dataIds=[ref.datasetRef.dataId
223 for ref in inputRefs.ref_cat],
224 refCats=butlerQC.get(inputRefs.ref_cat),
225 name=self.config.connections.ref_cat,
226 log=self.log,
227 config=ref_config)
228 self.makeSubtask('fgcmLoadReferenceCatalog',
229 refObjLoader=ref_obj_loader,
230 refCatName=self.config.connections.ref_cat)
231
232 # The visit summary handles for use with fgcmMakeVisitCatalog must be keyed with
233 # visit, and values are a list with the first value being the visit_summary_handle,
234 # the second is the source catalog (which is not used, but a place-holder is
235 # left for compatibility reasons.)
236 visit_summary_handle_dict = {handle.dataId['visit']: [handle, None] for
237 handle in input_ref_dict['visit_summaries']}
238
239 camera = input_ref_dict["camera"]
240
241 lookup_table_handle = input_ref_dict["fgcm_lookup_table"]
242
243 struct = self.run(
244 camera=camera,
245 visit_summary_handle_dict=visit_summary_handle_dict,
246 isolated_star_cat_handle_dict=isolated_star_cat_handle_dict,
247 isolated_star_source_handle_dict=isolated_star_source_handle_dict,
248 lookup_table_handle=lookup_table_handle,
249 )
250
251 butlerQC.put(struct.fgcm_visit_catalog, outputRefs.fgcm_visit_catalog)
252 butlerQC.put(struct.fgcm_star_observations, outputRefs.fgcm_star_observations)
253 butlerQC.put(struct.fgcm_star_ids, outputRefs.fgcm_star_ids)
254 if self.config.doReferenceMatches:
255 butlerQC.put(struct.fgcm_reference_stars, outputRefs.fgcm_reference_stars)
256
257 def run(self, *, camera, visit_summary_handle_dict, isolated_star_cat_handle_dict,
258 isolated_star_source_handle_dict, lookup_table_handle):
259 """Run the fgcmBuildFromIsolatedStarsTask.
260
261 Parameters
262 ----------
263 camera : `lsst.afw.cameraGeom.Camera`
264 Camera object.
265 visit_summary_handle_dict : `dict` [`int`, [`lsst.daf.butler.DeferredDatasetHandle`]]
266 Visit summary dataset handles, with the visit as key.
267 isolated_star_cat_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
268 Isolated star catalog dataset handles, with the tract as key.
269 isolated_star_source_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
270 Isolated star source dataset handles, with the tract as key.
271 lookup_table_handle : `lsst.daf.butler.DeferredDatasetHandle`
272 Data reference to fgcm look-up table.
273
274 Returns
275 -------
276 struct : `lsst.pipe.base.struct`
277 Catalogs for persistence, with attributes:
278
279 ``fgcm_visit_catalog``
280 Catalog of per-visit data (`lsst.afw.table.ExposureCatalog`).
281 ``fgcm_star_observations``
282 Catalog of individual star observations (`astropy.table.Table`).
283 ``fgcm_star_ids``
284 Catalog of per-star ids and positions (`astropy.table.Table`).
285 ``fgcm_reference_stars``
286 Catalog of reference stars matched to fgcm stars (`astropy.table.Table`).
287 """
288 # Compute aperture radius if necessary. This is useful to do now before
289 # any heave lifting has happened (fail early).
290 calib_flux_aperture_radius = None
291 if self.config.doSubtractLocalBackground:
292 try:
293 calib_flux_aperture_radius = computeApertureRadiusFromName(self.config.instFluxField)
294 except RuntimeError as e:
295 raise RuntimeError("Could not determine aperture radius from %s. "
296 "Cannot use doSubtractLocalBackground." %
297 (self.config.instFluxField)) from e
298
299 # Check that we have the lookup_table_handle if we are doing reference matches.
300 if self.config.doReferenceMatches:
301 if lookup_table_handle is None:
302 raise RuntimeError("Must supply lookup_table_handle if config.doReferenceMatches is True.")
303
304 lut_cat = lookup_table_handle.get()
305 if len(camera) == lut_cat[0]["nCcd"]:
306 use_science_detectors = False
307 else:
308 # If the LUT has a different number of detectors than
309 # the camera, then we only want to use science detectors
310 # in the focal plane projector.
311 # Note that in the LUT building code we either have
312 # all detectors or only science detectors.
313 use_science_detectors = True
314 del lut_cat
315
316 visit_cat = self.fgcmMakeVisitCatalog(
317 camera,
318 visit_summary_handle_dict,
319 useScienceDetectors=use_science_detectors,
320 )
321
322 # Select and concatenate the isolated stars and sources.
323 fgcm_obj, star_obs = self._make_all_star_obs_from_isolated_stars(
324 isolated_star_cat_handle_dict,
325 isolated_star_source_handle_dict,
326 visit_cat,
327 camera,
328 calib_flux_aperture_radius=calib_flux_aperture_radius,
329 use_science_detectors=use_science_detectors,
330 )
331
333 visit_cat,
334 star_obs,
335 )
336
337 # Do density down-sampling.
338 self._density_downsample(fgcm_obj, star_obs)
339
340 # Mark reserve stars
341 self._mark_reserve_stars(fgcm_obj)
342
343 # Do reference association.
344 if self.config.doReferenceMatches:
345 fgcm_ref = self._associate_reference_stars(lookup_table_handle, visit_cat, fgcm_obj)
346 else:
347 fgcm_ref = None
348
349 return pipeBase.Struct(
350 fgcm_visit_catalog=visit_cat,
351 fgcm_star_observations=star_obs,
352 fgcm_star_ids=fgcm_obj,
353 fgcm_reference_stars=fgcm_ref,
354 )
355
357 self,
358 isolated_star_cat_handle_dict,
359 isolated_star_source_handle_dict,
360 visit_cat,
361 camera,
362 calib_flux_aperture_radius=None,
363 use_science_detectors=False,
364 ):
365 """Make all star observations from isolated star catalogs.
366
367 Parameters
368 ----------
369 isolated_star_cat_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
370 Isolated star catalog dataset handles, with the tract as key.
371 isolated_star_source_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
372 Isolated star source dataset handles, with the tract as key.
373 visit_cat : `lsst.afw.table.ExposureCatalog`
374 Catalog of per-visit data.
375 camera : `lsst.afw.cameraGeom.Camera`
376 The camera object.
377 calib_flux_aperture_radius : `float`, optional
378 Radius for the calibration flux aperture.
379 use_science_detectors : `bool`, optional
380 Only use science detectors?
381
382 Returns
383 -------
384 fgcm_obj : `astropy.table.Table`
385 Catalog of ids and positions for each star.
386 star_obs : `astropy.table.Table`
387 Catalog of individual star observations.
388 """
389 source_columns = [
390 "sourceId",
391 "visit",
392 "detector",
393 "ra",
394 "dec",
395 "x",
396 "y",
397 "physical_filter",
398 "band",
399 "obj_index",
400 self.config.instFluxField,
401 self.config.instFluxField + "Err",
402 self.config.apertureInnerInstFluxField,
403 self.config.apertureInnerInstFluxField + "Err",
404 self.config.apertureOuterInstFluxField,
405 self.config.apertureOuterInstFluxField + "Err",
406 ]
407
408 if self.config.doSubtractLocalBackground:
409 source_columns.append(self.config.localBackgroundFluxField)
410 local_background_flag_name = self.config.localBackgroundFluxField[0: -len('instFlux')] + 'flag'
411 source_columns.append(local_background_flag_name)
412
413 if self.sourceSelector.config.doFlags:
414 source_columns.extend(self.sourceSelector.config.flags.bad)
415
416 if self.config.doSubtractLocalBackground:
417 local_background_area = np.pi*calib_flux_aperture_radius**2.
418
419 # Compute the approximate pixel area over the full focal plane
420 # from the WCS jacobian using the camera model.
421 approx_pixel_area_fields = computeApproxPixelAreaFields(camera)
422
423 # Construct mapping from detector number to index.
424 detector_mapping = {}
425 for detector_index, detector in enumerate(camera):
426 detector_mapping[detector.getId()] = detector_index
427
428 star_obs_dtype = [
429 ("ra", "f8"),
430 ("dec", "f8"),
431 ("x", "f8"),
432 ("y", "f8"),
433 ("visit", "i8"),
434 ("detector", "i4"),
435 ("inst_mag", "f4"),
436 ("inst_mag_err", "f4"),
437 ("jacobian", "f4"),
438 ("delta_mag_bkg", "f4"),
439 ("delta_mag_aper", "f4"),
440 ("delta_mag_err_aper", "f4"),
441 ]
442
443 fgcm_obj_dtype = [
444 ("fgcm_id", "i8"),
445 ("isolated_star_id", "i8"),
446 ("ra", "f8"),
447 ("dec", "f8"),
448 ("obs_arr_index", "i4"),
449 ("n_obs", "i4"),
450 ("obj_flag", "i4"),
451 ]
452
453 fgcm_objs = []
454 star_obs_cats = []
455 merge_source_counter = 0
456
457 k = 2.5/np.log(10.)
458
459 visit_cat_table = visit_cat.asAstropy()
460
461 for tract in sorted(isolated_star_cat_handle_dict):
462 stars = isolated_star_cat_handle_dict[tract].get()
463 sources = isolated_star_source_handle_dict[tract].get(parameters={"columns": source_columns})
464
465 # Down-select sources.
466 good_sources = self.sourceSelector.selectSources(sources).selected
467 if self.config.doSubtractLocalBackground:
468 good_sources &= (~sources[local_background_flag_name])
469 local_background = local_background_area*sources[self.config.localBackgroundFluxField]
470 good_sources &= ((sources[self.config.instFluxField] - local_background) > 0)
471
472 # We need to do this check as early as possible to
473 # ensure we can do the match test below.
474 if good_sources.sum() == 0:
475 self.log.info("No good sources found in tract %d", tract)
476 continue
477
478 # We also need to make sure that all sources are in the
479 # input list of visits and detectors.
480
481 detector_ids = []
482 for detector in camera:
483 if use_science_detectors:
484 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE:
485 continue
486 detector_ids.append(detector.getId())
487
488 missing_sources1 = np.ones(len(sources), dtype=np.bool_)
489 _, sources_match = esutil.numpy_util.match(visit_cat_table["visit"], sources["visit"])
490 missing_sources1[sources_match] = False
491 good_sources &= (~missing_sources1)
492
493 missing_sources2 = np.ones(len(sources), dtype=np.bool_)
494 _, sources_match = esutil.numpy_util.match(detector_ids, sources["detector"])
495 missing_sources2[sources_match] = False
496 good_sources &= (~missing_sources2)
497
498 # Check again after further tests.
499 if good_sources.sum() == 0:
500 self.log.info("No good sources found in tract %d", tract)
501 continue
502
503 self.log.info("Tract %d contains %d good sources.", tract, good_sources.sum())
504
505 # Need to count the observations of each star after cuts, per band.
506 # If we have requiredBands specified, we must make sure that each star
507 # has the minumum number of observations in each of thos bands.
508 # Otherwise, we must make sure that each star has at least the minimum
509 # number of observations in _any_ band.
510 if len(self.config.requiredBands) > 0:
511 loop_bands = self.config.requiredBands
512 else:
513 loop_bands = np.unique(sources["band"])
514
515 n_req = np.zeros((len(loop_bands), len(stars)), dtype=np.int32)
516 for i, band in enumerate(loop_bands):
517 (band_use,) = (sources[good_sources]["band"] == band).nonzero()
518 np.add.at(
519 n_req,
520 (i, sources[good_sources]["obj_index"][band_use]),
521 1,
522 )
523
524 if len(self.config.requiredBands) > 0:
525 # The min gives us the band with the fewest observations, which must be
526 # above the limit.
527 (good_stars,) = (n_req.min(axis=0) >= self.config.minPerBand).nonzero()
528 else:
529 # The max gives us the band with the most observations, which must be
530 # above the limit.
531 (good_stars,) = (n_req.max(axis=0) >= self.config.minPerBand).nonzero()
532
533 if len(good_stars) == 0:
534 self.log.info("No good stars found in tract %d", tract)
535 continue
536
537 # With the following matching:
538 # sources[good_sources][b] <-> stars[good_stars[a]]
539 obj_index = sources["obj_index"][good_sources]
540 a, b = esutil.numpy_util.match(good_stars, obj_index)
541
542 # Update indexes and cut to selected stars/sources
543 _, index_new = np.unique(a, return_index=True)
544 stars["source_cat_index"][good_stars] = index_new
545 sources = sources[good_sources][b]
546 sources["obj_index"][:] = a
547 stars = stars[good_stars]
548
549 nsource = np.zeros(len(stars), dtype=np.int32)
550 np.add.at(
551 nsource,
552 sources["obj_index"],
553 1,
554 )
555 stars["nsource"][:] = nsource
556
557 # After these cuts, the catalogs have the following properties:
558 # - ``stars`` only contains isolated stars that have the minimum number of good
559 # sources in the required bands.
560 # - ``sources`` has been cut to the good sources.
561 # - The slice [stars["source_cat_index"]: stars["source_cat_index"]
562 # + stars["nsource"]]
563 # applied to ``sources`` will give all the sources associated with the star.
564 # - For each source, sources["obj_index"] points to the index of the associated
565 # isolated star.
566
567 # We now reformat the sources and compute the ``observations`` that fgcm expects.
568 star_obs = Table(data=np.zeros(len(sources), dtype=star_obs_dtype))
569 star_obs["ra"] = sources["ra"]
570 star_obs["dec"] = sources["dec"]
571 star_obs["x"] = sources["x"]
572 star_obs["y"] = sources["y"]
573 star_obs["visit"] = sources["visit"]
574 star_obs["detector"] = sources["detector"]
575
576 visit_match, obs_match = esutil.numpy_util.match(visit_cat_table["visit"], sources["visit"])
577
578 exp_time = np.zeros(len(star_obs))
579 exp_time[obs_match] = visit_cat_table["exptime"][visit_match]
580
581 with warnings.catch_warnings():
582 # Ignore warnings, we will filter infinities and nans below.
583 warnings.simplefilter("ignore")
584
585 inst_mag_inner = -2.5*np.log10(sources[self.config.apertureInnerInstFluxField])
586 inst_mag_err_inner = k*(sources[self.config.apertureInnerInstFluxField + "Err"]
587 / sources[self.config.apertureInnerInstFluxField])
588 inst_mag_outer = -2.5*np.log10(sources[self.config.apertureOuterInstFluxField])
589 inst_mag_err_outer = k*(sources[self.config.apertureOuterInstFluxField + "Err"]
590 / sources[self.config.apertureOuterInstFluxField])
591 star_obs["delta_mag_aper"] = inst_mag_inner - inst_mag_outer
592 star_obs["delta_mag_err_aper"] = np.sqrt(inst_mag_err_inner**2. + inst_mag_err_outer**2.)
593 # Set bad values to sentinel values for fgcm.
594 bad = ~np.isfinite(star_obs["delta_mag_aper"])
595 star_obs["delta_mag_aper"][bad] = 99.0
596 star_obs["delta_mag_err_aper"][bad] = 99.0
597
598 if self.config.doSubtractLocalBackground:
599 # At the moment we only adjust the flux and not the flux
600 # error by the background because the error on
601 # base_LocalBackground_instFlux is the rms error in the
602 # background annulus, not the error on the mean in the
603 # background estimate (which is much smaller, by sqrt(n)
604 # pixels used to estimate the background, which we do not
605 # have access to in this task). In the default settings,
606 # the annulus is sufficiently large such that these
607 # additional errors are negligibly small (much less
608 # than a mmag in quadrature).
609
610 # This is the difference between the mag with local background correction
611 # and the mag without local background correction.
612 local_background = local_background_area*sources[self.config.localBackgroundFluxField]
613 star_obs["delta_mag_bkg"] = (-2.5*np.log10(sources[self.config.instFluxField]
614 - local_background) -
615 -2.5*np.log10(sources[self.config.instFluxField]))
616
617 # Need to loop over detectors here.
618 for detector in camera:
619 if use_science_detectors:
620 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE:
621 continue
622
623 detector_id = detector.getId()
624 # used index for all observations with a given detector
625 (use,) = (star_obs["detector"][obs_match] == detector_id).nonzero()
626 # Prior to running the calibration, we want to remove the effect
627 # of the jacobian of the WCS because that is a known quantity.
628 # Ideally, this would be done for each individual WCS, but this
629 # is extremely slow and makes small differences that are much
630 # smaller than the variation in the throughput due to other issues.
631 # Therefore, we use the approximate jacobian estimated from the
632 # camera model.
633 jac = approx_pixel_area_fields[detector_id].evaluate(
634 star_obs["x"][obs_match][use],
635 star_obs["y"][obs_match][use],
636 )
637 star_obs["jacobian"][obs_match[use]] = jac
638 scaled_inst_flux = (sources[self.config.instFluxField][obs_match[use]]
639 * visit_cat_table["scaling"][visit_match[use],
640 detector_mapping[detector_id]])
641 star_obs["inst_mag"][obs_match[use]] = (-2.5 * np.log10(scaled_inst_flux
642 / exp_time[use]))
643
644 # Compute instMagErr from inst_flux_err/inst_flux.
645 # This does not need to be computed per-detector because there
646 # is no need for the per-detector scaling which cancels out.
647 star_obs["inst_mag_err"] = k*(sources[self.config.instFluxField + "Err"]
648 / sources[self.config.instFluxField])
649
650 # Apply the jacobian if configured to do so.
651 if self.config.doApplyWcsJacobian:
652 star_obs["inst_mag"] -= 2.5*np.log10(star_obs["jacobian"])
653
654 # We now reformat the stars and compute the ''objects'' that fgcm expects.
655 fgcm_obj = Table(data=np.zeros(len(stars), dtype=fgcm_obj_dtype))
656 fgcm_obj["isolated_star_id"] = stars["isolated_star_id"]
657 fgcm_obj["ra"] = stars["ra"]
658 fgcm_obj["dec"] = stars["dec"]
659 fgcm_obj["obs_arr_index"] = stars["source_cat_index"]
660 fgcm_obj["n_obs"] = stars["nsource"]
661
662 # Offset indexes to account for tract merging
663 fgcm_obj["obs_arr_index"] += merge_source_counter
664
665 fgcm_objs.append(fgcm_obj)
666 star_obs_cats.append(star_obs)
667
668 merge_source_counter += len(star_obs)
669
670 fgcm_obj = vstack(fgcm_objs)
671
672 # Set the fgcm_id to a unique 64-bit integer for easier sorting.
673 fgcm_obj["fgcm_id"][:] = np.arange(len(fgcm_obj)) + 1
674
675 return fgcm_obj, vstack(star_obs_cats)
676
677 def _density_downsample(self, fgcm_obj, star_obs):
678 """Downsample stars according to density.
679
680 Catalogs are modified in-place.
681
682 Parameters
683 ----------
684 fgcm_obj : `astropy.table.Table`
685 Catalog of per-star ids and positions.
686 star_obs : `astropy.table.Table`
687 Catalog of individual star observations.
688 """
689 if self.config.randomSeed is not None:
690 np.random.seed(seed=self.config.randomSeed)
691
692 ipnest = hpg.angle_to_pixel(self.config.densityCutNside, fgcm_obj["ra"], fgcm_obj["dec"])
693 # Use the esutil.stat.histogram function to pull out the histogram of
694 # grouped pixels, and the rev_indices which describes which inputs
695 # are grouped together. The fgcm histogram_rev_sorted shim
696 # ensures that the indices are properly sorted.
697 hist, rev_indices = histogram_rev_sorted(ipnest)
698
699 obj_use = np.ones(len(fgcm_obj), dtype=bool)
700
701 (high,) = (hist > self.config.densityCutMaxPerPixel).nonzero()
702 (ok,) = (hist > 0).nonzero()
703 self.log.info("There are %d/%d pixels with high stellar density.", high.size, ok.size)
704 for i in range(high.size):
705 # The pix_indices are the indices of every star in the pixel.
706 pix_indices = rev_indices[rev_indices[high[i]]: rev_indices[high[i] + 1]]
707 # Cut down to the maximum number of stars in the pixel.
708 cut = np.random.choice(
709 pix_indices,
710 size=pix_indices.size - self.config.densityCutMaxPerPixel,
711 replace=False,
712 )
713 obj_use[cut] = False
714
715 fgcm_obj = fgcm_obj[obj_use]
716
717 obs_index = np.zeros(np.sum(fgcm_obj["n_obs"]), dtype=np.int32)
718 ctr = 0
719 for i in range(len(fgcm_obj)):
720 n_obs = fgcm_obj["n_obs"][i]
721 obs_index[ctr: ctr + n_obs] = (
722 np.arange(fgcm_obj["obs_arr_index"][i], fgcm_obj["obs_arr_index"][i] + n_obs)
723 )
724 fgcm_obj["obs_arr_index"][i] = ctr
725 ctr += n_obs
726
727 star_obs = star_obs[obs_index]
728
729 def _mark_reserve_stars(self, fgcm_obj):
730 """Run the star reservation task to flag reserved stars.
731
732 Parameters
733 ----------
734 fgcm_obj : `astropy.table.Table`
735 Catalog of per-star ids and positions.
736 """
737 reserved = self.reserve_selection.run(
738 len(fgcm_obj),
739 extra=','.join(self.config.requiredBands),
740 )
741 fgcm_obj["obj_flag"][reserved] |= objFlagDict["RESERVED"]
742
743 def _associate_reference_stars(self, lookup_table_handle, visit_cat, fgcm_obj):
744 """Associate fgcm isolated stars with reference stars.
745
746 Parameters
747 ----------
748 lookup_table_handle : `lsst.daf.butler.DeferredDatasetHandle`, optional
749 Data reference to fgcm look-up table (used if matching reference stars).
750 visit_cat : `lsst.afw.table.ExposureCatalog`
751 Catalog of per-visit data.
752 fgcm_obj : `astropy.table.Table`
753 Catalog of per-star ids and positions
754
755 Returns
756 -------
757 ref_cat : `astropy.table.Table`
758 Catalog of reference stars matched to fgcm stars.
759 """
760 # Figure out the correct bands/filters that we need based on the data.
761 lut_cat = lookup_table_handle.get()
762
763 std_filter_dict = {filter_name: std_filter for (filter_name, std_filter) in
764 zip(lut_cat[0]["physicalFilters"].split(","),
765 lut_cat[0]["stdPhysicalFilters"].split(","))}
766 std_lambda_dict = {std_filter: std_lambda for (std_filter, std_lambda) in
767 zip(lut_cat[0]["stdPhysicalFilters"].split(","),
768 lut_cat[0]["lambdaStdFilter"])}
769 del lut_cat
770
771 reference_filter_names = self._getReferenceFilterNames(
772 visit_cat,
773 std_filter_dict,
774 std_lambda_dict,
775 )
776 self.log.info("Using the following reference filters: %s", (", ".join(reference_filter_names)))
777
778 # Split into healpix pixels for more efficient matching.
779 ipnest = hpg.angle_to_pixel(self.config.coarseNside, fgcm_obj["ra"], fgcm_obj["dec"])
780 hpix, revpix = histogram_rev_sorted(ipnest)
781
782 pixel_cats = []
783
784 # Compute the dtype from the filter names.
785 dtype = [("fgcm_id", "i4"),
786 ("refMag", "f4", (len(reference_filter_names), )),
787 ("refMagErr", "f4", (len(reference_filter_names), ))]
788
789 (gdpix,) = (hpix > 0).nonzero()
790 for ii, gpix in enumerate(gdpix):
791 p1a = revpix[revpix[gpix]: revpix[gpix + 1]]
792
793 # We do a simple wrapping of RA if we need to.
794 ra_wrap = fgcm_obj["ra"][p1a]
795 if (ra_wrap.min() < 10.0) and (ra_wrap.max() > 350.0):
796 ra_wrap[ra_wrap > 180.0] -= 360.0
797 mean_ra = np.mean(ra_wrap) % 360.0
798 else:
799 mean_ra = np.mean(ra_wrap)
800 mean_dec = np.mean(fgcm_obj["dec"][p1a])
801
802 dist = esutil.coords.sphdist(mean_ra, mean_dec, fgcm_obj["ra"][p1a], fgcm_obj["dec"][p1a])
803 rad = dist.max()
804
805 if rad < hpg.nside_to_resolution(self.config.coarseNside)/2.:
806 # Small radius, read just the circle.
807 ref_cat = self.fgcmLoadReferenceCatalog.getFgcmReferenceStarsSkyCircle(
808 mean_ra,
809 mean_dec,
810 rad,
811 reference_filter_names,
812 )
813 else:
814 # Otherwise the less efficient but full coverage.
815 ref_cat = self.fgcmLoadReferenceCatalog.getFgcmReferenceStarsHealpix(
816 self.config.coarseNside,
817 hpg.nest_to_ring(self.config.coarseNside, ipnest[p1a[0]]),
818 reference_filter_names,
819 )
820 if ref_cat.size == 0:
821 # No stars in this pixel; that's okay.
822 continue
823
824 with Matcher(fgcm_obj["ra"][p1a], fgcm_obj["dec"][p1a]) as matcher:
825 idx, i1, i2, d = matcher.query_radius(
826 ref_cat["ra"],
827 ref_cat["dec"],
828 self.config.matchRadius/3600.,
829 return_indices=True,
830 )
831
832 if len(i1) == 0:
833 # No matched stars in this pixel; that's okay.
834 continue
835
836 pixel_cat = Table(data=np.zeros(i1.size, dtype=dtype))
837 pixel_cat["fgcm_id"] = fgcm_obj["fgcm_id"][p1a[i1]]
838 pixel_cat["refMag"][:, :] = ref_cat["refMag"][i2, :]
839 pixel_cat["refMagErr"][:, :] = ref_cat["refMagErr"][i2, :]
840
841 pixel_cats.append(pixel_cat)
842
843 self.log.info(
844 "Found %d reference matches in pixel %d (%d of %d).",
845 len(pixel_cat),
846 ipnest[p1a[0]],
847 ii,
848 gdpix.size - 1,
849 )
850
851 ref_cat_full = vstack(pixel_cats)
852 ref_cat_full.meta.update({'FILTERNAMES': reference_filter_names})
853
854 return ref_cat_full
855
856 def _compute_delta_aper_summary_statistics(self, visit_cat, star_obs):
857 """Compute delta aperture summary statistics (per visit).
858
859 Parameters
860 ----------
861 visit_cat : `lsst.afw.table.ExposureCatalog`
862 Catalog of per-visit data.
863 star_obs : `astropy.table.Table`
864 Catalog of individual star observations.
865 """
866 (ok,) = ((star_obs["delta_mag_aper"] < 99.0)
867 & (star_obs["delta_mag_err_aper"] < 99.0)).nonzero()
868
869 visit_index = np.zeros(len(star_obs[ok]), dtype=np.int32)
870 a, b = esutil.numpy_util.match(visit_cat["visit"], star_obs["visit"][ok])
871 visit_index[b] = a
872
873 h, rev = histogram_rev_sorted(visit_index)
874
875 visit_cat["deltaAper"] = -9999.0
876 h_use, = np.where(h >= 3)
877 for index in h_use:
878 i1a = rev[rev[index]: rev[index + 1]]
879 visit_cat["deltaAper"][index] = np.median(np.asarray(star_obs["delta_mag_aper"][ok[i1a]]))
880
881 n_detector = visit_cat.schema["deltaAperDetector"].asKey().getSize()
882
883 # This is a fast and clever way of doing a multi-dimensional grouping
884 # between visit + detector (hashed together).
885 visit_detector_hash = visit_index * (n_detector + 1) + star_obs["detector"][ok]
886
887 h, rev = histogram_rev_sorted(visit_detector_hash)
888
889 visit_cat["deltaAperDetector"][:] = -9999.0
890 h_use, = np.where(h >= 3)
891 for index in h_use:
892 i1a = rev[rev[index]: rev[index + 1]]
893 v_ind = visit_index[i1a[0]]
894 d_ind = visit_detector_hash[i1a[0]] % (n_detector + 1)
895 delta_mag = np.median(np.asarray(star_obs["delta_mag_aper"][ok[i1a]]))
896 visit_cat["deltaAperDetector"][v_ind, d_ind] = delta_mag
_make_all_star_obs_from_isolated_stars(self, isolated_star_cat_handle_dict, isolated_star_source_handle_dict, visit_cat, camera, calib_flux_aperture_radius=None, use_science_detectors=False)
run(self, *, camera, visit_summary_handle_dict, isolated_star_cat_handle_dict, isolated_star_source_handle_dict, lookup_table_handle)
_getReferenceFilterNames(self, visitCat, stdFilterDict, stdLambdaDict)
fgcmMakeVisitCatalog(self, camera, groupedHandles, useScienceDetectors=False)