Coverage for python/lsst/drp/tasks/gbdesAstrometricFit.py: 71%
1111 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-11-08 02:01 -0800
1# This file is part of drp_tasks.
2#
3# LSST Data Management System
4# This product includes software developed by the
5# LSST Project (http://www.lsst.org/).
6# See COPYRIGHT file at the top of the source tree.
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <https://www.lsstcorp.org/LegalNotices/>.
21#
22import dataclasses
23import itertools
24import re
26import astropy.coordinates
27import astropy.time
28import astropy.units as u
29import astshim
30import numpy as np
31import wcsfit
32import yaml
33from astropy.table import vstack
34from sklearn.cluster import AgglomerativeClustering
35from smatch.matcher import Matcher
37import lsst.afw.geom as afwgeom
38import lsst.afw.table
39import lsst.geom
40import lsst.pex.config as pexConfig
41import lsst.pipe.base as pipeBase
42import lsst.sphgeom
43from lsst.daf.butler import DataCoordinate
44from lsst.meas.algorithms import (
45 LoadReferenceObjectsConfig,
46 ReferenceObjectLoader,
47 ReferenceSourceSelectorTask,
48)
49from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry
51from .build_camera import BuildCameraFromAstrometryTask
53__all__ = [
54 "calculate_apparent_motion",
55 "GbdesAstrometricFitConnections",
56 "GbdesAstrometricFitConfig",
57 "GbdesAstrometricFitTask",
58 "GbdesAstrometricMultibandFitConnections",
59 "GbdesAstrometricMultibandFitTask",
60 "GbdesGlobalAstrometricFitConnections",
61 "GbdesGlobalAstrometricFitConfig",
62 "GbdesGlobalAstrometricFitTask",
63 "GbdesGlobalAstrometricMultibandFitConnections",
64 "GbdesGlobalAstrometricMultibandFitTask",
65]
68def calculate_apparent_motion(catalog, refEpoch):
69 """Calculate shift from reference epoch to the apparent observed position
70 at another date.
72 This function calculates the shift due to proper motion combined with the
73 apparent motion due to parallax. This is not used in the
74 `GbdesAstrometricFitTask` or related child tasks, but is useful for
75 assessing results.
77 Parameters
78 ----------
79 catalog : `astropy.table.Table`
80 Table containing position, proper motion, parallax, and epoch for each
81 source, labeled by columns 'ra', 'dec', 'pmRA', 'pmDec', 'parallax',
82 and 'MJD'.
83 refEpoch : `astropy.time.Time`
84 Epoch of the reference position.
86 Returns
87 -------
88 apparentMotionRA : `np.ndarray` [`astropy.units.Quantity`]
89 RA shift in degrees.
90 apparentMotionDec : `np.ndarray` [`astropy.units.Quantity`]
91 Dec shift in degrees.
92 """
93 ra_rad = catalog["ra"].to(u.rad)
94 dec_rad = catalog["dec"].to(u.rad)
96 dt = (catalog["MJD"] - refEpoch).to(u.yr)
97 properMotionRA = catalog["pmRA"].to(u.deg / u.yr) * dt
98 properMotionDec = catalog["pmDec"].to(u.deg / u.yr) * dt
100 # Just do calculations for unique mjds:
101 mjds = astropy.time.Time(
102 np.unique(catalog["MJD"].mjd), scale=catalog["MJD"][0].scale, format=catalog["MJD"][0].format
103 )
104 sun = astropy.coordinates.get_body("sun", time=mjds)
105 frame = astropy.coordinates.GeocentricTrueEcliptic(equinox=mjds)
106 tmpSunLongitudes = sun.transform_to(frame).lon.radian
108 # Project back to full table:
109 newTable = astropy.table.Table({"MJD": mjds, "SL": tmpSunLongitudes})
110 joined = astropy.table.join(catalog[["MJD"]], newTable, keys="MJD", keep_order=True)
111 sunLongitudes = joined["SL"]
113 # These equations for parallax come from Equations 5.2 in Van de Kamp's
114 # book Stellar Paths. They differ from the parallax calculated in gbdes by
115 # ~0.01 mas, which is acceptable for QA and plotting purposes.
116 parallaxFactorRA = np.cos(wcsfit.EclipticInclination) * np.cos(ra_rad) * np.sin(sunLongitudes) - np.sin(
117 ra_rad
118 ) * np.cos(sunLongitudes)
119 parallaxFactorDec = (
120 np.sin(wcsfit.EclipticInclination) * np.cos(dec_rad)
121 - np.cos(wcsfit.EclipticInclination) * np.sin(ra_rad) * np.sin(dec_rad)
122 ) * np.sin(sunLongitudes) - np.cos(ra_rad) * np.sin(dec_rad) * np.cos(sunLongitudes)
123 parallaxDegrees = catalog["parallax"].to(u.degree)
124 parallaxCorrectionRA = parallaxDegrees * parallaxFactorRA
125 parallaxCorrectionDec = parallaxDegrees * parallaxFactorDec
127 apparentMotionRA = properMotionRA + parallaxCorrectionRA
128 apparentMotionDec = properMotionDec + parallaxCorrectionDec
130 return apparentMotionRA, apparentMotionDec
133def _make_ref_covariance_matrix(
134 refCat, inputUnit=u.radian, outputCoordUnit=u.marcsec, outputPMUnit=u.marcsec, version=1
135):
136 """Make a covariance matrix for the reference catalog including proper
137 motion and parallax.
139 The output is flattened to one dimension to match the format expected by
140 `gbdes`.
142 Parameters
143 ----------
144 refCat : `lsst.afw.table.SimpleCatalog`
145 Catalog including proper motion and parallax measurements.
146 inputUnit : `astropy.unit.core.Unit`
147 Units of the input catalog
148 outputCoordUnit : `astropy.unit.core.Unit`
149 Units required for the coordinates in the covariance matrix. `gbdes`
150 expects milliarcseconds.
151 outputPMUnit : `astropy.unit.core.Unit`
152 Units required for the proper motion/parallax in the covariance matrix.
153 `gbdes` expects milliarcseconds.
154 version : `int`
155 Version of the reference catalog. Version 2 includes covariance
156 measurements.
157 Returns
158 -------
159 cov : `list` [`float`]
160 Flattened output covariance matrix.
161 """
162 cov = np.zeros((len(refCat), 25))
163 if version == 1:
164 # Here is the standard ordering of components in the cov matrix,
165 # to match the PM enumeration in C++ code of gbdes package's Match.
166 # Each tuple gives: the array holding the 1d error,
167 # the string in Gaia column names for this
168 # the ordering in the Gaia catalog
169 # and the ordering of the tuples is the order we want in our cov matrix
170 raErr = (refCat["coord_raErr"]).to(outputCoordUnit).to_value()
171 decErr = (refCat["coord_decErr"]).to(outputCoordUnit).to_value()
172 raPMErr = (refCat["pm_raErr"]).to(outputPMUnit).to_value()
173 decPMErr = (refCat["pm_decErr"]).to(outputPMUnit).to_value()
174 parallaxErr = (refCat["parallaxErr"]).to(outputPMUnit).to_value()
175 stdOrder = (
176 (raErr, "ra", 0),
177 (decErr, "dec", 1),
178 (raPMErr, "pmra", 3),
179 (decPMErr, "pmdec", 4),
180 (parallaxErr, "parallax", 2),
181 )
183 k = 0
184 for i, pr1 in enumerate(stdOrder):
185 for j, pr2 in enumerate(stdOrder):
186 if pr1[2] < pr2[2]:
187 cov[:, k] = 0
188 elif pr1[2] > pr2[2]:
189 cov[:, k] = 0
190 else:
191 # diagnonal element
192 cov[:, k] = pr1[0] * pr2[0]
193 k = k + 1
195 elif version == 2:
196 positionParameters = ["coord_ra", "coord_dec", "pm_ra", "pm_dec", "parallax"]
197 units = [outputCoordUnit, outputCoordUnit, outputPMUnit, outputPMUnit, outputPMUnit]
198 k = 0
199 for i, pi in enumerate(positionParameters):
200 for j, pj in enumerate(positionParameters):
201 if i == j:
202 cov[:, k] = ((refCat[f"{pi}Err"].value) ** 2 * inputUnit**2).to(units[j] * units[j]).value
203 elif i > j:
204 cov[:, k] = (refCat[f"{pj}_{pi}_Cov"].value * inputUnit**2).to_value(units[i] * units[j])
205 else:
206 cov[:, k] = (refCat[f"{pi}_{pj}_Cov"].value * inputUnit**2).to_value(units[i] * units[j])
207 k += 1
208 return cov
211def _nCoeffsFromDegree(degree):
212 """Get the number of coefficients for a polynomial of a certain degree with
213 two variables.
215 This uses the general formula that the number of coefficients for a
216 polynomial of degree d with n variables is (n + d) choose d, where in this
217 case n is fixed to 2.
219 Parameters
220 ----------
221 degree : `int`
222 Degree of the polynomial in question.
224 Returns
225 -------
226 nCoeffs : `int`
227 Number of coefficients for the polynomial in question.
228 """
229 nCoeffs = int((degree + 2) * (degree + 1) / 2)
230 return nCoeffs
233def _degreeFromNCoeffs(nCoeffs):
234 """Get the degree for a polynomial with two variables and a certain number
235 of coefficients.
237 This is done by applying the quadratic formula to the
238 formula for calculating the number of coefficients of the polynomial.
240 Parameters
241 ----------
242 nCoeffs : `int`
243 Number of coefficients for the polynomial in question.
245 Returns
246 -------
247 degree : `int`
248 Degree of the polynomial in question.
249 """
250 degree = int(-1.5 + 0.5 * (1 + 8 * nCoeffs) ** 0.5)
251 return degree
254def _convert_to_ast_polymap_coefficients(coefficients):
255 """Convert vector of polynomial coefficients from the format used in
256 `gbdes` into AST format (see Poly2d::vectorIndex(i, j) in
257 gbdes/gbutil/src/Poly2d.cpp). This assumes two input and two output
258 coordinates.
260 Parameters
261 ----------
262 coefficients : `list`
263 Coefficients of the polynomials.
264 degree : `int`
265 Degree of the polynomial.
267 Returns
268 -------
269 astPoly : `astshim.PolyMap`
270 Coefficients in AST polynomial format.
271 """
272 polyArray = np.zeros((len(coefficients), 4))
273 N = len(coefficients) / 2
274 degree = _degreeFromNCoeffs(N)
276 for outVar in [1, 2]:
277 for i in range(degree + 1):
278 for j in range(degree + 1):
279 if (i + j) > degree:
280 continue
281 vectorIndex = int(((i + j) * (i + j + 1)) / 2 + j + N * (outVar - 1))
282 polyArray[vectorIndex, 0] = coefficients[vectorIndex]
283 polyArray[vectorIndex, 1] = outVar
284 polyArray[vectorIndex, 2] = i
285 polyArray[vectorIndex, 3] = j
287 astPoly = astshim.PolyMap(polyArray, 2, options="IterInverse=1,NIterInverse=10,TolInverse=1e-7")
288 return astPoly
291def _get_instruments(inputVisitSummaries):
292 """Make `wcsfit.Instrument` objects for all of the instruments and filters
293 used for the input visits. This also returns the indices to match the
294 visits to the instrument/filter used.
296 Parameters
297 ----------
298 inputVisitSummaries: `list` [`lsst.afw.table.ExposureCatalog`]
299 List of catalogs with per-detector summary information.
301 Returns
302 -------
303 instruments : `list` [`wcsfit.Instrument`]
304 List of instrument objects.
305 instrumentIndices : `list` [`int`]
306 Indices matching each visit to the instrument/filter used.
307 """
308 instruments = []
309 filters = []
310 instrumentIndices = []
311 for visitSummary in inputVisitSummaries:
312 filter = visitSummary[0]["physical_filter"]
313 instrumentName = visitSummary[0].getVisitInfo().instrumentLabel
314 if filter not in filters:
315 filters.append(filter)
316 filter_instrument = wcsfit.Instrument(instrumentName)
317 filter_instrument.band = filter
318 instruments.append(filter_instrument)
319 instrumentIndices.append(filters.index(filter))
320 return instruments, instrumentIndices
323class CholeskyError(pipeBase.AlgorithmError):
324 """Raised if the Cholesky decomposition in the model fit fails."""
326 def __init__(self) -> None:
327 super().__init__(
328 "Cholesky decomposition failed, likely because data is not sufficient to constrain the model."
329 )
331 @property
332 def metadata(self) -> dict:
333 """There is no metadata associated with this error."""
334 return {}
337class GbdesAstrometricFitConnections(
338 pipeBase.PipelineTaskConnections,
339 dimensions=("skymap", "tract", "instrument", "physical_filter"),
340 defaultTemplates={
341 "outputName": "gbdesAstrometricFit",
342 },
343):
344 """Middleware input/output connections for task data."""
346 inputCatalogRefs = pipeBase.connectionTypes.Input(
347 doc="Source table in parquet format, per visit.",
348 name="preSourceTable_visit",
349 storageClass="DataFrame",
350 dimensions=("instrument", "visit"),
351 deferLoad=True,
352 multiple=True,
353 )
354 inputVisitSummaries = pipeBase.connectionTypes.Input(
355 doc=(
356 "Per-visit consolidated exposure metadata built from calexps. "
357 "These catalogs use detector id for the id and must be sorted for "
358 "fast lookups of a detector."
359 ),
360 name="visitSummary",
361 storageClass="ExposureCatalog",
362 dimensions=("instrument", "visit"),
363 multiple=True,
364 )
365 referenceCatalog = pipeBase.connectionTypes.PrerequisiteInput(
366 doc="The astrometry reference catalog to match to loaded input catalog sources.",
367 name="gaia_dr3_20230707",
368 storageClass="SimpleCatalog",
369 dimensions=("skypix",),
370 deferLoad=True,
371 multiple=True,
372 )
373 isolatedStarSources = pipeBase.connectionTypes.Input(
374 doc="Catalog of matched sources.",
375 name="isolated_star_presources",
376 storageClass="DataFrame",
377 dimensions=("instrument", "skymap", "tract"),
378 multiple=True,
379 deferLoad=True,
380 )
381 isolatedStarCatalogs = pipeBase.connectionTypes.Input(
382 doc="Catalog of objects corresponding to the isolatedStarSources.",
383 name="isolated_star_presource_associations",
384 storageClass="DataFrame",
385 dimensions=("instrument", "skymap", "tract"),
386 multiple=True,
387 deferLoad=True,
388 )
389 colorCatalog = pipeBase.connectionTypes.Input(
390 doc="The catalog of magnitudes to match to input sources.",
391 name="fgcm_Cycle4_StandardStars",
392 storageClass="SimpleCatalog",
393 dimensions=("instrument",),
394 )
395 inputCameraModel = pipeBase.connectionTypes.PrerequisiteInput(
396 doc="Camera parameters to use for 'device' part of model",
397 name="gbdesAstrometricFit_cameraModel",
398 storageClass="ArrowNumpyDict",
399 dimensions=("instrument", "physical_filter"),
400 )
401 inputCamera = pipeBase.connectionTypes.PrerequisiteInput(
402 doc="Input camera to use when constructing camera from astrometric model.",
403 name="camera",
404 storageClass="Camera",
405 dimensions=("instrument",),
406 isCalibration=True,
407 )
408 outputWcs = pipeBase.connectionTypes.Output(
409 doc=(
410 "Per-tract, per-visit world coordinate systems derived from the fitted model."
411 " These catalogs only contain entries for detectors with an output, and use"
412 " the detector id for the catalog id, sorted on id for fast lookups of a detector."
413 ),
414 name="{outputName}SkyWcsCatalog",
415 storageClass="ExposureCatalog",
416 dimensions=("instrument", "visit", "skymap", "tract"),
417 multiple=True,
418 )
419 outputCatalog = pipeBase.connectionTypes.Output(
420 doc=(
421 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
422 "plane coordinates and chisq values."
423 ),
424 name="{outputName}_fitStars",
425 storageClass="ArrowNumpyDict",
426 dimensions=("instrument", "skymap", "tract", "physical_filter"),
427 )
428 starCatalog = pipeBase.connectionTypes.Output(
429 doc=(
430 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
431 "fitProperMotion is True."
432 ),
433 name="{outputName}_starCatalog",
434 storageClass="ArrowNumpyDict",
435 dimensions=("instrument", "skymap", "tract", "physical_filter"),
436 )
437 modelParams = pipeBase.connectionTypes.Output(
438 doc="WCS parameters and covariance.",
439 name="{outputName}_modelParams",
440 storageClass="ArrowNumpyDict",
441 dimensions=("instrument", "skymap", "tract", "physical_filter"),
442 )
443 outputCameraModel = pipeBase.connectionTypes.Output(
444 doc="Camera parameters to use for 'device' part of model",
445 name="{outputName}_cameraModel",
446 storageClass="ArrowNumpyDict",
447 dimensions=("instrument", "skymap", "tract", "physical_filter"),
448 )
449 camera = pipeBase.connectionTypes.Output(
450 doc="Camera object constructed using the per-detector part of the astrometric model",
451 name="{outputName}Camera",
452 storageClass="Camera",
453 dimensions=("instrument", "skymap", "tract", "physical_filter"),
454 )
455 dcrCoefficients = pipeBase.connectionTypes.Output(
456 doc="Per-visit coefficients for DCR correction.",
457 name="{outputName}_dcrCoefficients",
458 storageClass="ArrowNumpyDict",
459 dimensions=("instrument", "skymap", "tract", "physical_filter"),
460 )
462 def getSpatialBoundsConnections(self):
463 return ("inputVisitSummaries",)
465 def adjust_all_quanta(self, adjuster: pipeBase.QuantaAdjuster) -> None:
466 if not self.config.useIsolatedStars:
467 return
468 # The goal of this implementation of the hook is to expand the set
469 # of isolatedStar* inputs to the tracts that overlap the visits that
470 # overlap the quantum, as opposed to just the tracts that overlap the
471 # quantum.
472 #
473 # We start by extracting mappings from quantum data ID (i.e. tract or
474 # healpix) to tract data ID and visit data ID from the naive-overlap
475 # QG we're starting from.
476 tracts_by_quantum: dict[DataCoordinate, set[DataCoordinate]] = {}
477 visits_by_quantum: dict[DataCoordinate, set[DataCoordinate]] = {}
478 for quantum_data_id in adjuster.iter_data_ids():
479 quantum_inputs = adjuster.get_inputs(quantum_data_id)
480 # We assume we have the exact same tracts for isolatedStarCatalogs
481 # and isolatedStar sources; something is very wrong if we do not.
482 tracts_by_quantum[quantum_data_id] = set(quantum_inputs["isolatedStarCatalogs"])
483 visits_by_quantum[quantum_data_id] = set(quantum_inputs["inputVisitSummaries"])
484 # In order to get from quantum data ID to visit ID to tract data ID,
485 # we next need a mapping from visit to tract. Here's a start, with
486 # just the keys populated so far.:
487 tracts_by_visit: dict[DataCoordinate, set[DataCoordinate]] = {
488 visit_data_id: set()
489 for visit_data_id in itertools.chain.from_iterable(visits_by_quantum.values())
490 }
491 if self.config.healpix is None:
492 # If quanta are per-tract, we just have to invert
493 # visits_by_quantum.
494 for tract_data_id, visit_data_ids_for_quantum in visits_by_quantum.items():
495 for visit_data_id in visit_data_ids_for_quantum:
496 tracts_by_visit[visit_data_id].add(tract_data_id)
497 else:
498 # If quanta are per-healpix, we have to query the butler for the
499 # tract/visit overlaps.
500 tracts_flat = set(itertools.chain.from_iterable(tracts_by_quantum.values()))
501 tract_dimensions = adjuster.butler.dimensions.conform(["tract", "instrument"])
502 visit_dimensions = adjuster.butler.dimensions.conform(["visit"])
503 with adjuster.butler.query() as query:
504 for joint_data_id in (
505 query.join_data_coordinates(tracts_flat)
506 .join_data_coordinates(tracts_by_visit.keys())
507 .data_ids(["tract", "visit"])
508 ):
509 tract_data_id = joint_data_id.subset(tract_dimensions)
510 visit_data_id = joint_data_id.subset(visit_dimensions)
511 tracts_by_visit[visit_data_id].add(tract_data_id)
512 # Loop over quantum data IDs and then visits and then tracts to extend
513 # the inputs of each quantum.
514 for quantum_data_id, visit_data_ids_for_quantum in visits_by_quantum.items():
515 # We'll use tracts_for_quantum to track which inputs are already
516 # present.
517 tracts_for_quantum = tracts_by_quantum[quantum_data_id]
518 for visit_data_id in visit_data_ids_for_quantum:
519 for tract_data_id in tracts_by_visit[visit_data_id]:
520 if tract_data_id not in tracts_for_quantum:
521 adjuster.add_input(quantum_data_id, "isolatedStarCatalogs", tract_data_id)
522 adjuster.add_input(quantum_data_id, "isolatedStarSources", tract_data_id)
523 tracts_for_quantum.add(tract_data_id)
525 def __init__(self, *, config=None):
526 super().__init__(config=config)
528 if self.config.healpix is not None:
529 self.dimensions.remove("tract")
530 self.dimensions.remove("skymap")
531 healpixName = f"healpix{self.config.healpix}"
532 self.dimensions.add(healpixName)
533 self.outputWcs = dataclasses.replace(
534 self.outputWcs, dimensions=("instrument", "visit", healpixName)
535 )
536 self.outputCatalog = dataclasses.replace(
537 self.outputCatalog, dimensions=("instrument", "physical_filter", healpixName)
538 )
539 self.starCatalog = dataclasses.replace(
540 self.starCatalog, dimensions=("instrument", "physical_filter", healpixName)
541 )
542 self.modelParams = dataclasses.replace(
543 self.modelParams, dimensions=("instrument", "physical_filter", healpixName)
544 )
545 self.outputCameraModel = dataclasses.replace(
546 self.outputCameraModel, dimensions=("instrument", "physical_filter", healpixName)
547 )
548 self.camera = dataclasses.replace(
549 self.camera, dimensions=("instrument", "physical_filter", healpixName)
550 )
551 self.dcrCoefficients = dataclasses.replace(
552 self.dcrCoefficients, dimensions=("instrument", "physical_filter", healpixName)
553 )
555 if not self.config.useColor:
556 self.inputs.remove("colorCatalog")
557 self.outputs.remove("dcrCoefficients")
558 if not self.config.saveModelParams:
559 self.outputs.remove("modelParams")
560 if not self.config.useInputCameraModel:
561 self.prerequisiteInputs.remove("inputCameraModel")
562 if not self.config.saveCameraModel:
563 self.outputs.remove("outputCameraModel")
564 if not self.config.saveCameraObject:
565 self.prerequisiteInputs.remove("inputCamera")
566 self.outputs.remove("camera")
567 if self.config.useIsolatedStars:
568 del self.inputCatalogRefs
569 else:
570 del self.isolatedStarCatalogs
571 del self.isolatedStarSources
574class GbdesAstrometricFitConfig(
575 pipeBase.PipelineTaskConfig, pipelineConnections=GbdesAstrometricFitConnections
576):
577 """Configuration for GbdesAstrometricFitTask"""
579 sourceSelector = sourceSelectorRegistry.makeField(
580 doc="How to select sources for cross-matching.", default="science"
581 )
582 referenceSelector = pexConfig.ConfigurableField(
583 target=ReferenceSourceSelectorTask,
584 doc="How to down-select the loaded astrometry reference catalog.",
585 )
586 referenceFilter = pexConfig.Field(
587 dtype=str,
588 doc="Name of filter to load from reference catalog. This is a required argument, although the values"
589 "returned are not used.",
590 default="phot_g_mean",
591 )
592 setRefEpoch = pexConfig.Field(
593 dtype=float,
594 doc="Set the reference epoch to a fixed value in MJD (if None, median observation date is used)",
595 default=None,
596 optional=True,
597 )
598 applyRefCatProperMotion = pexConfig.Field(
599 dtype=bool,
600 doc="Apply proper motion to shift reference catalog to epoch of observations.",
601 default=True,
602 )
603 matchRadius = pexConfig.Field(
604 doc="Matching tolerance between associated objects (arcseconds).", dtype=float, default=1.0
605 )
606 minMatches = pexConfig.Field(
607 doc="Number of matches required to keep a source object.", dtype=int, default=2
608 )
609 allowSelfMatches = pexConfig.Field(
610 doc="Allow multiple sources from the same visit to be associated with the same object.",
611 dtype=bool,
612 default=False,
613 )
614 sourceFluxType = pexConfig.Field(
615 dtype=str,
616 doc="Source flux field to use in source selection and to get fluxes from the catalog.",
617 default="apFlux_12_0",
618 )
619 systematicError = pexConfig.Field(
620 dtype=float,
621 doc=(
622 "Systematic error padding added in quadrature for the science catalogs (marcsec). The default"
623 "value is equivalent to 0.02 pixels for HSC."
624 ),
625 default=0.0034,
626 )
627 referenceSystematicError = pexConfig.Field(
628 dtype=float,
629 doc="Systematic error padding added in quadrature for the reference catalog (marcsec).",
630 default=0.0,
631 )
632 useColor = pexConfig.Field(
633 dtype=bool,
634 doc="Use color information to correct for differential chromatic refraction.",
635 default=False,
636 )
637 color = pexConfig.ListField(
638 dtype=str,
639 doc="The bands to use for calculating color.",
640 default=["g", "i"],
641 listCheck=(lambda x: (len(x) == 2) and (len(set(x)) == len(x))),
642 )
643 referenceColor = pexConfig.Field(
644 dtype=float,
645 doc="The color for which DCR is defined as zero.",
646 default=0.61,
647 )
648 modelComponents = pexConfig.ListField(
649 dtype=str,
650 doc=(
651 "List of mappings to apply to transform from pixels to sky, in order of their application."
652 "Supported options are 'INSTRUMENT/DEVICE' and 'EXPOSURE'."
653 ),
654 default=["INSTRUMENT/DEVICE", "EXPOSURE"],
655 )
656 deviceModel = pexConfig.ListField(
657 dtype=str,
658 doc=(
659 "List of mappings to apply to transform from detector pixels to intermediate frame. Map names"
660 "should match the format 'BAND/DEVICE/<map name>'."
661 ),
662 default=["BAND/DEVICE/poly"],
663 )
664 exposureModel = pexConfig.ListField(
665 dtype=str,
666 doc=(
667 "List of mappings to apply to transform from intermediate frame to sky coordinates. Map names"
668 "should match the format 'EXPOSURE/<map name>'."
669 ),
670 default=["EXPOSURE/poly"],
671 )
672 devicePolyOrder = pexConfig.Field(dtype=int, doc="Order of device polynomial model.", default=4)
673 exposurePolyOrder = pexConfig.Field(dtype=int, doc="Order of exposure polynomial model.", default=6)
674 fitProperMotion = pexConfig.Field(dtype=bool, doc="Fit the proper motions of the objects.", default=False)
675 excludeNonPMObjects = pexConfig.Field(
676 dtype=bool, doc="Exclude reference objects without proper motion/parallax information.", default=True
677 )
678 fitReserveFraction = pexConfig.Field(
679 dtype=float, default=0.2, doc="Fraction of objects to reserve from fit for validation."
680 )
681 fitReserveRandomSeed = pexConfig.Field(
682 dtype=int,
683 doc="Set the random seed for selecting data points to reserve from the fit for validation.",
684 default=1234,
685 )
686 saveModelParams = pexConfig.Field(
687 dtype=bool,
688 doc=(
689 "Save the parameters and covariance of the WCS model. Default to "
690 "false because this can be very large."
691 ),
692 default=False,
693 )
694 useInputCameraModel = pexConfig.Field(
695 dtype=bool,
696 doc=(
697 "Use a preexisting model for the 'device' part of the model. When true, the device part of the"
698 " model will be held fixed in the fitting process."
699 ),
700 default=False,
701 )
702 saveCameraModel = pexConfig.Field(
703 dtype=bool,
704 doc="Save the 'device' part of the model to be used as input in future runs.",
705 default=False,
706 )
707 buildCamera = pexConfig.ConfigurableField(
708 target=BuildCameraFromAstrometryTask, doc="Subtask to build camera from astrometric model."
709 )
710 saveCameraObject = pexConfig.Field(
711 dtype=bool,
712 doc="Build and output an lsst.afw.cameraGeom.Camera object using the fit per-detector model.",
713 default=False,
714 )
715 clipThresh = pexConfig.Field(
716 dtype=float,
717 doc="Threshold for clipping outliers in the fit (in standard deviations)",
718 default=5.0,
719 )
720 clipFraction = pexConfig.Field(
721 dtype=float,
722 doc="Minimum fraction of clipped sources that triggers a new fit iteration.",
723 default=0.0,
724 )
725 healpix = pexConfig.Field(
726 dtype=int,
727 doc="Run using all visits overlapping a healpix pixel with this order instead of a tract. Order 3 "
728 "corresponds to pixels with angular size of 7.329 degrees.",
729 optional=True,
730 default=None,
731 )
732 minDetectorFraction = pexConfig.Field(
733 dtype=float,
734 doc=(
735 "Minimum fraction of detectors with valid WCSs per visit required to include the visit in the "
736 "fit."
737 ),
738 default=0.25,
739 )
740 useIsolatedStars = pexConfig.Field(
741 dtype=bool,
742 default=False,
743 doc=(
744 "If True, use the pre-matched isolated star catalogs instead of loading and matching per-visit "
745 "input catalogs."
746 ),
747 )
749 def setDefaults(self):
750 # Use only stars because aperture fluxes of galaxies are biased and
751 # depend on seeing.
752 self.sourceSelector["science"].doUnresolved = True
753 self.sourceSelector["science"].unresolved.name = "sizeExtendedness"
755 # Use only isolated sources.
756 self.sourceSelector["science"].doIsolated = True
757 self.sourceSelector["science"].isolated.parentName = "parentSourceId"
758 self.sourceSelector["science"].isolated.nChildName = "deblend_nChild"
759 # Do not use either flux or centroid measurements with flags,
760 # chosen from the usual QA flags for stars.
761 self.sourceSelector["science"].doFlags = True
762 badFlags = [
763 "pixelFlags_edge",
764 "pixelFlags_saturated",
765 "pixelFlags_interpolatedCenter",
766 "pixelFlags_interpolated",
767 "pixelFlags_crCenter",
768 "pixelFlags_bad",
769 "hsmPsfMoments_flag",
770 f"{self.sourceFluxType}_flag",
771 ]
772 self.sourceSelector["science"].flags.bad = badFlags
774 # Use only primary sources.
775 self.sourceSelector["science"].doRequirePrimary = True
777 self.sourceSelector["science"].doSignalToNoise = True
778 self.sourceSelector["science"].signalToNoise.minimum = 8.0
779 self.sourceSelector["science"].signalToNoise.maximum = 1000.0
780 self.sourceSelector["science"].signalToNoise.fluxField = self.sourceFluxType + "_instFlux"
781 self.sourceSelector["science"].signalToNoise.errField = self.sourceFluxType + "_instFluxErr"
783 def validate(self):
784 super().validate()
786 # Check if all components of the device and exposure models are
787 # supported.
788 for component in self.deviceModel:
789 mapping = component.split("/")[-1]
790 if mapping not in ["poly", "identity"]: 790 ↛ 791line 790 didn't jump to line 791 because the condition on line 790 was never true
791 raise pexConfig.FieldValidationError(
792 GbdesAstrometricFitConfig.deviceModel,
793 self,
794 f"deviceModel component {component} is not supported.",
795 )
797 for component in self.exposureModel:
798 mapping = component.split("/")[-1]
799 if mapping not in ["poly", "identity", "dcr"]: 799 ↛ 800line 799 didn't jump to line 800 because the condition on line 799 was never true
800 raise pexConfig.FieldValidationError(
801 GbdesAstrometricFitConfig.exposureModel,
802 self,
803 f"exposureModel component {component} is not supported.",
804 )
806 if self.saveCameraModel and self.useInputCameraModel: 806 ↛ 807line 806 didn't jump to line 807 because the condition on line 806 was never true
807 raise pexConfig.FieldValidationError(
808 GbdesAstrometricFitConfig.saveCameraModel,
809 self,
810 "saveCameraModel and useInputCameraModel cannot both be true.",
811 )
813 if self.saveCameraObject and (self.exposurePolyOrder != 1): 813 ↛ 814line 813 didn't jump to line 814 because the condition on line 813 was never true
814 raise pexConfig.FieldValidationError(
815 GbdesAstrometricFitConfig.saveCameraObject,
816 self,
817 "If saveCameraObject is True, exposurePolyOrder must be set to 1.",
818 )
821class GbdesAstrometricFitTask(pipeBase.PipelineTask):
822 """Calibrate the WCS across multiple visits of the same field using the
823 GBDES package.
824 """
826 ConfigClass = GbdesAstrometricFitConfig
827 _DefaultName = "gbdesAstrometricFit"
829 def __init__(self, **kwargs):
830 super().__init__(**kwargs)
831 self.makeSubtask("sourceSelector")
832 self.makeSubtask("referenceSelector")
833 if self.config.saveCameraObject: 833 ↛ 834line 833 didn't jump to line 834 because the condition on line 833 was never true
834 self.makeSubtask("buildCamera")
836 def runQuantum(self, butlerQC, inputRefs, outputRefs):
837 # We override runQuantum to set up the refObjLoaders
838 inputs = butlerQC.get(inputRefs)
840 instrumentName = butlerQC.quantum.dataId["instrument"]
842 # Ensure the inputs are in a consistent and deterministic order
843 if self.config.useIsolatedStars:
844 isolatedStarCatalogs = inputs.pop("isolatedStarCatalogs")
845 inputCatTracts = np.array([inputCat.dataId["tract"] for inputCat in isolatedStarCatalogs])
846 isolatedStarCatalogs = [isolatedStarCatalogs[v] for v in inputCatTracts.argsort()]
848 isolatedStarSources = inputs.pop("isolatedStarSources")
849 inputSourceTracts = np.array([inputCat.dataId["tract"] for inputCat in isolatedStarSources])
850 isolatedStarSources = [isolatedStarSources[v] for v in inputSourceTracts.argsort()]
851 inputs["inputCatalogRefs"] = None
852 else:
853 inputCatVisits = np.array([inputCat.dataId["visit"] for inputCat in inputs["inputCatalogRefs"]])
854 inputs["inputCatalogRefs"] = [inputs["inputCatalogRefs"][v] for v in inputCatVisits.argsort()]
855 isolatedStarCatalogs = None
856 isolatedStarSources = None
857 inputSumVisits = np.array([inputSum[0]["visit"] for inputSum in inputs["inputVisitSummaries"]])
858 inputs["inputVisitSummaries"] = [inputs["inputVisitSummaries"][v] for v in inputSumVisits.argsort()]
859 inputRefHtm7s = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputRefs.referenceCatalog])
860 inputRefCatRefs = [inputRefs.referenceCatalog[htm7] for htm7 in inputRefHtm7s.argsort()]
861 inputRefCats = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputs["referenceCatalog"]])
862 inputs["referenceCatalog"] = [inputs["referenceCatalog"][v] for v in inputRefCats.argsort()]
864 refConfig = LoadReferenceObjectsConfig()
865 if self.config.applyRefCatProperMotion:
866 refConfig.requireProperMotion = True
867 refObjectLoader = ReferenceObjectLoader(
868 dataIds=[ref.datasetRef.dataId for ref in inputRefCatRefs],
869 refCats=inputs.pop("referenceCatalog"),
870 config=refConfig,
871 log=self.log,
872 )
874 nCores = butlerQC.resources.num_cores
875 self.log.info("Running with nCores = %d", nCores)
877 if self.config.useColor:
878 colorCatalog = inputs.pop("colorCatalog")
879 else:
880 colorCatalog = None
881 try:
882 output = self.run(
883 **inputs,
884 isolatedStarCatalogs=isolatedStarCatalogs,
885 isolatedStarSources=isolatedStarSources,
886 instrumentName=instrumentName,
887 refObjectLoader=refObjectLoader,
888 colorCatalog=colorCatalog,
889 nCores=nCores,
890 )
891 except pipeBase.AlgorithmError as e:
892 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
893 # No partial outputs for butler to put
894 raise error from e
896 wcsOutputRefDict = {outWcsRef.dataId["visit"]: outWcsRef for outWcsRef in outputRefs.outputWcs}
897 for visit, outputWcs in output.outputWcss.items():
898 butlerQC.put(outputWcs, wcsOutputRefDict[visit])
899 butlerQC.put(output.outputCatalog, outputRefs.outputCatalog)
900 butlerQC.put(output.starCatalog, outputRefs.starCatalog)
901 if self.config.saveModelParams:
902 butlerQC.put(output.modelParams, outputRefs.modelParams)
903 if self.config.saveCameraModel:
904 butlerQC.put(output.cameraModelParams, outputRefs.outputCameraModel)
905 if self.config.saveCameraObject:
906 butlerQC.put(output.camera, outputRefs.camera)
907 if self.config.useColor:
908 butlerQC.put(output.colorParams, outputRefs.dcrCoefficients)
909 if output.partialOutputs:
910 e = RuntimeError("Some visits were dropped because data was insufficient to fit model.")
911 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
912 raise error from e
914 def run(
915 self,
916 inputCatalogRefs,
917 inputVisitSummaries,
918 isolatedStarSources=None,
919 isolatedStarCatalogs=None,
920 instrumentName="",
921 refEpoch=None,
922 refObjectLoader=None,
923 inputCameraModel=None,
924 colorCatalog=None,
925 inputCamera=None,
926 nCores=1,
927 ):
928 """Run the WCS fit for a given set of visits
930 Parameters
931 ----------
932 inputCatalogRefs : `list` [`DeferredDatasetHandle`]
933 List of handles pointing to visit-level source
934 tables.
935 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
936 List of catalogs with per-detector summary information.
937 instrumentName : `str`, optional
938 Name of the instrument used. This is only used for labelling.
939 refEpoch : `float`
940 Epoch of the reference objects in MJD.
941 refObjectLoader : instance of
942 `lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader`
943 Reference object loader instance.
944 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
945 Parameters to use for the device part of the model.
946 colorCatalog : `lsst.afw.table.SimpleCatalog`
947 Catalog containing object coordinates and magnitudes.
948 inputCamera : `lsst.afw.cameraGeom.Camera`, optional
949 Camera to be used as template when constructing new camera.
950 nCores : `int`, optional
951 Number of cores to use during WCS fit.
953 Returns
954 -------
955 result : `lsst.pipe.base.Struct`
956 ``outputWcss`` : `list` [`lsst.afw.table.ExposureCatalog`]
957 List of exposure catalogs (one per visit) with the WCS for each
958 detector set by the new fitted WCS.
959 ``fitModel`` : `wcsfit.WCSFit`
960 Model-fitting object with final model parameters.
961 ``outputCatalog`` : `pyarrow.Table`
962 Catalog with fit residuals of all sources used.
963 ``starCatalog`` : `pyarrow.Table`
964 Catalog with best-fit positions of the objects fit.
965 ``modelParams`` : `dict`
966 Parameters and covariance of the best-fit WCS model.
967 ``cameraModelParams`` : `dict` [`str`, `np.ndarray`]
968 Parameters of the device part of the model, in the format
969 needed as input for future runs.
970 ``colorParams`` : `dict` [`int`, `np.ndarray`]
971 DCR parameters fit in RA and Dec directions for each visit.
972 ``camera`` : `lsst.afw.cameraGeom.Camera`
973 Camera object constructed from the per-detector model.
974 """
975 self.log.info("Gather instrument, exposure, and field info")
977 # Get RA, Dec, MJD, etc., for the input visits
978 exposureInfo, exposuresHelper, extensionInfo, instruments = self._get_exposure_info(
979 inputVisitSummaries
980 )
982 # Get information about the extent of the input visits
983 fields, fieldCenter, fieldRadius = self._prep_sky(inputVisitSummaries, exposureInfo.medianEpoch)
984 self.log.info("Field center set at %s with radius %s degrees", fieldCenter, fieldRadius.asDegrees())
986 self.log.info("Load catalogs and associate sources")
987 if not self.config.useIsolatedStars:
988 # Set up class to associate sources into matches using a
989 # friends-of-friends algorithm
990 associations = wcsfit.FoFClass(
991 fields,
992 instruments,
993 exposuresHelper,
994 [fieldRadius.asDegrees()],
995 (self.config.matchRadius * u.arcsec).to(u.degree).value,
996 )
997 else:
998 associations = None
1000 # Add the reference catalog to the associator
1001 medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format="jyear").mjd
1002 refObjects, refCovariance = self._load_refcat(
1003 refObjectLoader,
1004 extensionInfo,
1005 epoch=medianEpoch,
1006 center=fieldCenter,
1007 radius=fieldRadius,
1008 associations=associations,
1009 )
1011 if self.config.useIsolatedStars:
1012 allRefObjects = {0: refObjects}
1013 associations, sourceDict = self._associate_from_isolated_sources(
1014 isolatedStarSources, isolatedStarCatalogs, extensionInfo, allRefObjects
1015 )
1016 else:
1017 # Add the science catalogs and associate new sources as they are
1018 # added.
1019 sourceIndices, usedColumns = self._load_catalogs_and_associate(
1020 associations, inputCatalogRefs, extensionInfo
1021 )
1022 self._check_degeneracies(associations, extensionInfo)
1024 self.log.info("Fit the WCSs")
1025 # Set up a YAML-type string using the config variables and a sample
1026 # visit
1027 inputYaml, mapTemplate = self.make_yaml(
1028 inputVisitSummaries[0],
1029 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
1030 )
1032 # Set the verbosity level for WCSFit from the task log level.
1033 # TODO: DM-36850, Add lsst.log to gbdes so that log messages are
1034 # properly propagated.
1035 loglevel = self.log.getEffectiveLevel()
1036 if loglevel >= self.log.WARNING: 1036 ↛ 1037line 1036 didn't jump to line 1037 because the condition on line 1036 was never true
1037 verbose = 0
1038 elif loglevel == self.log.INFO: 1038 ↛ 1039line 1038 didn't jump to line 1039 because the condition on line 1038 was never true
1039 verbose = 1
1040 else:
1041 verbose = 2
1043 # Set up the WCS-fitting class using the results of the FOF associator
1044 if self.config.useInputCameraModel:
1045 fixMaps = ",".join(inputCameraModel.keys())
1046 else:
1047 fixMaps = ""
1048 wcsf = wcsfit.WCSFit(
1049 fields,
1050 instruments,
1051 exposuresHelper,
1052 extensionInfo.visitIndex,
1053 extensionInfo.detectorIndex,
1054 inputYaml,
1055 extensionInfo.wcs,
1056 associations.sequence,
1057 associations.extn,
1058 associations.obj,
1059 sysErr=self.config.systematicError,
1060 refSysErr=self.config.referenceSystematicError,
1061 usePM=self.config.fitProperMotion,
1062 verbose=verbose,
1063 fixMaps=fixMaps,
1064 num_threads=nCores,
1065 )
1066 # Add the science and reference sources
1067 if self.config.useIsolatedStars:
1068 self._add_objects_from_dict(wcsf, sourceDict, extensionInfo)
1069 else:
1070 self._add_objects(wcsf, inputCatalogRefs, sourceIndices, extensionInfo, usedColumns)
1071 self._add_ref_objects(wcsf, refObjects, refCovariance, extensionInfo)
1072 if self.config.useColor:
1073 self._add_color_objects(wcsf, colorCatalog)
1075 # There must be at least as many sources per visit as the number of
1076 # free parameters in the per-visit mapping. Set minFitExposures to be
1077 # the number of free parameters divided by the fraction of non-reserved
1078 # sources, so that visits with fewer sources are dropped.
1079 nCoeffVisitModel = _nCoeffsFromDegree(self.config.exposurePolyOrder)
1080 minFitExposures = int(np.ceil(nCoeffVisitModel / (1 - self.config.fitReserveFraction)))
1081 # Do the WCS fit
1082 try:
1083 wcsf.fit(
1084 reserveFraction=self.config.fitReserveFraction,
1085 randomNumberSeed=self.config.fitReserveRandomSeed,
1086 minFitExposures=minFitExposures,
1087 clipThresh=self.config.clipThresh,
1088 clipFraction=self.config.clipFraction,
1089 )
1090 except RuntimeError as e:
1091 if "Cholesky decomposition failed" in str(e):
1092 raise CholeskyError() from e
1093 else:
1094 raise
1096 self.log.info("WCS fitting done")
1098 outputWcss, cameraParams, colorParams, camera, partialOutputs = self._make_outputs(
1099 wcsf,
1100 inputVisitSummaries,
1101 exposureInfo,
1102 mapTemplate,
1103 extensionInfo,
1104 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
1105 inputCamera=(inputCamera if self.config.buildCamera else None),
1106 )
1107 outputCatalog = wcsf.getOutputCatalog()
1108 outputCatalog["exposureName"] = np.array(outputCatalog["exposureName"])
1109 outputCatalog["deviceName"] = np.array(outputCatalog["deviceName"])
1111 starCatalog = wcsf.getStarCatalog()
1112 modelParams = self._compute_model_params(wcsf) if self.config.saveModelParams else None
1114 return pipeBase.Struct(
1115 outputWcss=outputWcss,
1116 fitModel=wcsf,
1117 outputCatalog=outputCatalog,
1118 starCatalog=starCatalog,
1119 modelParams=modelParams,
1120 cameraModelParams=cameraParams,
1121 colorParams=colorParams,
1122 camera=camera,
1123 partialOutputs=partialOutputs,
1124 )
1126 def _prep_sky(self, inputVisitSummaries, epoch, fieldName="Field"):
1127 """Get center and radius of the input tract. This assumes that all
1128 visits will be put into the same `wcsfit.Field` and fit together.
1130 Paramaters
1131 ----------
1132 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
1133 List of catalogs with per-detector summary information.
1134 epoch : float
1135 Reference epoch.
1136 fieldName : str
1137 Name of the field, used internally.
1139 Returns
1140 -------
1141 fields : `wcsfit.Fields`
1142 Object with field information.
1143 center : `lsst.geom.SpherePoint`
1144 Center of the field.
1145 radius : `lsst.sphgeom._sphgeom.Angle`
1146 Radius of the bounding circle of the tract.
1147 """
1148 allDetectorCorners = []
1149 for visSum in inputVisitSummaries:
1150 detectorCorners = [
1151 lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees).getVector()
1152 for (ra, dec) in zip(visSum["raCorners"].ravel(), visSum["decCorners"].ravel())
1153 if (np.isfinite(ra) and (np.isfinite(dec)))
1154 ]
1155 allDetectorCorners.extend(detectorCorners)
1156 boundingCircle = lsst.sphgeom.ConvexPolygon.convexHull(allDetectorCorners).getBoundingCircle()
1157 center = lsst.geom.SpherePoint(boundingCircle.getCenter())
1158 ra = center.getRa().asDegrees()
1159 dec = center.getDec().asDegrees()
1160 radius = boundingCircle.getOpeningAngle()
1162 # wcsfit.Fields describes a list of fields, but we assume all
1163 # observations will be fit together in one field.
1164 fields = wcsfit.Fields([fieldName], [ra], [dec], [epoch])
1166 return fields, center, radius
1168 def _get_exposure_info(
1169 self,
1170 inputVisitSummaries,
1171 fieldNumber=0,
1172 refEpoch=None,
1173 fieldRegions=None,
1174 ):
1175 """Get various information about the input visits to feed to the
1176 fitting routines.
1178 Parameters
1179 ----------
1180 inputVisitSummaries : `list [`lsst.afw.table.ExposureCatalog`]
1181 Tables for each visit with information for detectors.
1182 fieldNumber : `int`, optional
1183 Index of the field for these visits. Should be zero if all data is
1184 being fit together. This is ignored if `fieldRegions` is not None.
1185 refEpoch : `float`, optional
1186 Epoch of the reference objects in MJD.
1187 fieldRegions : `dict` [`int`, `lsst.sphgeom.ConvexPolygon`], optional
1188 Dictionary of regions encompassing each group of input visits
1189 keyed by an arbitrary index.
1191 Returns
1192 -------
1193 exposureInfo : `lsst.pipe.base.Struct`
1194 Struct containing general properties for the visits:
1195 ``visits`` : `list`
1196 List of visit names.
1197 ``detectors`` : `list`
1198 List of all detectors in any visit.
1199 ``ras`` : `list` [`float`]
1200 List of boresight RAs for each visit.
1201 ``decs`` : `list` [`float`]
1202 List of borseight Decs for each visit.
1203 ``medianEpoch`` : float
1204 Median epoch of all visits in decimal-year format.
1205 exposuresHelper : `wcsfit.ExposuresHelper`
1206 Object containing information about the input visits.
1207 extensionInfo : `lsst.pipe.base.Struct`
1208 Struct containing properties for each extension (visit/detector):
1209 ``visit`` : `np.ndarray`
1210 Name of the visit for this extension.
1211 ``detector`` : `np.ndarray`
1212 Name of the detector for this extension.
1213 ``visitIndex` : `np.ndarray` [`int`]
1214 Index of visit for this extension.
1215 ``detectorIndex`` : `np.ndarray` [`int`]
1216 Index of the detector for this extension.
1217 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1218 Initial WCS for this extension.
1219 ``extensionType`` : `np.ndarray` [`str`]
1220 "SCIENCE" or "REFERENCE".
1221 instruments : `list` [`wcsfit.Instrument`]
1222 List of instrument objects used for the input visits.
1223 """
1224 exposureNames = []
1225 ras = []
1226 decs = []
1227 visits = []
1228 detectors = []
1229 airmasses = []
1230 exposureTimes = []
1231 mjds = []
1232 observatories = []
1233 wcss = []
1234 fieldNumbers = []
1236 extensionType = []
1237 extensionVisitIndices = []
1238 extensionDetectorIndices = []
1239 extensionVisits = []
1240 extensionDetectors = []
1242 instruments, instrumentNumbers = _get_instruments(inputVisitSummaries)
1244 # Get information for all the science visits
1245 for v, visitSummary in enumerate(inputVisitSummaries):
1246 visitInfo = visitSummary[0].getVisitInfo()
1247 visit = visitSummary[0]["visit"]
1248 visits.append(visit)
1249 exposureNames.append(str(visit))
1250 raDec = visitInfo.getBoresightRaDec()
1251 ras.append(raDec.getRa().asRadians())
1252 decs.append(raDec.getDec().asRadians())
1253 if fieldRegions is not None:
1254 inField = [r for r, region in fieldRegions.items() if region.contains(raDec.getVector())]
1255 if len(inField) != 1: 1255 ↛ 1256line 1255 didn't jump to line 1256 because the condition on line 1255 was never true
1256 raise RuntimeError(
1257 f"Visit should be in one and only one field, but {visit} is contained "
1258 f"in {len(inField)} fields."
1259 )
1260 fieldNumbers.append(inField[0])
1261 else:
1262 fieldNumbers.append(fieldNumber)
1263 airmasses.append(visitInfo.getBoresightAirmass())
1264 exposureTimes.append(visitInfo.getExposureTime())
1265 obsDate = visitInfo.getDate()
1266 obsMJD = obsDate.get(obsDate.MJD)
1267 mjds.append(obsMJD)
1268 # Get the observatory ICRS position for use in fitting parallax
1269 obsLon = visitInfo.observatory.getLongitude().asDegrees()
1270 obsLat = visitInfo.observatory.getLatitude().asDegrees()
1271 obsElev = visitInfo.observatory.getElevation()
1272 earthLocation = astropy.coordinates.EarthLocation.from_geodetic(obsLon, obsLat, obsElev)
1273 observatory_gcrs = earthLocation.get_gcrs(astropy.time.Time(obsMJD, format="mjd"))
1274 observatory_icrs = observatory_gcrs.transform_to(astropy.coordinates.ICRS())
1275 # We want the position in AU in Cartesian coordinates
1276 observatories.append(observatory_icrs.cartesian.xyz.to(u.AU).value)
1278 validDetectors = [row for row in visitSummary if row.wcs is not None]
1279 nDetectorFraction = len(validDetectors) / len(visitSummary)
1280 if nDetectorFraction < self.config.minDetectorFraction: 1280 ↛ 1281line 1280 didn't jump to line 1281 because the condition on line 1280 was never true
1281 self.log.warning(
1282 "Visit %d has only %d detectors with valid WCSs (%s of total) and will be dropped.",
1283 visit,
1284 len(validDetectors),
1285 nDetectorFraction,
1286 )
1287 continue
1289 for row in visitSummary:
1290 detector = row["id"]
1292 wcs = row.getWcs()
1293 if wcs is None:
1294 self.log.warning(
1295 "WCS is None for visit %d, detector %d: this extension (visit/detector) will be "
1296 "dropped.",
1297 visit,
1298 detector,
1299 )
1300 continue
1301 else:
1302 wcsRA = wcs.getSkyOrigin().getRa().asRadians()
1303 wcsDec = wcs.getSkyOrigin().getDec().asRadians()
1304 tangentPoint = wcsfit.Gnomonic(wcsRA, wcsDec)
1305 mapping = wcs.getFrameDict().getMapping("PIXELS", "IWC")
1306 gbdes_wcs = wcsfit.Wcs(wcsfit.ASTMap(mapping), tangentPoint)
1307 wcss.append(gbdes_wcs)
1309 if detector not in detectors:
1310 detectors.append(detector)
1311 if not instruments[instrumentNumbers[v]].hasDevice(str(detector)):
1312 detectorBounds = wcsfit.Bounds(
1313 row["bbox_min_x"], row["bbox_max_x"], row["bbox_min_y"], row["bbox_max_y"]
1314 )
1315 instruments[instrumentNumbers[v]].addDevice(str(detector), detectorBounds)
1317 detectorIndex = np.flatnonzero(detector == np.array(detectors))[0]
1318 extensionVisitIndices.append(v)
1319 extensionDetectorIndices.append(detectorIndex)
1320 extensionVisits.append(visit)
1321 extensionDetectors.append(detector)
1322 extensionType.append("SCIENCE")
1324 if len(wcss) == 0: 1324 ↛ 1325line 1324 didn't jump to line 1325 because the condition on line 1324 was never true
1325 raise pipeBase.NoWorkFound("No input extensions have valid WCSs.")
1327 # Set the reference epoch to be the median of the science visits.
1328 # The reference catalog will be shifted to this date.
1329 if self.config.setRefEpoch is None: 1329 ↛ 1333line 1329 didn't jump to line 1333 because the condition on line 1329 was always true
1330 medianMJD = np.median(mjds)
1331 self.log.info(f"Ref epoch set to median: {medianMJD}")
1332 else:
1333 medianMJD = self.config.setRefEpoch
1334 self.log.info(f"Ref epoch set by user: {medianMJD}")
1335 medianEpoch = astropy.time.Time(medianMJD, format="mjd").jyear
1337 # Add information for the reference catalog. Most of the values are
1338 # not used. There needs to be a separate catalog for each field.
1339 if fieldRegions is None:
1340 fieldRegions = {0: None}
1341 for f in fieldRegions:
1342 exposureNames.append("REFERENCE")
1343 # Make the "visit" number the field * -1 to disambiguate it from
1344 # any potential visit number:
1345 visits.append(-1 * f)
1346 fieldNumbers.append(f)
1347 if self.config.fitProperMotion: 1347 ↛ 1348line 1347 didn't jump to line 1348 because the condition on line 1347 was never true
1348 instrumentNumbers.append(-2)
1349 else:
1350 instrumentNumbers.append(-1)
1351 ras.append(0.0)
1352 decs.append(0.0)
1353 airmasses.append(0.0)
1354 exposureTimes.append(0)
1355 mjds.append((refEpoch if (refEpoch is not None) else medianMJD))
1356 observatories.append(np.array([0, 0, 0]))
1357 identity = wcsfit.IdentityMap()
1358 icrs = wcsfit.SphericalICRS()
1359 refWcs = wcsfit.Wcs(identity, icrs, "Identity", np.pi / 180.0)
1360 wcss.append(refWcs)
1362 extensionVisitIndices.append(len(exposureNames) - 1)
1363 extensionDetectorIndices.append(-1) # REFERENCE device must be -1
1364 extensionVisits.append(-1 * f)
1365 extensionDetectors.append(-1)
1366 extensionType.append("REFERENCE")
1368 # Make a table of information to use elsewhere in the class
1369 extensionInfo = pipeBase.Struct(
1370 visit=np.array(extensionVisits),
1371 detector=np.array(extensionDetectors),
1372 visitIndex=np.array(extensionVisitIndices),
1373 detectorIndex=np.array(extensionDetectorIndices),
1374 wcs=np.array(wcss),
1375 extensionType=np.array(extensionType),
1376 )
1378 # Make the exposureHelper object to use in the fitting routines
1379 exposuresHelper = wcsfit.ExposuresHelper(
1380 exposureNames,
1381 fieldNumbers,
1382 instrumentNumbers,
1383 ras,
1384 decs,
1385 airmasses,
1386 exposureTimes,
1387 mjds,
1388 observatories,
1389 )
1391 exposureInfo = pipeBase.Struct(
1392 visits=visits, detectors=detectors, ras=ras, decs=decs, medianEpoch=medianEpoch
1393 )
1395 return exposureInfo, exposuresHelper, extensionInfo, instruments
1397 def _load_refcat(
1398 self,
1399 refObjectLoader,
1400 extensionInfo,
1401 epoch=None,
1402 fieldIndex=0,
1403 associations=None,
1404 center=None,
1405 radius=None,
1406 region=None,
1407 ):
1408 """Load the reference catalog and add reference objects to the
1409 `wcsfit.FoFClass` object.
1411 Parameters
1412 ----------
1413 refObjectLoader :
1414 `lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader`
1415 Object set up to load reference catalog objects.
1416 extensionInfo : `lsst.pipe.base.Struct`
1417 Struct containing properties for each extension (visit/detector).
1418 ``visit`` : `np.ndarray`
1419 Name of the visit for this extension.
1420 ``detector`` : `np.ndarray`
1421 Name of the detector for this extension.
1422 ``visitIndex` : `np.ndarray` [`int`]
1423 Index of visit for this extension.
1424 ``detectorIndex`` : `np.ndarray` [`int`]
1425 Index of the detector for this extension.
1426 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1427 Initial WCS for this extension.
1428 ``extensionType`` : `np.ndarray` [`str`]
1429 "SCIENCE" or "REFERENCE".
1430 epoch : `float`, optional
1431 MJD to which to correct the object positions.
1432 fieldIndex : `int`, optional
1433 Index of the field. Should be zero if all the data is fit together.
1434 associations : `wcsfit.FoFClass`, optional
1435 Object to which to add the catalog of reference objects.
1436 center : `lsst.geom.SpherePoint`, optional
1437 Center of the circle in which to load reference objects. Ignored if
1438 `region` is set. If used, `radius` must also be set.
1439 radius : `lsst.sphgeom._sphgeom.Angle`, optional
1440 Radius of the circle in which to load reference objects. Ignored if
1441 `region` is set. If used, `center` must also be set.
1442 region : `lsst.sphgeom.ConvexPolygon`, optional
1443 Region in which to load reference objects.
1445 Returns
1446 -------
1447 refObjects : `dict`
1448 Position and error information of reference objects.
1449 refCovariance : `list` [`float`]
1450 Flattened output covariance matrix.
1451 """
1452 if self.config.applyRefCatProperMotion: 1452 ↛ 1453line 1452 didn't jump to line 1453 because the condition on line 1452 was never true
1453 formattedEpoch = astropy.time.Time(epoch, format="mjd")
1454 else:
1455 formattedEpoch = None
1457 neededColumns = ["coord_ra", "coord_dec", "coord_raErr", "coord_decErr"]
1458 if self.config.applyRefCatProperMotion: 1458 ↛ 1459line 1458 didn't jump to line 1459 because the condition on line 1458 was never true
1459 neededColumns += [
1460 "pm_ra",
1461 "pm_dec",
1462 "parallax",
1463 "pm_raErr",
1464 "pm_decErr",
1465 "parallaxErr",
1466 ]
1467 # Get refcat version from refcat metadata
1468 refCatMetadata = refObjectLoader.refCats[0].get().getMetadata()
1469 # DM-47181: Added this to work for LSSTComCam with
1470 # the_monster_20240904 catalog that does not have this key.
1471 if "REFCAT_FORMAT_VERSION" not in refCatMetadata: 1471 ↛ 1472line 1471 didn't jump to line 1472 because the condition on line 1471 was never true
1472 refCatVersion = 2
1473 else:
1474 refCatVersion = refCatMetadata["REFCAT_FORMAT_VERSION"]
1475 if refCatVersion == 2: 1475 ↛ 1476line 1475 didn't jump to line 1476 because the condition on line 1475 was never true
1476 neededColumns += [
1477 "coord_ra_coord_dec_Cov",
1478 "coord_ra_pm_ra_Cov",
1479 "coord_ra_pm_dec_Cov",
1480 "coord_ra_parallax_Cov",
1481 "coord_dec_pm_ra_Cov",
1482 "coord_dec_pm_dec_Cov",
1483 "coord_dec_parallax_Cov",
1484 "pm_ra_pm_dec_Cov",
1485 "pm_ra_parallax_Cov",
1486 "pm_dec_parallax_Cov",
1487 ]
1489 # Load each shard of the reference catalog separately to avoid a spike
1490 # in the memory load.
1491 refCatShards = []
1492 for dataId, cat in zip(refObjectLoader.dataIds, refObjectLoader.refCats):
1493 miniRefObjectLoader = ReferenceObjectLoader(
1494 dataIds=[dataId],
1495 refCats=[cat],
1496 config=refObjectLoader.config,
1497 log=self.log,
1498 )
1499 try:
1500 if region is not None:
1501 skyRegion = miniRefObjectLoader.loadRegion(
1502 region, self.config.referenceFilter, epoch=formattedEpoch
1503 )
1504 elif (center is not None) and (radius is not None): 1504 ↛ 1509line 1504 didn't jump to line 1509 because the condition on line 1504 was always true
1505 skyRegion = miniRefObjectLoader.loadSkyCircle(
1506 center, radius, self.config.referenceFilter, epoch=formattedEpoch
1507 )
1508 else:
1509 raise RuntimeError("Either `region` or `center` and `radius` must be set.")
1510 except RuntimeError:
1511 self.log.debug("Reference catalog shard has no objects inside the region.")
1512 continue
1513 selected = self.referenceSelector.run(skyRegion.refCat)
1514 # Need memory contiguity to get reference filters as a vector.
1515 if not selected.sourceCat.isContiguous(): 1515 ↛ 1516line 1515 didn't jump to line 1516 because the condition on line 1515 was never true
1516 refCatShard = selected.sourceCat.copy(deep=True)
1517 else:
1518 refCatShard = selected.sourceCat
1519 refCatShard = refCatShard.asAstropy()[neededColumns]
1521 # In Gaia DR3, missing values are denoted by NaNs.
1522 finiteInd = np.isfinite(refCatShard["coord_ra"]) & np.isfinite(refCatShard["coord_dec"])
1523 refCatShard = refCatShard[finiteInd]
1524 refCatShards.append(refCatShard)
1525 refCat = vstack(refCatShards)
1527 if self.config.excludeNonPMObjects and self.config.applyRefCatProperMotion: 1527 ↛ 1529line 1527 didn't jump to line 1529 because the condition on line 1527 was never true
1528 # Gaia DR2 has zeros for missing data, while Gaia DR3 has NaNs:
1529 hasPM = (
1530 (refCat["pm_raErr"] != 0) & np.isfinite(refCat["pm_raErr"]) & np.isfinite(refCat["pm_decErr"])
1531 )
1532 refCat = refCat[hasPM]
1534 ra = (refCat["coord_ra"]).to(u.degree).to_value().tolist()
1535 dec = (refCat["coord_dec"]).to(u.degree).to_value().tolist()
1536 raCov = ((refCat["coord_raErr"]).to(u.degree).to_value() ** 2).tolist()
1537 decCov = ((refCat["coord_decErr"]).to(u.degree).to_value() ** 2).tolist()
1539 if refCatVersion == 2: 1539 ↛ 1540line 1539 didn't jump to line 1540 because the condition on line 1539 was never true
1540 raDecCov = (refCat["coord_ra_coord_dec_Cov"]).to(u.degree**2).to_value().tolist()
1541 else:
1542 raDecCov = np.zeros(len(ra))
1544 refObjects = {"ra": ra, "dec": dec, "raCov": raCov, "decCov": decCov, "raDecCov": raDecCov}
1545 refCovariance = []
1547 if self.config.fitProperMotion: 1547 ↛ 1548line 1547 didn't jump to line 1548 because the condition on line 1547 was never true
1548 raPM = (refCat["pm_ra"]).to(u.marcsec).to_value().tolist()
1549 decPM = (refCat["pm_dec"]).to(u.marcsec).to_value().tolist()
1550 parallax = (refCat["parallax"]).to(u.marcsec).to_value().tolist()
1551 cov = _make_ref_covariance_matrix(refCat, version=refCatVersion)
1552 pmDict = {"raPM": raPM, "decPM": decPM, "parallax": parallax}
1553 refObjects.update(pmDict)
1554 refCovariance = cov
1556 if associations is not None:
1557 extensionIndex = np.flatnonzero(extensionInfo.extensionType == "REFERENCE")[0]
1558 visitIndex = extensionInfo.visitIndex[extensionIndex]
1559 detectorIndex = extensionInfo.detectorIndex[extensionIndex]
1560 instrumentIndex = -1 # -1 indicates the reference catalog
1561 refWcs = extensionInfo.wcs[extensionIndex]
1563 associations.addCatalog(
1564 refWcs,
1565 "STELLAR",
1566 visitIndex,
1567 fieldIndex,
1568 instrumentIndex,
1569 detectorIndex,
1570 extensionIndex,
1571 np.ones(len(refCat), dtype=bool),
1572 ra,
1573 dec,
1574 np.arange(len(ra)),
1575 )
1577 return refObjects, refCovariance
1579 @staticmethod
1580 def _find_extension_index(extensionInfo, visit, detector):
1581 """Find the index for a given extension from its visit and detector
1582 number.
1584 If no match is found, None is returned.
1586 Parameters
1587 ----------
1588 extensionInfo : `lsst.pipe.base.Struct`
1589 Struct containing properties for each extension.
1590 visit : `int`
1591 Visit number
1592 detector : `int`
1593 Detector number
1595 Returns
1596 -------
1597 extensionIndex : `int` or None
1598 Index of this extension
1599 """
1600 findExtension = np.flatnonzero((extensionInfo.visit == visit) & (extensionInfo.detector == detector))
1601 if len(findExtension) == 0:
1602 extensionIndex = None
1603 else:
1604 extensionIndex = findExtension[0]
1605 return extensionIndex
1607 def _load_catalogs_and_associate(
1608 self, associations, inputCatalogRefs, extensionInfo, fieldIndex=0, instrumentIndex=0
1609 ):
1610 """Load the science catalogs and add the sources to the associator
1611 class `wcsfit.FoFClass`, associating them into matches as you go.
1613 Parameters
1614 ----------
1615 associations : `wcsfit.FoFClass`
1616 Object to which to add the catalog of source and which performs
1617 the source association.
1618 inputCatalogRefs : `list`
1619 List of DeferredDatasetHandles pointing to visit-level source
1620 tables.
1621 extensionInfo : `lsst.pipe.base.Struct`
1622 Struct containing properties for each extension (visit/detector).
1623 ``visit`` : `np.ndarray`
1624 Name of the visit for this extension.
1625 ``detector`` : `np.ndarray`
1626 Name of the detector for this extension.
1627 ``visitIndex` : `np.ndarray` [`int`]
1628 Index of visit for this extension.
1629 ``detectorIndex`` : `np.ndarray` [`int`]
1630 Index of the detector for this extension.
1631 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1632 Initial WCS for this extension.
1633 ``extensionType`` : `np.ndarray` [`str`]
1634 "SCIENCE" or "REFERENCE".
1635 fieldIndex : `int`
1636 Index of the field for these catalogs. Should be zero assuming all
1637 data is being fit together.
1638 instrumentIndex : `int`
1639 Index of the instrument for these catalogs. Should be zero
1640 assuming all data comes from the same instrument.
1642 Returns
1643 -------
1644 sourceIndices : `list`
1645 List of boolean arrays used to select sources.
1646 columns : `list` [`str`]
1647 List of columns needed from source tables.
1648 """
1649 columns = [
1650 "detector",
1651 "sourceId",
1652 "x",
1653 "xErr",
1654 "y",
1655 "yErr",
1656 "ixx",
1657 "iyy",
1658 "ixy",
1659 f"{self.config.sourceFluxType}_instFlux",
1660 f"{self.config.sourceFluxType}_instFluxErr",
1661 ]
1662 if self.sourceSelector.config.doFlags: 1662 ↛ 1664line 1662 didn't jump to line 1664 because the condition on line 1662 was always true
1663 columns.extend(self.sourceSelector.config.flags.bad)
1664 if self.sourceSelector.config.doUnresolved: 1664 ↛ 1666line 1664 didn't jump to line 1666 because the condition on line 1664 was always true
1665 columns.append(self.sourceSelector.config.unresolved.name)
1666 if self.sourceSelector.config.doIsolated: 1666 ↛ 1669line 1666 didn't jump to line 1669 because the condition on line 1666 was always true
1667 columns.append(self.sourceSelector.config.isolated.parentName)
1668 columns.append(self.sourceSelector.config.isolated.nChildName)
1669 if self.sourceSelector.config.doRequirePrimary: 1669 ↛ 1672line 1669 didn't jump to line 1672 because the condition on line 1669 was always true
1670 columns.append(self.sourceSelector.config.requirePrimary.primaryColName)
1672 sourceIndices = [None] * len(extensionInfo.visit)
1673 for inputCatalogRef in inputCatalogRefs:
1674 visit = inputCatalogRef.dataId["visit"]
1675 inputCatalog = inputCatalogRef.get(parameters={"columns": columns})
1676 # Get a sorted array of detector names
1677 detectors = np.unique(inputCatalog["detector"])
1679 for detector in detectors:
1680 detectorSources = inputCatalog[inputCatalog["detector"] == detector]
1681 xCov = detectorSources["xErr"] ** 2
1682 yCov = detectorSources["yErr"] ** 2
1683 xyCov = (
1684 detectorSources["ixy"] * (xCov + yCov) / (detectorSources["ixx"] + detectorSources["iyy"])
1685 )
1686 # Remove sources with bad shape measurements
1687 goodShapes = xyCov**2 <= (xCov * yCov)
1688 selected = self.sourceSelector.run(detectorSources)
1689 goodInds = selected.selected & goodShapes
1690 isStar = np.ones(goodInds.sum())
1691 extensionIndex = self._find_extension_index(extensionInfo, visit, detector)
1692 if extensionIndex is None:
1693 # This extension does not have information necessary for
1694 # fit. Skip the detections from this detector for this
1695 # visit.
1696 continue
1697 detectorIndex = extensionInfo.detectorIndex[extensionIndex]
1698 visitIndex = extensionInfo.visitIndex[extensionIndex]
1700 sourceIndices[extensionIndex] = goodInds
1702 wcs = extensionInfo.wcs[extensionIndex]
1703 associations.reprojectWCS(wcs, fieldIndex)
1705 associations.addCatalog(
1706 wcs,
1707 "STELLAR",
1708 visitIndex,
1709 fieldIndex,
1710 instrumentIndex,
1711 detectorIndex,
1712 extensionIndex,
1713 isStar,
1714 detectorSources[goodInds]["x"].to_list(),
1715 detectorSources[goodInds]["y"].to_list(),
1716 np.arange(goodInds.sum()),
1717 )
1719 associations.sortMatches(
1720 fieldIndex, minMatches=self.config.minMatches, allowSelfMatches=self.config.allowSelfMatches
1721 )
1723 return sourceIndices, columns
1725 def _associate_from_isolated_sources(
1726 self, isolatedStarSourceRefs, isolatedStarCatalogRefs, extensionInfo, refObjects
1727 ):
1728 """Match the input catalog of isolated stars with the reference catalog
1729 and transform the combined isolated star sources and reference source
1730 into the format needed for gbdes.
1732 Parameters
1733 ----------
1734 isolatedStarSourceRefs : `list` [`DeferredDatasetHandle`]
1735 List of handles pointing to isolated star sources.
1736 isolatedStarCatalogRefs: `list` [`DeferredDatasetHandle`]
1737 List of handles pointing to isolated star catalogs.
1738 extensionInfo : `lsst.pipe.base.Struct`
1739 Struct containing properties for each extension (visit/detector).
1740 ``visit`` : `np.ndarray`
1741 Name of the visit for this extension.
1742 ``detector`` : `np.ndarray`
1743 Name of the detector for this extension.
1744 ``visitIndex` : `np.ndarray` [`int`]
1745 Index of visit for this extension.
1746 ``detectorIndex`` : `np.ndarray` [`int`]
1747 Index of the detector for this extension.
1748 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1749 Initial WCS for this extension.
1750 ``extensionType`` : `np.ndarray` [`str`]
1751 "SCIENCE" or "REFERENCE".
1752 refObjects : `dict`
1753 Dictionary of dictionaries containing the position and error
1754 information of reference objects.
1756 Returns
1757 -------
1758 associations : `lsst.pipe.base.Struct`
1759 Struct containing the associations of sources with objects.
1760 sourceDict : `dict` [`int`, [`int`, [`str`, `list` [`float`]]]]
1761 Dictionary containing the source centroids for each visit.
1762 """
1763 sequences = []
1764 extensions = []
1765 object_indices = []
1767 sourceColumns = ["x", "y", "xErr", "yErr", "ixx", "ixy", "iyy", "obj_index", "visit", "detector"]
1768 catalogColumns = ["ra", "dec"]
1770 sortedVisits = np.unique(extensionInfo.visit)
1772 sourceDict = dict([(visit, {}) for visit in np.unique(extensionInfo.visit)])
1773 for visit, detector in zip(extensionInfo.visit, extensionInfo.detector):
1774 sourceDict[visit][detector] = {"x": [], "y": [], "xCov": [], "yCov": [], "xyCov": []}
1776 for isolatedStarCatalogRef, isolatedStarSourceRef in zip(
1777 isolatedStarCatalogRefs, isolatedStarSourceRefs
1778 ):
1779 isolatedStarCatalog = isolatedStarCatalogRef.get(parameters={"columns": catalogColumns})
1780 isolatedStarSources = isolatedStarSourceRef.get(parameters={"columns": sourceColumns})
1781 if len(isolatedStarCatalog) == 0: 1781 ↛ 1784line 1781 didn't jump to line 1784 because the condition on line 1781 was never true
1782 # This is expected when only one visit overlaps with a given
1783 # tract, meaning that no sources can be associated.
1784 self.log.debug(
1785 "Skipping tract %d, which has no associated isolated stars",
1786 isolatedStarCatalogRef.dataId["tract"],
1787 )
1788 continue
1790 # Only use isolated sources that are from the visits in this fit.
1791 sub1 = np.clip(
1792 np.searchsorted(sortedVisits, isolatedStarSources["visit"]),
1793 0,
1794 len(sortedVisits) - 1,
1795 )
1796 visitsInFit = sortedVisits[sub1] == isolatedStarSources["visit"].to_numpy()
1797 if not visitsInFit.any(): 1797 ↛ 1798line 1797 didn't jump to line 1798 because the condition on line 1797 was never true
1798 continue
1799 isolatedStarSources = isolatedStarSources[visitsInFit]
1801 # Match the reference stars to the existing isolated stars, then
1802 # insert the reference stars into the isolated star sources.
1803 allVisits = np.copy(isolatedStarSources["visit"])
1804 allDetectors = np.copy(isolatedStarSources["detector"])
1805 allObjectIndices = np.copy(isolatedStarSources["obj_index"])
1806 issIndices = np.copy(isolatedStarSources.index)
1807 for f, regionRefObjects in refObjects.items():
1808 # Use the same matching technique that is done in
1809 # isolatedStarAssociation and fgcmBuildFromIsolatedStars.
1810 with Matcher(
1811 isolatedStarCatalog["ra"].to_numpy(), isolatedStarCatalog["dec"].to_numpy()
1812 ) as matcher:
1813 idx, idx_isoStarCat, idx_refObjects, d = matcher.query_radius(
1814 np.array(regionRefObjects["ra"]),
1815 np.array(regionRefObjects["dec"]),
1816 self.config.matchRadius / 3600.0,
1817 return_indices=True,
1818 )
1820 # Remove sources that were matched to multiple reference
1821 # objects.
1822 if len(idx_isoStarCat) != len(np.unique(idx_isoStarCat)): 1822 ↛ 1823line 1822 didn't jump to line 1823 because the condition on line 1822 was never true
1823 _, matchInverse, matchCount = np.unique(
1824 idx_isoStarCat, return_inverse=True, return_counts=True
1825 )
1826 idx_isoStarCat = idx_isoStarCat[matchCount[matchInverse] == 1]
1827 idx_refObjects = idx_refObjects[matchCount[matchInverse] == 1]
1829 refSort = np.searchsorted(isolatedStarSources["obj_index"], idx_isoStarCat)
1830 sub1 = np.clip(refSort, 0, len(isolatedStarSources) - 1)
1831 goodMatches = isolatedStarSources["obj_index"].iloc[sub1].to_numpy() == idx_isoStarCat
1833 refDetector = np.ones(len(idx_isoStarCat)) * -1
1834 # The "visit" for the reference catalogs is the field times -1.
1835 refVisit = np.ones(len(idx_isoStarCat)) * f * -1
1837 allVisits = np.insert(allVisits, refSort[goodMatches], refVisit[goodMatches])
1838 allDetectors = np.insert(allDetectors, refSort[goodMatches], refDetector[goodMatches])
1839 allObjectIndices = np.insert(
1840 allObjectIndices, refSort[goodMatches], idx_isoStarCat[goodMatches]
1841 )
1842 issIndices = np.insert(issIndices, refSort[goodMatches], idx_refObjects[goodMatches])
1844 # Loop through the associated sources to convert them to the gbdes
1845 # format, which requires the extension index, the source's index in
1846 # the input table, and a sequence number corresponding to the
1847 # object with which it is associated.
1848 tractExtensions = np.zeros(len(allObjectIndices), dtype=int)
1849 tractObjIndices = np.zeros(len(allObjectIndices), dtype=int)
1850 skippedSources = np.zeros(len(allObjectIndices), dtype=bool)
1851 for visit in np.unique(isolatedStarSources["visit"]):
1852 for detector in np.unique(
1853 isolatedStarSources[isolatedStarSources["visit"] == visit]["detector"]
1854 ):
1855 outputInds = (allVisits == visit) & (allDetectors == detector)
1856 extensionIndex = np.flatnonzero(
1857 (extensionInfo.visit == visit) & (extensionInfo.detector == detector)
1858 )
1859 if len(extensionIndex) == 0:
1860 # This happens for runs where you are not using all the
1861 # visits overlapping a given tract that were included
1862 # in the isolated star association task.
1863 skippedSources[outputInds] = True
1864 continue
1865 else:
1866 extensionIndex = extensionIndex[0]
1868 sourceInds = (isolatedStarSources["visit"] == visit) & (
1869 isolatedStarSources["detector"] == detector
1870 )
1872 tractExtensions[outputInds] = extensionIndex
1873 objIndices = np.arange(sourceInds.sum()) + len(sourceDict[visit][detector]["x"])
1874 tractObjIndices[outputInds] = objIndices
1876 source = isolatedStarSources[sourceInds]
1877 sourceDict[visit][detector]["x"].extend(list(source["x"]))
1878 sourceDict[visit][detector]["y"].extend(list(source["y"]))
1879 xCov = source["xErr"] ** 2
1880 yCov = source["yErr"] ** 2
1881 xyCov = source["ixy"] * (xCov + yCov) / (source["ixx"] + source["iyy"])
1882 # TODO: add correct xyErr if DM-7101 is ever done.
1883 sourceDict[visit][detector]["xCov"].extend(list(xCov))
1884 sourceDict[visit][detector]["yCov"].extend(list(yCov))
1885 sourceDict[visit][detector]["xyCov"].extend(list(xyCov))
1887 # Add ref objects:
1888 for refField in np.unique(allVisits[allVisits <= 0]):
1889 extensionIndex = np.flatnonzero(
1890 (extensionInfo.detector == -1) & (extensionInfo.visit == refField)
1891 )
1892 extensionIndex = extensionIndex[0]
1893 outputInds = (allVisits == refField) & (allDetectors == -1)
1894 tractExtensions[outputInds] = extensionIndex
1895 tractObjIndices[outputInds] = issIndices[outputInds]
1897 uniqueObjects = np.arange(np.array(allObjectIndices).max() + 2)
1898 sorted = np.searchsorted(np.array(allObjectIndices), uniqueObjects)
1899 objectCount = np.diff(sorted)
1900 tractSequence = np.concatenate([np.arange(count) for count in objectCount])
1902 extensions.extend(list(tractExtensions[~skippedSources]))
1903 object_indices.extend(list(tractObjIndices[~skippedSources]))
1904 sequences.extend(list(tractSequence[~skippedSources]))
1906 associations = pipeBase.Struct(extn=extensions, obj=object_indices, sequence=sequences)
1907 return associations, sourceDict
1909 def _check_degeneracies(self, associations, extensionInfo):
1910 """Check that the minimum number of visits and sources needed to
1911 constrain the model are present.
1913 This does not guarantee that the Hessian matrix of the chi-square,
1914 which is used to fit the model, will be positive-definite, but if the
1915 checks here do not pass, the matrix is certain to not be
1916 positive-definite and the model cannot be fit.
1918 Parameters
1919 ----------
1920 associations : `wcsfit.FoFClass`
1921 Object holding the source association information.
1922 extensionInfo : `lsst.pipe.base.Struct`
1923 Struct containing properties for each extension (visit/detector):
1924 ``visit`` : `np.ndarray`
1925 Name of the visit for this extension.
1926 ``detector`` : `np.ndarray`
1927 Name of the detector for this extension.
1928 ``visitIndex` : `np.ndarray` [`int`]
1929 Index of visit for this extension.
1930 ``detectorIndex`` : `np.ndarray` [`int`]
1931 Index of the detector for this extension.
1932 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1933 Initial WCS for this extension.
1934 ``extensionType`` : `np.ndarray` [`str`]
1935 "SCIENCE" or "REFERENCE".
1936 """
1937 # As a baseline, need to have more stars per detector than per-detector
1938 # parameters, and more stars per visit than per-visit parameters.
1939 whichExtension = np.array(associations.extn)
1940 whichDetector = np.zeros(len(whichExtension))
1941 whichVisit = np.zeros(len(whichExtension))
1943 for extension, (detector, visit) in enumerate(zip(extensionInfo.detector, extensionInfo.visit)):
1944 ex_ind = whichExtension == extension
1945 whichDetector[ex_ind] = detector
1946 whichVisit[ex_ind] = visit
1948 if (not self.config.useInputCameraModel) and ("BAND/DEVICE/poly" in self.config.deviceModel):
1949 nCoeffDetectorModel = _nCoeffsFromDegree(self.config.devicePolyOrder)
1950 unconstrainedDetectors = []
1951 for detector in np.unique(extensionInfo.detector):
1952 numSources = (whichDetector == detector).sum()
1953 if numSources < nCoeffDetectorModel: 1953 ↛ 1954line 1953 didn't jump to line 1954 because the condition on line 1953 was never true
1954 unconstrainedDetectors.append(str(detector))
1956 if unconstrainedDetectors: 1956 ↛ 1957line 1956 didn't jump to line 1957 because the condition on line 1956 was never true
1957 raise pipeBase.NoWorkFound(
1958 "The model is not constrained. The following detectors do not have enough "
1959 f"sources ({nCoeffDetectorModel} required): ",
1960 ", ".join(unconstrainedDetectors),
1961 )
1963 def make_yaml(self, inputVisitSummary, inputFile=None, inputCameraModel=None):
1964 """Make a YAML-type object that describes the parameters of the fit
1965 model.
1967 Parameters
1968 ----------
1969 inputVisitSummary : `lsst.afw.table.ExposureCatalog`
1970 Catalog with per-detector summary information.
1971 inputFile : `str`
1972 Path to a file that contains a basic model.
1973 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
1974 Parameters to use for the device part of the model.
1976 Returns
1977 -------
1978 inputYaml : `wcsfit.YAMLCollector`
1979 YAML object containing the model description.
1980 inputDict : `dict` [`str`, `str`]
1981 Dictionary containing the model description.
1982 """
1983 if inputFile is not None: 1983 ↛ 1984line 1983 didn't jump to line 1984 because the condition on line 1983 was never true
1984 inputYaml = wcsfit.YAMLCollector(inputFile, "PixelMapCollection")
1985 else:
1986 inputYaml = wcsfit.YAMLCollector("", "PixelMapCollection")
1987 inputDict = {}
1988 modelComponents = ["BAND/DEVICE", "EXPOSURE"]
1989 baseMap = {"Type": "Composite", "Elements": modelComponents}
1990 inputDict["EXPOSURE/DEVICE/base"] = baseMap
1992 xMin = str(inputVisitSummary["bbox_min_x"].min())
1993 xMax = str(inputVisitSummary["bbox_max_x"].max())
1994 yMin = str(inputVisitSummary["bbox_min_y"].min())
1995 yMax = str(inputVisitSummary["bbox_max_y"].max())
1997 deviceModel = {"Type": "Composite", "Elements": self.config.deviceModel.list()}
1998 inputDict["BAND/DEVICE"] = deviceModel
1999 for component in self.config.deviceModel:
2000 mapping = component.split("/")[-1]
2001 if mapping == "poly": 2001 ↛ 2011line 2001 didn't jump to line 2011 because the condition on line 2001 was always true
2002 componentDict = {
2003 "Type": "Poly",
2004 "XPoly": {"OrderX": self.config.devicePolyOrder, "SumOrder": True},
2005 "YPoly": {"OrderX": self.config.devicePolyOrder, "SumOrder": True},
2006 "XMin": xMin,
2007 "XMax": xMax,
2008 "YMin": yMin,
2009 "YMax": yMax,
2010 }
2011 elif mapping == "identity":
2012 componentDict = {"Type": "Identity"}
2014 inputDict[component] = componentDict
2016 if (inputCameraModel is not None) and self.config.useInputCameraModel:
2017 # This assumes that the input camera model is a 'poly' model
2018 nCoeffs = _nCoeffsFromDegree(self.config.devicePolyOrder)
2019 for key, coeffs in inputCameraModel.items():
2020 if len(coeffs) != nCoeffs * 2: 2020 ↛ 2021line 2020 didn't jump to line 2021 because the condition on line 2020 was never true
2021 raise RuntimeError(
2022 "Input camera model polynomial order does not match the devicePolyOrder"
2023 )
2024 mapDict = {
2025 "Type": "Poly",
2026 "XPoly": {
2027 "OrderX": self.config.devicePolyOrder,
2028 "SumOrder": True,
2029 "Coefficients": coeffs[:nCoeffs].tolist(),
2030 },
2031 "YPoly": {
2032 "OrderX": self.config.devicePolyOrder,
2033 "SumOrder": True,
2034 "Coefficients": coeffs[nCoeffs:].tolist(),
2035 },
2036 "XMin": xMin,
2037 "XMax": xMax,
2038 "YMin": yMin,
2039 "YMax": yMax,
2040 }
2041 inputDict[key] = mapDict
2043 exposureModelComponents = self.config.exposureModel.list()
2044 if self.config.useColor:
2045 exposureModelComponents.append("EXPOSURE/dcr")
2046 exposureModel = {"Type": "Composite", "Elements": exposureModelComponents}
2047 inputDict["EXPOSURE"] = exposureModel
2048 for component in exposureModelComponents:
2049 mapping = component.split("/")[-1]
2050 if mapping == "poly":
2051 componentDict = {
2052 "Type": "Poly",
2053 "XPoly": {"OrderX": self.config.exposurePolyOrder, "SumOrder": "true"},
2054 "YPoly": {"OrderX": self.config.exposurePolyOrder, "SumOrder": "true"},
2055 }
2056 elif mapping == "identity": 2056 ↛ 2057line 2056 didn't jump to line 2057 because the condition on line 2056 was never true
2057 componentDict = {"Type": "Identity"}
2058 elif mapping == "dcr": 2058 ↛ 2065line 2058 didn't jump to line 2065 because the condition on line 2058 was always true
2059 componentDict = {
2060 "Type": "Color",
2061 "Reference": self.config.referenceColor,
2062 "Function": {"Type": "Constant"},
2063 }
2065 inputDict[component] = componentDict
2067 inputYaml.addInput(yaml.dump(inputDict))
2068 inputYaml.addInput("Identity:\n Type: Identity\n")
2070 return inputYaml, inputDict
2072 def _add_objects(self, wcsf, inputCatalogRefs, sourceIndices, extensionInfo, columns):
2073 """Add science sources to the wcsfit.WCSFit object.
2075 Parameters
2076 ----------
2077 wcsf : `wcsfit.WCSFit`
2078 WCS-fitting object.
2079 inputCatalogRefs : `list`
2080 List of DeferredDatasetHandles pointing to visit-level source
2081 tables.
2082 sourceIndices : `list`
2083 List of boolean arrays used to select sources.
2084 extensionInfo : `lsst.pipe.base.Struct`
2085 Struct containing properties for each extension (visit/detector):
2086 ``visit`` : `np.ndarray`
2087 Name of the visit for this extension.
2088 ``detector`` : `np.ndarray`
2089 Name of the detector for this extension.
2090 ``visitIndex` : `np.ndarray` [`int`]
2091 Index of visit for this extension.
2092 ``detectorIndex`` : `np.ndarray` [`int`]
2093 Index of the detector for this extension.
2094 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
2095 Initial WCS for this extension.
2096 ``extensionType`` : `np.ndarray` [`str`]
2097 "SCIENCE" or "REFERENCE".
2098 columns : `list` [`str`]
2099 List of columns needed from source tables.
2100 """
2101 for inputCatalogRef in inputCatalogRefs:
2102 visit = inputCatalogRef.dataId["visit"]
2103 inputCatalog = inputCatalogRef.get(parameters={"columns": columns})
2104 detectors = np.unique(inputCatalog["detector"])
2106 for detector in detectors:
2107 detectorSources = inputCatalog[inputCatalog["detector"] == detector]
2109 extensionIndex = self._find_extension_index(extensionInfo, visit, detector)
2110 if extensionIndex is None:
2111 # This extension does not have information necessary for
2112 # fit. Skip the detections from this detector for this
2113 # visit.
2114 continue
2116 sourceCat = detectorSources[sourceIndices[extensionIndex]]
2118 xCov = sourceCat["xErr"] ** 2
2119 yCov = sourceCat["yErr"] ** 2
2120 xyCov = sourceCat["ixy"] * (xCov + yCov) / (sourceCat["ixx"] + sourceCat["iyy"])
2121 # TODO: add correct xyErr if DM-7101 is ever done.
2123 d = {
2124 "x": sourceCat["x"].to_numpy(),
2125 "y": sourceCat["y"].to_numpy(),
2126 "xCov": xCov.to_numpy(),
2127 "yCov": yCov.to_numpy(),
2128 "xyCov": xyCov.to_numpy(),
2129 }
2131 wcsf.setObjects(
2132 extensionIndex,
2133 d,
2134 "x",
2135 "y",
2136 ["xCov", "yCov", "xyCov"],
2137 defaultColor=self.config.referenceColor,
2138 )
2140 def _add_objects_from_dict(self, wcsf, sourceDict, extensionInfo):
2141 """Add science sources to the wcsfit.WCSFit object.
2143 Parameters
2144 ----------
2145 wcsf : `wcsfit.WCSFit`
2146 WCS-fitting object.
2147 sourceDict : `dict`
2148 Dictionary containing the source centroids for each visit.
2149 extensionInfo : `lsst.pipe.base.Struct`
2150 Struct containing properties for each extension (visit/detector).
2151 ``visit`` : `np.ndarray`
2152 Name of the visit for this extension.
2153 ``detector`` : `np.ndarray`
2154 Name of the detector for this extension.
2155 ``visitIndex` : `np.ndarray` [`int`]
2156 Index of visit for this extension.
2157 ``detectorIndex`` : `np.ndarray` [`int`]
2158 Index of the detector for this extension.
2159 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
2160 Initial WCS for this extension.
2161 ``extensionType`` : `np.ndarray` [`str`]
2162 "SCIENCE" or "REFERENCE".
2163 """
2164 for visit, visitSources in sourceDict.items():
2165 # Visit numbers equal or below zero connote the reference catalog.
2166 if visit <= 0:
2167 # This "visit" number corresponds to a reference catalog.
2168 continue
2170 for detector, sourceCat in visitSources.items():
2171 extensionIndex = np.flatnonzero(
2172 (extensionInfo.visit == visit) & (extensionInfo.detector == detector)
2173 )[0]
2175 d = {
2176 "x": np.array(sourceCat["x"]),
2177 "y": np.array(sourceCat["y"]),
2178 "xCov": np.array(sourceCat["xCov"]),
2179 "yCov": np.array(sourceCat["yCov"]),
2180 "xyCov": np.array(sourceCat["xyCov"]),
2181 }
2182 wcsf.setObjects(
2183 extensionIndex,
2184 d,
2185 "x",
2186 "y",
2187 ["xCov", "yCov", "xyCov"],
2188 defaultColor=self.config.referenceColor,
2189 )
2191 def _add_ref_objects(self, wcsf, refObjects, refCovariance, extensionInfo, fieldIndex=0):
2192 """Add reference sources to the wcsfit.WCSFit object.
2194 Parameters
2195 ----------
2196 wcsf : `wcsfit.WCSFit`
2197 WCS-fitting object.
2198 refObjects : `dict`
2199 Position and error information of reference objects.
2200 refCovariance : `list` [`float`]
2201 Flattened output covariance matrix.
2202 extensionInfo : `lsst.pipe.base.Struct`
2203 Struct containing properties for each extension (visit/detector):
2204 ``visit`` : `np.ndarray`
2205 Name of the visit for this extension.
2206 ``detector`` : `np.ndarray`
2207 Name of the detector for this extension.
2208 ``visitIndex` : `np.ndarray` [`int`]
2209 Index of visit for this extension.
2210 ``detectorIndex`` : `np.ndarray` [`int`]
2211 Index of the detector for this extension.
2212 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
2213 Initial WCS for this extension.
2214 ``extensionType`` : `np.ndarray` [`str`]
2215 "SCIENCE" or "REFERENCE".
2216 fieldIndex : `int`, optional
2217 Index of the field to which these sources correspond.
2218 """
2219 extensionIndex = np.flatnonzero(
2220 (extensionInfo.extensionType == "REFERENCE") & (extensionInfo.visit == fieldIndex)
2221 )[0]
2222 if self.config.fitProperMotion: 2222 ↛ 2223line 2222 didn't jump to line 2223 because the condition on line 2222 was never true
2223 wcsf.setObjects(
2224 extensionIndex,
2225 refObjects,
2226 "ra",
2227 "dec",
2228 ["raCov", "decCov", "raDecCov"],
2229 pmDecKey="decPM",
2230 pmRaKey="raPM",
2231 parallaxKey="parallax",
2232 pmCovKey="fullCov",
2233 pmCov=refCovariance,
2234 )
2235 else:
2236 wcsf.setObjects(extensionIndex, refObjects, "ra", "dec", ["raCov", "decCov", "raDecCov"])
2238 def _add_color_objects(self, wcsf, colorCatalog):
2239 """Associate input matches with objects in color catalog and set their
2240 color value.
2242 Parameters
2243 ----------
2244 wcsf : `wcsfit.WCSFit`
2245 WCSFit object, assumed to have fit model.
2246 colorCatalog : `lsst.afw.table.SimpleCatalog`
2247 Catalog containing object coordinates and magnitudes.
2248 """
2250 # Get current best position for matches
2251 starCat = wcsf.getStarCatalog()
2253 # TODO: DM-45650, update how the colors are read in here.
2254 catalogBands = colorCatalog.metadata.getArray("BANDS")
2255 colorInd1 = catalogBands.index(self.config.color[0])
2256 colorInd2 = catalogBands.index(self.config.color[1])
2257 colors = colorCatalog["mag_std_noabs"][:, colorInd1] - colorCatalog["mag_std_noabs"][:, colorInd2]
2258 goodInd = (colorCatalog["mag_std_noabs"][:, colorInd1] != 99.0) & (
2259 colorCatalog["mag_std_noabs"][:, colorInd2] != 99.0
2260 )
2262 with Matcher(np.array(starCat["starX"]), np.array(starCat["starY"])) as matcher:
2263 idx, idx_starCat, idx_colorCat, d = matcher.query_radius(
2264 (colorCatalog[goodInd]["coord_ra"] * u.radian).to(u.degree).value,
2265 (colorCatalog[goodInd]["coord_dec"] * u.radian).to(u.degree).value,
2266 self.config.matchRadius / 3600.0,
2267 return_indices=True,
2268 )
2270 matchesWithColor = starCat["starMatchID"][idx_starCat]
2271 matchColors = np.ones(len(matchesWithColor)) * self.config.referenceColor
2272 matchColors = colors[goodInd][idx_colorCat]
2273 wcsf.setColors(matchesWithColor, matchColors)
2275 def _make_afw_wcs(self, mapDict, centerRA, centerDec, doNormalizePixels=False, xScale=1, yScale=1):
2276 """Make an `lsst.afw.geom.SkyWcs` from a dictionary of mappings.
2278 Parameters
2279 ----------
2280 mapDict : `dict`
2281 Dictionary of mapping parameters.
2282 centerRA : `lsst.geom.Angle`
2283 RA of the tangent point.
2284 centerDec : `lsst.geom.Angle`
2285 Declination of the tangent point.
2286 doNormalizePixels : `bool`
2287 Whether to normalize pixels so that range is [-1,1].
2288 xScale : `float`
2289 Factor by which to normalize x-dimension. Corresponds to width of
2290 detector.
2291 yScale : `float`
2292 Factor by which to normalize y-dimension. Corresponds to height of
2293 detector.
2295 Returns
2296 -------
2297 outWCS : `lsst.afw.geom.SkyWcs`
2298 WCS constructed from the input mappings
2299 """
2300 # Set up pixel frames
2301 pixelFrame = astshim.Frame(2, "Domain=PIXELS")
2302 normedPixelFrame = astshim.Frame(2, "Domain=NORMEDPIXELS")
2304 if doNormalizePixels: 2304 ↛ 2309line 2304 didn't jump to line 2309 because the condition on line 2304 was always true
2305 # Pixels will need to be rescaled before going into the mappings
2306 normCoefficients = [-1.0, 2.0 / xScale, 0, -1.0, 0, 2.0 / yScale]
2307 normMap = _convert_to_ast_polymap_coefficients(normCoefficients)
2308 else:
2309 normMap = astshim.UnitMap(2)
2311 # All of the detectors for one visit map to the same tangent plane
2312 tangentPoint = lsst.geom.SpherePoint(centerRA, centerDec)
2313 cdMatrix = afwgeom.makeCdMatrix(1.0 * lsst.geom.degrees, 0 * lsst.geom.degrees, True)
2314 iwcToSkyWcs = afwgeom.makeSkyWcs(lsst.geom.Point2D(0, 0), tangentPoint, cdMatrix)
2315 iwcToSkyMap = iwcToSkyWcs.getFrameDict().getMapping("PIXELS", "SKY")
2316 skyFrame = iwcToSkyWcs.getFrameDict().getFrame("SKY")
2318 frameDict = astshim.FrameDict(pixelFrame)
2319 frameDict.addFrame("PIXELS", normMap, normedPixelFrame)
2321 currentFrameName = "NORMEDPIXELS"
2323 # Dictionary values are ordered according to the maps' application.
2324 for m, mapElement in enumerate(mapDict.values()):
2325 mapType = mapElement["Type"]
2327 if mapType == "Poly":
2328 mapCoefficients = mapElement["Coefficients"]
2329 astMap = _convert_to_ast_polymap_coefficients(mapCoefficients)
2330 elif mapType == "Identity": 2330 ↛ 2333line 2330 didn't jump to line 2333 because the condition on line 2330 was always true
2331 astMap = astshim.UnitMap(2)
2332 else:
2333 raise ValueError(f"Converting map type {mapType} to WCS is not supported")
2335 if m == len(mapDict) - 1:
2336 newFrameName = "IWC"
2337 else:
2338 newFrameName = "INTERMEDIATE" + str(m)
2339 newFrame = astshim.Frame(2, f"Domain={newFrameName}")
2340 frameDict.addFrame(currentFrameName, astMap, newFrame)
2341 currentFrameName = newFrameName
2342 frameDict.addFrame("IWC", iwcToSkyMap, skyFrame)
2344 outWCS = afwgeom.SkyWcs(frameDict)
2345 return outWCS
2347 def _make_outputs(
2348 self,
2349 wcsf,
2350 visitSummaryTables,
2351 exposureInfo,
2352 mapTemplate,
2353 extensionInfo,
2354 inputCameraModel=None,
2355 inputCamera=None,
2356 ):
2357 """Make a WCS object out of the WCS models.
2359 Parameters
2360 ----------
2361 wcsf : `wcsfit.WCSFit`
2362 WCSFit object, assumed to have fit model.
2363 visitSummaryTables : `list` [`lsst.afw.table.ExposureCatalog`]
2364 Catalogs with per-detector summary information from which to grab
2365 detector information.
2366 extensionInfo : `lsst.pipe.base.Struct`
2367 Struct containing properties for each extension (visit/detector):
2368 ``visit`` : `np.ndarray`
2369 Name of the visit for this extension.
2370 ``detector`` : `np.ndarray`
2371 Name of the detector for this extension.
2372 ``visitIndex` : `np.ndarray` [`int`]
2373 Index of visit for this extension.
2374 ``detectorIndex`` : `np.ndarray` [`int`]
2375 Index of the detector for this extension.
2376 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
2377 Initial WCS for this extension.
2378 ``extensionType`` : `np.ndarray` [`str`]
2379 "SCIENCE" or "REFERENCE".
2380 mapTemplate : `dict` [`str`, `str`]
2381 Dictionary containing the model description.
2382 extensionInfo : `lsst.pipe.base.Struct`
2383 Struct containing properties for each extension (visit/detector).
2384 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
2385 Parameters to use for the device part of the model. This must be
2386 provided if an input camera model was used.
2388 Returns
2389 -------
2390 catalogs : `dict` [`str`, `lsst.afw.table.ExposureCatalog`]
2391 Dictionary of `lsst.afw.table.ExposureCatalog` objects with the WCS
2392 set to the WCS fit in wcsf, keyed by the visit name.
2393 cameraParams : `dict` [`str`, `np.ndarray`], optional
2394 Parameters for the device part of the model in the format needed
2395 when used as input for future runs.
2396 colorFits : `dict` [`int`, `np.ndarray`], optional
2397 DCR parameters fit in RA and Dec directions for each visit.
2398 camera : `lsst.afw.cameraGeom.Camera`, optional
2399 Camera object constructed from the per-detector model.
2400 """
2401 # Get the parameters of the fit models
2402 mapParams = wcsf.mapCollection.getParamDict()
2403 cameraParams = {}
2404 if self.config.saveCameraModel:
2405 for element in mapTemplate["BAND/DEVICE"]["Elements"]:
2406 for detector in exposureInfo.detectors:
2407 detectorTemplate = element.replace("DEVICE", str(detector))
2408 detectorTemplate = detectorTemplate.replace("BAND", ".+")
2409 for k, params in mapParams.items():
2410 if re.fullmatch(detectorTemplate, k):
2411 cameraParams[k] = params
2412 if self.config.saveCameraObject: 2412 ↛ 2413line 2412 didn't jump to line 2413 because the condition on line 2412 was never true
2413 if inputCamera is None:
2414 raise RuntimeError(
2415 "inputCamera must be provided to _make_outputs in order to build output camera."
2416 )
2417 camera = self.buildCamera.run(
2418 mapParams,
2419 mapTemplate,
2420 exposureInfo.detectors,
2421 exposureInfo.visits,
2422 inputCamera,
2423 )
2424 else:
2425 camera = None
2426 if self.config.useInputCameraModel:
2427 if inputCameraModel is None: 2427 ↛ 2428line 2427 didn't jump to line 2428 because the condition on line 2427 was never true
2428 raise RuntimeError(
2429 "inputCameraModel must be provided to _make_outputs in order to build output WCS."
2430 )
2431 mapParams.update(inputCameraModel)
2433 # Set up the schema for the output catalogs
2434 schema = lsst.afw.table.ExposureTable.makeMinimalSchema()
2435 schema.addField("visit", type="L", doc="Visit number")
2436 schema.addField(
2437 "recoveredWcs",
2438 type="Flag",
2439 doc="Input WCS missing, output recovered from other input visit/detectors.",
2440 )
2442 # Pixels will need to be rescaled before going into the mappings
2443 xscale = int(mapTemplate["BAND/DEVICE/poly"]["XMax"]) - int(mapTemplate["BAND/DEVICE/poly"]["XMin"])
2444 yscale = int(mapTemplate["BAND/DEVICE/poly"]["YMax"]) - int(mapTemplate["BAND/DEVICE/poly"]["YMin"])
2446 # Make dictionary of bboxes for each detector.
2447 detectorBBoxes = {}
2448 for detector in exposureInfo.detectors:
2449 for visitSummary in visitSummaryTables: 2449 ↛ 2448line 2449 didn't jump to line 2448 because the loop on line 2449 didn't complete
2450 if detInfo := visitSummary.find(detector): 2450 ↛ 2449line 2450 didn't jump to line 2449 because the condition on line 2450 was always true
2451 detectorBBoxes[detector] = detInfo.getBBox()
2452 break
2454 catalogs = {}
2455 colorFits = {}
2456 partialOutputs = False
2457 for v, visitSummary in enumerate(visitSummaryTables):
2458 visit = visitSummary[0]["visit"]
2459 if visit not in extensionInfo.visit: 2459 ↛ 2460line 2459 didn't jump to line 2460 because the condition on line 2459 was never true
2460 self.log.warning("Visit %d was dropped because no detectors had valid WCSs.", visit)
2461 partialOutputs = True
2462 continue
2464 visitMaps = wcsf.mapCollection.orderAtoms(f"{visit}")
2465 visitMap = visitMaps[0]
2466 visitMapType = wcsf.mapCollection.getMapType(visitMap)
2467 if (visitMap not in mapParams) and (visitMapType != "Identity"): 2467 ↛ 2468line 2467 didn't jump to line 2468 because the condition on line 2467 was never true
2468 self.log.warning("Visit %d was dropped because of an insufficient amount of data.", visit)
2469 partialOutputs = True
2470 continue
2472 if self.config.useColor:
2473 if f"{visit}/dcr" not in visitMaps: 2473 ↛ 2474line 2473 didn't jump to line 2474 because the condition on line 2473 was never true
2474 raise RuntimeError("DCR map missing for visit %s", visit)
2475 colorMap = visitMaps.pop(visitMaps.index(f"{visit}/dcr"))
2476 colorFits[visit] = mapParams[colorMap]
2478 catalog = lsst.afw.table.ExposureCatalog(schema)
2479 catalog.resize(len(exposureInfo.detectors))
2480 catalog["visit"] = visit
2482 for d, detector in enumerate(exposureInfo.detectors):
2483 mapName = f"{visit}/{detector}"
2484 if mapName in wcsf.mapCollection.allMapNames():
2485 mapElements = wcsf.mapCollection.orderAtoms(f"{mapName}/base")
2486 catalog[d]["recoveredWcs"] = False
2487 else:
2488 # This extension was not fit, but the WCS can be recovered
2489 # using the maps fit from sources on other visits but the
2490 # same detector and from sources on other detectors from
2491 # this visit.
2492 genericElements = mapTemplate["EXPOSURE/DEVICE/base"]["Elements"]
2493 mapElements = []
2494 band = visitSummary[0]["physical_filter"]
2496 # Go through the generic map components to build the names
2497 # of the specific maps for this extension.
2498 for component in genericElements:
2499 elements = mapTemplate[component]["Elements"]
2500 for element in elements:
2501 # TODO: DM-42519, gbdes sets the "BAND" to the
2502 # instrument name currently. This will need to be
2503 # disambiguated if we run on multiple bands at
2504 # once.
2505 element = element.replace("BAND", str(band))
2506 element = element.replace("EXPOSURE", str(visit))
2507 element = element.replace("DEVICE", str(detector))
2508 mapElements.append(element)
2509 catalog[d]["recoveredWcs"] = True
2510 mapDict = {}
2511 for m, mapElement in enumerate(mapElements):
2512 mapType = wcsf.mapCollection.getMapType(mapElement)
2513 if mapType == "Color":
2514 # DCR fit should not go into the generic WCS.
2515 continue
2516 mapDict[mapElement] = {"Type": mapType}
2518 if mapType == "Poly":
2519 mapCoefficients = mapParams[mapElement]
2520 mapDict[mapElement]["Coefficients"] = mapCoefficients
2522 # The RA and Dec of the visit are needed for the last step of
2523 # the mapping from the visit tangent plane to RA and Dec
2524 outWCS = self._make_afw_wcs(
2525 mapDict,
2526 exposureInfo.ras[v] * lsst.geom.radians,
2527 exposureInfo.decs[v] * lsst.geom.radians,
2528 doNormalizePixels=True,
2529 xScale=xscale,
2530 yScale=yscale,
2531 )
2533 catalog[d].setId(detector)
2534 catalog[d].setBBox(detectorBBoxes[detector])
2535 catalog[d].setWcs(outWCS)
2536 catalog.sort()
2537 catalogs[visit] = catalog
2538 if self.config.useColor:
2539 colorVisits = np.array(list(colorFits.keys()))
2540 colorRA = np.array([colorFits[vis][0] for vis in colorVisits])
2541 colorDec = np.array([colorFits[vis][1] for vis in colorVisits])
2542 colorFits = {"visit": colorVisits, "raCoefficient": colorRA, "decCoefficient": colorDec}
2544 return catalogs, cameraParams, colorFits, camera, partialOutputs
2546 def _compute_model_params(self, wcsf):
2547 """Get the WCS model parameters and covariance and convert to a
2548 dictionary that will be readable as a pandas dataframe or other table.
2550 Parameters
2551 ----------
2552 wcsf : `wcsfit.WCSFit`
2553 WCSFit object, assumed to have fit model.
2555 Returns
2556 -------
2557 modelParams : `dict`
2558 Parameters and covariance of the best-fit WCS model.
2559 """
2560 modelParamDict = wcsf.mapCollection.getParamDict()
2561 modelCovariance = wcsf.getModelCovariance()
2563 modelParams = {k: [] for k in ["mapName", "coordinate", "parameter", "coefficientNumber"]}
2564 i = 0
2565 for mapName, params in modelParamDict.items():
2566 nCoeffs = len(params)
2567 # There are an equal number of x and y coordinate parameters
2568 nCoordCoeffs = nCoeffs // 2
2569 modelParams["mapName"].extend([mapName] * nCoeffs)
2570 modelParams["coordinate"].extend(["x"] * nCoordCoeffs)
2571 modelParams["coordinate"].extend(["y"] * nCoordCoeffs)
2572 modelParams["parameter"].extend(params)
2573 modelParams["coefficientNumber"].extend(np.arange(nCoordCoeffs))
2574 modelParams["coefficientNumber"].extend(np.arange(nCoordCoeffs))
2576 for p in range(nCoeffs):
2577 if p < nCoordCoeffs:
2578 coord = "x"
2579 else:
2580 coord = "y"
2581 modelParams[f"{mapName}_{coord}_{p}_cov"] = modelCovariance[i]
2582 i += 1
2584 # Convert the dictionary values from lists to numpy arrays.
2585 for key, value in modelParams.items():
2586 modelParams[key] = np.array(value)
2588 return modelParams
2591class GbdesAstrometricMultibandFitConnections(
2592 GbdesAstrometricFitConnections, dimensions=("skymap", "tract", "instrument")
2593):
2594 outputCatalog = pipeBase.connectionTypes.Output(
2595 doc=(
2596 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2597 "plane coordinates and chisq values."
2598 ),
2599 name="gbdesAstrometricMultibandFit_fitStars",
2600 storageClass="ArrowNumpyDict",
2601 dimensions=("instrument", "skymap", "tract"),
2602 )
2603 starCatalog = pipeBase.connectionTypes.Output(
2604 doc=(
2605 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2606 "fitProperMotion is True."
2607 ),
2608 name="gbdesAstrometricMultibandFit_starCatalog",
2609 storageClass="ArrowNumpyDict",
2610 dimensions=("instrument", "skymap", "tract"),
2611 )
2612 modelParams = pipeBase.connectionTypes.Output(
2613 doc="WCS parameters and covariance.",
2614 name="gbdesAstrometricMultibandFit_modelParams",
2615 storageClass="ArrowNumpyDict",
2616 dimensions=("instrument", "skymap", "tract"),
2617 )
2620class GbdesAstrometricMultibandFitConfig(
2621 GbdesAstrometricFitConfig, pipelineConnections=GbdesAstrometricMultibandFitConnections
2622):
2623 pass
2626class GbdesAstrometricMultibandFitTask(GbdesAstrometricFitTask):
2627 """Calibrate the WCS across multiple visits in multiple filters of the same
2628 field using the GBDES package.
2629 """
2631 ConfigClass = GbdesAstrometricMultibandFitConfig
2632 _DefaultName = "gbdesAstrometricMultibandFit"
2635class GbdesGlobalAstrometricFitConnections(
2636 pipeBase.PipelineTaskConnections, dimensions=("instrument", "physical_filter")
2637):
2638 inputVisitSummaries = pipeBase.connectionTypes.Input(
2639 doc=(
2640 "Per-visit consolidated exposure metadata built from calexps. "
2641 "These catalogs use detector id for the id and must be sorted for "
2642 "fast lookups of a detector."
2643 ),
2644 name="visitSummary",
2645 storageClass="ExposureCatalog",
2646 dimensions=("instrument", "visit"),
2647 multiple=True,
2648 )
2649 referenceCatalog = pipeBase.connectionTypes.PrerequisiteInput(
2650 doc="The astrometry reference catalog to match to loaded input catalog sources.",
2651 name="gaia_dr3_20230707",
2652 storageClass="SimpleCatalog",
2653 dimensions=("skypix",),
2654 deferLoad=True,
2655 multiple=True,
2656 )
2657 colorCatalog = pipeBase.connectionTypes.Input(
2658 doc="The catalog of magnitudes to match to input sources.",
2659 name="fgcm_Cycle4_StandardStars",
2660 storageClass="SimpleCatalog",
2661 dimensions=("instrument",),
2662 )
2663 isolatedStarSources = pipeBase.connectionTypes.Input(
2664 doc="Catalog of matched sources.",
2665 name="isolated_star_presources",
2666 storageClass="DataFrame",
2667 dimensions=(
2668 "instrument",
2669 "skymap",
2670 "tract",
2671 ),
2672 multiple=True,
2673 deferLoad=True,
2674 )
2675 isolatedStarCatalogs = pipeBase.connectionTypes.Input(
2676 doc="Catalog of objects corresponding to the isolatedStarSources.",
2677 name="isolated_star_presource_associations",
2678 storageClass="DataFrame",
2679 dimensions=(
2680 "instrument",
2681 "skymap",
2682 "tract",
2683 ),
2684 multiple=True,
2685 deferLoad=True,
2686 )
2687 inputCameraModel = pipeBase.connectionTypes.PrerequisiteInput(
2688 doc="Camera parameters to use for 'device' part of model",
2689 name="gbdesAstrometricFit_cameraModel",
2690 storageClass="ArrowNumpyDict",
2691 dimensions=("instrument", "physical_filter"),
2692 )
2693 inputCamera = pipeBase.connectionTypes.PrerequisiteInput(
2694 doc="Input camera to use when constructing camera from astrometric model.",
2695 name="camera",
2696 storageClass="Camera",
2697 dimensions=("instrument",),
2698 isCalibration=True,
2699 )
2700 outputWcs = pipeBase.connectionTypes.Output(
2701 doc=(
2702 "Per-visit world coordinate systems derived from the fitted model. These catalogs only contain "
2703 "entries for detectors with an output, and use the detector id for the catalog id, sorted on id "
2704 "for fast lookups of a detector."
2705 ),
2706 name="gbdesGlobalAstrometricFitSkyWcsCatalog",
2707 storageClass="ExposureCatalog",
2708 dimensions=("instrument", "visit"),
2709 multiple=True,
2710 )
2711 outputCatalog = pipeBase.connectionTypes.Output(
2712 doc=(
2713 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2714 "plane coordinates and chisq values."
2715 ),
2716 name="gbdesGlobalAstrometricFit_fitStars",
2717 storageClass="ArrowNumpyDict",
2718 dimensions=("instrument", "physical_filter"),
2719 )
2720 starCatalog = pipeBase.connectionTypes.Output(
2721 doc=(
2722 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2723 "fitProperMotion is True."
2724 ),
2725 name="gbdesGlobalAstrometricFit_starCatalog",
2726 storageClass="ArrowNumpyDict",
2727 dimensions=("instrument", "physical_filter"),
2728 )
2729 modelParams = pipeBase.connectionTypes.Output(
2730 doc="WCS parameters and covariance.",
2731 name="gbdesGlobalAstrometricFit_modelParams",
2732 storageClass="ArrowNumpyDict",
2733 dimensions=("instrument", "physical_filter"),
2734 )
2735 outputCameraModel = pipeBase.connectionTypes.Output(
2736 doc="Camera parameters to use for 'device' part of model",
2737 name="gbdesAstrometricFit_cameraModel",
2738 storageClass="ArrowNumpyDict",
2739 dimensions=("instrument", "physical_filter"),
2740 )
2741 dcrCoefficients = pipeBase.connectionTypes.Output(
2742 doc="Per-visit coefficients for DCR correction.",
2743 name="gbdesGlobalAstrometricFit_dcrCoefficients",
2744 storageClass="ArrowNumpyDict",
2745 )
2746 camera = pipeBase.connectionTypes.Output(
2747 doc="Camera object constructed using the per-detector part of the astrometric model",
2748 name="astrometry_camera",
2749 storageClass="Camera",
2750 isCalibration=True,
2751 dimensions=("instrument", "physical_filter"),
2752 )
2754 def getSpatialBoundsConnections(self):
2755 return ("inputVisitSummaries",)
2757 def __init__(self, *, config=None):
2758 super().__init__(config=config)
2760 if not self.config.useColor:
2761 self.inputs.remove("colorCatalog")
2762 self.outputs.remove("dcrCoefficients")
2763 if not self.config.saveModelParams:
2764 self.outputs.remove("modelParams")
2765 if not self.config.useInputCameraModel:
2766 self.prerequisiteInputs.remove("inputCameraModel")
2767 if not self.config.saveCameraModel:
2768 self.outputs.remove("outputCameraModel")
2769 if not self.config.saveCameraObject:
2770 self.prerequisiteInputs.remove("inputCamera")
2771 self.outputs.remove("camera")
2774class GbdesGlobalAstrometricFitConfig(
2775 GbdesAstrometricFitConfig, pipelineConnections=GbdesGlobalAstrometricFitConnections
2776):
2777 visitOverlap = pexConfig.Field(
2778 dtype=float,
2779 default=1.0,
2780 doc=(
2781 "The linkage distance threshold above which clustered groups of visits will not be merged "
2782 "together in an agglomerative clustering algorithm. The linkage distance is calculated using the "
2783 "minimum distance between the field-of-view centers of a given visit and all other visits in a "
2784 "group, and is in units of the field-of-view radius. The resulting groups of visits define the "
2785 "fields for the astrometric fit."
2786 ),
2787 )
2790class GbdesGlobalAstrometricFitTask(GbdesAstrometricFitTask):
2791 """Calibrate the WCS across multiple visits and multiple fields using the
2792 GBDES package.
2794 This class assumes that the input visits can be separated into contiguous
2795 groups, for which an individual group covers an area of less than a
2796 hemisphere.
2797 """
2799 ConfigClass = GbdesGlobalAstrometricFitConfig
2800 _DefaultName = "gbdesAstrometricFit"
2802 def runQuantum(self, butlerQC, inputRefs, outputRefs):
2803 # We override runQuantum to set up the refObjLoaders
2804 inputs = butlerQC.get(inputRefs)
2806 instrumentName = butlerQC.quantum.dataId["instrument"]
2808 # Ensure the inputs are in a consistent and deterministic order
2809 inputSumVisits = np.array([inputSum[0]["visit"] for inputSum in inputs["inputVisitSummaries"]])
2810 inputs["inputVisitSummaries"] = [inputs["inputVisitSummaries"][v] for v in inputSumVisits.argsort()]
2811 inputRefHtm7s = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputRefs.referenceCatalog])
2812 inputRefCatRefs = [inputRefs.referenceCatalog[htm7] for htm7 in inputRefHtm7s.argsort()]
2813 inputRefCats = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputs["referenceCatalog"]])
2814 inputs["referenceCatalog"] = [inputs["referenceCatalog"][v] for v in inputRefCats.argsort()]
2815 inputIsolatedStarSourceTracts = np.array(
2816 [isolatedStarSource.dataId["tract"] for isolatedStarSource in inputs["isolatedStarSources"]]
2817 )
2818 inputIsolatedStarCatalogTracts = np.array(
2819 [isolatedStarCatalog.dataId["tract"] for isolatedStarCatalog in inputs["isolatedStarCatalogs"]]
2820 )
2821 for tract in inputIsolatedStarCatalogTracts:
2822 if tract not in inputIsolatedStarSourceTracts:
2823 raise RuntimeError(f"tract {tract} in isolated_star_cats but not isolated_star_sources")
2824 inputs["isolatedStarSources"] = np.array(
2825 [inputs["isolatedStarSources"][t] for t in inputIsolatedStarSourceTracts.argsort()]
2826 )
2827 inputs["isolatedStarCatalogs"] = np.array(
2828 [inputs["isolatedStarCatalogs"][t] for t in inputIsolatedStarSourceTracts.argsort()]
2829 )
2831 refConfig = LoadReferenceObjectsConfig()
2832 if self.config.applyRefCatProperMotion:
2833 refConfig.requireProperMotion = True
2834 refObjectLoader = ReferenceObjectLoader(
2835 dataIds=[ref.datasetRef.dataId for ref in inputRefCatRefs],
2836 refCats=inputs.pop("referenceCatalog"),
2837 config=refConfig,
2838 log=self.log,
2839 )
2841 nCores = butlerQC.resources.num_cores
2842 self.log.info("Running with nCores = %d", nCores)
2844 if self.config.useColor:
2845 colorCatalog = inputs.pop("colorCatalog")
2846 else:
2847 colorCatalog = None
2849 try:
2850 output = self.run(
2851 **inputs,
2852 instrumentName=instrumentName,
2853 refObjectLoader=refObjectLoader,
2854 colorCatalog=colorCatalog,
2855 nCores=nCores,
2856 )
2857 except pipeBase.AlgorithmError as e:
2858 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
2859 # No partial outputs for butler to put
2860 raise error from e
2862 wcsOutputRefDict = {outWcsRef.dataId["visit"]: outWcsRef for outWcsRef in outputRefs.outputWcs}
2863 for visit, outputWcs in output.outputWcss.items():
2864 butlerQC.put(outputWcs, wcsOutputRefDict[visit])
2865 butlerQC.put(output.outputCatalog, outputRefs.outputCatalog)
2866 butlerQC.put(output.starCatalog, outputRefs.starCatalog)
2867 if self.config.saveModelParams:
2868 butlerQC.put(output.modelParams, outputRefs.modelParams)
2869 if self.config.saveCameraModel:
2870 butlerQC.put(output.cameraModelParams, outputRefs.outputCameraModel)
2871 if self.config.saveCameraObject:
2872 butlerQC.put(output.camera, outputRefs.camera)
2873 if self.config.useColor:
2874 butlerQC.put(output.colorParams, outputRefs.dcrCoefficients)
2875 if output.partialOutputs:
2876 e = RuntimeError("Some visits were dropped because data was insufficient to fit model.")
2877 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
2878 raise error from e
2880 def run(
2881 self,
2882 inputVisitSummaries,
2883 isolatedStarSources,
2884 isolatedStarCatalogs,
2885 instrumentName="",
2886 refEpoch=None,
2887 refObjectLoader=None,
2888 inputCameraModel=None,
2889 colorCatalog=None,
2890 inputCamera=None,
2891 nCores=1,
2892 ):
2893 """Run the WCS fit for a given set of visits
2895 Parameters
2896 ----------
2897 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
2898 List of catalogs with per-detector summary information.
2899 isolatedStarSources : `list` [`DeferredDatasetHandle`]
2900 List of handles pointing to isolated star sources.
2901 isolatedStarCatalog: `list` [`DeferredDatasetHandle`]
2902 List of handles pointing to isolated star catalogs.
2903 instrumentName : `str`, optional
2904 Name of the instrument used. This is only used for labelling.
2905 refEpoch : `float`, optional
2906 Epoch of the reference objects in MJD.
2907 refObjectLoader : instance of
2908 `lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader`,
2909 optional
2910 Reference object loader instance.
2911 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
2912 Parameters to use for the device part of the model.
2913 colorCatalog : `lsst.afw.table.SimpleCatalog`
2914 Catalog containing object coordinates and magnitudes.
2915 inputCamera : `lsst.afw.cameraGeom.Camera`, optional
2916 Camera to be used as template when constructing new camera.
2917 nCores : `int`, optional
2918 Number of cores to use during WCS fit.
2920 Returns
2921 -------
2922 result : `lsst.pipe.base.Struct`
2923 ``outputWcss`` : `list` [`lsst.afw.table.ExposureCatalog`]
2924 List of exposure catalogs (one per visit) with the WCS for each
2925 detector set by the new fitted WCS.
2926 ``fitModel`` : `wcsfit.WCSFit`
2927 Model-fitting object with final model parameters.
2928 ``outputCatalog`` : `pyarrow.Table`
2929 Catalog with fit residuals of all sources used.
2930 ``starCatalog`` : `pyarrow.Table`
2931 Catalog with best-fit positions of the objects fit.
2932 ``modelParams`` : `dict`
2933 Parameters and covariance of the best-fit WCS model.
2934 ``cameraModelParams`` : `dict` [`str`, `np.ndarray`]
2935 Parameters of the device part of the model, in the format
2936 needed as input for future runs.
2937 ``colorParams`` : `dict` [`int`, `np.ndarray`]
2938 DCR parameters fit in RA and Dec directions for each visit.
2939 ``camera`` : `lsst.afw.cameraGeom.Camera`
2940 Camera object constructed from the per-detector model.
2941 """
2942 self.log.info("Gather instrument, exposure, and field info")
2944 # Get information about the extent of the input visits
2945 fields, fieldRegions = self._prep_sky(inputVisitSummaries)
2947 # Get RA, Dec, MJD, etc., for the input visits
2948 exposureInfo, exposuresHelper, extensionInfo, instruments = self._get_exposure_info(
2949 inputVisitSummaries, fieldRegions=fieldRegions
2950 )
2952 self.log.info("Load associated sources")
2953 medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format="jyear").mjd
2954 allRefObjects, allRefCovariances = {}, {}
2955 for f, fieldRegion in fieldRegions.items():
2956 refObjects, refCovariance = self._load_refcat(
2957 refObjectLoader, extensionInfo, epoch=medianEpoch, region=fieldRegion
2958 )
2959 allRefObjects[f] = refObjects
2960 allRefCovariances[f] = refCovariance
2962 associations, sourceDict = self._associate_from_isolated_sources(
2963 isolatedStarSources, isolatedStarCatalogs, extensionInfo, allRefObjects
2964 )
2966 self.log.info("Fit the WCSs")
2967 # Set up a YAML-type string using the config variables and a sample
2968 # visit
2969 inputYaml, mapTemplate = self.make_yaml(
2970 inputVisitSummaries[0],
2971 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
2972 )
2974 # Set the verbosity level for WCSFit from the task log level.
2975 # TODO: DM-36850, Add lsst.log to gbdes so that log messages are
2976 # properly propagated.
2977 loglevel = self.log.getEffectiveLevel()
2978 if loglevel >= self.log.WARNING: 2978 ↛ 2979line 2978 didn't jump to line 2979 because the condition on line 2978 was never true
2979 verbose = 0
2980 elif loglevel == self.log.INFO: 2980 ↛ 2981line 2980 didn't jump to line 2981 because the condition on line 2980 was never true
2981 verbose = 1
2982 else:
2983 verbose = 2
2985 # Set up the WCS-fitting class using the source matches from the
2986 # isolated star sources plus the reference catalog.
2987 wcsf = wcsfit.WCSFit(
2988 fields,
2989 instruments,
2990 exposuresHelper,
2991 extensionInfo.visitIndex,
2992 extensionInfo.detectorIndex,
2993 inputYaml,
2994 extensionInfo.wcs,
2995 associations.sequence,
2996 associations.extn,
2997 associations.obj,
2998 sysErr=self.config.systematicError,
2999 refSysErr=self.config.referenceSystematicError,
3000 usePM=self.config.fitProperMotion,
3001 verbose=verbose,
3002 num_threads=nCores,
3003 )
3005 # Add the science and reference sources
3006 self._add_objects_from_dict(wcsf, sourceDict, extensionInfo)
3007 for f in fieldRegions.keys():
3008 self._add_ref_objects(
3009 wcsf, allRefObjects[f], allRefCovariances[f], extensionInfo, fieldIndex=-1 * f
3010 )
3011 if self.config.useColor:
3012 self._add_color_objects(wcsf, colorCatalog)
3014 # Do the WCS fit
3015 try:
3016 wcsf.fit(
3017 reserveFraction=self.config.fitReserveFraction,
3018 randomNumberSeed=self.config.fitReserveRandomSeed,
3019 clipThresh=self.config.clipThresh,
3020 clipFraction=self.config.clipFraction,
3021 )
3022 except RuntimeError as e:
3023 if "Cholesky decomposition failed" in str(e):
3024 raise CholeskyError() from e
3025 else:
3026 raise
3028 self.log.info("WCS fitting done")
3030 outputWcss, cameraParams, colorParams, camera, partialOutputs = self._make_outputs(
3031 wcsf,
3032 inputVisitSummaries,
3033 exposureInfo,
3034 mapTemplate,
3035 extensionInfo,
3036 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
3037 inputCamera=(inputCamera if self.config.buildCamera else None),
3038 )
3039 outputCatalog = wcsf.getOutputCatalog()
3040 outputCatalog["exposureName"] = np.array(outputCatalog["exposureName"])
3041 outputCatalog["deviceName"] = np.array(outputCatalog["deviceName"])
3043 starCatalog = wcsf.getStarCatalog()
3044 modelParams = self._compute_model_params(wcsf) if self.config.saveModelParams else None
3046 return pipeBase.Struct(
3047 outputWcss=outputWcss,
3048 fitModel=wcsf,
3049 outputCatalog=outputCatalog,
3050 starCatalog=starCatalog,
3051 modelParams=modelParams,
3052 cameraModelParams=cameraParams,
3053 colorParams=colorParams,
3054 camera=camera,
3055 partialOutputs=partialOutputs,
3056 )
3058 def _prep_sky(self, inputVisitSummaries):
3059 """Cluster the input visits into disjoint groups that will define
3060 separate fields in the astrometric fit, and, for each group, get the
3061 convex hull around all of its component visits.
3063 The groups are created such that each visit overlaps with at least one
3064 other visit in the same group by the `visitOverlap` amount, which is
3065 calculated as a fraction of the field-of-view radius, and no visits
3066 from separate groups overlap by more than this amount.
3068 Paramaters
3069 ----------
3070 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
3071 List of catalogs with per-detector summary information.
3073 Returns
3074 -------
3075 fields : `wcsfit.Fields`
3076 Object with field information.
3077 fieldRegions : `dict` [`int`, `lsst.sphgeom.ConvexPolygon`]
3078 Dictionary of regions encompassing each group of input visits,
3079 keyed by an arbitrary index.
3080 """
3081 allDetectorCorners = []
3082 mjds = []
3083 radecs = []
3084 radii = []
3085 for visSum in inputVisitSummaries:
3086 detectorCorners = [
3087 lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees).getVector()
3088 for (ra, dec) in zip(visSum["raCorners"].ravel(), visSum["decCorners"].ravel())
3089 if (np.isfinite(ra) and (np.isfinite(dec)))
3090 ]
3091 if len(detectorCorners) == 0: 3091 ↛ 3095line 3091 didn't jump to line 3095 because the condition on line 3091 was never true
3092 # Skip this visit if none of the detectors have finite ra/dec
3093 # corners, which happens when the WCSs are missing. The visit
3094 # will get formally dropped elsewhere.
3095 continue
3096 allDetectorCorners.append(detectorCorners)
3098 # Get center and approximate radius of field of view
3099 boundingCircle = lsst.sphgeom.ConvexPolygon.convexHull(detectorCorners).getBoundingCircle()
3100 center = lsst.geom.SpherePoint(boundingCircle.getCenter())
3101 ra = center.getRa().asDegrees()
3102 dec = center.getDec().asDegrees()
3103 radecs.append([ra, dec])
3104 radius = boundingCircle.getOpeningAngle()
3105 radii.append(radius)
3107 obsDate = visSum[0].getVisitInfo().getDate()
3108 obsMJD = obsDate.get(obsDate.MJD)
3109 mjds.append(obsMJD)
3111 # Find groups of visits where any one of the visits overlaps another by
3112 # a given fraction of the field-of-view radius.
3113 distance = self.config.visitOverlap * np.median(radii)
3114 clustering = AgglomerativeClustering(
3115 distance_threshold=distance.asDegrees(), n_clusters=None, linkage="single"
3116 )
3117 clusters = clustering.fit(np.array(radecs))
3119 medianMJD = np.median(mjds)
3120 medianEpoch = astropy.time.Time(medianMJD, format="mjd").jyear
3122 fieldNames = []
3123 fieldRAs = []
3124 fieldDecs = []
3125 epochs = []
3126 fieldRegions = {}
3128 for i in range(clusters.n_clusters_):
3129 fieldInd = clusters.labels_ == i
3130 # Concatenate the lists of all detector corners that are in this
3131 # field
3132 fieldDetectors = sum([allDetectorCorners[f] for f, fInd in enumerate(fieldInd) if fInd], [])
3133 hull = lsst.sphgeom.ConvexPolygon.convexHull(fieldDetectors)
3134 center = lsst.geom.SpherePoint(hull.getCentroid())
3135 ra = center.getRa().asDegrees()
3136 dec = center.getDec().asDegrees()
3138 fieldRegions[i] = hull
3139 fieldNames.append(str(i))
3140 fieldRAs.append(ra)
3141 fieldDecs.append(dec)
3142 # Use the same median epoch for all fields so that the final object
3143 # positions are calculated for the same epoch.
3144 epochs.append(medianEpoch)
3146 fields = wcsfit.Fields(fieldNames, fieldRAs, fieldDecs, epochs)
3148 return fields, fieldRegions
3151class GbdesGlobalAstrometricMultibandFitConnections(
3152 GbdesGlobalAstrometricFitConnections,
3153 dimensions=("instrument",),
3154):
3155 outputCatalog = pipeBase.connectionTypes.Output(
3156 doc=(
3157 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
3158 "plane coordinates and chisq values."
3159 ),
3160 name="gbdesGlobalAstrometricMultibandFit_fitStars",
3161 storageClass="ArrowNumpyDict",
3162 dimensions=("instrument",),
3163 )
3164 starCatalog = pipeBase.connectionTypes.Output(
3165 doc=(
3166 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
3167 "fitProperMotion is True."
3168 ),
3169 name="gbdesGlobalAstrometricMultibandFit_starCatalog",
3170 storageClass="ArrowNumpyDict",
3171 dimensions=("instrument",),
3172 )
3173 modelParams = pipeBase.connectionTypes.Output(
3174 doc="WCS parameters and covariance.",
3175 name="gbdesGlobalAstrometricMultibandFit_modelParams",
3176 storageClass="ArrowNumpyDict",
3177 dimensions=("instrument",),
3178 )
3181class GbdesGlobalAstrometricMultibandFitConfig(
3182 GbdesAstrometricFitConfig,
3183 pipelineConnections=GbdesGlobalAstrometricMultibandFitConnections,
3184):
3185 """Configuration for the GbdesGlobalAstrometricMultibandFitTask"""
3187 pass
3190class GbdesGlobalAstrometricMultibandFitTask(GbdesGlobalAstrometricFitTask):
3191 """Calibrate the WCS across multiple visits in multiple filters and
3192 multiple fields using the GBDES package.
3193 """
3195 ConfigClass = GbdesGlobalAstrometricMultibandFitConfig
3196 _DefaultName = "gbdesAstrometricMultibandFit"