LSST Applications g1635faa6d4+215bc75b8c,g1653933729+a8ce1bb630,g22ce9dc20b+d972d8df89,g28da252d5a+0fcf840c6d,g29321ee8c0+e558be0e74,g2bbee38e9b+9634bc57db,g2bc492864f+9634bc57db,g2cdde0e794+c2c89b37c4,g3156d2b45e+41e33cbcdc,g347aa1857d+9634bc57db,g35bb328faa+a8ce1bb630,g3a166c0a6a+9634bc57db,g3e281a1b8c+9f2c4e2fc3,g414038480c+077ccc18e7,g41af890bb2+2a6f257a1d,g5fbc88fb19+17cd334064,g781aacb6e4+a8ce1bb630,g7ab3e175f3+59ce30aec6,g80478fca09+f8b2ab54e1,g82479be7b0+ba9d578ff8,g858d7b2824+d972d8df89,g9125e01d80+a8ce1bb630,g9726552aa6+10f999ec6a,ga5288a1d22+630363936d,gae0086650b+a8ce1bb630,gb58c049af0+d64f4d3760,gb9c6c11c1e+9553554aa7,gbd46683f8f+0c4209622a,gc28159a63d+9634bc57db,gcf0d15dbbd+2db122af0a,gda3e153d99+d972d8df89,gda6a2b7d83+2db122af0a,gdaeeff99f8+1711a396fd,ge2409df99d+d1dc2f3b25,ge33fd446bb+d972d8df89,ge79ae78c31+9634bc57db,gf0baf85859+147a0692ba,gf3967379c6+02b11634a5,w.2024.45
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.pex.config as pexConfig
41import lsst.pipe.base as pipeBase
42from lsst.pipe.base import connectionTypes
43from lsst.meas.algorithms import ReferenceObjectLoader, LoadReferenceObjectsConfig
44from lsst.pipe.tasks.reserveIsolatedStars import ReserveIsolatedStarsTask
45
46from .fgcmBuildStarsBase import FgcmBuildStarsConfigBase, FgcmBuildStarsBaseTask
47from .utilities import computeApproxPixelAreaFields, computeApertureRadiusFromName
48
49__all__ = ["FgcmBuildFromIsolatedStarsConfig", "FgcmBuildFromIsolatedStarsTask"]
50
51
52class FgcmBuildFromIsolatedStarsConnections(pipeBase.PipelineTaskConnections,
53 dimensions=("instrument",),
54 defaultTemplates={}):
55 camera = connectionTypes.PrerequisiteInput(
56 doc="Camera instrument",
57 name="camera",
58 storageClass="Camera",
59 dimensions=("instrument",),
60 isCalibration=True,
61 )
62 fgcm_lookup_table = connectionTypes.PrerequisiteInput(
63 doc=("Atmosphere + instrument look-up-table for FGCM throughput and "
64 "chromatic corrections."),
65 name="fgcmLookUpTable",
66 storageClass="Catalog",
67 dimensions=("instrument",),
68 deferLoad=True,
69 )
70 ref_cat = connectionTypes.PrerequisiteInput(
71 doc="Reference catalog to use for photometric calibration.",
72 name="cal_ref_cat",
73 storageClass="SimpleCatalog",
74 dimensions=("skypix",),
75 deferLoad=True,
76 multiple=True,
77 )
78 isolated_star_cats = pipeBase.connectionTypes.Input(
79 doc=("Catalog of isolated stars with average positions, number of associated "
80 "sources, and indexes to the isolated_star_sources catalogs."),
81 name="isolated_star_presource_associations",
82 storageClass="ArrowAstropy",
83 dimensions=("instrument", "tract", "skymap"),
84 deferLoad=True,
85 multiple=True,
86 )
87 isolated_star_sources = pipeBase.connectionTypes.Input(
88 doc=("Catalog of isolated star sources with sourceIds, and indexes to the "
89 "isolated_star_cats catalogs."),
90 name="isolated_star_presources",
91 storageClass="ArrowAstropy",
92 dimensions=("instrument", "tract", "skymap"),
93 deferLoad=True,
94 multiple=True,
95 )
96 visit_summaries = connectionTypes.Input(
97 doc=("Per-visit consolidated exposure metadata. These catalogs use "
98 "detector id for the id and must be sorted for fast lookups of a "
99 "detector."),
100 name="visitSummary",
101 storageClass="ExposureCatalog",
102 dimensions=("instrument", "visit"),
103 deferLoad=True,
104 multiple=True,
105 )
106 fgcm_visit_catalog = connectionTypes.Output(
107 doc="Catalog of visit information for fgcm.",
108 name="fgcmVisitCatalog",
109 storageClass="Catalog",
110 dimensions=("instrument",),
111 )
112 fgcm_star_observations = connectionTypes.Output(
113 doc="Catalog of star observations for fgcm.",
114 name="fgcm_star_observations",
115 storageClass="ArrowAstropy",
116 dimensions=("instrument",),
117 )
118 fgcm_star_ids = connectionTypes.Output(
119 doc="Catalog of fgcm calibration star IDs.",
120 name="fgcm_star_ids",
121 storageClass="ArrowAstropy",
122 dimensions=("instrument",),
123 )
124 fgcm_reference_stars = connectionTypes.Output(
125 doc="Catalog of fgcm-matched reference stars.",
126 name="fgcm_reference_stars",
127 storageClass="ArrowAstropy",
128 dimensions=("instrument",),
129 )
130
131 def __init__(self, *, config=None):
132 super().__init__(config=config)
133
134 if not config.doReferenceMatches:
135 self.prerequisiteInputs.remove("ref_cat")
136 self.prerequisiteInputs.remove("fgcm_lookup_table")
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.instFluxFieldinstFluxField = 'normCompTophatFlux_instFlux'
164 self.localBackgroundFluxFieldlocalBackgroundFluxField = 'localBackground_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.instFluxFieldinstFluxField
182 source_selector.signalToNoise.errField = self.instFluxFieldinstFluxField + '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 lookup_table_handle = input_ref_dict["fgcm_lookup_table"]
219
220 # Prepare the reference catalog loader
221 ref_config = LoadReferenceObjectsConfig()
222 ref_config.filterMap = self.config.fgcmLoadReferenceCatalog.filterMap
223 ref_obj_loader = ReferenceObjectLoader(dataIds=[ref.datasetRef.dataId
224 for ref in inputRefs.ref_cat],
225 refCats=butlerQC.get(inputRefs.ref_cat),
226 name=self.config.connections.ref_cat,
227 log=self.log,
228 config=ref_config)
229 self.makeSubtask('fgcmLoadReferenceCatalog',
230 refObjLoader=ref_obj_loader,
231 refCatName=self.config.connections.ref_cat)
232 else:
233 lookup_table_handle = None
234
235 # The visit summary handles for use with fgcmMakeVisitCatalog must be keyed with
236 # visit, and values are a list with the first value being the visit_summary_handle,
237 # the second is the source catalog (which is not used, but a place-holder is
238 # left for compatibility reasons.)
239 visit_summary_handle_dict = {handle.dataId['visit']: [handle, None] for
240 handle in input_ref_dict['visit_summaries']}
241
242 camera = input_ref_dict["camera"]
243
244 struct = self.run(
245 camera=camera,
246 visit_summary_handle_dict=visit_summary_handle_dict,
247 isolated_star_cat_handle_dict=isolated_star_cat_handle_dict,
248 isolated_star_source_handle_dict=isolated_star_source_handle_dict,
249 lookup_table_handle=lookup_table_handle,
250 )
251
252 butlerQC.put(struct.fgcm_visit_catalog, outputRefs.fgcm_visit_catalog)
253 butlerQC.put(struct.fgcm_star_observations, outputRefs.fgcm_star_observations)
254 butlerQC.put(struct.fgcm_star_ids, outputRefs.fgcm_star_ids)
255 if self.config.doReferenceMatches:
256 butlerQC.put(struct.fgcm_reference_stars, outputRefs.fgcm_reference_stars)
257
258 def run(self, *, camera, visit_summary_handle_dict, isolated_star_cat_handle_dict,
259 isolated_star_source_handle_dict, lookup_table_handle=None):
260 """Run the fgcmBuildFromIsolatedStarsTask.
261
262 Parameters
263 ----------
264 camera : `lsst.afw.cameraGeom.Camera`
265 Camera object.
266 visit_summary_handle_dict : `dict` [`int`, [`lsst.daf.butler.DeferredDatasetHandle`]]
267 Visit summary dataset handles, with the visit as key.
268 isolated_star_cat_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
269 Isolated star catalog dataset handles, with the tract as key.
270 isolated_star_source_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
271 Isolated star source dataset handles, with the tract as key.
272 lookup_table_handle : `lsst.daf.butler.DeferredDatasetHandle`, optional
273 Data reference to fgcm look-up table (used if matching reference stars).
274
275 Returns
276 -------
277 struct : `lsst.pipe.base.struct`
278 Catalogs for persistence, with attributes:
279
280 ``fgcm_visit_catalog``
281 Catalog of per-visit data (`lsst.afw.table.ExposureCatalog`).
282 ``fgcm_star_observations``
283 Catalog of individual star observations (`astropy.table.Table`).
284 ``fgcm_star_ids``
285 Catalog of per-star ids and positions (`astropy.table.Table`).
286 ``fgcm_reference_stars``
287 Catalog of reference stars matched to fgcm stars (`astropy.table.Table`).
288 """
289 # Compute aperture radius if necessary. This is useful to do now before
290 # any heave lifting has happened (fail early).
291 calib_flux_aperture_radius = None
292 if self.config.doSubtractLocalBackground:
293 try:
294 calib_flux_aperture_radius = computeApertureRadiusFromName(self.config.instFluxField)
295 except RuntimeError as e:
296 raise RuntimeError("Could not determine aperture radius from %s. "
297 "Cannot use doSubtractLocalBackground." %
298 (self.config.instFluxField)) from e
299
300 # Check that we have the lookup_table_handle if we are doing reference matches.
301 if self.config.doReferenceMatches:
302 if lookup_table_handle is None:
303 raise RuntimeError("Must supply lookup_table_handle if config.doReferenceMatches is True.")
304
305 visit_cat = self.fgcmMakeVisitCatalog(camera, visit_summary_handle_dict)
306
307 # Select and concatenate the isolated stars and sources.
308 fgcm_obj, star_obs = self._make_all_star_obs_from_isolated_stars(
309 isolated_star_cat_handle_dict,
310 isolated_star_source_handle_dict,
311 visit_cat,
312 camera,
313 calib_flux_aperture_radius=calib_flux_aperture_radius,
314 )
315
317 visit_cat,
318 star_obs,
319 )
320
321 # Do density down-sampling.
322 self._density_downsample(fgcm_obj, star_obs)
323
324 # Mark reserve stars
325 self._mark_reserve_stars(fgcm_obj)
326
327 # Do reference association.
328 if self.config.doReferenceMatches:
329 fgcm_ref = self._associate_reference_stars(lookup_table_handle, visit_cat, fgcm_obj)
330 else:
331 fgcm_ref = None
332
333 return pipeBase.Struct(
334 fgcm_visit_catalog=visit_cat,
335 fgcm_star_observations=star_obs,
336 fgcm_star_ids=fgcm_obj,
337 fgcm_reference_stars=fgcm_ref,
338 )
339
341 self,
342 isolated_star_cat_handle_dict,
343 isolated_star_source_handle_dict,
344 visit_cat,
345 camera,
346 calib_flux_aperture_radius=None,
347 ):
348 """Make all star observations from isolated star catalogs.
349
350 Parameters
351 ----------
352 isolated_star_cat_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
353 Isolated star catalog dataset handles, with the tract as key.
354 isolated_star_source_handle_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
355 Isolated star source dataset handles, with the tract as key.
356 visit_cat : `lsst.afw.table.ExposureCatalog`
357 Catalog of per-visit data.
358 camera : `lsst.afw.cameraGeom.Camera`
359 The camera object.
360 calib_flux_aperture_radius : `float`, optional
361 Radius for the calibration flux aperture.
362
363 Returns
364 -------
365 fgcm_obj : `astropy.table.Table`
366 Catalog of ids and positions for each star.
367 star_obs : `astropy.table.Table`
368 Catalog of individual star observations.
369 """
370 source_columns = [
371 "sourceId",
372 "visit",
373 "detector",
374 "ra",
375 "dec",
376 "x",
377 "y",
378 "physical_filter",
379 "band",
380 "obj_index",
381 self.config.instFluxField,
382 self.config.instFluxField + "Err",
383 self.config.apertureInnerInstFluxField,
384 self.config.apertureInnerInstFluxField + "Err",
385 self.config.apertureOuterInstFluxField,
386 self.config.apertureOuterInstFluxField + "Err",
387 ]
388
389 if self.config.doSubtractLocalBackground:
390 source_columns.append(self.config.localBackgroundFluxField)
391 local_background_flag_name = self.config.localBackgroundFluxField[0: -len('instFlux')] + 'flag'
392 source_columns.append(local_background_flag_name)
393
394 if self.sourceSelector.config.doFlags:
395 source_columns.extend(self.sourceSelector.config.flags.bad)
396
397 if self.config.doSubtractLocalBackground:
398 local_background_area = np.pi*calib_flux_aperture_radius**2.
399
400 # Compute the approximate pixel area over the full focal plane
401 # from the WCS jacobian using the camera model.
402 approx_pixel_area_fields = computeApproxPixelAreaFields(camera)
403
404 # Construct mapping from detector number to index.
405 detector_mapping = {}
406 for detector_index, detector in enumerate(camera):
407 detector_mapping[detector.getId()] = detector_index
408
409 star_obs_dtype = [
410 ("ra", "f8"),
411 ("dec", "f8"),
412 ("x", "f8"),
413 ("y", "f8"),
414 ("visit", "i8"),
415 ("detector", "i4"),
416 ("inst_mag", "f4"),
417 ("inst_mag_err", "f4"),
418 ("jacobian", "f4"),
419 ("delta_mag_bkg", "f4"),
420 ("delta_mag_aper", "f4"),
421 ("delta_mag_err_aper", "f4"),
422 ]
423
424 fgcm_obj_dtype = [
425 ("fgcm_id", "i8"),
426 ("isolated_star_id", "i8"),
427 ("ra", "f8"),
428 ("dec", "f8"),
429 ("obs_arr_index", "i4"),
430 ("n_obs", "i4"),
431 ("obj_flag", "i4"),
432 ]
433
434 fgcm_objs = []
435 star_obs_cats = []
436 merge_source_counter = 0
437
438 k = 2.5/np.log(10.)
439
440 visit_cat_table = visit_cat.asAstropy()
441
442 for tract in sorted(isolated_star_cat_handle_dict):
443 stars = isolated_star_cat_handle_dict[tract].get()
444 sources = isolated_star_source_handle_dict[tract].get(parameters={"columns": source_columns})
445
446 # Down-select sources.
447 good_sources = self.sourceSelector.selectSources(sources).selected
448 if self.config.doSubtractLocalBackground:
449 good_sources &= (~sources[local_background_flag_name])
450 local_background = local_background_area*sources[self.config.localBackgroundFluxField]
451 good_sources &= ((sources[self.config.instFluxField] - local_background) > 0)
452
453 if good_sources.sum() == 0:
454 self.log.info("No good sources found in tract %d", tract)
455 continue
456
457 # Need to count the observations of each star after cuts, per band.
458 # If we have requiredBands specified, we must make sure that each star
459 # has the minumum number of observations in each of thos bands.
460 # Otherwise, we must make sure that each star has at least the minimum
461 # number of observations in _any_ band.
462 if len(self.config.requiredBands) > 0:
463 loop_bands = self.config.requiredBands
464 else:
465 loop_bands = np.unique(sources["band"])
466
467 n_req = np.zeros((len(loop_bands), len(stars)), dtype=np.int32)
468 for i, band in enumerate(loop_bands):
469 (band_use,) = (sources[good_sources]["band"] == band).nonzero()
470 np.add.at(
471 n_req,
472 (i, sources[good_sources]["obj_index"][band_use]),
473 1,
474 )
475
476 if len(self.config.requiredBands) > 0:
477 # The min gives us the band with the fewest observations, which must be
478 # above the limit.
479 (good_stars,) = (n_req.min(axis=0) >= self.config.minPerBand).nonzero()
480 else:
481 # The max gives us the band with the most observations, which must be
482 # above the limit.
483 (good_stars,) = (n_req.max(axis=0) >= self.config.minPerBand).nonzero()
484
485 if len(good_stars) == 0:
486 self.log.info("No good stars found in tract %d", tract)
487 continue
488
489 # With the following matching:
490 # sources[good_sources][b] <-> stars[good_stars[a]]
491 obj_index = sources["obj_index"][good_sources]
492 a, b = esutil.numpy_util.match(good_stars, obj_index)
493
494 # Update indexes and cut to selected stars/sources
495 _, index_new = np.unique(a, return_index=True)
496 stars["source_cat_index"][good_stars] = index_new
497 sources = sources[good_sources][b]
498 sources["obj_index"][:] = a
499 stars = stars[good_stars]
500
501 nsource = np.zeros(len(stars), dtype=np.int32)
502 np.add.at(
503 nsource,
504 sources["obj_index"],
505 1,
506 )
507 stars["nsource"][:] = nsource
508
509 # After these cuts, the catalogs have the following properties:
510 # - ``stars`` only contains isolated stars that have the minimum number of good
511 # sources in the required bands.
512 # - ``sources`` has been cut to the good sources.
513 # - The slice [stars["source_cat_index"]: stars["source_cat_index"]
514 # + stars["nsource"]]
515 # applied to ``sources`` will give all the sources associated with the star.
516 # - For each source, sources["obj_index"] points to the index of the associated
517 # isolated star.
518
519 # We now reformat the sources and compute the ``observations`` that fgcm expects.
520 star_obs = Table(data=np.zeros(len(sources), dtype=star_obs_dtype))
521 star_obs["ra"] = sources["ra"]
522 star_obs["dec"] = sources["dec"]
523 star_obs["x"] = sources["x"]
524 star_obs["y"] = sources["y"]
525 star_obs["visit"] = sources["visit"]
526 star_obs["detector"] = sources["detector"]
527
528 visit_match, obs_match = esutil.numpy_util.match(visit_cat_table["visit"], sources["visit"])
529
530 exp_time = np.zeros(len(star_obs))
531 exp_time[obs_match] = visit_cat_table["exptime"][visit_match]
532
533 with warnings.catch_warnings():
534 # Ignore warnings, we will filter infinities and nans below.
535 warnings.simplefilter("ignore")
536
537 inst_mag_inner = -2.5*np.log10(sources[self.config.apertureInnerInstFluxField])
538 inst_mag_err_inner = k*(sources[self.config.apertureInnerInstFluxField + "Err"]
539 / sources[self.config.apertureInnerInstFluxField])
540 inst_mag_outer = -2.5*np.log10(sources[self.config.apertureOuterInstFluxField])
541 inst_mag_err_outer = k*(sources[self.config.apertureOuterInstFluxField + "Err"]
542 / sources[self.config.apertureOuterInstFluxField])
543 star_obs["delta_mag_aper"] = inst_mag_inner - inst_mag_outer
544 star_obs["delta_mag_err_aper"] = np.sqrt(inst_mag_err_inner**2. + inst_mag_err_outer**2.)
545 # Set bad values to sentinel values for fgcm.
546 bad = ~np.isfinite(star_obs["delta_mag_aper"])
547 star_obs["delta_mag_aper"][bad] = 99.0
548 star_obs["delta_mag_err_aper"][bad] = 99.0
549
550 if self.config.doSubtractLocalBackground:
551 # At the moment we only adjust the flux and not the flux
552 # error by the background because the error on
553 # base_LocalBackground_instFlux is the rms error in the
554 # background annulus, not the error on the mean in the
555 # background estimate (which is much smaller, by sqrt(n)
556 # pixels used to estimate the background, which we do not
557 # have access to in this task). In the default settings,
558 # the annulus is sufficiently large such that these
559 # additional errors are negligibly small (much less
560 # than a mmag in quadrature).
561
562 # This is the difference between the mag with local background correction
563 # and the mag without local background correction.
564 local_background = local_background_area*sources[self.config.localBackgroundFluxField]
565 star_obs["delta_mag_bkg"] = (-2.5*np.log10(sources[self.config.instFluxField]
566 - local_background) -
567 -2.5*np.log10(sources[self.config.instFluxField]))
568
569 # Need to loop over detectors here.
570 for detector in camera:
571 detector_id = detector.getId()
572 # used index for all observations with a given detector
573 (use,) = (star_obs["detector"][obs_match] == detector_id).nonzero()
574 # Prior to running the calibration, we want to remove the effect
575 # of the jacobian of the WCS because that is a known quantity.
576 # Ideally, this would be done for each individual WCS, but this
577 # is extremely slow and makes small differences that are much
578 # smaller than the variation in the throughput due to other issues.
579 # Therefore, we use the approximate jacobian estimated from the
580 # camera model.
581 jac = approx_pixel_area_fields[detector_id].evaluate(
582 star_obs["x"][obs_match][use],
583 star_obs["y"][obs_match][use],
584 )
585 star_obs["jacobian"][obs_match[use]] = jac
586 scaled_inst_flux = (sources[self.config.instFluxField][obs_match[use]]
587 * visit_cat_table["scaling"][visit_match[use],
588 detector_mapping[detector_id]])
589 star_obs["inst_mag"][obs_match[use]] = (-2.5 * np.log10(scaled_inst_flux
590 / exp_time[use]))
591
592 # Compute instMagErr from inst_flux_err/inst_flux; scaling will cancel out.
593 star_obs["inst_mag_err"] = k*(sources[self.config.instFluxField + "Err"]
594 / sources[self.config.instFluxField])
595
596 # Apply the jacobian if configured to do so.
597 if self.config.doApplyWcsJacobian:
598 star_obs["inst_mag"] -= 2.5*np.log10(star_obs["jacobian"])
599
600 # We now reformat the stars and compute the ''objects'' that fgcm expects.
601 fgcm_obj = Table(data=np.zeros(len(stars), dtype=fgcm_obj_dtype))
602 fgcm_obj["isolated_star_id"] = stars["isolated_star_id"]
603 fgcm_obj["ra"] = stars["ra"]
604 fgcm_obj["dec"] = stars["dec"]
605 fgcm_obj["obs_arr_index"] = stars["source_cat_index"]
606 fgcm_obj["n_obs"] = stars["nsource"]
607
608 # Offset indexes to account for tract merging
609 fgcm_obj["obs_arr_index"] += merge_source_counter
610
611 fgcm_objs.append(fgcm_obj)
612 star_obs_cats.append(star_obs)
613
614 merge_source_counter += len(star_obs)
615
616 fgcm_obj = vstack(fgcm_objs)
617
618 # Set the fgcm_id to a unique 64-bit integer for easier sorting.
619 fgcm_obj["fgcm_id"][:] = np.arange(len(fgcm_obj)) + 1
620
621 return fgcm_obj, vstack(star_obs_cats)
622
623 def _density_downsample(self, fgcm_obj, star_obs):
624 """Downsample stars according to density.
625
626 Catalogs are modified in-place.
627
628 Parameters
629 ----------
630 fgcm_obj : `astropy.table.Table`
631 Catalog of per-star ids and positions.
632 star_obs : `astropy.table.Table`
633 Catalog of individual star observations.
634 """
635 if self.config.randomSeed is not None:
636 np.random.seed(seed=self.config.randomSeed)
637
638 ipnest = hpg.angle_to_pixel(self.config.densityCutNside, fgcm_obj["ra"], fgcm_obj["dec"])
639 # Use the esutil.stat.histogram function to pull out the histogram of
640 # grouped pixels, and the rev_indices which describes which inputs
641 # are grouped together. The fgcm histogram_rev_sorted shim
642 # ensures that the indices are properly sorted.
643 hist, rev_indices = histogram_rev_sorted(ipnest)
644
645 obj_use = np.ones(len(fgcm_obj), dtype=bool)
646
647 (high,) = (hist > self.config.densityCutMaxPerPixel).nonzero()
648 (ok,) = (hist > 0).nonzero()
649 self.log.info("There are %d/%d pixels with high stellar density.", high.size, ok.size)
650 for i in range(high.size):
651 # The pix_indices are the indices of every star in the pixel.
652 pix_indices = rev_indices[rev_indices[high[i]]: rev_indices[high[i] + 1]]
653 # Cut down to the maximum number of stars in the pixel.
654 cut = np.random.choice(
655 pix_indices,
656 size=pix_indices.size - self.config.densityCutMaxPerPixel,
657 replace=False,
658 )
659 obj_use[cut] = False
660
661 fgcm_obj = fgcm_obj[obj_use]
662
663 obs_index = np.zeros(np.sum(fgcm_obj["n_obs"]), dtype=np.int32)
664 ctr = 0
665 for i in range(len(fgcm_obj)):
666 n_obs = fgcm_obj["n_obs"][i]
667 obs_index[ctr: ctr + n_obs] = (
668 np.arange(fgcm_obj["obs_arr_index"][i], fgcm_obj["obs_arr_index"][i] + n_obs)
669 )
670 fgcm_obj["obs_arr_index"][i] = ctr
671 ctr += n_obs
672
673 star_obs = star_obs[obs_index]
674
675 def _mark_reserve_stars(self, fgcm_obj):
676 """Run the star reservation task to flag reserved stars.
677
678 Parameters
679 ----------
680 fgcm_obj : `astropy.table.Table`
681 Catalog of per-star ids and positions.
682 """
683 reserved = self.reserve_selection.run(
684 len(fgcm_obj),
685 extra=','.join(self.config.requiredBands),
686 )
687 fgcm_obj["obj_flag"][reserved] |= objFlagDict["RESERVED"]
688
689 def _associate_reference_stars(self, lookup_table_handle, visit_cat, fgcm_obj):
690 """Associate fgcm isolated stars with reference stars.
691
692 Parameters
693 ----------
694 lookup_table_handle : `lsst.daf.butler.DeferredDatasetHandle`, optional
695 Data reference to fgcm look-up table (used if matching reference stars).
696 visit_cat : `lsst.afw.table.ExposureCatalog`
697 Catalog of per-visit data.
698 fgcm_obj : `astropy.table.Table`
699 Catalog of per-star ids and positions
700
701 Returns
702 -------
703 ref_cat : `astropy.table.Table`
704 Catalog of reference stars matched to fgcm stars.
705 """
706 # Figure out the correct bands/filters that we need based on the data.
707 lut_cat = lookup_table_handle.get()
708
709 std_filter_dict = {filter_name: std_filter for (filter_name, std_filter) in
710 zip(lut_cat[0]["physicalFilters"].split(","),
711 lut_cat[0]["stdPhysicalFilters"].split(","))}
712 std_lambda_dict = {std_filter: std_lambda for (std_filter, std_lambda) in
713 zip(lut_cat[0]["stdPhysicalFilters"].split(","),
714 lut_cat[0]["lambdaStdFilter"])}
715 del lut_cat
716
717 reference_filter_names = self._getReferenceFilterNames(
718 visit_cat,
719 std_filter_dict,
720 std_lambda_dict,
721 )
722 self.log.info("Using the following reference filters: %s", (", ".join(reference_filter_names)))
723
724 # Split into healpix pixels for more efficient matching.
725 ipnest = hpg.angle_to_pixel(self.config.coarseNside, fgcm_obj["ra"], fgcm_obj["dec"])
726 hpix, revpix = histogram_rev_sorted(ipnest)
727
728 pixel_cats = []
729
730 # Compute the dtype from the filter names.
731 dtype = [("fgcm_id", "i4"),
732 ("refMag", "f4", (len(reference_filter_names), )),
733 ("refMagErr", "f4", (len(reference_filter_names), ))]
734
735 (gdpix,) = (hpix > 0).nonzero()
736 for ii, gpix in enumerate(gdpix):
737 p1a = revpix[revpix[gpix]: revpix[gpix + 1]]
738
739 # We do a simple wrapping of RA if we need to.
740 ra_wrap = fgcm_obj["ra"][p1a]
741 if (ra_wrap.min() < 10.0) and (ra_wrap.max() > 350.0):
742 ra_wrap[ra_wrap > 180.0] -= 360.0
743 mean_ra = np.mean(ra_wrap) % 360.0
744 else:
745 mean_ra = np.mean(ra_wrap)
746 mean_dec = np.mean(fgcm_obj["dec"][p1a])
747
748 dist = esutil.coords.sphdist(mean_ra, mean_dec, fgcm_obj["ra"][p1a], fgcm_obj["dec"][p1a])
749 rad = dist.max()
750
751 if rad < hpg.nside_to_resolution(self.config.coarseNside)/2.:
752 # Small radius, read just the circle.
753 ref_cat = self.fgcmLoadReferenceCatalog.getFgcmReferenceStarsSkyCircle(
754 mean_ra,
755 mean_dec,
756 rad,
757 reference_filter_names,
758 )
759 else:
760 # Otherwise the less efficient but full coverage.
761 ref_cat = self.fgcmLoadReferenceCatalog.getFgcmReferenceStarsHealpix(
762 self.config.coarseNside,
763 hpg.nest_to_ring(self.config.coarseNside, ipnest[p1a[0]]),
764 reference_filter_names,
765 )
766 if ref_cat.size == 0:
767 # No stars in this pixel; that's okay.
768 continue
769
770 with Matcher(fgcm_obj["ra"][p1a], fgcm_obj["dec"][p1a]) as matcher:
771 idx, i1, i2, d = matcher.query_radius(
772 ref_cat["ra"],
773 ref_cat["dec"],
774 self.config.matchRadius/3600.,
775 return_indices=True,
776 )
777
778 if len(i1) == 0:
779 # No matched stars in this pixel; that's okay.
780 continue
781
782 pixel_cat = Table(data=np.zeros(i1.size, dtype=dtype))
783 pixel_cat["fgcm_id"] = fgcm_obj["fgcm_id"][p1a[i1]]
784 pixel_cat["refMag"][:, :] = ref_cat["refMag"][i2, :]
785 pixel_cat["refMagErr"][:, :] = ref_cat["refMagErr"][i2, :]
786
787 pixel_cats.append(pixel_cat)
788
789 self.log.info(
790 "Found %d reference matches in pixel %d (%d of %d).",
791 len(pixel_cat),
792 ipnest[p1a[0]],
793 ii,
794 gdpix.size - 1,
795 )
796
797 ref_cat_full = vstack(pixel_cats)
798 ref_cat_full.meta.update({'FILTERNAMES': reference_filter_names})
799
800 return ref_cat_full
801
802 def _compute_delta_aper_summary_statistics(self, visit_cat, star_obs):
803 """Compute delta aperture summary statistics (per visit).
804
805 Parameters
806 ----------
807 visit_cat : `lsst.afw.table.ExposureCatalog`
808 Catalog of per-visit data.
809 star_obs : `astropy.table.Table`
810 Catalog of individual star observations.
811 """
812 (ok,) = ((star_obs["delta_mag_aper"] < 99.0)
813 & (star_obs["delta_mag_err_aper"] < 99.0)).nonzero()
814
815 visit_index = np.zeros(len(star_obs[ok]), dtype=np.int32)
816 a, b = esutil.numpy_util.match(visit_cat["visit"], star_obs["visit"][ok])
817 visit_index[b] = a
818
819 h, rev = histogram_rev_sorted(visit_index)
820
821 visit_cat["deltaAper"] = -9999.0
822 h_use, = np.where(h >= 3)
823 for index in h_use:
824 i1a = rev[rev[index]: rev[index + 1]]
825 visit_cat["deltaAper"][index] = np.median(np.asarray(star_obs["delta_mag_aper"][ok[i1a]]))
run(self, *camera, visit_summary_handle_dict, isolated_star_cat_handle_dict, isolated_star_source_handle_dict, lookup_table_handle=None)
_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)
_getReferenceFilterNames(self, visitCat, stdFilterDict, stdLambdaDict)