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