Coverage for python/lsst/drp/tasks/gbdesAstrometricFit.py: 11%
1006 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-12 01:19 -0800
1# This file is part of drp_tasks.
2#
3# LSST Data Management System
4# This product includes software developed by the
5# LSST Project (http://www.lsst.org/).
6# See COPYRIGHT file at the top of the source tree.
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <https://www.lsstcorp.org/LegalNotices/>.
21#
22import dataclasses
23import 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 )
639 minDetectorFraction = pexConfig.Field(
640 dtype=float,
641 doc=(
642 "Minimum fraction of detectors with valid WCSs per visit required to include the visit in the "
643 "fit."
644 ),
645 default=0.25,
646 )
648 def setDefaults(self):
649 # Use only stars because aperture fluxes of galaxies are biased and
650 # depend on seeing.
651 self.sourceSelector["science"].doUnresolved = True
652 self.sourceSelector["science"].unresolved.name = "sizeExtendedness"
654 # Use only isolated sources.
655 self.sourceSelector["science"].doIsolated = True
656 self.sourceSelector["science"].isolated.parentName = "parentSourceId"
657 self.sourceSelector["science"].isolated.nChildName = "deblend_nChild"
658 # Do not use either flux or centroid measurements with flags,
659 # chosen from the usual QA flags for stars.
660 self.sourceSelector["science"].doFlags = True
661 badFlags = [
662 "pixelFlags_edge",
663 "pixelFlags_saturated",
664 "pixelFlags_interpolatedCenter",
665 "pixelFlags_interpolated",
666 "pixelFlags_crCenter",
667 "pixelFlags_bad",
668 "hsmPsfMoments_flag",
669 f"{self.sourceFluxType}_flag",
670 ]
671 self.sourceSelector["science"].flags.bad = badFlags
673 # Use only primary sources.
674 self.sourceSelector["science"].doRequirePrimary = True
676 self.sourceSelector["science"].doSignalToNoise = True
677 self.sourceSelector["science"].signalToNoise.minimum = 8.0
678 self.sourceSelector["science"].signalToNoise.maximum = 1000.0
679 self.sourceSelector["science"].signalToNoise.fluxField = self.sourceFluxType + "_instFlux"
680 self.sourceSelector["science"].signalToNoise.errField = self.sourceFluxType + "_instFluxErr"
682 def validate(self):
683 super().validate()
685 # Check if all components of the device and exposure models are
686 # supported.
687 for component in self.deviceModel:
688 mapping = component.split("/")[-1]
689 if mapping not in ["poly", "identity"]:
690 raise pexConfig.FieldValidationError(
691 GbdesAstrometricFitConfig.deviceModel,
692 self,
693 f"deviceModel component {component} is not supported.",
694 )
696 for component in self.exposureModel:
697 mapping = component.split("/")[-1]
698 if mapping not in ["poly", "identity", "dcr"]:
699 raise pexConfig.FieldValidationError(
700 GbdesAstrometricFitConfig.exposureModel,
701 self,
702 f"exposureModel component {component} is not supported.",
703 )
705 if self.saveCameraModel and self.useInputCameraModel:
706 raise pexConfig.FieldValidationError(
707 GbdesAstrometricFitConfig.saveCameraModel,
708 self,
709 "saveCameraModel and useInputCameraModel cannot both be true.",
710 )
712 if self.saveCameraObject and (self.exposurePolyOrder != 1):
713 raise pexConfig.FieldValidationError(
714 GbdesAstrometricFitConfig.saveCameraObject,
715 self,
716 "If saveCameraObject is True, exposurePolyOrder must be set to 1.",
717 )
720class GbdesAstrometricFitTask(pipeBase.PipelineTask):
721 """Calibrate the WCS across multiple visits of the same field using the
722 GBDES package.
723 """
725 ConfigClass = GbdesAstrometricFitConfig
726 _DefaultName = "gbdesAstrometricFit"
728 def __init__(self, **kwargs):
729 super().__init__(**kwargs)
730 self.makeSubtask("sourceSelector")
731 self.makeSubtask("referenceSelector")
732 if self.config.saveCameraObject:
733 self.makeSubtask("buildCamera")
735 def runQuantum(self, butlerQC, inputRefs, outputRefs):
736 # We override runQuantum to set up the refObjLoaders
737 inputs = butlerQC.get(inputRefs)
739 instrumentName = butlerQC.quantum.dataId["instrument"]
741 # Ensure the inputs are in a consistent and deterministic order
742 inputCatVisits = np.array([inputCat.dataId["visit"] for inputCat in inputs["inputCatalogRefs"]])
743 inputs["inputCatalogRefs"] = [inputs["inputCatalogRefs"][v] for v in inputCatVisits.argsort()]
744 inputSumVisits = np.array([inputSum[0]["visit"] for inputSum in inputs["inputVisitSummaries"]])
745 inputs["inputVisitSummaries"] = [inputs["inputVisitSummaries"][v] for v in inputSumVisits.argsort()]
746 inputRefHtm7s = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputRefs.referenceCatalog])
747 inputRefCatRefs = [inputRefs.referenceCatalog[htm7] for htm7 in inputRefHtm7s.argsort()]
748 inputRefCats = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputs["referenceCatalog"]])
749 inputs["referenceCatalog"] = [inputs["referenceCatalog"][v] for v in inputRefCats.argsort()]
751 refConfig = LoadReferenceObjectsConfig()
752 if self.config.applyRefCatProperMotion:
753 refConfig.requireProperMotion = True
754 refObjectLoader = ReferenceObjectLoader(
755 dataIds=[ref.datasetRef.dataId for ref in inputRefCatRefs],
756 refCats=inputs.pop("referenceCatalog"),
757 config=refConfig,
758 log=self.log,
759 )
761 nCores = butlerQC.resources.num_cores
762 self.log.info("Running with nCores = %d", nCores)
764 if self.config.useColor:
765 colorCatalog = inputs.pop("colorCatalog")
766 else:
767 colorCatalog = None
768 try:
769 output = self.run(
770 **inputs,
771 instrumentName=instrumentName,
772 refObjectLoader=refObjectLoader,
773 colorCatalog=colorCatalog,
774 nCores=nCores,
775 )
776 except pipeBase.AlgorithmError as e:
777 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
778 # No partial outputs for butler to put
779 raise error from e
781 wcsOutputRefDict = {outWcsRef.dataId["visit"]: outWcsRef for outWcsRef in outputRefs.outputWcs}
782 for visit, outputWcs in output.outputWcss.items():
783 butlerQC.put(outputWcs, wcsOutputRefDict[visit])
784 butlerQC.put(output.outputCatalog, outputRefs.outputCatalog)
785 butlerQC.put(output.starCatalog, outputRefs.starCatalog)
786 if self.config.saveModelParams:
787 butlerQC.put(output.modelParams, outputRefs.modelParams)
788 if self.config.saveCameraModel:
789 butlerQC.put(output.cameraModelParams, outputRefs.outputCameraModel)
790 if self.config.saveCameraObject:
791 butlerQC.put(output.camera, outputRefs.camera)
792 if self.config.useColor:
793 butlerQC.put(output.colorParams, outputRefs.dcrCoefficients)
794 if output.partialOutputs:
795 e = RuntimeError("Some visits were dropped because data was insufficient to fit model.")
796 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
797 raise error from e
799 def run(
800 self,
801 inputCatalogRefs,
802 inputVisitSummaries,
803 instrumentName="",
804 refEpoch=None,
805 refObjectLoader=None,
806 inputCameraModel=None,
807 colorCatalog=None,
808 inputCamera=None,
809 nCores=1,
810 ):
811 """Run the WCS fit for a given set of visits
813 Parameters
814 ----------
815 inputCatalogRefs : `list` [`DeferredDatasetHandle`]
816 List of handles pointing to visit-level source
817 tables.
818 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
819 List of catalogs with per-detector summary information.
820 instrumentName : `str`, optional
821 Name of the instrument used. This is only used for labelling.
822 refEpoch : `float`
823 Epoch of the reference objects in MJD.
824 refObjectLoader : instance of
825 `lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader`
826 Reference object loader instance.
827 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
828 Parameters to use for the device part of the model.
829 colorCatalog : `lsst.afw.table.SimpleCatalog`
830 Catalog containing object coordinates and magnitudes.
831 inputCamera : `lsst.afw.cameraGeom.Camera`, optional
832 Camera to be used as template when constructing new camera.
833 nCores : `int`, optional
834 Number of cores to use during WCS fit.
836 Returns
837 -------
838 result : `lsst.pipe.base.Struct`
839 ``outputWcss`` : `list` [`lsst.afw.table.ExposureCatalog`]
840 List of exposure catalogs (one per visit) with the WCS for each
841 detector set by the new fitted WCS.
842 ``fitModel`` : `wcsfit.WCSFit`
843 Model-fitting object with final model parameters.
844 ``outputCatalog`` : `pyarrow.Table`
845 Catalog with fit residuals of all sources used.
846 ``starCatalog`` : `pyarrow.Table`
847 Catalog with best-fit positions of the objects fit.
848 ``modelParams`` : `dict`
849 Parameters and covariance of the best-fit WCS model.
850 ``cameraModelParams`` : `dict` [`str`, `np.ndarray`]
851 Parameters of the device part of the model, in the format
852 needed as input for future runs.
853 ``colorParams`` : `dict` [`int`, `np.ndarray`]
854 DCR parameters fit in RA and Dec directions for each visit.
855 ``camera`` : `lsst.afw.cameraGeom.Camera`
856 Camera object constructed from the per-detector model.
857 """
858 self.log.info("Gather instrument, exposure, and field info")
860 # Get RA, Dec, MJD, etc., for the input visits
861 exposureInfo, exposuresHelper, extensionInfo, instruments = self._get_exposure_info(
862 inputVisitSummaries
863 )
865 # Get information about the extent of the input visits
866 fields, fieldCenter, fieldRadius = self._prep_sky(inputVisitSummaries, exposureInfo.medianEpoch)
867 self.log.info("Field center set at %s with radius %s degrees", fieldCenter, fieldRadius.asDegrees())
869 self.log.info("Load catalogs and associate sources")
870 # Set up class to associate sources into matches using a
871 # friends-of-friends algorithm
872 associations = wcsfit.FoFClass(
873 fields,
874 instruments,
875 exposuresHelper,
876 [fieldRadius.asDegrees()],
877 (self.config.matchRadius * u.arcsec).to(u.degree).value,
878 )
880 # Add the reference catalog to the associator
881 medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format="decimalyear").mjd
882 refObjects, refCovariance = self._load_refcat(
883 refObjectLoader,
884 extensionInfo,
885 epoch=medianEpoch,
886 center=fieldCenter,
887 radius=fieldRadius,
888 associations=associations,
889 )
891 # Add the science catalogs and associate new sources as they are added
892 sourceIndices, usedColumns = self._load_catalogs_and_associate(
893 associations, inputCatalogRefs, extensionInfo
894 )
895 self._check_degeneracies(associations, extensionInfo)
897 self.log.info("Fit the WCSs")
898 # Set up a YAML-type string using the config variables and a sample
899 # visit
900 inputYaml, mapTemplate = self.make_yaml(
901 inputVisitSummaries[0],
902 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
903 )
905 # Set the verbosity level for WCSFit from the task log level.
906 # TODO: DM-36850, Add lsst.log to gbdes so that log messages are
907 # properly propagated.
908 loglevel = self.log.getEffectiveLevel()
909 if loglevel >= self.log.WARNING:
910 verbose = 0
911 elif loglevel == self.log.INFO:
912 verbose = 1
913 else:
914 verbose = 2
916 # Set up the WCS-fitting class using the results of the FOF associator
917 if self.config.useInputCameraModel:
918 fixMaps = ",".join(inputCameraModel.keys())
919 else:
920 fixMaps = ""
921 wcsf = wcsfit.WCSFit(
922 fields,
923 instruments,
924 exposuresHelper,
925 extensionInfo.visitIndex,
926 extensionInfo.detectorIndex,
927 inputYaml,
928 extensionInfo.wcs,
929 associations.sequence,
930 associations.extn,
931 associations.obj,
932 sysErr=self.config.systematicError,
933 refSysErr=self.config.referenceSystematicError,
934 usePM=self.config.fitProperMotion,
935 verbose=verbose,
936 fixMaps=fixMaps,
937 num_threads=nCores,
938 )
940 # Add the science and reference sources
941 self._add_objects(wcsf, inputCatalogRefs, sourceIndices, extensionInfo, usedColumns)
942 self._add_ref_objects(wcsf, refObjects, refCovariance, extensionInfo)
943 if self.config.useColor:
944 self._add_color_objects(wcsf, colorCatalog)
946 # There must be at least as many sources per visit as the number of
947 # free parameters in the per-visit mapping. Set minFitExposures to be
948 # the number of free parameters divided by the fraction of non-reserved
949 # sources, so that visits with fewer sources are dropped.
950 nCoeffVisitModel = _nCoeffsFromDegree(self.config.exposurePolyOrder)
951 minFitExposures = int(np.ceil(nCoeffVisitModel / (1 - self.config.fitReserveFraction)))
952 # Do the WCS fit
953 try:
954 wcsf.fit(
955 reserveFraction=self.config.fitReserveFraction,
956 randomNumberSeed=self.config.fitReserveRandomSeed,
957 minFitExposures=minFitExposures,
958 clipThresh=self.config.clipThresh,
959 clipFraction=self.config.clipFraction,
960 )
961 except RuntimeError as e:
962 if "Cholesky decomposition failed" in str(e):
963 raise CholeskyError() from e
964 else:
965 raise
967 self.log.info("WCS fitting done")
969 outputWcss, cameraParams, colorParams, camera, partialOutputs = self._make_outputs(
970 wcsf,
971 inputVisitSummaries,
972 exposureInfo,
973 mapTemplate,
974 extensionInfo,
975 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
976 inputCamera=(inputCamera if self.config.buildCamera else None),
977 )
978 outputCatalog = wcsf.getOutputCatalog()
979 outputCatalog["exposureName"] = np.array(outputCatalog["exposureName"])
980 outputCatalog["deviceName"] = np.array(outputCatalog["deviceName"])
982 starCatalog = wcsf.getStarCatalog()
983 modelParams = self._compute_model_params(wcsf) if self.config.saveModelParams else None
985 return pipeBase.Struct(
986 outputWcss=outputWcss,
987 fitModel=wcsf,
988 outputCatalog=outputCatalog,
989 starCatalog=starCatalog,
990 modelParams=modelParams,
991 cameraModelParams=cameraParams,
992 colorParams=colorParams,
993 camera=camera,
994 partialOutputs=partialOutputs,
995 )
997 def _prep_sky(self, inputVisitSummaries, epoch, fieldName="Field"):
998 """Get center and radius of the input tract. This assumes that all
999 visits will be put into the same `wcsfit.Field` and fit together.
1001 Paramaters
1002 ----------
1003 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
1004 List of catalogs with per-detector summary information.
1005 epoch : float
1006 Reference epoch.
1007 fieldName : str
1008 Name of the field, used internally.
1010 Returns
1011 -------
1012 fields : `wcsfit.Fields`
1013 Object with field information.
1014 center : `lsst.geom.SpherePoint`
1015 Center of the field.
1016 radius : `lsst.sphgeom._sphgeom.Angle`
1017 Radius of the bounding circle of the tract.
1018 """
1019 allDetectorCorners = []
1020 for visSum in inputVisitSummaries:
1021 detectorCorners = [
1022 lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees).getVector()
1023 for (ra, dec) in zip(visSum["raCorners"].ravel(), visSum["decCorners"].ravel())
1024 if (np.isfinite(ra) and (np.isfinite(dec)))
1025 ]
1026 allDetectorCorners.extend(detectorCorners)
1027 boundingCircle = lsst.sphgeom.ConvexPolygon.convexHull(allDetectorCorners).getBoundingCircle()
1028 center = lsst.geom.SpherePoint(boundingCircle.getCenter())
1029 ra = center.getRa().asDegrees()
1030 dec = center.getDec().asDegrees()
1031 radius = boundingCircle.getOpeningAngle()
1033 # wcsfit.Fields describes a list of fields, but we assume all
1034 # observations will be fit together in one field.
1035 fields = wcsfit.Fields([fieldName], [ra], [dec], [epoch])
1037 return fields, center, radius
1039 def _get_exposure_info(
1040 self,
1041 inputVisitSummaries,
1042 fieldNumber=0,
1043 refEpoch=None,
1044 fieldRegions=None,
1045 ):
1046 """Get various information about the input visits to feed to the
1047 fitting routines.
1049 Parameters
1050 ----------
1051 inputVisitSummaries : `list [`lsst.afw.table.ExposureCatalog`]
1052 Tables for each visit with information for detectors.
1053 fieldNumber : `int`, optional
1054 Index of the field for these visits. Should be zero if all data is
1055 being fit together. This is ignored if `fieldRegions` is not None.
1056 refEpoch : `float`, optional
1057 Epoch of the reference objects in MJD.
1058 fieldRegions : `dict` [`int`, `lsst.sphgeom.ConvexPolygon`], optional
1059 Dictionary of regions encompassing each group of input visits
1060 keyed by an arbitrary index.
1062 Returns
1063 -------
1064 exposureInfo : `lsst.pipe.base.Struct`
1065 Struct containing general properties for the visits:
1066 ``visits`` : `list`
1067 List of visit names.
1068 ``detectors`` : `list`
1069 List of all detectors in any visit.
1070 ``ras`` : `list` [`float`]
1071 List of boresight RAs for each visit.
1072 ``decs`` : `list` [`float`]
1073 List of borseight Decs for each visit.
1074 ``medianEpoch`` : float
1075 Median epoch of all visits in decimal-year format.
1076 exposuresHelper : `wcsfit.ExposuresHelper`
1077 Object containing information about the input visits.
1078 extensionInfo : `lsst.pipe.base.Struct`
1079 Struct containing properties for each extension (visit/detector):
1080 ``visit`` : `np.ndarray`
1081 Name of the visit for this extension.
1082 ``detector`` : `np.ndarray`
1083 Name of the detector for this extension.
1084 ``visitIndex` : `np.ndarray` [`int`]
1085 Index of visit for this extension.
1086 ``detectorIndex`` : `np.ndarray` [`int`]
1087 Index of the detector for this extension.
1088 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1089 Initial WCS for this extension.
1090 ``extensionType`` : `np.ndarray` [`str`]
1091 "SCIENCE" or "REFERENCE".
1092 instruments : `list` [`wcsfit.Instrument`]
1093 List of instrument objects used for the input visits.
1094 """
1095 exposureNames = []
1096 ras = []
1097 decs = []
1098 visits = []
1099 detectors = []
1100 airmasses = []
1101 exposureTimes = []
1102 mjds = []
1103 observatories = []
1104 wcss = []
1105 fieldNumbers = []
1107 extensionType = []
1108 extensionVisitIndices = []
1109 extensionDetectorIndices = []
1110 extensionVisits = []
1111 extensionDetectors = []
1113 instruments, instrumentNumbers = _get_instruments(inputVisitSummaries)
1115 # Get information for all the science visits
1116 for v, visitSummary in enumerate(inputVisitSummaries):
1117 visitInfo = visitSummary[0].getVisitInfo()
1118 visit = visitSummary[0]["visit"]
1119 visits.append(visit)
1120 exposureNames.append(str(visit))
1121 raDec = visitInfo.getBoresightRaDec()
1122 ras.append(raDec.getRa().asRadians())
1123 decs.append(raDec.getDec().asRadians())
1124 if fieldRegions is not None:
1125 inField = [r for r, region in fieldRegions.items() if region.contains(raDec.getVector())]
1126 if len(inField) != 1:
1127 raise RuntimeError(
1128 f"Visit should be in one and only one field, but {visit} is contained "
1129 f"in {len(inField)} fields."
1130 )
1131 fieldNumbers.append(inField[0])
1132 else:
1133 fieldNumbers.append(fieldNumber)
1134 airmasses.append(visitInfo.getBoresightAirmass())
1135 exposureTimes.append(visitInfo.getExposureTime())
1136 obsDate = visitInfo.getDate()
1137 obsMJD = obsDate.get(obsDate.MJD)
1138 mjds.append(obsMJD)
1139 # Get the observatory ICRS position for use in fitting parallax
1140 obsLon = visitInfo.observatory.getLongitude().asDegrees()
1141 obsLat = visitInfo.observatory.getLatitude().asDegrees()
1142 obsElev = visitInfo.observatory.getElevation()
1143 earthLocation = astropy.coordinates.EarthLocation.from_geodetic(obsLon, obsLat, obsElev)
1144 observatory_gcrs = earthLocation.get_gcrs(astropy.time.Time(obsMJD, format="mjd"))
1145 observatory_icrs = observatory_gcrs.transform_to(astropy.coordinates.ICRS())
1146 # We want the position in AU in Cartesian coordinates
1147 observatories.append(observatory_icrs.cartesian.xyz.to(u.AU).value)
1149 validDetectors = [row for row in visitSummary if row.wcs is not None]
1150 nDetectorFraction = len(validDetectors) / len(visitSummary)
1151 if nDetectorFraction < self.config.minDetectorFraction:
1152 self.log.warning(
1153 "Visit %d has only %d detectors with valid WCSs (%s of total) and will be dropped.",
1154 visit,
1155 len(validDetectors),
1156 nDetectorFraction,
1157 )
1158 continue
1160 for row in visitSummary:
1161 detector = row["id"]
1163 wcs = row.getWcs()
1164 if wcs is None:
1165 self.log.warning(
1166 "WCS is None for visit %d, detector %d: this extension (visit/detector) will be "
1167 "dropped.",
1168 visit,
1169 detector,
1170 )
1171 continue
1172 else:
1173 wcsRA = wcs.getSkyOrigin().getRa().asRadians()
1174 wcsDec = wcs.getSkyOrigin().getDec().asRadians()
1175 tangentPoint = wcsfit.Gnomonic(wcsRA, wcsDec)
1176 mapping = wcs.getFrameDict().getMapping("PIXELS", "IWC")
1177 gbdes_wcs = wcsfit.Wcs(wcsfit.ASTMap(mapping), tangentPoint)
1178 wcss.append(gbdes_wcs)
1180 if detector not in detectors:
1181 detectors.append(detector)
1182 if not instruments[instrumentNumbers[v]].hasDevice(str(detector)):
1183 detectorBounds = wcsfit.Bounds(
1184 row["bbox_min_x"], row["bbox_max_x"], row["bbox_min_y"], row["bbox_max_y"]
1185 )
1186 instruments[instrumentNumbers[v]].addDevice(str(detector), detectorBounds)
1188 detectorIndex = np.flatnonzero(detector == np.array(detectors))[0]
1189 extensionVisitIndices.append(v)
1190 extensionDetectorIndices.append(detectorIndex)
1191 extensionVisits.append(visit)
1192 extensionDetectors.append(detector)
1193 extensionType.append("SCIENCE")
1195 if len(wcss) == 0:
1196 raise pipeBase.NoWorkFound("No input extensions have valid WCSs.")
1198 # Set the reference epoch to be the median of the science visits.
1199 # The reference catalog will be shifted to this date.
1200 if self.config.setRefEpoch is None:
1201 medianMJD = np.median(mjds)
1202 self.log.info(f"Ref epoch set to median: {medianMJD}")
1203 else:
1204 medianMJD = self.config.setRefEpoch
1205 self.log.info(f"Ref epoch set by user: {medianMJD}")
1206 medianEpoch = astropy.time.Time(medianMJD, format="mjd").decimalyear
1208 # Add information for the reference catalog. Most of the values are
1209 # not used. There needs to be a separate catalog for each field.
1210 if fieldRegions is None:
1211 fieldRegions = {0: None}
1212 for f in fieldRegions:
1213 exposureNames.append("REFERENCE")
1214 # Make the "visit" number the field * -1 to disambiguate it from
1215 # any potential visit number:
1216 visits.append(-1 * f)
1217 fieldNumbers.append(f)
1218 if self.config.fitProperMotion:
1219 instrumentNumbers.append(-2)
1220 else:
1221 instrumentNumbers.append(-1)
1222 ras.append(0.0)
1223 decs.append(0.0)
1224 airmasses.append(0.0)
1225 exposureTimes.append(0)
1226 mjds.append((refEpoch if (refEpoch is not None) else medianMJD))
1227 observatories.append(np.array([0, 0, 0]))
1228 identity = wcsfit.IdentityMap()
1229 icrs = wcsfit.SphericalICRS()
1230 refWcs = wcsfit.Wcs(identity, icrs, "Identity", np.pi / 180.0)
1231 wcss.append(refWcs)
1233 extensionVisitIndices.append(len(exposureNames) - 1)
1234 extensionDetectorIndices.append(-1) # REFERENCE device must be -1
1235 extensionVisits.append(-1 * f)
1236 extensionDetectors.append(-1)
1237 extensionType.append("REFERENCE")
1239 # Make a table of information to use elsewhere in the class
1240 extensionInfo = pipeBase.Struct(
1241 visit=np.array(extensionVisits),
1242 detector=np.array(extensionDetectors),
1243 visitIndex=np.array(extensionVisitIndices),
1244 detectorIndex=np.array(extensionDetectorIndices),
1245 wcs=np.array(wcss),
1246 extensionType=np.array(extensionType),
1247 )
1249 # Make the exposureHelper object to use in the fitting routines
1250 exposuresHelper = wcsfit.ExposuresHelper(
1251 exposureNames,
1252 fieldNumbers,
1253 instrumentNumbers,
1254 ras,
1255 decs,
1256 airmasses,
1257 exposureTimes,
1258 mjds,
1259 observatories,
1260 )
1262 exposureInfo = pipeBase.Struct(
1263 visits=visits, detectors=detectors, ras=ras, decs=decs, medianEpoch=medianEpoch
1264 )
1266 return exposureInfo, exposuresHelper, extensionInfo, instruments
1268 def _load_refcat(
1269 self,
1270 refObjectLoader,
1271 extensionInfo,
1272 epoch=None,
1273 fieldIndex=0,
1274 associations=None,
1275 center=None,
1276 radius=None,
1277 region=None,
1278 ):
1279 """Load the reference catalog and add reference objects to the
1280 `wcsfit.FoFClass` object.
1282 Parameters
1283 ----------
1284 refObjectLoader :
1285 `lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader`
1286 Object set up to load reference catalog objects.
1287 extensionInfo : `lsst.pipe.base.Struct`
1288 Struct containing properties for each extension (visit/detector).
1289 ``visit`` : `np.ndarray`
1290 Name of the visit for this extension.
1291 ``detector`` : `np.ndarray`
1292 Name of the detector for this extension.
1293 ``visitIndex` : `np.ndarray` [`int`]
1294 Index of visit for this extension.
1295 ``detectorIndex`` : `np.ndarray` [`int`]
1296 Index of the detector for this extension.
1297 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1298 Initial WCS for this extension.
1299 ``extensionType`` : `np.ndarray` [`str`]
1300 "SCIENCE" or "REFERENCE".
1301 epoch : `float`, optional
1302 MJD to which to correct the object positions.
1303 fieldIndex : `int`, optional
1304 Index of the field. Should be zero if all the data is fit together.
1305 associations : `wcsfit.FoFClass`, optional
1306 Object to which to add the catalog of reference objects.
1307 center : `lsst.geom.SpherePoint`, optional
1308 Center of the circle in which to load reference objects. Ignored if
1309 `region` is set. If used, `radius` must also be set.
1310 radius : `lsst.sphgeom._sphgeom.Angle`, optional
1311 Radius of the circle in which to load reference objects. Ignored if
1312 `region` is set. If used, `center` must also be set.
1313 region : `lsst.sphgeom.ConvexPolygon`, optional
1314 Region in which to load reference objects.
1316 Returns
1317 -------
1318 refObjects : `dict`
1319 Position and error information of reference objects.
1320 refCovariance : `list` [`float`]
1321 Flattened output covariance matrix.
1322 """
1323 if self.config.applyRefCatProperMotion:
1324 formattedEpoch = astropy.time.Time(epoch, format="mjd")
1325 else:
1326 formattedEpoch = None
1328 if region is not None:
1329 skyRegion = refObjectLoader.loadRegion(region, self.config.referenceFilter, epoch=formattedEpoch)
1330 elif (center is not None) and (radius is not None):
1331 skyRegion = refObjectLoader.loadSkyCircle(
1332 center, radius, self.config.referenceFilter, epoch=formattedEpoch
1333 )
1334 else:
1335 raise RuntimeError("Either `region` or `center` and `radius` must be set.")
1337 selected = self.referenceSelector.run(skyRegion.refCat)
1338 # Need memory contiguity to get reference filters as a vector.
1339 if not selected.sourceCat.isContiguous():
1340 refCat = selected.sourceCat.copy(deep=True)
1341 else:
1342 refCat = selected.sourceCat
1343 refCat = refCat.asAstropy()
1345 # In Gaia DR3, missing values are denoted by NaNs.
1346 finiteInd = np.isfinite(refCat["coord_ra"]) & np.isfinite(refCat["coord_dec"])
1347 refCat = refCat[finiteInd]
1349 if self.config.excludeNonPMObjects and self.config.applyRefCatProperMotion:
1350 # Gaia DR2 has zeros for missing data, while Gaia DR3 has NaNs:
1351 hasPM = (
1352 (refCat["pm_raErr"] != 0) & np.isfinite(refCat["pm_raErr"]) & np.isfinite(refCat["pm_decErr"])
1353 )
1354 refCat = refCat[hasPM]
1356 ra = (refCat["coord_ra"]).to(u.degree).to_value().tolist()
1357 dec = (refCat["coord_dec"]).to(u.degree).to_value().tolist()
1358 raCov = ((refCat["coord_raErr"]).to(u.degree).to_value() ** 2).tolist()
1359 decCov = ((refCat["coord_decErr"]).to(u.degree).to_value() ** 2).tolist()
1361 # Get refcat version from refcat metadata
1362 refCatMetadata = refObjectLoader.refCats[0].get().getMetadata()
1364 # DM-47181: Added this to work for LSSTComCam with
1365 # the_monster_20240904 catalog that has not this key.
1366 if "REFCAT_FORMAT_VERSION" not in refCatMetadata:
1367 refCatVersion = 2
1368 raDecCov = (refCat["coord_ra_coord_dec_Cov"]).to(u.degree**2).to_value().tolist()
1369 else:
1370 refCatVersion = refCatMetadata["REFCAT_FORMAT_VERSION"]
1371 if refCatVersion == 2:
1372 raDecCov = (refCat["coord_ra_coord_dec_Cov"]).to(u.degree**2).to_value().tolist()
1373 else:
1374 raDecCov = np.zeros(len(ra))
1376 refObjects = {"ra": ra, "dec": dec, "raCov": raCov, "decCov": decCov, "raDecCov": raDecCov}
1377 refCovariance = []
1379 if self.config.fitProperMotion:
1380 raPM = (refCat["pm_ra"]).to(u.marcsec).to_value().tolist()
1381 decPM = (refCat["pm_dec"]).to(u.marcsec).to_value().tolist()
1382 parallax = (refCat["parallax"]).to(u.marcsec).to_value().tolist()
1383 cov = _make_ref_covariance_matrix(refCat, version=refCatVersion)
1384 pmDict = {"raPM": raPM, "decPM": decPM, "parallax": parallax}
1385 refObjects.update(pmDict)
1386 refCovariance = cov
1388 if associations is not None:
1389 extensionIndex = np.flatnonzero(extensionInfo.extensionType == "REFERENCE")[0]
1390 visitIndex = extensionInfo.visitIndex[extensionIndex]
1391 detectorIndex = extensionInfo.detectorIndex[extensionIndex]
1392 instrumentIndex = -1 # -1 indicates the reference catalog
1393 refWcs = extensionInfo.wcs[extensionIndex]
1395 associations.addCatalog(
1396 refWcs,
1397 "STELLAR",
1398 visitIndex,
1399 fieldIndex,
1400 instrumentIndex,
1401 detectorIndex,
1402 extensionIndex,
1403 np.ones(len(refCat), dtype=bool),
1404 ra,
1405 dec,
1406 np.arange(len(ra)),
1407 )
1409 return refObjects, refCovariance
1411 @staticmethod
1412 def _find_extension_index(extensionInfo, visit, detector):
1413 """Find the index for a given extension from its visit and detector
1414 number.
1416 If no match is found, None is returned.
1418 Parameters
1419 ----------
1420 extensionInfo : `lsst.pipe.base.Struct`
1421 Struct containing properties for each extension.
1422 visit : `int`
1423 Visit number
1424 detector : `int`
1425 Detector number
1427 Returns
1428 -------
1429 extensionIndex : `int` or None
1430 Index of this extension
1431 """
1432 findExtension = np.flatnonzero((extensionInfo.visit == visit) & (extensionInfo.detector == detector))
1433 if len(findExtension) == 0:
1434 extensionIndex = None
1435 else:
1436 extensionIndex = findExtension[0]
1437 return extensionIndex
1439 def _load_catalogs_and_associate(
1440 self, associations, inputCatalogRefs, extensionInfo, fieldIndex=0, instrumentIndex=0
1441 ):
1442 """Load the science catalogs and add the sources to the associator
1443 class `wcsfit.FoFClass`, associating them into matches as you go.
1445 Parameters
1446 ----------
1447 associations : `wcsfit.FoFClass`
1448 Object to which to add the catalog of source and which performs
1449 the source association.
1450 inputCatalogRefs : `list`
1451 List of DeferredDatasetHandles pointing to visit-level source
1452 tables.
1453 extensionInfo : `lsst.pipe.base.Struct`
1454 Struct containing properties for each extension (visit/detector).
1455 ``visit`` : `np.ndarray`
1456 Name of the visit for this extension.
1457 ``detector`` : `np.ndarray`
1458 Name of the detector for this extension.
1459 ``visitIndex` : `np.ndarray` [`int`]
1460 Index of visit for this extension.
1461 ``detectorIndex`` : `np.ndarray` [`int`]
1462 Index of the detector for this extension.
1463 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1464 Initial WCS for this extension.
1465 ``extensionType`` : `np.ndarray` [`str`]
1466 "SCIENCE" or "REFERENCE".
1467 fieldIndex : `int`
1468 Index of the field for these catalogs. Should be zero assuming all
1469 data is being fit together.
1470 instrumentIndex : `int`
1471 Index of the instrument for these catalogs. Should be zero
1472 assuming all data comes from the same instrument.
1474 Returns
1475 -------
1476 sourceIndices : `list`
1477 List of boolean arrays used to select sources.
1478 columns : `list` [`str`]
1479 List of columns needed from source tables.
1480 """
1481 columns = [
1482 "detector",
1483 "sourceId",
1484 "x",
1485 "xErr",
1486 "y",
1487 "yErr",
1488 "ixx",
1489 "iyy",
1490 "ixy",
1491 f"{self.config.sourceFluxType}_instFlux",
1492 f"{self.config.sourceFluxType}_instFluxErr",
1493 ]
1494 if self.sourceSelector.config.doFlags:
1495 columns.extend(self.sourceSelector.config.flags.bad)
1496 if self.sourceSelector.config.doUnresolved:
1497 columns.append(self.sourceSelector.config.unresolved.name)
1498 if self.sourceSelector.config.doIsolated:
1499 columns.append(self.sourceSelector.config.isolated.parentName)
1500 columns.append(self.sourceSelector.config.isolated.nChildName)
1501 if self.sourceSelector.config.doRequirePrimary:
1502 columns.append(self.sourceSelector.config.requirePrimary.primaryColName)
1504 sourceIndices = [None] * len(extensionInfo.visit)
1505 for inputCatalogRef in inputCatalogRefs:
1506 visit = inputCatalogRef.dataId["visit"]
1507 inputCatalog = inputCatalogRef.get(parameters={"columns": columns})
1508 # Get a sorted array of detector names
1509 detectors = np.unique(inputCatalog["detector"])
1511 for detector in detectors:
1512 detectorSources = inputCatalog[inputCatalog["detector"] == detector]
1513 xCov = detectorSources["xErr"] ** 2
1514 yCov = detectorSources["yErr"] ** 2
1515 xyCov = (
1516 detectorSources["ixy"] * (xCov + yCov) / (detectorSources["ixx"] + detectorSources["iyy"])
1517 )
1518 # Remove sources with bad shape measurements
1519 goodShapes = xyCov**2 <= (xCov * yCov)
1520 selected = self.sourceSelector.run(detectorSources)
1521 goodInds = selected.selected & goodShapes
1522 isStar = np.ones(goodInds.sum())
1523 extensionIndex = self._find_extension_index(extensionInfo, visit, detector)
1524 if extensionIndex is None:
1525 # This extension does not have information necessary for
1526 # fit. Skip the detections from this detector for this
1527 # visit.
1528 continue
1529 detectorIndex = extensionInfo.detectorIndex[extensionIndex]
1530 visitIndex = extensionInfo.visitIndex[extensionIndex]
1532 sourceIndices[extensionIndex] = goodInds
1534 wcs = extensionInfo.wcs[extensionIndex]
1535 associations.reprojectWCS(wcs, fieldIndex)
1537 associations.addCatalog(
1538 wcs,
1539 "STELLAR",
1540 visitIndex,
1541 fieldIndex,
1542 instrumentIndex,
1543 detectorIndex,
1544 extensionIndex,
1545 isStar,
1546 detectorSources[goodInds]["x"].to_list(),
1547 detectorSources[goodInds]["y"].to_list(),
1548 np.arange(goodInds.sum()),
1549 )
1551 associations.sortMatches(
1552 fieldIndex, minMatches=self.config.minMatches, allowSelfMatches=self.config.allowSelfMatches
1553 )
1555 return sourceIndices, columns
1557 def _check_degeneracies(self, associations, extensionInfo):
1558 """Check that the minimum number of visits and sources needed to
1559 constrain the model are present.
1561 This does not guarantee that the Hessian matrix of the chi-square,
1562 which is used to fit the model, will be positive-definite, but if the
1563 checks here do not pass, the matrix is certain to not be
1564 positive-definite and the model cannot be fit.
1566 Parameters
1567 ----------
1568 associations : `wcsfit.FoFClass`
1569 Object holding the source association information.
1570 extensionInfo : `lsst.pipe.base.Struct`
1571 Struct containing properties for each extension (visit/detector):
1572 ``visit`` : `np.ndarray`
1573 Name of the visit for this extension.
1574 ``detector`` : `np.ndarray`
1575 Name of the detector for this extension.
1576 ``visitIndex` : `np.ndarray` [`int`]
1577 Index of visit for this extension.
1578 ``detectorIndex`` : `np.ndarray` [`int`]
1579 Index of the detector for this extension.
1580 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1581 Initial WCS for this extension.
1582 ``extensionType`` : `np.ndarray` [`str`]
1583 "SCIENCE" or "REFERENCE".
1584 """
1585 # As a baseline, need to have more stars per detector than per-detector
1586 # parameters, and more stars per visit than per-visit parameters.
1587 whichExtension = np.array(associations.extn)
1588 whichDetector = np.zeros(len(whichExtension))
1589 whichVisit = np.zeros(len(whichExtension))
1591 for extension, (detector, visit) in enumerate(zip(extensionInfo.detector, extensionInfo.visit)):
1592 ex_ind = whichExtension == extension
1593 whichDetector[ex_ind] = detector
1594 whichVisit[ex_ind] = visit
1596 if (not self.config.useInputCameraModel) and ("BAND/DEVICE/poly" in self.config.deviceModel):
1597 nCoeffDetectorModel = _nCoeffsFromDegree(self.config.devicePolyOrder)
1598 unconstrainedDetectors = []
1599 for detector in np.unique(extensionInfo.detector):
1600 numSources = (whichDetector == detector).sum()
1601 if numSources < nCoeffDetectorModel:
1602 unconstrainedDetectors.append(str(detector))
1604 if unconstrainedDetectors:
1605 raise pipeBase.NoWorkFound(
1606 "The model is not constrained. The following detectors do not have enough "
1607 f"sources ({nCoeffDetectorModel} required): ",
1608 ", ".join(unconstrainedDetectors),
1609 )
1611 def make_yaml(self, inputVisitSummary, inputFile=None, inputCameraModel=None):
1612 """Make a YAML-type object that describes the parameters of the fit
1613 model.
1615 Parameters
1616 ----------
1617 inputVisitSummary : `lsst.afw.table.ExposureCatalog`
1618 Catalog with per-detector summary information.
1619 inputFile : `str`
1620 Path to a file that contains a basic model.
1621 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
1622 Parameters to use for the device part of the model.
1624 Returns
1625 -------
1626 inputYaml : `wcsfit.YAMLCollector`
1627 YAML object containing the model description.
1628 inputDict : `dict` [`str`, `str`]
1629 Dictionary containing the model description.
1630 """
1631 if inputFile is not None:
1632 inputYaml = wcsfit.YAMLCollector(inputFile, "PixelMapCollection")
1633 else:
1634 inputYaml = wcsfit.YAMLCollector("", "PixelMapCollection")
1635 inputDict = {}
1636 modelComponents = ["BAND/DEVICE", "EXPOSURE"]
1637 baseMap = {"Type": "Composite", "Elements": modelComponents}
1638 inputDict["EXPOSURE/DEVICE/base"] = baseMap
1640 xMin = str(inputVisitSummary["bbox_min_x"].min())
1641 xMax = str(inputVisitSummary["bbox_max_x"].max())
1642 yMin = str(inputVisitSummary["bbox_min_y"].min())
1643 yMax = str(inputVisitSummary["bbox_max_y"].max())
1645 deviceModel = {"Type": "Composite", "Elements": self.config.deviceModel.list()}
1646 inputDict["BAND/DEVICE"] = deviceModel
1647 for component in self.config.deviceModel:
1648 mapping = component.split("/")[-1]
1649 if mapping == "poly":
1650 componentDict = {
1651 "Type": "Poly",
1652 "XPoly": {"OrderX": self.config.devicePolyOrder, "SumOrder": True},
1653 "YPoly": {"OrderX": self.config.devicePolyOrder, "SumOrder": True},
1654 "XMin": xMin,
1655 "XMax": xMax,
1656 "YMin": yMin,
1657 "YMax": yMax,
1658 }
1659 elif mapping == "identity":
1660 componentDict = {"Type": "Identity"}
1662 inputDict[component] = componentDict
1664 if (inputCameraModel is not None) and self.config.useInputCameraModel:
1665 # This assumes that the input camera model is a 'poly' model
1666 nCoeffs = _nCoeffsFromDegree(self.config.devicePolyOrder)
1667 for key, coeffs in inputCameraModel.items():
1668 if len(coeffs) != nCoeffs * 2:
1669 raise RuntimeError(
1670 "Input camera model polynomial order does not match the devicePolyOrder"
1671 )
1672 mapDict = {
1673 "Type": "Poly",
1674 "XPoly": {
1675 "OrderX": self.config.devicePolyOrder,
1676 "SumOrder": True,
1677 "Coefficients": coeffs[:nCoeffs].tolist(),
1678 },
1679 "YPoly": {
1680 "OrderX": self.config.devicePolyOrder,
1681 "SumOrder": True,
1682 "Coefficients": coeffs[nCoeffs:].tolist(),
1683 },
1684 "XMin": xMin,
1685 "XMax": xMax,
1686 "YMin": yMin,
1687 "YMax": yMax,
1688 }
1689 inputDict[key] = mapDict
1691 exposureModelComponents = self.config.exposureModel.list()
1692 if self.config.useColor:
1693 exposureModelComponents.append("EXPOSURE/dcr")
1694 exposureModel = {"Type": "Composite", "Elements": exposureModelComponents}
1695 inputDict["EXPOSURE"] = exposureModel
1696 for component in exposureModelComponents:
1697 mapping = component.split("/")[-1]
1698 if mapping == "poly":
1699 componentDict = {
1700 "Type": "Poly",
1701 "XPoly": {"OrderX": self.config.exposurePolyOrder, "SumOrder": "true"},
1702 "YPoly": {"OrderX": self.config.exposurePolyOrder, "SumOrder": "true"},
1703 }
1704 elif mapping == "identity":
1705 componentDict = {"Type": "Identity"}
1706 elif mapping == "dcr":
1707 componentDict = {
1708 "Type": "Color",
1709 "Reference": self.config.referenceColor,
1710 "Function": {"Type": "Constant"},
1711 }
1713 inputDict[component] = componentDict
1715 inputYaml.addInput(yaml.dump(inputDict))
1716 inputYaml.addInput("Identity:\n Type: Identity\n")
1718 return inputYaml, inputDict
1720 def _add_objects(self, wcsf, inputCatalogRefs, sourceIndices, extensionInfo, columns):
1721 """Add science sources to the wcsfit.WCSFit object.
1723 Parameters
1724 ----------
1725 wcsf : `wcsfit.WCSFit`
1726 WCS-fitting object.
1727 inputCatalogRefs : `list`
1728 List of DeferredDatasetHandles pointing to visit-level source
1729 tables.
1730 sourceIndices : `list`
1731 List of boolean arrays used to select sources.
1732 extensionInfo : `lsst.pipe.base.Struct`
1733 Struct containing properties for each extension (visit/detector):
1734 ``visit`` : `np.ndarray`
1735 Name of the visit for this extension.
1736 ``detector`` : `np.ndarray`
1737 Name of the detector for this extension.
1738 ``visitIndex` : `np.ndarray` [`int`]
1739 Index of visit for this extension.
1740 ``detectorIndex`` : `np.ndarray` [`int`]
1741 Index of the detector for this extension.
1742 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1743 Initial WCS for this extension.
1744 ``extensionType`` : `np.ndarray` [`str`]
1745 "SCIENCE" or "REFERENCE".
1746 columns : `list` [`str`]
1747 List of columns needed from source tables.
1748 """
1749 for inputCatalogRef in inputCatalogRefs:
1750 visit = inputCatalogRef.dataId["visit"]
1751 inputCatalog = inputCatalogRef.get(parameters={"columns": columns})
1752 detectors = np.unique(inputCatalog["detector"])
1754 for detector in detectors:
1755 detectorSources = inputCatalog[inputCatalog["detector"] == detector]
1757 extensionIndex = self._find_extension_index(extensionInfo, visit, detector)
1758 if extensionIndex is None:
1759 # This extension does not have information necessary for
1760 # fit. Skip the detections from this detector for this
1761 # visit.
1762 continue
1764 sourceCat = detectorSources[sourceIndices[extensionIndex]]
1766 xCov = sourceCat["xErr"] ** 2
1767 yCov = sourceCat["yErr"] ** 2
1768 xyCov = sourceCat["ixy"] * (xCov + yCov) / (sourceCat["ixx"] + sourceCat["iyy"])
1769 # TODO: add correct xyErr if DM-7101 is ever done.
1771 d = {
1772 "x": sourceCat["x"].to_numpy(),
1773 "y": sourceCat["y"].to_numpy(),
1774 "xCov": xCov.to_numpy(),
1775 "yCov": yCov.to_numpy(),
1776 "xyCov": xyCov.to_numpy(),
1777 }
1779 wcsf.setObjects(
1780 extensionIndex,
1781 d,
1782 "x",
1783 "y",
1784 ["xCov", "yCov", "xyCov"],
1785 defaultColor=self.config.referenceColor,
1786 )
1788 def _add_ref_objects(self, wcsf, refObjects, refCovariance, extensionInfo, fieldIndex=0):
1789 """Add reference sources to the wcsfit.WCSFit object.
1791 Parameters
1792 ----------
1793 wcsf : `wcsfit.WCSFit`
1794 WCS-fitting object.
1795 refObjects : `dict`
1796 Position and error information of reference objects.
1797 refCovariance : `list` [`float`]
1798 Flattened output covariance matrix.
1799 extensionInfo : `lsst.pipe.base.Struct`
1800 Struct containing properties for each extension (visit/detector):
1801 ``visit`` : `np.ndarray`
1802 Name of the visit for this extension.
1803 ``detector`` : `np.ndarray`
1804 Name of the detector for this extension.
1805 ``visitIndex` : `np.ndarray` [`int`]
1806 Index of visit for this extension.
1807 ``detectorIndex`` : `np.ndarray` [`int`]
1808 Index of the detector for this extension.
1809 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1810 Initial WCS for this extension.
1811 ``extensionType`` : `np.ndarray` [`str`]
1812 "SCIENCE" or "REFERENCE".
1813 fieldIndex : `int`, optional
1814 Index of the field to which these sources correspond.
1815 """
1816 extensionIndex = np.flatnonzero(
1817 (extensionInfo.extensionType == "REFERENCE") & (extensionInfo.visit == fieldIndex)
1818 )[0]
1819 if self.config.fitProperMotion:
1820 wcsf.setObjects(
1821 extensionIndex,
1822 refObjects,
1823 "ra",
1824 "dec",
1825 ["raCov", "decCov", "raDecCov"],
1826 pmDecKey="decPM",
1827 pmRaKey="raPM",
1828 parallaxKey="parallax",
1829 pmCovKey="fullCov",
1830 pmCov=refCovariance,
1831 )
1832 else:
1833 wcsf.setObjects(extensionIndex, refObjects, "ra", "dec", ["raCov", "decCov", "raDecCov"])
1835 def _add_color_objects(self, wcsf, colorCatalog):
1836 """Associate input matches with objects in color catalog and set their
1837 color value.
1839 Parameters
1840 ----------
1841 wcsf : `wcsfit.WCSFit`
1842 WCSFit object, assumed to have fit model.
1843 colorCatalog : `lsst.afw.table.SimpleCatalog`
1844 Catalog containing object coordinates and magnitudes.
1845 """
1847 # Get current best position for matches
1848 starCat = wcsf.getStarCatalog()
1850 # TODO: DM-45650, update how the colors are read in here.
1851 catalogBands = colorCatalog.metadata.getArray("BANDS")
1852 colorInd1 = catalogBands.index(self.config.color[0])
1853 colorInd2 = catalogBands.index(self.config.color[1])
1854 colors = colorCatalog["mag_std_noabs"][:, colorInd1] - colorCatalog["mag_std_noabs"][:, colorInd2]
1855 goodInd = (colorCatalog["mag_std_noabs"][:, colorInd1] != 99.0) & (
1856 colorCatalog["mag_std_noabs"][:, colorInd2] != 99.0
1857 )
1859 with Matcher(np.array(starCat["starX"]), np.array(starCat["starY"])) as matcher:
1860 idx, idx_starCat, idx_colorCat, d = matcher.query_radius(
1861 (colorCatalog[goodInd]["coord_ra"] * u.radian).to(u.degree).value,
1862 (colorCatalog[goodInd]["coord_dec"] * u.radian).to(u.degree).value,
1863 self.config.matchRadius / 3600.0,
1864 return_indices=True,
1865 )
1867 matchesWithColor = starCat["starMatchID"][idx_starCat]
1868 matchColors = np.ones(len(matchesWithColor)) * self.config.referenceColor
1869 matchColors = colors[goodInd][idx_colorCat]
1870 wcsf.setColors(matchesWithColor, matchColors)
1872 def _make_afw_wcs(self, mapDict, centerRA, centerDec, doNormalizePixels=False, xScale=1, yScale=1):
1873 """Make an `lsst.afw.geom.SkyWcs` from a dictionary of mappings.
1875 Parameters
1876 ----------
1877 mapDict : `dict`
1878 Dictionary of mapping parameters.
1879 centerRA : `lsst.geom.Angle`
1880 RA of the tangent point.
1881 centerDec : `lsst.geom.Angle`
1882 Declination of the tangent point.
1883 doNormalizePixels : `bool`
1884 Whether to normalize pixels so that range is [-1,1].
1885 xScale : `float`
1886 Factor by which to normalize x-dimension. Corresponds to width of
1887 detector.
1888 yScale : `float`
1889 Factor by which to normalize y-dimension. Corresponds to height of
1890 detector.
1892 Returns
1893 -------
1894 outWCS : `lsst.afw.geom.SkyWcs`
1895 WCS constructed from the input mappings
1896 """
1897 # Set up pixel frames
1898 pixelFrame = astshim.Frame(2, "Domain=PIXELS")
1899 normedPixelFrame = astshim.Frame(2, "Domain=NORMEDPIXELS")
1901 if doNormalizePixels:
1902 # Pixels will need to be rescaled before going into the mappings
1903 normCoefficients = [-1.0, 2.0 / xScale, 0, -1.0, 0, 2.0 / yScale]
1904 normMap = _convert_to_ast_polymap_coefficients(normCoefficients)
1905 else:
1906 normMap = astshim.UnitMap(2)
1908 # All of the detectors for one visit map to the same tangent plane
1909 tangentPoint = lsst.geom.SpherePoint(centerRA, centerDec)
1910 cdMatrix = afwgeom.makeCdMatrix(1.0 * lsst.geom.degrees, 0 * lsst.geom.degrees, True)
1911 iwcToSkyWcs = afwgeom.makeSkyWcs(lsst.geom.Point2D(0, 0), tangentPoint, cdMatrix)
1912 iwcToSkyMap = iwcToSkyWcs.getFrameDict().getMapping("PIXELS", "SKY")
1913 skyFrame = iwcToSkyWcs.getFrameDict().getFrame("SKY")
1915 frameDict = astshim.FrameDict(pixelFrame)
1916 frameDict.addFrame("PIXELS", normMap, normedPixelFrame)
1918 currentFrameName = "NORMEDPIXELS"
1920 # Dictionary values are ordered according to the maps' application.
1921 for m, mapElement in enumerate(mapDict.values()):
1922 mapType = mapElement["Type"]
1924 if mapType == "Poly":
1925 mapCoefficients = mapElement["Coefficients"]
1926 astMap = _convert_to_ast_polymap_coefficients(mapCoefficients)
1927 elif mapType == "Identity":
1928 astMap = astshim.UnitMap(2)
1929 else:
1930 raise ValueError(f"Converting map type {mapType} to WCS is not supported")
1932 if m == len(mapDict) - 1:
1933 newFrameName = "IWC"
1934 else:
1935 newFrameName = "INTERMEDIATE" + str(m)
1936 newFrame = astshim.Frame(2, f"Domain={newFrameName}")
1937 frameDict.addFrame(currentFrameName, astMap, newFrame)
1938 currentFrameName = newFrameName
1939 frameDict.addFrame("IWC", iwcToSkyMap, skyFrame)
1941 outWCS = afwgeom.SkyWcs(frameDict)
1942 return outWCS
1944 def _make_outputs(
1945 self,
1946 wcsf,
1947 visitSummaryTables,
1948 exposureInfo,
1949 mapTemplate,
1950 extensionInfo,
1951 inputCameraModel=None,
1952 inputCamera=None,
1953 ):
1954 """Make a WCS object out of the WCS models.
1956 Parameters
1957 ----------
1958 wcsf : `wcsfit.WCSFit`
1959 WCSFit object, assumed to have fit model.
1960 visitSummaryTables : `list` [`lsst.afw.table.ExposureCatalog`]
1961 Catalogs with per-detector summary information from which to grab
1962 detector information.
1963 extensionInfo : `lsst.pipe.base.Struct`
1964 Struct containing properties for each extension (visit/detector):
1965 ``visit`` : `np.ndarray`
1966 Name of the visit for this extension.
1967 ``detector`` : `np.ndarray`
1968 Name of the detector for this extension.
1969 ``visitIndex` : `np.ndarray` [`int`]
1970 Index of visit for this extension.
1971 ``detectorIndex`` : `np.ndarray` [`int`]
1972 Index of the detector for this extension.
1973 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1974 Initial WCS for this extension.
1975 ``extensionType`` : `np.ndarray` [`str`]
1976 "SCIENCE" or "REFERENCE".
1977 mapTemplate : `dict` [`str`, `str`]
1978 Dictionary containing the model description.
1979 extensionInfo : `lsst.pipe.base.Struct`
1980 Struct containing properties for each extension (visit/detector).
1981 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
1982 Parameters to use for the device part of the model. This must be
1983 provided if an input camera model was used.
1985 Returns
1986 -------
1987 catalogs : `dict` [`str`, `lsst.afw.table.ExposureCatalog`]
1988 Dictionary of `lsst.afw.table.ExposureCatalog` objects with the WCS
1989 set to the WCS fit in wcsf, keyed by the visit name.
1990 cameraParams : `dict` [`str`, `np.ndarray`], optional
1991 Parameters for the device part of the model in the format needed
1992 when used as input for future runs.
1993 colorFits : `dict` [`int`, `np.ndarray`], optional
1994 DCR parameters fit in RA and Dec directions for each visit.
1995 camera : `lsst.afw.cameraGeom.Camera`, optional
1996 Camera object constructed from the per-detector model.
1997 """
1998 # Get the parameters of the fit models
1999 mapParams = wcsf.mapCollection.getParamDict()
2000 cameraParams = {}
2001 if self.config.saveCameraModel:
2002 for element in mapTemplate["BAND/DEVICE"]["Elements"]:
2003 for detector in exposureInfo.detectors:
2004 detectorTemplate = element.replace("DEVICE", str(detector))
2005 detectorTemplate = detectorTemplate.replace("BAND", ".+")
2006 for k, params in mapParams.items():
2007 if re.fullmatch(detectorTemplate, k):
2008 cameraParams[k] = params
2009 if self.config.saveCameraObject:
2010 # Get the average rotation angle of the input visits.
2011 rotations = [
2012 visTable[0].visitInfo.boresightRotAngle.asRadians() for visTable in visitSummaryTables
2013 ]
2014 rotationAngle = np.mean(rotations)
2015 if inputCamera is None:
2016 raise RuntimeError(
2017 "inputCamera must be provided to _make_outputs in order to build output camera."
2018 )
2019 camera = self.buildCamera.run(
2020 mapParams,
2021 mapTemplate,
2022 exposureInfo.detectors,
2023 exposureInfo.visits,
2024 inputCamera,
2025 rotationAngle,
2026 )
2027 else:
2028 camera = None
2029 if self.config.useInputCameraModel:
2030 if inputCameraModel is None:
2031 raise RuntimeError(
2032 "inputCameraModel must be provided to _make_outputs in order to build output WCS."
2033 )
2034 mapParams.update(inputCameraModel)
2036 # Set up the schema for the output catalogs
2037 schema = lsst.afw.table.ExposureTable.makeMinimalSchema()
2038 schema.addField("visit", type="L", doc="Visit number")
2039 schema.addField(
2040 "recoveredWcs",
2041 type="Flag",
2042 doc="Input WCS missing, output recovered from other input visit/detectors.",
2043 )
2045 # Pixels will need to be rescaled before going into the mappings
2046 xscale = int(mapTemplate["BAND/DEVICE/poly"]["XMax"]) - int(mapTemplate["BAND/DEVICE/poly"]["XMin"])
2047 yscale = int(mapTemplate["BAND/DEVICE/poly"]["YMax"]) - int(mapTemplate["BAND/DEVICE/poly"]["YMin"])
2049 catalogs = {}
2050 colorFits = {}
2051 partialOutputs = False
2052 for v, visitSummary in enumerate(visitSummaryTables):
2053 visit = visitSummary[0]["visit"]
2054 if visit not in extensionInfo.visit:
2055 self.log.warning("Visit %d was dropped because no detectors had valid WCSs.", visit)
2056 partialOutputs = True
2057 continue
2059 visitMaps = wcsf.mapCollection.orderAtoms(f"{visit}")
2060 if self.config.useColor:
2061 colorMap = visitMaps.pop(visitMaps.index(f"{visit}/dcr"))
2062 colorFits[visit] = mapParams[colorMap]
2063 visitMap = visitMaps[0]
2064 visitMapType = wcsf.mapCollection.getMapType(visitMap)
2065 if (visitMap not in mapParams) and (visitMapType != "Identity"):
2066 self.log.warning("Visit %d was dropped because of an insufficient amount of data.", visit)
2067 partialOutputs = True
2068 continue
2070 catalog = lsst.afw.table.ExposureCatalog(schema)
2071 catalog.resize(len(exposureInfo.detectors))
2072 catalog["visit"] = visit
2074 for d, detector in enumerate(exposureInfo.detectors):
2075 mapName = f"{visit}/{detector}"
2076 if mapName in wcsf.mapCollection.allMapNames():
2077 mapElements = wcsf.mapCollection.orderAtoms(f"{mapName}/base")
2078 catalog[d]["recoveredWcs"] = False
2079 else:
2080 # This extension was not fit, but the WCS can be recovered
2081 # using the maps fit from sources on other visits but the
2082 # same detector and from sources on other detectors from
2083 # this visit.
2084 genericElements = mapTemplate["EXPOSURE/DEVICE/base"]["Elements"]
2085 mapElements = []
2086 band = visitSummary[0]["physical_filter"]
2088 # Go through the generic map components to build the names
2089 # of the specific maps for this extension.
2090 for component in genericElements:
2091 elements = mapTemplate[component]["Elements"]
2092 for element in elements:
2093 # TODO: DM-42519, gbdes sets the "BAND" to the
2094 # instrument name currently. This will need to be
2095 # disambiguated if we run on multiple bands at
2096 # once.
2097 element = element.replace("BAND", str(band))
2098 element = element.replace("EXPOSURE", str(visit))
2099 element = element.replace("DEVICE", str(detector))
2100 mapElements.append(element)
2101 catalog[d]["recoveredWcs"] = True
2102 mapDict = {}
2103 for m, mapElement in enumerate(mapElements):
2104 mapType = wcsf.mapCollection.getMapType(mapElement)
2105 if mapType == "Color":
2106 # DCR fit should not go into the generic WCS.
2107 continue
2108 mapDict[mapElement] = {"Type": mapType}
2110 if mapType == "Poly":
2111 mapCoefficients = mapParams[mapElement]
2112 mapDict[mapElement]["Coefficients"] = mapCoefficients
2114 # The RA and Dec of the visit are needed for the last step of
2115 # the mapping from the visit tangent plane to RA and Dec
2116 outWCS = self._make_afw_wcs(
2117 mapDict,
2118 exposureInfo.ras[v] * lsst.geom.radians,
2119 exposureInfo.decs[v] * lsst.geom.radians,
2120 doNormalizePixels=True,
2121 xScale=xscale,
2122 yScale=yscale,
2123 )
2125 catalog[d].setId(detector)
2126 catalog[d].setWcs(outWCS)
2127 catalog.sort()
2128 catalogs[visit] = catalog
2129 if self.config.useColor:
2130 colorVisits = np.array(list(colorFits.keys()))
2131 colorRA = np.array([colorFits[vis][0] for vis in colorVisits])
2132 colorDec = np.array([colorFits[vis][1] for vis in colorVisits])
2133 colorFits = {"visit": colorVisits, "raCoefficient": colorRA, "decCoefficient": colorDec}
2135 return catalogs, cameraParams, colorFits, camera, partialOutputs
2137 def _compute_model_params(self, wcsf):
2138 """Get the WCS model parameters and covariance and convert to a
2139 dictionary that will be readable as a pandas dataframe or other table.
2141 Parameters
2142 ----------
2143 wcsf : `wcsfit.WCSFit`
2144 WCSFit object, assumed to have fit model.
2146 Returns
2147 -------
2148 modelParams : `dict`
2149 Parameters and covariance of the best-fit WCS model.
2150 """
2151 modelParamDict = wcsf.mapCollection.getParamDict()
2152 modelCovariance = wcsf.getModelCovariance()
2154 modelParams = {k: [] for k in ["mapName", "coordinate", "parameter", "coefficientNumber"]}
2155 i = 0
2156 for mapName, params in modelParamDict.items():
2157 nCoeffs = len(params)
2158 # There are an equal number of x and y coordinate parameters
2159 nCoordCoeffs = nCoeffs // 2
2160 modelParams["mapName"].extend([mapName] * nCoeffs)
2161 modelParams["coordinate"].extend(["x"] * nCoordCoeffs)
2162 modelParams["coordinate"].extend(["y"] * nCoordCoeffs)
2163 modelParams["parameter"].extend(params)
2164 modelParams["coefficientNumber"].extend(np.arange(nCoordCoeffs))
2165 modelParams["coefficientNumber"].extend(np.arange(nCoordCoeffs))
2167 for p in range(nCoeffs):
2168 if p < nCoordCoeffs:
2169 coord = "x"
2170 else:
2171 coord = "y"
2172 modelParams[f"{mapName}_{coord}_{p}_cov"] = modelCovariance[i]
2173 i += 1
2175 # Convert the dictionary values from lists to numpy arrays.
2176 for key, value in modelParams.items():
2177 modelParams[key] = np.array(value)
2179 return modelParams
2182class GbdesAstrometricMultibandFitConnections(
2183 GbdesAstrometricFitConnections, dimensions=("skymap", "tract", "instrument")
2184):
2185 outputCatalog = pipeBase.connectionTypes.Output(
2186 doc=(
2187 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2188 "plane coordinates and chisq values."
2189 ),
2190 name="gbdesAstrometricMultibandFit_fitStars",
2191 storageClass="ArrowNumpyDict",
2192 dimensions=("instrument", "skymap", "tract"),
2193 )
2194 starCatalog = pipeBase.connectionTypes.Output(
2195 doc=(
2196 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2197 "fitProperMotion is True."
2198 ),
2199 name="gbdesAstrometricMultibandFit_starCatalog",
2200 storageClass="ArrowNumpyDict",
2201 dimensions=("instrument", "skymap", "tract"),
2202 )
2203 modelParams = pipeBase.connectionTypes.Output(
2204 doc="WCS parameters and covariance.",
2205 name="gbdesAstrometricMultibandFit_modelParams",
2206 storageClass="ArrowNumpyDict",
2207 dimensions=("instrument", "skymap", "tract"),
2208 )
2211class GbdesAstrometricMultibandFitConfig(
2212 GbdesAstrometricFitConfig, pipelineConnections=GbdesAstrometricMultibandFitConnections
2213):
2214 pass
2217class GbdesAstrometricMultibandFitTask(GbdesAstrometricFitTask):
2218 """Calibrate the WCS across multiple visits in multiple filters of the same
2219 field using the GBDES package.
2220 """
2222 ConfigClass = GbdesAstrometricMultibandFitConfig
2223 _DefaultName = "gbdesAstrometricMultibandFit"
2226class GbdesGlobalAstrometricFitConnections(
2227 pipeBase.PipelineTaskConnections, dimensions=("instrument", "physical_filter")
2228):
2229 inputVisitSummaries = pipeBase.connectionTypes.Input(
2230 doc=(
2231 "Per-visit consolidated exposure metadata built from calexps. "
2232 "These catalogs use detector id for the id and must be sorted for "
2233 "fast lookups of a detector."
2234 ),
2235 name="visitSummary",
2236 storageClass="ExposureCatalog",
2237 dimensions=("instrument", "visit"),
2238 multiple=True,
2239 )
2240 referenceCatalog = pipeBase.connectionTypes.PrerequisiteInput(
2241 doc="The astrometry reference catalog to match to loaded input catalog sources.",
2242 name="gaia_dr3_20230707",
2243 storageClass="SimpleCatalog",
2244 dimensions=("skypix",),
2245 deferLoad=True,
2246 multiple=True,
2247 )
2248 colorCatalog = pipeBase.connectionTypes.Input(
2249 doc="The catalog of magnitudes to match to input sources.",
2250 name="fgcm_Cycle4_StandardStars",
2251 storageClass="SimpleCatalog",
2252 dimensions=("instrument",),
2253 )
2254 isolatedStarSources = pipeBase.connectionTypes.Input(
2255 doc="Catalog of matched sources.",
2256 name="isolated_star_presources",
2257 storageClass="DataFrame",
2258 dimensions=(
2259 "instrument",
2260 "skymap",
2261 "tract",
2262 ),
2263 multiple=True,
2264 deferLoad=True,
2265 )
2266 isolatedStarCatalogs = pipeBase.connectionTypes.Input(
2267 doc="Catalog of objects corresponding to the isolatedStarSources.",
2268 name="isolated_star_presource_associations",
2269 storageClass="DataFrame",
2270 dimensions=(
2271 "instrument",
2272 "skymap",
2273 "tract",
2274 ),
2275 multiple=True,
2276 deferLoad=True,
2277 )
2278 inputCameraModel = pipeBase.connectionTypes.PrerequisiteInput(
2279 doc="Camera parameters to use for 'device' part of model",
2280 name="gbdesAstrometricFit_cameraModel",
2281 storageClass="ArrowNumpyDict",
2282 dimensions=("instrument", "physical_filter"),
2283 )
2284 inputCamera = pipeBase.connectionTypes.PrerequisiteInput(
2285 doc="Input camera to use when constructing camera from astrometric model.",
2286 name="camera",
2287 storageClass="Camera",
2288 dimensions=("instrument",),
2289 isCalibration=True,
2290 )
2291 outputWcs = pipeBase.connectionTypes.Output(
2292 doc=(
2293 "Per-visit world coordinate systems derived from the fitted model. These catalogs only contain "
2294 "entries for detectors with an output, and use the detector id for the catalog id, sorted on id "
2295 "for fast lookups of a detector."
2296 ),
2297 name="gbdesGlobalAstrometricFitSkyWcsCatalog",
2298 storageClass="ExposureCatalog",
2299 dimensions=("instrument", "visit"),
2300 multiple=True,
2301 )
2302 outputCatalog = pipeBase.connectionTypes.Output(
2303 doc=(
2304 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2305 "plane coordinates and chisq values."
2306 ),
2307 name="gbdesGlobalAstrometricFit_fitStars",
2308 storageClass="ArrowNumpyDict",
2309 dimensions=("instrument", "physical_filter"),
2310 )
2311 starCatalog = pipeBase.connectionTypes.Output(
2312 doc=(
2313 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2314 "fitProperMotion is True."
2315 ),
2316 name="gbdesGlobalAstrometricFit_starCatalog",
2317 storageClass="ArrowNumpyDict",
2318 dimensions=("instrument", "physical_filter"),
2319 )
2320 modelParams = pipeBase.connectionTypes.Output(
2321 doc="WCS parameters and covariance.",
2322 name="gbdesGlobalAstrometricFit_modelParams",
2323 storageClass="ArrowNumpyDict",
2324 dimensions=("instrument", "physical_filter"),
2325 )
2326 outputCameraModel = pipeBase.connectionTypes.Output(
2327 doc="Camera parameters to use for 'device' part of model",
2328 name="gbdesAstrometricFit_cameraModel",
2329 storageClass="ArrowNumpyDict",
2330 dimensions=("instrument", "physical_filter"),
2331 )
2332 dcrCoefficients = pipeBase.connectionTypes.Output(
2333 doc="Per-visit coefficients for DCR correction.",
2334 name="gbdesGlobalAstrometricFit_dcrCoefficients",
2335 storageClass="ArrowNumpyDict",
2336 )
2337 camera = pipeBase.connectionTypes.Output(
2338 doc="Camera object constructed using the per-detector part of the astrometric model",
2339 name="gbdesGlobalAstrometricFitCamera",
2340 storageClass="Camera",
2341 dimensions=("instrument", "physical_filter"),
2342 )
2344 def getSpatialBoundsConnections(self):
2345 return ("inputVisitSummaries",)
2347 def __init__(self, *, config=None):
2348 super().__init__(config=config)
2350 if not self.config.useColor:
2351 self.inputs.remove("colorCatalog")
2352 self.outputs.remove("dcrCoefficients")
2353 if not self.config.saveModelParams:
2354 self.outputs.remove("modelParams")
2355 if not self.config.useInputCameraModel:
2356 self.prerequisiteInputs.remove("inputCameraModel")
2357 if not self.config.saveCameraModel:
2358 self.outputs.remove("outputCameraModel")
2359 if not self.config.saveCameraObject:
2360 self.prerequisiteInputs.remove("inputCamera")
2361 self.outputs.remove("camera")
2364class GbdesGlobalAstrometricFitConfig(
2365 GbdesAstrometricFitConfig, pipelineConnections=GbdesGlobalAstrometricFitConnections
2366):
2367 visitOverlap = pexConfig.Field(
2368 dtype=float,
2369 default=1.0,
2370 doc=(
2371 "The linkage distance threshold above which clustered groups of visits will not be merged "
2372 "together in an agglomerative clustering algorithm. The linkage distance is calculated using the "
2373 "minimum distance between the field-of-view centers of a given visit and all other visits in a "
2374 "group, and is in units of the field-of-view radius. The resulting groups of visits define the "
2375 "fields for the astrometric fit."
2376 ),
2377 )
2380class GbdesGlobalAstrometricFitTask(GbdesAstrometricFitTask):
2381 """Calibrate the WCS across multiple visits and multiple fields using the
2382 GBDES package.
2384 This class assumes that the input visits can be separated into contiguous
2385 groups, for which an individual group covers an area of less than a
2386 hemisphere.
2387 """
2389 ConfigClass = GbdesGlobalAstrometricFitConfig
2390 _DefaultName = "gbdesAstrometricFit"
2392 def runQuantum(self, butlerQC, inputRefs, outputRefs):
2393 # We override runQuantum to set up the refObjLoaders
2394 inputs = butlerQC.get(inputRefs)
2396 instrumentName = butlerQC.quantum.dataId["instrument"]
2398 # Ensure the inputs are in a consistent and deterministic order
2399 inputSumVisits = np.array([inputSum[0]["visit"] for inputSum in inputs["inputVisitSummaries"]])
2400 inputs["inputVisitSummaries"] = [inputs["inputVisitSummaries"][v] for v in inputSumVisits.argsort()]
2401 inputRefHtm7s = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputRefs.referenceCatalog])
2402 inputRefCatRefs = [inputRefs.referenceCatalog[htm7] for htm7 in inputRefHtm7s.argsort()]
2403 inputRefCats = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputs["referenceCatalog"]])
2404 inputs["referenceCatalog"] = [inputs["referenceCatalog"][v] for v in inputRefCats.argsort()]
2405 inputIsolatedStarSourceTracts = np.array(
2406 [isolatedStarSource.dataId["tract"] for isolatedStarSource in inputs["isolatedStarSources"]]
2407 )
2408 inputIsolatedStarCatalogTracts = np.array(
2409 [isolatedStarCatalog.dataId["tract"] for isolatedStarCatalog in inputs["isolatedStarCatalogs"]]
2410 )
2411 for tract in inputIsolatedStarCatalogTracts:
2412 if tract not in inputIsolatedStarSourceTracts:
2413 raise RuntimeError(f"tract {tract} in isolated_star_cats but not isolated_star_sources")
2414 inputs["isolatedStarSources"] = np.array(
2415 [inputs["isolatedStarSources"][t] for t in inputIsolatedStarSourceTracts.argsort()]
2416 )
2417 inputs["isolatedStarCatalogs"] = np.array(
2418 [inputs["isolatedStarCatalogs"][t] for t in inputIsolatedStarSourceTracts.argsort()]
2419 )
2421 refConfig = LoadReferenceObjectsConfig()
2422 if self.config.applyRefCatProperMotion:
2423 refConfig.requireProperMotion = True
2424 refObjectLoader = ReferenceObjectLoader(
2425 dataIds=[ref.datasetRef.dataId for ref in inputRefCatRefs],
2426 refCats=inputs.pop("referenceCatalog"),
2427 config=refConfig,
2428 log=self.log,
2429 )
2430 if self.config.useColor:
2431 colorCatalog = inputs.pop("colorCatalog")
2432 else:
2433 colorCatalog = None
2435 try:
2436 output = self.run(
2437 **inputs,
2438 instrumentName=instrumentName,
2439 refObjectLoader=refObjectLoader,
2440 colorCatalog=colorCatalog,
2441 )
2442 except pipeBase.AlgorithmError as e:
2443 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
2444 # No partial outputs for butler to put
2445 raise error from e
2447 for outputRef in outputRefs.outputWcs:
2448 visit = outputRef.dataId["visit"]
2449 butlerQC.put(output.outputWcss[visit], outputRef)
2450 butlerQC.put(output.outputCatalog, outputRefs.outputCatalog)
2451 butlerQC.put(output.starCatalog, outputRefs.starCatalog)
2452 if self.config.saveModelParams:
2453 butlerQC.put(output.modelParams, outputRefs.modelParams)
2454 if self.config.saveCameraModel:
2455 butlerQC.put(output.cameraModelParams, outputRefs.outputCameraModel)
2456 if self.config.saveCameraObject:
2457 butlerQC.put(output.camera, outputRefs.camera)
2458 if self.config.useColor:
2459 butlerQC.put(output.colorParams, outputRefs.dcrCoefficients)
2460 if output.partialOutputs:
2461 e = RuntimeError("Some visits were dropped because data was insufficient to fit model.")
2462 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
2463 raise error from e
2465 def run(
2466 self,
2467 inputVisitSummaries,
2468 isolatedStarSources,
2469 isolatedStarCatalogs,
2470 instrumentName="",
2471 refEpoch=None,
2472 refObjectLoader=None,
2473 inputCameraModel=None,
2474 colorCatalog=None,
2475 inputCamera=None,
2476 ):
2477 """Run the WCS fit for a given set of visits
2479 Parameters
2480 ----------
2481 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
2482 List of catalogs with per-detector summary information.
2483 isolatedStarSources : `list` [`DeferredDatasetHandle`]
2484 List of handles pointing to isolated star sources.
2485 isolatedStarCatalog: `list` [`DeferredDatasetHandle`]
2486 List of handles pointing to isolated star catalogs.
2487 instrumentName : `str`, optional
2488 Name of the instrument used. This is only used for labelling.
2489 refEpoch : `float`, optional
2490 Epoch of the reference objects in MJD.
2491 refObjectLoader : instance of
2492 `lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader`,
2493 optional
2494 Reference object loader instance.
2495 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
2496 Parameters to use for the device part of the model.
2497 colorCatalog : `lsst.afw.table.SimpleCatalog`
2498 Catalog containing object coordinates and magnitudes.
2499 inputCamera : `lsst.afw.cameraGeom.Camera`, optional
2500 Camera to be used as template when constructing new camera.
2502 Returns
2503 -------
2504 result : `lsst.pipe.base.Struct`
2505 ``outputWcss`` : `list` [`lsst.afw.table.ExposureCatalog`]
2506 List of exposure catalogs (one per visit) with the WCS for each
2507 detector set by the new fitted WCS.
2508 ``fitModel`` : `wcsfit.WCSFit`
2509 Model-fitting object with final model parameters.
2510 ``outputCatalog`` : `pyarrow.Table`
2511 Catalog with fit residuals of all sources used.
2512 ``starCatalog`` : `pyarrow.Table`
2513 Catalog with best-fit positions of the objects fit.
2514 ``modelParams`` : `dict`
2515 Parameters and covariance of the best-fit WCS model.
2516 ``cameraModelParams`` : `dict` [`str`, `np.ndarray`]
2517 Parameters of the device part of the model, in the format
2518 needed as input for future runs.
2519 ``colorParams`` : `dict` [`int`, `np.ndarray`]
2520 DCR parameters fit in RA and Dec directions for each visit.
2521 ``camera`` : `lsst.afw.cameraGeom.Camera`
2522 Camera object constructed from the per-detector model.
2523 """
2524 self.log.info("Gather instrument, exposure, and field info")
2526 # Get information about the extent of the input visits
2527 fields, fieldRegions = self._prep_sky(inputVisitSummaries)
2529 # Get RA, Dec, MJD, etc., for the input visits
2530 exposureInfo, exposuresHelper, extensionInfo, instruments = self._get_exposure_info(
2531 inputVisitSummaries, fieldRegions=fieldRegions
2532 )
2534 self.log.info("Load associated sources")
2535 medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format="decimalyear").mjd
2536 allRefObjects, allRefCovariances = {}, {}
2537 for f, fieldRegion in fieldRegions.items():
2538 refObjects, refCovariance = self._load_refcat(
2539 refObjectLoader, extensionInfo, epoch=medianEpoch, region=fieldRegion
2540 )
2541 allRefObjects[f] = refObjects
2542 allRefCovariances[f] = refCovariance
2544 associations, sourceDict = self._associate_from_isolated_sources(
2545 isolatedStarSources, isolatedStarCatalogs, extensionInfo, allRefObjects
2546 )
2548 self.log.info("Fit the WCSs")
2549 # Set up a YAML-type string using the config variables and a sample
2550 # visit
2551 inputYaml, mapTemplate = self.make_yaml(
2552 inputVisitSummaries[0],
2553 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
2554 )
2556 # Set the verbosity level for WCSFit from the task log level.
2557 # TODO: DM-36850, Add lsst.log to gbdes so that log messages are
2558 # properly propagated.
2559 loglevel = self.log.getEffectiveLevel()
2560 if loglevel >= self.log.WARNING:
2561 verbose = 0
2562 elif loglevel == self.log.INFO:
2563 verbose = 1
2564 else:
2565 verbose = 2
2567 # Set up the WCS-fitting class using the source matches from the
2568 # isolated star sources plus the reference catalog.
2569 wcsf = wcsfit.WCSFit(
2570 fields,
2571 instruments,
2572 exposuresHelper,
2573 extensionInfo.visitIndex,
2574 extensionInfo.detectorIndex,
2575 inputYaml,
2576 extensionInfo.wcs,
2577 associations.sequence,
2578 associations.extn,
2579 associations.obj,
2580 sysErr=self.config.systematicError,
2581 refSysErr=self.config.referenceSystematicError,
2582 usePM=self.config.fitProperMotion,
2583 verbose=verbose,
2584 )
2586 # Add the science and reference sources
2587 self._add_objects(wcsf, sourceDict, extensionInfo)
2588 for f in fieldRegions.keys():
2589 self._add_ref_objects(
2590 wcsf, allRefObjects[f], allRefCovariances[f], extensionInfo, fieldIndex=-1 * f
2591 )
2592 if self.config.useColor:
2593 self._add_color_objects(wcsf, colorCatalog)
2595 # Do the WCS fit
2596 try:
2597 wcsf.fit(
2598 reserveFraction=self.config.fitReserveFraction,
2599 randomNumberSeed=self.config.fitReserveRandomSeed,
2600 clipThresh=self.config.clipThresh,
2601 clipFraction=self.config.clipFraction,
2602 )
2603 except RuntimeError as e:
2604 if "Cholesky decomposition failed" in str(e):
2605 raise CholeskyError() from e
2606 else:
2607 raise
2609 self.log.info("WCS fitting done")
2611 outputWcss, cameraParams, colorParams, camera, partialOutputs = self._make_outputs(
2612 wcsf,
2613 inputVisitSummaries,
2614 exposureInfo,
2615 mapTemplate,
2616 extensionInfo,
2617 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
2618 inputCamera=(inputCamera if self.config.buildCamera else None),
2619 )
2620 outputCatalog = wcsf.getOutputCatalog()
2621 outputCatalog["exposureName"] = np.array(outputCatalog["exposureName"])
2622 outputCatalog["deviceName"] = np.array(outputCatalog["deviceName"])
2624 starCatalog = wcsf.getStarCatalog()
2625 modelParams = self._compute_model_params(wcsf) if self.config.saveModelParams else None
2627 return pipeBase.Struct(
2628 outputWcss=outputWcss,
2629 fitModel=wcsf,
2630 outputCatalog=outputCatalog,
2631 starCatalog=starCatalog,
2632 modelParams=modelParams,
2633 cameraModelParams=cameraParams,
2634 colorParams=colorParams,
2635 camera=camera,
2636 partialOutputs=partialOutputs,
2637 )
2639 def _prep_sky(self, inputVisitSummaries):
2640 """Cluster the input visits into disjoint groups that will define
2641 separate fields in the astrometric fit, and, for each group, get the
2642 convex hull around all of its component visits.
2644 The groups are created such that each visit overlaps with at least one
2645 other visit in the same group by the `visitOverlap` amount, which is
2646 calculated as a fraction of the field-of-view radius, and no visits
2647 from separate groups overlap by more than this amount.
2649 Paramaters
2650 ----------
2651 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
2652 List of catalogs with per-detector summary information.
2654 Returns
2655 -------
2656 fields : `wcsfit.Fields`
2657 Object with field information.
2658 fieldRegions : `dict` [`int`, `lsst.sphgeom.ConvexPolygon`]
2659 Dictionary of regions encompassing each group of input visits,
2660 keyed by an arbitrary index.
2661 """
2662 allDetectorCorners = []
2663 mjds = []
2664 radecs = []
2665 radii = []
2666 for visSum in inputVisitSummaries:
2667 detectorCorners = [
2668 lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees).getVector()
2669 for (ra, dec) in zip(visSum["raCorners"].ravel(), visSum["decCorners"].ravel())
2670 if (np.isfinite(ra) and (np.isfinite(dec)))
2671 ]
2672 if len(detectorCorners) == 0:
2673 # Skip this visit if none of the detectors have finite ra/dec
2674 # corners, which happens when the WCSs are missing. The visit
2675 # will get formally dropped elsewhere.
2676 continue
2677 allDetectorCorners.append(detectorCorners)
2679 # Get center and approximate radius of field of view
2680 boundingCircle = lsst.sphgeom.ConvexPolygon.convexHull(detectorCorners).getBoundingCircle()
2681 center = lsst.geom.SpherePoint(boundingCircle.getCenter())
2682 ra = center.getRa().asDegrees()
2683 dec = center.getDec().asDegrees()
2684 radecs.append([ra, dec])
2685 radius = boundingCircle.getOpeningAngle()
2686 radii.append(radius)
2688 obsDate = visSum[0].getVisitInfo().getDate()
2689 obsMJD = obsDate.get(obsDate.MJD)
2690 mjds.append(obsMJD)
2692 # Find groups of visits where any one of the visits overlaps another by
2693 # a given fraction of the field-of-view radius.
2694 distance = self.config.visitOverlap * np.median(radii)
2695 clustering = AgglomerativeClustering(
2696 distance_threshold=distance.asDegrees(), n_clusters=None, linkage="single"
2697 )
2698 clusters = clustering.fit(np.array(radecs))
2700 medianMJD = np.median(mjds)
2701 medianEpoch = astropy.time.Time(medianMJD, format="mjd").decimalyear
2703 fieldNames = []
2704 fieldRAs = []
2705 fieldDecs = []
2706 epochs = []
2707 fieldRegions = {}
2709 for i in range(clusters.n_clusters_):
2710 fieldInd = clusters.labels_ == i
2711 # Concatenate the lists of all detector corners that are in this
2712 # field
2713 fieldDetectors = sum([allDetectorCorners[f] for f, fInd in enumerate(fieldInd) if fInd], [])
2714 hull = lsst.sphgeom.ConvexPolygon.convexHull(fieldDetectors)
2715 center = lsst.geom.SpherePoint(hull.getCentroid())
2716 ra = center.getRa().asDegrees()
2717 dec = center.getDec().asDegrees()
2719 fieldRegions[i] = hull
2720 fieldNames.append(str(i))
2721 fieldRAs.append(ra)
2722 fieldDecs.append(dec)
2723 # Use the same median epoch for all fields so that the final object
2724 # positions are calculated for the same epoch.
2725 epochs.append(medianEpoch)
2727 fields = wcsfit.Fields(fieldNames, fieldRAs, fieldDecs, epochs)
2729 return fields, fieldRegions
2731 def _associate_from_isolated_sources(
2732 self, isolatedStarSourceRefs, isolatedStarCatalogRefs, extensionInfo, refObjects
2733 ):
2734 """Match the input catalog of isolated stars with the reference catalog
2735 and transform the combined isolated star sources and reference source
2736 into the format needed for gbdes.
2738 Parameters
2739 ----------
2740 isolatedStarSourceRefs : `list` [`DeferredDatasetHandle`]
2741 List of handles pointing to isolated star sources.
2742 isolatedStarCatalogRefs: `list` [`DeferredDatasetHandle`]
2743 List of handles pointing to isolated star catalogs.
2744 extensionInfo : `lsst.pipe.base.Struct`
2745 Struct containing properties for each extension (visit/detector).
2746 ``visit`` : `np.ndarray`
2747 Name of the visit for this extension.
2748 ``detector`` : `np.ndarray`
2749 Name of the detector for this extension.
2750 ``visitIndex` : `np.ndarray` [`int`]
2751 Index of visit for this extension.
2752 ``detectorIndex`` : `np.ndarray` [`int`]
2753 Index of the detector for this extension.
2754 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
2755 Initial WCS for this extension.
2756 ``extensionType`` : `np.ndarray` [`str`]
2757 "SCIENCE" or "REFERENCE".
2758 refObjects : `dict`
2759 Dictionary of dictionaries containing the position and error
2760 information of reference objects.
2762 Returns
2763 -------
2764 associations : `lsst.pipe.base.Struct`
2765 Struct containing the associations of sources with objects.
2766 sourceDict : `dict` [`int`, [`int`, [`str`, `list` [`float`]]]]
2767 Dictionary containing the source centroids for each visit.
2768 """
2769 sequences = []
2770 extensions = []
2771 object_indices = []
2773 sourceColumns = ["x", "y", "xErr", "yErr", "ixx", "ixy", "iyy", "obj_index", "visit", "detector"]
2774 catalogColumns = ["ra", "dec"]
2776 sourceDict = dict([(visit, {}) for visit in np.unique(extensionInfo.visit)])
2777 for visit, detector in zip(extensionInfo.visit, extensionInfo.detector):
2778 sourceDict[visit][detector] = {"x": [], "y": [], "xCov": [], "yCov": [], "xyCov": []}
2780 for isolatedStarCatalogRef, isolatedStarSourceRef in zip(
2781 isolatedStarCatalogRefs, isolatedStarSourceRefs
2782 ):
2783 isolatedStarCatalog = isolatedStarCatalogRef.get(parameters={"columns": catalogColumns})
2784 isolatedStarSources = isolatedStarSourceRef.get(parameters={"columns": sourceColumns})
2785 if len(isolatedStarCatalog) == 0:
2786 # This is expected when only one visit overlaps with a given
2787 # tract, meaning that no sources can be associated.
2788 self.log.debug(
2789 "Skipping tract %d, which has no associated isolated stars",
2790 isolatedStarCatalogRef.dataId["tract"],
2791 )
2792 continue
2794 # Match the reference stars to the existing isolated stars, then
2795 # insert the reference stars into the isolated star sources.
2796 allVisits = np.copy(isolatedStarSources["visit"])
2797 allDetectors = np.copy(isolatedStarSources["detector"])
2798 allObjectIndices = np.copy(isolatedStarSources["obj_index"])
2799 issIndices = np.copy(isolatedStarSources.index)
2800 for f, regionRefObjects in refObjects.items():
2801 # Use the same matching technique that is done in
2802 # isolatedStarAssociation and fgcmBuildFromIsolatedStars.
2803 with Matcher(
2804 isolatedStarCatalog["ra"].to_numpy(), isolatedStarCatalog["dec"].to_numpy()
2805 ) as matcher:
2806 idx, idx_isoStarCat, idx_refObjects, d = matcher.query_radius(
2807 np.array(regionRefObjects["ra"]),
2808 np.array(regionRefObjects["dec"]),
2809 self.config.matchRadius / 3600.0,
2810 return_indices=True,
2811 )
2813 refSort = np.searchsorted(isolatedStarSources["obj_index"], idx_isoStarCat)
2814 refDetector = np.ones(len(idx_isoStarCat)) * -1
2815 # The "visit" for the reference catalogs is the field times -1.
2816 refVisit = np.ones(len(idx_isoStarCat)) * f * -1
2818 allVisits = np.insert(allVisits, refSort, refVisit)
2819 allDetectors = np.insert(allDetectors, refSort, refDetector)
2820 allObjectIndices = np.insert(allObjectIndices, refSort, idx_isoStarCat)
2821 issIndices = np.insert(issIndices, refSort, idx_refObjects)
2823 # Loop through the associated sources to convert them to the gbdes
2824 # format, which requires the extension index, the source's index in
2825 # the input table, and a sequence number corresponding to the
2826 # object with which it is associated.
2827 sequence = 0
2828 obj_index = allObjectIndices[0]
2829 for visit, detector, row, obj_ind in zip(allVisits, allDetectors, issIndices, allObjectIndices):
2830 extensionIndex = np.flatnonzero(
2831 (extensionInfo.visit == visit) & (extensionInfo.detector == detector)
2832 )
2833 if len(extensionIndex) == 0:
2834 # This happens for runs where you are not using all the
2835 # visits overlapping a given tract that were included in
2836 # the isolated star association task."
2837 continue
2838 else:
2839 extensionIndex = extensionIndex[0]
2841 extensions.append(extensionIndex)
2842 if visit <= 0:
2843 object_indices.append(row)
2844 else:
2845 object_indices.append(len(sourceDict[visit][detector]["x"]))
2846 source = isolatedStarSources.loc[row]
2847 sourceDict[visit][detector]["x"].append(source["x"])
2848 sourceDict[visit][detector]["y"].append(source["y"])
2849 xCov = source["xErr"] ** 2
2850 yCov = source["yErr"] ** 2
2851 xyCov = source["ixy"] * (xCov + yCov) / (source["ixx"] + source["iyy"])
2852 # TODO: add correct xyErr if DM-7101 is ever done.
2853 sourceDict[visit][detector]["xCov"].append(xCov)
2854 sourceDict[visit][detector]["yCov"].append(yCov)
2855 sourceDict[visit][detector]["xyCov"].append(xyCov)
2856 if obj_ind != obj_index:
2857 sequence = 0
2858 sequences.append(sequence)
2859 obj_index = obj_ind
2860 sequence += 1
2861 else:
2862 sequences.append(sequence)
2863 sequence += 1
2865 associations = pipeBase.Struct(extn=extensions, obj=object_indices, sequence=sequences)
2866 return associations, sourceDict
2868 def _add_objects(self, wcsf, sourceDict, extensionInfo):
2869 """Add science sources to the wcsfit.WCSFit object.
2871 Parameters
2872 ----------
2873 wcsf : `wcsfit.WCSFit`
2874 WCS-fitting object.
2875 sourceDict : `dict`
2876 Dictionary containing the source centroids for each visit.
2877 extensionInfo : `lsst.pipe.base.Struct`
2878 Struct containing properties for each extension (visit/detector).
2879 ``visit`` : `np.ndarray`
2880 Name of the visit for this extension.
2881 ``detector`` : `np.ndarray`
2882 Name of the detector for this extension.
2883 ``visitIndex` : `np.ndarray` [`int`]
2884 Index of visit for this extension.
2885 ``detectorIndex`` : `np.ndarray` [`int`]
2886 Index of the detector for this extension.
2887 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
2888 Initial WCS for this extension.
2889 ``extensionType`` : `np.ndarray` [`str`]
2890 "SCIENCE" or "REFERENCE".
2891 """
2892 for visit, visitSources in sourceDict.items():
2893 # Visit numbers equal or below zero connote the reference catalog.
2894 if visit <= 0:
2895 # This "visit" number corresponds to a reference catalog.
2896 continue
2898 for detector, sourceCat in visitSources.items():
2899 extensionIndex = np.flatnonzero(
2900 (extensionInfo.visit == visit) & (extensionInfo.detector == detector)
2901 )[0]
2903 d = {
2904 "x": np.array(sourceCat["x"]),
2905 "y": np.array(sourceCat["y"]),
2906 "xCov": np.array(sourceCat["xCov"]),
2907 "yCov": np.array(sourceCat["yCov"]),
2908 "xyCov": np.array(sourceCat["xyCov"]),
2909 }
2910 wcsf.setObjects(
2911 extensionIndex,
2912 d,
2913 "x",
2914 "y",
2915 ["xCov", "yCov", "xyCov"],
2916 defaultColor=self.config.referenceColor,
2917 )
2920class GbdesGlobalAstrometricMultibandFitConnections(
2921 GbdesGlobalAstrometricFitConnections,
2922 dimensions=("instrument",),
2923):
2924 outputCatalog = pipeBase.connectionTypes.Output(
2925 doc=(
2926 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2927 "plane coordinates and chisq values."
2928 ),
2929 name="gbdesGlobalAstrometricMultibandFit_fitStars",
2930 storageClass="ArrowNumpyDict",
2931 dimensions=("instrument",),
2932 )
2933 starCatalog = pipeBase.connectionTypes.Output(
2934 doc=(
2935 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2936 "fitProperMotion is True."
2937 ),
2938 name="gbdesGlobalAstrometricMultibandFit_starCatalog",
2939 storageClass="ArrowNumpyDict",
2940 dimensions=("instrument",),
2941 )
2942 modelParams = pipeBase.connectionTypes.Output(
2943 doc="WCS parameters and covariance.",
2944 name="gbdesGlobalAstrometricMultibandFit_modelParams",
2945 storageClass="ArrowNumpyDict",
2946 dimensions=("instrument",),
2947 )
2950class GbdesGlobalAstrometricMultibandFitConfig(
2951 GbdesAstrometricFitConfig,
2952 pipelineConnections=GbdesGlobalAstrometricMultibandFitConnections,
2953):
2954 """Configuration for the GbdesGlobalAstrometricMultibandFitTask"""
2956 pass
2959class GbdesGlobalAstrometricMultibandFitTask(GbdesGlobalAstrometricFitTask):
2960 """Calibrate the WCS across multiple visits in multiple filters and
2961 multiple fields using the GBDES package.
2962 """
2964 ConfigClass = GbdesGlobalAstrometricMultibandFitConfig
2965 _DefaultName = "gbdesAstrometricMultibandFit"